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

УДК 539.3

А.Б. Ефимов, О,В. Аксененко, А.В. Цвелих

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

ABSTRACT

Quadrilateral element by T.H.&.Pi an and Eel Iinger-Re issuer variational principle are used to design algorithm and software for solving axisymmetric problems of incompressible solid states. All the necessary matrices were derived and a computer program is described. The problem of stresses and rigidity analysis of the spherical damper is solved. It is shown that numerical solution is close to the exact one. The error of numerical result is about 8 %.

Широко распространенный метод конечных элементов в Форме метода перемещений [1] непременим для прочностного анализа деталей из несжимаемых материалов (коэффициент Пуассона равен

0 5 >- Причина заключается в том, что невозможно использовать закон Гука в виде:

о^[Т)]£(1 )

:;з-за неограниченного роста элементов матрицы [В] при стремлении коэффициента Пуассона к 0.5- Альтернативным подходом является использование гибридного метода конечных элементов СГМКЭ} [3], позволяющего избежать указанной неприятности. В данной работе предлагается версия ГМКЭ, дающая возможность решать осесимметричную задачу теории упругости без ограничений на величину коэффициента Пуассона. Подход основан на применении четырехугольного изопараметрического элемента, в пределах которого поля напряжений и перемещений аппроксимируются независимо друг от друга. Преимуществом данного подхода является также и возможность задания узлов конечно элементной модели непосредственно на оси симметрии, что невозможно для обычного метода конечных элементов в форме перемещений из-за возникновения неинтегрируемой особенности при вычислении матрицы жесткости [ 1 ].

Применим для решения указанной задачи Функционал Рейссне

ра. Для одного конечного элемента Спри условии согласованности

поля перемещений и на границах с соседними элементами) Функ-

ч

ционал имеет вид:

П = [(--=-аС[3]о> + сгТ[Ь]и ) ЙУ , к 2 я е

(2)

где [Ь] -дифференциально-алгебраическая матрица, задающая геометрические соотношения: е = [Ь] и ,

ч

[3] - матрица в законе Гука е — [8] о.

Для осесимметричной задачи интегрирование по объему можно заменить интегрированием по площади:

П = 2п[г(- ~сТ[3]сг + [Ь]и ) <1А к •* £ я в

(3)

Сиспользуем обычную цилиндрическую систему координат). Тогда векторы и матрицы в П имеют следующий вид:

^ = [а*с'ъ*-г.с'е]Т ; е = [е-£~Гш~ео]Г : «_= •*

Н 2 НИ 0'

Ч ч

Рис. 1. Конечный элемент.

Для четырехугольных изопараметрических элементов рис. 1 сог-эванное поле перемещен

ными базисными Функциями-.

ласованное поле перемещений и будем аппроксимировать билиней-

ч

ЯСЯОНОЯО

12 3 4

О N О Ш О N О к

12 3-

и 1

С|£

[Я зи ч д

де- « - вектор узловых перемещений; [V ] - матрица согласован-

ч ч

!Ь1х базис ных Функций, ее компаненты имеют вид:

N

(6)

I номер узла. Выборе*! ^ - насно [3]> для напряжений ап-жсимаиикт, независимую от V

к 1 0 о 0 1) 0 0 0 ? 0 0 0

и 0 1 о 0 0 ? 0 0 0 Г) 0 0

0 0 і 0 ■у 0 0 о 0 0 к V

'е 0 0 0 1 0 0 V ? 0 0 0 0

э

р

12

= Р + СФЖ

(7)

где /3

С

В этом выражении /3

[р Р Р . . .р ]

5 в' 7 1 12

вектор внутренних параметров напряжений,

задающих постоянную аппроксимацию на элементе, /3 - вектор па-

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

' “х ' г-н 4-у о 1

•V. о о И' -Н1!

= [Шх]\ .

(в)

где X - вектор внутренних параметров перемещений; [Ш^] ~ матрица несогласованных базисных Функций.

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

П

= 2п|г<:-А

[3]<у+суГ [Ь1ы ) )а.Ал~2п$1'ТТикаГ9 , (9)

где Г - граница элемента; X - усилия, действующие на границе в

со стороны соседних элементов.

Вновь введенное слагаемое имеет простой Физический смысл, поскольку - не что иное, как невязка в перемещениях, а поле

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

т т . ,т Т = о [п]

(10)

где [п] - матрица направляющих косинусов нормали участка границы (ЗГ к осям ОД и 02- Матрица [п] имеет вид:

мт =

п г 0

О п

п п

2 г

0 о

(11 )

Представим вектор а в виде суммы постоянной и линейной компонент:

а = а + о.

(12)

Z661 ‘COT м "Посій и wexono ■dHireti'ow xєнэхєн

.'WHhAlfOU

rw иь зхєжохкинЛ инєхєнс ииннєєсі о Eifed аэхни хняомнниОо ear? ох

(91) * + xn[U*0 = ^[Ч]

:ЛииЛэ а онєжоігєесі яхнэ хэжон eoweeJHiro вчхвсіх (Єї) а охь ‘яхихвиеє иігов е'иВЄйтоіги ои ired -лвхни ниро а вoxeAecdgoedи elwHedj ои ігегїлехни ‘ончігехєаоОеіо

Є

У

(ді) ■ о = °ТР f3of«rj;aJ

і л* J

охь ‘вэяхиетеэЛ оноЛсІьвн квинєаосіисиеіни ютннваі.оїУвсіоои©н

а ар

- + ~т &

ze

О

_± +

І Є

= C1J

et/j

(PI) ' 9VP + вГР (*<>[%] )*f =°JP

:eHHdj eirAwdocfc ou weAeedgoedu a ©oweejeiro евнРвігооц

J J

• *JP ''‘nfuj Члі Гі<г - ejP ^nruj^Mjue -

X X •!> Ju

(ЄІ)

-*yp + bnfTJx" + x>fsJx-o-|-= *и

etfaox Ч£/ f®J -4<o э&л

иьєОєє ионкисіхоїчмиоооо ohhoihoj

4 V т т

j г i -p ' fSjfo -f о- fb?« -f &^[L]u^) d4.e-

A« nr;

- 2я JVo1* f « jTti^ dT& .

Г

Нэдл' <’>мым условием стационарности Функционала является с ледуюг^° условие'

Jro-^fn jrii^ £®Р jp-ti * fnj<^h dtT^ . f I8J

Г Г'

Данному условию «удеудовлетворять априори за счет соответствующего рывора аппг‘'-к-'ж-ации для о». Подставляя в полученное

П

выражение а л п рок с имэ*; уал для и <у , получим;

ь {г[к У ы][ф.:фжг1 or-ft. = о

** А. I 11 © Ь

(19)

'Pfi

:я„ i-jr£»x] in][Фх:фнjare

Поскольку з оощем случае Х^О, то [U]ft =О- Тогда

ь

р

п|

(20)

где /З ~[в р р р /, Р =[р р р р ]Т .

1 1 1- 51 б' 7 8 ' II Р* Ю 11 12

Преобразуя последнее выражение, получим

- <21> ш

если только сЗеИ;[М ]&0. Тогда, обозначив а вектор напряжений

11

в случае» если он аппроксимируется не двенадцатью, а восемью

параметрами напряжений р , получим;

* * * * * *

о = сг + о- = р +[ф.]Р=[1:ф.]

С П С П I п

* *

*[Ф № .

(22)

где [I] - единичная матрица; [Ф*]=[Ф ]-[Ф ] [Ы ] * [М ].

Г» Т. II II I

С учетом замены о на » Функционал Рейсснера примет вид:

Пн = 2п$г(- -£-о [3]СГ +<? [Ь]и +0^ [Ъ]^) ОА&. (23)

А

&

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

и В таком случае имеем следующий Функционал: ч

*

П = 2пГг(--^-о'*Г[3]о*+<у*т[1,]и ) <Ы . (24)

к ^ до

А

в

СЗаметим, что удаленное слагаемое должно иметь более высокий порядок малости по сравнению с другими слагаемыми и тождественно равняться нулю без принятия каких-либо допущений в слу“ чае,если конечный элемент представляет собой параллелограм, пара сторон которого параллельна оси ОД- Для элемента, произвольно ориентированного по отношению к глобальным осям координат, этот Факт равенства указанного слагаемого нулю подтвердить или опровергнуть не удалось ввиду исключительной сложности аналитических выкладок).

Условия устойчивости Бабушки-Бреззи для элемента, основанного на с и и , будут следующими: ч

сИт(о' ) > (Цт(и )-г , (25)

ч

где г - число степеней свободы элемента как твердого целого, а

<Ит() ~ размерность величины (), т.е. число параметров этой

величины В случае осевой симметрии г=1, сИт^ ) = В и

<Ит(и )=8, т.е. введенный элемент удовлетворяет условию устой-ч

чивости.

Для задачи определения осесимметричного напряженно-

деформированного состояния некоторого тела Функционал энергии запишется так:

П У П - V

а, к!

(26)

где суммирование производится по всем конечным элементам, а V обозначает потенциал внешней нагрузки:

V = 2п Г грТи ДГ . (27)

■> я Р

Г

р

Здесь Г часть границы, на которой заданы силовые граничные р

условия (кинематические граничные условия будут добавлены в

т

качестве дополнительных). Вектор о = [р р 1 - это распреде-

гх

ленная Сна единицу боковой поверхности) осесимметричная нагрузка.

Условие стационарности Функционала П имеет вид:

<5П - О , (28)

*

где варьирование следует провести отдельно по (3 и и -

я

Преобразуем выражение (24) , подставив в него аппроксима-

*

ции для а и и - Тогда: ч

1 * *т *

пл =-£-/3 [ЮР + Р [С]р . (29)

где [Н] = 2п£г[ Ф* ]Т [в] [Ф* ] &Ав А

[с] = 2пГг[Ф* ]Т[Ь][Я ] ОА

ч е

А

В локальных координатах С?,г?> матрицы [Н1 и [с] примут вид:

[КЗ = т-[Ф*]г [3][ФШ ]\(1еЫ\ <ХОп :

(30)

[с] = 2njj г[ф*]Т[ь][я^]\аеы\ акая .

-1 - 1

Здесь: г = а +а ?+а 7)+а2£т7 ; |еіеі.7| =д.г+6.^ +д.±~0 . (31)

-де

Л1=а2ЬЭ~аЭЬ2: 'г2=а1Ьэ'аЗЬ1: <13=аіЬ2_а2Ь1 *

а1ь1

а2Ь2

а Ь Э 3

а Ь 4 4

-1 1 1-1 Г-*

-*-і ; 1 і 1 і 1

каждая пара Сг >г > дает координаты соответствующей вершины ь * * цемента в глобальной системе координат. Матрица [Ф ]—[11Ф.]

. ь

состоит из двух подматриц, каждая размером 4x4. [I] ная подматрица, а £ф ] имеет вид:

Г|

единич-

[ф*ъ] =

-та г? 2 і

-т V 2 1

-,В31?-те417?

ТЇ-ТО1І?

-шя±п

Г)-те ?

11

О

О

О О

. (32;

гаесь т. ~ этс элемент матрицы [Н ] [Ы ], где матрица [и ]

>.1 II I II

преаставлкма в тасэляце I, а элементы [Н ] таковы;

Таблица 1

8 ъ. 16 . — -- Ь а - . (3. 3 з 4 45 з Ь а 45 2 1 0 8 ~За За4 + ■ + 45 а1а2 16 —а а 9 2 3

б ^ в ^ Т Ь1а2+"9-Ь2а1 0 16 тг а, а 9 12 8 —а а - 3 14 88 - 45 агаз

0 «б —-я- а а 9 2 3 в, ?6 . 3 за4 45 з 88 , - 45 Ь2а1 - — Ь а - 9 з 2 8 . — о а 9 2 3

; /"1 ' 8 -^а1а4- вв - . „ а, а 45 2 э я , , в . Ъ а + о а 9 1 2 9 2 1 в . »б , 3 1а4 45 1 08 . т^“ Ь а 45 2 з

" дъ За2 9 Ь2аз: ш21" 3 “1“4 45 “г“з 45 “1

в

г- а а +

вв

32 з 3 4 45 1 2

а а ; т.

16

88 . >6 ,

Ь а-----------;■=- (2 ;

9 а*.а2 *

Г33>

'остальные равны нулюХ Введем обозначение [В]= г [Ь] [Я ]■

я

Матрица [В] имеет размерность 4x8, при этом

■ц- -32- 40^ЕТ [(ъ±-ьэ)+(ь3+ь2)п-(ь^ъ2)!;]:

Ь3- =

Ь

13- “34- 4ШТ 1(Ь*+Ьа)+(Ь*-Ь*>1>+<\+ЬгХ]:

ь^= ь*«= -4ШТ [ (Ьэ -Ь1 >+ < ‘Ьэ -ь2 »>+ ^ь2 -ь* * ]:

ь =

17

ьзв- да ^-ГЬ1+ЬЭ ;-ГЬэ + Ь2 )Т)+(Ъ1-Ь2 )К1; ь

'31- “22- 4а^Т[(а3-*±>-(*г+*3>*+(\+*гХ 3:

99- “24- 4й^^Га»+а9}+^а9-а2:>Т?-('а1+а2^:г;

ь = ь =

99 24

ь = ъ =

-^Г^Х [(а.1-аэ)+(а2-аз)т)+(а1-а2 К 1:

(34)

ьэ7 = ь28 = -4^ТТ 1Га1 +аэ ■) + Га2 +а3 * (а2 -а, X ^;

70

Ма'гемат. моделир. сис'гем и проц. N КО, 1992

Решение < ч »симмет-ричной Зс'.плчч Ь,4=-~ 04 )а~п);

■I

Ь С | ,1£ > { >’ ■}■« ) ;

а остальные элементы матрицы [В] равны нулю.

Интегрирование при вычислении матрицы [с] и [Н] целееооб{ зно провести по квадратурной Формуле Гаусса [I]:

11 эз

Я/С? .г)) дК&п = £ £ А.*. , (35)

-1-1

где

?2=Т>2= О ;

кг=8/9

После варьирования Функционала П отдельно по « и Г< мот/

я

чим систему линеиных алгебраических уравнений (СЛАУ) следующего вида:

Г [Н] [СП Г ь 1 Г о ]

I СС]Г о ] [ и ] [р]

(36)

Рассмотрим каждую компоненту СЛАУ в отдельности, [и] пред-' с тавляет собой блочно-диагональную матрицу размерностью (8я$3)* (8*113) Элементов, где |?8 ~ общее число конечных элементов. Она имеет вид:

[Е]

О

О

(37)

Матрица [с] может быть представлена в виде:

Т Т Т

[с] =[с :с,: 1 , 2 .

:с т]

. М!3

(38)

где [С ]Т получается из [с ]* следующим образом. Если в г.'

V

••бальной нумерации узлов конечно-элементной модели локалы-

номера узлов 1,2,3,4 Сдля конечного элемента с номером О равны соответственно к,і,т,п, то первая и вторая строки [о ]т пе-

І

Т

реносятся в строки 2У -1 и 2к матрицы [С ] , третья и четвертая

- в строки 2И--1 и 21 и т.д. Таким образом, в [С ]Т из 2 * яд

(где ЯД - общее число узлов в модели) лишь восемь отличны от нулевых; остальные строки - нулевые. Вектор Ь имеет вид:

т *т *т *т

b =[ftt ,Р2...............Рыз] . (39)

* * *

где p. ......Рві]. i.= t............NЗ

Вектор и выглядит так:

U=[U ,W ,V ,W ..................U .W f. (40)

ql ql q2‘ q2 ' qMR qNR 1 '

Вычисление вектора p индентично обычному методу конечных элементов в Форме метода перемещений [1]. Отметим, что СЛАУ (37) может быть очень плохо обусловлена из-за того, что элементы матриц [Н] и [С] имеют разный порядок. Поэтому перед решением системы С37) следует применить масштабирование.

Представленная версия ГТ1КЭ реализована в виде прикладной программы для ЭВМ IBM PC/AT (транслятор Турбо-Си 2.0). В качестве тестового примера рассмотрим задачу анализа НДС шарового подпятника (рис.2) (решение этой задачи, представленное в [2], является ошибочным). В сферической системе координат уравнения в перемещениях и условие несжимаемости выглядит так [21:

1 а ~ 1 д ,

—-------------(Bin Р )------------------[ain Р-----] +

г в in <р д<р г a in &<р д г

д2 и д2 (т- w) &S

= О ;

д т &<р а т*2 а$>

а (г2 и) 1 d(w sin <р)

(41)

г2 а г г ain <р дір

= О .

где и - переменная вдоль координаты г ; ш- переменная угловой координаты р ; в - гидростатическое давление.

*7777777

Рис. 2. Расчетная схема.

Кинематические граничные условия имеют вид:

и(Иг)=0 ; и(111)=-& совр ;

т(Я_)=0 ; и>(й' )- А згпр .

Аналитическое решение задачи таково:

Д д(г /(г))

■и,- Д/(г) ооар ; и>=----- ------------з<тг р

2г

в г

о>~=ОА(бС г+ЗС г 2 +6С г *) соа<р

г 2 Э 4

а =?а =СА (12С г-ЗС г ) ооар ;

С7 2 4

т~ =• ОД ҐЗС г+ЗС г ) вітир ;

Гр 2 4

С >ОС2г+Сдг ) совр .

ВДОЛЬ

(42)

(43)

где в ~ модуль сдвига материала;

Рис. 3. Конечно-элементная аппроксимация

С± = (5-да2 +4хх~3 )ТГ% ; =3(а2 -» №~2*2>~‘ :

С =6Н4 (а3 -сГ 2 )В~ 1 ; С =21? (аГ 3 ~) )В~± :

(44)

а=И2/Я±; 1)=10-9(<я~2 +аа )+4(ы3 -сх~а } ;

/(г)=С±+С2г*+Сзг %С4ГЭ .

В качестве примера рассмотрим задачу со следующими исходными данными: О* 1 МПа, V>-0.5, А —1 мм, Д - ЮО мм, Д =В0 мм. Ко-

1 2

нечно-элементная модель Срис.З) содержит 100 узлов и 76 элементов. Результаты таковы: максимальное вертикальное перемещение имеет узел 41', при этом аналитические Формулы дают иг*3.89 мм, а численное решение иг-3.92 мм. Напряжения вычислялись в центрах конечных элементов, при этом интенсивность напряжений напряжений достигает максимума в элементах- 58, 59, 60, 61

Стабл.2)-Как можно видеть, погрешность численного решения не превышает 8% .

Таблица 2

Элемент а аналитическое, МПа t а численное, МПа

58 0.988 0.999

59 0.981 1 .06

60 0.968 1 .02

61 0.948 1 .01

Литература

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

2.Лавендел Э.Э. Расчет резино-технических изделий. М:Маши-ностроение, 1976. 232 с.

3.Piati X.ЕГ-H..Wu Chun. A rational approach for choosing stress terms for hybrid finite element formulations. Int. J. for Sum. Ueth. in Eng. Vo I 26. 1988. P. 2331-2343.

Московский институт электронного машиностроения