Постановка задачи численного моделирования конечных деформаций МКЭ Текст научной статьи по специальности «Физика»

М. К. Сагдатуллин

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

Ключевые слова: метод конечных элементов, конечные деформации, метрический тензор.

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

Keywords: a method of finite elements, finite strains, metric tensor.

In paper the fundamentals of a technique of numerical research of finite strains of the isotropic hyperelastic bodies, oriented on application FEM are stated. The first section is devoted kinematics of finite strains in the Lagrangian frame, tensors and measures of the strain, defining are entered speeds, various tensors of stress are reduced. In the second section physical parities of the hyperelastic isotropic environment are formulated, using the thermodynamics equations. In the third section it is reduced resolving the equation in a current configuration and the parities defining speed of change of a stress tensor of Koshi-Euler as linear function from a tensor of a spatial gradient of speed are output.

Введение

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

современных вычислительных технологий на базе МКЭ. Варианты постановки задачи и возможных алгоритмов решения представлены в работах [1-4, 17]. Они основаны на фундаментальных результатах нелинейной механики упругих сред, изложенных в известных монографиях и учебных пособиях [6,7,9,11,13,14,15,16]. Как правило, используется Лагранжевое описание кинематики среды и

трехмерная постановка. В качестве Лагранжевых

координат в [2-4] используются декартовые

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

криволинейная система координат, так как каждая

материальная точка среды однозначно

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

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

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

нагружения, и строятся выражения для скорости тензора напряжений как линейно зависящие от тензора пространственного градиента скорости.

Кинематика конечных деформаций

Процесс деформирования будем рассматривать в некоторой инерциальной системе отсчета, в которую введем декартовую систему координат с ортами е1, е2, е3.

Пусть в исходной конфигурации радиус -вектор материальной точки имеет вид

я = Xі (, #2, ..,

где £,х, І;2, - криволинейные Лагранжевы

координаты.

Определим

- базисные векторы

Як = дек Єі = ЯкЄі,

Щ

- сопряженные базисные векторы

- дРк 1 - -

~>к -* ктп ^ т~) т-)к,і-

Rm x Rn = R ’■5..

m n

дХ‘ ' 2^0

- элементарный объем

40=Д •[ Д х Яз ],

- метрический тензор

(с ) = 0() = 0 (Щ ) = ОД^),

где

а. = я. • я. =х дх:дх:=Е КтКт,

,] ' 1 т д? д? т 1

О« = Я1 • Я1 = У-д?_ Л? = у я,’тЯ1 т, т дХт дХт т

О. = £0тпЯ‘тЯ1.

т,п

По аналогии в актуальном состоянии

введем

■ радиус - вектор материальной точки

г = У (?,??> ),

■ базисные векторы

_ дУ _ _,■_

гк =д7Те,- = гке,, д?

■ сопряженные базисные векторы

~к д? _ 1 ктп ^ ^ _к ,■_

Гк ^ ^ Х Гп = Гк 4,

■ элементарный объем

= Г •[ Г2 Х Г3 ] ,

■ метрический тензор

(Ё) = Ё,. Ё (г1?1) = X Ё. (),

(1)

где

Р^т р)у.т

* * Х^ ил ил Х"' т т

о,, = г. • Г. = 7------------------=7 Г Г. ,

^ 1 3 ^ д^ д^ 3

О Ї = Гі . Г3 = -^1 ^ = У Гі-тГ3т дхт дхт т ’

=7

утпГ'Г3

5 і] ^ & т п '

В дальнейшем нам понадобится вектор скорости и , который определим в следующем виде

0 = | Г = и ),

где ґ - время.

В задачах статики под вектором скорости будет пониматься приращение радиус - вектора Г по мере нарастания деформаций, то есть, фактически, и есть

вектор приращений перемещений АП, где П -вектор перемещений.

и = г - Я = П (, #2, #3).

Введем в рассмотрение группу тензоров, описывающих кинематику среды при ее деформации.

- Тензор градиента деформации

(р) = (Г) = ои (Якг) = Ёк‘ (я/к) = У р (ее),

и j

где

~ дХу

Р. = —— = отпЯ1г' = О Ят’’гп’3.

] дХ 3 п т тп

- Тензор градиента места

(р )Т = (Я) = о* (гкЯі) = о" (г-Як) =

=у р (ее Ь-Х (еіе.),

- Обратный тензор градиента деформации

(р-1 ) = (Я ) = Окі (гкЯ ) = ок (Як ) =

=У Р*(ее )=^ХТ(ее),

- Обратный тензор градиента места

(Р-1 )Т =(#г ) = Окі ((кГі ) = Окі (яігк )

(2)

Ур (іЄз) = дХ_(ее.)-

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

(р-1 )= 1 є є — дХ^(ее)

' ' 2аег|р.| тп ]Ы дХк дХ1 ' ' ’’.

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

Далее при описании деформации будем использовать

- Правый тензор Коши - Грина (мера деформации Коши - Грина)

(с) = (р) ‘(О) -(р) = 8]{Я'Я) = У Ё . (еЄ),

где

= ОтпЯт,і Яп

■ Левый тензор Пиолы (мера деформации Альманси)

(В-) = (р-1 )Т ‘(О)‘(р-1 ) = О(ГГ) =

=Ув* (ее). (3)

где

- Тензор деформации Коши - Грина

(Е )=1 [(с )-(О )] =1 [ О, - О ](* ) = У 4(^3),

где

Ё-. =-\о - О 1Ят ''Яп

у ^ тп тп J

Тензор деформации Альманси

^^) = 2[(О) -(в-1 )^ = -2[Оі3 -О3](Г 1'г-<^

= У 4 (ее),

(4)

і 3

где

А.. =1 [ - О \гт’‘гп’1.

1 2 тп тп J

Отметим, что компоненты тензоров деформаций Коши - Грина и Альманси в криволинейных базисах совпадают между собой. Справедливо соотношение, связывающее тензоры деформации Коши - Грина и Альманси

(Р) = (Р) •(А^ (5)

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

(А) = (Р-1 )Т •(ЕУ(Р-1),

преобразует метрический тензор в пространственный.

Теперь рассмотрим тензоры, используемые для описания течения среды. Базовыми тензорами здесь являются

- Тензор пространственного градиента скорости

RO^M^ ) = (ОГ ) = У Г; (ee ) (6)

где

до дО

доя dxm

0 = =--------- e = U^e = У---------:------Г

1 д? д? 1 ' 1 — д? д?

7~ i J до

h.. = о г =---------,

1 я dxJ

Тензор деформации скорости

(d T = j[(h T+(h )T ] = У dij ((ej),

(7)

(В)

где

дО дО Ixj Ixl

d.. =1 Ur-1 + oJr-J ]= 1 ■j 2 [ я я ] 2

Тензор скорости поворота

И = 1 [(h )-(h )T ] = У®;

где

Iol до1

IxJ Ixl

Введем в рассмотрение материальную производную (полную производную по времени) тензора деформации Коши - Грина

(E T=1«, (RiRJ )=У Ё> (eie;T

(9)

где

g; =огГ; +UJ • r

=У(

г- T

E = 2 g -„Rm,iR"’1.

Далее воспользуемся преобразованиями. Из (1) следует

_ IxJ ?k

следующими

Из (7) следует

(d ) = 2 У

e; =--------- Г

1 I?

доо Ixя доо Ixm

+ -

I? I? I? I?

2 У s..(r r).

(rlrj) =

(10)

., 1

Таким образом, скорость изменения тензора деформаций Коши - Грина (9) и тензора деформаций скорости (8) в криволинейных базисах имеют одинаковые значения компонент. Следовательно, справедливо

(а ) = (р-1 )Т-(р)-(р-1). (її)

Запишем (11) с учетом (5), (6) в виде

7-1\т . а аґ \-

= (р-1 )Т{(р)Т '(Л)‘(рИр)Т-(ЛУ(р) +

+ (р )Т.(Л).(р )].(р1 ) =

= (ЛЖТ -(^Ж^К). (12)

(d ) = (F -1 )T • [(F )T •^A)•(F )KF -1)

Выражение в правой части (12) называют производной Ли (Lie Rate), для которой, с учетом (10), справедливо выражение

M=(d)=). аз)

2 i, j

По аналогии с операцией

дифференцирования по времени определим вариации основных тензоров.

Итак, справедливо:

- из (9) следует

1 с--

где

и )=2 s. (RiR; t=xse, (, т,

2 i’J

sg; = SO, • rj +SUj • Г. =У(-Г- + So-r-),

я

SEy =1 Sg^R-^”^1.

- из (11) следует

Sd) = (F-1 )T •(SE?T(F-1) = 2S (rJ) = У15~d; (e)

2 i,j

где

Sdy = 2 У(rk +Sir— )r^r^1,

2 k

- из (6) имеем

(Sh) = ySokrk (rr) = У Sh; (ee,),

где

Shi; =yso;

krkrя’irn’ j я n

Физическая модель гиперупругого тела

Для построения определяющих

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

деформирования упругой изотропной среды в виде рц-(ст)--(а) = 0. (14)

Здесь р - плотность материала, ц - функция свободной энергии, (ст) - тензор напряжений

Коши, (а) - тензор деформации скорости, который

определяется различными формами (8), (10-13).

Построим общие структурные

соотношения, принимая в качестве базового выражения - тензор деформаций Альманси (4). Выбор базового тензора деформаций предполагает, что функция свободной энергии зависит от компонент этого тензора. Таким образом, будем считать заданной функцию

Ц = Ц(А). (15)

Согласно правилу дифференцирования скалярной функции по тензору [5,11] имеем

*=(f К4),

где появляется тензор второго ранга

(16)

I*

и

Воспользуемся соотношением (12), из которого следует, что

k

m

(Л) = (а)-(л)-(й)-(*)Т (Л) =

- [(а)-(л)( а)-(а).(л)] - [(л) »-(»)■ (л)] .(17)

Подставим (17) в (16) и далее в уравнение (14) и

получим

Р

ду

л([(а )-(л)- (а )-(а ) (л)]-

-[(л).(»)-(»).(л)]}-(^)..(а )=о.

Преобразуем уравнение (18) к форме

-И}..(а)-

(18)

р

ГНИ ).(л)-(л).(1

-'Ш И-М-ЙУ)}» = 0. (19)

Так как тензор деформации скорости (а) и тензор скорости вращения (о) являются совершенно

независимыми тензорами, то из (19) следуют два уравнения

И = р

Зг НЇ )<лнлН£

(20)

Ц ]'(А) = (А)-(! ). (21)

Соотношения типа (20) называются физическими (определяющими) соотношениями или уравнением состояния и, по сути, определяют тензор напряжения в виде функции от тензора деформации. Уравнение (21) является ограничением на выбор функции свободной энергии (15).

С учетом (21) определяющие соотношения (20) можно записать либо в виде

И-Й )[(* )-2 (л)] = р(§ )(в-1Т,

либо

И = Р )-2 )-(**)]% р

Оба выражения эквивалентны, вследствие

справедливости (21).

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

Ц Ц (11А ,12 А , Л А ) .

Здесь, и далее, обозначение скаляров упростим, а именно

11 = 11А , 12 = 12 А , А = А А .

В результате, после несложных преобразований, получим

(<г) = [ц1 + Ц Л ](ё )-Ц (А) + Ц Л (А-1), (22)

где

дц дц дц

Ц = Р^Г, ^2 = р^-, Ц = р^-.

д/1 д/2 д/3

ду)

Для получения представления тензора напряжений (и) в виде суммы диад по соответствующему

базису, необходимо определить, в каком базисе определен тензор деформаций Альманси. Наиболее простой вид получаем при использовании базиса еі . В этом случае

И) = Уи(,ёiёj)

и

= [у + Г2Л ] ^ - Г2Л + УзЛл*,

где |Аи| образуется как матрица, обратная .

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

{<?) = &! (Гги) = ст (г/.)

= [ Г + Г211 ] О3 - Г2л + Гз734*3,

где

л = О3 - О

л* - компоненты обратного тензора в базисе (т/),

которые имеют достаточно сложный вид [6,11].

В практической реализации, вычисление компонент л* значительно проще, чем л* . Поэтому

целесообразно вычислять л* по л*. Для этого из тождества

л (ҐҐ) = У л (еіе3) = (л-1)

i.3

получим

4 = Г .у £п (ЄтЄп ) . Г = ^тп^^ .

т,п

Контравариантные компоненты тензора напряжений определяются в виде

і ^_ті пі

И = и О О .

тп

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

Разрешающее уравнение на шаге нагружения

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

|(ст)..(^а)аК = |/*.дисіУ + | ґ*8ис18, (23)

V V

где /* - вектор заданных внешних объемных сил, —**

ґп - вектор заданных напряжений на части

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

последовательность равновесных состояний Vk,

и

и

которые реализуются при заданных значениях

к у'»* к“** /"V ч_/

внешних сил / , /п . Определим в качестве основной

ч_/ к ^

неизвестной величины вектор скорости и, который можно трактовать как вектор приращения текущей конфигурации при переходе к состоянию V+1,

(24) где

ко = Дки = к+1г - кг,

кг = кУ (?,?2,?).

Разрешающее уравнение на текущем шаге строится путем линеаризации исходного уравнения

дки‘

(23) в предположении ^ и П 1. Детали построения

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

|{(к„ ^(л)-2 (V)-(Зк)•(кк) + (кк)Т • (Зк)Т

д кО

~дкУ‘

•[(kст)••((5kd)- к/*Зо]1 йУк

+1К*•(* )т -

д кО

• % \3odS, =

д кУ

= | 7 * •8ойУк + | % •8оёБк -

- {| [(кст) • • (Зка) -к?*0 ^ - [ % • О|.

Чтобы полностью определить это уравнение необходимо построить выражение скорости напряжений (кст) для известной конфигурации кг

через неизвестный вектор скорости (24) в виде линейной функции.

Рассмотрим общий случай изотропного материала с определяющим уравнением (22). Индекс « к » для сокращения записи опустим. Справедливо

(СТ) = [Ц1 + Ц2/1 + Ц/1 ] (Ё) + Ц1 + ЦА ] (Ё) -

-ц>2 (А)- ц (А) + Ц13 (А-1) +

+цД (А-1 ) + ц/3 (.А-1). (25)

Далее распишем каждое из слагаемых

д2ц ( д/1'

= Р

д/,. д/ ^ дА д/,. д/2 ^ дА ^

^ № \(,а ) =

д/, д/3 ^ дА ^ 1 ’

Нрад (ё )+р__а [ '■ (ё Н'4)]+

д2ц

&_Г3

+РЩ. /3 ).

Здесь использовались следующие представления

А = (£ ^ ) = (ё Иа),

^ =(д& )"(А )=[ '■ (* ИА)]"(А ),

;3 ^-(-'1 )=/ 3 (а-‘ ).

Справедливо тождество

Ж4-1 )=(ё). (26)

Если продифференцировать по времени тождество (26), то получим

(А-‘ ) = ( а-)•[(* )-(-4 Х)].

Итак, удалось построить выражения всех слагаемых в соотношении (25) через два тензора скоростей (А) и ( ё ). Причем все тензоры, которые свертываются с ними, однозначно определяются текущей конфигурацией кг и могут быть вычислены по соответствующим формулам. Теперь

необходимо выразить тензоры (А), (ё ) через

пространственный тензор градиента деформаций (6) в виде линейной зависимости.

Из (12) имеем

(кА) = (кё)-(кк) •(кА) + (кА)-(кН), (27)

то есть с учетом соотношения (8) (кА) есть линейная функция от (кН).

Из соотношения (4) получаем

(Ё ) = 2 (А ) + (В-1). (28)

Из (3) следует

(В-1 ) = (Р-11 •(0)•(Р-1 ) + (Р-1) <0)•(Р-1). Дифференцируя по времени тождество

(Р1 )•(Р )=(с),

получим

(Р " ) = -(Р-‘ )•(Р )'(Р ■' ) = -(Р-‘ )•( к)

и подставим его

(В1 )=-(к )т •(Р1 )т ^0 )•(Р-1 )-

-(Р 1 )т ^0 )•(Р 1 )•(h )=

= -(к)т •(В-1 )-(В-1 \(к). (29)

Собирая вместе представления (26), (27) и (28),

получаем выражение

(Ё) = 2[(а)-(к)т •(А)-(Ау(к)]-

-(к)Т •(В-1)-(В-1 )•(к) =

= 2(а)-(к)Т •[2(А) + (В 1 )]-[2^А) + (В-1 )^^(к) = = 2 (а )-(к )т •(ё )-(ё )^(к),

то есть соотношение полностью аналогичное (27).

Таким образом, представление для скорости изменения тензора напряжений (25) в глобальном базисе, то есть в виде

(кст )=Х кст,1(ее)

допускает представление

к к^тп к 1

ст.. = О к .

1 1 тп

(30)

Выражение для к0’тп из (30) легко строится с

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

Заключение

Полученные соотношения представляют собой теоретическую основу конечно - элементного алгоритма исследования конечных деформаций нелинейно-упругих тел при силовом их нагружении. Необходимо лишь добавить конкретную физическую модель в виде выражения функционала свободной энергии, справедливого для соответствующего материала. Некоторые рекомендации по их выбору даны в работах [8,10-14].

Литература

1. Голованов А.И. Метод конечных элементов в механике деформируемых твердых тел / А.И. Голованов, Д.В.

Бережной. - Казань: «ДАС», 2001. - 300с.

2. Голованов А.И. Численное исследование больших

деформаций гиперупругих материалов. Часть 1.

Кинематика и вариационные уравнения / А.И. Голованов, Ю.Г. Коноплев, Л.У. Султанов // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. — 2008. — Т.150.- №1. — С.25-37.

3. Голованов А.И. Численное исследование больших

деформаций гиперупругих материалов. Часть 2.

Физические соотношения / А.И. Голованов, Ю.Г. Коноплев, Л.У. Султанов // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. — 2008. — Т.150.- №1. — С.25-37.

4. Голованов А.И. Численное исследование больших

деформаций гиперупругих материалов. Часть 3.

Постановка задачи и алгоритмы решения / А.И. Голованов, Ю.Г. Коноплев, Л.У. Султанов // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. — 2009. — Т.151.- №1.

5. Голованов А.И. Математические модели вычислительной нелинейной механики деформируемых сред / А.И. Голованов, Л.У. Султанов. - Казань: Казан. гос. ун-т., 2009.

- 465 с.

6. Грин А. Большие упругие деформации и нелинейная механика сплошной среды / А. Грин, Д. Адкинс - М.: Мир. - 1965. - 455 с.

7. Елисеев В.В. Механика упругих тел / В.В. Елисеев. - С.-Петербург, СПбГТУ, 1999. - 341 с.

8. Корнеев С. А. Термодинамически согласованные уравнения состояния нелинейной теории упругости / С .А. Корнеев // Изв. РАН. МТТ. - 2003, №2, С. 71-82.

9. Перелыгин, О.А. Исследование прочности цилиндрических оболочек при наличии увода или смещения кромок сварных швов / О.А. Перелыгин, Н.М. Туйкин, Д.В. Бережной, М.Н. Серазутдинов // Вестник КГТУ, Казань, изд-во КГТУ. - 2001. - №1-2. - С.77-79.

10. Коробейников С.Н. Нелинейное деформирование твердых тел / С.Н. Коробейников // Новосибирск СО РАН, 2000. - 262 с.

11. Кузнецова В.Г. Эффект учета слабой сжимаемости материала в упругих задачах с конечными деформациями / В.Г. Кузнецова, А.А. Роговой // Изв. РАН. МТТ. - 1999, №4, С. 64-76.

12. Валиуллин, А.Х. Большие деформации и перемещения композитной цилиндрической оболочки / А.Х. Валиуллин // Вестник КГТУ, Казань, изд-во КГТУ. -2011. - №9. - С.109-117.

13. Лурье А.И. Нелинейная теория упругости / А.И. Лурье.

- М.: Наука, 1980. - 512 с.

14. Мальков В.М. Нелинейный закон упругости для тензора условных напряжений и градиента деформации / В.М. Мальков // Изв. РАН. МТТ. - 1998, №1, С. 91-98.

15. Оден Д. Конечные элементы в нелинейной механике сплошных сред / Д. Оден. - М.: Мир, 1976. - 464 с.

16. Черных К.Ф. Нелинейная теория упругости в машиностроительных расчетах / К.Ф. Черных. - Л.: Машиностроение, 1986. - 336 с.

17. Голованов А.И. «Нелинейная задача о гиперупругом деформировании полилинейного конечного элемента оболочки средней толщины» / А.И. Голованов, М.К. Сагдатуллин // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. -2010. - №4 (16). - С. 39-49.

© М. К. Сагдатуллин - канд. физ.-мат. наук, доц. каф. теоретической механики и сопротивления материалов КНИТУ, ssmarat@mail.ru.