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

МАТЕМАТИКА

УДК 510: 53.072:621.1.016.4(03)

ИССЛЕДОВАНИЕ ИНТЕГРОДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ФИЛЬТРАЦИИ

Р. В. Арутюнян

Арутюнян Роберт Владимирович, кандидат физико-математических наук, доцент кафедры математического анализа, Московский технический университет связи и информатики, rob57@mail.ru

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

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

DOI: 10.18500/1816-9791 -2016-16-1-5-12

1. ИНТЕГРОДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ ФИЛЬТРАЦИИ

В статье для исследования процесса зарастания отверстий в решетчатой структуре, играющей роль фильтра, в отличие от [1-4], использован стохастический подход. Рассмотрим одномерную периодическую структуру типа «решета» (рис. 1).

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

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

2

распределена по некоторому закону с плотностью р(г), не зависящей от времени г, на полуинтервале (0, 2]. Примем следующие допущения:

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

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

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

Со временем размер проницаемой части «решета» уменьшается из-за налипания палок. Искомые — плотности распределения вероятностей размеров отверстий и палок на выходе из «решета» С(х,г) и г). Исходное уравнение баланса, описывающее процесс зарастания отверстий, имеет следующий вид:

С(х, г + Дг)Дх - С(х, г)Дх = (/т(х, г)Дх - /ои4(х, г)Дх)--Ъ о(1),

т (1)

Дг ^ о, Дх ^ о, V х е [0,1], г> о,

где левая часть описывает с точностью до бесконечно малых приращение за интервал (г, г + Дг) вероятности существования отверстий с размерами от х до х + Дх; /¿п(х,г) — плотность вероятности образования, а (х, г) — соответственно исчезновение отверстий шириной х в момент г в результате падения палки. Очевидно, /¿п(1,г) = 0, для всех г > 0 (0, г) = 0, С(х, 0) = 5(х — 1) — дельта-функция; т — частота падения палок на период решета. Без ограничения общности примем т = 1. Для ненулевых отверстий (х > 0):

1

2Дх

7т(х, г)Дх = — С (у, г)Р (у — х)ф, V х е (0,1), г > 0,

х

2

где Р(у—х) = / — вероятность существования палок с длинами, способными образовать из

2(у-х)

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

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

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

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

Таким образом,

/ои4(х,г)Дх = С (х,г)

х 2х 2

[ 3 Р(г) ^ + / ^ Р(г) ^ + / ^Р(г) ^

.0 х 2х

Дх,

нуль-отверстие может образовываться из всех ненулевых отверстий (0 < х ^ 1) при попадании центра тяжести палок с размерами г ^ х на интервал длиной г — х в пределах периода «решетки», поэтому

1 2

/т (0,г)Дх = У С (у, г) У р(г) ^у.

ху

Устремляя Дх и Дг к нулю, получаем из (1) с учетом выражений для /¿п(х,г) и (х, г) , а также начальных условий:

1

дС Г

—(х,г) = — ^(х)С(х,г)+ / Р(у — х)С(у,г) V х е (0,1], г> 0, (2)

q(x) =

1 2

de0 (t) = У с (y, t)j ^ p(z) dz dy, V t > 0,

о y

C(x, 0) = 5 (x - 1 + 0), V x G (0,1], Co(0) = 0, 2

P(w) = 2 y p(z) dz, 0 ^ w ^ 1,

2w

x 2x 2

i 2P(z) dz + /^^TpP(z) dz + i P(z) dz , V x G (0,1),

(3)

(4)

.0 х 2х

С0(Ь) — вероятность существования на периоде решета нуль-отверстия. Представим С(х,Ь) в виде

С(х, *) = Сх(ж, Ь) + ¿(х - 1 + 0)е-9(1)*, (5)

где С1 (х,Ь) — ограниченная составляющая С(х,Ь). После подстановки (5) в (2) и преобразований получаем задачу относительно С1(х,^) для любых х £ (0,1], Ь > 0:

dCi

(x, t) = -q(x)Ci (x, t) W P (y - x)Ci (y, t) dy + P (1 - x)e-q(1)t

(6)

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

Ci (x, 0) = 0, V x G (0,1].

(7)

Выражение для плотности условного распределения вероятностей длин палок на выходе из отверстий ^(z,t) получается аналогично соотношениям для Iin(x, t) и Iout(x,t) и имеет вид

p(z,t) = p(z) у R(y,z)C(y,t) dy, V z G [0, 2), t> 0,

z/2

где R(y, z) есть вероятность преодолеть отверстие шириной y палке длиной z:

R(yz) = J(2y - z)/3, z/2 ^ y ^ min(z, 1), [y/3, min(z, 1) < y ^ 1.

Подстановка C1 (x, t) = A(x, t)e-q(1)t преобразует (6)-(7) для любых x G (0,1] и t > 0 к виду

i

dA Г

— (x, t) = -Q(x)A(x, t) + P (y - x)A(y, t) dy + P (1 - x),

x

A(x, 0) = 0, V x G (0,1],

(8)

(9) (10)

где ф(х) = д(1) — д(х). Функция А(х,Ь) для любых х £ (0,1] и Ь ^ 0 может быть разложена в степенной ряд:

A(x,t) = £ Dj (x)tj,

(11)

j=1

коэффициенты которого определяются из соотношений для j = 1, 2,...,

(j + 1)Dj+1 (x)= Q(x)Dj (x)+ / P (y - x)Dj (y) dy, D1 (x) = P (1 - x), V x G [0,1]. (12)

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

д к С

Co(t) G Cœ(0, œ),

■ (x,t) G C1 (0,1), V t> 0, k = 0,1,... .

1

1

1

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

Промажорируем коэффициенты ряда (11):

/х 2 /х 3 х \ 0 < о(х) < - шах -, —, 1 + - < 1,

чу } 3 2х' 2) '

2

0 ^ Р(х) ^ з,

¿х

<

1О(х)| ^ 3(1 — х),

V х е [0,1],

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

311с(0,1)

3 + 1

з3 (3 — 1)!'

3 = 1, 2,

(13)

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

1

(3 + 1) (х) = — ¿х(х)^з (х) + О(х)(х) — Р(х)Ю (х) — / ¿Р (у — х)Ю (у) ¿у,

V х е [0,1], 3 = 1, 2,...,

откуда следуют оценки

(3 + 1)

3+1

¿х

С(0,1)

5 1

^ 3 ||ю3||С(0,1) +

3

¿х

С(0,1)

3 = 1, 2,...,

на основании которых с учетом (13) и начальных условий имеем:

4||р|С(0,2) + 2

ю

1Нс!(0,1) =

3

п ц _ ||р|с(0,2) +43'

ю Ус1 (0,1) = 33-1 (3 — 2)! ,

3 = 2, 3, ...

(14)

(15)

Из (13))-(15) следует, что ряд (11) и ему соответствующие, полученные почленным дифференцированием (11) по х и г, сходятся равномерно, причем при всех г ^ 0, к = 0,1,...

д к С1

дгк

Так как

^^(г) < 2

1

<

С1(0,1)

¿к С1

(3) (Л + 2)3(г + ||р||с(0,2) +8)3в(1/3-я(1))* .

¿гк

С(0,1)

¿гк+1 ^ ^ 3 то утверждение доказано.

Свойство 2. Имеют место оценки снизу

+ 15к(г)в-я(1)*, V г ^ 0, к = 0,1,..., 3

□

в-я(х)г _ е—д(1)* е—^(1)* _ в Яшах *

С1 (х,г) ^ Р(1 — х) —- ^ Р(1 — х) ■

5(1) — 5(х)

5тах — 5(1)

V х е [0, 1], г ^ 0, 5тах = |Ы|с(0,1). Доказательство. Так как Р(х) ^ 0 V х е (0,1), то из (9)-(10): А(х,г) ^ 0, А(х,г) ^

^ Р(х — 1)-

(х)* — 1

. С учетом монотонности

е х — в у (у — х)

ство.

Заметим, что С1(1,г) = Зтв -я(1)* .

по обоим аргументам, следует доказатель-

□

1

3

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

С (х Ь) < 2 л /_—_/Л 2а/2М - х)Ь | е-дт1п* < (3/2)3/4_—_е-дт1п2(1-х)^

С1(х ,Ь) < 3\/2(1 — х) 'М ^3(1 х)М6 < (1 — х)з/2 е '

V х £ [0,1], Ь ^ 0, дт1П = :го.1г1 д(х). Доказательство. Поскольку д(х) — Отт ^ 0, 0 < Р(х) < 2/3, V х £ [0, 1], то С1(х, Ь) <

1

-А+ , 2 [ , 2

< А+(х,Ь)е , где ——— (х,Ь) = - А(у,Ь) ¿у +—, V х £ (0,1], Ь ^ 0, при начальном условии

дЬ 3 3

х

А+ (х, 0) = 0, V х £ (0,1]. Аналитическое решение данной задачи Коши имеет вид

А+(х,Ь) = 2у2(13—х /1 (2^3(1 — х)^, V х £ (0,1), А+(1,Ь) = |, V Ь > 0,

что вследствие свойств функций Бесселя [5] доказывает данное свойство. □

Свойство 4. Имеет место условие нормировки:

1

Со(Ь) + У С1 (х,Ь) ¿х = 1, V Ь ^ 0. 0

Для доказательства достаточно доказать

1

(Ь)+ / дС1 (х,Ь)¿х = 0, V Ь ^ 0, дЬ 7 дЬ ' '

0

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

2. ИССЛЕДОВАНИЕ АСИМПТОТИЧЕСКИХ СВОЙСТВ (6)-(7)

Исследуем важный случай, когда д(х) < 0, V х £ (0,1), вследствие чего на полуинтервале [0,1) существует не более одной точки с абсциссой х*, для которой выполняется равенство д(х*) = д(1), если же д(х*) ^ о(1), V х £ [0,1), то положим х* = 1. Будем искать коэффициенты асимптотического разложения А(х,Ь) на (х*, 1] в ряд Лорана (предполагаем, что х* < 1):

А(х,Ь)=^ аэ (х)Ь—, г> 0, Ь ^ го (16)

3 = -г

Подставляя (16) в (9), находим соотношения его коэффициентов:

1

(1 — 3К-1(х) = ^(хК(х) + У Р(У — хК(У) ^ 3 = —^ •••, —1,1,•••;

х

о^-(х) = 0, (—3) = г, г + 1,--- ;

1

0-1 (х) = ^(х)о0(х) ^ У Р(у — х)о0(у) ¿у + Р(1 — х), ^(х) > 0, V х £ (х*, 1]-

х

Для определения степени г главного члена асимптотики (16) сделаем оценку снизу и сверху для функции С1 (х,Ь). Поскольку Р(х) < 2/3, то С1(х,Ь) < А0(х,Ь), где А0(х,Ь) есть решение задачи Коши для каждого х £ (х*, 1] и Ь > 0:

1

о л О Г

(х, Ь) + о(х)А0(х, Ь) = 2 А0(у, Ь) ¿у + Р(1 — х)в^(1)*, (17)

Ао(ж,0) = 0, V х е (ж*, 1]. (18)

1

Рассмотрим эд0(ж,Ь) = / А0(у,Ь) ¿у при произвольных ж е (ж*, 1] и Ь > 0. Точное решение (17), (18):

X

1 2 У

^о(ж,*)= /( + рр1(1)7+) ( )) е3 Х ^^ ¿у, V ж е (ж*, 1], (19)

} (5 + „(1))(в + „(у))

X

где г&0(ж,з) — образ функции (ж, Ь), 5 — параметр преобразования Лапласа.

Из (19) получаем, что гЪ0(ж, 5 — „(1)) = £(ж)з-п(1)-1 (1 + о(1)) для любого ж е (ж*, 1], 5 ^ го, где

£(ж) = [„(ж) — д(1)]п(1)е X ч(-)-ч(1) *, п(ж) =----„'(ж) = V ж е (ж, 1],

3д'(ж) аж

откуда

*„0(ж ь) = Ьп(1)е-ч(1)^

Г(п(1) + 1)

и с учетом

^0(ж, Ь) = Ьга(1)е_ч(1)* ^(ж) (1 + о(1)), V ж е (ж*, 1), Ь ^го

А(ж,Ь) = д'(1)Ьп(1)е-ч(1)*Г(П|1щ(1 + о(1)), V ж е (ж*, 1), Ь ^ го. (20)

Из непрерывности р(г) и условия д"(ж) < 0 при любом ж е (0,1) следует, что „'(ж) < 0, поэтом п(1) конечно и ж* ^ 1. В силу справедливости неравенства Р(ж—у) > Р(1—ж) для всех у е (ж, 1), ж е (0,1) А1 (ж, Ь), являющаяся решением задачи Коши:

1

д А С

-1 (ж, Ь) + д(ж)А1(ж,Ь)= Р(1 — ж) / А1 (у, Ь) ¿у + Р(1 — ж)е-ч(1)*, V ж е (ж, 1], Ь> 0, дЬ ]

X

А1 (ж, 0) = 0, V ж е (ж*, 1],

есть оценка снизу для С1 (ж, Ь) . Асимптотика для функции А1(ж,Ь) находится также при помощи преобразования Лапласа и имеет вид

А1(ж,Ь) = — „(1)Р-У^'1'-' ^*_""_*) ^(1 + о(1)), (21)

V ж е (ж*, 1), Ь ^ го.

Из (20) и (21) следует, что максимальный порядок степеней Ь в слагаемых асимптотического разложения (16) г = п(1), а также выводится краевое условие для коэффициента апс) (ж) главного члена асимптотики

^ «п(1)(ж) _ [ — „(1)]п(1)

х^1-0 (1 — ж)п(1)-1 Г(п(1)) '

(22)

являющееся условием однозначности для соответствующего уравнения системы соотношений для коэффициентов а,](ж), з = — п(1), — п(1) + 1,...

3. КОНЕЧНО-РАЗНОСТНАЯ СХЕМА ДЛЯ РЕШЕНИЯ ЗАДАЧИ КОШИ (6)-(7)

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

В+1 - В

= — а В] + 5] В] ^ + Р], г = 1,М, (23)

к=1

В] = Р0 Ь] е_Ч0 ^, а = „(1 — ^), Р = Р (^), ^ = гй, Р] = Р е_Чс^,

C B^-i (1 + o(1)), т, h ^ 0, i = 0,M,

j = jT j = 0,7, J =[t/t ],

h = M'

где квадратные скобки означают вычисление целой части. Для невязок (23), получающихся после подстановки сеточной функции В^ = {В^}^=00М точного решения задачи (6)-(7), имеют место неравенства

ej ^ /it + /2h,

d 2 C

dt2

C(0,1)

(24)

; f2 = f3 SUp ||Ci Ус1 (0,1);

i>0

где т ^ 0, h ^ 0, i = 1, M, j = 1, J; коэффициенты /1 = /3 sup 2

/3 = -(1 + 2||p||C(0,2)) в силу ранее доказанных в пп. 1 и 2 свойств C1 (x,t) являются конечными. 3

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

IBh1' - Bh

(2) II < Y|R(1) Т?(2Ь

h ^ 7||Rh - Rh I|h,

„шТ

1

Y =

w = IIPI

w

C (0,1) — qmin

(25)

сеточная норма равна: ||Bh||h = max |B||, О = {(i, j) : i = 0, M, j = 0, J}.

Таким образом [6], схема (23) является корректной, а приближенное решение Bh сходится к точному, причем из (24)-(25) следуют оценки:

C1(z<, tj) — Bm-i| ^ /4T + /5h, V (i, j) e О,

где /4 = const /1, /5 = const /2. Аналогичная схема первого порядка точности для задачи (6)-(7)

получается из (23) при q = Q(1 — z), Rj = , тогда A(z, tj) = Bm^^! + о(т + h)), т, h

0,

г = 0, М, 3 = 0, J. Значения коэффициентов в мажоранте (24) для этой схемы становятся зависящими степенным образом от Т. Из упомянутых схем предпочтительнее вторая, так как в ней свободный член не зависит от времени. Результаты моделирования даны на рис. 2 (задача с коэффициентами д(х) = 1 + (1 — х)/6, Р(х) = 2/3, и квазиполиномиальным решением, параметры к = 0^025, т = 0^5) и рис. 3 (данные п. 2.2, параметры к = 0^025, т = 0^(3)).

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

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

По оси ординат откладывались: на рис. 2 — a1(x,t) = t) (приближенное решение),

IA

a2(x,t) = .. (точное решение); на рис. 3 ||Ah У

соответствующий вектору Ah.

h| h

аз(х, t) = At), Ah(x,t)

линейный сплайн,

h| h

ЗАКЛЮЧЕНИЕ

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

Библиографический список

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

2. Гавич И. К., Зекцер И. С., Ковалевский В. С., Язвин Л. С., Пиннекер Е. В., Бондаренко С. С., Боревский Л. В., Дзюба А. А. Основы гидрогеологии. Гидрогеодинамика. Новосибирск : Наука, Сиб. отд-ние, 1983. 241 с.

3. Гуревич А. Е. Практическое руководство по изучению движения подземных вод при поисках полезных ископаемых. Л. : Недра, 1980. 216 с.

4. Колесников А. В. Математическое моделирование фильтрации жидкости в неоднородных и периодических пористых телах методом однородно-анизотропного эквивалентирования : автореф. дис. ... канд. техн. наук / Северо-Кавказский федеральный университет. Ставрополь, 2014.

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

6. Треногин В. А. Функциональный анализ. М. : Наука, 1980. 496 с.

Stochastic Simulation of Diffusion Filtering R. V. Arutyunyan

Arutyunyan Robert Vladimirovich, Moscow Technical University of Communications and Informatics, 8a, Aviamotornaya st., Moscow, Russia, 111024, rob57@mail.ru

Formulated and investigated is the system of kinetic equations describing the process of diffusion filtering based on a stochastic approach. The theorem of existence and uniqueness of the solution for the case of a continuous density is prove. We obtain the representation of solution in the form of a uniformly convergent and asymptotic series, and explore the nature of its behavior at infinity. The concrete particular cases such as the density of the delta function and a uniform distribution are considered. The finite-difference scheme for the solution of the corresponding Cauchy problem on finite intervals of time is constructed and justified. The results of computer simulation are given.

Key words: filtration, diffusion, kinetics, stochastic equation, existence, uniqueness, numerical method.

References

1. Reznikov G. D., Zhikhar' A. S. Chislenno-analiti-cheskii podkhod k modelirovaniiu perenosa chastits v fil'truiushchem sloe [Numerical and analytical approach to the modeling of the transport of particles in the filter bed]. Matematicheskoe mod-elirovanie [Mathematical modeling], 1995, vol. 7, no. 6, pp. 118-125 (in Russian).

2. Gavich I. K., Zektser I. S., Kovalevskii V. S., Iazvin L. S., Pinneker E. V., Bondarenko S. S., Borevskii L. V., Dziuba A. A. Osnovy gidroge-ologii. Gidrogeodinamika [Fundamentals of hydro-geology. Hidrogeodinamika]. Novosibirsk, Nauka, 1983, 241 p. (in Russian).

3. Gurevich A. E. Prakticheskoe rukovodstvo po izu-cheniiu dvizheniia podzemnykh vod pri poiskakh poleznykh iskopaemykh [A Practical Guide to the study of groundwater movement in the search for

minerals]. Leningrad, Nedra, 1980, 216 p. (in Russian).

4. Kolesnikov A. V. Matematicheskoe modelirovanie fil'tratsii zhidkosti v neodnorodnykh i periodich-eskikh poristykh telakh metodom odnorodno-anizotropnogo ekvivalentirovaniia [Mathematical modeling of fluid flow in heterogeneous and periodic porous bodies by homogeneous anisotropic equivalenting] : avtoref. dis. ... kand. tekhn. nauk / Severo-Kavkazskii federal'nyi universitet. Stavropol', 2014 (in Russian).

5. Spravochnik po spetsial'nym funktsiiam [Handbook of special functions] / eds. M. Abramovitsa, I. Stigan. Moscow, Nauka, 1979, 832 p. (in Russian).

6. Trenogin V. A. Funktsional'nyi analiz [Function analysis]. Moscow, Nauka, 1980, 496 p. (in Russian).