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

МЕХАНИКА

J

УДК 539.3

Ю. И. Димитриенко, С. М. Царев, А. В. Веретенников

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

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

Резино-технические изделия (РТИ) — амортизаторы, уплотнения, клеевые соединения и др. — находят большое применение в различных отраслях современной техники. Однако надежных численных методов расчета эксплуатационных характеристик РТИ, учитывающих реальную сложную геометрию изделий, несжимаемость резин, большие деформации и нелинейное поведение РТИ, практически нет, в то время как теоретические исследования в области нелинейной теории упругости ведутся уже достаточно давно [1-5]. Эти эффекты в должной мере не учитываются и широко известными конечно-элементными программными комплексами, такими как А^УБ, МАБТЕАМ и другими. Цель настоящей работы — разработка конечно-элементного метода расчета напряженно-деформированного состояния осесимметричных элементов РТИ из эластомерных несжимаемых материалов на основе предложенных ранее [6-9] моделей материалов с конечными упругими деформациями.

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

Примем следующие допущения: 1) рассматриваемые материалы — идеально-упругие, обладающие конечными упругими деформациями (необратимые деформации вязкоупругости, пластичности и т.п. отсутствуют); 2) материалы несжимаемые и изотропные; 4) определяющие соотношения для материалов соответствуют полулинейной модели Лу, предложенной в работе [4].

Тогда математическая постановка задачи о расчете напряженно-деформированного состояния (НДС) элементов конструкций из рас-

сматриваемых материалов, формулируемая в терминах отсчетной (ла-

о

гранжевой) конфигурации К в тензорной форме имеет следующий вид [4]:

V •Р = 0, Xг е V; (1)

det F = 1, Хг G V; (2)

V

T = -p G-i +h I1(C)E + /2C; (3)

< P = T ; (4)

G-1 = F-1 • F , Xг G V U £; (5)

1 ° ° ° ° ° ° C = -(V®u + (V®u)T + V®u • (V®u)T), Xг G V U £; (6) 2

F = E + (V ®u)T, Xг G V U £; (7)

°

P • n = s, s = -pe n •F-1, Xг G £p; (8)

u • n = ue, P • ba = 0, Xг G £u, (9)

где (1) — уравнение равновесия; (2) — уравнение несжимаемости; (3) — уравнения состояния эластомерного материала; (4) — соотношения между первым тензором напряжений Пиолы-Кирхгофа P и пятым

энергетическим тензором напряжений T [4]; (5) — соотношение между правой обратной мерой деформации Коши-Грина G-1 и градиентом деформаций F; (6) и (7) — кинематические соотношения между правым тензором деформаций Коши-Грина C и градиентом деформаций

°

с одной стороны и градиентом вектора перемещений V ®u с другой стороны.

Уравнения (8) и (9) представляют собой граничные условия на по°

верхности области £. Рассматривается случай, когда на части границы

Ер в актуальной конфигурации задано давление pe, а на части границы

°

Е заданы смешанные граничные условия: нормальное перемещение ue

и

и условие свободного скольжения; n — вектор нормали к поверхности

°

£р, а ba, а = 1, 2 — векторы в касательной плоскости, ортогональные

°

к n, s — вектор поверхностных усилий. Предполагаем, что перехода материальных точек с одной части границы на другую в актуальной конфигурации не происходит.

В системе (1)-(9) также обозначено: p — гидростатическое давле-

° °. д

ние; р = const — плотность материала; V = r' qj^, — набла-оператор

в отсчетной конфигурации; • и ® — знаки скалярного и тензорного

о • °

произведений [10]; гг — локальные векторы взаимного базиса в К; Хг

оо

— лагранжевы координаты; V — область в К, занятая рассматривае-

оо

мым элементом конструкции; £ — граница области V; Е — метрический тензор; и — вектор перемещений; 11,12 — константы материала; /1(С) = С"Е — первый инвариант тензора [10]. Методика определения материальных констант /1,/2, а также примеры ее реализации рассмотрены в [5, 6]

Подставляя соотношения (6)-(7) в определяющие соотношения (3)-(5), а те в свою очередь — в уравнения (1), (2) и граничные условия (8), (9), получаем постановку задачи нелинейной теории упругости с большими деформациями, состоящую из 4 скалярных уравнений относительно 4 неизвестных функций: и и р.

Вариационная формулировка. Вариационная формулировка задачи (1)-(9) состоит из двух вариационных уравнений и имеет следующий вид:

V

T • • V dV = s • 8u d £,

(10)

V

(det G - 1 )SpdV = 0.

(11)

V

Здесь и — вектор кинематически допустимых перемещений, т.е. удовлетворяющих граничным условиям (9) для перемещений, а 0и — вектор возможных перемещений, удовлетворяющий нулевым граничным

° • о

условиям (9): о и • п = 0, Xг € £и.

Запись в физических компонентах. Запишем уравнения (1)-(9) в физическом базисе еа = га/На, где На — параметры Ламе [10]. Компоненты вектора перемещений и, градиента вектора перемещений

о

V ®и, тензора деформаций С, меры деформаций С и обратной меры деформаций С-1 в базисе еа имеют вид:

3

3

3

с = ^ £aß еа ® ее, G = ^ daß еа ® eß, G 1 = 9°'в ea ® ee,

a,ß= 1 3

a,ß=1

a,ß= 1

3

u = ua.ea, V ®u = uaßea ® e^,

a,ß=1 a,ß= 1

(12)

где обозначены ua — физические компоненты вектора перемещений u:

_ - 1 ÖUß Ua дИа

uaß — uaß + "aß u a, uaß —

u a —

Иа dXa HaHß дХ

1 3 ли (13)

u^ dHa

Ha И дХY '

7=1 '

а также физические компоненты тензора деформаций и метрической и обратной метрической матриц:

1 3

£ав = 2 (иав + ива + иа7 ^7), (14)

а,в,7=1

gae — ¿aß + 2eaß, (15)

¿aß (1 + 2eii) — 2eaß

gaß — ¿ae(1 + 2eM) - 2£aß + 2 eaij Eßkl eikeß. (16)

Здесь Еаг^ — символы Леви-Чивиты (по повторяющимся латинским индексам идет суммирование).

V

Компоненты Тав тензора Т в базисе еа совпадают с компонентами истинного тензора напряжений Т в физическом базисе актуальной конфигурации [4] и на основании уравнений (3) имеют вид:

V 3

Т = ^ Тав еа ® ее, (17)

а, в=1

Тав = ®ав - рдав, (18)

= Ь^гг^ав + 212£ав. (19)

Подставляя выражения (13)-(19) в соотношения (10), (11), получаем запись вариационных уравнений в физических компонентах:

3 р 3

У, / Тав$£ав<IV = ^ / 8а5иа<1 £; (20)

а,в=1 о а=1 {

V

(det(5ae + 2eae) - 1) ¿pdV — 0. (21)

V

Метод МКЭ для решения нелинейной осесимметричной задачи. Рассмотрим далее случай осесимметричных конструкций, тогда в качестве Лагранжевых координат Хг удобно выбрать цилиндрические координаты: X1 = г, X2 = г, X3 = в, причем Н1 = 1, Н2 = 1, Н3 = г. Полагаем, что имеет место осесимметричное нагружение, т.е.

Рис. 1. Шестиузловой КЭ с 15 степенями свободы

и3е = 0, з3 = 0 на В этом случае окружные перемещения отсутствуют (и3 = щ = 0), а радиальное и осевое перемещения щ = иг, и2 = щ зависят только от г, г. Компоненты тензоров деформаций , егв) и напряжений (Тг^, Тг$) также в осесимметричном случае равны нулю, а остальные их компоненты зависят только от г, г.

Для решения вариационных уравнений (20), (21) в осесимметричном случае, был применен метод конечных элементов со специально разработанным новым шестиузловым треугольным конечным элементом. Этот элемент имеет 15 степеней свободы: по два перемещения в каждом узле и по одному значению гидростатического давления в каждой вершине (рис. 1).

Аппроксимация в каждом конечном элементе (КЭ) по перемещениям была квадратичной, а по гидростатическому давлению линейной:

{и} = [Ф] {q}, p = (ФР|т {y},

(22)

2 2x12 12 1 3 3

где {и} = (и1, и2) = (щ, щ) — координатный столбец перемещений в КЭ; {д} — координатный столбец перемещений в узлах; {у} — координатный столбец гидростатических давлений в вершинах КЭ. Матрица [Ф] и столбец {Фр} имеют следующий вид:

[Ф] = [Фав ] =

2x12

Ф1 0 0 Ф1

Ф2 0

0

Ф2

Фз 0

0

Фз

Ф4 0

0

Ф4

Ф5 0

0

Ф5

Фб 0

0

Фб

(23)

{Фр}т = ( Li L2 L3 ),

где Ф1 = Li(2Li -1); Ф1 = L2(2L2-1); Фз = Ьз(2Ьз-1); Ф4 = 4LiL2; Ф5 = 4L2L3; Фб = 4L1L3;

и = + Ь(г)Г + ;

а(1) = г(2)г(3) - г(3)^(2); Ь(1) = г(2) - г(3); (24)

с(1) = г(3) - г(2); 2^е = Ь(1)С(2) - Ь(2)с(1),

а Т({), г а) — координаты вершин треугольника. Остальные коэффициенты а а), Ь(ь), О(ь) получаются с помощью круговой перестановки индексов; Бе — площадь КЭ.

Координатный столбец деформаций |е}т = (£11, £22, £33, £12) = = (£гг ,£гг ,£$0, £гг) представим в виде суммы линейной и нелинейной частей:

{£} = {£} + {£} = ([В ] + [В ] + [ВШ, (25)

где матрица [В] имеет вид: [В ] = [{В(1)},..., {В(х)},..., {В(12)}],

4x12

{В (х)}т =(и*,и*,иХ3, (^ + Ь%1 )/2),

4

их =±_ фах дНа +

ав На дХа НаЩ дХв

3

ая

+наНт,дХ, х=1...12,»,/з=1,2,(2б)

7=1 '

а элементы матриц [В] и [В] вычисляются по следующей формуле:

12 12 1 12

В шх = £ НХЧ; Вшх = £ нрхдр; НХР =2 £ ¿а7^7, (27)

р=1 р=1 7=1

где х, р =1... 12, ш = 1... 4 (значениям индекса ш = 1... 4 соответ-свуют сочетания пар индексов (а, в) = (1, 1)(2, 2)(3, 3)(1, 2)).

Вариации компонент вектора перемещений и тензора деформаций, а также гидростатического давления р с помощью соотношений (22) и (25) представляем в виде

¿{и} = [Ф] 5{д}, 8р = {Фр}т ¿{у}, ^{£} = ([В] + [В] + [В]Щд}. (28)

Вводим координатный столбец {а}т = (а11, а22, а33, а12) = = (агг,а22,о00,оГ2), который на основании определяющего соотношения (19) связан с координатным столбцом деформаций линейным соотношением

М = [А]{£}, (29)

4 4x4 4

где [А] — матрица модулей упругости изотропного материала [10]:

[A] =

44

(30)

/ /1 + 2/2 /1 /1 0 \

/1 /1 + 2/2 /1 0 /1 /1 /1 + 2/2 0 0 0 0 2/2 Вводим также столбец компонент обратной метрической матрицы {д}т = (у11, у22, у33, у12) = (#ь у2, у3, #4), тогда соотношения (16) можно записать в матричной форме:

где

{/о}т = (1,1,1, 0);

{g} = {/о} + 2 [I] {f} ,

[I ] =[Ео ] + [R]; [Ео] =

4 4 4 4 4 4 4 4

(31)

0110 1010 110 0 0000

(32)

Элементы Raß матрицы [R] зависят от деформаций и вычисляются следующим образом:

4

Raß = Ro

V F

aß FY,

7=1

1122зз

r23 = R31 = R13 = R31 = R12 = R21 = 1,

44 r34 = r43 = 2-

(33)

Все остальные элементы матриц Л^в, кроме указанных в (33), равны нулю.

С помощью выражений (29) и (31) соотношение (18) можно также записать в матричной форме:

{Т} = {а}-р {у} = [А]([В] + [ВВ]){д}- р({1о} + 2 [I] ([В] + [В]){д}).

(34)

Здесь {Т}т = (Т11, Т22, Т33, Т12) — координатный столбец напряжений.

Раскрывая детерминант в выражении (21), представим его следующим образом:

ае^ав + 2е«в) - 1 = 2 {1о}т {е} + 2к;

(35)

к = 2(е11е22 + е33е11 + е33е22 — е12) + 4е33(е11е22 — е12).

Введем координатный столбец внешних усилий {5} = (51,52) = = (зг, Зг), где вг, Зг — компоненты вектора усилий б (см. (8)) на поверхности тела.

Подставляя теперь (28), (34) и (35) в формулы (20), (21), получаем разрешающую систему для КЭ в матричной форме:

[K(q)] {q}-{Np(q)}{y} = {P};

12x12 12 12x3 3 12

[Bp(q)]T {q} = {hp(q)},

3 12 12 3

(36)

где обозначено:

[K(q)] = / ([B]T + [B]T + [B]T)[A]([B] + [B])d V;

12x12 ◦

V

[Bp(q)]= /" ([B]T + [B]Tpp]dV;

(37)

[Np(q)] = I ([B]T + [B]T + [B]T)(^p] + [I]([B] + [B])[#q]{y})dV;

12x3 J

О

V

12x3

О

V

{Р} = / [Ф№ = / •

о V

Е

Здесь также обозначены матрицы, полученные тензорным умножением столца на строку:

[Фр] = {/о}{ФР}т; ф]= {q}{Фp}т • (38)

4x3 4 3 15x3 15 3

Численный алгоритм. На основе локальной системы алгебраических уравнений (37), записанной для одного КЭ, далее была построе-

о

на глобальная система для всей области V. Обе эти системы являются

нелинейными, поскольку матрица [В] (см. (27)), а следовательно, и матрицы [К], [Вр], [^р] зависят от вектора неизвестных

Для решения системы нелинейных уравнений применяли итерационный метод, при котором на каждом текущем шаге итерации значения компонент вектора неизвестных в матрицах [К], [Вр], [^р] и векторах {Р} и {Л,р} выбирали с предыдущего шага итерации. Решение линеаризованных систем было проведено методом Гаусса.

Генерацию конечно-элементной (КЭ) треугольной сетки осуществляли "вручную", с помощью преобразования прямоугольника в

о

двумерную область решения задачи V в координатах г, г, при этом

о

полагали, что V ограничена двумя торцевыми плоскостями г = 0 и г = гтах (см. далее рис. 4, а). В результате была получена симметричная относительно плоскости г = ^гтах КЭ сетка, которая, как показали

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

£ = -¿тах и устойчивости итерационного процесса по сравнению с типовой нерегуляризиро-ванной КЭ сеткой, характерной для стандартных коммерческих КЭ-пакетов.

Результаты расчетов. Тестирование разработанного метода проводилось путем сравнения численных результатов с аналитическим решением нелинейной задачи Ламе [1] о трубе, нагруженной внутренним давлением, торцы которой свободно скользят по плоскостям £ = 0 и £ = £тах, внешняя поверхность цилиндра свободная. На рис. 2 показана КЭ сетка для цилиндра в отсчетной и актуальной конфигурациях, а на рис. 3 — распределение радиальных и окружных напряжений в зависимости от радиальной координаты г. Получено хорошее совпадение аналитических и численных результатов, относительная ошибка не превышает 0,1 %.

Далее разработанный метод применен для расчета напряжений в осесимметричном амортизаторе, работающем при осевом сжатии. Граничные условия заданы в следующем виде: на поверхности £ = гтах

Рис. 2. Конечно-элементная сетка для цилиндра в отсчетной и актуальной конфигурациях под действием внутреннего давления (задача Ламе)

Рис. 3. Распределение радиальных и окружных напряжений в задаче Ламе, полученных конечно-элементным расчетом и по аналитическому решению

задано осевое смещение и2е = £тахгтах, где £тах — средняя деформация сжатия амортизатора, на поверхности г = 0 задано свободное скольжение, поверхность г = / (г) свободна от нагрузок, а на оси симметрии заданы условия симметрии

г = 0 : Тгг = 0, иг = 0; ¿тах = 0 : Тгг = 0, иг = и2е;

(39)

г = 0 : Тгг = 0, иг = 0; г = /(г) : Тгг = Тгг = 0.

На рис.4,а показана конечно-элементная сетка для амортизатора в отсчетной и актуальной (деформированной) конфигурациях. На рис. 4, б, в и 5-7 приведены результаты расчетов напряжений в амортизаторе при следующих значениях параметров: /1 = 1 МПа, /2 = 10 МПа, высота гтах = 20 мм, максимальный и минимальный радиусы гтах = 40 мм, гтш = 10 мм. Расчет проведен для значений £тах, равных -0,1, -0,5 —5 и -10 %.

Рис. 4. Конечно-элементная сетка для амортизатора в отсчетной и актуальной конфигурациях при осевом сжатии (а), распределение напряжений по осевой координате г на свободной поверхности амортизатора (б) и по радиальной координате г в зоне утонения амортизатора (в) при деформации осевого сжатия 2,5 %

Рис.5. Распределение касательных напряжений в амортизаторе (|Тгг|, МПа) при деформациях осевого сжатия 0,5 % (а), 5 % (б) и 10 % (в)

На рис. 5-7 показаны распределения напряжений Тгг, Тгг и Т$в в амортизаторе для трех значений деформации етах, а на рис. 4, б, в приведены графики зависимости напряжений от координаты £ на свободной поверхности амортизатора и от координаты г в зоне утонения амортизатора.

На основании проведенных расчетов сделано заключение о наличии следующих особенностей распределения напряжений в амортизаторе: значения осевых напряжений Тгг достигают максимума в зоне

утонения амортизатора (г = - £тах), а в его широких частях (г = 0 и

£ = £тах) они снижаются почти до нуля. По всей зоне утонения при 0 < г < гт1п осевые напряжения распределены практически равномерно (см. рис. 6).

Рис.6. Распределение осевых напряжений в амортизаторе (|Тгг|, МПа) при деформациях осевого сжатия 0,5 % (а), 5 % (б) и 10 % (в)

Касательные напряжения Тгг и окружные Т$$ имеют локальный максимум в зоне утонения, в окрестности свободной поверхности г = гтщ (см. рис. 5 и 7).

На распределения напряжений, полученные при значении деформации етах = —0,1 % (рис. 5, а, 6, а и 7, а), нелинейность задачи практически не влияет, поэтому это решение фактически совпадает с решением линейной задачи теории упругости несжимаемой среды с малыми деформациями. При значениях деформации етах, равных —5 и — 10 % влияние конечности деформаций оказывается уже существенным (рис. 5, в, 6, в и 7, в); линейный закон увеличения значений напряжений с ростом значений деформации етах уже не соблюдается — в нелинейной области в целом напряжения имеют меньшие значения, чем получаемые по линейной модели. Кроме того, меняется характер

Рис.7. Распределение окружных напряжений в амортизаторе (Т$$, МПа) при деформациях осевого сжатия 0,5 % (а), 5 % (б) и 10 % (в)

распределения окружного напряжения Т$$ (см. рис. 7) — при конечных деформациях появляется еще одна зона максимальных значений напряжений на оси симметрии тела.

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

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

Работа выполнена при поддержке грантов РФФИ№ 06-08-01448a и 07-08-00574-a.

СПИСОК ЛИТЕРАТУРЫ

1. Л у р ь е А. И. Нелинейная теория упругости. - М.: Наука, 1980. - 572 с.

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

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

4. Dimitrienko Yu. I. Tensor analysis and nonlinear tensor functions. - Kluwer Academic Publishers. - Dordrecht, 2004.- 675 p.

5. Dimitrienko Yu. I. Novel viscoelastic models for elastomers under finite stains // European Journal of Mechanics. A/Solids. - 2001. - N 1. - P. 38-56.

6. Димитриенко Ю. И., Даштиев И. З. Модели вязкоупругого поведения эластомеров при конечных деформациях // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2000. - № 1. - С. 21-41.

7. Димитриенко Ю. И., Даштиев И. З. Прогнозирование ползучести и релаксации эластомеров при больших деформациях // Вопросы оборонной техники. Сер. Композиционные и неметаллические материалы. - 2001. - № 1. -С. 51-55.

8. Димитриенко Ю. И., Царев С. М. Разработка конечно-элементного метода расчета эластомерных несжимаемых амортизаторов, работающих в области больших упругих деформаций // Труды конф., посвященной 90-летию В.И. Феодосьева. - Май, 2006.

9. Димитриенко Ю. И., Царев С. М. Расчет напряженно-деформированного состояния резино-технических элементов конструкций при больших деформациях // В сб. Аэрокосмические технологии. - М.: Изд-во МГТУ им. Н.Э. Баумана. - 2002. - C. 85-89.

10. Д и м и т р и е н к о Ю. И. Тензорное исчисление. - М.: Высшая школа, 2001. - 575 с.

Статья поступила в редакцию 25.12.2006

Юрий Иванович Димитриенко родился в 1962 г., окончил в 1984 г. МГУ им. М.В. Ломоносова. Д-р физ.-мат. наук, профессор, заведующий кафедрой "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана, действительный член академии инженерных наук. Автор более 100 научных работ в области вычислительной механики, нелинейного тензорного анализа, термомеханики композитов, математического моделирования в материаловедении.

Yu.I. Dimitrienko (b.1962) graduated from the Lomonosov Moscow State University in 1984. D. Sc. (Phys.-Math.), professor, head of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University, Full member of the Russian Academy of Engineering Sciences. Author of more than 100 publications in the field of computational mechanics, nonlinear tensor analysis, thermomechanics of composite materials, mathematical simulations in science of materials.