Конечно-элементное моделирование больших деформаций нелинейно-упругих материалов с использованием модели АV Текст научной статьи по специальности «Физика»

УДК 539.3

Конечно-элементное моделирование больших деформаций нелинейно-упругих материалов с использованием модели Ау

© Ю.И. Димитриенко, А. А. Веретенников МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: конечные деформации, большие деформации, нелинейно-упругие материалы, модель Ар, неогуковская модель, метод конечного элемента.

Введение. В последнее время значительно возросло число исследований, направленных на разработку конечно-элементных методов решения задач нелинейной механики с большими деформациями, отметим лишь некоторые работы в этой области [1-12]. Наиболее известные методы конечно-элементного решения задач нелинейной теории упругости при конечных деформациях основаны на использовании вариационной формулировки в недеформированной конфигурации [12, 13]. В работе [14] Боне был предложен метод, основанный на использовании вариационной слабой формулировки задачи в приращениях и актуальной (деформированной) конфигурации. Для сжимаемой неогуковской модели нелинейно-упругого поведения материала метод достаточно эффективен. Известно, что реализация того или иного численного метода существенно зависит от выбранной модели нелинейно-упругого материала. Целью настоящей работы является применение метода Боне для решения задач нелинейной теории упругости для модели, отличной от неогуковской, — для полулинейной модели Лу, согласно классификации моделей, предложенной в [15-17].

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

ми деформациями для изотропной сжимаемой гиперупругой среды. В актуальной конфигурации уравнение равновесия в дифференциальной формулировке имеет вид [17]

V-Т + рГ = 0, (1)

где Т — тензор напряжений Коши; Г — плотность внешних массовых сил; р — плотность материала; V — набла-оператор. Соответствующий уравнениям (1) принцип виртуальной работы (вариационное уравнение) имеет вид [17]

5Ж(и,5у) = |Т • ЬШУ - |рГ • ЫёУ -11 п • ЫёЪ = 0, (2)

V V I

где и — вектор перемещений; 5у — произвольная виртуальная скорость; V — объем тела в актуальной конфигурации; Б = ® V + V ® уг) —

тензор скоростей деформации; 1п — плотность внешних поверхностных сил (вектор напряжений). Последний интеграл берется по границе Е объема V.

Для задачи (1) выберем два типа определяющих соотношений нелинейной гиперупругой изотропной сжимаемой среды: соответствующие модели А5 по классификации, введенной в [17], и неогуковской модели. Определяющие соотношения модели Ау имеют вид

Т = 1Г •(Я/! (С) + 2|С)), (3)

здесь Г — градиент деформаций; 3 = ёй Г; С — тензор деформаций Коши—Грина; I и | — параметры Ламе материала. Соотношения между вектором перемещений и градиентом деформации и тензором деформации Коши—Грина имеют следующий вид:

Г = Е + ^ иг, (4)

С = 1 [V® и + ^ иг и •V® и^, (5)

- ди

где Е — единичный тензор, а V® и = г' ® ^—- — градиент вектора

перемещений, заданный в отсчетной конфигурации; г' — локальные векторы взаимного базиса.

Определяющие соотношения неогуковской модели выберем в варианте Боне из [14]:

Т = у (В - Е) + ~ (1п У)Е, (6)

где В = Г • Гт — мера деформаций Фингера.

Граничные условия для уравнений (1) выберем в виде заданных перемещений и0 на части Еи поверхности и свободной от нагрузок оставшейся части ЕТ = Е \ Еи поверхности в актуальной конфигурации

и| ^ = ио, п • Т| ^ = 0. (7)

Линеаризация задачи. Задача решалась с применением метода конечных элементов в актуальной конфигурации. В силу постановки задачи (1) в конечных деформациях результирующая система в методе конечного элемента является нелинейной. Для решения системы нелинейных уравнений применялась следующая модификация метода Ньютона: пусть Г(х) = 0 — система нелинейных уравнений, х0 — приближенное решение, а х = х0 + и — точное решение системы, тогда алгоритм метода Ньютона имеет вид

Г(х.) + БГ(х,.) [и] = 0, х,. + 1 = х,. + и, (8)

где БГ(х.) [и] — производная по направлению вектора и, определяемая как слабый дифференциал (дифференциал Гато [18]):

ПГ(х,,и) = ±Г(хо +И5=0 = ¡¡5 Г(х0 + ^-Г(х0).

(9)

В случае, если функционал БГ(х0, и) линеен по то БГ(х0, и) = = БГ(х0)^и. Если же при этом и — вектор из конечномерного пространства, то обычно используют такое обозначение [14]: БГ(х0)^и = = Ш(хо)[и].

Линеаризация принципа виртуальной работы (2) согласно алгоритму (8) метода Ньютона имеет вид

5Ж(ф„ 5у) + £5Ж(фк, 5у) [и] = 0, (10)

здесь фк - к-е приближение решения в перемещениях; и — приращение решения; Б — символ производной по направлению неизвестного вектора перемещений и.

В рассматриваемой задаче внешние поверхностные и массовые силы отсутствовали и граничные условия задавались только в виде перемещений (7). С учетом этого из (9) и (2) можно получить выражение для слабой производной [14]:

БбЖ (, 6у )[и ] = {бБ ••4 С -(¡¿Г + { Т -(®и )-(У®5у )) (11)

V V

Тогда, перенеся в (10) бЩф, 5у) в правую часть со знаком минус, линеаризованный принцип виртуальной работы принимает вид

|бБ--4С• •¡¡¿Г +1Т••( ® и)Т • (V ®5у))) = -[Т• •бDdV. (12)

V V V

Здесь е = 1 (V® и + V® иТ) — тензор деформаций линейной теории

упругости; 4С — тензор упругих постоянных четвертого ранга в пространственном описании, зависящий от материала. Для моделей Лу и неогуковской компоненты тензора 4С соответственно имеют вид

с№ = ((б + Ц5*б + ЦббД ), С™ = -1 ( би + 2ц'б* б ,),

ц' = ц-Я1п 3. (13)

Материальную составляющую линеаризованного уравнения виртуальной работы обозначим следующим образом:

БбШС (,6у)[и] = |бБ ••4С ••(¿¡У, (14)

V

а геометрическую составляющую, соответственно, так

БбШа (,6у)[и] = |Т • •["(V ® и)Т • V ® 6у]¿V. (15)

V

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

и = (X ) , (16)

1=1

где N (Х]) — функции формы; Х] — лагранжевы координаты; qг. — неизвестные перемещения в узлах конечного элемента.

После серии преобразований вклад дискретизированной материальной составляющей (14) в линеаризованное уравнение виртуальной работы на элементе е и узлов а и Ь имеет вид

В5ЖСе (, ^5у в )[ьиь ] =

= |2(- ®™а + ™а ®5уа)4С--2(и, <8^ )ё¥е. (17)

V. 2 2

Вклад геометрической составляющей (15) дискретизируется следующим образом:

Б5Жае (,^5уа) [ьиь] = (5у • и) { • Т• Шь • ЕЖе. (18)

Ve

Вклад эквивалентных узловых усилий имеет вид

5Же(фк,а) = 5уа • |Т(19)

Ve

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

[КА + К]{и} = -{Л}. (20)

Здесь {и} — вектор приращения решения; [КА] и [Ко] — материальная и геометрическая компоненты матрицы жесткости конечного элемента:

[Ка ] = | [В]т • [С] • [В]ОУ, [Кс ] = { [B]т{o}dV, (21)

V V

где [В] — матрица производных функций формы конечного элемента; [С] — матричная запись симметричного тензора 4С упругих постоянных; {а} — столбец компонент тензора истинных напряжений Коши в текущий момент времени (рассчитывается в результате итерационной процедуры); {Л} — столбец эквивалентных узловых сил:

{Л} = | [ В]т •{о}ау.

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

В качестве первого приближения использовались напряжения, полученные при применении к дискретизированной геометрии задачи граничных условий в виде перемещений: полагая решением на нулевой итерации исходную геометрию, измененную на граничные условия, т. е. если х0 — вектор координат узлов в исходной конфигурации, и и0 — вектор приращений координат узлов, задаваемый граничными условиями в виде перемещений, то х01 = х0 + и0 — вектор координат узлов в нулевой итерации. Разница между двумя этими векторами

(х01 - х0 = и0) задает деформацию узлов, по которой можно рассчитать градиент деформации и тензор деформации, а значит, используя (3) или (5), можно получить и компоненты тензора напряжений Коши в выбранной точке конечного элемента.

Алгоритм решения задачи. Алгоритм состоит из двух циклов: цикл по итерациям нагрузки и цикл по итерациям метода Ньютона.

Схема алгоритма следующая.

1. Ввод геометрии, граничных условий и параметров материала.

2. Инициализация х = х0 — глобальный вектор координат узлов.

3. Цикл по итерациям нагрузки:

3.1. обновить глобальный вектор узлов х решением (перемещениями), полученным на предыдущей итерации нагрузки и граничными условиями (хи = х. + и0);

3.2. найти {Я};

3.3. построить глобальную матрицу жесткости;

3.4. реализовать цикл по итерациям метода Ньютона:

3.4.1. построить глобальную матрицу жесткости;

3.4.2. найти вектор Я;

3.4.3. исключить из матрицы жесткости предписанные

перемещения;

3.4.4. решить СЛАУ (20);

3.4.5. вычислить невязку Т • и;

3.4.6. вычислить смещение координат точек х = х + и.

Для решения СЛАУ применялась модификация метода Холецкого с использованием дерева исключения (еНштайоШтее) [19], а также метод сопряженных градиентов с предобуславливателем в виде неполного ЬЦ-разложения (ГШ).

Пример численного решения. Для тестирования предложенного алгоритма была рассмотрена задача об одноосном растяжении бруса, для которой известно аналитическое решение [17]. Были выбраны следующие значения констант задачи:

X = 100 МПа, ц = 100 МПа, Их = 1, И2 = 6, Иъ = 1,

где hi — геометрические размеры бруса в отсчетной конфигурации.

На каждой итерации брусу придавалось приращение по оси растяжения бруса 0,8333 % за 1 итерацию.

На рис. 1 приведен пример решения задачи о растяжении бруса для модели Ау. Показано распределение эквивалентного напряжения в брусе при деформации растяжения 20 %. В виде сетки представлена отсчетная конфигурация бруса до деформации. Цветом обозначено эквивалентное напряжение фон Мизеса (второй инвариант тензора

12(Т)

39.3

78.5

Рис. 1. Распределение интенсивности напряжений 12(Т), МПа, в брусе модели Лу при растяжении до деформации 20 %. Сеткой показана отсчетная

конфигурация бруса

Рис. 2. Распределение интенсивности напряжений 12(Т), МПа, в брусе неогуковской модели при растяжении до деформации 100 %. Сеткой показана отсчетная конфигурация бруса

напряжений 12(Т) [20]. Полученное решение согласуется с известным аналитическим решением [17], относительная ошибка не превышает 0,01 %. На рис. 3 показана диаграмма напряжения — деформации для численного и аналитических решений. Аналогичные графики приведены для неогуковской модели (рис. 2 и 4). Относительная ошибка неогуковской модели не превышает 0,016 %.

Малое количество итераций для модели Лу по сравнению с неогуковской моделью объясняется потерей сходимости метода Ньютона при деформациях выше 20 %. Это явление можно объяснить сильным влиянием нелинейности модели на сходимость метода Ньютона.

Рис. 3. Диаграмма напряжения (фон Мизеса) 12(Т), МПа, — относительное удлинение 5Х бруса при растяжении для численного (точки) и аналитического (сплошная линия) решений в модели Ау

Рис. 4. Диаграмма напряжения (фон Мизеса) 12(Т), МПа, — относительное удлинение 5Х бруса при растяжении для численного (точки) и аналитического (сплошная линия) решений для неогуковской модели

Все исходные коды конечно-элементного комплекса, использовавшегося для данной работы, вместе с исходными данными, доступны по адресу https://github.com/fourier/fea-large. Авторская реализация метода Холецкого с использованием дерева исключения доступна по адресу https://github.com/fourier/libspmatrix.

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

Исследование выполнено при финансовой поддержке Министерства образования и науки Российской Федерации (номер НИР 1.5433.2011) и РФФИ (грант № 12-08-00998-а).

ЛИТЕРАТУРА

[1] Трусов П.В., Швейкин А.И. Теория определяющих соотношений: Ч. II. Теория пластичности. Пермь, ПГТУ, 2008, 243 с.

[2] Коробейников С.Н. Нелинейное деформирование твердых тел. Новосибирск, 2000, 262 с.

[3] Левитас В.И. Большиеупругопластические деформации материалов при высоком давлении. Киев, Наукова думка, 1987, 232 с.

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

[5] Simo J.S., Meschke G. A new class of algorithms for classical plasticity extended to finite strains. Application to geomaterials. Comput. Mech., 1993, vol. 11, pp. 253-278.

[6] Betsch P., Stein E. Numerical implementation of multiplicative elasto-plasticity into assumed strain elements with application to shells at large strains. Comput. Meth. Appl. Mech. Eng., 1999, vol. 179, pp. 215-245.

[7] Ibrahimbegovic A., Gharzeddin F. Finite deformation plasticity in principal axes: from a manifold to the Euclidean setting. Comput. Meth. Appl. Mech. Eng., 1999, vol. 171, pp. 341-369.

[8] Rosati L., Varloso N. A return map algorithm for general isotropic elasto/visco-plastic materials in principal space. Int. J. Numer. Meth. Eng., 2004, vol. 60, pp. 461-498.

[9] Idesman A.V. Comparison of different isotropic elastoplastic models at finite strains used in numerical analysis. Comput. Meth. Appl. Mech. Eng., 2003, vol. 192, pp. 4659—4674.

[10] Simo J.S., Ortiz M. A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations. Comput. Meth. Appl. Mech. Eng., 1985, vol. 49, pp. 221—245.

[11] Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упругопластические деформации: теория, алгоритм, приложения. Москва, Наука, 1986, 232 с.

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

[13] Димитриенко Ю.И., Царев С.М., Веретенников А.В. Разработка метода конечных элементов для расчета элементов конструкций из несжимаемых материалов с большими деформациями. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2007, № 3, с. 69—83.

[14] Bonet J., Wood R.D. Nonlinear Continuum Mechanics for Finite Element Analysis. 2nd edition. Cambridge, Cambridge University Press, 2008.

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

[16] Dimitrienko Yu.I. Novel viscoelastic models for elastomwrs under finite strains. European Journal of Mechanics, A /Solids, 2002, vol. 21 (2), pp. 133—150.

[17] Димитриенко Ю.И. Нелинейная механика сплошной среды. Москва, ФИЗ-МАТЛИТ, 2009, 610 c.

[18] Колмогоров А.Н., Фомин С.В. Элементы теории функций и функционального анализа. Москва, Наука, 1976, 542 с.

[19] Davis T.A. Direct Methods for Sparse Linear Systems. Philadelphia, SIAM, 2006.

[20] Димитриенко Ю.И. Механика сплошной среды. Т. 1: Тензорный анализ. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2011, 463 с.

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

Ссылку на эту статью просим оформлять следующим образом:

Димитриенко Ю.И., Веретенников А.А. Конечно-элементное моделирование больших деформаций нелинейно-упругих материалов с использованием модели А^ Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engjournal.ru/catalog/mathmodel/material/1118.html

Димитриенко Юрий Иванович родился в 1962 г., окончил МГУ им. М.В. Ломоносова в 1984 г. Д-р физ.-мат. наук, профессор, заведующий кафедрой «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана, директор Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана (НОЦ «СИМПЛЕКС»), действительный член академии инженерных наук. Автор более 250 научных работ в области механики сплошной среды, вычислительной механики, нелинейной механики деформируемых сред, термомеханики композитов, математического моделирования в науке о материалах, моделирования в экономике. e-mail: dimit.bmstu@gmail.com; txm.fourier@gmail.com

Веретенников Алексей Анатольевич родился в 1983 г., окончил МГТУ им. Н.Э. Баумана в 2005 г. Аспирант кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор научных работ, посвященных применению метода конечных элементов для решения задач механики при конечных деформациях.