Учет смещения конечного элемента как жесткого целого в смешанной формулировке МКЭ Текст научной статьи по специальности «Физика»

УЧЕТ СМЕЩЕНИЯ КОНЕЧНОГО ЭЛЕМЕНТА КАК ЖЕСТКОГО ЦЕЛОГО В СМЕШАННОЙ ФОРМУЛИРОВКЕ МКЭ

Н.А. ГУРЕЕВА, канд. техн. наук, доцент

Ю.В. КЛОЧКОВ, д-р техн. наук, профессор

А.П. НИКОЛАЕВ, д-р техн. наук, профессор

Волгоградская государственная сельскохозяйственная академия,

400002, г. Волгоград, пр. Университетский, 26.

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

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

КЛЮЧЕВЫЕ СЛОВА: метод конечного элемента, оболочка, изотропная среда.

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

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

1. Геометрия оболочки. Радиус - вектор произвольной точки срединной поверхности оболочки (рис. 1, а) определяется выражением

R = xi + rsindj + rcosdk , (1-1)

^ ^ —

где i, j, k - орты декартовой системы координат; х - осевая координата; в -угловая координата; r - радиус вращения произвольной точки.

Данная работа выполнена при финансовой поддержке администрации Волгоградской области, предоставившей государственный научный грант победителям конкурса в соответствии с Законом Волгоградской области от 4 марта 2005 г. № 1020-ОД "О государственных научных грантах Волгоградской области", постановлением Волгоградской областной Думы от 9 октября 2008 г. N° 8/653 "О выделении государственных научных грантов Волгоградской области", согласно договору от 27 октября 2008г.

Дифференцированием (1.1) определяются базисные векторы в произвольной точке срединной поверхности, касательные к ней

а = R,s = x,s i + r,s sinej + r,s cos вк; (1.2)

52 = R,e = x,e i + r cose j - r sin вк , (1.3)

где 5 - меридиональная координата. Вектор, нормальный к срединной

Рис. 1

поверхности, определяется векторным произведением

A = а х а = -r,s ri + rx,s sine j + rx,s cos вк . (1.4)

Нормаль к срединной поверхности определяется выражением

a = = -ir,s + jx,s sine + rx,s cose. (1.5)

Базис векторов произвольной точки серединной поверхности {a}T = {ax a2 a} можно представить в матричном виде

{a} = [M ]{j} , (1.6)

3х1 3х3 3х1

где {i } = {i j k} - строка ортов декартовой системы координат; [M ] - матрица коэффициентов.

Определим в матричном виде производные базисных векторов произвольной точки срединной поверхности дифференцированием (1.2), (1.3) и (1.5)

Rs } = [N ]{i}; {a,e} = [N2 ]{j}, (1.7)

где {a,s} = {a1 ,s a2 ,s a,s} - вектор-строка производных базисных векторов по

координате 5; {a,e} = {a

2 'в ae

} - вектор-строка производных базисных век-

торов по координате в; [N1], [N2] - матрицы коэффициентов.

С учетом (1.6) выразим производные базисных векторов (1.7) через векторы этого же базиса

{а„} = [*!][М]1 {а} = [ш]{а} ; {а,в} = [Н2][М]1 {а} = [п]{а} . (1.8) 2. Перемещения и деформации. Вектор перемещения произвольной точки серединной поверхности представим компонентами в локальном базисе

{а а2 а} в виде V = о1 ах + о2а2 + оа . (2.1)

Производные вектора перемещения по координатам 5, в, ^ с учетом (1.8) представим выражениями

V = (и,1 +о1шп + о2ш21 + ош31 )а + (оШ12 +°,2 +о2ш22 + от32 )а2 + +(о1ш13 + о2т23 + о,5 +ит33)а = /а + /2а2 + /За ; (2.2)

11 о ) 1 9 9 —*

Vв = (о,в +о пп + о п21 +ип31 )а + (о п12 +и,в +о п22 +оп32)а2 +

12 12 3 ""12 +(рпи + v «23 +V,Q +vri33)a = f2a + f a2 + f>a ; V,t = u,ta + V a + v,ta .

k

б

h

m

I

n

Радиус-вектор точки М1, отстоящей на расстоянии t от срединной поверхности оболочки, определяется выражением

R' = R + ta , (2.3)

а радиус - вектор этой точки в деформированном состоянии оболочки можно

представить в виде R' = R' + V. (2.4)

Векторы, касательные к произвольной поверхности в исходном и деформированном состояниях, определяются дифференцированием (2.3), (2.4)

g = R,i = a + ta,i;g Г = R > = g + V ,í. (2.5)

где под символом i понимаются координаты s, в, t.

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

е„ = gi • V,1 ; = g2 'V,2 ; £33 = g3 ■ V,3 ; e12 = e21 = ^ (g1 V>2 +g2 V1) ;

e13 = e31 = 1 (g1 ■V,3 +g3 ■V,1 ) ; e23 = e32 = 1 (g 2 ■V,3 +g3 ^2 ) , (2.6)

где g1 = (1 + tm31)a1 + tm32a2 + tm33a ; g2 = tn31a1 + (1 + tn32)a2 + tn33a ; g3 = a .

Учитывая (2.2), деформации (2.6) можно представить в матричном виде

м=[ ад. (27)

6x1 6x3 3x1

где {e}T = {en е22 езз 2ei2 2е1з 2егз} - вектор-строка деформаций точки М1; {u}T = {и1 и2 и} - вектор-строка контравариантных компонент перемещений

точки М1; [L] - матрица дифференциальных операторов.

3. Соотношения между деформациями и напряжениями. Закон Гука для изотропной среды имеет вид [1]

М ,-11 . 22 . . М „11.

e11 =-£ g11(CT g11 + ст g 22 + ст g33) + + -)g!1g11CT ;

e22 = - М g22 1 g11 + g 22 + g33 ) + ^ + M) g22 g22^ ; Е E E

M 11 22 33 \ /- 1 M\ 33 /1 M\ 12

e33 = -TT g33(ff g11 + ст g22 + ст g33) + +^g33g33^ ; e12 = e21 = + Т;Х?1^22ст ; Е E E E E

S13 = e31 = O1 +M)g11 g33^13 ; e23 = e32 = O1 + M)g22g33CT23 , (31)

или матричной форме {e} = [ D] {ст}, (3.2)

6x1 6x6 6x1

где {ст}Т = {ст11 ст22 ст33 ст12 ст13 ст23} - вектор - строка контравариантных компонент тензора напряжений в произвольной точке срединной поверхности; [ D ] -

матрица податливости материала.

Функционал Рейснера для оболочки вращения имеет вид [2]

ПК =${ст}Г [L]{u}dV -1 |{ст}г [я]{ст} dv -jj>}T {q}dF , (3.3)

V 2 V F

V - объем оболочки; F - поверхность оболочки с заданными внешними силами.

4. Матрица деформирования конечного элемента. Конечный элемент выбран в форме шестигранника с узлами i, j, k, l, m, n, p, h (рис. 1, б). Для выполнения численного интегрирования глобальные координаты s, в , t аппроксимировались через локальные координаты куба - 1 — ц, Z — 1 (рис.1, в) трилинейными соотношениями 2 = {f {Ау}, (4.1)

где {Лу} ={л' & 1к 1ш 1" 1Р 1И} - вектор-строка узловых значений величины X. Под символом X понимаются координаты 5, в,

Производные глобальных координат в локальной системе 8,п, 8,$ в в ,п, в 1,п, t,t; и локальные координаты в глобальной системе %,в, ц,5, т],в, ц,ь С,5-, С,в, определяются дифференцированием (4.1).

Векторные поля перемещений и тензорные поля напряжений внутренней точки конечного элемента аппроксимируются через узловые неизвестные так же трилинейными соотношениями

V = {/(^,0}г V}; & = {=,0}г {&у}, (4.2)

где V = о151 +о2а2 +оа - вектор перемещения внутренней точки конечного элемента; { Vy } = {V' Vk V1 Vш Vп УР VИ | - строка узловых векторов перемещений; & = а11а1а1 + а22а2а2 + а33аа + 2а12а1а2 + 2сг13а1а + 2а23а2а - тензор напряжений внутренней точки конечного элемента; {&у }т ={ &' &1 &к &1 &ш &п &Р &И}

- строка узловых тензоров напряжений.

С учетом (1.6) выразим базисные векторы узловых точек через базисные векторы произвольной внутренней точки конечного элемента в виде

{аш } = [мш][М\ 1 {5} = [жш]{5},(® = ',к, 1,ш,п, р,И) , (4.3)

3x1 3x3 3x1

где {аш} = {аШ аШ аш }- базисные векторы узловой точки.

Выразим диадные произведения базисных векторов узловой точки

ш)Т (—ш—ш — ш—ш —а—а — ш — ш п.—т-т — а—а\

а } = {а1 а1 а2 а2 а а 2а1 а2 2а1 а 2а2 а }

1x6

через диадные произведения произвольной точки конечного элемента {а} = {а1а1 а2а2 аа 2а1а2 2а1а 2а2а}

1x6

в матричной форме, используя соотношения (4.3),

\аш} = [ош]{а}; (ш = ',],k, 1, ш, п, p, и) , (4.4)

6x1 6x6 6x1

где ^- матрица коэффициентов, получаемая на основе (4.3).

6x6

Выражения тензоров напряжений внутренней и узловых точек конечного элемента в их базисах имеют вид [1]

& = ЩТ {а};& ш = ^ {<г ш }, (ш= ', ],к, 1,ш,п, р,И) . (4.5)

1x6 6x1 1x6 6x1

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

V ЖК }, (46)

8x1 8x24 24x1

где }Т ={о1г о2' о' о1 о2л ол... о1 о2И оИ} - строка скалярных узловых неиз-

1x24

вестных; [О ] - матрица, элементами которой являются базисные векторы узловых точек {аш}. После замены векторов {аш} в матрице [О] по (4.3) вектор перемещения внутренней точки конечного элемента и его производные по координатам 8, в, t можно записать в виде

V=Н7

1x3

/ (£чХ)[*''] / ] -Л (£ч,Х)[*

3x3

3x3

V!;

v ,д-(г!'

1x3

(v, }'

24x1

(4.7)

А,л[ ] ^ ] .../8,л[ n ]

3x3 3x3 3x3

где^ - функции являющиеся элементами аппроксимирующей матрицы (/(%,ч,()}Т.

С учетом (2.1) и (2.2) из соотношений (4.7) получаются матричные выражения для аппроксимации компонент вектора перемещений внутренней точки конечного элемента и их производных:

М = [ а](Уу }; = [/1 ]т (уу }; ^,2 = [/2 Г V}; = [/3 Г У};

(4.8)

^ = [/4 ] т (Уу} ;о,1 = [/5 ] т (Уу}; ^ = [/6 ] т У};»,1 = [/7 ] т У};

°,2 = [/8 ] т (Уу};»„= [/9]т (Уу}. Узловые тензоры напряжений, используя (4.5) представим в матричном виде

{&} = р]{<Гу }, (4.9)

8x48

48x1

где (оу} = ](ог} (о} ■■■(ок} ?- строка узловых неизвестных напряжений;

1x48

1x6 1x6

матрица, элементами которой являются строки диадных произведений

Г-®]г Г3!

а > . После замены строк < а > по (4.4) в матрице т (4.8) тензор напряжений (4.2) внутренней точки конечного элемента можно представить соотношением

= (а}

1x6

/1 (#,ч,Х)[ег]т/2 (£чХ)|У]т(^,ч,С)\ок]т ]}. (4.10)

6x6 6x6 6x6 ] 48x1

Учитывая диадное представление (4.5) тензора 0 , компоненты тензора напряжений внутренней точки конечного элемента на основании (4.9) представим

в виде (о}=[Б](оу} . (4.11)

6x1 6x48 48^

Принимая во внимание соотношения (4.11) и (4.8), запишем функционал (3.3) в виде

Пя = (оу } | [^]Т [^^(оу}-1 (оу }Г _[ [5] ^ [51 &}-

\°у. - ...... - -

1x48 У 48x6 6x24 24x1 2 1x48 У 48x6 6x6 6x48 48x1

-(°у Г ш а] М .

1x24 я 24x3 3x1

(4.12)

Дифференцируя функционал (4.12) по узловым неизвестным (оу} и (^}Г и приравнивая производные нулю, получим систему уравнений

^Пн = в ь}- [н ] (о, }=0; ^ = [в]т (о,}-/ }=0, (4.13)

0(°у} 48x24 24x1 48x48 48x1 С(иу} 24x48 48x1 24x1

где [в]=|[^]Г [В] ¿У ; [Н]=|[^]г [Б][Б] ¿У; [/] = Ц[А]Г (д} <Ь.

; 24x3 3x1

24x1

Представим систему (4.13) в традиционной для МКЭ форме

[H] [Q]

[к] {Z/ } = {fR }, где [k] =

72x1

72x1

48x48 48x24

[Q]" [о]

24x48 24x24.

- матрица деформирования конечного

элемента;

мента

; {z / f =к }Ъ г!

1x72 I 1х 48 1х 24 J

; f }г=\{о}гшг\

72x1 Ux48 1x24 J

вектор узловых неизвестных конечного эле-

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

Для формирования матрицы деформирования всей конструкции используется традиционная процедура МКЭ [3].

Пример 1. Определено напряженно-деформированное состояние кольца, нагруженного сосредоточенной силой P, приложенной в верхней точке вертикального диаметра, при расположении напружиненной опоры в нижней точке диаметра. Ввиду симметрии конструкции рассматривалась половина кольца (рис. 2) при следующих исходных данных:

R = 30см; h = 0,5см; Р = 2Н; /л = 0,3; Е = 2-107Н/см2 .

Рассматриваемая часть кольца разбивалась на 60 элементов вдоль дуги и на два элемента вдоль ее толщины. При изменении жесткости опорной пружины точка 5 получала различные перемещения vn (см). В таблице приведены кольцевые физические напряжения в точках 1, 2, 3, 4 (рис. 2) при различных жесткостях опорной пружины и получаемых соответственно смещениях vn кольца как жесткого целого. В числителе приведены значения напряжений при использовании тензорно-векторной аппроксимации искомых величин, а в знаменателе -при их аппроксимации как скалярных величин.

Анализ результатов, представленных в табл. 1, показывает, что напряжения (в числителе) меняются незначительно даже при смещении кольца как жесткого целого на 111,7 см. Результаты, представленные в знаменателе, оказываются неприемлемыми уже при смещении кольца как жесткого целого на величину 1 см.

Таблица

Рис. 2

72x72

Vn (см) 0,0 0,1 1,0 2,0 10,1 111,7

№ точки Напряжения ст2^ (H/см2)

1 228,46 228,46 228,43 228,40 228,16 225,15

229,48 222,99 164,53 99,43 -427,15 -7001,6

2 -225,97 -225,97 -225,94 -225,90 -225,61 -222,14

-228,27 -221,45 -160,07 -91,72 461,15 7363,9

3 -133,69 -133,69 -133,66 -133,62 -133,32 -129,99

-137,62 -137,23 -169,72 -205,90 -498,580 -4152,8

4 130,29 130,29 130,26 130,22 129,91 126,46

130,87 134,5 162,66 194,53 452,31 3670,7

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

Конечный элемент на основе разработанного способа аппроксимации векторов перемещений и тензоров напряжений может найти применение в про-

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

Л и т е р а т у р а

1. Голованов А.И., Тюленева О.Н., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. - М.: Физматлит, 2006. - 391с.

2. Постное В.А., Хархурим И.Я. Метод конечных элементов в расчетах судовых конструкций. - Л.: Судостроение, 1974. - 344с.

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

4. Рикардс Р.Б. Метод конечных элементов в теории оболочек и пластин. - Рига: Зинанте, 1988. - 284с.

5. Седов Л.И. Механика сплошной среды, т.1. - М.: Наука, 1976. - 536 с.

THE ACCOUNT OF DISPLACEMENT OF THE FINAL ELEMENT AS RIGID WHOLE IN MIXED FORMULATION FEM

Gureeva N.A., Klochkov J.V., Nikolaev A.P.

The volume finite element is developed for calculation of covers of rotation in shape hexahedron with central unknown persons in a kind a component of a vector of moving and a component tensor strains. For expression the component of a vector of moving of an internal point of a finite element and its component tensor strains through central unknown persons was used a way of interpolation vector and tensor fields three linear by form functions that has allowed to consider implicitly displacement of a finite element as the rigid whole. For formation of a matrix of deformation hexahedron the variation principle in the mixed formulation was used.

On an example efficiency of the offered way of approximation of required sizes as vector and tensor fields in comparison with traditional way of approximation of required sizes as scalar fields is shown.

KEY WORDS: finite element method, shell.

Hb -0- -0-