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

УДК 510:53.072:621.1.016.4

DOI 10.18698/2309-3684-2017-4-1730

Моделирование стохастических процессов фильтрации в решетчатых системах

© Р.В. Арутюнян МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: фильтрация, диффузия, кинетика, стохастическое уравнение, существование, единственность, численный метод

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

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

В работах [15, 16] реализуется подход Бринкмана в случае фильтрации жидкости сквозь твердый недеформируемый пористый материал, рассматриваемый как периодическая среда, в которой выделяется ряд характерных масштабов и применяется асимптотическое усреднение. Для каждого из масштабов формулируется система уравнений. Нахождение усредненных фильтрационных характеристик среды, а также скорости и давления жидкости сводится к решению соответствующих периодических задач на ячейке.

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

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

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

Математическая модель. Рассмотрим одномерную периодическую структуру типа «решето» (рис. 1). Длина непроницаемой части равна двум, проницаемой (отверстие) — единице. Сквозь это «решето» просеивается поток одномерных частиц (стержней) случайных размеров. Длина произвольного стержня 2 распределена по некоторому закону с плотностью р(г), не зависящей от времени I, на полуинтервале (0, 2].

Рис. 1. Схема процесса фильтрации в одномерной решетчатой структуре

Примем следующие допущения.

1. Поток частиц однороден во времени, на период «решета» в единицу времени падает один стержень.

2. Падение частиц равновероятно в любую точку на периоде.

3. Стержень проходит через отверстие, если его центр тяжести попадает в створ отверстия, в противном случае — прилипает к решетке.

Со временем размер проницаемой части «решета» уменьшается из-за налипания стержней.

Задача заключается в отыскании плотностей распределения вероятностей размеров отверстий и стержней на выходе из «решета» С(х, г) и ф(г, г) соответственно.

Исходное уравнение баланса, описывающее процесс зарастания отверстий, имеет следующий вид:

С(х, г + Лг)Лх - С(х, г)Лх = (( (х, г)Дх - 1Ш (х, г)Ах— + о(1),

X (1)

Лг ^ 0, Ах ^ 0, Ух е [0, 1], г > 0.

Здесь левая часть равенства описывает с точностью до бесконечно малых значений приращение за интервал времени (г, г + Дг) вероятности существования отверстий с размерами в диапазоне х...х + Дх; С(х, 0) = 5(х -1) — дельта-функция; 1Ы (х, г) — плотность вероятности образования; 1оШ (х, г) — заполнение отверстий шириной х в момент времени г в результате падения одного стержня; (1, г) = 0, 1оШ (0, г) = 0 У г > 0; х — частота падения стержней на период «решета».

Без ограничения общности результатов исследования для конкретности примем х = 1, а для ненулевых отверстий х > 0. Тогда

2Д 1

1Ы (х, г )Дх = — | С (у, г )Р(у - х)ф, Ух е (0,1), г > 0.

х

2

Здесь Р(у - х) = | р(г)dz представляет собой вероятность суще-

2( У - х)

ствования стержней с длинами, способными образовать из отверстия шириной у отверстие размером х (см. рис. 1).

Множитель Ах /3 — это вероятность попадания центра тяжести стержня на интервал длиной Ах. Коэффициент 2 учитывает налипание стержней с обеих сторон отверстий. Для вычисления уходного члена 1оШ (х, г) рассмотрим три случая возможных соотношений между длиной падающего стержня г и шириной отверстия х:

1) для случая 0 < г < х для изменения размера отверстия центр тяжести стержня должен попасть внутрь интервала на периоде решетки длиной г;

2) для случая 0 < г < 2х суммарная длина интервалов на периоде, попадание на которые приводит к уменьшению размера отверстия, составляет 2г - х, ненулевое отверстие создается с помощью г, а (г - х) / 2 на нуль-отверстие;

3) в случае 2х < г < 2 длина соответствующего интервала равна 2 + х, где 2х идет на создание ненулевого отверстия, а 2 - х — на нуль-отверстия. Таким образом,

Кш (х г)Ах = С (х, г)

Ах,

} гР( г )ёг + { ^ Р( г )ёг + | ^ Р( г)ёг

_ 0 х 2 х

Ух е (0, 1), г > 0.

Нуль-отверстие может образовываться из всех ненулевых отверстиях (0 < х < 1) при попадании центра тяжести стержня с размерами

г < х на интервал длиной г - х в пределах периода «решетки», по-

1 2

этому (0, г )Ах = | С (у, г)| ^^ р( г)ёгёу, У г > 0.

х у

В итоге, устремляя Ах и Аг к нулю, с учетом выражений для (х, г) и 1ш (х, г), а также соответствующих начальных условий

получаем из формулы (1) задачу Коши для системы линейных инте-гродифференциальных уравнений

дС (х, г) = -д( х)С (х, г) + {Р(у-х)С(у, г)ёу, Ух е( 0,1], Уг > 0; (2)

ас

дг

(г) = | С (у, г)} ^ р( г )ёгёу, У г > 0;

С(х, 0) = 5(х -1 + 0), Ух е(0,1]; С0(0) = 0,

(3)

(4)

где

2

Р(^) = - Г р(г)ёг, 0 < w < 1;

2^

д( х) = ■

2 х ,

12 Р( г)ёг + Р( г)ёг + | ^ Р( г )ёг

, Ух е (0,1).

Здесь также С0(г) — вероятность существования на периоде «решета» нуль-отверстия.

Для удобства представим функцию С(х, г) в виде суммы

С(х, г) = С1(х, г) + 5(х -1 + 0)в-д(1)г, (5)

где С1(х, г) — ограниченная составляющая С(х, г).

0

После подстановки суммы (5) в уравнение (2) и преобразований получаем задачу Коши для интегродифференциального уравнения относительно С1(х, г):

х, г) = -(х)С(х, г) + ГР(у-х)С(у, г^у + Р(1 -х)в-*(1)г,

дг х (6)

Ух е( 0, 1], г > 0

с начальным условием

С1(х, 0) = 0, Ух е( 0,1]. (7)

Выражение для плотности условного распределения вероятностей длин стержней на выходе из отверстия ф(г, г) находим аналогично соотношениям для 1Ы (х, г) и 1оШ (х, г), оно имеет вид

1

ф( г, г) = р( г) | Я( у, г )С (у, г ^у, У г е [0, 2), г > 0. (8)

г/2

Здесь Я(у, г) — вероятность того, что стержень длиной г преодолеет отверстие шириной у. Определим ее по формуле

Я(у, г) =

2у - г г . . . 3 ,2 < у <т1п(z, 1);

у, т1п(г, 1) < у < 1.

Решение системы (3)-(4), (6)-(7) и его свойства. Посредством подстановки С1(х, г) = А(х, г)в-д(1)' задача Коши (6)-(7) преобразуется к виду

ВА 1

— (х, г) = -ё(х)А(х, г) +ГР(у -х)А(у, г)dy+Р(1 -х), Вг х (9)

Ух е (0, 1], г > 0; А(х, 0) = 0, Ух е (0,1], (10)

где 0( х) = д(1) - д( х).

Функция А(х, г), Ух е(0,1], г > 0 может быть разложена в степенной ряд

А(х, г) = ^ В] (х)г], (11)

]=1

коэффициенты которого определяются из рекуррентных соотношений

1

(. + 1)£>,.+1(х) = 0(х)Б.(х) +1Р(у-х)Бу(у)ёу, . = 1, 2, ... (12)

х

с начальным условием Б1 (х) = Р(1 - х), Ух е [0, 1].

Свойство 1. Если плотность распределения р(г) непрерывна на отрезке [0, 2], то решение задачи Коши (3)-(4), (6)-(7) существует

д кС

единственно, причем С0(г) е С(0, да), —г- (х, г) е С1(0, 1), У г > 0,

дг

к = 0,1, ....

Доказательство. Достаточно показать, что ряд (12) и соответствующие ему слагаемые ряда, полученные почленным дифференцированием по переменным х и г, сходятся равномерно и являются

непрерывными функциями Уг > 0, х е [0, 1].

Поскольку р(г) непрерывна на отрезке [0, 2], то д(х) е С2(0, 1), Р(х) е С1(0, 1), Б. (х) е С1(0, 1), . = 1, 2, ....

Промажорируем коэффициенты ряда (11) по норме пространств функций С (0, 1) и С1(0, 1) [3]. Запишем неравенства

0 < а(х) < 2тах | х, —, 1 + х)< 1, 3 ^ 2 2 х 2)

-Х < -3, \<2(х)\ < 3(1 - х), 0 < Р(х) < |, Ух е[0,1],

с учетом которых из (2) получаем мажоранты:

. +1

С С0,1) 3 . (. -1)!

, . = 1,2,.... (13)

Посредством дифференцирования (12) по х получаем рекуррентные соотношения для производных коэффициентов Б. (х):

-Б.+, -а -Б.

(. + 1)-хЦ х) = --а (х) Б. (х) + 0( х)-х- ( х) - Р( х) Б. (х) -ах ах ах

1 -Р

-{—(у-х)Б.(у)-у, Ухе[0,1], . = 1,2,...,

из которых следуют оценки

а+1)

+1

йх

с(0, 1)

5 < —

3

В,

1

с(0,1) з

йх

] = 1, 2,....

с(0,1)

На основании этих оценок с учетом (13) и соответствующих начальных условий имеем

Ч сх(0,1)

41Р\с (0,2) + 2.

В

с(0, 2)

3

+ 4 ]

с(°Л) 31 -1(] - 2)!

-, ] = 2,3,....

(14)

(15)

Из соотношений (13)—(15) следует, что ряд (11) и соответствующие ему слагаемые ряда, полученные почленным дифференцированием (11) по переменным х и (, сходятся равномерно, причем для функции С1 (х, г) справедливы оценки

а кс1

&к

<[ 4 1 (Я + 2)3 (

С(0, 2)

+ 8) е

е(1/3-? (1)))

С'(0,1) 4 у

Уг >0, к = 0, 1, .... Учитывая элементарным способом получаемые неравенства

йк+1С 2

й с0-(г) < 2

йг

к+1

йкС1

йг

с (0,1)

+1 дк (Ое-Я(1)г, Уг > 0, к = 0,1,...,

констатируем справедливость утверждения. Свойство 2. Имеют место оценки снизу

С1(х, г) > Р(1 - х)

е-я(х)г - е-Я(1)г

Я (1) - Я (х) Ух е[0, 1], г > 0

> Р(1 - х)

е-Я(1)г - е-Ятахг

- я(1)

где Ятах = Ь 11с(0, 1) •

Доказательство. В силу неотрицательности функции Р(х) на интервале (0, 1) из (9), (10) запишем неравенства

А(х, г) > 0, А(х, г) > Р(х -1)-

М (х )г

-1

0( х)

—x — V

e — e у

Из этого с учетом монотонного убывания функции - по

(У - х)

обоим аргументам следует требуемое доказательство.

Заметим, что С (1, г) = — е~ч(1)'.

1 3г

Свойство 3. Справедливы оценки сверху

С1 (x, t) < 1 3\ 2(1 — x) 1

f [2-1

1 2J-a — x)t

e qmin^ <

(3 / 2) t1/4 —qmintx)t

<V r- --"3Tf e , Vx e[0,1], t > 0,

ЗуЫ (1 — x)3/2 L J

где qmin = m\n q( x)-

2

Доказательство. Поскольку q(x) — qmin > 0, 0 < P(x) < —, Vx e

e[0,1], то C1(x, t) < A+ (x, t)e—qmint, где A+ (x, t) — решение уравне-

CA 2 1 2

ния (x, t) = 31A(y, t)dy + —, Vxe(0,1], t > 0 при начальном

условии Л+ (х, 0) = 0, Ух е (0, 1].

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

w ч 2 3t Т

A+ (x, t) = --L

+ 3 V 2(1 — x) 1

f 12-1

2J-(1 — x)t , Vx e (0,1),

v ' 3

Л. (1, г) =—, Vг > о, + зг

что ввиду свойств функций Бесселя [17] доказывает данное свойство. Свойство 4. Имеет место условие нормировки С0(г) +

1

+1С1(х, г)ёх = 1, Уг > 0.

0

Доказательство. Достаточно показать выполнение тождества —С°(г) + [~С1(х, г)йХ = 0, Уг > 0, что достигается путем подстановки

дг 0 —г

вместо производных соответствующих выражений правых частей системы кинетических уравнений задачи Коши (2)-(4).

Конечно-разностная схема для приближенного решения задачи Коши (6)-(7). На конечных интервалах (0, Т), где Т > 0 — величина порядка постоянной времени процесса, одним из эффективных способов решения задачи Коши (6)-(7) является применение конечно-разностных методов, из которых для краткости изложения рассмотрим схему первого порядка

Б/+1 - В/

= -а/В/ + £Р_кВ/И + Я, / = 1, М;

(16)

к=1

в0 = Р,г/

а = я(1 - ъ); р = р( ъ); ъ = гк;

Я = Ре~%'1; с1 (, г/ ) = (1 + о(1)), т, к ^ 0, г = 0, М; к = -1;

г/ = л; / = 0 1; 1 =

Т

т

Здесь квадратные скобки означают операцию вычисления целой части. Для невязок системы уравнений (16), получающихся после под-

( .у=0~1

становки сеточной функции Вь = |В/} г _ точного решения задачи

(6)-(7), имеют место неравенства

е/ < /1Т + /2к,

(17)

где т > 0; к > 0; / = 1, М; / = 1, 1; коэффициенты /1 =

= Л 8иР

г >0

а 2с1

аг2

; /2 = /3 8ир||сЦс1(0,1); /3 = -((^1 р||с(0,2)) в си-

с (0,1) г >0 3

лу ранее доказанных выше свойств функции С^х, г) являются конечными числами.

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

В'? - В£\ <у 4" - ЯР

(18)

Здесь у =

ешТ -1

ш

ш=НР1с(0 1) -Ятш; сеточная норма определяется

по правилу

В

Ни

= тах

В]

, □ = {(/, /): / = 0, М; / = 0, 1}.

т

Таким образом [3], схема (16) является корректной, а приближенное решение Bh сходится к точному решению, причем из уравнений (17, 18) следуют оценки

Ci (, tj) - Bj< f,x + fh V (i, j) g Q,

где f4 = const fl, f5 = const f2.

Аналогичная конечно-разностная схема первого порядка точности, соответствующая задаче Коши (6)-(7), получается из схемы (16) при

q = Q(1 -z), R = р. Тогда A(, tj) = Bj- (1 + o(x + h)), x, h ^ 0,

i=oTm, j=07.

Значения коэффициентов в мажоранте (17) для этой схемы становятся зависимыми степенным образом от длины интервала интегрирования T.

Из упомянутых схем предпочтительнее вторая, так как в ней свободный член не зависит от времени. Конечно-разностная схема (16) является практически приемлемой при вычислениях на ЭВМ.

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

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

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

циентами q(х) = 1 н--, P(х) = — и параметрами схемы h = 0,025,

63

x = 0,5, представлены на рис. 2. На рис. 3 приведены результаты моделирования той же задачи с исходными данными (см. предыдущий раздел) и параметрами схемы h = 0,025, x = 0.

Ah (х, t)

По оси ординат отложены значения функций a (х, t) = -

IIA

hllh

А(х г)

(приближенное решение), х, г) = ^—(точное решение) (см.

1К11

А (х г)

рис. 2); а3(х, г) = —р— (см. рис. 3), где Ак(х, г) — линейный

\АЛи

сплайн, соответствующий вектору Ак.

1,0 х

Рис. 2. Результаты решения МКР модельной задачи при г, равном 10 (1), 30 (2) и 60 (3): --МКР;---точное решение

1,0 я:

Рис. 3. Графики численного решения МКР при равномерном распределении размеров частиц при г, равном 3 (1), 13 (2) и 60 (3)

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

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

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

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

литература

[1] Alexiades V., Solomon A.D. Mathematical modeling of melting and freezing. processes. Washington DC, Hemisphere Publ. Co, 1993, 323 p.

[2] Baiocchi C. Sur une probleme a frontiere libre traduisant le filtrage de liquids a travers des milieux poreux. Comptes Rendus de l'Académie des Sciences, Série A, 1971, vol. 273, pp. 1215-1217.

[3] Caffarelli L.A. The smoothness of the free surface in a filtration problem. Archive for Rational Mechanics and Analysis, 1976, vol. 63, iss. 1, pp. 77-86.

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

[5] Dimitrienko Yu.I., Bogdanov I.O. Multiscale modeling of filtration liquid binding processes in composite designs at RTM production method. Mathematical Modeling and Computational Methods, 2017, № 2 (14), pp. 32-54.

[6] Димитриенко Ю.И., Шпакова Ю.В., Богданов И.О., Сборщиков С.В. Моделирование процесса многоуровневой фильтрации жидкого связующего в тканевом композите при RTM-методе изготовления. Инженерный журнал: наука и инновации — Engineering Journal: Science and Innovation, 2015, no. 12.

URL: http://engjournal.ru/catalog/msm/pmcm/1454.html (дата обращения 20.10.2017).

[7] Арутюнян Р.В. Моделирование воздействия сильноточного импульса на электрод с учетом нелинейностей характеристик материала и фазовых переходов. Наука и образование, 2016, № 4, с. 139-155. DOI 10.7463/0416.0837747

[8] Арутюнян Р.В. Расчет оттаивания грунтов посредством СВЧ-нагрева на основе метода интегральных уравнений. Известия высших учебных заведений. Электромеханика, 2015, № 6 (542), с. 32-38.

[9] Арутюнян Р.В., Некрасов С.А. Асимптотические и численные методы моделирования диффузной фильтрации. Сибирские электронные математические известия, 2016, т. 13, с. 525-540.

[10] Бондаренко С.С., Боревский Л.В., Гавич И.К. Дзюба А.А., Зекцер И.С. и др. Основы гидрогеологии. Гидрогеодинамика. Новосибирск, Наука, 1983, 242 с.

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

[12] Димитриенко Ю.И., Богданов И.О. Многомасштабное моделирование процессов фильтрации жидкого связующего в композитных конструкциях, из-

готавливаемых методом RTM. Математическое моделирование и численные методы, 2017, № 2, с. 3-27.

[13] Мельникова Ю.С. Математическое моделирование управления нестационарным температурным полем в двухфазных средах. Наука и образование, 2012, № 2, с. 1-19.

[14] Резников Г.Д., Жихарь А.С. Численно-аналитический подход к моделированию переноса частиц в фильтрующем слое. Математическое моделирование, 1995, т. 7, № 6, с. 118-125.

[15] Саваторова В.Л. Моделирование фильтрации жидкости сквозь пористую среду с периодической структурой. Москва, Горная книга, 2010, 62 с.

[16] Саваторова В.Л. Математическое моделирование процессов теплопроводности и фильтрации в неоднородных средах со структурой, близкой к периодической. Дисс. ... д-ра физ.-мат. наук. Москва, 2010, 304 с.

[17] Абрамовиц М., Стиган И., ред. Справочник по специальным функциям. Москва, Наука, 1979, 832 с.

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

Ссылку на эту статью просим оформлять следующим образом:

Арутюнян Р.В. Моделирование стохастических процессов фильтрации в решетчатых системах. Математическое моделирование и численные методы, 2017, № 4, с. 17-30.

Арутюнян Роберт Владимирович — канд. физ.-мат. наук, доцент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор более 60 научных работ. Область научных интересов: математическое моделирование, ортогональные многочлены, теория приближения функций в действительной и комплексной областях, методика преподавания высшей математики в технических университетах. e-mail: rob57@mail.ru

Modeling of stochastic filtration processes in lattice systems

© R.V. Arutyunyan Bauman Moscow State Technical University, Moscow, 105005, Russia

The purpose of the paper was to formulate and study the system of kinetic equations modeling the process of diffusion filtration based on a stochastic approach. Within the research we proved the theorem of existence and uniqueness of the solution with respect to the case of continuous density, obtained the solutions in uniformly convergent and asymptotic series and examined its behavior at infinity. Moreover, we considered the specific cases of density of the Delta-function type and uniform distribution. As a result, the finite-difference scheme for solving the corresponding Cauchy problem on finite time intervals is built and justified. The results of computer simulation are also given.

Keywords: filtration, diffusion, kinetics, stochastic equation, existence, uniqueness, numerical method

references

[1] Alexiades V., Solomon A.D. Mathematical modeling of melting and freezing. processes. Washington DC, Hemisphere Publ. Co, 1993, 323 p.

[2] Baiocchi C. Comptes Rendus de l'Académie des Sciences, Série A, 1971, vol. 273, pp. 1215-1217.

[3] Caffarelli L.A. Archive for Rational Mechanics and Analysis, 1976, vol. 63, iss. 1, pp. 77-86.

[4] Dimitrienko Yu.I., Dimitrienko I.D. European Journal of Mechanics, B/Fluids, 2013, vol. 37, pp. 174-179.

[5] Dimitrienko Yu.I., Bogdanov I.O. Mathematical Modeling and Computational Methods, 2017, no. 2 (14), pp. 32-54.

[6] Dimitrienko Yu.I., Shpakova Yu.V., Bogdanov I.O., Sborschikov S.V. Inzhe-nernyy zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2015, no. 12. Available at:

http://engjournal.ru/catalog/msm/pmcm/1454.html (accessed October 20, 2017).

[7] Arutyunyan R.V. Nauka i obrazovanie — Science and Education, 2016, no. 4, pp. 139-155. DOI 10.7463/0416.0837747

[8] Arutyunyan R.V. Izvestiya vysshikh uchebnykh zavedeniy. Elektromekhanika — Russian Electromechanics, 2015, no. 6 (542), pp. 32-38.

[9] Arutyunyan R.V., Nekrasov S.A. Sibirskie elektronnye matematicheskie izvestiya — Siberian Electronic Mathematical Reports, 2016, vol. 13, pp. 525-540.

[10] Bondarenko S.S., Borevskiy L.V., Gavich I.K., Dzyuba A.A., Zektser I.S. et al. Osnovy gidrogeologii. Gidrogeodinamika [Fundamentals of hydrogeology. Hy-drogeodynamics]. Novosibirsk, Nauka Publ., 1983, 242 p.

[11] Dimitrienko Yu.I., Glazikov M.L. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki — Herald of the Bauman Moscow State Technical University. Series Natural Sciences, 2003, no. 1, pp. 59-71.

[12] Dimitrienko Yu.I., Bogdanov I.O. Matematicheskoe modelirovanie i chislennye metody — Mathematical Modeling and Computational Methods, 2017, no. 2, pp. 3-27.

[13] Melnikova Yu.S. Nauka i obrazovanie — Science and Education, 2012, no. 2, pp. 1-19.

[14] Reznikov G.D., Zhikhar A.S. Matematicheskoe modelirovanie — Mathematical Models and Computer Simulations, 1995, vol. 7, no. 6, pp. 118-125.

[15] Savatorova V.L. Modelirovanie filtratsii zhidkosti skvoz poristuyu sredu s peri-odicheskoy strukturoy [Modeling of fluid filtration through a porous medium with a periodic structure]. Moscow, Gornaya kniga Publ., 2010, 62 p.

[16] Savatorova V.L. Matematicheskoe modelirovanie protsessov teploprovodnosti i filtratsii v neodnorodnykh sredakh so strukturoy, blizkoy k periodicheskoy. Diss. dokt. fiz.-mat. nauk [Mathematical modeling of heat conduction and filtration processes in inhomogeneous media with a structure close to periodic. Dr. phys. and math. sc. diss.]. Moscow, 2010, 304 p.

[17] Abramowitz M., Stegun I.A., ed. Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables. National Bureau of Standards, 1972, 1060 p. [In Russ.: Abramowitz M., Stegun I., ed. Spravochnikpo spetsial-nym funktsiyam. Moscow, Nauka Publ., 1979, 832 p.].

Arutyunyan R.V., Cand. Sc. (Phys.-Math.), Assoc. Professor, Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. Author of over 60 scientific papers. Science research interests: mathematical modeling, orthogonal polynomials, theory of approximation of functions in real and complex fields, methods of teaching higher mathematics at technical universities. e-mail: rob57@mail.ru