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

УДК 517.958:531.72, 517.958:539.3(4)

ЧИСЛЕННОЕ РЕШЕНИЕ УСРЕДНЕННОЙ МОДЕЛИ СОВМЕСТНОГО ДВИЖЕНИЯ ЖИДКОСТЕЙ РАЗЛИЧНОЙ ПЛОТНОСТИ В ПОРИСТЫХ СРЕДАХ 13)

О.В. Гальцев

Белгородский государственный национальный исследовательский университет,

ул. Студенческая, 14, Белгород, 308007, Россия, e-mail: oleggaltsev@yandex.ru

Аннотация. В работе исследуется задача со свободной границей, описывающая на макроскопическом уровне совместную фильтрацию двух несмешивающихся несжимаемых жидкостей. Основными уравнениями математической модели на микроскопическом уровне являются стационарная система уравнений Стокса для несжимаемой неоднородной вязкой жидкости в порах и стационарная система уравнений Ламе для несжимаемого упругого скелета, дополненные граничными условиями непрерывности перемещений и нормальных напряжений на общей границе «твердый скелет - жидкость» в порах, и транспортного уравнения для неизвестной плотности жидкости. Методом двухмасштабного разложения выводится макроскопическая математическая модель совместного движения жидкости и упругого скелета, и обсуждается численная аппроксимация усредненной модели для специальной структуры порового пространства.

Ключевые слова: задача Маскета, задача со свободной границей, фильтрация жидкости, усреднение.

1. Введение

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

Среди математических моделей, описывающих такое движение широко известна задача Маскета [1]. Эта феноменологическая модель описывает фильтрацию двух несме-шивающихся несжимаемых жидкостей различной плотности, разделенных некоторой подвижной (свободной) границей. Движение первой жидкости под воздействием силы тяжести в области Q+(t) с постоянной плотностью р+ описывается системой уравнений Дарси

v+ = -№р+ + , V- v+ = 0, х е П+(Ь), (1)

13 Работа выполнена в рамках ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы (госконтракт №02.740.11.0613).

для микроскопической скорости "У+ и давления р+ жидкости. Единичный вектор ^ совпадает с направлением силы тяжести, а заданная постоянная к называется проницаемостью грунта.

Соответственно, движение второй жидкости под воздействием силы тяжести в области П-(£) с постоянной плотностью р- описывается той же самой системой уравнений Дарси

V- = -Ур- + р-^, V- V- = 0, х е П-СО (2)

для микроскопической скорости V- и давления р-. На общей свободной границе Г(£) = дП+(£) Р| дП-(¿) давления и нормальные скорости непрерывны

р+ = р-, х е г(г), (3)

v+ ■ п = V- ■ п = к, х е г(г). (4)

Условие (4) означает, что граница Г(£) есть материальная поверхность. То есть во время движения она состоит из одних и тех же материальных точек. Этот факт позволяет сформулировать обобщенную постановку задачи Маскета. Определим давление р неоднородной жидкости как

р = р+ если х е П+(£), р = р- если х е П-(£),

плотность р/ неоднородной жидкости как

р/ = р+ если х е П+(£), р/ = р- если х е П-(£),

и скорость V неоднородной жидкости как

V = v+ если х е П+(£), V = V- если х е П-(¿).

Неизвестные функции V, р, и р/ будут решениями системы уравнений фильтрации Дар-си

к

« = —-Ур + р/ V • V = 0, ж € п, I > 0, (5)

и уравнения переноса

dpf _ dpf

+ v •V pf = 0, x G Q, t > 0. (6)

^ д*

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

V ■ Ур = V ■ (vp), V -Vр/ = V ■ (vр/).

Проблема дополняется однородным граничным условием

V ■ п = 0, х е 5 = дП, ¿> 0, (7)

где п есть нормальный вектор к границе Б, и начальным условием

р/ (х, 0) = Р0(х), х е П (8)

с разрывными начальными значениями

р0(х) = р+, если х е П+ и р0(х) = р-, если х е П-.

В настоящее время только в нескольких работах рассматривается классическая разрешимость задачи Маскета локальной по времени или глобальной по времени, но вблизи известного решения (см. [2], [3], [4]). Какие-либо результаты обобщенного решения автору неизвестны.

В работе А.М. Мейрманова [5], был предложен новый подход для математического моделирования указанного типа фильтрации. А именно, физический процесс, в первую очередь, описывается на микроскопическом уровне системой уравнений, состоящей из уравнения Стокса для движения вязкой несжимаемой жидкости в порах, уравнения Ламэ для движения упругого несжимаемого скелета

д'wе

V • (Х£^оЩх, —) + (1 - оВ{х, ю£) - р£1) + р£ Р = 0, (9)

V- wе = 0, (10)

we(x,t) = 0 при х е Б = дП, (11)

X£w£(x, 0) при х е П, (12)

/ р£ = 0 (13)

ип

и транспортного уравнения для плотности

с1р£ _ др£ ди1£ £ £ £

в области П, где П - единичный квадрат (0 < ж < 1) х (0 < ж2 < 1), П = П/ и и П^, П/ поровое пространство (область, занятая жидкостью), П - твердый скелет и Бе есть поверхность «твердый скелет - поровое пространство», wе(x,t) = w|(x,t)) -

перемещение сплошной среды, ре(х,^ - давление в сплошной среде, D(ж,дw/дt) есть симметричная часть градиента дw/дt (тензор напряжения), I - единичная матрица, Xе - характеристическая функция порового пространства, р/ - плотность жидкости, соотнесенная к плотности воды р0, р8 есть плотность упругого скелета, р0 - безразмерная вязкость жидкости, А0 - безразмерная постоянная Ламэ, Р - известная вектор-функция.

Далее эта система путем усреднения аппроксимируется задачей Маскета для вязко-упругой фильтрации

V- Р — ^ + р Р = 0, (15)

V ■ w = 0, р = тр/ + (1 — т)р5,

/ дм \ Г*

F = ^^osЯ1:B^x,—J+Xo^^2^Щx,w) + J УХ3{1 - т) : Щх, ю{х, т))йт , (16)

где N1, N2 и N3 - тензора четвертого порядка, которые описаны ниже.

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

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

1 = 50

1 = 860

1 = 1200

1 = 2600

1 = 4800

Рис.1 Сравнение результатов численного моделирования усредненной задачи движения

жидкости

в поровом пространстве (снизу) и задачи совместного движения двух жидкостей в упругом

скелете на микроскопическом уровне (сверху).

§2. Вывод усредненных уравнений

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

Математическая модель, описывающая совместное движение вязкой несжимаемой жидкости и упругого несжимаемого скелета на микроскопическом уровне, имеет вид (9)-(13). Эта система дополняется транспортным уравнением (14). Будем искать решение задачи в виде

дм£ дм д^

^£(ж,*) = + е^(ж, ж/е,*) + ... ,

р£ = р(ж,*), р£ = Р(ж, ж/е,*) + ... ,

где точками обозначены члены более высокого порядка малости.

При выводе усредненных уравнений нам понадобятся хорошо известное свойство

/ ^(ж, ж/е,*) ^ж ^ ( / <^(ж, у, ¿) ^у) ^ж ./Пт ./пт -/г

для всякой гладкой 1-периодической по переменной у функции <^(ж, у) ( [5], [8]). Здесь и далее ^ж означает дифференциал лебеговой меры в М3.

Переходить к пределу в интегральном тождестве (1) при е д 0 будем с двумя различными видами пробных функций. Прежде всего с пробной функцией ^ = ^(ж, *), которая является произвольной гладкой вектор-функцией, равной нулю на границе Б. Умножим (9) на функцию ^(ж,*) и проинтегрируем по частям по области Пу

/ Р : Р(ж, ^(ж,*)) ¿ж^ = р£Р ■ ^(ж,*) ¿ж^. (17)

Пт Пт

Далее подставим и р£ в Р, распишем каждое слагаемое.

О (ж,— ) =В ж,—

= + О (у, ^(ж, ж/е, *) ) + ... , (18)

D(x, w£) = D(x, w(x, t) + eW (x, x/e, t)) =

= D(x, w(x,t)) + eD (x, W(x, x/e,t)) + D(y, W(x, x/e,t)) + ... . (19)

Первое слагаемое в (9) при e g 0 примет вид Г ( dw£ \

J j : Р(ж, <р(х, t)) dxdt =

= jf х(х/е)(ло + ) : Р(ж, <p(x, t)) dxdt , (20)

f ( dw \

lim / I ж, -7— (x,t) ) : B>(x,(p(x,t)) dxdt

£go JnT V dt /

dw

mjUoD ( ж, —— (ж,i) ) : , (21)

'nT V dt

Г ( dW х \

lim / роЬ —:Ю>(x,<p(x,t)) dxdt = s\°JnT V dt e /

= J £ x(y)ß0° (у, : Ю)(ж' *)) dxdVdt =

где m = J x(y)dy есть пористость среды и eD(x, dW/dt(x, x/e, t)) \ 0 при e \ 0. Второе слагаемое при e \ 0 будет

/ (1 - x(x/e))A0D(x, w£) : D(x, p(x,t)) dxdt =

JnT

= / (1 - x(x/e))A0[D(x, w(x,t)) + D(y, W (x, y,t))] : D(x, p(x,t)) dxdt, (23)

JnT

lim / A0(1 - x(x/e))D(x, w(x,t)) : D(x, ^(x,t)) dxdt = sgoJnT

= / A0(1 - x(y))D(x, w(x,t)): D(x, p(x,t)) dxdydt =

JnT Jy

= (1 - m)A0D(x, w(x,t)): D(x, p(x,t)) dxdt, (24)

JnT

lim / A0(1 - x(x/e))D(y, W(x, x/e,t)) : D(x, p(x,t)) dxdt =

eg0JnT

= i i (1 - x(y))A0D(y, W(x, y, t)) : D(x, p(x, t)) dxdydt

./nT ./ Y

= ^ А^ ^(1 - х(у))0(у, ^(х, у,*))^ : Ю>(®, р(®,*)) (25)

где (1 — т) = J (1 — х(у))^У есть пористость среды и еО(ж, ^(х, ¿)) \ 0 при е \ 0. Третье слагаемое при е д 0

Pdy) div^dxdt = pdiv^dxdt, (26)

П v JY y JnT

так как (Р)у = р.

Наконец, четвертое слагаемое

/ р£Р ■ ^ ^ж^* = р(ж,*)Р ■ ^ ^ж^*

'пт ./пт

Собрав все слагаемые вместе, получим

J тц0В : Р(ж, </?(ж, ;£)) сМ +

+ / (1 - т)Л0О(ж, ад(ж,*)): »(ж, р(ж,*)) ^ж^ +

+ / М / (1 - хЫЭ(у, ^(ж, у,*))^у) : »(ж, р(ж,*)) М -./пт -ЗУ

- / р(ж,*) а1у ^(ж,*) ^ж^ + / р(ж,*)Р ■ ^(ж,*) = 0 . (27)

пт пт

Снова интегрируя (3), приходим к макроскопическому уравнению сохранения импульса

V

/ дад N / / д^ \\ /

роСтВ ( ж,—(ж,^ 1 ( у,—(х,у,г) ) ^ ) + Ло^1-?в)0(ж,ю(ж,^) +

+ <»(у, ^(ж, у,*))>у) - р1] + рР = 0 . (28) Проделаем аналогичные действия с уравнением неразрывности

/ div м£п ^ж^* = div(w£n) ^ж^* - / ■ Уп ^ж^ = ./пт ./пт .Упт

= (м£ ■ - / ■ Уп ^ж^ = 0 , (29)

•/ дПт ./Пт

где п = п(ж, *) - произвольная гладкая функция. Первое слагаемое равно нулю в силу граничного условия (11) и, подставив в (13) = м(ж,*) + е^(ж, у,*), получим

Иш / (м(ж,*) + е^(ж,у,*)) -Уп^ж^* = м(ж,*) -Уп^ж^*. (30)

£до ,/пт ,/пт

Снова интегрируя (14) получим макроскопическое уравнению неразрывности

У- м = 0 . (31)

В качестве пробной функции в (3) возьмем функцию вида ^ = еЛ,(х, ¿)^0(х/е) и проведем аналогичные рассуждения. В результате получим следующее микроскопическое уравнение сохранения импульса

V,,

х и О у

дШ

~дГ

+ Ъ) + Ао (1 — х)О(у, Ш) — РI

0

(32)

/ й \ 2 Z(ж,í) = р0Ю) (х,-^) - ЛоЮ>(а5, го) = ^ У х (0,Т) .

^ ' 1,3=1

Аналогично, уравнение неразрывности на микроскопическом уровне примет вид

V, ■ Ш = 0, у е У, г е (о,т). (33)

Для нахождения тензоров N1, N2 и N3 в (15)-(16) необходимо решить систему (16)-(33), найдя О(у, дШ/дг) и О(у, Ш) как операторы от О(х, д^/дг) и О(х, ад) и, заменив эти выражение в (4).

Пусть {Ш (*'3)(у,г), Р(^'3)(у,г)} и (Ш0^'3)(у), Р(?'з)(у)} г,^ = 1, 2 -решение периодической задачи, где в первую очередь в У^ решается

V,

(х (Ро О(у, ш(3))+ 1(3) — Р0(3) I)) =0

х(у)Ш ^(у^у = 0

(34)

'у

Систему (18) можно представить в виде

Д, Ш03 — V,Р0'3 = 0, у е Уf ,

0 0, у е Yf ,

(¿3)

где

(ро О(у, Ш0гз)) + 1(3) — Р0(3) I) ■ п = 0 , 7 = дУf \дУ

1

Г3 = - (вг О е,- + е,- О ег) . Затем решается система уравнений во всей элементарной ячейке У

V,

хР0 О у,

д\¥{гЛ' дЬ

+ А0 (1 — х)О(у, Ш(3)) — Р(3) ^ =0 ,

V, ■ Ш(3) = 0 , х(у)Ш(3)(у, 0) = Ш0гз)(у) . J

(35)

Тогда

2 /•*

Ш(х,у,г) = V / %(х,г)Ш(гз)(у,г — т)^т

¿,3=1

0

0

и

2

Р (ж, у,*) = х(у)£ ^ (ж,*)Ро(7)(у)+^ / %(ж,т)Р(7)(у,* - т ¿,7=1 ¿,7=1 ^

С учетом вышесказанного уравнение (4) примет вид (16), где тензора

N1 = ро т 1 + ро (Ао)у/ , N2 = Ло (1 - т)1 + Ло (Ао>Ув + ро (А1(у, 0))У/ '5А1

лм,

т3(г) = ро (-^(зМ)) + ло {Ыу^))га ,

где

2

Ао(у) = Ро »(у, ^?7)(у)) 0 1(7), (37)

¿,7=1

2

(36)

/ / д^ (¿7) \ \

211(3/,*) = (й)Ю> -Л0О(у,^)(у,*))) (38)

¿,7=1 V /

2

1 = ^ 1г3 0 17 • ¿,7=1

Аналогично поступим с уравнением переноса (14)

С д£ дм£ /*

J (р£^ + р£^- • УС) с1хсЫ = - у ро £(ж, 0) с?ж , (39)

где £ - произвольная гладкая функция, равная нулю при * = Т. Переходя к пределу при е \ 0 получим

С д£ с д£ с д£

Ит / р£^с1хсИ = / р^ с1хсЫ = / р(х,1)^с1хсЫ , (40)

( дм£

Г (дм д^ \

= Ит / р£ — (ж, г) + - (ж, у, • УС с1хсИ = £до ./пт ^ д* д* /

Г дм

= / р— (ж,*) • У£сМ , (41)

Иш / ро£(ж, 0) ^ж = ро(ж)£(ж, 0) ^ж (42)

£до Уп Уп

или записывая все слагаемые вместе

Г д£ Г дм Г

/ р(х^)^с1хсИ+ I р—(ж,£) • У£сЫ = - / ро(ж)£(ж,0)сЫ . (43)

После реинтегрирования уравнение (12) переходим к усредненному уравнению переноса

^р др дад „ . „ „.

§3. Вычислительный алгоритм

В качестве численного метода решения системы уравнений (15) - (16), был выбран один из наиболее универсальных и эффективных разностных методов - модифицированный метод конечных разностей.

Построение решения с использованием этого метода следующее:

1 этап. Область непрерывного изменения аргумента заменяется некоторым конечным дискретным множеством точек (узлами).

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

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

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

Типичное эйлерово представление, состоящее из множества (г + 1) точек (х) описывает одномерную пространственную область < х < х+1). Размеры ячеек определяются разностями (х — ж^-1).

Построение подходящей сетки в двух или трех измерениях часто является трудной задачей. Часть трудностей возникает из-за геометрической сложности расчетной области, часть же трудностей возникает из-за того, что сама структура волновых процессов в многомерном случае более сложна, чем в одномерном. Существует 6 типов многомерных эйлеровых сеток: регулярная прямоугольная сетка; прямоугольная адаптивная сетка; регулярная треугольная сетка; треугольная сетка с локальным измельчением; прямоугольная сетка с вложенной мелкой сеткой; ортогональная сетка с вложенной мелкой сеткой.

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

Для численного решения задач (15)-(16), (18)-(17) была применена двухуровневая сеточная модель, в основу которой положены следующие принципы:

• сетка является равномерной по всем осям координат. Данное условие позволяет в дальнейшем упростить аппроксимирующие выражения;

• все граничные узлы являются регулярными по всем осям координат, т.е. каждый граничный узел совпадает с узлом сетки, а не ложится между двумя соседними узлами сетки. Данное условие зависит от величины шага сетки. Чем больше дискретность сетки, т.е. чем она «мельче», тем более точным будет геометрическая границ;

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

На первом уровне хранятся значения искомой сеточной функции на трех временных слоях. На втором - признаки расположения узла, а именно:

• основные геометрические признаки расположения узла (внутри, вне или на границе расчетной области); если узел лежит на границе, то конкретизируется тип границы;

• вспомогательные геометрические признаки расположения границы относительно центра расчетной области.

Для построения разностных схем данной сеточной модели введем в области П = [0 < х < х [0 < х2 < N^[0 < г < Т] следующую сетку (Рис. 2):

х(/+1) = гйь й( > 0; г = 0,1,...,ЛТ(;

ж23+1) = ^2, ^ > 0; з = 0,

где Л,(, Л,2 - размер сетки, N1, N - количество ячеек сетки, соответственно, в направлении Х( и х2.

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

_ /,;.,- . /,; /.;,- Ь) „ /,;,- . /,; 2/.;,-;. - /■;,• /,; I (х) 2Ь . ./ (х; - ь;2

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

Рис. 2. Шаблон разностной сетки.

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

Этап 1: Находим значения №к = ^К02'Ы), решая систему уравнений (18) для

элементарной ячейки У f. Основная трудность при численном решении связана с вычислением поля давления. Первый значительный успех в преодолении этих трудностей был достигнут благодаря идее искусственной сжимаемости [9]. Эта идея очень важна для вычисления неизвестного давления Р0 из уравнения.

др к1

^ + = 0. (45)

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

Так как к = 1, 2 и / = 1, 2, то первое уравнение можно записать в удобном виде для дальнейшего представления в виде разностной схемы:

д2ж1'ы 1 /д2ж1'ы д2ж2'ы\ дРо

_|__. и _|___0

ду\ 2 у ду\ дугду-2 ) дуг 1 /д2ж1'ы . д2Ж2'ы\ д2Ж2'ы дРо

2 + ., .', 1 + /'••

2 V ду| ду2 ду2

дЖ

дУ1

+

дЖ

дУ2

0 •

Разностные схемы для вычисления и Р^1 на следующем временном слое

будут иметь вид

Ж

1М

Ш1М -I--I--I- __— < ры — ры \

уу0Л+1^ ' УУ0Л~1^ ' ^ СМ,7+1 ' уу 04^-1 \г0Л+1^ г0Л-1^>

лм

г1М

к

Ж

Т/Г/ ' -У ТЛ/2'Ы -У ТЛ/2'Ы -У ТД/2Д'г__—(Рш — рш \

см,;/ — ^0,г+1,7 "Т" ^0,г-1,7 "Г ^о,^-! \г0,г,7+1 М),г,7-17

р0

Г2М

р0

к

с»г I 2Н I +

Р

о,г,7

ир <

(46)

(47)

(48)

где п - шаг по времени при расчете искусственной сжимаемости.

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

• для 111, то есть при к =1 и I = 1 имеем Ж1,'11 = 2(Жо2^1|11,7 - Ж2'11) + Жо1,'^171+1 и

Ж2,11 = Ж2,11

у о,г,7 = У о,г,7+1

для 122 имеем Жо

Р11 к ро

ол7 = + уо;»+1,7-- уо2й2+1 и =

Р22 к ро

12 1,12 1,12 2,12 2,12 2,12 1,12 12 для 112 имеем Уо/г,7 = Уо;г,7+1 + Уо/г+1,7 - и Уо/г,7 = Уо/г,7+1 - к(Ро,2,7 - 1) ;

для 121 имеем = + <¿+1,7 - и = Ж^ц - к(Р02,1,7 - 1) .

1,21

2,21

2,21

2,21

1,21

Этап 2: Численно находим решение системы уравнений (17). Для этого решаем систему уравнений в жидкой части У/

роАсШи _ к1 у.^./

2 дЬ

(49)

и систему в упругом скелете У

ты,

2 дЬ

(50)

где на общей границе «жидкость - упругий скелет» выполняется равенство нормальных напряжений

( у, 1 - Р*1

■ п = (Ло»(у, ^м1) - Рм11) ■ п •

(51)

Пусть V = д^/д£, V = (VI, У2). Тогда, учитывая равенство нормальных напряжений и перемещений, условие на общей границе будет иметь вид

V-

1,к1

»,3

лЛМ , Т./2М , т,/2М _ ^^г,/ 4.7+1 + 4+1,3 + Уг,3 ..

Ро

V-

2,к1 г,3

, -2.1,1

Уг,з+1 ..

Ро

(52)

где

Щ

1,к1

г,3

р к1 г,

г+1,3

Л,

о

Щ2,к1 = Щ 1,.к,' Щ 1,.кг + Щ2;кг ■

(53)

Этап 3: Находим тензора (36), предварительно рассчитав А0 и А1. Для более удобного преобразования тензора в разностном виде проведем процедуру тензорного произведения. Например тензор А0 примет вид

где

ад

«11 «12 «13 «14

«21 «22 «23 «24

«31 «32 «33 «34

\«41 «42 «43 «44

(54)

«11

5ж1

«12 = т:

дЩ12

0,1

+

дЩ21

0,1

дх1 дх1

«13=2

дХ2

+

дЩ11

,2

дх1

«14 = «23 = «32 = «41

+ + + + Ш^Х дх2 дх1 дх2 дх1 дх1 У

\( дЩ1,! дЩ?!Л

«21 = - ( —--V

2 1

«31=2

дх

2

дх1 У '

дх

2

дх

1

«22

«33

<9Щ22

дхт

дх2

1 / дЩ22 дЩ^

«24 = "Г I —--Ь

2 V дх2 1

дх1 У '

2 V дх2

1 /<9Щ22 <9Щ22\

"42= 2 (^Г+ )

«43

| <9Ж212\ дх2 дх2 у

«44

дх2

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

1

Литература

1. Muskat M. Two-fluid system in porous media. The encroachment of water into an oil sand // Physics. - 1934. - 5. - P.250-264.

2. Fahuai Yu. Global classical solution of Muskat free boundary problem //J. Math. Anal. Appl. - 2003. - 288. - P.442-461.

3. Radkevich E. On the spectrum of the pencil in the Verigin-Muskat problem // Sbornik: Mathematics. - 1995. - 80;1. - P.33-74.

4. Siegel M., Caflish R.E., Howison S. Global existence, singular solutions, and ill-posedness for the Muskat problem // Common Pure and Appl. Math. - 2004. - LVII. - P.1-38.

5. Meirmanov A. The Muskat problem for a viscoelastic filtration // Interfaces and Free Boundaries. - 2011. - 13. - P.463-484.

6. Гальцев О. В., Мейрманов А. М. Численное усреднение в задаче Рэлея-Тейлора при фильтрации двух несмешивающихся несжимаемых жидкостей // Математическое моделирование. - 2011. - 23. - P.33 — 43.

7. Nguetseng G. A general convergence result for a functional related to the theory of homoge-nization // SIAM J. Math. Anal. - 1989. - 20. - P.608-623.

8. Nguetseng G. Asymptotic analysis for a stiff variational problem arising in mechanics // SIAM J. Math. Anal. - 1990. - 21. - P.1394-1414.

9. Владимирова Н.Н., Кузнецова Б.Г., Яненко Н.Н. Численные рассчеты симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости // Некоторые вопросы вычислительной и прикладной математики. Новосибирск. - 1996. - P.186-192.

NUMERICAL SOLUTION OF AVERAGED MODEL FOR THE JOIN MOTION OF LIQUIDS WITH DIFFERENT DENSITIES IN POROUS MEDIA

O.V. Galtsev

Belgord State University, Studencheskaya St., 14, Belgorod, 308015, Russia, email: oleggaltsev@yandex.ru

Abstract. The free boundary problem describing the joint filtration of two immiscible incompressible liquids on macroscopic level is studied. Basic equations of mathematical model at microscopic level represent stationary system of Stokes' equations for inhomogeneous incompressible viscous liquid in pores and stationary Lame's system of equations for incompressible elastic solid skeleton. The are supplemented by boundary conditions of continuity of displacements and normal stresses on the common border «solid skeleton-liquid in pores» and transport equations for unknown liquid density. It is built the macroscopic mathematical model of joint motion of fluid and elastic solid skeleton by the method of two-scale expansion. Numerical approximation of the averaged model for the special structure of the pore space is discussed.

Key words: Musket problem, the Rayleigh-Taylor instability, Stokes and Lame's equations.