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

УДК 539.3

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

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

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

South-Russian State Technical University (Novocherkassk Polytechnic Institute)

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

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

It's developed flat four-nodal finite element with a multilinear approximation of geometry and displacements, which is intended for the decision ofplane problem of the theory of elasticity. Covariant components of a tensor of deformations in a direction of local axes are presented in the form of pieces of numbers of Maklorena for the purpose of improvement stiffness properties of an element. Shift deformation is necessary to a constant. Convergence and accuracy of an offered final element is investigated on test examples.

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

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

Целью настоящей работы является разработка математического и программного обеспечения, основанного на идее моментной схемы конечных элементов [1], применительно к решению задач о плоской деформации и плоском напряженном состоянии.

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

местную X1 систему координат, такую, что значения

х1 = + 1 в узлах и на сторонах КЭ.

Z.

4

3

2

Х1= 1

Рис. 1

Полагаем, что оси X1 рассматриваемого КЭ совпадают с осями Zl во всех точках. Зададим поле перемещений и1 в глобальных осях в виде полилинейного полинома

(1)

_ 00 10 01 ,11 иi — ai + ai x 1 + ai x. + a j x 1 x.

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

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

Уравнение для определения компонент тензора деформаций в глобальном базисе имеет вид

1

sö 2

д и, д и. - + - '

д x,

V j

д x.

i, j = 1, 2.

(2)

Подставив аппроксимирующий полином (1) в уравнение (2), получим:

_ (0) , (1) 811 — 811 + 8 11 x 2

— (0) + (1) 8 22 — 8 22 + 8 22 x1

— (0) + (1) + (2) 12 — 12 + 12 x 1 + 12 x

(3)

где обозначено: е^-1 = aj0 ; е1(11) = a"; е^^1 = a 21

(0) 1 / 01 , 10 \ е12 = 21 a1 + a 2 j;

(1) — н 8 22 — a 2

(1) 11 812 — 2 a1

(2) —± 11 8 12 — 2 a 2

2

Z

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

й I = ^ + Д и1. (4)

о л и 20 2 . , 02 2 . , 21 2 , , 12 2 ,

Здесь Д и1 = Ь1 х1 + Ь1 х 2 + Ь1 х1 х 2 + Ь1 х1 х 2 +

+Ь22 х12 х2 , где - дополнительные коэффициенты.

Подставив расширенный аппроксимирующий полином (4) в уравнение (2), получим следующие выражения для компонент тензора деформаций:

~ (0) , ~ (1) (2) (3) ~ (4) 2 ~ (5) 2.

611 — 611 + 611 X 2 + 611 Х1 + 611 Х1 X 2 + 611 X 2 + 611 Х1 X 2 ;

6 22 = 6 221 + 6 22 Х1 + 6 22) х2 + 6 221 х1 + 6 22) х1 х2 + 6 У х1 х2;

(5)

- =-(0) + ~ (1) + ~ (2) + ~ (3) + ~ (4) 2 + ~ (5) 2 +

612 = 612 +612 Х1 +612 Х2 +612 Х1 Х2 +6 22 Х1 +612 Х1 Х2 +

+~ (6) 2 +- (7) 2 +6 12 Х 2 + 6 12 Х1 Х 2 ,

где обозначено: е{1; = аг 21

(0) 10 ~(1) 11 w - " е 11 = a1

: (2) _ л L20

: (3) _

= 2 b:

~ (1) = н . ~ 22 = a2 '

~ (5) = 2Ь 22 ■

22 2

:(4) _А12 11

12 = b1

;(5) = 2ь 22

~ 22 = 2 b 2

02

: (2) = 22

~ ^ = 2( a+ a 1

(3) 21

~ 22 = Ь2

Л L 20 .

~ (0) = 01 . ~ 22 = a 2 '

(4) 12

~22 = 2b

2 (a101 + a 20 ):

~i°2) = 2a" + b2

.02.

'1

(22) = 2 a 21 + b (5) _ и 22. ~ (6) ,12

(3) 12 21 ~ 12 = b1 + b 2

: (4) = I b 21 . '12 = 2b1 '

'12 = b1

= , 12 . ~ (7) = ь 22 12 =b 2 ' ~ 12 =b 2 .

~11( x1,x 2) = ~11(0,0) +

5 ~п(0,0) -x

5 x 2

5~22 (0,0)

~22(x1, x 2) =~ 22(0,0) + -

5 x1

~ 12 (X1, x 2) = ~12(0,0),

(7)

где (0,0) - компоненты тензора деформаций в центре КЭ.

Для вычисления коэффициентов рядов 6i]■(0,0) и 56и (0,0)

производных - воспользуемся аппроксима-

5 х1

цией перемещений КЭ в виде произведения одномерных полиномов Лагранжа

1 4 ,

й1 = тЕ й1 (1 + Р1кх1)(1 + Р 2кХ 2 ) , 1 = 1, 2, (8)

4 к—1

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

[ Р ] :

(2x4)

-111 -1 -1 -11 1

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

{И — {{й1(1)}{й1(2)}^{и(4)}}Г ,

~(1) = ~(1) = Ii. ~(0) = ~(0) =_. a

22 22 "2 ' 12 12 2

2( a?1 + a20 ).

Оставляя в разложении (3) только совпадающие коэффициенты, представим выражения для компонент тензора деформаций в окончательном виде:

= (0) + (1) . ~ 11 = ~ 11 +~ 11 Х2 .

• =~ (0) +~ (1) Х . ~ = ~(°) (6) 22 22 22 1 12 12

{u(k)}={u 1k) u ¥)}, к = 1, 2, ..., 4

и деформаций {е}={е11622 2е12} . Здесь символ «Т»

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

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

(9)

Рис. 2

В соответствующих выражениях (3) и (5) совпадают только следующие (подчеркнутые) коэффициен-

: (0) = ~ (0) = 10. (1) = ~ (1) = 11. (0) = ~ (0) = 01. ты: 611 = 611 = а1 ; 611 = 611 = а1 ; 6 22 = 6 22 = а2 ;

{6} = [ D ]{ w},

где [В] = [[ВЫВ] 2 • (3x8) ••[ D ] 4 ] - блочная

субматрица - [ В ]к =1 (3x2) 4 d 1к 0

0 d 2 к d 0 d 0 " 2 к 1к ;

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

d 1к ={ 1 0 }[ J ]

d 2 к ={ 0 1}[ J ]

-1 | Р1к (1 + Р 2 кХ 2)] { Р 2к (1 + Р1кХ1) J

-1 | Р1к (1 + Р 2к Х 2 )| I Р 2 к (1 + Р1кХ1) [

x

6

Z?

5

7

4

1

3

Z

\р1к I ,

[Р 2к ] '

Г0П -1 \ Р 1к | I Р 2 к ] .

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

д Z1 д Z 2

d 0 ={ 1 о}[ j 0] 1

d 0k ={ 0 1}[ J 0] -

[ J ] =

д x 1 д x 1 д Z1 д Z 2 д x 2 д x 2

[J ] - матрица Якоби в центре КЭ.

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

1 4 к

21 =~гТ 21 (1 + Р 1кХ 1)(1 + Р 2 кХ 2 ) , 1 = 1 2,

4

k=1

Zi

Рис. 3

д x,

д м >

координат; ии =-. В выражении (10) принято

д Xi

Производные в центре КЭ, входящие в выражение (7), определяем по формулам:

де п(0,0) 14

д x 2

=т Z мк (,12 + zu p2k)Pik; (11) 4 k=1

д8 22 (0, 0) 14 kt ^ \

-Z-= тЕ Ml ( zl,12 + zl,2 P1k ) P

д x1 4 k=1

2k

д z.

где смешанная частная производная - {■ =-

д х, д

x, д x,

1 J

На основании выражений (7) и (11) субматрица [ D ] к получает вид

[ D ] k =1

(3x2)

d 11 d 12 d 21 d22 d31 d32

к гу

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

Отметим, что полученное соотношение (9), устанавливающее связь между деформациями и перемещениями КЭ в глобальных координатах, является основой для формирования матрицы жесткости КЭ. Рассмотрим случай, когда перемещения и1 и и 2

задаются в глобальных осях Zl (I = 1, 2), а деформации еi]■ (I, ] = 1, 2) определяются в местных в общем случае неортогональных осях Х1 (рис. 3).

3

4

12

d 11 = Р1к [ г 1,1 + ( г 1,12 + г 1,1 Р 2 к ) х 2 ] , d 12 = Р1к [ Г 2,1 + ( Г 2,12 + Г 2,1 Р 2к ) Х 2] , d 21 = Р 2 к [ г 1,2 + ( г 1,12 + г 1,2 Р1к ) х1] , d22 = Р 2к [ Г 2,2 + ( Г 2,12 + Г 2,2 Р1к ) х1] , d 31 = г 1,1 Р 2к + г 1,2 Р 1к , d 32 = г 2,1 Р 2к + г 2,2 Р1к .

Зависимость, связывающая ковариантные еi]■ и физические е(^) компоненты тензора деформаций,

имеет вид е^) -

V g ii g J

г. В этом выражении gii и

gУ] - элементы метрического тензора, образующего

g 11 g 12 g 21 g 22

где

Ковариантные компоненты тензора деформаций еi]■ в осях XI в случае малых перемещений определяем по формуле

еи = 2 ( , ]иИ + 21,'и1,] ), (10)

д

где у = —--компоненты тензора преобразования

симметричную матрицу [ g ] =

(2x2)

gij = 2ц21у (суммирование по повторяющемуся индексу). Если оси Х1 ортогональны, то

[g ] = [g 11 g 22 ].

В соответствии с вариационным принципом метода возможных перемещений формируем матрицу жесткости КЭ

1 1

[ h ] = t II [ D ]т [ Е ][ D ]| J|dx 1 dx 2, (12)

(8x8) -1 -1

где t - толщина пластины; | 31 - определитель матрицы Якоби; [ Е ] - матрица упругости материала

(структура зависит от класса задачи). Интегрирование в выражении (12) выполняем численно с использованием квадратурной формулы Гаусса. При решении задач о плоской деформации считаем t = 1.

Компоненты тензора напряжений в глобальном базисе определяем по формуле

{а} = [ E ][ D ]{ w},

(13)

правило суммирования по повторяющемуся индексу I. где {ст} ={стп ст22 ст12 }т - вектор-столбец напряжений.

1

2

Процедура вычисления напряжений в узлах КЭ представляет следующую последовательность:

1. Формируем массив [ г ] - статически эквива-

(4x3)

лентной нагрузки в узлах КЭ:

1 4

= 7 Е kij (1 + Plnx k)(1 + Р 2nx

4 k=1

n = 1,2,...,4;

при m = i = 1, j = 1; при m = 2^ i = 2, j = 2; при m = 3^ i = 1, j = 2,

1 1

[V ] = J J (ф } {Ф } dx1 dx2 =-

(4x4)

-1 -1

4 2 2 4

1

2

24

1

2 1 2

24

Пример 1. Консольная балка, нагруженная сосредоточенной (рис. 4 а) и равномерно распределенной нагрузкой (рис. 4 б). Упругие константы материала: модуль упругости Е = 2-1011 Н/м2; коэффициент Пуассона V = 0,25. Длина консоли I = 0,1 м. Рассматриваем задачу о плоском напряженном состоянии. Балку моделируем одним слоем КЭ по толщине.

Р =10 кН

к к

где х 1 , х 2 - значения локальных координат в точ-

к

ках интегрирования; ст ^ - компоненты тензора напряжений в точках интегрирования, вычисленные по формуле (13).

2. Формируем массив узловых напряжений [ О ] =[ у ][ г ] , где фундаментальная матрица

(4x3) (4x4) (4x3)

плоского четырехузлового КЭ

!

q = 50 кн/м2

ант

0,01 м

{ф} = {Ф1 ф 2"'ф 4} - вектор-столбец функций формы

КЭ.

Для сглаживания скачков поля напряжений выполняем осреднение напряжений в смежных узлах конечноэлементной сетки.

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

б

Рис. 4

Результаты конечноэлементного решения задачи в

виде значений прогиба на конце консоли для двух видов нагружения, полученные с использованием разработанного билинейного КЭ, классического изо-параметрического КЭ и КЭ PLANE42 комплекса ANSYS, представлены в табл. 1. Значения прогиба приведены в метрах.

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

(нижнее) значение - включенной (include) опции К2.

Таблица 1

r

l

а

l

Расчетная схема Сетка Тип КЭ

Билинейный Изопараметрический ANSYS Аналитическое решение

Рис. 4 а 2 КЭ 0,01770 0,001819 0,001824 0,01894 0,02

4 КЭ 0,01858 0,005645 0,005661 0,01987

8 КЭ 0,01880 0,01190 0,01193 0,02010

16 КЭ 0,01886 0,01646 0,01651 0,02016

Рис. 4 б 2 КЭ 0,003547 0,0003701 0,000351 0,00379 0,00375

4 КЭ 0,003547 0,001083 0,00109 0,00379

8 КЭ 0,003547 0,002248 0,00225 0,00379

16 КЭ 0,003547 0,003098 0,00310 0,00379

Анализируя полученные результаты, приходим к выводу:

- изопараметрический КЭ и элемент PLANE42 с отключенной опцией К2 дают практически одинаковую точность на всех сетках;

- билинейный КЭ превосходит по точности изопараметрический КЭ, но уступает элементу PLANE42 с включенной опцией К2.

С целью дальнейшего исследования возможностей разработанного билинейного КЭ рассмотрим так называемый patch test. Суть данного вычислительного эксперимента состоит в анализе точности на «плохой» конечноэлементной сетке. Схема «плохой» однослойной разбивки консоли на КЭ показана на рис. 5.

Z2, м 0,01 г-

0,02

0,04 0,06 Рис. 5

0,08

Z1, м

fm

P п R EJ

2

В данном тестовом примере, так же как и в примере 1, использована однослойная схема разбивки кольца на КЭ по толщине. В процессе численного эксперимента установлено, что на сетках 36 КЭ и 72 КЭ с использованием билинейного и изопараметриче-ского КЭ матрица жесткости ансамбля КЭ как бы «запирается изнутри», т. е. при отсутствии опорных связей не происходит смещения кольца как жесткого целого, а оно деформируется. На более густых сетках 144 КЭ и 288 КЭ этот эффект исчезает.

P = 100 Н

А-А

Для схемы нагружения, приведенной на рис. 4 а, получены следующие значения максимального прогиба: билинейный КЭ - 0,004455 м; PLANE42 с отключенной опцией К2 - 0, 002357 м; PLANE42 с включенной опцией К2 - 0, 01352 м. Таким образом, на «плохой» сетке билинейный КЭ и элемент PLANE42 с отключенной опцией К2 значительно уступают по точности элементу PLANE42 с включенной опцией К2.

Пример 2. Разрезное кольцо средним радиусом Я = 0,2 м, жестко закрепленное на одном конце и нагруженное сосредоточенной силой на свободном конце (рис. 6). Упругие постоянные материала кольца: Е = 1-1011 Н/м2, V = 0,3. Рассматриваем задачу о плоском напряженном состоянии.

Выражение для точного значения перемещения в вертикальном направлении под силой имеет вид

А

0,01 м

В табл. 2 представлены результаты расчетов /max

в метрах, полученные с помощью билинейного КЭ, изопараметрического КЭ и элемента PLANE42 с отключенной (exclude) и включенной (include) опцией К2.

Рис. 6

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

Пример 3. Растяжение прямоугольной пластины, ослабленной круговым отверстием (задача Кирша). Упругие константы материала пластины: Е = 2-105 МПа; V = 0,3. Учитывая осевую симметрию геометрии и нагрузки, рассмотрим 1/4 часть пластины с соответствующими граничными условиями (рис. 7).

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

Таблица 2

0

Сетка Тип КЭ

Билинейный Изопараметрический ANSYS Аналитическое решение

36 КЭ 0,0009924 0,0006112 0,005645 0,03016

0,03265

72 КЭ 0,004086 0,007331 0,01453

0,03296

144 КЭ 0,02866 0,02167 0,02374

0,03302

288 КЭ 0,03243 0,02574 0,02819

0,03304

q = 180 МПа >

Z2, м 0,12

0,10

0,08

0,06

0,04

0,02

0

Z2, м 0,12

0,10

0,08

0,06

0,04

0,02

0

t = 0,01 м

Рис. 7

!////// /

/////// / У /

/ У

-

0,04 0,08 а

0,04

0,12

0,16 Z1, м

0,08 б

Рис. 8

0,12

0,16 Z1, м

Экстремальные значения суммарного перемещения и ъ тах, нормальных напряжений а х тах и а у т1п ,

касательного напряжения т т1п, интенсивности напряжений а ^ тах , полученные с использованием билинейного КЭ и КЭ PLANE42 с включенной опцией К2, представлены в табл. 3.

Из данных табл. 1 следует, что билинейный КЭ и КЭ PLANE42 на рассматриваемых сетках дают практически одинаковые значения и Е тах . Однако значения напряжений совпадают только в первой значащей цифре.

В данном примере отношение полуширины пластины к радиусу отверстия составляет 3,25. Для классической задачи Кирша указанное отношение обычно принимается больше 10. Поэтому для принятой геометрии полученные значения напряжений а х тах и

а у тт значительно превышают аналитическое решение

а х тах = 3 q = 540 МПа и а у т1п =-д = -180 МПа.

Пример. 4. Растяжение мембраны эллиптической формы, ослабленной эллиптическим отверстием. На рис. 9 показана 1/4 часть мембраны. Упругие константы материала: Е = 2,1105 МПа; V = 0,3. Размеры мембраны: а = 1,75 м; Ь = 1,0 м; с = 2,0 м; < = 1,25 м. Толщина мембраны t = 0,1 м.

q=10 МПа

х 2 У 2 i^c+d ) + (a+b )

=1

Рис. 9

Таблица 3

a

b

Расчетная схема u 2 max•103, м а X max , МПа а Y min , МПа Т min , МПа а i max , МПа

Билинейный КЭ PLANE42 Билинейный КЭ PLANE42 Билинейный КЭ PLANE42 Билинейный КЭ PLANE42 Билинейный КЭ PLANE42

Рис. 8 а 0,2087 0,209 644,3 636 -262,1 -257 -127,6 -147 622,4 637

Рис. 8 б 0,2090 0,209 642,7 638 -263,7 -260 -150,9 -163 629,8 639

Z2, м

1,5 2,0 2,5 Рис. 10

Z1, м

а , МПа

-1,339 -8,20 17,74 27,28 36,81 46,36 55,89 165,43 74,97 84,51 194,05

Задачу о плоском напряженном состоянии мембраны решаем с использованием изопараметрическо-го и билинейного КЭ.

Точное значение напряжения а в точке D (рис. 9)

составляет 92,7 МПа [2]. Численные значения а в

указанной точке, полученные на регулярной сетке с помощью изопараметрического и билинейного КЭ, принимают значения соответственно 93,72 МПа и 94,05 МПа. На рис. 10 показана визуализация поля напряжений а для билинейного КЭ.

Литература

1. Метод конечных элементов в механике твердых тел / под общ. ред. А.С. Сахарова, И. Альтенбаха. Киев, 1982. 480 с.

2 Barlow. J., Davis G.A.O. Selected FE Benchmarks in Structural and Thermal Analysis // NAFEMS Report FEBSTA, Rev. 1, October 1986, Test No. LE1 (modified).

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

21 марта 2011 г.

Гайджуров Петр Павлович - д-р техн. наук, профессор, заведующий кафедрой, Южно-Российский государственный технический университет (Новочеркасский политехнический институт). Тел. (8635)255-4-14.

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

Guyjurov Peter Pavlovich - Doctor of Technical Sciences, professor, head of department, South-Russia State Technical University (Novocherkassk Polytechnic Institute). Ph. (8635)255-4-14.

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