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

Том 152, кн. 4

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Физико-математические пауки

2010

УДК 539.3

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ КОНЕЧНЫХ

ДЕФОРМАЦИЙ ГИПЕРУПРУГИХ ТЕЛ. IV. КОНЕЧНОЭЛЕМЕНТНАЯ РЕАЛИЗАЦИЯ. ПРИМЕРЫ РЕШЕНИЯ ЗАДАЧ

Ю.Г. Копоплев, Л. У. Султанов

Аннотация

Настоящая работа является продолжением цикла статей «Численное исследование конечных деформаций гиперупругих тел» и посвящена копечпоэлемептпой реализации алгоритма, который описан в первых трех частях. Приведены различные потенциалы упругих деформаций. На примере одного из mix проведено численное численное исследование конечных деформаций.

Ключевые слова: конечные деформации, гиперупругость, метод конечных элемеп-

А.И. Голованов

Введение

Первая часть цикла статей была посвящена общим вопросам нелинейной механики твердого деформируемого тела, рассмотрены основные положения кинематики конечных деформаций, приведены основные тензоры и соотношения, описывающие движение сплошной среды. Изложены основные виды вариационных уравнений, использующихся в численных методиках решения [1]. Вторая часть посвящена вопросу построения определяющих соотношений при конечных деформациях для гиперупругого изотропного материала, а также для слабосжимаемых материалов [2]. Третья часть посвящена описанию методики численной реализации исследования гиперупругих тел па основе метода последовательных нагружения. приводятся различные постановки и способы реализации алгоритмов решения [3].

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

1. Слабосжимаемый изотропный материал

Рассмотрим слабосжимаемый изотропный материал, кинематика среды будет определяться левым тензором деформаций Коши-Грина (В). Тогда упругий потенциал запишется в следующем виде:

W = W {11в,/2в) + (7).

В результате имеем физические соотношения для тензора истинных напряжений в виде суммы шарового тензора и девиатора [2]:

(Е) = фо (I) + 2 (<?') ,

где

3W,

Для описания сжимаемости часто используют величину истинной сжимаемости в виде

с J d¿To d¿To ^ ^

Если эту величину считать постоянной, из (1) следует закон сжимаемости при умеренных давлениях:

ст0 = K ln J, (2)

где K = 1/kc - модуль объемного расширения. Для малых давлений (малых изменениях объема) из (1) получаем упрощенный закон сжимаемости:

В качестве других возможных законов сжимаемости в [4] предлагается использовать:

закон Мурнагана (F. Murnaghari)

(J - l) + (kc - 1)-1 JJ1-kc - l) ; (3)

закон умеренных давлений

W0 = K [J ln J - (J - 1)] ;

упрощенный закон

K 2

W0 = -(J-iy. (4)

В работе [5] предлагается закон, обобщающий неогуковский материал и материал Муни Ривлина:

k=i k

В работе [6] приводится выражение

" J2 - 1

W0 = К--ln J

(6)

В работах [7 9] используется потенциал

Wo = |(lnJ)2,

который приводит к закону (2). В работе [11] дано выражение

А 2

W0 = - (ln jy -¿tln J.

Отмстим также потенциал [10]

2yU, 3

«-¿Шл— -О-

k=1

Сдвиговые напряжения описывает функционал Ш (/1§, ) > простейший вид которого [9]:

* = I (11ё ~ 3) • (7)

Отсюда следует

(Е') = м(в').

Дополнительные сведения о возможных формах потенциала упругой энергии изменения объема можно найти в [14]. Для справки приведем некоторые выражения потенциала Ш (/^, /2в) • В паре с выражением (5) используется

W = Ё Ck (I1B - 3)

k=1

с выражением ( 6)

где

k=1 k

r l. г 1 . г 11 ■ Г 19 ■ г 519

<-1 — С'2 — —, Сз — ———, С'4 — , С'5 —

2' 20 1050' 7000' 673750

2. Физические соотношения и разрешающее уравнение для неогуковского материала

Рассмотрим пример построения физических соотношений для следующего потенциала упругих деформаций (здесь используются выражения (4) и (7)):

(/1ё-3) + |(7-1)2, (8)

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

которое можно преобразовать к виду:

(¿) = (Л) • • (d) + (h) • (Е) + (Е) • (h)T - (Е) I1d, (9)

(Л) = (ЛЕ') + (ЛСТо),

4 / d2W' N

2J

-5/3

\hB (Си) ~ ¡ (В) (/) - i (/) (В) + h1B (/) (/)

2К

±(.7-1)(/) (1) + \.7(1) (I)-(.7-1) (СП)

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

(£Т^) = (кЛ • • (¿) .

(Ю)

Для решения задачи применяется метод последовательных нагружеиий [3]. В качестве основного используется вариационное уравнение принципа виртуальных мощностей, записанное в актуальном состоянии, которое для к-го состояния будет иметь вид [3]:

Ук

(к Ё) •• (¿к ¿) + [к У у • ки (к Е) •• (¿к ¿) -

_ I . . ^Н) • (к1г) + (к1г)Т • (6к1г)Т] - [кЧу • кг>] кГ • <5г>} с1Ук +

I {к • (кН)Т - [кУу • ки] к• ¿и ¿Бк = Iкг • ¿V ¿Ук + у Н; • ¿V ¿Бк

1

У (к Ё) • • (¿к¿) ¿Ук -1кг* • ¿и ¿Ук -1 к С •¿и ¿БЛ .

I Ук Ук «г

Подставляя в последнее уравнение физические соотношения (10), получим разрешающее уравнение

[ {(кй) • • (кЛ) • • (<М) + 1 (к£) • • [(51г)Т • (к/г) + (к1г)Т ■ (51г)

Ук

- [к У у • к и кг • ¿и} ¿Ук + {к < • Т - [к У у •к ик • ¿и ¿Бк

к Г • ¿и ¿Ук + к 1 п • ¿и ¿Бк -

Ук

" ~Ь. 1 / ' '^ сШк ~ I "Г ' сШк ~ I ЛБк \ (П)

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

к

уместно принять производную по времени как отношение приращения соответствующих величин, получаемые при переходе с к-го состояния в (к + 1)-е. Например,

к Дки

Аг

Теперь, в силу произвольности параметра £ вполне допустимо припять Д£ = 1. В результате из (11) получаем уравнения для приращений перемещений Дки.

3. Конечноэлементная дискретизация

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

8

к у (£0 = Ек у№ (£0,

'Г1 (12)

k и j = Еk un (£j),

г=1

где ^ {?) = 1 (1 + &1) (1 + (1 + функция формы, -1 <

< 1, ку\ - координаты узлов (к - номер шага нагружения, £ - номер узла в элементе), ££ = ±1 — координаты соответствующих узлов в локальной системе координат, кскорости узлов.

Приведенные выше функции формы обладают свойствами полноты и конформности, то есть образуют полные полиномы и при стыковке двух элементов обеспе-

дки дки

чивают непрерывность функций кг/ и кК'. Для перехода от -^т к . строим матрицу Якоби [А]:

З'1 Qij Z^ •"' Qej •

1;: ^rr ^и,

t=i

Затем вычисляем матрицу, обратную к [A] , и с ее помощью находим пско-

д д

мыс производные ifcCl = [fcAl . ——г = кС-ц —— . Так. например, вычисляется

д k yl o£j

производная от функции формы:

m_ m

('— Ht,i- (13)

к

Теперь с помощью (12) и (13) определим компоненты тензоров (Б), (К), (¿):

8

Б* =]Т М(14) г=1

■■ диг ■■ vlNt (£j) А . dNt (ej

h%3 = v""3 = ' V ' =VG„,ui *} дуо дуз ^ 3 1 dim

¿Cjm v\Ntm = £ utHtj, (15) t=i t=i

1 ° 1

(16)

t= 1

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

(И):

(d) • • (Л) • • (öd) = d^Ajiklödlk = ± [(v^ + v**) -] Ajm\ {5vl'k + 5vk'1)

(Лjiik + Ajjik + Kjiki + Ajjki) Su ' —

— vtHt,jT (Ajjifc + Ajjifc + Ajjfc; + Ajjfc;) Hgjöv

(E) • • (Sh)T • (h) + (h)T • (Sh)

aij un'jSunj — aij vnHtjHj Su?,

[kVy • kU kf* • Sv — umH/NqNM, tn • (h)T - [Vy • u] tn — (ujtftj- t;mumHt,m) Suj.

Таким образом. (11) после аппроксимации принимает вид

| {(d) •• (A) •• (Sd) +

V

+ i (fcE) • • [(Shf • (fc/>) + (khf • - [Vy • v] f* • <h>} dl/

E E E

V q,t,s=1 k,l=1 i,j = 1

t Ht ,j'7 (Ajük + A.yifc + Ajiki + Aijki) Hsjövs +

+ ajjutHt,iHa,jSuls -utHtjjNqNsSul} det [J] d^3d£3. (17)

/ {^П • - [Уу • и] • ¿и dS =

л л 3 8

= // Е Е [и?^ - ¿и^ ^] de1 de3. (18)

При вычислении интегралов используется схема численного интегрирования интегрирование заменяется суммированием значений подынтегральных выражений в квадратурных точках, умноженных на весовые коэффициенты. Тогда после интегрирования (17), 18) получим матрицу левых частей {¿и}Т [К] {и}. Используя соотношения аппроксимации (12). получаем вектор скорости узловых нагрузок, вектор невязки:

jf * • ¿и dV + ! {П • ¿и dS ={Р} {¿и} , ~~Кг ) / ^ " ^ с1У ~ / Г'6л> с1У ~ /^ [ = ^'

V )

1

— „,j,j

—u

Рис. 1. Растяжение полосы

Рис. 2. Деформированное состояние полосы

Окончательно получается система линейных алгебраических уравнений на к-м шаге нагружения:

[кК] {Дки} = {ДкР} + {кн}.

4. Численные примеры

Приведем результаты расчетов для упругого потенциала, заданного в виде (8).

Рассмотрена задача о плоской деформации квадратной полосы со сторонами 20 х 20 мм, / = 0.4225 Н/мм2, К = 5 Н/мм2 (рис. 1). По вертикальным краям полосы заданы перемещения. Так как полоса имеет две оси симметрии, то была рассмотрена четверть полосы с заданием соответствующих условий симметрии с разбиением на сетку конечных элементов 16 х 16.

Полоса в горизонтальном направлении была растянута в 3 раза, а в вертикальном в 0.3712 (для среднего сечения) (рис. 2). Полученные результаты хорошо согласуются с решением этой задачи другими авторами [13] (0.3711).

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

д = 0.1128 Н/мм2. Радиус плиты 50 мм, толщина 10 мм, / = 0.4225 Н/мм2 , К = = 5 2

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

На рис. 7 представлена зависимость перемещения центральной точки А, лежащей на верхней грани, на рис. 8 — зависимость нормальных напряжении ахх точки А и В центральной точки, лежащей на нижней грани. Из рис. 8 видно, что в верхней точке А возникают растягивающие усилия, а в точке В до определенного

Рис. 4. Деформирование плиты. Промежуточное состояние

Рис. 5. Деформирование плиты. Промежуточное состояние

Рис. 6. Деформированное состояние плиты

! 20-

-0.04 -0.08

нагрузка

-0.04 -0.08

нагрузка

Рис. 7. Перемещения точки А

Рис. 8. Нормальные напряжения ах

Рис. 9. Начальное состояние плиты

0.8

0.6 -

0.4 -

« 0.2

0 -

0

-0.2

-0.12

-0.12

Рис. 10. Промежуточное деформированное состояние плиты

момента наблюдаются сжимающие усилия, а затем и растягивающие, то есть плита сначала изгибается, а затем начинает постепенно раздуваться.

Решена задача об упругом деформировании квадратной плиты под действием равномерного давления д = 0.12625. Граничные условия и свойства материала задавались так же, как в предыдущем примере. Плита со сторонами а =100 мм

Рис. 11. Промежуточное деформированное состояние плиты

Рис. 12. Деформированное состояние плиты

I 20

0 -0.04 -0.08 -0.12

нагрузка

а 0.4 -

-1 -0.4

-0.16 0

-0.04 -0.08 -0.12 -0.16

нагрузка

Рис. 13. Перемещения точки А и точки С Рис. 14. Нормальные напряжения ах

точки А и В

1.6

1.2

0.8

0

-20

и толщиной Н =10 мм (рис. 9). На рис. 10, 11 изображены промежуточные состояния, на 12 деформированное состояние плиты. На рис. 13 представлена зависимость от нагрузки перемещений: вертикального перемещения центральной точки Л, лежащей па верхней грани плиты, и поперечного горизонтального перемещения точки С. На рис. 14 показана зависимость от нагрузки нормальных напряжений ахх точки Л и центральной точки В, лежащей на нижней грани плиты. Плита сначала изгибается, а затем постепенно начинает раздуваться, то есть края плиты сначала смещаются в стороны, а затем внутрь.

Заключение

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

Работа выполнена при финансовой поддержке РФФИ (проект Л*1' 08-01-00546а).

Summary

A.I. Golovanov

Yu. G. Konoplev, L. V. Sultanov. Numerical Investigation of Large Deformations of Hyperelast.ic Solids. IV. Finite Element Realization.

The article is the fourth part of a series of papers called "Numerical investigation of large deformations of hyperelast.ic solids" and deals with the finite element realization of an algorithm described in the first three parts. Some energy potentials of elastic deformations are presented, and for one of them the finite element realization is derived. The examples of problems and their solutions are given.

Key words: large deformations, hyperelast.icit.y, finite element method.

Литература

1. Голованов А.И., Коиопле.в Ю.Г., Султанов Л.У. Числешюе исследование конечных деформаций гиперупругих тел. I. Кинематика и вариационные уравнения // Учен, зап. Казап. уп-та. Сер. физ.-матем. пауки. 2008. Т. 150, кп. 1. С. 29 37.

2. Голованов А.И., Конопле.в Ю.Г., Султанов Л.У. Числешюе исследование конечных деформаций гиперупругих тел. II. Физические соотношения // Учен. зап. Казап. уп-та. Сер. физ.-матем. пауки. 2008. Т. 150, кп. 3. С. 122 132.

3. Голованов А.И., Конопле.в Ю.Г., Султанов Л.У. Численное исследование конечных деформаций гиперупругих тел. III. Постановки задачи и алгоритмы решения // Учен. зап. Казап. уп-та. Сер. физ.-матем. пауки. 2009. Т. 151, кп. 3. С. 86 98.

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

5. Yeuh О.Н. Some forms of strain energy function for rubber // Rubber Cliem. Teclin. 1994. V. 66. P. 754 771.

6. Arruda E.M., Buyce M.C. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials // J. Mecli. Physics Solids. 1993. V. 41. P. 389 412.

7. Auricchio F. A robust integration-algorithm for a finite-strain shape- memory-alloy superelastic model // Int. J. Plasticity. 2001. V. 17. P. 971 990.

8. Auricchio F., Taylor R.L. A return-map algorithm for general associative isotropic elast.o-plastic materials in large deformation regimes // Int. J. Plasticity. 1999. V. 15. P. 1359 1378.

9. Simo J.S., Taylor R.L., Pinter K.S. Variational and projection methods for the volume constraint in finite deformation elast.o-plastisit.y // Comput. Motli. Appl. Mecli. Eng. 1985. V. 51. P. 177 208.

10. Simo J.S., Pinter K.S. Remarks on rate constitutive equations for finite deformation problems: computational implications // Comput. Moth. Appl. Mecli. Eng. 1984. V. 46. P. 201 215.

11. Storakers B. On material representation and constitutive branching in finite compressible elasticity /,/ J. Mecli. Pliys. Solids. 1986. V. 34. P. 125 145.

12. Bonet J. Large strain viscoelastic constitutive models // Int.. J. Solids Struct. 2001. V. 38. P. 2953 2968.

13. Bonet J., Wood R.D. Nonlinear continuum mechanics for finite element, analysis. 1997. 283 p.

14. Hartmann S., Naff P. Polyconvexit.y of general polynomial-type hyperplastic strain energy functions for liear-incompressibility // Int.. J. Solids Struct. 2003. V. 40. P. 2767 2791.

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

Голованов Александр Иванович доктор физико-математических паук, профессор кафедры теоретической механики Казанского (Приволжского) федерального университета.

Коноплев Юрий Геннадьевич доктор физико-математических паук, профессор, заведующий кафедрой теоретической механики Казанского (Приволжского) федерального университета.

Султанов Ленар Усманович кандидат физико-математических паук, старший научный сотрудник НИИММ им. Н.Г. Чеботарева Казанского (Приволжского) федерального университета.

E-mail: Leñar. Sultanov Qksu.ru