Об одном способе построения двумерных 4-узловых и трехмерных 8-узловых конечных элементов для решения задач теории упругости Текст научной статьи по специальности «Физика»

____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Том 155, кн. 3 Физико-математические науки

2013

УДК 519.6:539.3

ОБ ОДНОМ СПОСОБЕ ПОСТРОЕНИЯ ДВУМЕРНЫХ 4-УЗЛОВЫХ И ТРЕХМЕРНЫХ 8-УЗЛОВЫХ КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ ЗАДАЧ ТЕОРИИ УПРУГОСТИ

Д. Т. Чекмарев

Аннотация

Предложен новый способ построения численных схем МКЭ решения двумерных и трехмерных задач механики сплошных сред. Метод основан на проектировании ажурных схем линейных конечных элементов большей размерности на двумерную или трехмерную конечноэлементную сетку. На примере линейной задачи теории упругости рассмотрен процесс построения двумерных 4-узловых и трехмерных 8-узловых конечных элементов. Полученные конечные элементы близки к известным полилинейным элементам и превосходят их по эффективности. Схемы содержат параметры, позволяющие регулировать сходимость численных решений. Указано на возможность применения данного подхода к построению численных схем решения других задач математической физики.

Ключевые слова: метод конечных элементов, ажурная схема, эффективность численных схем, аппроксимация, сходимость численных решений, теория упругости.

Введение

Проблема эффективности численных методов решения многомерных задач механики сплошных сред остается актуальной. Известен эффект экспоненциального роста объема вычислений при увеличении размерности задач. Один из возможных подходов к повышению эффективности - применение конечноэлементных (КЭ) схем на ажурных сетках [1-4]. В ажурных схемах МКЭ расчетные элементы покрывают область задачи с регулярными промежутками, что позволяет на сетках одинаковых размеров заметно уменьшить число элементов и тем самым сократить объем вычислений. Наиболее удачной ажурной схемой МКЭ решения трехмерных задач теории упругости и пластичности представляется схема на основе гексаэд-ральной сетки, в каждой ячейке (гексаэдре) которой располагается один линейный элемент в виде центрального тетраэдра (рис. 1), а остальные 4 тетраэдра в расчетах не участвуют. При этом по сравнению с традиционной схемой линейного КЭ на одинаковых сетках число элементов сокращается в 5 раз, а число узлов - в 2 раза. При лучшей экономичности ажурная схема обладает вторым порядком точности (не хуже, чем у схем на основе линейных и билинейных КЭ). На ряде тестовых задач она также показала лучшую сходимость по сравнению с указанными схемами.

Как отмечалось в [1], удачные ажурные схемы получаются при симметричном расположении линейного элемента (симплекса) внутри n-мерного куба. Так, при n = 2 (треугольник в квадрате) расположение является несимметричным (рис. 2), а при n = 3 (тетраэдр в кубе) - симметричным (центр тетраэдра совпадает с центром куба (рис. 1)). Следующий случай удачной размерности - n = 7 (правильный 8-узловой симплекс может быть симметрично вписан в 7-мерный куб). При решении

150

ОБ ОДНОМ СПОСОБЕ ПОСТРОЕНИЯ ДВУМЕРНЫХ 4-УЗЛОВЫХ...

151

задач «неудачной» размерности эффективные схемы МКЭ могут быть получены путем проектирования ажурной схемы более высокой «удачной» размерности. При этом естественно предположить, что высокое качество схемы при проектировании сохраняется. Пример проектирования трехмерной ажурной схемы МКЭ решения задачи теории упругости в двумерное пространство рассмотрен в [5]. В настоящей работе рассматривается развитие данного подхода на задачу большей размерности.

1. Двумерная схема МКЭ решения задач теории упругости на основе трехмерной ажурной схемы

В данном пункте кратко излагаются результаты, полученные в [5].

Ажурная схема линейного КЭ решения трехмерной статической задачи теории упругости на равномерной ортогональной сетке х1 = xQ + hii, xj = x0 + h-2j, x| = = x0 + hsk может быть записана в следующем конечно-разностном виде (см. [6]):

D11 ui + D12U2 + D13U3 U1 F1

(А + p) D21U1 + D22U2 + D23U3 + р Dд U2 + P F2

D31U1 + D32U2 + D33U3 U3 F3

(1)

аналогичном по форме системе уравнений Ламе

(А + р) grad div u + р Дu + pF = 0. (2)

Здесь А, р - параметры Ламе, p - плотность среды, F - внешняя массовая сила, м - поле перемещений, D ij , D д - сеточные операторы, аппроксимирующие вторые производные и оператор Лапласа. Операторы строятся следующим образом. Базисные операторы, аппроксимирующие первые частные производные в элементе (d+ « д/дхт), записываются в виде

(d1 f )ijk

(d2 f )ijk

(d3 f )ijk Операторы

(d1 f )ijk -(d2 f )ijk : (d3 f )ijk :

1

2hQ 1

2hj 1

2h33

(fi+1,j + 1,k+1 + fi+1,j,k — fi,j + 1,k — fi,j,k + 1)j

(fi+1,j + 1,k+1 + fi,j+1,k — fi+1,j,k — fi,j,k + 1)j

(fi+1,j + 1,k+1 + fi,j,k +1 — fi,j + 1,k — fi+1,j,k )•

1 (-f

2hQ Fi-1,j- -1,k

1 (-f

2M Fi-1,j- -1,k

1 (-f

2h3 Fi-1,j- -1,k

(3)

(4)

152

Д.Т. ЧЕКМАРЕВ

равны сопряженным к (3), взятым со знаком «минус». Операторы Dj определим через суперпозиции операторов (3) и (4): Dj = (d+ d- + d+d-) /2, Da = D11 + + D22 + D33 .

Численную схему решения двумерной задачи теории упругости (плоская деформация) получим путем проектирования трехмерной на плоскость X1OX2 . Полагая, что трехмерная расчетная область имеет вид И х [0, Н3], где И - область в Д2, возьмем один ряд ячеек по координате хз и наложим на решение ограничение из = 0. При этом разностная схема (1) примет вид

(Л + д) (D11U1 + D12U2) + д (D11U1 + D22U1 + D33U1) + pFi — 0 (Л + д) (D21U1 + D22U2) + д (D11U2 + D22U2 + D33U2) + pF2 = 0,

В (5) операторы D j получаются проектированием рассмотренных выше операторов на двумерную сетку: Dj = (d+ d- + d+d-) /2, где

(d+ f )ij

(d+ f )ij

(d+f )ij

(d1 f )ij

(d2 f )ij

(d3 f )ij

1

2Н1

1

2Н2

1

2Н3

1

2Н1

1

2Н2

1

2Н3

(fi+1,j + 1 + fi+1,j — fi,j+1 — fi,j), (fi+1,j + 1 + fi,j+1 — fi+1,j — fi,j), (fi+1,j + 1 + fi,j — fi,j+1 — fi+1,j), (— fi- 1,j-1 — fi-1,j + fi,j-1 + fi,j ), (— fi- 1,j-1 — fi,j-1 + fi-1,j + fi,j ), (— fi-1,j-1 — fi,j + fi,j-1 + fi-1,j)-

(6)

Отметим, что операторы d+ и d- в двумерном случае аппроксимируют с точностью до положительного множителя оператор второй смешанной производной,

д2 —д2 )

дх1дх2 дх1дх2

Размер элемента Н3 по координате х3 в построенной численной схеме превращается в регулируемый параметр схемы. При Н3 = то схема (5) переходит в вариационно-разностную схему [6], не учитывающую моментные составляющие в элементе. Известная схема билинейного КЭ совпадает со схемой (5) в случае квадратной сетки( Н1 = Н2 ) при значении параметра

взятой с разными знаками (соответственно

Н3

N

3

2

1 - 2v

+ 3

Н1,

где v - коэффициент Пуассона.

Таким образом, в результате проектирования в Д2 ажурной схемы получено однопараметрическое семейство численных схем МКЭ решения плоской задачи теории упругости, близких к схеме билинейного КЭ. Настраивая параметр Н3, можно изменять влияние моментной составляющей в элементе и тем самым регулировать сходимость численных решений.

ОБ ОДНОМ СПОСОБЕ ПОСТРОЕНИЯ ДВУМЕРНЫХ 4-УЗЛОВЫХ...

153

2. Трехмерная схема МКЭ решения задач теории упругости на основе 7-мерной ажурной схемы

Обобщая описанный выше подход, построим семейство 8-узловых схем МКЭ решения трехмерных задач теории упругости, близких к схеме 8-узлового полилинейного КЭ.

Рассмотрим 7-мерную ажурную схему.

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

(7)

Все ребра данного симплекса имеют одинаковую длину. Отметим также, что наборы из первых трех координат данной системы векторов образуют множество вершин единичного куба в Д3.

Для построения 7-мерной задачи рассмотрим расширение задачи теории упругости на 7-мерное пространство. Такое расширение очевидно. Определим тензор деформаций в Д7 как обобщение соотношений Коши

0 0 0 0 1 1 1 1

0 0 1 1 0 0 1 1

0 1 0 1 0 1 0 1

0 ■ 0 ■ 1 ■ 1 ■ 1 ■ 1 ■ 0 ■ 0

0 1 0 1 1 0 1 0

0 1 1 0 0 1 1 0

0 1 1 0 1 0 0 1

0 (ui,j + Uj,i') ■ i,j 1, • • • ■ 7

£

Связь напряжений и деформаций установим на основе «закона Гука»

@ij A$ij£kk + 2^£ij ■ i,j 1> • • • ■ 7.

В результате 7-мерные уравнения равновесия теории упругости в перемещениях запишутся в виде (2), где градиент, дивергенция и оператор Лапласа определены в Д7. Ажурная схема МКЭ на основе линейного конечного элемента на равномерной ортогональной сетке ж1 = xj + hiii,..., x\7 = x7 + h.717 примет вид, аналогичный (1):

(А + и) Dnui + • • + D17U7 + и Dд Ui + Р Fi

D71U1 + • • + D77U7 U7 F7

0.

(8)

Для построения трехмерной схемы поступим аналогично рассмотренному выше случаю двумерной схемы. Полагая во всех узлах U4 = U5 = и = щ = 0, получим после проектирования в Д3

схему МКЭ в виде

D ii ui + D12U2 + D13U3 Ui

(А + и) D21 Ui + D22U2 + D23U3 + и Dд U2

D3i Ui + D32U2 + D33U3 U3

Ui Fi

+ и (D44 + D55 + D66 + D77) U2 + Р F2

U3 F3

0. (9)

154

Д.Т. ЧЕКМАРЕВ

В (9) операторы Dj — (d+d^- + d+ di ) /2, Da — Du + D22 + D33 выражаются

через базисные операторы

(d1 f )ijk = 4^ (fi+1,j,k + fi+1,j + 1,k + fi+1!j,k+1 + fi+1!j + 1,k+1 —

fi,j,k fi,j+1,k fi,j,k +1 fi!j+1,k+1) ?

(d2 f )ij k

, 7 (fi,j+1,k + fi+1,j+1,k + fi,j + 1,k+1 + fi+1!j+1,k+1

4h2

fi,j,k fi+1,j,k fi,j,k + 1

A+0 ,

j

(d3 f )ijk — 4^ (fi,jjk+1 + fi+1,j,k+1 + fiJj+1Jk+1 + fi+1,j + 1,k+1 —

fi,j,k fi+1,j,k fi,j + 1,k fi+1!j+1,k) ,

(d4 f )ijk

1

4h4

(fi,j,k + fi+1,j + 1,k + fi!j,k+1 + fi+1,j + 1,k + 1 —

fi+ 1,j,k — fi,j+1,k — fi+1,j,k+1

fijj+1,k + 1) , (10)

(d5 f )ijk

TT~ (fi,j,k + fi+1!j,k+1 + fi,j+1,k + fi+1,j + 1,k + 1— 4h5

fi+1,j,k fi!j,k+1 fi+1,j + 1,k fi!j+1,k+1) ?

(d6f )ijk — -ГГ- (fi,j,k + fi,j + 1,k + 1 + fi+1,j,k + fi+1,j + 1,k + 1— 4h6

fi,j+1,k fi!j,k+1

,k

,jjk+1) ,

(d7 f )ijk 4^ (fi+1,j,k + fi,j + 1,k + fi!j,k+1 + fi+1!j + 1,k + 1

fi,j,k

fi+1!j,k+1 fi!j+1,k + 1)

k

и

(d1 f )ijk

— 1 4hi

(fi—1,j,k + fi

1,j-1,k + fi— 1,j,k — 1 + fi—1,j —1,k

1

fi,j,k fi,j— 1,k fi,j,k—1 fi,j— 1,k —1) ,

(d2 f )ijk — (fi,j—1,k + fi—1,j —1,k + fi,j —1,k—1 + fi—1,j—1,k—1 —

4h2

fi,j,k fi — 1,j,k fi,j,k— 1 fi— 1,j,k —1) ?

(d3 f )ijk

— 1

4h3

(fi,j,k —1 + fi—1,j,k—1 + fi,j—1,k —1 + fi—1,j —1,k

1

'fi,j,k fi— 1,j,k fi,j— 1,k fi— 1,j — 1,k) ?

ОБ ОДНОМ СПОСОБЕ ПОСТРОЕНИЯ ДВУМЕРНЫХ 4-УЗЛОВЫХ...

155

-1

(d4 f'jijk 4^ (fi,j,k + fi — 1 — 1,k + fi,j,k — 1 + fi—1,j — 1,k — 1

fi—1,j,k fi,j—1,k fi — 1,j,k — 1 fi,j—1,k—1) ? (11)

-1

(d5 f)ijk — 7T” (fi,j,k + fi—1,j,k—1 + fi,j — 1,k + fi—1,j — 1,k — 1 —

4h5

fi—1,j,k fi,j,k— 1 fi — 1,j—1,k fi,j—1,k — 1) ?

-1

(d6 f)ijk 4^ (fi,j,k + fi,j—1,k — 1 + fi—1,j,k + fi—1,j — 1,k — 1

fi,j — 1,k fi,j,k—1 fi—1,j—1,k fi—1,j,k — 1) ?

(d7 f)ijk 4^ (fi—1,j,k + fi,j — 1,k + fi,j,k — 1 + fi—1,j — 1,k — 1

fi,j,k fi—1,j—1,k fi—1,j,k—1 fi,j — 1,k—1)

Операторы d+, d+, d+, d+ в данном случае аппроксимируют с точностью до положительного множителя следующие операторы дифференцирования:

д2 д2 д2 д3 ;

дх1дх2’ дх1дх3’ дх2дх3’ дх1дх2дх3’

аналогично операторы d4— , d5— , d6— , d7— аппроксимируют операторы дифференцирования

—д2 — д2 — д2 д3

дх1дх2’ дх1дх3 ’ дх2дх3 ’ дх1дх2дх3

Подобные аппроксимации моментных составляющих имеют место и в традиционной схеме полилинейного 8-узлового конечного элемента (см., например, [7]). Последняя совпадает со схемой (9) при h1 — h2 — Ъз и соответствующем наборе параметров h4, h5j, Ъб, h7.

Таким образом, в результате проектирования в R3 ажурной схемы получено 4-параметрическое семейство численных схем МКЭ решения трехмерной задачи теории упругости. Параметры h4, h5, h6, h7 можно настраивать, изменяя влияние моментных составляющих в элементе.

3. Схемы МКЭ на неравномерных и нерегулярных сетках

Проведенный выше конечно-разностный анализ нового класса двумерных 4-узловых и трехмерных 8-узловых схем МКЭ показывает возможность их применения при решении задач теории упругости. Используя технику метода конечных элементов (по ней имеется обширная литература: см., например, [7-11]), данный подход можно распространить на случай использования неравномерных и нерегулярных сеток, а также на геометрически и физически нелинейные задачи теории упругости и пластичности.

Рассмотрим один из аспектов техники МКЭ. Матрица жесткости конечного элемента строится для элемента стандартной формы и размера (единичного квадрата или куба), а для определения узловых сил произвольного элемента используется его отображение на стандартный элемент. При вычислении интегралов с учетом замены переменных необходимо учитывать якобиан преобразования. Как правило, применяется численное интегрирование с использованием гауссовых точек

156

Д.Т. ЧЕКМАРЕВ

(4 точки в двумерном элементе и 8 - в трехмерном). При описанном выше подходе для отображения произвольного 4-угольника на квадрат или гексаэдра на куб естественно использовать линейное отображение тетраэдра на правильный тетраэдр или 7-мерного симплекса на правильный 7-мерный симплекс, вписанный в единичный 7-мерный куб. При этом полученные отображения получаются линейными, а якобиан - постоянным в пределах элемента. Учитывая, что напряжения и деформации (а также их моментные составляющие) в элементе также являются постоянными, достаточно ограничиться одной точкой интегрирования. Таким образом, получается большая экономия в объеме вычислений. В трехмерной задаче вместо 48 значений компонент тензора напряжений (по 6 компонент в 8 точках) нужно вычислить 10 (6 компонент тензора напряжений и 4 момента). В итоге получается почти пятикратная экономия вычислений в элементе. Точность при решении линейных задач не ухудшится. При решении нелинейных (например, упругопластических) задач вопрос остается открытым, поскольку применение гауссовых квадратурных формул обеспечивает высокую точность для функций с высокой гладкостью, чего нельзя гарантировать для поля напряжений в упругопластическом теле.

Двумерные 4-узловые и трехмерные 8-узловые конечные элементы с одной точкой интегрирования широко применяются при решении динамических задач механики деформируемого твердого тела в сочетании с явной схемой «крест» интегрирования по времени. При этом, как правило, возникает проблема неустойчивости типа «песочные часы» (hourglass instability) [6, 12, 13], связанная с неполнотой системы базисных операторов и приводящая к сильным нефизическим искажениям расчетной сетки. Для борьбы с ней используются специальные средства, требующие дополнительных вычислительных затрат. Отметим, что схемы (5) и (9) свободны от данного недостатка, поскольку (например, для трехмерной схемы) система базисных операторов d+, d+, d+, d+, d+, d+, d+, d+ является полной

в том смысле, что пересечение ядер этих операторов является нулевой сеточной функцией. Здесь оператор

(d0f)ijk = ^ (fi+l,j,k + fi+l,j+l,k + fi+l,j,k + l +

+fi+l,j + l,k + l + fi,j,k + fi,j + 1,k + fi,j,k+1 + fi,j+1,k + 1)

аппроксимирует функцию f элементе. Заметим, что поскольку данный базисный оператор в схеме (9) не используется, то при отсутствии кинематических граничных условий задача может иметь бесконечное множество решений, отличающихся постоянным вектором перемещений (параллельным переносом).

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

Заключение

Рассмотренный подход к построению двумерных и трехмерных схем МКЭ на основе ажурных схем более высокой размерности позволяет строить эффективные численные схемы решения задач механики деформируемого твердого тела. Он также может быть применен и к другим задачам математической физики. Описанная в работе идеология может быть применена и к задачам более высокой размерности, чем 3, хотя актуальность численного решения таких задач пока представляется неопределенной.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 13-01-97052-р_поволжье_а).

ОБ ОДНОМ СПОСОБЕ ПОСТРОЕНИЯ ДВУМЕРНЫХ 4-УЗЛОВЫХ...

157

Summary

D. T. Chekmarev. A Method of Constructing 2D 4-Node and 3D 8-Node Finite Elements for Solving Elasticity Problems.

A new method of constructing FEM schemes for solving 2D and 3D problems of continuum mechanics is proposed. The method is based on the projection of the rare mesh schemes of higher dimension linear finite elements on a 2D or 3D finite element mesh. Based on the example of a linear elasticity problem, the construction of 2D 4-node and 3D 8-node FEM schemes is considered. The obtained finite elements are similar to the known multilinear elements and are more efficient. The schemes contain parameters that make it possible to adjust the convergence of the numerical solutions. The possibility of applying this approach to the construction of numerical schemes for other problems of mathematical physics is shown.

Keywords: finite element method, rare mesh scheme, efficiency of numerical schemes, approximation, convergence of numerical solutions, elasticity theory.

Литература

1. Чекмарев Д. Т. Численные схемы метода конечного элемента на «ажурных» сетках // Вопр. атомной науки и техники. Сер. Матем. моделирование физ. процессов. - 2009. -Вып. 2. - С. 49-54.

2. Жидков А.В., Зефиров С.В., Кастальская К.А., Спирин С.В., Чекмарев Д.Т. Ажурная схема численного решения трехмерных динамических задач теории упругости и пластичности // Вестн. Нижегор. гос. ун-та им. Н.И. Лобачевского. - 2011. -№ 4, Ч. 4. - С. 1480-1482.

3. Чекмарев Д. Т. Об эффективности использования узловой информации при решении трехмерных задач методом конечных элементов // Сеточные методы для краевых задач и приложения: Материалы Девятой Всерос. конф. - Казань: Отечество, 2012. -С. 415-420.

4. Жидков А.В., Спирин С.В., Чекмарев Д.Т. Ажурная схема метода конечных элементов решения статических задач теории упругости // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2012. - Т. 154, кн. 4. - С. 26-32.

5. Чекмарев Д. Т, Гладильщикова К.М. Об одном классе двумерных схем МКЭ // Проблемы прочности и пластичности: Межвуз. сб. - Н. Новгород: Нижегор. гос. ун-т, 2006. - Вып. 68. - С. 236-242.

6. Баженов В.Г., Чекмарев Д.Т. Решение задач нестационарной динамики пластин и оболочек вариационно-разностным методом. - Н. Новгород: Изд-во Нижегор. ун-та, 2000. - 118 с.

7. Голованов А.И., Бережной Д.В. Метод конечных элементов в механике деформируемых твердых тел. - Казань: Дас, 2001. - 300 с.

8. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 541 с.

9. Зенкевич О., Морган К. Конечные элементы и аппроксимация. - М.: Мир, 1986. -318 с.

10. Галлагер Р. Метод конечных элементов. Основы. - М.: Мир, 1984. - 428 с.

11. Каплун А.Б., Морозов Е.М., Олферьева М.А. ANSYS в руках инженера. - М.: Еди-ториал УРСС, 2003. - 272 с.

12. Jacquotte O.-P., Oden J. T., Becker E.B. Numerical control of the hourglass instability // Int. J. Numer. Methods Eng. - 1986. - V. 22, No 1. - P. 219-228.

158

Д.Т. ЧЕКМАРЕВ

13. Баженов В.Г., Кибец А.И., Тулинцев О.В. Применение моментной схемы МКЭ для анализа нелинейных трехмерных задач динамики массивных и оболочечных элементов конструкций // Прикладные проблемы прочности и пластичности. Методы решения. - Н. Новгород: Изд-во Нижегор. ун-та, 1991. - Вып. 47. - С. 46-53.

Поступила в редакцию 12.09.13

Чекмарев Дмитрий Тимофеевич - доктор физико-математических наук, профессор кафедры численного моделирования физико-механических процессов, Нижегородский государственный университет им. Н.И. Лобачевского, г. Нижний Новгород, Россия. E-mail: 4ekm@mm.unn.ru