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

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

Том 149, кл. 4

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

2007

УДК 539.3

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ СИЛЬНОГО НЕЛИНЕЙНОГО ДЕФОРМИРОВАНИЯ В КООРДИНАТАХ ЭЙЛЕРА

М.С. Агапов, Е.Б. Кузнецов, В.И. Шалашилип

Аннотация

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

Введение

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

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

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

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

1. Постановка задачи

Рассмотрим напряженно деформированное состояние среды, занимающей область П. Уравнения равновесия среды, отнесенные к криволинейным координатам Эйлера (х1,х2,х3), в случае отсутствия внешних массовых и поверхностных моментов взаимодействия, будут иметь вид [4. 5]

У оц + р^ = 0, = а(1.1)

Здесь р - плотность среды, ¥® — компоненты вектора массовых сил, а13 - кон-травариаптпые компоненты тензора напряжений, УЦ — символы ковариантного дифференцирования

датп

У^тп | ^.шгт^п . „.пз рш

-к

рического тензора д при помощи формул

где Г3 - символы Кристофеля второго рода, вычисляемые через компоненты мет-

Г« = ^ + ~ ^ , т, п = 1, 2, 3.

13 2 \дхЗ дх1 дхп

Если в системе криволинейных координат заданы ковариантные компоненты вектора перемещения и точек среды при деформации, то составляющие тензора конечной деформации вц в координатах Эйлера будут вычисляться по формулам [6, 7] 1

еЧ = + У]иг - У.гД/.пУ^ы"), (1.2)

где символы ковариантного дифференцирования определены соотношениями

__ ди з _„ __ д:иР _о ^

= д£~ Г>п' УгП = ^ +г'»;ы •

Ограничимся физическими соотношениям, которые могут быть заданы в дифференциальной форме

а = Л1зтп 1втп. (1.3)

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

Полагаем, что все неизвестные, входящие в уравнения равновесия (1.1). геометрические соотношения (1.2) и физические соотношения (1.3). зависят от некоторого параметра А, смысл которого определим позже, то есть агз = агз (А), е^ = е^ (А), и = и(А), г,з = 1, 2, 3.

Введем следующие обозначения

¿аг (1вг! . ¿и

—— = а1'3, -г- = е.И, — = и. (1.4

¿X ' ¿X 3 ¿X у '

Чтобы получить уравнения, определяющие поле скоростей изменения неизвестных а, в, и, продифференцируем соотношения (1.1), (1.2) по А, а физические соотношения (1.3) разделим па ¿А. В результате, считая плотность среды неизменной, получим следующую систему уравнений, линейную относительно скоростей напряжений, деформаций и перемещений:

д¥г

У^+р—= 0, о»=сг>%, (1.5)

^(У;,«^ + Vз'Щ - ЧгйпЧ^и" - У^ЫпУ^м"), (1.6)

2

аго = лг^тпетп. (1.7)

Обратимся теперь к граничным условиям. Обозначим через ^1,^2, ^з компоненты единичного вектора нормали к поверхности Я = Я1 + Я2 - границе области В. Тогда граничные условия могут быть записаны в виде

и = и(0) на Я1, агз = / г на Я2, (1.8)

где ]г — компоненты полного вектора внешней силы.

А

¿и- д{г

и = 0 на Й'ь сг%3щ■■ + а*1=-¿г "а (1.9)

¿А дА

Выберем параметр А так, чтобы при А = 0 было известно напряженно-деформированное состояние среды а = а(0\ в = е(0), и = и(0). Обычно это неде-формпрованное состояние а(0) = 0, е(0) = 0, и(0) = 0. Тогда обозначения (1.4), дополненные начальными условиями, можно рассматривать как начальную задачу для системы обыкновенных дифференциальных уравнений

¿а ¿в ¿и .

Тх=^ Тх=е' Ж = и (1Л0)

а|л=0 = а(0), е|л=0 = е(0), и|л=0 = и(0).

Входящие в правые части этих уравнений величины скоростей <г, е, и по пара-А

(1.7) при граничных условиях (1.9) на деформированной поверхности среды. Эти

А

система (1.5)—(1.7) линейна относительно функций <г, е, и, то ее удобно свести к уравнениям для скоростей и изменения перемещений и. Эти уравнения примут

Ь(а,е, и)и = Ф, (1.11)

где Ь - линейный оператор, определенный в деформированной среде, то есть в области В.

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

Г(а,в, и)и = 7, (1.12)

где Г - линейный оператор, определенный на поверхности деформированной среды, то есть на границе S области D.

Тогда параметр Л можно определить, как в [9], функционалом

d,X2 = ~jj J(du> du)dD или J(ii, u)dD = 1, (1.13)

D D

где оператор (,) задает скалярное произведение.

Функционал (1.13) является обобщением на гильбертово пространство наилучшего параметра продолжения [2, 3] в евклидовом пространстве.

2. Одномерная задача

Уравнения равновесия (1.1) и выражения для компонентов тензора конечной деформации (1.2), записанные в ортогональной декартовой системе координат Эйлера (xi,x2,x3), примут вид

daij ■ ....

— +PF>=0, а*=сг>\ (2.1)

- I (, _ ди^диЛ

€ij ~ 2 \дхЗ дхi дхi dxi) ' 1 ' '

Ось ж1 надраим вдоль оси стержпя x. Примем квадратичное приближение в форме, предложенной в [1], т.е. будем считать относительное удлинение £i рав-

в11

но вполне корректно, но оно допускает точное решение, что позволяет использовать его для тестирования численных методов. Тогда для простейшего закона Гука a ii = Евц после введения обозначений fi = pFi, u = ui уравнения (1.5)— (1.7) можно записать в форме

¿¿а . dfx . ( du\ dii .

—,--Ь —= 0, ец = 1 — —— ——, сг ц= Sen. (2.3)

dx dЛ \ dx J dx

Граничные условия для стержня, один конец которого закреплен, а другой нагружен напряжением a, можно задать в виде

U(0) = 0, aii (xr) = a. (2.4)

Здесь a = da/dЛ xr - координата незакрепленного конца стержня, нагруженного a

Для случая Fi =0 получим точное решение этой задачи. Из первого уравнения (2.3) следует, что функция an не изменяется вдоль координаты x, то есть an = = const.

С учетом второго граничного условия (2.4) имеем a ii = a. xi = x

du 1 f du\ 2

611 = dx ~ 2 \dx) ' (2'5)

А

условий принимает вид

<711 = Ее11.

¿и

Разрешив уравнение (2.5) относительно — . получим

¿х

^ = (2.6)

В равенстве (2.6) следует взять знак минус. Выбор знака обусловлен физическим смыслом: перемещение должно возрастать при увеличении деформации е11 а

Отсюда, имея в виду, что и(0) = 0 и ец = — = е, заключаем:

Е

и= ( 1_ </1 ж.

Для стержня, начальная длина которого была равна I, получим перемещение правого конца в виде

¿(хг) = (1 - \/1 - 2е) (/ + м(.г'г)).

хг

ства следует, что

Для получения численного решения запишем задачу (2.3), (2.4) в перемещениях. В этом случае начальная задача (1.10) примет вид

^ = Ы|А=О = «(0). (2.8)

Краевая задача (1.11), (1.12) может быть записана в следующем виде

с1и\ сРй ¿Риски, ¿х

¿х2 ¿х2 ¿х ¿и а

"1*=о = °> = 7-л--ч—• (2-9)

1-^1 Х=ХГ)Е

¿х

Рассмотрим алгоритм численного решения этой задачи при помощи конечно-разностного метода.

Пусть х0 - недеформпрованные координаты точек стержня, (I = 0,1, ■■■ ,к); х? - координаты тех же точек, то в деформированном состоянии при А = А^, то есть х3 = хгА), это лаграпжевы (вмороженные) координаты; хг - координаты х

Пусть также < = аг(А^) - напряжение в точке с координатой х^и А = А^.

А

< = <11 (хг, Аз), и\ = и(хг, А^), изГ = и(хг, А^),

иг = и(хг, Аз), иГ = и(хг, Аз).

Примем также следующие обозначения: Н - шаг сетки по переменной ж, ДА -шаг интегрирования задачи Коши по переменной А —з - количество точек сетки на ]-м шаге.

На нулевом шаге алгоритма примем

ж 0 — ж^, ж г — ж^, а% — 0, и% — 0,

то есть стержень находится в покое.

Применим метод сеток к задаче (2.9). аппроксимируя первые и вторые производные на трехточечном шаблоне с точностью до 0(Н2) по формулам

и'0

< и'г ■■

и'ы

Поскольку граничная точка совпадает с последней точкой в сетке, то приращения и производные их равны:

иМ0 — uг, и1%0 — uг, и1Ч0 — иг •

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

Используя найденное поле скоростей перемещений по параметру и интегрируя начальную задачу (2.8) методом Эйлера, находим поле перемещений и1. Отметим, что на первом шаге вычислений приходится использовать метод Эйлера с точностью О(ДА), а далее появляется возможность применять более точные методы.

Координата конца стержня определялась экстраполированием аппроксимируемых по методу наименьших квадратов восьми уже найденных эйлеровых координат точек сетки. Рассматривалась линейная аппроксимация.

На г-м шаге алгоритма имеем (нагрузка на правом конце равна аГ): а1 (ж) - поле напряжений, и3 (ж) - поле перемещений,

3

ж\ - координаты вмороженных точек, иГ — перемещение правой граничной точки, жг

Используя найденное на предыдущем шаге значение перемещения правой граничной точки, проверяется, требуется ли добавление еще одной точки сетки. Для этого сравниваются два значения жГ и (N3-1 + 1)Н. Если первое больше второго, то добавляется еще одна точка, значение в которой вычисляется с помощью интерполяции кубическими сплайнами.

Далее, аналогично начальному шагу алгоритма применяется метод сеток к задаче (2.9). Производная в граничной точке аппроксимируется с неравномерным

—-(-3«0 + 4«1 - и2), и" = — (—Зм'-, + 4Ц - и'0),

2Н 2Н

= 7^7 Ы-2 - 4'Ыдг_1 + З'Ыдг), = Г-2 ~ 4ыДГ-1 + Зыдг)-

Рис. 1

шагом по формуле

1

Кг

К + Кг

1г

"¿к

К { КГ

где Кг — расстояние между конечной точкой и последней точкой сетки. Полученная система N3 + 1 уравнений решается методом прогонки. По найденному полю скоростей перемещения и3 находим поле перемещений и\+1 для следующего значения параметра с помощью метода Эйлера-Кошп [10]

= и>-1 +2ДАи/.

иг

алгоритма применяется аппроксимация по методу наименьших квадратов.

По найденному полю перемещений находим координату правой граничной точки. используя формулу

„¿+1 _ „,0

3 + 1

= хг + и

Далее процесс повторяется.

В ходе реализации алгоритма были рассмотрены различные подходы.

Было показано, что при добавлении новой точки в сетку одинаково хороши как линейные, так и нелинейные методы аппроксимации. Влияние выбора метода добавления точки па конечный результат не превосходит 0(К2).

Экстраполирование для нахождения перемещения конца стержня также было реализовано разными методами: методом наименьших квадратов для кривой первого и второго порядков и кубическими сплайнами. Наименьшее отличие от точного решения (2.7) при сильных деформациях (более 3/4 от начальной длины) показал метод наименьших квадратов с кривой первого порядка, построенной по 8-ми предшествующим точкам.

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

2

3

3

и

и

г

г

Рис. 2

Рис. 3

Табл. 1

Возможность уточнения результата уменьшением шага по параметру, шаг сетки стержня 10-2 м. Начальное количество просчитываемых точек па стержне 100. конечное 170

ДА Ошибка (б) Относительное время счета

0.6 5.631 • 10~2 0.5

0.3 3.679 • 10"2 0.7

0.15 2.482 • 10~2 1.2

0.06 1.246 • 10~2 3.0

0.04 6.500 • 10~3 5.0

0.03 3.596 • 10~3 8.5

Увеличение точности можно достигнуть не только уменьшением шага разбиения сетки стержня, но и уменьшением шага интегрирования. Табл. 1 демонстрирует влияние изменения шага интегрирования на среднеквадратичную накопленную ошибку и время, затраченное на вычисления.

Рис. 4

3. Плоская задача

Исследуем напряженно деформированное состояние прямоугольной упругой пластииы. два смежных края которой закреплены по нормальным перемещениям н свободны по тангенциальным. Два других края свободны, причём к правому приложена равномерно распределенная нагрузка интенсивностью а (см. рис. 4). Для этой задачи уравнения равновесия (1.1) и выражения для компонентов ец тензора конечной деформации (1.2). записанные в ортогональной декартовой системе координат Эйлера (х1,х2), выглядят следующим образом:

еи

е22

ди 1 дх1

ди2

дх2

дан да12

дх1 дх2

да 22 ! да и

дх2 дх1

1 ( ¿>г/1 \

~ 2 ^¿>.гл )

1 ( диЛ

2 у дх2 )

+

+ Р^2

(3.1)

+

+

V дх1

дх2 )

(3.2)

1 / ди1 ди2 612 = 2 (¿Ь2 + д^

Закон Гука примем в следующем виде:

Е

а 11

а22

(1 + М)(1 - 2М) Е

(1 + М)(1 - 2М)

а12 = Се12,

ди1ди1 ди2 ди2 дх1 дх2 дх1 дх2

[(1 - ц)еп + ^22]; \рец + (1 - м)е22] :

(3.3)

где ^ - коэффициент Пуассона.

Следуя методу продолжения решения по параметру, полагаем, что все неизвестные, входящие в уравнения, зависят от некоторого параметра А. Продифференцировав по этому параметру уравнения (3.1), (3.2) и (3.3), получаем следующую систему уравнений линейных относительно скоростей изменения неизвестных функций:

да 12

11

дд_ дх1

да22

+

дх2

дай дх2 дх1

+

+ рЛ = 0,

+ рЛ2 = 0,

(3.4)

2

2

2

2

eii

dui\ dui du2 dÚ2

дх1 I дх1 дх1 <9.гл

/ \ dúo дгч дгч

дх2 / 9ж2 9ж2 9ж2 '

(3.5)

ei2

^ дгч А дгч Л \ дг/,2 дгч дгч дг/,2 дг/,2

дх2 ) дх2 \ дх1 ) дх1 дх2 дх1 дх1 дх2

E

& 11

&22 =

(1 + М)(1 - 2М) Е

(1 + м)(1-2м)

[(1 - ц)ёп + ^¿22] , \pén + (1 - H)¿22] ,

(3.6)

&12 = Gé

12-

А

цирования по этому параметру граничные условия примут вид (см. рис. 4.):

х1 =0, U1 = 0; é12 = 0. &(1 + м)(1 - 2м)

х\ = а, (1 - yttjen + ytte'22 = ---; ei2 = 0.

E

х2 =0, U2 = 0; é12 = 0-

х1 = b, péu + (1 - n)¿22 = 0; é12

0-

(3.7)

(3.8)

(3.9) (3.10)

Еслн данную краевую задачу решать в перемещениях, то приходим к линейной

u1 u2

жения по параметру Л определяет правые части начальной задачи (1.10). Найдём аналитическое решение задачи (3.1) (3.3).

Точное решение будем искать в виде: U1 = 0^1, U2 = С2х2, где C1, C2 — некоторые константы. Полагаем, что axx = a, ayy = 0. Из уравнения (3.1) получаем:

ей = Gi - -С2,

€22 — С2 ~ -СЦ.

Добавляя в эту систему уравнение (3.2) и выражая оттуда С11 С2 ? получаем следующую систему квадратных уравнений:

Решая полученную систему и подставляя найденные значения С1 и С2 в вы-и1 и2

«f = (1-,/1-2(1- S)-)X1

„2 = ( 1-А/1 + 2(М + М2)- |Х2

(З.И)

1

где X1 и X2 - начальные лагранжевы координаты точки, а иг и и^ — перемещения соответствующей граничной точки.

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

а,

12 х1 х2

сторону.

Чтобы получить численное решение задачи, нужно преобразовать уравнения

и1 и2 и1 и2 х1 х2

нужно, и подставим полученные выражения в уравнения (3.4). Получим два громоздких уравнения, в которых будут присутствовать частные производные первого и второго порядка:

/1 (иьи2,иьи2) = 0,

(3.12)

/ (иьи2,иьи2) = 0.

Дальнейшее численное решение этого уравнения происходит аналогично одномерной задаче. Будем использовать конечно-разностный метод. Пусть разбиение сетки по осям х^ и х2 будет состоять из N и М элементов соответственно. Про-

и1 и2 и1 и2 х1 х2

уравнениях (3.12) и граничных условиях (3.7) (3.11). В начальном состоянии нам и1 и2

и1 и2

используются во всех внутренних точках сетки, а граничные условия (3.7) (3.11) на соответствующих границах. Заметим, что в каждой точке сетки присутствуют две неизвестные иЦ и иЦ и определены два уравнения. Итак, получена система линейных уравнений с сильно разряженной квадратной матрицей А 2MN-гo порядка, которую схематично можно представить так:

«Л ( 0

" ) = ( i /1 , .л/1 п.,Л (3-13)

где в правых частях уравнении, соответствующих правому граничному условию (3.8), стоит некая функция а, зависящая от параметра Л, а вектор неизвестных состоит из ¿i и «2 во всех точках начальной сетки. Функция а выбирается произвольно. В описанной реализации она берётся в виде а = Л2/2 , тогда а = Л.

Полагая нагрузку на правом краю пластины равной а | Х=1, решаем систему линейных уравнений. Далее используем найденные значения функций ¿i и ¿2 для интегрирования задачи Коши (1.10). Полученные значения ui, «2 являются перемещениями в точках сетки при нагрузке на правом конце пластины, вычисленными

Л

интегрирования краевые точки сетки и пластины совпадали, то очевидно, что после деформирования тело стало больше начальной сетки по оси x1 и меньше по оси x2. В этом случае требуется добавить точки на правой границе. Это можно сделать, например, линейной экстраполяцией:

Хг|Л=1 = -=Ц1| ' (ЗЛ4)

хг1а=о и11 А=0

где Хр - х1-координата правой границы пластины. Эта формула следует из предположения, что начальные лагранжевы координаты линейно деформировались.

(Хр11

^2

Рис. 5

Точки верхней границы деформированного тела можно найти, используя координаты Лагранжа до деформации:

_£

( X2 I

^к -,= Г|л=.°/, . (3-15)

т1а=0 "2|л=0

где ХГ - х2 -координата верхней границы пластины. После вычисления следующего значения функции ¿г|л=2 повторяется процесс построения матрицы Л, решение системы линейных уравнения и интегрирование задачи Коши. При значении Л больше единицы процесс экстраполяции, описанный формулами (3.14). (3.15). изменится следующим образом:

X11 X11 X 21 X 21

У1| _ ЛГ 1л=к-1 ЛГ |л=0 х^2 1 _ ЛГ 1л=к-1 ЛГ |л=0

Лг , , — ~—г,-;-, Аг

"г|л=й -у-II _?/1| ' Г|л=й X21 - по\

ЛГ|л=0 и11л=к-1 ЛГ|л=0 и2|л=к-1

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

и1

в каждой точке конечно-разностной сетки. Этот график соответствует увеличению пластины более чем па 1/2 по оси, к которой приложили нагрузку.

4. Выводы

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

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

Summary

M.S. Agapov, Е.В. Kuznetsov, V.I. Shalashilin. Numerical solution of strong nonlinear deformation problems in Euler's coordinates.

Numerical solution of problem of strong nonlinear deformation using method of solution continuation with respect to a parameter is considered. Euler's coordinates are used. Numerical results demonstrate efficiency of the approach.

Литература

1. Новожилов В.В. Теория упругости. Л.: Судпромгиз, 1958. 370 с.

2. Шалашилии В.И., Кузнецов Е.В. Метод продолжения решения по параметру и пай-лучшая параметризация в прикладной математике и механике. М.: Эдиториал УР-СС, 1999. 222 с.

3. Shalashilin V.I., Kuznetsov Е.В. Parametric Continuation and Optimal Parametrization in Applied Mathematics and Mechanics. Dordrecht/Bostou/Loudou: Kluwer Academic Publishers, 2003. 236 p.

4. Седов JI.H. Механика сплошной среды. Т. 1. М.: Наука, 1970. 492 с.

5. Ильюшин А.А. Механика сплошной среды. М.: Изд-во Моск. ун-та, 1971. 248 с.

6. Гольденбла/т И.И. Нелинейные проблемы теории упругости. М.: НаукаД969. 336 с.

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

8. Данилин А.Н., Шалашилии В.И. О параметризации нелинейных уравнений деформирования твердого тела // Изв. РАН. Механика твердого тела. 2000. Л' 1. С. 82 92.

9. Паймушии В.Н., Шалашилии В. И. Непротиворечивый вариант теории деформаций сплошных сред в квадратичном приближении // Докл. РАН. 2004. Т. 396, Л' 4. С. 492 495.

10. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987. 600 с.

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

Шалашилии Владимир Иванович доктор физико-математических паук, профессор Московского авиационного института.

Кузнецов Евгений Борисович доктор физико-математических паук, профессор Московского авиационного института.

E-mail: kuznetsov Qmail. ru

Агапов Максим Сергеевич аспирант Московского авиационного института.

E-mail: agapov.rn.axim.0gm.ail.com.