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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 539.3

ТЕСТИРОВАНИЕ ДВУМЕРНЫХ КВАДРАТИЧНЫХ КОНЕЧНЫХ ЭЛЕМЕНТОВ ПРИ МОДЕЛИРОВАНИИ ИЗГИБА ПЛОСКИХ СТЕРЖНЕЙ

БОЛЬШОЙ КРИВИЗНЫ

© 2011 г. П.П. Гайджуров, Э.Р. Исхакова

Южно-Российский государственный технический университет (Новочеркасский политехнический институт)

South-Russian State Technical University (Novocherkassk Polytechnic Institute)

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

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

Two two-dimensional square-law final elements, accordingly based on classical isoparametrical technology and a method of double approximation of deformations are developed. On test examples about a flat cross-section bend of a thin-walled curvilinear core convergence and accuracy of offered final elements is numerically investigated.

Keywords: the finite element method; flat problem of the theory of elasticity; components of a tensor of deformations; row Maklorena; stiffness matrix; accuracy and convergence.

В работе [1] предложен новый билинейный конечный элемент (КЭ) с полилинейным восполнением геометрии и перемещений, имеющий более высокую скорость сходимости к точному решению, чем элемент аналогичного типа, построенный по классической изопараметрической технологии (ИТ). Вместе с тем в данной работе указывается на то, что разработанный билинейный КЭ имеет определенные ограничения по точности при моделировании изгибных деформаций двумерных объектов с отношением длины к высоте поперечного сечения I /h > 10. Причем для вытянутых в осевом направлении криволинейных тел скорость сходимости оказывается ниже, чем для прямолинейных тел с аналогичным значением параметра I / h .

Целью настоящей работы является разработка универсальных двумерных квадратичных КЭ повышенной точности, позволяющих рассчитывать как объекты с концентрациями напряжений, так и тонкостенные конструкции произвольной геометрической формы с I / h > 10, работающими на изгиб.

Рассмотрим плоский восьмиузловой конечный элемент (рис. 1), отнесенный к базисной (глобальной) декартовой системе осей Zl, I =1, 2. Введем также

местную (локальную) Х1 систему координат, оси

которой в частном случае совпадают с соответствующими направлениями глобальных координат. Локаль-

ные координаты в узлах и на сторонах КЭ принимают значения х, = + 1.

Z".

x„ = +1

(- »-(1 X2

1- X1 " •-1»

x2=

3

6 + I I

2

>

21

Рис. 1

Получим выражения для компонент тензора деформаций е ^ , /, у =1, 2 в глобальном базисе на основании принципов, изложенных в работе [1].

Уравнение для определения компонент е ^ имеет

вид

1

8 2

д и{ д u - + -

д x,

V j

д Xi

(1)

Зададим поле перемещений и 1 в глобальных осях в виде неполного квадратичного полинома

_ 10 ,01 ,11 20 2 Uf — + Xj + X 2 + X1 X 2 + X1 +

02 2 21 2 12 2 + a i X 2 + a^ X1 X 2 + a i x 1 X 2 ,

где а- - коэффициенты, зависящие от принятого

порядка локальной нумерации узлов КЭ. Количество членов полинома (1) соответствует числу узлов рассматриваемого КЭ.

Достроим выражение (2) до полного бикубического полинома, соответствующего шестнадцатиузлово-му КЭ, показанному на рис. 2,

ü i = ul +Д ut

(3)

Здесь

а г. 30 3 . , 03 3 . , 13 3 , , 31 3

Д и1 = Ь1 х1 + Ь1 х 2 + Ь1 х1 х 2 + Ь1 х1 х 2 +

+ Ь 23 2 3 + Ь 32 3 2 + ь 33 3 3 + Ь1 х1 х 2 + Ь1 Х1 Х 2 + Ь1 Х1 Х 2 ,

где Ь-~ - дополнительные коэффициенты.

4

I

121

11

I

1

10 9

.?15"'?16

i I

6 12(x1 , x 2 ) = 6 12 +

5 е (°) (0) . 5612

12 + ^ Х1 + ^ x 2 5 x 1 5 x

5б (°) 12

x+

1

+ — 2

( 5 2 6 (°) ^ 12

5 x 1 5 x 2

12

где е г-°'1 - компоненты тензора деформаций в центре

КЭ (x1 = 0, x 2 = 0).

Для вычисления коэффициентов рядов е

(0)

56

(0)

5 x,

... воспользуемся аппроксимацией перемеще-

ний КЭ в виде

8 , , ч

ui = Е uixv k ( xJ, x 2), l = 1 2,

k=1

(6)

где ик - узловые перемещения КЭ в глобальном базисе; функции формы восьмиузлового квадратичного КЭ

Vк (хьх2) = ^(3р\кр2к -2) х

Рис. 2

Подставив в уравнение (1) полиномы (2) и (3), получим два соответствующих набора компонент тензора деформаций

е 11, е22, е 12 и ец, е22, е 12 в виде отрезков степенных рядов. Оставляя в разложениях е 11, е 22, е 12 и е 11, е 22, е 12 только совпадающие

члены, запишем выражения для деформаций в следующем виде:

_ (0) , (1) , (2) (3) (4) 2 .

е 11 = е 11 +е 11 Х1 + е 11 Х2 + е 11 х1 х2 + е 11 х2 ;

Х[ PlkP 22 k (1 + P1kx1) (1 + P 2 kx 2 ) "

" P 2 k (1" x 2)(1 + P 2 kx 2 ) " P fk (1" x 2)(1 + P 1kx1)] :

р1 к и р 2 к - элементы матрицы локальных координат узлов КЭ в соответствии с принятой на рис. 1 нумерацией

[ P ]■■

(2x8)

-111 -10 10 -1 -1 -11 1 -10 1 0

= (0) + (1) + (2) + (3) + (4) 2 . (4)

6 22 = 6 22 +6 22 x1 + 6 22 x2 + 6 22 x1 x2 + 6 22 x1 ' (4)

(0) (1) (2) (3)

6 12 = 6 12 +6 12 x1 + 6 12 x2 + 6 12 x1 x2

Здесь введены обозначения для коэффициентов:

^ = а|0; е® = 2а;20; е® = а;1; е® = 2а;21; е™ = а|2 ;

6(0) = 01. „(1) = a 11. 6(2) = 2 02. 6(3) = о„12. 6(4) = 21. 6 22 = a 2 ; 6 22 =a 2 ; 6 22 = 2 a 2 ' 6 22 =2a 2 ' 622 = a2 :

61(02) = 2 (a01 + a20).

(4) = 12 + 21 6 12 = a1 + a 2 '

6(2) = 2 a" + a 220 ; 6$ = a02 +1 a ";

6ii (XJ, x2 ) = 40) +

56 (0) Ss^

x1+ _ x^ +

5 x 1

5x

2

2

1

+ — 2

( 52 6 (0) ^ 5 x 1 5 x 2

( 52 6 (0) ^ 5 x 2

i =1, 2; (5)

Введем векторы-столбцы узловых перемещений

{ w} ={ {u«^}...^} j Г , {üi(k)}={

(16x1) '

k =1, 2,

(k ь={u (k) и 2k)},

, 8 и деформаций {е}={е 11 е22 2е 12} в

глобальных осях Zl. Здесь символ « Т » обозначает

операцию транспонирования.

На основании уравнения (2) и соотношений (4) -(6) выражения для компонент тензора деформаций представим в матричной форме

{ 6} =[ D ]{ w},

(7)

В дальнейшем выражения (4) будем рассматривать как отрезки рядов Маклорена:

где [D] =[[D]1[D]2...[D]8] - блочная матрица;

(3х16)

субматрица - [D]k =-

(3x2)

d1k 0

0 d 2k ; d1k ={ 1 0j[ J ]

d 0 2k d0k

d 2 k ={ 0 1j[ J ] -1{ r k}; d 0k ={ 1 0 j[ J ] -1{ Sk }:

^2k ={ 0 1j[ J ] -1{ Sk}.

2

6

5

1

x 1 x 2 +

x

Компоненты векторов { гк } = { г1к г2 к }

r1k = s1k = c [ PlkP 2k (1 + Р 2 kx 2) +

{ } = { s1k s 2 к } определяем по следующим формулам:

г1к = с[(Р 1кР2к " Р13к) + 2Р2к Х1 +

, 3 3 ,3 ,3 2-.

+Pik Р 2k Х2 + Р 2k X1 Х2 + Pik X 2 ];

2 3 3 3 3

r2k = c [( P 1kP 2 k " P 2k ) + PikP 2k X1 +

, ") 2 ,3 ,32-.

+2 P 1kX 2 + PikX1 X 2 + P 2 kX 2];

s1k = c [( PfkP 22k " P 13k ) + 2 P 2k X1 +

3 3 3

+P 1kP 2kX2 + P 2kX1 X2];

2 3 3 3 3

?2k = c [( P 1kP 2k " P 2k ) + P 1kP 2 k X1 +

+2 p 22kxi(1+P 2 kx 2) - Pi3k(1 - x 2)];

r2k = s 2k = c [ P IkP 2k (1 + Plkx1) +

+2 Pl2k X2 (1 + Pik x1 ) - P 2>k (1 - Xl2 )].

С целью исследования сходимости и точности разработанных квадратичных КЭ решим ряд тестовых примеров.

Пример 1. Плоский поперечный изгиб консольного криволинейного стержня постоянного поперечного сечения (рис. 3). Ось стержня имеет круговое очертание радиусом R . К свободному концу стержня приложена сосредоточенная сила Р . Упругие константы материала: модуль упругости Е = 2^105 МПа; коэффициент Пуассона V = 0,25. Размеры поперечного сечения стержня 0,01^0,01 м.

+2 P12kx 2 + P13kX1 Х 2]

где множитель с = 1/4(3р12кр2к -2).

Здесь матрица Якоби в произвольной точке КЭ

[ J ] =

дz, дz-

д x1 д x1 д z 1 д z 2 д x 2 д x 2

Для вычисления элементов матрицы Якоби используем аппроксимацию геометрии КЭ, аналогичную той, которая была принята для аппроксимации перемещений,

8 к , ч

=Е к ( х l, х 2 ), 1 = 1, 2

к=1

к

где г 1 - координаты узлов в осях 21.

Соотношение (7), устанавливающее связь между деформациями и перемещениями КЭ в глобальных координатах, является основой для формирования матрицы жесткости квадратичного КЭ. При этом направление локальных X1 и глобальных осей может не совпадать, т. е. вполне допустимо, когда д / д х{ Ф 1, | д / д Ху | Ф 0, I Ф у .

Отметим, что в литературе, посвященной реализации метода конечных элементов в перемещениях применительно к задачам строительной механики, рассмотренный подход называется методом двойной аппроксимации (МДА). Такое название объясняется тем, что перемещения и деформации в пределах КЭ аппроксимируются раздельно.

При построении матрицы жесткости квадратичного КЭ по ИТ векторы {г к } и {5 к } имеют структуру:

Рис. 3

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

EJ

, а sin2a

P RI---

2 4

Рассмотрим два стержня с длиной дуги окружности I = 0,5 м и I = 1,0 м. Соответствующие значения радиусов равны 0,6366 м и 0,3183 м. Приняв для стержня с I = 0,5 м Р = 100 Н, получаем аналитическое значение перемещения в исследуемой точке V = = 0,01520 м. Аналогично для стержня с I = 1,0 м принимаем Р = 10 Н. Тогда получаем V = 0,01216 м. Для рассматриваемых стержней отношение длины к высоте поперечного сечения I / h составляет 50 и 100, т. е.

I / h > 10, что соответствует относительно тонким стержням.

В процессе тестирования квадратичных КЭ, построенных по ИТ и МДА, установлено, что значительное влияние на скорость сходимости численного решения оказывают применяемая схема интегрирования по Гауссу и способ решения результирующей системы уравнений (РСУ). Вычислительные эксперименты выполнялись только на однослойных схемах дискретизации стержня с постоянным шагом сетки вдоль дуги окружности.

т

и

v

v =

На рис. 4 представлены графики зависимости отношения v / v точ (v - численное значение прогиба,

v точ - точная величина прогиба) от числа КЭ n для стержня l /h = 50. На этих графиках приняты следующие обозначения для номеров кусочно ломаных линий: 1 и 2 - квадратичные КЭ с порядком интегрирования 3 х 3, построенные соответственно по МДА и ИТ, с последующим решением РСУ предобусловлен-ным методом сопряженных градиентов (МСГ); 3 и 4 -квадратичные КЭ с порядком интегрирования 3 х 3, построенные соответственно по МДА и ИТ, с последующим решением РСУ методом Холецкого (МХ), основанном на LDL1 разложении; 5 и 6 - квадратичные КЭ с порядком интегрирования 2 х 2, построенные по ИТ, с последующим решением РСУ соответственно МСГ и МХ. Как известно, МСГ является итерационным методом, а МХ - прямым методом.

Из приведенных данных видно, что численное решение, полученное на базе МСГ, имеет значительно меньшую скорость сходимости, чем решение, основанное на использовании МХ. Кроме этого на начальных участках линий 1, 2, 5 наблюдается своеобразный «провал», что указывает на неустойчивость вычислительного процесса МСГ на грубых сетках. Причина такого поведения МСГ кроется в структуре глобальной матрицы жесткости, для которой при использовании квадратичных КЭ характерно появление недиагональных элементов, сопоставимых по модулю с диагональными элементами.

0,90,8 0,7 Н 0,6 0,5 i 0,4 0,3 0,2 0,1 i

3 / \1, 2, 5

мости в виде графиков V/Vточ ~п для данного стержня показаны на рис. 5. На этом рисунке линия 1 отображает решение для КЭ, построенного на базе ИТ с порядком интегрирования 3 х 3, линия 2 - для КЭ, построенного с помощью МДА с порядком интегрирования 3 х 3, линия 3 - для КЭ, построенного на базе ИТ с порядком интегрирования 2 х 2.

0,90,8 0,7 0,6 0,5 Н 0,4 0,3 0,2 0,1

20 40 60 80 100 120 п Рис. 4

Квадратичный КЭ, построенный по ИТ в сочетании с усеченным порядком численного интегрирования (линия 6) обладает так называемой «суперсходимостью». Данное явление общеизвестно и объясняется тем, что в точках Гаусса х1 = + , I = 1, 2 точность

вычисления деформаций для квадратичных КЭ максимальна.

В процессе числовых экспериментов для стержня с отношением I / h =100, моделируемого квадратичными КЭ, было установлено, что удовлетворительная монотонная сходимость обеспечивается только при решении РСУ МХ. Результаты исследования сходи-

20 40 60 80 100 120 n

Рис. 5

Сравнивая приведенные на рис. 4 и 5 графики, приходим к выводу, что квадратичный КЭ, базирующийся на МДА, с порядком интегрирования 3 х 3 практически не имеет преимущества по скорости сходимости перед квадратичным КЭ, построенным на основе ИТ, с аналогичным порядком интегрирования.

Для криволинейного стержня с отношением I /h =50 помимо поликвадратичных 8-узловых КЭ применялась также однослойная конечноэлементная модель, образованная из полилинейных 4-узловых КЭ. Результаты исследования сходимости для полилинейных КЭ в виде графиков V /V точ ~ п представлены на

рис. 6. Здесь линия 1 соответствует билинейному КЭ [1], линия 2 - классическому изопараметрическому КЭ, линия 3 - несовместному КЭ. Выражение для аппроксимации перемещений 4-узлового несовмест-

6 , , . ного КЭ имеет вид и1 = ^ и1 у к I х1, х2 I, I = 1, 2,

к=1

где ик - узловые и неузловые перемещения КЭ; функции формы ук(Х1,х2) = -4(1+^1кХ1)(1+^2кХ2) ,

k =1,2,.

,4; у 5( x1, x2) =1 - Xj2; у 6( x1, x2) =1 - x 22.

Номера узлов 4-узлового КЭ совпадают с принятой нумерацией угловых узлов 8-узлового КЭ (рис. 1).

Неузловые перемещения и5 и и 6 (внутренние степени свободы) рассматриваемого КЭ исключаем с помощью известной процедуры статической «конденсации» соответствующих элементов матрицы жесткости.

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

0,80,60,40,2

20 40 60 80 100 120 n

Рис. 6

для криволинеиных стержней v = -

P п R

EJ

3

Рассмотрим кольцо с длиной окружности I = 1 м. Соответствующее значение радиуса равно 0,1592 м. Аналитическое значение перемещения в точке приложения сосредоточенной силы при Р =100 Н составляет V =0,76064 0-2 м.

На рис. 8 представлены графики зависимости отношения V / V точ ~ п для однослойной схемы разбивки

с равномерным шагом по окружности при использо-

вании различных типов КЭ. На этом рисунке приняты следующие обозначения для номеров линий: 1 - квадратичный 8-узловой КЭ PLANE82 комплекса ANSYS; 2 и 3 - тестируемые в настоящей работе квадратичные КЭ, построенные по ИТ, с порядком интегрирования соответственно 3 х 3 и 2 х 2 (РСУ решалась МХ); 4 -билинейный 4-узловой КЭ [1].

Напротив, несовместный КЭ позволяет получить достаточно высокую точность даже на грубых сетках. Отметим, что эффективность несовместного КЭ обусловлена добавлением квадратичных функций формы у 5 и у 6. Практически аналогичные значения прогиба v были получены при использовании элемента PLANE42 комплекса ANSYS с активизированной опцией К2 «Extra displacement shapes» («внешние функции формы для перемещений»).

В процессе вычислительных экспериментов установлено, что при использовании полилинейных 4-узловых КЭ численное решение не зависит от способа решения РСУ.

Пример 2. Плоский поперечный изгиб консольно закрепленного разрезного кольца радиуса R (рис. 7). К свободному концу кольца приложена сосредоточенная сила P . Упругие константы материала: E =2-105 МПа; v = 0,25. Размеры поперечного сечения кольца 0,01x0,01 м.

Аналитическое значение перемещения на свободном конце кольца определяем по известной формуле

Рис. 7

0,9 0,8 0,7 0,6 0,5 0,4 0,3 0,2 0,1

120 n

Рис. 8

Из сравнения представленных данных следует, что элементы PLANE82 и изопараметрический КЭ с порядком интегрирования 2 х 2 являются наиболее эффективными по скорости сходимости.

Литература

1. Гайджуров П.П., Исхакова Э.Р. Билинейный четырехуз-ловой конечный элемент для решения двумерных задач теории упругости // Изв. вузов. Сев.-Кавк. регион. Техн. науки. 2011. № 4. С. 7 - 13.

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

1 сентября 2011 г.

Гайджуров Петр Павлович - д-р техн. наук, профессор, Южно-Российский государственный технический университет (Новочеркасский политехнический институт). Тел. 8-908-180-83-67.

Исхакова Эльвира Рашидовна - студентка, Южно-Российский государственный технический университет (Новочеркасский политехнический институт).

Guyjurov Peter Pavlovich - Doctor of Technical Sciences, professor, South-Russia State Technical University (Novocherkassk Polytechnic Institute). Ph. 8-908-180-83-67.

Iskhakova Elvira Rashidovna - student, South-Russia State Technical University (Novocherkassk Polytechnic Institute).