Сравнительный анализ методов интегрирования уравнений нелинейной теории упругости Текст научной статьи по специальности «Математика»

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

1. Болотин Ю.В., Голован А.А., Парусников Н.А. Уравнения аэрогравиметрии. Алгоритмы и результаты испытаний. М.: Изд-во ЦПИ при мех.-мат. ф-те МГУ, 2002.

2. Fengler M.J., Freeden W., Gutting M. Multiscale modeling from EIGEN-1S, EIGEN-2, EIGEN-GRACE01S, GGM01S, UCPH2002_0,5, EGM96: wavelet coefficients, variances and reconstruction // Proc. 2nd CHAMP Science Meeting. Berlin; Heidelberg; N. Y.: Springer, 2004. 145-150.

3. Moritz H. Advanced physical geodesy. Karlsruhe: Herbert Wichmann Verlag, 1980.

4. Freeden W., Schneider F. Wavelet approximation on closed surfaces and their application to boundary-value problems of potential theory // Math. Meth. Appl. Sci. 1998. 21. 129-163.

5. Матасов А.И. Метод гарантирующего оценивания. M.: Изд-во МГУ, 2009.

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

УДК 511

СРАВНИТЕЛЬНЫЙ АНАЛИЗ МЕТОДОВ ИНТЕГРИРОВАНИЯ УРАВНЕНИЙ НЕЛИНЕЙНОЙ ТЕОРИИ УПРУГОСТИ

М. В. Козлов1, С. В. Шешенин2

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

Ключевые слова: резинокорд, композит, нелинейность, вариационное уравнение, линеаризация, метод конечных элементов, метод Эйлера, метод Ньютона, метод Рунге-Кутты, многошаговый метод.

In this paper we describe the way of practical solution of the variational equation in geometrically and physically nonlinear problems of deformable body mechanics. A system of a large number of nonlinear ordinary differential equations usually appears in such problems. The Euler method is typically used in this approach. We propose to use a Runge-Kutta method and multistep methods and consider the solution complexity in terms of computing cost to find a method that provides a more efficient solving procedure of nonlinear problems.

Key words: rubber-cord, composite, nonlinearity, variational equation, linearization, finite element method, Euler method, Newton method, Runge-Kutta method, multistep method.

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

В настоящей работе рассматривается общий вид квазистатической краевой задачи в начальной области в случае, когда нагрузка S0 не зависит от вектора перемещения (случай "мертвой нагрузки"):

( 0 0 -V P(и) + Р f = 0,

< P(и) ■ N(X)| 0 = S0(X,t), (1)

~ £2

и| 0 = 0,

£1

1 Козлов Михаил Владимирович — студ. каф. механики композитов мех.-мат. ф-та МГУ, e-mail: my_skylineQmail.ru.

2 Шешенин Сергей Владимирович — доктор физ.-мат. наук, проф. каф. механики композитов мех.-мат. ф-та МГУ, e-mail: sergey.shesheninQmail.ru.

0 0 где V — градиент в начальной конфигурации, р — плотность в начальной конфигурации, Р — первый

тензор Пиолы-Кирхгофа, / — массовые силы. Определяющие соотношения уже учтены, и тензор Пио-

0

лы Р есть оператор от вектора перемещения и. На части поверхности £1 задано нулевое перемещение.

0 0 0 0 Граничные условия на части поверхности £2 заданы в напряжениях, причем £ = £1 + £2- При этом

N — нормаль в начальной конфигурации, X = {XI, Х2, Х3} — радиус-вектор в начальной конфигурации.

Задаче (1) соответствует вариационное уравнение

У р(и): V = у Р /{Х,г) ■ + у 50(X,Í) ■ £ , (2)

0 0 0 V V £2

0

где и — пробная функция, (V = dXldX2dXз.

Линеаризованное по параметру нагружения Ь вариационное уравнение имеет вид

{ 0 0 Г 0 - - 0 Г _ 0

/ (р(и) : V и (V = Р (/(X, Ь) ■ wdV + (Б0(Х,Ь) ■ £ .

0 0 0 V V £2

Последнее соотношение удобно записать с использованием тензора Кирхгофа 5 [1] (второй тензор

Пиолы-Кирхгофа). Он связан с градиентом деформации Р [1] соотношением Р = Р ■ Б. В декартовых

координатах линеаризованное вариационное уравнение с использованием тензора Кирхгофа Б принимает вид

У С™ (иРд( V + У (иРд( V = ^ Р (/г ' ( V + ! (Б,0 ' (£, (3)

0 0 0 0 V V V £2

где

дБ

Г,ТАК_Г,Е ГС Г. с. ГЕ _

гзрд ~ — 01'р0]Ч1 ^ — —

и Е — тензор деформаций Лагранжа-Грина [1]. В уравнении (3) все дифференциалы взяты относительно параметра Ь.

К уравнению (3) применяется метод конечных элементов [2, 3]. Обозначив С= СТ^ + перепишем (3) следующим образом:

[ , , 0 Г 0 ■ 0 Г . п 0

/ (С$а(и) иЦ (V = Р 1к ик (V + Б0 'к (£ , (4)

0 0 0 V V £2

где точка обозначает производную по Ь.

Пусть N — число узлов сетки. В используемом подходе функции и, и и г аппроксимируются в виде сумм глобальных базисных функций ф:(X) с коэффициентами, равными соответствующим узловым значениям и/, и)1. Необходимые представления имеют вид [1-3]

N N N

йг(x,t) = Y,Ф1 (X) и (ь), щ^^) = ^2ф:(X) и(Ь), иг(X) = Ф1 (X)

1=1 1=1 1=1

Для частных производных справедливы представления

N N

(X, = Фи (X) щ (Ь), иг,3 (X) = £ фи (X) Щ, 1=1 1=1

подставляя которые в уравнение (4) получим

Пусть

Киы(и) = J Ск$й(и) ф1>г фз^ йУ, доз -0 V

С учетом последних обозначений имеем соотношение

N 3 N 3

^,J=1 г,к=1 3=1 0=1

которое должно выполняться для любой пробной функции. Поэтому справедлива система из ЗУ линейных

-1

уравнений относительно и :

N 3

щ = доз, J = 1,...,У, к = 1,2,3. (5)

1=1i=1

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

К (и)Ц7 + ^ = 0. (6)

Матрица К в (6) называется глобальной матрицей жесткости.

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

N (и) + ^ = 0, (7)

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

Уравнение (7) всегда может быть получено для определяющего соотношения 5 = Б(Е). Для решения

уравнения (7) часто используют метод Ньютона [2, 4], метод поиска вдоль линии (усовершенствованный метод Ньютона) [4] или метод длин дуг [5].

Однако, если существуют производные Гато определяющего соотношения, система алгебраических уравнений (7) может быть приведена к виду (6), тогда можно использовать методы, разработанные для решения систем обыкновенных дифференциальных уравнений (ОДУ). Далее речь пойдет о таких методах, применяемых к (6).

Приведем (6) к форме и = —К-1(иобратив матрицу К. В вычислительном алгоритме нет необ-

К

соотношение в теории ОДУ принято записывать кратко:

17 = {(г, и). (8)

Явный метод Эйлера [6], примененный к (8), выглядит следующим образом: ип+1 = ип+ДгГ(гп, ип), где Дг — шаг параметра нагружения.

Неявный метод Эйлера [6] для решения задачи имеет вид ип+1 = ип + ДИ(гп+1, ип+1). В общем случае вектор ип+1 не может быть явно найден из уравнений, поэтому на практике используется схема с несколькими итерациями по параметру з, состоящая го соотношений ип+1'° = ип + Дг{(гп, ип) и и«.+М+1 = ип + ДгГ(гп+1, ип+1'5).

Метод трапеций [6] является неявным и имеет вид ип+1 = ип + Дг Г(гп, ип) + {(гп+1, ип+1)) /2. При его реализации применяется итерационная схема го соотношений ип+1'° = ип + ДгГ(гп, ип) и ип+1'5+1 =

ип + Дг Г (гп, ип) + г (ьп+1, ип+1>*)) /2.

00 Р ¡к Фз й У +

5 0 Фз й Е

0 V

0 £2

Обозначим Г = Гиг). Явный четырехшаговый метод Адамса-Башфорта [6] записывается в форме ип+1 = ип + § (55£та - 59£га_1 + 37^-2 ~

Метод Адамса-Мултона [6] того же порядка точности является неявным и имеет вид ип+1 = ип + (9{п-\-\ + 19£п — 5£п-1 + ¡п-2)- На практике можно применять схему, где в первой итерации используется метод Адамса-Башфорта

Г ип+1,о = ип + ^ (55£та - 59^-1 + 37(п.2 - 9(п.3), | ^п+1,,+1 = ип + ^ (9?*+1 + Жп - Ып-1 + ,

здесь £П+1 = Г(Ьп+1, ип+1'°).

Метод Рунге-Кутты четвертого порядка [6] имеет вид ип+1 = IIй + ^ (Ау\ + 2Ау2 + 2Дуз + Ау^),

где

Ауг = ип), Ау2 = Aíf + ип + ,

Ауз = А1( + + , Ду4 = Д^п + Д*, С/га + Ау3).

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

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

В таблице приводятся характеристики описанных методов решения систем ОДУ. Для неявного метода Эйлера и метода трапеций используются реализации с тремя итерациями. В стоимости решения для многошаговых методов учитываются начальные шаги, на которых применяется метод Рунге-Кутты.

Характеристика метода Метод решения систем ОДУ

Явный метод Эйлера Неявный метод Эйлера Метод трапеций Метод Адамса- Вашфорта Метод Адамса- Мултона Метод Рунге-Кутты

Стоимость одного шага 1 3 3 1 2 4

Стоимость решения из М шагов М 3 м 3 М М + 9 2М + 6 4М

Порядок локальной ошибки 0(Д£2) 0(Д£2) 0(Д£3) 0(Д£б) 0(Д£б) 0(Д£б)

Порядок глобальной ошибки 0(Д£) 0(Д£) 0(Д£2) 0(Д£4) 0(Д£4) 0(Д£4)

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

Целью исследования являлось отыскание численного метода с самой низкой стоимостью решения при наименьшем уклонении от точного, иными словами, с самой низкой стоимостью глобальной ошибки. В качестве "точного" использовалось решение, полученное с шагом, в 40 раз меньшим, чем типичный шаг в описанных выше методах.

В ходе эксперимента пластина нагружалась одной и той же силой. При этом решение задачи отыскивалось с различным числом шагов по времени (1, 2, 3,...) до тех пор, пока его стоимость не превышала 400 единиц. Для каждого метода вычислялись абсолютные величины относительных ошибок по завершении нагружения. Величина, которая подвергалась сравнению во всех случаях, — прогиб деформированной пластины. По осям графика откладывались модуль относительной ошибки и общая стоимость решения.

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

На рис. 2 приведены кривые для всех вышеупомянутых методов. По ним можно получить представление о стоимости определенной величины погрешности для каждого метода на примере реальной задачи.

са Ы VO К

3 о

4

сз

а л

ч

10"

10"

10"

Й 10

-4

К о О X н

О

10

-5

10"6 -

10

га 10": и ю в

,-1

[Т

\\ N

\Ч " \ ч

10"3

; ю-4

_L

_L

_L

3

я л

4

о н К о о аа

О 10"5

10"

Метод

■ Адамса-Мултона

■ Рунге-Кутты

■ трапеций

■ Адамса-Башфорта неявный Эйлера

■ явный Эйлера

\\

_L

_L

100 200 300

Стоимость

400

100 200

Стоимость

300

400

Рис. 1. Степенная аппроксимация экспериментального графика

Рис. 2. Стоимость относительной ошибки

Методы высоких порядков оказались значительно эффективнее. Лучшую скорость убывания ошибки продемонстрировал мнш'ошах'овый метод Адамса Башфорта. Однако так как метод является явным, то мш'ут возникать проблемы, связанные с устойчивостью решения. Поэтому для быстршх) и эффективншх) решения нелинейных задач больше вемх) подходят метод Руш'е Кутты четвертшх) порядка и мнш'ошах'о-вый метод Адамса Мултона.

В экспериментах первоначально решались линейные системы из 945 уравнений. Затем та же самая задача была решена на более мелкой сетке с использованием 4305 уравнений. Оказалось, что картина в этих случаях качественно одна и та же.

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

Работа выполнена при финансовой поддержке РФФИ, грант № 13 01 00688.

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

1. Bonet J., Wood R. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge: Cambridge Univ. Press, 1997.

2. ПобеОря Б.Е. Численные методы в теории упругости и пластичности. М.: Илд-во МГУ, 1995.

3. Зенкевич О. Метод конечных элементов в технике. М.: Мир, 1975.

4. Kelley С. Т. Iterative Methods for Linear and Nonlinear Equations. Philadelphia: SIAM, 1995.

5. Метоп B.A., Su X.Z. Arc-length technique for nonlinear finite element analysis // J. Zhejiang Univ. SCI. 2004. 5, N 5. 618 628.

6. Butcher J.CI Numerical Methods for Ordinary Differential Equations. Chichester: John Wiley & Sons, 2008.

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