УДК 539.2:544.2:678.01:519.7:539.3:517.958
ЭВОЛЮЦИОННЫЕ ИЗОТРОПНЫЕ СХЕМЫ ЧИСЛЕННОГО РЕШЕНИЯ НЕЛИНЕЙНЫХ КРАЕВЫХ ЗАДАЧ КВАЗИСТАТИЧЕСКОГО ДЕФОРМИРОВАНИЯ. Часть 5. Особенности пошагового решения проекционно-сеточных систем при конечных деформациях
АЛЬЕС М.Ю.
Институт механики Уральского отделения РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. Исследуются эволюционные схемы численного решения систем нелинейных проекционно-сеточных уравнений неупругого поведения полимерных изделий в условиях конечных деформаций. Анализируются пошаговые методы.
КЛЮЧЕВЫЕ СЛОВА: эластомерные композиты, нелинейное механическое поведение, конечные деформации, краевые задачи, метод конечных элементов, эволюционные численные схемы.
Трудности, возникающие при решении полученных в [1] систем нелинейных проекционно-сеточных уравнений большого порядка идентичны тем, что подробно обсуждались в [2]. К проблемам, связанным с существенно нелинейным механическим поведением рассматриваемых материалов, добавились существенные трудности, обусловленные учетом конечности деформаций. Из сравнения проекционно-сеточных систем [1, 2] видно, что уравнения не только существенно усложнились, но произошли структурные изменения. В слагаемом Ь(и)и [1] отсутствует вклад от шаровой составляющей тензора напряжении. При рассмотрении задачи в рамках малых деформаций [2, 3] относительное изменение объема 9 определяется первым инвариантом 11 (е), что позволяет внести вклад члена K9gij в матрицу Ь(и) системы и по возможности предусмотреть
уменьшение нормы || ¿ (и)||. При конечных деформациях относительное изменение объема 9 определяется третьим инвариантом тензора деформации Коши-Грина 13 (О) [4] и возможно несколько вариантов выделения линейности части в виде произведения некоторой функции на дивергенцию вектора перемещений. Эта процедура, с одной стороны отнюдь не однозначна с точки зрения влияния на сходимость алгоритма в условиях, когда объемные изменения могут быть существенными, существенно нелинейными и неоднородными по области. С другой стороны - по существу во многом предопределяет схему разрабатываемого вычислительного алгоритма и сужает возможности маневра при его отладке, или в ситуации, когда, например, при расчетах возникает необходимость изменения модели объемного поведения. Отмеченное обстоятельство характерно и по отношению к сдвиговым характеристикам (если член ЦЭЭц в модели явно не выделяется). "Лобовая атака"
краевых задач неупругого поведения полимерных изделий сложной формы в условиях больших деформаций и при сложных программах нагружения, как правило, безнадежна. Предварительные же теоретические исследования сходимости алгоритма при учете конечности деформации еще более затруднены по сравнению с задачами в геометрически линейном приближении. Все это, с одной стороны, сильно затрудняет получение численных результатов и использование сложных нелинейных моделей деформирования в практических расчетах, а с другой, на основе анализа полученных результатов - разработку новых, более совершенных и физически содержательных моделей деформирования.
Существует достаточно большое число методов решения уравнений, получающихся при применении метода конечных элементов к задачам с учетом геометрической нелинейности (см., например, [5 - 14]). Всестороннее исследование применения МКЭ к нелинейным упругим задачам при больших деформациях дано Оденом [5]. По-видимому,
наибольшее применение в нелинейных задачах механики деформируемого упругого и упруго-пластического тела нашли пошаговые методы, которые позволяют исследовать процесс деформирования при изменении внешней нагрузки и тем самым контролировать (путем сравнения с решением аналогичной линейной задачи, например) достоверность (физичность) результатов расчета, точки бифуркации решения.
В основе значительной части пошаговых вычислительных алгоритмов решения нелинейных задач лежит метод продолжения по параметру [10,11], заключающийся в "погружении" исходной системы
Н (и) = 0 (1)
в однопараметрическое семейство уравнений
Н (и,£ ) = 0 (2)
таким образом, что при £ = 0 система (2) имеет заведомо известное решение и(£ = 0) = и0, при X = 1 переходит в исходное уравнение (1) Н (и, 1) ° Н (и). Далее считая, что и непрерывно дифференцируемо на отрезке [0,1] а отображение Н(и, X имеет непрерывные частные производные диН(и, £), д^ Н(и, £), исходная задача (1) путем дифференцирования (2) по X сводится к последовательности задач Коши
д^и + д~и Н (и, X д^Н (и, X ): : 0, и
V
0
0
"X е [0,1]
(3)
с начальными условиями: X
Разнообразие имеющихся пошаговых алгоритмов обусловлено возможностью широкого выбора различных схем численного интегрирования уравнения (3). Используя, например, простейшую явную схему, получим широко применяемый в нелинейных задачах упругости и упруго-пластичности, теории оболочек, метод последовательных нагружений
[5].
Применим метод последовательных нагружений для решения системы Ь(и)и + Г (и) — О — р(и) = 0 [1]. С целью парирования накапливающейся погрешности вычислений в разностный аналог задачи Коши (3) на шаге интегрирования (е +1) введем дополнительный член, равный невязке системы (2) на предыдущем шаге с некоторым весовым коэффициентом [9, 12]. Имеем
е е+1 е ее е+1
А А и = / —деН (и, X )Л X, е = 0, Е — 1
А = д„
Ь(и ) и
— дuP(u)X+дuF (и),
I =
О + Р(и)
АX,
Н (u,X ) = Ци ) и —X
О + Р(и)
+ Г (и),
е+1
£=£+А£
е+1 е+1
и
и + А и
(4)
где Е - число интервалов разбиения области изменения X. Поскольку тело в
ненагруженном состоянии находится в р -конфигурации, то и(£ = 0) = и = 0 . Аналитическое отыскание выражений для частных производных ди [Ь(и )и], диГ (и) затруднено (см. (25),(26)
в [1]). При расчетах, как правило, используется простейшая явная конечно-разностная аппроксимация.
е
ее
е
е
е
Здесь отметим следующие важные моменты. Популярность метода последовательных нагружений в практических приложениях в немалой степени обусловлена ясной физической интерпретацией каждого шага решения в задачах нелинейного деформирования упругих тел. "Погружение" системы (1) в семейство (2) естественно произвести, умножив вектор нагрузки на параметр X/. Приращение параметра X/ соответствует тогда приращению нагрузки.
Продолжая этот процесс до момента X/ = 1, имеем полную картину нелинейного упругого поведения тела при переходе из конфигурации р в р. Очевидно, однако, что при рассмотрении процесса деформирования, обусловленного воздействием потенциального силового поля (например, при действии сил тяжести) физический смысл имеет решение только при X/ = 1. Требует некоторого внимания также расчет напряженно-
деформированного состояния этим методом в условиях немонотонного нагружения и зависимости свойств среды от предыстории деформирования. Программа нагружения разбивается на временные интервалы t0 ..Лх, ^ ..Л2, t2 ...t3, ..., в пределах которых нагрузка
изменяется монотонно. Задача (3) решается для каждого отдельного интервала, при этом, начальным условием для следующего интервала является решение, полученное в последней точке (X/ = 1) нагружения на предыдущем.
Не следует переоценивать простоту и универсальность подобных приемов для любых задач деформирования. Расчеты показали, что применение метода последовательных нагружений с коррекцией (4) для решения системы уравнений Ь(и)и + Г (и) - G - р(и) = 0 [1] дает хороший результат только на участке нагрузки в первом цикле деформирования и только до момента достижения некоторого критического уровня деформаций ес » 0,6, при превышении которого накопление погрешности начинает очень быстро возрастать. Некоторые иллюстрирующие сходимость результаты для модельной задачи растяжения-сжатия стержня показаны на рис. 1, 2.
0.5 '
Рис. 1. Режим нагружения и изменение относительной евклидовой нормы невязки
£ =
Н (иЛ)
Н (u„Xf)
, * - решение, полученное сходящимся методом (5)
При деформациях меньших, чем ес наилучшая сходимость достигалась при
^Xf » 0,8... 1,2 (рис. 2).
2
2
Рис. 2. Влияние корректирующего множителя на норму невязки
На участке разгрузки процесс (4) начинает расходиться сразу, даже в условиях малых деформаций (е < 0,1). Не удается получить удовлетворительного результата и при повторном нагружении.
С целью парирования трудностей, обусловленных возможным превышением деформациями вышеуказанного критического уровня ес » 0,6, в [13] предложен алгоритм
промежуточных отсчетный конфигураций, заключающийся в организации процесса относительно приращений перемещений. Алгоритм показал хорошую сходимость при простых нагружениях, однако, для немонотонных режимов удовлетворительной сходимости добиться не удается.
Теоретические оценки условий сходимости схемы (4) получить затруднительно. Наблюдаемые же счетные эффекты сходимости объясняются следующим образом. Как только нелинейность системы (1) начинает проявляться наиболее сильно, так представление решения в виде последовательности кусочно-линейных шагов (4) приводит к недопустимому его отклонению от траектории, описываемой системой параметрических уравнений (2). На этапе активного процесса нагружения (в первом цикле) нелинейные эффекты системы уравнений (2) определяются достигнутым уровнем деформаций. При разгрузке и повторном нагружении до момента достижения деформациями значений || е ||¥ сразу же, даже в
условиях малых деформаций, очень сильно проявляется нелинейность, связанная с механическим поведением.
Для уменьшения погрешности, накапливающейся при интегрировании уравнения (3),
можно воспользоваться каким-либо итерационным методом, использующим точку и , определяемую системой (4), в качестве начального приближения. Однако, при этом, неминуемо возникают те же трудности, которые подробно обсуждались в [2]. При этом к проблемам, связанным с нелинейностью механического поведения материала, добавляются проблемы геометрической нелинейности. Предложенные в [15 - 17] подходы к построению алгоритмов решения физически нелинейных задач можно эффективно распространить на задачи с учетом конечности деформаций. Записав модель деформирования в следующем,
г
достаточно общем виде Р'] = { ек 1 (X), Л(%) } и введя вместо данного соотношения
Х=о
следующее эволюционное уравнение
Р* = Щыдх еы + Я* { ек,(£),?(£)}
(5)
Х=о
Х= о
-ы
0
можно построить эффективный вычислительный алгоритм, сходящийся при немонотонных режимах нагружения, нелинейной связи между напряжениями и деформациями и конечных деформациях. В (6): Рг, еК1 - симметричный тензор напряжений Пиола-Киргоффа и тензор конечных деформаций Грина; Л(%) - совокупность скалярных параметров, характеризующих процесс; - некоторый функционал; Щы = 2т1Ег]к1 + (К/ -1 | )&1 ЕЫ, | > 0, Кг > 0, Кг - у | > 0, |, Кг - параметры алгоритма, размерность которых совпадает с размерностью модулей упругости; Е']ке = 0,5(glkg1 1 + g11 glk) - единичный тензор четвертого ранга; X/ - фиктивное безразмерное время. Уравнениями (5) вводится некоторый
переходный процесс «изменения» связи между напряжениями и деформациями. Физический смысл имеет только стационарное решение. Некоторые иллюстрирующие сходимость результаты для модельной задачи растяжения-сжатия стержня показаны на рис. 3.
а)
8а = 1\И(и)||¥ : 1) зе = 0 ; 2) зеДХ/ = 1; 3) - линеализированный метод наложения малой деформации на предварительно нагруженное тело [14]; 4) - алгоритм (5)
г* I >
б.. -
ол
о, л
1
к
г иЧ
\\ ' V*; Ч _ * ,___
б) изменение относительной ошибки при итерациях по алгоритму (5)
1Л
Рис. 3. Сходимость алгоритмов при заданном режиме деформирования
Детальное изложение вычислительных алгоритмов на основе (5) будет дано в следующей статье.
СПИСОК ЛИТЕРАТУРЫ
1. Альес М.Ю. Эволюционные изотропные схемы численного решения нелинейных краевых задач квазистатического деформирования. Часть 4. Проекционно-сеточная постановка при конечных деформациях // Химическая физика и мезоскопия. 2014. Т. 16, № 2. С. 195-201.
2. Альес М.Ю. Эволюционные изотропные схемы численного решения нелинейных краевых задач квазистатического деформирования. Часть 1. Особенности сходимости нелинейного численного решения в условиях малых деформаций // Химическая физика и мезоскопия. 2013. Т. 15, № 3. С. 337-342.
3. Альес М.Ю. Феноменологические описания нелинейного сопротивления деформированию высоконаполненных эластомерных (нано) композитов. Часть 1. Малые деформации // Химическая физика и мезоскопия. 2010. Т. 12, № 1. С. 69-77.
4. Альес М.Ю. Феноменологические описания нелинейного сопротивления деформированию высоконаполненных эластомерных (нано) композитов. Часть 2. Конечные деформации // Химическая физика и мезоскопия. 2010. Т. 12, № 2. С. 219-223.
5. Оден Дж. Конечные элементы в нелинейной механике сплошных сред. М. : Мир, 1976. 464 с.
6. Поздеев А.А., Трусов В.П., Няшин Ю.И. Большие упругопластические деформации: теория, алгоритмы, приложения. М. : Наука, 1986. 231 с.
7. Григолюк Г.И., Шалашилин В.И. Проблемы нелинейного деформирования: метод продолжения по параметру в нелинейных задачах механики твердого тела. М. : Наука, 1988. 232 с.
8. Хофмейстер Л., Гринбаум Г., Ивенсен Д. Упругопластический расчет больших деформаций методом конечных элементов // Ракетная техника и космонавтика. 1971. Т.9, № 7. С.42-51.
9. Стриклин Д., Хейслер В., Риземанн В. Оценка методов решения задач строительной механики, нелинейность которых связана со свойствами материала и (или) геометрией // Ракетная техника и космонавтика. 1978. Т. 11, № 3. С. 45-56.
10. Ортега Дж., Рейнболдт В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. М. : Мир, 1975. 560 с.
11. Давиденко Д.Ф. Об одном новом методе численного решения систем нелинейных уравнений // Доклады академии наук СССР. 1953. Т. 88, № 4. С. 601-602
12. Темис Ю. М. Метод последовательных нагружений с коррекцией погрешности в геометрически нелинейных упругих задачах // Всесоюзный межвуз. сборник «Прикладные проблемы прочности и пластичности. Алгоритмизация решения задач упругости и пластичности». Горький : Изд-во Горьковского государственного университета, 1980. Вып. 16. С. 3-10.
13. Альес М.Ю., Булгаков В.К., Липанов А.М. Об одном алгоритме решения геометрически нелинейной задачи о напряженно-деформированного состоянии полых цилиндров сложной формы на основе метода конечных элементов // Известия Академии наук СССР. Механика твердого тела. 1985. № 2. С. 106-112.
14. Прикладные методы расчета изделий из высокоэластичных материалов // под ред. Э.Э. Лавендела. Рига : Зинатне, 1980. 238 с.
15. Альес М. Ю. Эволюционные схемы численного решения нелинейных краевых задач неньютоновских ползущих течений // Химическая физика и мезоскопия. 2013. Т. 15, № 2. С. 193-200.
16. Альес М.Ю. Эволюционные изотропные схемы численного решения нелинейных краевых задач квазистатического деформирования. Часть 2. Метод асимптотической сходимости для нелинейных сжимаемых сред в условиях малых деформаций // Химическая физика и мезоскопия. 2013. Т. 15, № 4. С. 502-507.
17. Альес М.Ю. Эволюционные изотропные схемы численного решения нелинейных краевых задач квазистатического деформирования. Часть 3. Метод асимптотической сходимости для нелинейных несжимаемых сред в условиях малых деформаций // Химическая физика и мезоскопия. 2014. Т. 16, № 1. С. 25-30.
EVOLUTIONARY ISOTROPIC NUMERICAL SCHEMES FOR ACTION BOUNDARY VALUE PROBLEMS OF NONLINEAR QUASI-STATIC DEFORMATION. PART 5. FEATURES STEP SOLUTIONS PROJECTION-GRID SYSTEMS AT FINITE DEFORMATIONS
Alies M.Yu.
Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. We investigate the evolutionary scheme of numerical solutions of systems of nonlinear equations projection-grid inelastic behavior of polymeric products under finite deformations. Analyzed step by step methods.
KEYWORDS: elastomer composites, nonlinear mechanical behavior, finite deformation, boundary value problems, finite element method, numerical evolutionary scheme.
Альес Михаил Юрьевич, доктор физико-математических наук, профессор, главный научный сотрудник ИМ УрО РАН, тел. 89128563824, e-mail: [email protected]