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

УДК 532.51

Многомасштабное моделирование процессов фильтрации жидкого связующего в композитных конструкциях, изготавливаемых методом RTM

© Ю.И. Димитриенко, И.О. Богданов МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: многомасштабное моделирование, фильтрация, метод RTM, композиционные материалы, метод конечных элементов, метод асимптотического осреднения

Введение. Среди большого числа технологий изготовления композитов в последнее время все большее распространение получает метод пропитки армирующего наполнителя связующим в оснастке [1-5]. Существует множество методов, основанных на данном подходе, но большая их часть являются вариациями метода RTM (Resin Transfer Moulding — метод инжекции смолы в закрытую форму), в рамках которого используется жесткая оснастка, состоящая из матрицы и пуансона. Сначала на матрицу выкладывается сухой раскроенный материал, например стеклоткань, стекломат или иной вид армирующего материала, который прижимается пуансоном к матрице. Затем под необходимым давлением в оснастку подается жидкое связующее. Инжекция продолжается до полной пропитки наполнителя, после чего композит отверждается при нормальной или повышенной температуре.

Качество изделия, полученного в результате применения метода RTM, зависит от многих факторов: параметров оснастки, характеристик наполнителя, свойств связующего, рабочего давления пропитки, размеров изделия и др. В связи с этим большое значение имеет про-

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

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

Наиболее эффективным методом математического моделирования процесса фильтрации в пористой системе является метод асимптотического осреднения [11-16]. Методу решения задач нестационарной фильтрации в многофазных композитных структурах посвящены работы [17-20]. Цель настоящей работы — развитие метода асимптотического осреднения для случая многомасштабного моделирования процесса фильтрации с учетом движения фаз как в отдельных порах, так и в пористой среде в целом, для слабосжимаемых жидкостей применительно к задаче течения жидкого связующего в порах сухой композитной конструкции.

Геометрическая модель движения жидкого связующего в пористом каркасе композитного материала. Для простоты будем считать, что композиционный материал имеет форму параллелепипеда и состоит из трех фаз: армирующего наполнителя, газовой фазы, заполняющей поры наполнителя в отсутствии жидкого связующего, и самого жидкого связующего, движущегося по порам наполнителя, вытесняя при этом газовую фазу. Предполагается, что существует подвижная граница раздела £^ (рис. 1) между областями У1 и V,

причем в области У1 есть только наполнитель и жидкая фаза, а в области V — только наполнитель и газовая фаза. Граница раздела £ ^

представляет собой плоскость, параллельную координатной плоскости Ох1 х2 (см. рис. 1) во время ее движения.

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

xi

Рис. 1. Модель процесса пропитки армирующего наполнителя жидким связующим:

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

сти областей V и V соответственно i g

НI

Щ i [--■■■■■ ....................................

II' :: iE Ш [ ' "'Г' •• •• kjÄ::" *

* И:" » §|(

pjiiii;:::»}" "Ё Vg

f. 1 i k - А-/

' ' i 1 1

^ \ ^ 1 Pc

1 1 1 ®

О

У

*2

Ре!

Системы уравнений Навье — Стокса и Стокса для моделирования движения жидкости и газа в порах композита. Под действием давления, создаваемого в пористой среде при методе RTM, происходит движение жидкой фазы по порам, при котором вытесняется газовая фаза. Общая постановка задачи движения слабосжимае-мой жидкости и газа основана на системе уравнений Навье — Стокса, которая в безразмерной форме имеет вид [21]

Pmt + Vi (PmVmi ) = 0

ЕЦ- ((mi )t + Vj (PmVmiVmj )) = "ViPm + (^'m + Km )ViVjVmj + KmVjVjVmi , (1) Pm = AmPm + 5ml (öl " KPOl ),

i = 1,3, j = 1,3, m = l, g.

Здесь Pm — плотность жидкости (газа); Vi — ковариантная производная; vmi — компоненты вектора скорости; X'm, д'т — коэффициенты вязкости жидкости (газа); pm — давление в жидкости (газе); Am = R5mg + K5ml — константа упругости газа и жидкости; Kg = R9 — коэффициент упругости газа; R — безразмерная газовая постоянная; 9 — температура; 5mg и 5ml — символы Кронекера; K = KPl / рг —

безразмерная константа упругости жидкости; K — коэффициент упругости жидкости; POl — гидростатическая часть давления в жидкости, отвечающая за изменение давления без изменения плотности

жидкости, в общем случае является неизвестной переменной величиной; p0l — плотность жидкости в начальный момент времен; Eu m —

число Эйлера, вычисляемое по формуле Eum = J'l ; величины

р v

г m m

с символом л являются характерными значениями соответствующих функций.

Индекс m = g соответствует газовой фазе, а m = l — слабосжи-маемой жидкости. Уравнение состояния в системе (1) при m = g представляет собой уравнение состояния Менделеева — Клапейрона в безразмерной форме для газа, а при m = l — уравнение состояния Маскета [6, 7] для слабосжимаемой жидкости. Модель, для которой p01 = const, назовем классической моделью слабосжимаемой жидкости (КСЖ), а модель, в которой p01 является переменой величиной, — обобщенной. Для обобщенной модели слабосжимаемой жидкости (ОСЖ) систему уравнений (1) необходимо дополнить, поскольку p0l является самостоятельной неизвестной величиной.

Уравнение состояния в системе (1) позволяет рассматривать в единой форме все три модели: модель совершенного сжимаемого газа, КСЖ и ОСЖ.

Система (1) дополняется начальными и граничными условиями на гранях параллелепипеда и границе раздела фаз. На поверхностях раздела «твердое тело — жидкость» £sl и «твердое тело — газ» £ sg ставятся граничные условия прилипания

v4 = vgi

= 0. (2)

На нижней £l и верхней £ гранях композита считается заданным перепад внешних давлений

Pi lll = Pel, Pg £ = Peg . (3)

На боковых гранях £ и £g ставится условие непротекания

ntVu L = a niVgi I = (4)

l l£g

На границе раздела фаз £lg формулируются условия для скачков давления и нормальной скорости

[p] = Pc = , n (vu - Vgi ) = 0. (5)

£sg

Здесь рс — безразмерное капиллярное давление; аЕ — безразмерный коэффициент поверхностного натяжения; Я — безразмерный радиус кривизны поверхности раздела фаз.

В начальный момент времени « = «0 выполняются условия

4=«0 = р01, РХ =, = р0я. (6)

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

дх3-

3Z

dt

= Vl3 ( (tXt)

где х3Е («) — координата поверхности раздела фаз.

Безразмерные коэффициенты вязкости в системе (1) можно выразить через числа Эйлера и Рейнольдса

Л, Л ^

. _ Мт"^т _ _ Л' _ т^т

£ * -г^ т-к 5 т

Lpm Euт Rem m Lpm

где Rem = ^^mVmL — число Рейнольдса; L — характерный линейный Mm

размер композиционного материала.

Если ввести малый безразмерный параметр к = // L ^ 1, равный отношению характерных размеров / ячейки периодичности (ЯП) композита, и характерный размер всей пористой среды L, можно построить соотношения между безразмерными комплексами вида [21]

Eu m к* = Eum = 0(1), Rem к" = Rem = 0(1), (7)

где к, п е Z — некоторые целые числа, соответствующие разным режимам течения фаз. Тогда безразмерные коэффициенты вязкости можно записать как

,,' = MmVm = К = кк+",,0 V = VmVm = Kk+nv 0 (8)

Mm = = ^ 0 „ 0 =К Mm, Vm = = К Vm. (8)

LPm EUm Rem LPm

Здесь M0m = У (Eum Rem ) = 0(1), V m = 0(1) — главные части безразмерных вязкостей. В настоящей статье рассмотрен режим течения к = 2 и п = 0. В этом случае

Eu » 0(1), Re = 0(1), ,m =K2,m, vm = k2Vm. (9)

Наличие второго порядка малости в коэффициентах вязкости д' m

и X'm является допущением рассматриваемой задачи и физически

обосновано малой вязкостью жидкости. Из формулы (9) следует, что имеет место медленное Стоксово (безынерционное) движение сла-босжимаемой жидкости и газа в пористой среде. В этом случае нелинейное уравнение Навье — Стокса преобразуется в линейное уравнение Стокса, и система (1) принимает вид

Pmt + Vi (Pm0 = 0

-Vp +К2 (( +K2MmV;V;Vm! = 0,

Pm = AmPm + ^ml (Pol - KPol )

i = 1,3, j = 1,3, m = l, g.

Для случая обобщенной модели слабосжимаемой жидкости система уравнений Стокса включает в качестве неизвестной гидростатическое давление p0l, входящее в уравнение состояния, поэтому

она является незамкнутой. Доопределим ее следующим образом. Неизвестные функции представим в аддитивном виде:

Pm = Pom + Plm , Pm = Pom +Plm , Vm. = V0m. + V1m., i = 1,3 ™ = /, g. (11)

Причем полагаем, что

Po g = 0, Po g = 0, Vogi = 0, i = 1,3. (12)

После подстановки уравнения (11) в систему уравнений Стокса (10) и разделения ее на две системы получаем:

VjVo/i = 0, x е V,

-V.Pol +к2(0 +^0)V.VjVo/j +KV0Vj.VjVo. = 0, x, eV, (13)

Voli y = 0,

Plmt + V. (PomV1mi + P1mV0m. + P1mV1m.) = 0 Xi eVl,

-VjPim +K2 ( +^m ) Vi V jV1m. + K ^V. V^ = 0, X, G V ,

P1m = AmP1m , = 0.

(14)

s/

Если выполняются обе системы (13) и (14), то автоматически выполняется и исходная система (10). Системы (13) и (14) являются замкнутыми: система (13) представляет собой постановку задачи

о медленном движении несжимаемой жидкости относительно функций p0/, v0/i (для газа эта система отсутствует в силу условия (12)), а система (14) является постановкой задачи движения газа и ОСЖ относительно функций p1m, v1mi, в ней p0m, v0mi рассматриваются

как входные данные.

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

Асимптотический анализ уравнений движений газа и сла-босжимаемой вязкой жидкости в пористой среде. Следуя общей концепции метода асимптотического осреднения [8-12, 22], помимо

глобальных безразмерных координат xi = xjL введем локальные координаты i = xjl = xj к. В этом случае все функции (обозначим их О), описывающие течение газа и жидкости в порах, можно считать квазипериодическими, т. е. зависящими от локальных i координат, глобальных xi координат и времени t. Также введем операцию осреднения функций по области поры V^m в рамках одной ЯП V:

(О m) = Jo mdV, m = g, /. (15)

Vim

Тогда 9m = J dV — пористость среды. Пористость композита

Vm

считаем известной и постоянной величиной.

Решение систем (13) и (14) будем искать в виде асимптотических разложений по малому геометрическому параметру

Vo/i= v002} + 0(к2), p()L = p00} +КР0;} + O( к2),

V1m.=vQ +К«+0(к2), Plm =p(m) +кр® +0(к2), (16)

p.m = pm +Kpm + 0(к2). ( )

Все члены в асимптотических разложениях (16) являются функциями глобальных и локальных координат

0(m) =О( m )(ii, Xj, t), О = (v. 0, p0, vn, p„ P1}, m = 0, 1, ..., (17)

причем по i эти функции полагаются периодическими.

Подставляя разложения (15) в задачу (13) и применяя правило дифференцирования квазипериодических функций [16, 22], получаем асимптотические разложения системы уравнений движения несжи-

маемой жидкости. Члены при отрицательных степенях малого параметра к"1 в этой системе и остальных степенях к", п = 0, 1, ... приравниваем к нулю и получаем:

к"1: V,^ = 0, } = 0, ^ = 0; (18)

к0: V ^ = 0,

-^.р00}" V,гp01l}+(Мн+^ V, ^+^ V,; V, ^ = 0, (19)

у0й = 0.

Из первой системы (18) определяем, что р0°) не зависит от локальных координат , : р(^1>(х\). Тогда, присоединяя к уравнению неразрывности и граничному условию системы (18) уравнение равновесия жидкости из системы (19), получаем локальную задачу на ЯП для несжимаемой жидкости в нулевом приближении:

V,гv0ю) = 0,

-V *гр00} - ^р01}+КV, V, ^ = 0,

(20)

х, = 0,

у(0) 0Н

< р0? >= 0.

Неизвестными в этой системе являются функции у0<°^), р0,), а градиент давления V^р^) рассматривается как входные данные. К системе (20) присоединяются условия периодичности функций, которые обозначаются как

[[О]] = 0 о О(,,2,-1/2) = О((,,2,1/2),

1,-1/2,,з) = О(,1,1/2,, ), О(-1/2,,2,,з) = О(1/2,,2,,з), (21)

О = (у(0) р0)) 1 г \ У0И > У01 ) •

Здесь -1/2 < , < 1/2 — область ячейки периодичности У^.

Условие нормировки < р^ >= 0 в системе (20) введено для единственности решения задач (20), (21).

Общее решение локальной задачи (20) можно представить как формальную зависимость скорости и давления у0°г), р01 от входных

данных — градиента V ^р00}

С =—Я 1 )У Л>, ^ = р(1, )У ^>. (22)

М21

Здесь Я. 1), р(1) — структурные функции, которые зависят

только от локальных координат и вычисляются на основе решения задачи (20) на ЯП.

Осредняя уравнение неразрывности в системе (19), с учетом того, что для любой периодической функции < V,.. О >= 0, получаем

осредненное макроскопическое уравнение движения несжимаемой жидкости

V,< v00 >= 0. (23)

Подставляя (22) в (23), получаем осредненное уравнение для вычисления давления р0О)

V .(К^ „Р00У) = 0, (24)

где

Кщ = —— < Я1ЧЬ) >, < С >= —Кщ(Ь^^ (25)

являются компонентами тензора фильтрации несжимаемой части жидкости в пористой среде.

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

к—1 : V (п у(°) + п(%№ + п(°У°)) = 0

к : У ь (П0шЧшг + П1ш У0ш. + П1ш У1шг ) = 0,

Р(0) = А п(°). V Р(°) = 0 V(0М = 0

Пш 1ш > У Ьг/Лш и> У1шг '

к0: п(0) + ^ (2 + п^ш5 (vО0Уi +vШ)) +

+ п^ос + )+п(ш)(е + V®) )=0,

—?.Рш —V,Рш +(М1ш +М2ш ) + М = 0,

(26)

(27)

л(1) = А п(1), V!1

1ш 1ш' 1шг

Е = 0.

Из системы (26) получаем, что р) не зависит от локальных координат , следовательно, и р^ не зависит от : р00)(х}), р^Чх}).

Тогда, присоединяя к уравнению неразрывности и граничному условию системы (26) уравнение равновесия сжимаемой части жидкости из системы (27), получаем локальную задачу на ЯП для сжимаемой части жидкости в нулевом приближении:

V,. (р0 +р(0С +р(0Ч) ) = О,

- V,.р) + (М- + М-) V V,;у1(0). + м-V,;V,;<) = 0, (28)

V(0) = 0 < Р(1) >= 0

Неизвестными в этой системе являются функции V0, р,-. Величины , р(Я> и градиент давления Vхгр1(°) рассматриваются как входные данные.

Преобразуем уравнение неразрывности в систему (28). Поскольку

р(Я-( X), р0я не зависят от , получаем:

V (р V(0) + Р(0)V(0) +Р(0)V(0)) = р V У(0) +Р(0)V V(0) +Р(0)V V(0) = 0 (29)

у 0т у1тг^ г 1т у0т. ^ Г1т у1тО К0я у Ц у1тг ^ Г1т у 0т. ^ Ит у , у1тг )

Скорость v00)i удовлетворяет уравнению несжимаемости (20) V,^0т)г. = 0, поэтому (29) сводится к уравнению V, у) = 0, а задача (28) принимает вид

v,Vl(m)=0,

-V „р) +(М1т + м2тV,V) + ;V,;<) = 0, (30)

V?)

1тг

=0, < рот >=0.

Задача (30) формально в точности совпадает с задачей (20), поэтому ее решение имеет такую же структуру, как и (22):

У(0) = — Ж(;)(, )V Р(0) Р(1) = Р( 1)(, )V Р(0)

г1тг ''1т. \±>к; у х]У1т ' У1т 1 1т у х}¥1т '

м

2т

< V) >=-ктуVхР?, кту =-—< Ж1(т.)(,к) >. (31)

М2т

Осредняя уравнение неразрывности в системе (27), с учетом того, что для всякой периодической функции < V,г.fi >= 0, получаем

осредненное уравнение движения для газа и сжимаемой части жидкости

© о(0) +о V <V(0) >+У (о(0)(< У(0) > + <у(0) >)) = 0 (32)

УтГ \тг И0шу хг ^ х^К\т \ ^ ^0тг ^ ^ ^ У\тг ^^

С учетом системы (23) это уравнение можно записать как

© о(0) +V ((о +о(0)) <у(0) >)+<у(0) ^ о(0) = 0 (33)

УтК\тг^ у хг И0т ^ К\т / ^ Мм у0тг ^ у хгУ\т и' )

Подставляя в выражение (33) уравнения (25) и (32), получаем искомое осредненное уравнение для вычисления давления р\П :

© о(0) ^ ((о +о(0))К\ V Р(0)) + К V Р(0)V о(0) т = I я (34)

Ч^тКт у хг ^ЧК0т ^ к\т )1^тг] у ху//\т ^ ^ ту у здУ0ту хгУ\т ' ^^

Здесь введено формальное определение тензора фильтрации для газа К%ц = 0, тензор фильтрации К1Ц для жидкости определяется по формулам (25), тензор фильтрации К1т11- вычисляется по формуле (31).

Решение локальных задач. Далее рассмотрим многоканальную пористую структуру (см. рис. \), у которой поры образуют продольные каналы вдоль каждого координатного направления 0,г. Такая

структура возникает у многих пористых материалов, в том числе у не-пропитанных композиционных материалов, армированных тканями.

Подставляя формулу (22) в задачу (20), после исключения градиента V х;р0°} получаем набор локальных задач:

vW1 ) = 0,

-V?1 = ^, , ^, (35)

, .,= 0, , ^,

< р(1) >= 0, ) ]] = 0, [[ р(1) ]] = 0

для определения функций р(1) (,г), ), которые, в отличие от

задачи (20), не содержат никаких свойств жидкости и не зависят от входных данных. Здесь А, = V,: . V, . — оператор Лапласа.

Решение задач (35) определяется только внутренней геометрией пор. Функции И( 1) вычисляются следующим образом:

.) | 0, г * 1 или (г = 1 и V^00) = 0), ^) = 1 _ (0) - 1 (36)

т1;

|\, г = 1 и V х,р00) * 0.

Введем компоненты тензора напряжений Т^) и тензора скоро

стей деформаций ) в жидкости:

1гк

1гк 13 'ш"т;гттъ"г>'',гт,т;г

1%) = -(1 ^ + 2Б»\ Бк) = 1 ))- (37)

Тогда уравнение равновесия жидкости в системе (35) перепишем в виде

^) = 1). (38)

Применяя теорему о симметричном и антисимметричном продолжении, аналогичную приведенной в работах [16, 22], можно сформулировать локальную задачу (35) на 1/8 ЯП с граничными условиями на противоположных гранях 1/8 ЯП:

^ = 0; Ъ = 2

(■) -Ж(1) г и

Ж1) (5]к + 5И - 25]к5и) + [1 - (5]к + 5и - 25]к5И) ] = 0;

дР(1) _

Р(1) 5 к +-Р- (1 -5 ]к) = 0; 1, и к = 1,3.

(39)

Каждая из задач системы (35) представляет собой стационарную задачу течения фиктивной линейно-вязкой несжимаемой среды. Решение задач (35) зависит только от внутренней геометрии пор, поэтому их постановка применима для расчетов фильтрации любых газов и жидкостей, которые могут быть рассмотрены в рамках выбранных моделей.

Аналогичным образом в силу условия (31) можно представить решение задач (30). Поскольку эти задачи подобны, после исключения коэффициентов вязкости структурные функции ^1(/), )

и Ж. 1) совпадают, поэтому

К = Кги = ---- < Ж1 )£к) >, К: =--^ < Ж) >,

—2, —2 g

< Сг >=-Кт11 VX), т = I, g. (40)

Методы численного решения локальных задач (35) разработаны в работах [22-25].

Постановка макроскопических задач фильтрации газов и жидкостей. Второе соотношение (25) представляет собой закон фильтрации Дарси для макроскоспического движения жидкости [22], а осредненное уравнение (24) является уравнением фильтрации несжимаемой части жидкости в порах композита

V^(К^^}) = 0, X; е¥1. (41)

С учетом того, что р^ не зависит от локальных координат и имеет место уравнение состояния р^ = Атр^, после его подстановки в систему (34) получаем второе уравнение фильтрации для газовой фазы и сжимаемой части жидкости в порах композита

=М(РотАт +Р!т)Жту^А(т) ^^РтЧ^, т = I, g, (42)

где К'ту = К'т,] / Фт , К т,} = Кт,} 1 Фт •

Поскольку для газа v0g) = О, К^ = 0 и р0 = 0, из формулы (42) получаем уравнение фильтрации газовой фазы:

Р(2 = ^ (Р^К^Р^ х}. ^, (43)

а для сжимаемой части жидкости

р(0} (( +р1(°))К1 )+Кй;ух;р00)У„.Р1(0), х;. (44)

Здесь градиент давления Vхз-р0)0) является известной величиной,

его можно найти из решения задачи (41).

Метод введения динамических независимых координат для макроскопических задач фильтрации. Рассмотрим макроскопические задачи (41), (43) и (44) для случая, когда области У1 и V представляют собой параллелепипед, а граница раздела между ними

— плоскость, движущаяся вдоль координаты х3 в направлении

области V. Эта задача более похожа на классическую задачу Стефана [26] с неизвестной заранее границей раздела взаимодействующих сред, положение и форма которых меняется с течением времени, но является существенно более сложной ввиду того, что сама система уравнений (43) и (41), (44) различна по разные стороны от поверхности раздела, кроме того, эти системы уравнений нелинейны, и уравнение движения поверхности раздела также нелинейное.

Применим метод введения динамических независимых координат, в которых положение границы раздела и области газовой и жидкой фаз не меняются в течение всего рассматриваемого времени. Для этого введем новые координаты

Х1 = х1, Х2 = х2, Х3 = —^ + Ь. (45)

а + кх3

Здесь коэффициенты а, Ь и к определяются из условий

Х3 = 0 при х3 = 0,

Х3 = Ь при х3 = Ь, (46)

Х3 = х3Е (0) пРи х3 = х3Е ^).

Тогда производные по времени и пространственным координатам можно представить в виде

О, + Ь5,зУлО,

V „о^ , Х] о,

(47)

где V х О =

ао

ах,.

а также

ь = --а^кхг+Ь(, а =5,-5, г = 1+

(а + кх3)

(а + кх3)

(48)

После подстановки (47) в уравнения (41), (43) и (44) получаем формулировку макроскопической задачи неустановившейся совместной фильтрации слабосжимаемой жидкости и газа, которая состоит из уравнений

ЖьО*УхкК0 = о, х, ¿V; = 0,пУхп(н^кР™)-И5,зУх,р(0\ х, ¿V;

(49)

(50)

р(0} = (Н1,кVхкРи}) + (ЯЛр00} -Ь5,з)Ух,р(0), х, ¿Г1. (51)

На нижней Е1 и верхней Ег гранях композита считается задан-

(52)

"I----^--------я

ным перепад внешних давлений

Р (°) У о 1

= Ре1, Р1(0

= 0,

р£05

= Рея.

На боковых гранях Ег' и Е^ задаются условия непротекания

п, КОVхкР00}| = о, п, Н,кVхкР^!Е, = о, (5з)

п1Н11к V хкР^ |Е,= 0.

На границе раздела фаз Е1я задаются условия для скачков давления и нормальной скорости

(( - Р10 ) = ^, Р(0) = 0, п, (КщЯ,кVхкР00) - Н^Р^ ) = 0. (54) В начальный момент времени I = выполняются условия

Р (0) = Р Р(0) Р01 ¥01? ¥11

М—1с

= 0,

I = Р0 я.

(55)

Многомасштабное моделирование процессов фильтрации... Здесь обозначены

я*=М+р(0))), я1г]=, =р1(0)кк;уб;к. (56)

Вариационные постановки макроскопической задачи неустановившейся фильтрации. Все три уравнения фильтрации (49)-(51) имеют общую структуру

Яр! = 0V, ЯVхкр ) + Л. VХ.р, (57)

где Я, Нк и Я. определяются самими формулами (49)—(51), а

р = {р00}, А?, р(0)}.

Рассмотрим уравнение (57) и построим для него вариационную формулировку. Введем специальные классы:

1) скалярных переменных полей $(Х,г), определенных в V и Е,

V! > 0, непрерывно-дифференцируемых один раз по X. иг и удовлетворяющих нулевому граничному условию на частях Е1 и Е поверхности Е области V:

^ UЕg = 0 (58)

2) стационарных скалярных полей Р(Х), определенных в области V.

Тогда, умножая (57) на $ и интегрируя по области V, получим: К - % V* НVхкр ) - Я, VХ.р) = 0. (59)

V

Для начального условия

р' = р0,Х е V, г = 0 (60)

можно записать

|(р'-р, )pdV = 0. (61)

V

Здесь р' = р (Х. ,0).

Вычитая равенство (61) из формулы (59), получаем:

|бь^-м-V„(НкVхкр -VХ^ -|(р'-р0)№ = 0.

V с V V V

(62)

Преобразуем второй интеграл в формуле (62). В силу того что на поверхностях Е1 и Е выполняется условие (53), на боковых поверхностях образца — условие непроницаемости, а на границе раздела — условие равенства нормальных скоростей (5), получаем:

¡ЩV,НVкР = -{V,тн)Ухкр йУ. (63)

V У

Будем называть нестационарное поле р(X,I) (X еУ,I > 0) допустимым, если оно является один раз непрерывно-дифференцируемым по Xг в V и Е и I и удовлетворяет граничному условию

(3). Допустимым стационарным полем р'(X) будем называть произвольное скалярное поле в V.

Выберем $ и Р в виде $ = 5р, Р = 5р' . Тогда из формул (62) и (63) получаем:

{Жр др-йУ -¡V,Ьр(ан к^лр йУ -

У У (64)

- ¡5рЯг VХгрйУ -¡(р '-ро)5р' йУ = 0.

У У

В силу произвольности вариаций 5р и 5р ' получаем искомую вариационную формулировку задачи:

{£5рдр-йУ - ¡V,5р(<2уДк)Vлр йУ - ¡5рЯг VхрйУ = 0; (65)

у У У

¡(р' - ро)5р'йУ = 0. (66)

У

Метод конечных элементов для решения глобальной задачи фильтрации. Для численного решения вариационной задачи (65), (66) глобальной задачи будем использовать метод конечного элемента [16, 24] с тетраэдральным симплекс-элементом. Введем координатный столбец значений нестационарного поля давления р в узлах конечного элемента и матрицу функций формы

{Р} = (р Р2 Ръ РА); (67)

1x4

[ N ] = (( N2 N3 N4) = Ц, г = 1,4. (68)

1x4

Здесь Ц — Ц -координаты [27, 28]. Тогда давление и вариация давления на конечном элементе могут быть представлены как

р = [N]{Р} 5р = [N]{5Р}.

(69)

Аналогично для стационарного поля давления р' и его вариации 5р' можно получить:

р' = [ N ]{Р'}, 5р' = [ N ]{5Р'}. (70)

1x4 4x1 1x4 4x1

Здесь введены вектор значений стационарного поля давления в узлах конечного элемента

{Р'}Т =( р р' Р Р).

(71)

Далее введем векторы градиентов стационарного и нестационарного полей давлений:

Ыг =

(

др др др

Л

дХ1 дХ2 дХ

, {g Г =

(

3

др' др' др' ~дХ1 ~дХ~2 дХ

Л

(72)

3 У

После подстановки соотношений (70) в формулы (72) получаем выражения для градиентов:

{^ = [В ]{Р}, {g '} = [ В ]{Р'}, (73)

3x1 3x4 4x1

где введена матрица производных функций формы [ В ].

3x4

Подставляя полученные соотношения в вариационную постановку (65), получаем:

д{р}

[С + К ]{Р} = К }.

4x4 дг 4x4 4x1 4x1

Здесь введены матрицы и координатный столбец правой части

Л

(74)

т

т 2 I 4x4 У

[С] = ]М [N]dV, [Кт ] = [Кт1 уШтг >[К

4x4 V 4x1 1x4 4x4 4x4 2 ^ 4x4

[Кт1 ] = |[В ] [Ат ][В [ Кт 2 ] = ] {Ьт } [В ] dV, (75)

4x4 V 4x3 3x3 3x4 4^ V 4x1 1x3 3x4

{т } = //т [N]TdV.

4x1 V 4x1

Решение системы (75) может быть осуществлено методом конечных разностей. Пусть гп — значение времени на п -м временном

слое, гп+1 — значение времени на следующем временном слое,

т = гп+1 - гп — шаг по времени. Тогда конечно-разностный аналог си-

стемы (75) можно получить с помощью неявной схемы на основе правила трапеций [29]. Для производной по времени получаем:

(76)

& т

Предполагается, что производная вычисляется в средней точке временного интервала, поэтому все функции времени также должны быть вычислены в этой точке. Для этого можно использовать при, {Р' , + {Р'

ближенное соотношение {Р} =-^^--. Тогда с учетом выражения (76) уравнение (74) записываем в виде

[с]{Р}"+'~{Р}" +[кт]' {Р}п+'+{Р}" ={Гт}', т = ¡,Е. (77)

Т 2 4x1

Здесь символ ' означает принадлежность средней точке временного интервала. Выполняя группировку в последнем уравнении, приходим к искомой СЛАУ для каждого отдельного конечного элемента:

[ Ат ]>т}„+1 =[К ]>т}„ + 2{т }' , (78)

где введены матрицы

[ Ат ]' = 2 [С] + [Кт ]' , [Кт ]' = ^ [С] " [Кт ]' . (79)

Отметим, что данная конечно-разностная схема построена на основе центрально-разностной схемы Кранка — Николсона, которая является безусловно устойчивой [23].

Моделирование процесса пропитки армирующего материала композита жидким связующим. Рассмотрим моделирование процесса пропитки армирующего материала жидким связующим на основе разработанного численного метода. Свойства сред и значения основных параметров приведены в таблице, результаты расчетов отображены на рис. 2 и 3.

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

На рис. 2 и 3 вертикальными прерывистыми линиями показано положение поверхности раздела фаз (фронта пропитки жидким связующим) в разные моменты времени. Распределение давления в композиционном материале в процессе пропитки при методе ЯТМ изменяется монотонным образом: от максимального значения на поверхности £ подачи связующего под давлением до минимального значения на поверхности £ , где происходит откачка воздуха. По мере распростра-

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

Значения основных параметров расчета

Параметр Значение

Размеры образца, м 2 х 0,2 х 0,05

Характерный размер ячейки пористого материала (предполагается одинаковым для всех фаз), м 5 -10-6

Безразмерные коэффициенты проницаемости пористой среды К11 = 0,003687676 К22 = 0,002496272 К33 = 0,002496272

Вязкость ц жидкости, Па с 0,2

Вязкость ц'^ газа (воздуха), Па-с 1,81 -10-5

Коэффициент упругости жидкости К, м2/с2 106

Коэффициент упругости газа К , м2/с2 105

Давление р01, р0 в фазах в начальный момент времени, Па 105

Давление ре1 жидкости на внешней границе, Па 106

Давление р газа на внешней границе, Па 105

0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0 Продольная координата, м

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

--0,1 (3 с);-- 1,005022 (48 с);-- 1,714717 (96 с);-- 0,340197 (12 с);

— 1,191137 (60 с);--1,884458 (108 с);--0,603497 (24 с);--1,369582 (72 с);

-- 1,96899 (114 с);--0,810645 (36 с);-- 1,543524 (84 с)

о £

о о.

о §

я

0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0 Продольная координата, м

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

= — 0,1 (3 с);-- 1,005022 (48 с);-- 1,714717 (96 с);-- 0,340197 (12 с);

■ — 1,191137 (60 с);--1,884458 (108 с);--0,603497 (24 с);--1,369582 (72 с);

-- 1,96899 (114 с);--0,810645 (36 с);-- 1,543524 (84 с)

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

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

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

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

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

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

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

ЛИТЕРАТУРА

[1] Джоган О.М., Костенко О.П. Методы изготовления деталей из композиционных материалов пропиткой в оснастке. В 2 ч. Ч. 1. Методы пропитки под давлением. Вопросы проектирования и производства конструкций летательных аппаратов. Харьков, ХАИ, 2011, № 4 (68), с. 111-125.

[2] Джоган О.М., Костенко О.П. Методы изготовления деталей из композиционных материалов пропиткой в оснастке. В 2 ч. Ч. 2. Методы вакуумной пропитки. Вопросы проектирования и производства конструкций летательных аппаратов. Харьков, ХАИ, 2012, № 1 (69), с. 80-92.

[3] Джоган О.М., Костенко О.П. Практическая классификация методов изготовления деталей из полимерных композиционных материалов пропиткой в оснастке. Вопросы проектирования и производства конструкций летательных аппаратов. Харьков, ХАИ, 2013, № 1 (73), с. 21-32.

[4] Chen X., Zhang Y., Shilin Yan S. Two-dimensional simulations of resin flow in dual-scale fibrous porous medium under constant pressure. Journal of Reinforced Plastics and Composites, 2013, vol. 32 (22), pp. 1757-1766.

[5] Yang B., Tang Q., Wang S., Jin T., Bi F. Three-dimensional numerical simulation of the filling stage in resin infusion process. Journal of Composite Materials, 2016, vol. 50 (29), pp. 4171-4186.

[6] Маскет М. Течение однородных жидкостей в пористой среде. Москва, Ижевск, Институт компьютерных исследований, 2004, 628 с.

[7] Ho C.K., Webb S.W., eds. Gas transport in porous media. Springer, 2006, 440 p.

[8] Баренблатт Г.И., Ентов В.М., Рыжик В.М. Движение жидкостей и газов в природных пластах. Москва, Недра, 1984, 211 с.

[9] Jenny P., Lee S.H., Tchelepi H.A. Adaptive multiscale finite-volume method for multiphase flow and transport in porous media. Multiscale Modeling and Simulation, 2004, 3 (1), pp. 50-64.

[10] Nield D.A., Bejan A. Convection in porous media. New York, Springer-Verlag, 2006. DOI 10.1007/978-1-4614-5541-7

[11] Бардзокас Д.И., Зобнин А.И. Математическое моделирование физических процессов в композиционных материалах периодической структуры. Москва, Едиториал УРСС, 2003, 376 с.

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

[13] Беляев А.Ю. Усреднение в задачах теории фильтрации. Москва, Наука, 2004, 200 с.

[14] Победря Б.Е. Механика композиционных материалов. Москва, Изд-во Моск. ун-та, 1984, 336 с.

[15] Санчес-Паленсия Э. Неоднородные среды и теория колебаний. Москва, Мир, 1984, 472 с.

[16] Димитриенко Ю.И., Кашкаров А.И., Макашов А.А. Конечно-элементный расчет эффективных упругопластических характеристик композитов на основе метода асимптотического осреднения. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2007, № 1, с. 26-46.

[17] Dimitrienko Yu.I., Dimitrienko I.D. Effect of thermomechanical erosion on heterogeneous combustion of composite materials in high-speed flows. Combustion and Flame, 2000, vol. 122, no. 3, pp. 211-226.

[18] Dimitrienko Yu.I. Thermal stresses in ablative composite thin-walled structures under intensive heat flows. International Journal of Engineering Science, 1997, vol. 35, no. 1, pp. 15-31.

[19] Димитриенко Ю.И., Яковлев Д.О. Асимптотическая теория термоупругости многослойных композитных пластин. Механика композиционных материалов и конструкций, 2014, т. 20, № 2, с. 260-282 с.

[20] Dimitrienko Yu.I. A structural thermomechanical model of textile composite materials at high temperatures. Composite science and technologies, 1999, vol. 59, pp. 1041-1053.

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

[22] Dimitrienko Yu.I., Dimitrienko I.D. Simulation of local transfer in periodic porous media. European Journal of Mechanics/B-Fluids, 2013, no. 1, pp. 174-179.

[23] Димитриенко Ю.И., Захарова Ю.В., Богданов И.О. Математическое и численное моделирование процесса фильтрации связующего в тканевом композите при RTM методе изготовления. Университетский научный журнал, 2016, № 19, с. 33-43.

[24] Димитриенко Ю.И., Левина А.И., Боженик П. Конечно-элементное моделирование локальных процессов переноса в пористых средах. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2008, № 3, с. 90-103.

[25] Димитриенко Ю.И., Шпакова Ю.В., Богданов И.О., Сборщиков С.В. Моделирование процесса многоуровневой фильтрации жидкого связующего в тканевом композите при RTM-методе изготовления. Инженерный журнал: Наука и инновации, 2015, вып. 12. DOI 10.18698/2308-6033-2015-12-1454

[26] Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. Москва, Едиториал УРСС, 2003, 784 с.

[27] Сегерлинд Л. Применение метода конечных элементов. Москва, Мир, 1979, 392 с.

[28] Zienkiewicz O.C, Taylor R.L., Zhu J.Z. The finite element method. Its Basis and fundamentals. Butterworth-Heinemann, 2013, 756 p.

[29] Амосов А.А, Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. Москва, Высшая школа, 1994, 544 с.

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

Ссылку на эту статью просим оформлять следующим образом: Димитриенко Ю.И., Богданов И.О. Многомасштабное моделирование процессов фильтрации жидкого связующего в композитных конструкциях, изготавливаемых методом RTM. Математическое моделирование и численные методы, 2017, № 2, с. 3-27.

Димитриенко Юрий Иванович — д-р физ.-мат. наук, профессор, директор Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана, заведующий кафедрой «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор более 350 научных работ в области механики сплошных сред, вычислительной механики, механики и термомеханики композитов, математического моделирования в науке о материалах, вычислительной газодинамики. e-mail: [email protected]

Богданов Илья Олегович — аспирант кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Специалист в области численного моделирования процессов фильтрации. e-mail: [email protected]

Multiscale modeling of liquid binder filtration processes in composite structures manufactured by RTM

© Yu.I. Dimitrienko, I.O. Bogdanov Bauman Moscow State Technical University, Moscow, 105005, Russia

A mathematical model for the multiscale process of filtration of weakly compressible liquids and gases in periodic porous media is proposed with reference to the process of composite material production based on the RTM method. Using the method of asymptotic averaging made it possible to formulate the so-called local filtration problems for a single pore and the global problem of unsteady filtration of weakly compressible liquids. Two models of a weakly compressible fluid are considered: classical and generalized. The classical model is based on the Musket's equation of the state, which requires initial constant values for fluid pressure and density to be preset. The generalized model is based on the same equation, but requires presetting only the initial fluid density, using the unknown hydrostatic pressure instead of the initial constant liquid pressure. The results of simulation of the impregnation process of a of filler material sample by a binder are presented using the two models of a weakly compressible liquid.

Keywords: multiscale modeling, filtration, RTM method, composite materials, finite element method, asymptotic averaging method

REFERENCES

[1] Dzhogan O.M., Kostenko O.P. Metody izgotovleniya detaley iz kompo-zitsionnykh materialov propitkoy v osnastke. V 2 chastyakh. Chast 1. Metody pro-pitki pod davleniem [Methods of manufacturing parts of composite materials by impregnating in a mounting. In 2 parts. Part 1. Methods of impregnation under pressure]. Sbornik trudov "Voprosyproektirovaniya iproizvodstva konstruktsiy le-tatelnykh apparatov" [Proc. of the Problems of design and production of aircraft structures]. Kharkov, Kharkov Aviation Institute, 2011, no. 4 (68), pp. 111-125.

[2] Dzhogan O.M., Kostenko O.P. Metody izgotovleniya detaley iz kompo-zitsionnykh materialov propitkoy v osnastke. V 2 chastyakh. Chast 2. Metody vakuumnoy propitki [Methods of manufacturing parts of composite materials by impregnating in a mounting. In 2 parts. Part 2. Vacuum impregnation methods]. Sbornik trudov "Voprosy proektirovaniya i proizvodstva konstruktsiy letatelnykh apparatov" [Proc. of the Problems of design and production of aircraft structures]. Kharkov, Kharkov Aviation Institute, 2012, no. 1 (69), pp. 80-92.

[3] Dzhogan O.M., Kostenko O.P. Prakticheskaya klassifikatsiya metodov izgotovleniya detaley iz polimernykh kompozitsionnykh materialov propitkoy v osnastke [Practical classification of methods for manufacturing parts of polymer composite materials by impregnating in a mounting]. Sbornik trudov "Voprosy proektirovaniya i proizvodstva konstruktsiy letatelnykh apparatov" [Proc. of the Problems of design and production of aircraft structures]. Kharkov, Kharkov Aviation Institute, 2013, no. 1 (73), pp. 21-32.

[4] Chen X., Zhang Y., Shilin Yan S. Journal of Reinforced Plastics and Composites, 2013, vol. 32 (22), pp. 1757-1766.

[5] Yang B., Tang Q., Wang S., Jin T., Bi F. Journal of Composite Materials, 2016, vol. 50 (29), pp. 4171-4186.

[6] Muskat M. The flow of homogeneous fluids through porous media. McGraw-hill book company, Inc. Publ., 1937 [In Russ.: Muskat M. Techenie odnorodnykh zhidkostey v poristoy srede. Moscow, Izhevsk, Institut kompyuternykh issledo-vaniy Publ., 2004, 628 p.].

[7] Ho C.K., Webb S.W., eds. Gas transport in porous media. Springer, 2006, 440 p.

[8] Barenblatt G.I., Yentov V.M., Ryzhik V.M. Dvizhenie zhidkostey i gazov v pri-rodnykh plastakh [Movement of liquids and gases in natural seams]. Moscow, Nedra Publ., 1984, 211 p.

[9] Jenny P., Lee S.H., Tchelepi H.A. Multiscale Modeling and Simulation, 2004, 3 (1), pp. 50-64.

[10] Nield D.A., Bejan A. Convection in porous media. New York, Springer-Verlag Publ., 2006. DOI 10.1007/978-1-4614-5541-7

[11] Bardzokas D.I., Zobnin A.I. Matematicheskoe modelirovanie fizicheskikh protsessov v kompozitsionnykh materialakh periodicheskoy struktury [Mathematical modeling of physical processes in composite materials of periodic structure]. Moscow, Editorial URSS Publ., 2003, 376 p.

[12] Bakhvalov N.S., Panasenko G.P. Osrednenie protsessov v periodicheskikh sredakh. Matematicheskie zadachi mekhaniki kompozitsionnykh materialov [Averaging processes in periodic media. Mathematical problems of the composite material mechanics]. Moscow, Nauka Publ., 1984, 352 p.

[13] Belyaev A.Yu. Usrednenie v zadachakh teorii filtratsii [Averaging in problems of the filtration theory]. Moscow, Nauka Publ., 2004, 200 p.

[14] Pobedrya B.E. Mekhanika kompozitsionnykh materialov [Mechanics of composite materials]. Moscow, Lomonosov Moscow State University Publ., 1984, 336 p.

[15] Sanchez-Palencia E. Non-Homogeneous media and Vibration Theory. Berlin Heil-deberg New York, Springer-Verlag Publ., 1980 [In Russ.: Sanchez-Palencia E. Neodnorodnye sredy i teoriya kolebaniy. Moscow, Mir Publ., 1984, 472 p.].

[16] Dimitrienko Yu.I., Kashkarov A.I., Makashov A.A. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki — Herald of the Bauman Moscow State Technical University. Series Natural Sciences, 2007, no. 1, pp. 26-46.

[17] Dimitrienko Yu.I., Dimitrienko I.D. Combustion and Flame, 2000, vol. 122, no. 3, pp. 211-226.

[18] Dimitrienko Yu.I. International Journal of Engineering Science, 1997, vol. 35, no. 1, pp. 15-31.

[19] Dimitrienko Yu.I., Yakovlev N.O. Mekhanika kompozitsionnykh materialov i konstruktsiy — Mathematical modeling and Computational Methods, 2014, vol. 20, no. 2, pp. 260-282.

[20] Dimitrienko Yu.I. Composite science and technologies, 1999, vol. 59, pp. 1041-1053.

[21] Dimitrienko Yu.I., Ivanov M.Yu. Vestnik MGTU im. N.E. Baumana. Seria Estestvennye nauki — Herald of the Bauman Moscow State Technical University. Series: Natural Sciences, 2008, no. 1, pp. 39-56.

[22] Dimitrienko Yu.I., Dimitrienko I.D. European Journal of Mechanics/B-Fluids, 2013, no. 1, pp. 174-179.

[23] Dimitrienko Yu.I., Zakharova Yu.V., Bogdanov I.O. Universitetskiy nauchnyy zhurnal — Humanities and Science University Journal, 2016, no. 19, pp. 33-43.

[24] Dimitrienko Yu.I., Levina A.I., Bozhenik P. Vestnik MGTU im. N.E. Baumana. Seria Estestvennye nauki — Herald of the Bauman Moscow State Technical University. Series: Natural Sciences, 2008, no. 3, pp. 90-103.

[25] Dimitrienko Yu.I., Shpakova Yu.V., Bogdanov I.O., Sborshchikov S.V. Inzhe-nernyy zhurnal: Nauka i innovatsii — Engineering Journal: Science and Innovation, 2015, no. 12. DOI 10.18698/2308-6033-2015-12-1454

[26] Dimitrienko Yu.I., Samarsky A.A., Vabishchevich P.N. Vychislitelnaya teplop-eredacha [Computational heat transfer]. Moscow, Editorial URSS Publ., 2003, 784 p.

[27] Segerlind L.J. Applied Finite Element Analysis. John Wiley & Sons Publ., 1976 [In Russ.: Segerlind L. Primenenie metoda konechnykh elementov. Moscow, Mir Publ., 1979, 392 p.].

[28] Zienkiewicz O.C, Taylor R.L., Zhu J.Z. The finite element method. Its basis and fundamentals. Butterworth-Heinemann Publ., 2013, 756 p.

[29] Amosov A.A., Dubinsky Yu.A. Kopchenova N.V. Vychislitelnye metody dlya inzhenerov [Computational techniques for engineers]. Moscow, Vysshaya shko-la Publ., 1994, 544 p.

Dimitrienko Yu.I., Dr. Sc. (Phys.-Math.), Professor, Head of the Scientific and Educational Center for Supercomputer Engineering Modeling and Development of Software Complexes, Bauman Moscow State Technical University. Head of Computational Mathematics and Mathematical Physics Department at the Bauman Moscow State Technical University. Author of over 350 research papers in the field of continuum mechanics, computational mechanics, mechanics and thermomechanics of composites, mathematical modeling in the science of materials, computational gas dynamics. e-mail: [email protected]

Bogdanov I.O., Post graduate (Ph.D.) student, Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. Science research interests: numerical modeling filtration processes. e-mail: [email protected].