Применение метода двойной аппроксимации для построения матриц жесткости объемных конечных элементов

Введение. При конечно-элементном моделировании напряженно-деформированного состояния массивных тел используются объемные конечные элементы (КЭ) в форме параллелепипедов (гексаэдров), призм и тетраэдров, построение матриц жесткости которых выполнено, как правило, по изопараметрической технологии [1–5]. Вместе с тем известно, что полилинейные изопараметрические КЭ при использовании однослойной схемы неудовлетворительно моделируют изгибные деформации даже при существенном сгущении сетки [6][7]. Суть данной проблемы состоит в эффекте «locking» («запирания») элемента вследствие так называемой деформации «ложного сдвига» [8][9]. Для «улучшения» изопараметрических КЭ используют аппарат несовместных элементов, созданных путем введения дополнительных неузловых степеней свободы или вспомогательных аппроксимирующих полиномов [8]. Вместе с тем наиболее эффективным способом решения проблемы «заклинивания» КЭ является применение моментной схемы метода конечных элементов, теоретические основы которой были разработаны А.С. Сахаровым [7]. В последующем данный подход получил название метод двойной аппроксимации (МДА) [6]. Концептуально МДА базируется на раздельном представлении функций распределения перемещений и деформаций внутри элемента. Целью настоящего исследования является построение на базе МДА и тестирование новых объемных полилинейных КЭ, позволяющих моделировать поведение различных конструкций при различных видах внешнего воздействия.

Материалы и методы. Рассмотрим семейство объемных КЭ, состоящее из восьмиузлового и шестиузлового элементов, в глобальных декартовых осях , = 1, 2, 3 (рис. 1). Геометрию и перемещения КЭ представим в следующем виде:

;

,

где

, — узловые координаты и перемещения;

— «функции формы», представляющие собой произведение одномерных полиномов Лагранжа первой степени;

, , — локальные в общем случае неортогональные координаты КЭ;

— число узлов элемента.

Для базового восьмиузлового элемента (рис. 1 а) = 8 «функции формы» определяем по формуле:

, (1)

здесь — координаты узлов в локальных осях. Значения задаем в форме матрицы:

.

Зависимость между ковариантными компонентами тензора деформаций в локальном базисе и перемещениями в глобальных осях имеет вид [7]:

где

;

.

Связь между вектором деформаций и вектором узловых перемещений представим в матричной форме:

,

где блочная матрица

;

субматрица

,

.

Выражения для вектор-столбцов рассматриваемых КЭ имеют следующий вид [10]: восьмиузловой элемент (рис. 1 а):

(2)

шестиузловой элемент (рис. 1 б):

. (3)

Здесь введены обозначения:

,

;

.

;

;

.

Выражения для «функций формы» шестиузлового КЭ, полученные на основании полинома (1) с помощью принципа «вырождения», имеют следующий вид:

;

;

;

;

;

.

Формулы (2) и (3) являются основой для построения матриц жесткости рассматриваемых КЭ. Соответствующее программное обеспечение разработано на базе вычислительной платформы Microsoft Visual Studio и компилятора Intel® Parallel Studio XE с встроенным текстовым редактором Intel® Visual Fortran Composer XE. Процессы хранения и обработки глобальной матрицы жесткости реализованы в терминах разреженных матриц [11]. Для визуализации результатов расчетов использована дескрипторная графика компьютерной системы Matlab.

Результаты исследования. Исследование точности и сходимости разработанного конечно-элементного алгоритма выполнили на тестовых примерах, имеющих аналитическое решение. В тестовых примерах приведены численные решения, полученные с использованием разработанных элементов и аналогичного по размерности элемента SOLID185 программного комплекса ANSYS Mechanical [5][11]. Ниже приведены примеры, подобранные таким образом, чтобы в них присутствовало сочетание изгибных деформаций и жестких смещений КЭ.

Результаты сходимости представлены в таблице 1.

В таблице 1 в числителе приведено значение прогиба f, в знаменателе — относительная погрешность .

Визуализация перемещений , полученная на базе МДА и ANSYS на сетке 2×2×64, представлена на рис. 3 и 4. Отметим, что поле вертикальных перемещений, полученное с помощью ANSYS, не отражает зоны с экстремальным значением = –0,00314 м, показанной на рис. 3 стрелкой.

В данном примере особенно отчетливо прослеживается важность учета жестких смещений.

,

где

= 0,00126;

— цилиндрическая жесткость;

В данном примере рассматривалась ¼ часть пластины с учетом условий симметрии. Результаты сходимости в виде графиков для однослойной и двухслойной моделей представлены на рис. 5 и 6. Здесь и далее под слоями подразумевается разбивка пластины на КЭ по толщине. На этих рисунках значения параметра 1, 2, 3, 4 соответствуют сеткам: 4×4, 8×8, 16×16, 32×32. Приведенные графики отображают результаты решения, полученные с помощью ANSYS (линия 1) и МДА (линия 2). Горизонтальная линия, обозначенная цифрой 3, соответствует точному решению.

Из рис. 5 следует, что при однослойной схеме разбивки значения относительной погрешности на сетке 32×32 составляют: SOLID185 — 16 %; МДА — 10,5 %. При использовании двухслойной схемы (рис. 6) на сетке 32×32 имеем: SOLID185 — 36 %; МДА — 2,8 %.

Патч-тестирование КЭ выполнено для схемы разбивки 16×16×2 с искажением сетки (рис. 7). Результаты патч-теста в виде картин распределения прогибов приведены на рис. 8 и 9.

Как видно из рисунков, искажение сетки при использовании SOLID185 приводит к более заметной асимметрии поля , чем при использовании МДА. При этом значение максимального прогиба для КЭ МДА = 0,001263 м совпадает с точным решением.

Пример 3. Круглая пластина, жестко защемленная по контуру и нагруженная равномерно распределенной нагрузкой. Радиус и толщина пластины: = 1 м; = 0,01 м. Механические константы аналогичны данным примера 2.

Точное значение прогиба в центре пластины определяется по формуле [12]:

.

Результаты сходимости в виде графиков для дискретизации ¼ части пластины по варианту 1 при однослойной и двухслойной схемах разбивки представлены соответственно на рис. 11, 12 и рис. 13, 14.

На этих рисунках значения параметра 1, 2, 3, 4 соответствуют сеткам 4×4, 8×8, 16×16, 32×32. Как и в примере 1, линия 1 соответствует решению на базе SOLID185, а линия 2 — МДА. Горизонтальная линия, обозначенная цифрой 3, соответствует точному решению.

Из приведенных графиков следует, что элемент, построенный по моментной схеме на сетке 32×32×2 вариант 2, имеет относительную погрешность 4 %.

Картины визуализации поля распределения вертикальных перемещений для МДА и SOLID185 при радиальной схеме разбивки ¼ пластины (сетка 32×32×2) показаны соответственно на рис. 15 и 16.

Из приведенных рисунков видно, что несмотря на качественное совпадение картин , относительные погрешности для максимального прогиба составляют: МДА — 3,8 %; SOLID185 — 71 %. Такая значительная погрешность при использовании SOLID185 объясняется тем, что разработчики использовали для аппроксимации геометрии и перемещений шестиузлового КЭ «функции формы» аналогичные функциям, примененным для восьмиузлового элемента, т. е. без принципа «вырождения» [7].

Обсуждение и заключение. Построенные на базе метода двойной аппроксимации матрицы жесткостей объемных полилинейный конечных элементов позволяют моделировать напряженно-деформированное состояние строительных конструкций произвольной геометрии при различных видах внешнего воздействия. Принципиальным отличием предлагаемой концепции от ранее известных конечно-элементных технологий является то, что перемещения в данном случае задаются в глобальных координатах, а компоненты тензора деформаций определяются в местных, в общем случае, неортогональных осях.

На тестовых примерах показано, что объемные конечные элементы, построенные по методу двойной аппроксимации, обладают устойчивой сходимостью и успешно конкурируют с элементом аналогичного типа SOLID185 вычислительного комплекса ANSYS Mechanical.

Разработанное математическое обеспечение может быть внедрено в отечественные импортозамещающие программные комплексы, реализующие метод конечных элементов в форме метода перемещений.