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

НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА

НАУКА и ОБРАЗОВАНИЕ

Эл № ФС77 • 48211. Государственная регистрация №0421200025. ISSN 1994-0408

электронный научно-технический журнал

Моделирование высокоскоростных процессов в демпфирующих

системах с фазовыми превращениями

# 11, ноябрь 2012

Б01: 10.7463/1112.0493560

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

УДК 539.3

Россия, МГТУ им. Н.Э. Баумана IvanovM-fn11@yandex.ru dimit.bmstu@gmail.com

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

Физическая модель гетерогенного рабочего тела (ГРТ). Пористое тело имеет два масштаба пор: макропоры и нанопоры. В ходе проведенных исследований установлено, что поглощение механической энергии обусловлено проникновением жидкости в нанопоры, если давление в ГРТ достигло некоторого критического значения [3]. Такое затекание жидкости рассматривается как фазовое превращение, которое сопровождается образованием новой фазы (§-фазы). В рамках рассматриваемой модели считается, что §-фаза заполняет макропоры и нанопоры ГРТ (рис. 1).

Рис. 1. Схема течения §-фазы в ГРТ.

При моделировании рассматривалось ГРТ с периодической структурой макропор. Периодичность означает, что в гетерогенной системе можно выделить характерный элемент (ЯП), повторяющийся с некоторым периодом. Кроме того, в расчетной модели рабочего тела амортизатора выбрана одноканальная геометрия макропор (рис. 2).

Рис. 2. Расчетная модель микроструктуры ГРТ с одноканальной пористостью (показан разрез

ЯП плоскостью симметрии).

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

дрд „ дрд V „ _

^ +УХ рVg = 0, + УХ- ряví

дл

дл

IV„ = Vх аг

др (\+ '2/2} + VхРдVд (ед + vg/2) = Vя ■(аV,),

^ = ^ 0 + СК (^ -дд 0 ) , Рд = ЯРдвд , О д = - РдЕ + ^ Т д

Т д = Лд EV х • V д + кд (V, ^ V д + ^, 8 V д ^ ) , X ^д '

(1)

д2

Р* дцГ ^ ха*, £ = П--а,, 2£ ^ ® и* ® и*) , х ,

где рд - плотность, Vд, - вектор и модуль скорости, ад - тензор напряжений,

ед - внутренняя энергия, ск - теплоемкость при постоянном объеме, рд - давление,

Я - газовая постоянная (коэффициент сжимаемости, если в качестве §-фазы рассматривается жидкость), вд - температура, тд - тензор вязких напряжений, Лд, кд - параметры Ламе

§-фазы в области ; р5 - плотность, о5 - тензор напряжений, и5 - вектор перемещений, 8 - тензор деформаций, П - тензор упругих податливостей твердой фазы в области у; E - метрический тензор; Vх - оператор Гамильтона по декартовым координатам х]; t - время.

Граница раздела фаз является подвижной и на ней заданы следующие

соотношения: [р]Б + ргvg • п = 0, -р5vgD + [о]• п = 0, где [ • ]-оператор скачка функции

при переходе через поверхность Е, для которого [О] = О^" -ОЕ ; Б - скорость движения

поверхности Е^ (скорость фазового превращения); п - вектор нормали к поверхности Е^ .

Система (1) дополняется граничными условиями на внешней поверхности ГРТ и начальными условиями.

Математическое моделирование физических процессов в демпфирующей системе осуществляется с привлечением теории гетерогенных пористых систем с фазовыми превращениями [2, 4], в основе которой лежит метод асимптотического осреднения процессов в периодических структурах [5, 6]. Для ГРТ можно ввести малый параметр к = 1/Ь 1, где I - линейный размер макропор, Ь - линейный размер всей демпфирующей

системы, газодинамические и механические параметры О = {р^?,v3,...} которой выражаются

асимптотическими рядами относительно к: О = О(0) + к- О(1) + 0 (к2), где

О (г, х1, t)- периодические функции по локальным координатам %1 = х1 / к,

х1 - безразмерные декартовы (глобальные) координаты.

Локальные задачи на ЯП. С целью определения микропараметров нулевого

о2 (г, х], t) и первого О^ (г, х], t) порядков фаз ГРТ формулируются так называемые

локальные задачи газовой динамики и механики твердого тела на ЯП [1, 2].

Так, локальная задача газовой динамики нулевого уровня, описывающая микродвижения §-фазы в макропоре, имеет следующий вид:

Vг р^ = 0, Vг {р^ ® vg0) + р^Е) = 0,

Vг(рgO^{cg0gO^+ ^°')2/2) + рР)vg0) = 0, р(0)= К^ер, % *Уг; (2)

vg0)•n(0) = 0, % ^; |рр)| = | у?) = ео)

= 0;

(x,t), (р^0)vg0Л =Pg (Х^) Vg (х,0, т=в, (х,0

/„ „

где Vг - оператор Гамильтона по локальным координатам %, cg - теплоемкость §-фазы при постоянном давлении, У г, 2 г - части области Уg и поверхности 2^ в ЯП соответственно, [О] = 0 - условие периодичности функции О по локальным координатам % на границах ЯП, ( • - оператор среднего значения по /-ой фазе от функции О, для которого О = (О). = (1/|Уг |) ^ ШУг , Уг | - объем /-ой фазы в ЯП, / = 5, g.

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

соотношения на искомые функции (средние значения микропараметров в ЯП), б) в систему уравнений (2) включены условия периодичности на искомые функции. В работе [7] предложен приближенно-аналитический метод решения данной задачи, позволяющий вычислять параметры локального потока в отдельной макропоре в зависимости от скорости движения осредненного потока (макропотока) и геометрической формы макропор. Установлено существование периодических дозвуковых и сверхзвуковых режимов движения локального потока. Показано, что периодического трансзвукового режима течения для локальной задачи (2) не существует.

Метод предполагает создание банка данных решений задачи (2) для различных вариантов геометрии макропоры, физических свойств g-фазы и осредненного числа Маха

M = !а, где vg0 - осевая компонента вектора скорости v(g0) во введенной на ЯП

локальной цилиндрической системы координат (ЛЦСК), а - осредненная скорость звука,

а = ^R6g . В частных случаях, если форма поверхности Zsg в процессе фазовых

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

Qg0) = {pf, vg0>, 6f} задачи (2) из анализа сформированного банка данных решений

локальной задачи относительно среднего значения Qg и числа Маха M по ЯП:

Q g°) (z, M, Qg ) = Qg ([q^ ,0 (M ) - Qg% д (M)] exp {- z2/ ( 2<Q )} + Qg°i д (M )), (3)

где z - осевая координата ЛЦСК; Qg^o (M), Qg°m x (M)- интерполяционные полиномы m-ой степени, отражающие поведение функции Qg0) в начале ЛЦСК и на границе ЯП

соответственно при различных значениях осредненного числа Маха M; < - параметр, выбираемый из условия наилучшего приближения аппроксимационного решения к численному путем минимизации максимальной относительной погрешности аппроксимации вдоль оси симметрии макропоры при фиксированной геометрии и различных числах Маха M, т.е. <Q зависит только от геометрической формы области течения и свойств g-фазы. Установлено, что функция вида (3) дает наилучшее приближение для сверхзвуковых периодических решений с относительной погрешностью не превышающей 3 %.

При моделировании функция формы f области V^g макропоры была выбрана в

следующем виде: /z(z) = (A0 + B0-(A0 -B0)cos(2nz))/2 , где z e[-l/2;l/2]; A0 = const,

B0 = const, для которых 0 < B0 < A0 < 1/2. В качестве g-фазы рассматривался воздух. Рис. 3 и 4 отображают характер поведения плотности и температуры в макропоре при различных значениях (дозвуковых и сверхзвуковых) осредненного числа Маха M. На рис. 5 и 6 показаны графики изменения радиальной и осевой компонент скорости в ЛЦСК для макропоры вдоль одиннадцати линий тока при фиксированном значении M = 5.55. Каждая линия тока пропорциональна функции формы /. На приведенных графиках геометрия

области V4g определяется следующими параметрами: A0 = 0.45, B0 = 0.1.

¡•••••••а НИГГ»»**« 1* **••*••)

°990< V 0 97

\ 0 95

0 91

0 89 • Плотн -Плотн ость (IV = ость (IV = 0.01) 0.05) -

-0гв7- х Температура (II = 0.01) о Температура (IV! = 0.05)

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2

Безразмерная осевая локальная координата

0.3

0.4

0.5

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

ж Плотность (IV = 5.55) о Плотность (IV = 7.55) а Плотность (II = 8.55) д Температура (II = 5.55) х Температура (II = 7.55) - Температура (II = 8.55) —

ж2 ж £ 8 ¡2Х 2ж 2 ж 2

ж & ж £ ж Ш ж £ ж

ъ 1 £ £

м ж к ш Я 1 Ж X к ■

........ [ -0- 1 *ш

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4

Безразмерная осевая локальная координата

0.5

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

I

ф

т „ го -0

I

(О

Безразмерная осевая локальная координата

5

Рис. 5. Графики изменения радиальной компоненты вектора скорости сверхзвукового локального потока для различных линий тока §-фазы в макропоре.

ЧИН*! \+о-< -и* ' Xх V*. уХХ ** « яунг •%-оу/

\\+< • жж ; ■о"- *. О - \+ " - X V хХхх, :ж х> Жжж»„ ЖЖ} 'ЖЖЖЖЖЖЖ) \\+ !|1ххххх> ¡ЖЖЖЖЖЖЖ) уХ' сххх ^ ;жж* жж • -:ж •*_ .• - о- ,• - о ' - п +/ ■0+1

\+ о \\ + ° \\+ °° \\ + \\ + ...... 1....... ж + // ......... > ...... о +/, ■ О + /У о° + // 1 + // + / /

\\ + °ооск ++++++н ' + / \ + >ооо° ■++++++ + / /

■ Линия № 0 о Линия № 1 * Линия № 2 х Линия № 3 ж Линия № 4 • Линия № 5 - Линия № 6 о Линия № 7 + Линия № 8 —Линия № 9 —Линия № 10

-3-

го

о о ср

о

^

0 ф

^

1

ф

т го

I

ОЭ ф

0

1

ср ф

го ср

(О ф

ш

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3

Безразмерная осевая локальная координата

0.4

0.5

Рис. 6. Графики изменения осевой компоненты вектора скорости сверхзвукового локального потока для различных линий тока §-фазы в макропоре.

Микронапряжения с!0 и перемещения и*1 твердой фазы в области У^ (часть области

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

V, • с(0) = 0, Е(0) = П•С0);

2е(0) = V* ® ц(°) + (V* ® и^ + V, ® и!1) + (V ® и!1))7 , % е ; с )-п(0 , % ; | и!1) | = 0; (и«) ^ = 0,

(4)

где Vх - оператор Гамильтона по глобальным координатам. Асимптотический ряд для

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

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

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

Глобальная задача. Осредненные по соответствующим фазам ГРТ функции (макропараметры): плотность р, , скорость V,, температура в§-фазы; тензор напряжений

С, скорость V, вектор перемещений и твердой фазы являются исходными данными для

локальных задач (2), (4) и зависят только от глобальных координат и времени г. Для вычисления макропараметров получена так называемая глобальная задача [1, 2], в которую входят осредненные уравнения фаз ГРТ. Система глобальных уравнений в декартовой системе координат Охк, определяющей решение для всего объема демпфирующей системы,

имеет вид: ^ + [С] = {Р}, в которой введены обозначения:

{А} = р р,рг, ^ ^ ад^ Р,РУ2, Р,рЛ3, р,ре)Т;

{д}; = (0, р,р,у;, р;1, рЯ2, р,а?, &1, &2, &3, Ь;)Т;

{Р} = (-J|р, 3, -РУ(рр), -Р2/(рр), -Р3/(рр), Р1, Р2, Р3, Е,)Т ;

=^+22

а=1

(а +Р рГК )'

г^е

.(0и/з(0)\ •

+ р Ав Р

0™ = Ре РУ + Р^ + <Р°ЧАу;

Ь; = Ре

(ре + р )Г + /р(0)АугАуЛ у,

Я Я Я Я / Я

- Е;

Е, = Р0)Ау;

Ч + 222(Ауа)2] >¡04) ,

, У^ - декартовы компоненты скорости твердой фазы и §-фазы соответственно; Ау, - компоненты вектора пульсаций скорости §-фазы относительно среднего значения V,; Е - компоненты вектора энергии относительного движения §-фазы; Спт - компоненты матрицы [С], для которых Спт ^ Спт8пт, если п, т = 1,9, причем С11 = 0,

С33 = С44 = С55 = - 1 (РР, ) , С22 = С66 = С77 = С88 = С99 = 1 ; Ре - пористость; Р, = 1 \

3 - интенсивность объемного фазового превращения; Рк - компоненты вектора сил межфазового взаимодействия; Е- интенсивность фазового энергетического обмена.

Макропараметры ГРТ вычислялись с применением модифицированной конечно-разностной схемы Мак-Кормака, состоящей из шагов 1) «предиктор», 2) «корректор», 3) «учет правой части» и 4) «сглаживание». Модификация общеизвестной схемы Мак-Кормака [9, 10] заключается в наличии 4-го этапа «сглаживание», который предусматривает использование оператора специального вида, предназначенного для устранения нефизических осцилляций решения, возникающих при прохождении шагов 1)-3) разностной схемы.

Выполнены расчеты для одномерной демпфирующей системы (пористой пластины) с цилиндрическими макропорами. Рассмотрен случай распространения импульса давления, возникшего в результате действия импульсной динамической нагрузки на внешнюю сторону пластины, у которой тыльная поверхность соприкасается с неподвижной опорой. На рис. 7 и 9 показаны волновые картины сжимающих напряжений, возникающие в твердой фазе ГРТ в процессе его деформирования в случае отсутствия и появления фазового превращения. На рис. 8 и 10 показаны графики заданного внешнего воздействия на демпфирующую систему и зависимость модуля сжимающих напряжений на ее тыльной поверхности от безразмерного времени в случае отсутствия и появления фазового превращения.

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

-0.02

к

I

ф

го

I

ф .0 I

о. ф

-0.04

-0.06

го

.

ф ш

-0.08

-0.1

, 0.04

Безразмерное время о 0.025 д 0.05 □ х 0.1 ■ 0.125 а ж 0.175 +0.2 - - -

0.25

0.275

06 ^Хл* к / £ * 4 1*2

Ч/\ / АЛ/ % # ж X У /* / * ж * у • ж * * ж * ж

^ ж ь. Ж » ¡Ик^^наиМЯ ж ^^ ж ^

ж ж ж ж ж ж ж ж ж

Безразмерная глобальная координата

Рис. 7. Волны напряжения в твердой фазе ГРТ (вдоль пористой пластины) без фазового превращения для различных значений безразмерного времени.

0.1

0.08

£ 0.06

0.04

го р

ф ш

0.02

-•-Внешнее воздействие на демпфирующую систему -Напряжения на тыльной поверхности пластины

/ \

\

1 \

0.05

0.1 0.15

Безразмерное время

0.2

0.25

0.3

0

0

Рис. 8. Зависимость внешнего воздействия и модуля сжимающих напряжений на тыльной поверхности пористой пластины от безразмерного времени (без фазового превращения).

Безразмерная глобальная координата

Рис. 9. Волны напряжения в твердой фазе ГРТ (вдоль пористой пластины) при наличии фазового превращения для различных значений безразмерного времени.

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

Безразмерное время

Рис. 10. Зависимость внешнего воздействия и модуля сжимающих напряжений на тыльной поверхности пористой пластины от безразмерного времени (при наличии фазового

превращения).

Если же в демпфирующей системе сложились условия для возникновения фазового превращения, то наблюдается значительное снижение напряжений на тыльной стороне пластины (рис. 10) за счет заполнения нанометровых пор твердого тела g-фазой.

Выводы. Разработана математическая модель, описывающая высокоскоростные нелинейные нестационарные процессы в гетерогенных деформируемых системах с фазовыми превращениями, показаны ее особенности. Модель может быть использована в процессе проектирования перспективных гетерогенных демпфирующих систем. На основе разработанных численных и численно-аналитических методов выполнены расчеты микро- и макропараметров такой динамической системы. Выявлены особенности локальных процессов переноса в макропорах ГРТ с криволинейной границей межфазной поверхности. Результаты расчетов макропараметров показали качественно правильное их соответствие имеющимся экспериментальным данным.

Список литературы

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

2. Dimitrienko Yu.I. Dynamic Transport Phenomena in Porouse Polymer Materials Under Impuls Thermal Effects // Transport in Porouse Media.- 1999.- Vol. 35, no. 3.- P. 299-326.- DOI: 10.1023/A:1006508227447

3. Конюков В.В., Крайнюков В.И., Реш Г.Ф., Иванов М.Ю. Использование гидрокапиллярных лиофобных систем в демпфирующих устройствах и аккумуляторах механической энергии // Аэрокосмические технологии: научные материалы Второй международной научно-технической конференции, посвященной 95-летию со дня рождения академика В.Н. Челомея (Российская Федерация, Реутов - Москва, 19-20 мая 2009 г.) / под ред. Р.П. Симоньянца. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2009.- С. 68.

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

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

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

7. Димитриенко Ю.И., Иванов М.Ю. Моделирование нелинейных динамических процессов переноса в пористых средах // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки.-2008. - № 1.- С. 24-38.

8. David V. Hutton. Fundamentals of Finite Element Analysis.- The McGraw-Hill Companies, 2004.- 494 p.

9. Роуч П. Вычислительная гидродинамика : пер. с англ. В.А. Гущина и В.Я. Митницкого.-М.: Изд-во «Мир», 1980.- 618 с.

10. Anderson John David. Computational fluid dynamics: the basics with applications. - McGraw-Hill, 1995.- 547 р.- (McGraw-Hill series in mechanical engineering).

SCIENTIFIC PERIODICAL OF THE RAIJMAN MS TU

SCIENCE and EDUCATION

EL № FS77 - 48211. №0421200025. ISSN 1994-0408

electronic scientific and technical journal

Modeling of high-speed processes in damping systems with phase transitions # 11, November 2012 DOI: 10.7463/1112.0493560 Dimitrienko Yu.I., Ivanov M.Yu.

Russia, Bauman Moscow State Technical University

IvanovM-fn11@yandex.ru dimit.bmstu@gmail.com

The authors consider problems of modeling processes in new types of damping devices, in which the working fluid is a two-phase system: a porous body - a fluid which does not wet it. On the basis of the theory of heterogeneous porous systems with phase changes they propose a mathematical model of non-stationary processes in a heterogeneous working fluid damper. Specific tasks for calculating micro-and macro parameters of the damping system are formulated. Numerical methods for solving these problems are developed.

Publications with keywords:mathematical modeling, heterogeneous system, porous body, quasi-periodic functions, asymptotic averaging method, cell of periodicity, phase transformation, difference scheme

Publications with words:mathematical modeling, heterogeneous system, porous body, quasi-periodic functions, asymptotic averaging method, cell of periodicity, phase transformation, difference scheme

References

1. Dimitrienko Iu.I. Mekhanika kompozitsionnykh materialovpri vysokikh temperaturakh [Mechanics of composite materials at high temperatures]. Moscow, Mashinostroenie, 1997. 367 p.

2. Dimitrienko Yu.I. Dynamic Transport Phenomena in Porouse Polymer Materials Under Impuls Thermal Effects. Transport in Porouse Media, 1999, vol. 35, no. 3, pp. 299-326. DOI: 10.1023/A:1006508227447

3. Koniukov V.V., Krainiukov V.I., Resh G.F., Ivanov M.Iu. Ispol'zovanie gidrokapilliarnykh liofobnykh sistem v dempfiruiushchikh ustroistvakh i akkumuliatorakh mekhanicheskoi energii [The use of hydro capillary lyophobic systems in the damping devices and batteries of mechanical energy]. Aerokosmicheskie tekhnologii: nauchnye materialy Vtoroi mezhdunarodnoi nauchno-tekhnicheskoi konferentsii, posviashchennoi 95-letiiu so dnia rozhdeniia akademika V.N. Chelomeia [Aerospace Technology: Scientific Proceedings of the Second International Scientific and Technical Conference dedicated to the 95th anniversary of the birth of Academician VN Chelomei], RF, Reutov - Moscow, 19-20 May 2009, Moscow, Bauman MSTU Publ., 2009, p. 68.

4. Nigmatulin R.I. Dinamika mnogofaznykh sred. Ch. 1. [Dynamics of multiphase media. Pt. 1.]. Moscow, Nauka, 1987. 467 p.

5. Sanchez-Palencia E. Nonhomogeneous Media and Vibration Theory. Berlin- New York, Springer- Verlag, 1980. 398 p. (Lecture notes in physics, vol. 127). (Russ. ed.: Sanches-Palensiia E. Neodnorodnye sredy i teoriia kolebanii. Moscow, Mir, 1984. 472 p.).

6. Bakhvalov N.S., Panasenko G.P. Osrednenieprotsessov vperiodicheskikh sredakh. Matematicheskie zadachi mekhaniki kompozitsionnykh materialov [Averaging of processes in periodic media. Mathematical problems of mechanics of composite materials]. Moscow, Nauka, 1984. 352 p.

7. Dimitrienko Iu.I., Ivanov M.Iu. Modelirovanie nelineinykh dinamicheskikh protsessov perenosa v poristykh sredakh [Modeling of nonlinear dynamical processes of transfer in porous media]. VestnikMGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2008, no. 1, pp. 24-38.

8. David V. Hutton. Fundamentals of Finite Element Analysis. The McGraw-Hill Companies, 2004. 494 p.

9. Roache P.J. Computational Fluid Dynamics. Hermosa Publishers, New Mexico, 1976. 446 p. (Russ. ed.: Rouch P. Vychislitel'naiagidrodinamika. Moscow, Mir, 1980. 618 p.).

10. Anderson John David. Computational fluid dynamics: the basics with applications. McGraw-Hill, 1995. 547 p. (McGraw-Hill series in mechanical engineering).