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

УДК 539.3

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

МОДЕЛИРОВАНИЕ НЕЛИНЕЙНЫХ ДИНАМИЧЕСКИХ ПРОЦЕССОВ ПЕРЕНОСА В ПОРИСТЫХ СРЕДАХ1

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

В настоящее время особый интерес представляют нелинейные высокоскоростные процессы переноса в пористых системах. Так, в связи с созданием динамических амортизирующих систем, основанных на принципе поглощения энергии за счет фазовых превращений, возникла проблема моделирования нелинейных высокоскоростных течений жидкости в пористых системах. Нелинейные процессы фильтрации в настоящее время практически еще не изучены, обычно применяются лишь эмпирические нелинейные модели, связывающие скорость с градиентом давления в пористой среде. В работах [1-4] предложен новый подход к нелинейным задачам газовой динамики в пористых системах, основанный на методе асимптотических разложений для уравнений в частных производных, заданных в областях с быстроосциллирующи-ми границами [1, 5, 6]. Показано, что в отличие от линейных задач в пористых системах [7, 8], характерных для медленных процессов переноса, описание высокоскоростных течений основано на решении нелинейных связанных локальных и глобальных задач переноса. Для их решения необходима разработка специального нового метода, поскольку существующие методы расчета параметров локальных потоков [9] достаточно приближенны и, как правило, не учитывают таких важных особенностей течения, как периодичность потоков по границам ячеек периодичности, влияние геометрической формы поры и др. Цель настоящей работы — разработка нового метода решения, основанного на модели асимптотических приближений для нелинейных процессов переноса в пористых средах.

Математическая постановка задачи. Рассмотрим пористый материал периодической структуры (рис. 1), поры которого заполнены вязкой сжимаемой совершенной жидкостью (газом).

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

Введем следующие обозначения: 10 — характерный линейный размер ячейки периодичности (ЯП) V; х0 — характерный глобальный размер всей пористой области V; £0 — характерное время;

= х0/Ь0 — характерная скорость; к = 10/х0 ^ 1 — малый параметр; £ = х/к — безразмерные локальные координаты, изменяющееся в пределах ячейки периодичности V; х = х/х0 — глобальные координаты.

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

Рис. 1. Схема периодической структуры пористого материала

dp ^ £+v

pv = 0, x е V;

dpv

~3t

д

9ip{e + 2) "у 2 p av = ¡i(Vx • v)E + ¡12(Vx 0 v + Vx

+ vx • (pv 0 v + pE) = vx • av;

P

+ vx • pv l e + — + - = vx • (av • v);

(1)

0 vJ

^ р = Ярв; е = су в; V =0,

где р, v,p,e,в - плотность, вектор скорости, давление, внутренняя энергия и температура жидкости; а у — тензор вязких напряжений;

— коэффициенты вязкости, которые полагаются "малыми", т.е представимыми в виде ца = к2; су,Я — соответственно теплоемкость при постоянном объеме и газовая постоянная; уж - набла-оператор Гамильтона, который в декартовом базисе имеет следующий вид [10]:уж = д/дх = егд/дхг.

Все функции П = {р,\,9,...}, входящие в систему (1), полагаются квазипериодическими, т.е. они зависят от трех аргументов: П = П (х, £, г) [1], X € V, £ € V "медленно" изменяются относительно аргумента X и являются однопериодическими по второму аргументу £: П ^Х, — ^, = П^х, £1 , £2, ^, и т.д. для всех

Условия периодичности далее обозначаются как [[П]] = 0. Дифференцирование квазипериодических функций выполняется в соответствии с правилом дифференцирования сложной функции:

ужп (х, £, г) ^ Ужп (х, £, г) +1 vп (х, £, г), (2)

к

где v, v^ — операторы Гамильтона по координатам х и £ соответственно.

Введем оператор среднего значения (•) по ЯП V для квазипериодических по локальным координатам £ € V функций П [1]:

(П) = ш / ^, (3)

1 г ~

где % = —— ¿/^ — объемная доля жидкости в ЯП (пористость), а

№ i 3

Vi — объем ЯП V.

В соответствии с методом асимптотического осреднения ([4]) функции системы уравнений (1) представляются асимптотическими рядами по малому параметру к:

П (х, £, г) = П(0) (х, £, г) + кП(1) (х, £, г) + к2п(2) (х, £, г) +... . (4)

Используя правила дифференцирования сложной функции (3), путем подстановки выражений (4) в систему уравнений (1) и приведения слагаемых при одинаковых степенях малого параметра к получаем последовательность так называемых локальных задач газовой динамики на ЯП V. Собирая члены при степени к-1, получаем локальную задачу газовой динамики нулевого уровня на ячейке периодичности:

f v • p(0)v(0) = 0;

v • (p(0)v(0) 0 v(0) + p(0)E) = 0;

v • (V0) (cvв(0) + + P(0)) v(0) = 0;

v(0)2' 2

p(0) = Rp(0)e(0), $ е Viff; v(0)-n(0) = 0, $ е £isg; [[П]] = 0;

<р(0М0)) = рv, <р(0)) = р, <в(0)) = 0,

(5)

в которой неизвестными функциями являются плотность р|°, вектор скорости V0 и температура в^0) в нулевом приближении. Задача (5) похожа на исходную (1), однако отличается от нее тем, что: 1) решение ищется только на ячейке периодичности, поэтому кроме условий на твердой поверхности в задачу (5) входят условия периодичности;

2) вследствие малости вязкости задача (5) не содержит вязких напряжений, т.е. она представляет собой систему уравнений Эйлера;

3) задача (5) относится к установившимся процессам, так как не содержит производных по 1 Кроме того, к задаче (5) присоединяются дополнительные условия нормировки — условия на средние значения функций, где величины р(хV (х,1) и 0 (х, ¿) представляют собой средние по ЯП плотность, вектор скорости и темепература жидкости. Они являются внешними данными задачи (5), т.е. полагаются заданными и зависят только от глобальных координат х и 1 Вследствие условий нормировки задача (5) относится к интегродифференциаль-ному типу. Это обстоятельство, а также наличие условий периодичности затрудняет применение большинства широко распространенных численных методов решения локальной задачи газовой динамики (5). Ее решение к настоящему времени удалось получить только для очень простой геометрической формы поры — цилиндра [2].

Далее рассмотрен новый метод нахождения численно-аналитического решения задачи (5) для ЯП с криволинейной геометрией пор.

Система уравнений Эйлера (5) допускает два первых интеграла — адиабату Пуассона и интеграл Бернулли вдоль линии тока. Следовательно, в нулевом приближении локальный процесс переноса в поре является адиабатическим, а задача (5) путем стандартных преобразований [11] может быть представлена в эквивалентом виде:

f • p(0)v(0) = 0; v(0) x (V¿ x v(0)) = V?i*: v2 v2

"2 + ^ = "2 + Cp0 * = i*

(0)

(0)

p(0) _ /P^V 0(o) _ p* V P* J ' 0* v(0) • n(0) = 0, £ G £iSg;

(0)

Y-1

(6)

p(0)v(0)) = Pv

p(0)> = p-,

р

[[П]] = 0; (в(0)> = в.

Первое уравнение системы (6) — уравнение неразрывности, второе — векторное уравнение установившегося движения газа в форме Громеки-Лемба, третье — интеграл Бернулли, четвертое и пятое — адиабата Пуассона и уравнение баротропии (оба являются следствием уравнения энергии и уравнения состояния). В системе (5) независимы только шесть скалярных уравнений. Здесь введены обозначения: р*, р *, в *, V *, г* — постоянные интегрирования, входящие в число неизвестных (независимы из них только три), и зависящие от линии тока; V = |у(0)| — модуль вектора скорости. Для постоянных интегрировния справедливы соотношения

Y-1

* * л --

p = p Cp0 Y .

Y = —, R = Cp — cv. cv

(7)

Осесимметричная локальная задача. В ЯП V кроме указанных выше локальных декартовых координат £г введем цилиндрические локальные координаты r, z, связанные с £г соотношениями [10]

Í1 = r cos £2 = r sin £3 = z (er, e^, ez — физический базис цилиндрической системы координат). Положим далее, что пора в ЯП имеет осесиммметричную форму с осью симметрии Oz и обозначим r = /s(z) функцию формы поверхности контакта жидкости с твердым телом в ЯП. В цилиндрической системе координат формула (3) примет вид

/s 1/2

(П) = / / П (r, z) rdrdz. (8)

0 -1/2

Положим, что входные данные задачи (5) согласованы с осесимме-тричной одноканальной структурой пор, т.е. являются одномерными и соответствуют течению жидкости в направлении оси Oz: р(Х3, t), V=vz (x3,t) ez и 0(x3,t). В этом случае решение задачи (6) также будет обладать осевой симметрией: р(0) (r, z), v(0) = ví0) (r, z) er + + vZ0) (r, z) ez и 0(0) (r, z), т.е. будет зависеть от двух координат: r и z.

Введем безразмерные неизвестные функции в задаче (6):

р(0)

р = —; р

0(0)

vr =

vf

— ) a

vz =

vZ0)

где а = \J~yR6 — осредненная скорость звука жидкости. Тогда задача (6) в цилиндрической системе координат принимает вид

' д {руг) + д (рУг) + ря = 0_ дг дг г '

/ dvr

\dZ

dvz dr

vz =

dvr dz

dvz

dr

dr ' dl*

vr =

dz'

^ v2 v2 ~

- + ъв = 2 + 7i0* =i*;

p рв в f р4 7-1

(10)

p*

р*

р*

(pvr. > = 0, ^ > = M, (р> = 1;

<^0 ^ = 1; (vr nr + vz n

z nz )|S

= 0,

где 7i = 7/(7 — 1); v2 = V2 + vf; MM = Vz/а — усредненное число Маха. Интегральное условие (pvr) = 0 задачи (10) связано с одноканальной структурой макропор, так как Vr = 0.

Задача (10) содержит только один внешний параметр — безразмерное число Маха M и одну константу j1, следовательно, ее решение можно представить в следующем обобщенном виде:

= (e,Yi ,m) ,

(11)

где Цд = |р,уг,уг,91. Предположим теперь, что получено решение

задачи (10) для нескольких значений чисел Маха М: М0 < М < < ... < Ма < ... < Мп, т.е. определен п +1 набор функций вида (11), а именно = (£г,/у1,М0). Тогда для произвольного значения Ма, а = 0, п, можно использовать сплайн-интерполяцию, например, третьей степени

Qß (f,7i,M) = £ (f,7i) Maш,

ш=0

где

^ (?,Ъ) = (?,Ъ) , еСЛИ Ma-i <M <Ma, ш = 07З, ß = М, а = 1"П.

(12)

(13)

*

При фиксированных значениях и y1 коэффициенты Q^a (£\ Y1) находятся стандартным способом вычисления сплайн-функций. Такой метод решения локальной задачи предполагает, что задача (10) решается n раз, затем в памяти ЭВМ формируются и хранятся 16n двумерных массивов 0,рша (£\ Y1) для каждого значения параметра y1 .

Отметим, что для нахождения осредненной скорости V (x, t) жидкости в методе формулируется специальная осредненная задача [2, 3]. В силу нелинейности локальной (10) и осредненной задач, они, вообще говоря, связаны и формально должны быть решены совместно. В предложенном методе оказывается возможным "развязывание" указанных задач за счет создания банков данных с решением локальной задачи (10), в которых накапливаются массивы Qg^ для раз-

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

Метод решения локальной задачи. Для решения задачи (10) применим следующий приближенный метод, который основан на следующих допущениях:

1. Все трубки тока в поре в ЯП пропорциональны поверхности контакта поры с твердой фазой r = /s(z), т.е. имеют вид r = /r(z) = /s(z)r/A0, где A0 — радиус поры при z = -1/2.

2. Модуль вектора скорости v зависит только от осевой координаты z: v = v(z) (здесь и далее, для простоты, символ "~"опущен).

3. Вместо уравнения неразрывности рассматриваем интегральное уравнение закона сохранения массы для произвольной подобласти

V с V :

J V • pvdV = 0.

Компоненты вектора скорости имеют вид

vr = v(z) sin <^(r, z); vz = v(z) cos <^(r, z). (14)

Здесь <^(r, z) — угол между касательной к линии тока и осью Oz, который в силу допущения 1 является известной величиной и определяется только функцией /r (z) по формулам [10]:

1 /' (z)

sin w(r,z) = —. -; cos w(r,z) = —. r -. (15)

, ) Vl + (/r(z))2' ^ , ) Vl + (/r(z))2

Модуль вектора скорости v = v(z) является неизвестной функцией и подлежит определению. Выбирая в качестве V часть всей ЯП, ограниченную плоскостями z = const, с учетом граничных условий на твердой стенке и допущений 1 и 3 получаем, как и в стандартной одномерной теории, условие постоянства скоростного напора жидкости

в каждом сечении ЯП:

/s

J p(z)vz(r,z)rdr = Q = const, vz. (16)

0

Проинтегрируем уравнение (16) по z от -1/2 до 1/2 и умножим получившееся выражение на 1/| Vg i, тогда с учетом интегрального условия (pVz) = MM системы (10) получим

1/2 /s

2п f f Q

MM = (pVz) = 1J p (z) V (z) J cos <^(r, z)rdrdz = -щ. (17)

-1/2 0 Отсюда находим, что Q = M | Vg |.

Пусть L — линия тока в ЯП V (см. рис. 2), начинающаяся от плоскости z = -1/2. Выберем два сечения ЯП, для которых z = -1/2 и z = const. Выберем в качестве значений p*,p*, e*,v*, г* газодинамические параметры, соответсвующие плоскости z = -1/2, тогда вдоль линии тока L из интеграла Бернулли, уравнения баротропии и адиабаты Пуассона системы (10) получим следующие соотношения:

в V2 - V2

w _ umax w . (18)

9* v2 — V2

w "max

p /v2 - V^1/(7-1)

P_ = Vfx-_ ; (19)

P* V^ax - V*V

p ( p \Y (V2 _ V2 \y/(y-1)

4 = 4 = Vfx-* . (20)

P* VPV V^ax - V2/ Здесь vmax обозначено максимальное значение модуля вектора скорости, которое достигается на линии тока

Vmax = л/2Т1в* + V2. (21)

Из формул (18)-(20) следует, что температура, плотность и давление в поре не зависят от координаты r, они зависят только от z, поскольку определяются модулем вектора скорости. Из формул (18)-(20) также получаем

в = e*G (v); p = p*H (v);

(22)

p = p*F(v); p* = pF (v),

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

F (v) =

2 2\ 1/(7-1)

V2 — v2 4

max *

чУтах - (23)

Ё (у) = 1/Ё(у); С (у) = [ё (у)]1-7 ; Н (у) = [Ё (у)р .

Выражая утах через 9* по формуле (21), приходим к следующему представлению функции Ё (у):

2719* 4 1/(7-11

Ё (у) = Ё (у,9«,у*)=^719* +'у2 - ■ (24)

Применяя оператор осреднения (8) к первому и третьему выражениям в (22), с учетом интегральных условий из (10) для р, 9 получаем, что должны выполняться следующие соотношения:

9*<С (У)> = 1; (25)

Р*(ё(у)> = 1,

поскольку 9* и р* не зависят от г. Кроме того, имеет место следующее соотношение, вытекающее из (16) и (22):

5(г) = М\Ъд\Ё (у) /(р*у), (26)

М*)

где Б (г) = 2-к J еоэ р(т, г )тйт. Таким образом, имеем систему урав-

0

нений (25) и (26) относительно трех неизвестных констант (р*, 9*, у*) и функции у (г). Исключая из этой системы р* с помощью второго соотношения (25), приходим к следующей системе:

' 1 - 9*<в(у)> =0;

1 {д \ <ё(у)> ; (27)

УБ

М\У,д \Ё (У) - = 0.

(Ё(У))

Здесь второе уравнение получено из третьего для г = -1/2, при котором Ё (у) = Ё (у*) = 1 согласно (23). После определения величин 9*, у* и у плотность р* вычисляем по второй из формул, а р* по формуле (7). Плотность р, температура 9 и давление р определяем по формулам (22).

Анализ решения локальной задачи. В зависимости от геометрии поровой области У^д, величины осредненного числа Маха ММ и коэффициента Пуассона 7 можно получать различные решения локальной

Рис. 3. Решения для дозвукового (а) и сверхзвукового (б) течения

задачи (10) — для дозвукового и сверхзвукового течения. Для рассмотрения этих решений рассмотрим третье уравнение в (27), из которого следует, что

5 (V) = М |%| ^ (V)/Д. (28)

Функция зависимости площади сечения Б поровой области У^д от модуля скорости V имеет локальный минимум в точке V = Vкр = = [(7 — 1)/(7 + 1)]1/2^тах (левее и правее которой ветви кривой направлены вверх, пересечений с осью абсцисс нет) и две асимптоты: V = 0

и V = Vmax.

Если в поре реализуется дозвуковой режим течения, то скорости частиц жидкости лежат на левой ветви кривой Б (V), газодинамические параметры — периодические функции локальной координаты г (рис.3,а). Если же в поре течение становится сверхзвуковым, то периодического решения локальной задачи (10) не существует. Иначе говоря, трансзвуковое решение является переходным. Наличие перехода определяется значением массового расхода 0 и функцией Б (г).

Численный метод решения локальной задачи. Для численного решения системы уравнений (27), фактически представляющей собой нелинейное интегральное уравнение относительно функции v(z), введем на отрезке [—1/2; 1/2] оси О г сетку узлов гт = —1/2 + + т/п,т = 0,..., п, п — целое число. Тогда разностная аппроксимация оператора осреднения (8) для случая функции О(г), зависящей только от г, будет иметь следующий вид:

(О) « ^^ ^ [От-1 (гт-1) + От(гт)] , (29)

21% i m=

m=1

ГДе Н — Zm zm—1.

9

Подставляя (29) в (27) и записывая эти уравнения в узлах сетки, получаем следующую нелинейную систему алгебраических уравнений относительно неизвестных значений функций в узлах сетки П = (п1,...,Пп+2 )т = (9*,у*,у1,...,ут,...,уп)т:

{ ф1 (т,...,Пп+2) = 0; Ф (п) = 0 ^ I ......................................................................(30)

[ фп+2 (п1, ...,Пп+2) = 0,

где ф = 1 - 9* (С (у)>;

- У3-2Б (г,3_2) _

ф = М\У^д\ Ё (у3-2)--/ \ , 3 = 2,п + 2, если Уо = у*.

(Р(у))

Система (30) решается численно методом Ньютона:

n(s+1) = n(s)

lp 4p

{K)(s) (31)

1

,(s)

где р = 1, п + 2 — номер неизвестной в координатном столбце п; в — номер итерации; I — индекс суммирования (номер строки в матрице Якоби); Ар — элементы матрицы Якоби, из которых ненулевые

= - (С (у)>д; А] = МШЁ' (У]-2) ^ - ^

д ] 7- 1 719* /ё(у)

3 = 2,n + 2;

M\Vig |

v2 v2

(32)

A = Fy (vj-2) -^Tjf, 3 = 3,n + 2;

Y - 1 2yi(в*)

A2 = M1^FY (vj-2)-^,3 =3ТПГ2.

М \У(, \

1 - 7 ~ 2/ 719*:

В методе Ньютона при вычислении элементов Ар матрицы Якоби на з-й итерации введена модификация, заключающаяся в том, что средние значения (С (у)> и (ё(у)> получаются путем подстановки в формулу (29) для П = {С (у) , Ё (у)} сеточных значений модуля скорости у частиц жидкости, взятых с итерации в - 1.

Пример численной реализации метода. В качестве примера численной реализации метода была рассмотрена пористая среда, поры которой заполнены воздухом (ср = 1,006 кДж/(кг-К), 7 = 1,2, Я = 0,168 кДж/(кг-К)). Характерные геометрические и газодинамические параметры были выбраны следующими: х0 = 1 м, ¿0 = 1 с, у0 = х0Д0, плотность р0 = 10 кг/м3, температура 90 = 293,15 К, давление р0 = р0у0. Функция формы (см. рис.2) поровой обла-

Рис. 4. Зависимость безразмерного модуля скорости V от осевой координаты поры г для различных чисел М и значений геометрических параметров А0, В0

Рис. 5. Зависимость безразмерной плотности р и температуры в от осевой координаты поры г для различных чисел М и значений геометрических параметров Ао = 0,45, Во = 0,15

сти У^д была выбрана в виде fs (z) =

Ao + Bo Ao - Bo

cos (2nz),

2 2

где z E [-1/2; 1/2]; A0 = const, B0 = const, для которых 0 < B0 ^ ^ A0 ^ 1/2.

На рис. 4 и 5 приведены результаты решения локальной задачи (10) для различных значений осредненных чисел Маха M E (0; Mmax], где Mmax — верхняя грань допустимого множества (0; Mmax] чисел Маха M, при которых периодическое решение существует, если геометрическая форма поровой области фиксирована. Так, для A0 = 0,45 и Bo = 0,15 Mmax = 0,10788, для Ао = 0,45 и Bo = 0,3 Mmax = 0,33489.

На рис. 6 на плоскости значений (А0, B0) геометрических параметров поры показаны области существования периодического решения

Рис. 6. Области существования периодического решения (показаны темным цветом) при различных значениях числа М и геометрических параметров

(Ао,Во):

а-д — М равно соответственно 0; 0,2; 0,3; 0,4; 0,5

задачи (10) при y1 = 6 и осредненных числах Маха M E {0,1; 0,2; 0,3; 0,4; 0, 5}. Из этих рисунков видно, что уменьшение площади критического сечения поровой области V^g при увеличении числа Маха M приводит к уменьшению (исчезновению) периодического решения локальной задачи (10).

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

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

1. Димитриенко Ю. И. Механика композиционных материалов при высоких температурах. - М.: Машиностроение, 1997. - 368 с.

2. Dimitrienko Yu. I. Dynamic Transport Phenomena in Porous Polymer Materials Under Impulse Thermal Effects. - Transport in Porous Media. - V. 35.

- 1999. - P. 299-326.

3.Димитриенко Ю. И., Иванов М. Ю. Разработка численного метода решения локальной задачи нелинейной фильтрации в периодических пористых средах. - В c6.: Современные естественно-научные и гуманитарные проблемы.

- М.: Логос, 2005. - С. 469-478.

4. Д и м и т р и е н к о Ю. И., И в а н о в М. Ю. Разработка метода асимптотического осреднения для решения нелинейных задач фильтрации в периодических пористых средах. - В c6.: Математика в современном мире / Под ред. Ю.А. Дробышева. - Калуга.: Изд-во КГПУ - 2004. - С. 155-163.

5. Санчес-Паленсия Э. Неоднородные среды и теория колебаний: Пер. с англ. - М.: Мир, 1984. - 472 с.

6. Бахвалов Н. С., Панасенко Г. П. Осреднение процессов в периодических средах. - М.: Наука. - 1984.

7. Димитриенко Ю. И., Глазиков М. Л. Моделирование процессов фильтрации в периодических пористых средах. - Вестник МГТУ им. Н.Э. Баумана. Сер. "Естественные науки". - № 1. - 2003. - C. 59-71.

8. Д и м и т р и е н к о Ю. И., Г л а з и к о в М. Л. Разработка метода асимптотического осреднения для решения задач газовой динамики в пористых средах.

- В сб.: Математика в современном мире / Под ред. Ю.А. Дробышева. - Калуга.: Изд-во КГПУ. - 2004. - С. 163-177.

9. Нигматулин Р. И. Динамика многофазных сред. Ч. I. - М.: Наука, 1987.

- 464 с.

10. Д и м и т р и е н к о Ю. И. Тензорное исчисление. - М.: Наука. - 2001. - 575 с.

11. С е д о в Л. И. Механика сплошной среды. T. 2. - М.: Наука. - 1976. - 552 с.

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