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

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Том 148, кн. 2

Физико-математические пауки

2006

УДК 519.63^532.5

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ФИЛЬТРАЦИИ ДВУХ НЕСМЕШИВАЮЩИХСЯ ЖИДКОСТЕЙ

Р.Ф. Кадыров, A.B. Лапин

Аннотация

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

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

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

Введение

Исследование процесса проникновения органического, не смешивающегося с водой, загрязнителя во влажную почву представляет интерес для физики почв и экологии [1]. Классические модели фильтрации [2] недостаточны для адекватного описания происходящих процессов. В частности, они не описывают явление остановки фронта проникновения загрязнителя в почву при непрерывно продолжающемся контакте поверхности почвы с загрязнителем [3]. Достаточно общие трехмерные математические модели процессов фильтрации в почвах и почвогрун-тах не смешивающихся с водой жидкостей приведены в [4. 5]. Применимость этих моделей к исследованию реальных процессов во многих случаях ограничена, поскольку широкий охват большого количества сопутствующих процессу механизмов влечет за собой необходимость определения большого количества параметров модели. Численное решение этих задач также представляет большую трудность.

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

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

Статья организована следующим образом.

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

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

Неявная сеточная аппроксимация и итерационный алгоритм решения предложены в п. 3.

Результаты численных экспериментов включены в последний пункт статьи.

Будем рассматривать одномерную задачу фильтрации на временном промежутке (0, Т). Пусть область фильтрации Q = (0, Ь) х (0, Т) разбита гладкой кривой Г = {(г,г) € Q : г = £(г)} та да подобласти: Q- = {(г, г) € Q : 0 < г < £(г)} и Q+ = {(г, г) € Q : £(£) < г < Ь}. В физической модели Q+ соответствует не загрязненной подобласти, насыщенной водой и воздухом, в то время как во вторую подобласть Q- проник загрязнитель, и она находится в условиях трехфазного насыщения - водой, воздухом и углеводородом. Кривая г = £(£) является неизвестной подвижной границей фронтом продвижения загрязнителя. Таким образом, если Б Б л„ Б - ияритпоннартп ипттт.т газа и углеводорода соответственно, то

где , Я0 - скорости фильтрации, а Н,, Н0 - гидродинамические напоры для воды и углеводорода соответственно.

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

1. Постановка задачи

Бш + Бд = 1, Бо = 0 при (г,г) € Q+,

Б/ш + Бд + Бо = 1, Бо > 0 при (г, г) € Q . Уравнения баланса масс воды и углеводорода имеют вид

да а* 1Й = ~

т

а

при (г, г) € Q

(3)

а

Все коэффициенты в (1) (3) являются функциями насыщонностой (эти зависимости приведены ниже) и имеют следующий физический смысл: m,k пористость и проницаемость среды,

fi (индексы принимают значения w, o) — вязкость и относительная фазовая проницаемость г-й фазы,

а* - максимально возможное количество впитываемого в почву загрязнителя («емкость»), т = const > 0.

Зависимости гидронапоров от насыщонностой даются уравнениями

Hw = H+ (z) = -Pc (0W) - pw gz ири (z,t) G Q+, (4)

Hw = H- (z) = ~^PC (1 - вд) - —Pc (вш) - Pwgz при (z,t) G Q-, (5)

®wg ®wg

H0(z) = -^Pc(l-eg)-p0gz при (z,t)eQ~, (6)

®wg

где (Jij > 0 - поверхностное натяжение на поверхности раздела фаз inj (индексы принимают значения w, о, д), ei = (Si — Sir)/(1 — Sir) - эффективная насыщенность пористой среды i-й фазой, Sir - величина остаточной насыщенности данной фазой, pi - плотноеть i-й фазы, д - величина гравитационного ускорения, капиллярное давление Pc определено равенством (уравнение Ван-Генухтена ( [7])):

= 4 («"'<'-"»-l)"",

где натуральное число n G N - константа, характеризующая геометрию поровой структуры.

Далее под S—, S— понимаем следы фупкций So, Sw на Г со стороны области Q- , под S+, S+ - следы функций So, Sw на Г со стороны области Q+. Аналогичный смысл имеют функции q—, q— и q+, q— . Из балансовых соотношений наряду с уравнениями в подобластях Q- и Q+ можно получить следующие уравнения на фронте:

т§ №(*)) "ЫШ = q+W)) ~q„W)), (')

™§S-m) = q^(№)- (8)

Sw So

на фронте z = £(t), система дифференциальных уравнений дополняется двумя условиями. Первое это условие непрерывности напора воды:

н— ш) = h+ т). (9)

Второе условие связывает насыщенности воды и углеводорода на фронте:

S— = ф^—), (10)

где ф : [0,1] ^ [0,1] - заданная гладкая функция.

Зависимость относительных фазовых проницаомостой от насыщонностой принимается в виде

( f // -,а("-1)/"\ 2 i it -n\2( n-1)/n

fw = ei (^1 — (i — ew/(n-1)) J , fo = ey (i — (1 — eo)n/(n-1))

с положительными постоянными £ и 7. Кроме того, m = m0 — a, где m0 - начальная пористость среды (m (z, 0) = m0), k (m) = Km3, K = const > 0.

Уравнения (1) (3) дополняются начальными и граничными условиями:

So (z, 0) = 0, Sw (z, 0) = SW (z);

(H)

Sw (0, t) + So (0, t) = 1, Sw (L, t) = 1;

a (z, 0) = 0,

(12) (13)

2. Математическая модель

Основная задача данного раздела это объединить дифференциальные уравнения (1). (2) и условия на фронте (7). (8) в так называемую сквозную дифференциальную постановку рассматриваемой задачи.

Пусть функция Ьа = Ь0(£) > 0 при каждом £ € [0, Т] равна следу (£) на Г. Очевидно, что Ьа - это скачок насыщенности Б0 на фронте Г. Подобным образом определим функцию = = Б+ — Б—, равную щи каждом £ € [0, Т] скачку-

насыщенности та фрон те Г.

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

k

Оо = О,о оj — f0(60) —

Mo ^wg(1 — Sgr ) a

1 p' ( 50 + Su

1 — Sg

К EE b0(S0) = —f0{S0) [Род-^Щрс

2 1 — Sg

bw =

к

Mo \ Wwg

= b° + -

^w g ^w g

(Sw) = bw + bw HPH(z,t) e Q+

gr

'b- (SW,S0) = b°w + ^bi + ^bl при (z, t) G g

— / c< с \ 1 °wo 2 / ГЛ —

aw (bw, b0) = -aw H--aw при (z,t) G Q ,

^wg ^wg

a

+ (Sw,So) = aw npn(z,t) e Q+.

Здесь

k

Pc(0) = aPc(0), b°w = —fw(Sw)Pwg

b<W fw(Sw ) n Pc

Mw

So + Sw

1 — Sgr

fw (Su

fw (Su

Mw a

1 p' ( SQ Sw

Sw Swr

1 - Sw

1 — Sgr

rPc

1 — Sgr 'I Sw Swr

(1 — Swr) c V 1 — S

Из (2) и (8) с учетом введенных обозначений следует:

d(mSo) dq

dt

+

dz

da d{S0 + Sw) r)_ -—, q0 = -a0-—--h b0{b0) в g

dt

dz

mL0(t)— = q0 на Г,

a

w

i

a

w

w

2

a

w

где Яо = -а0 {Ба , Бш) д(5о + Ъ0(Б0 ).

Аналогично, из (1) и (7) получаем

¿Ктб'ш) д(Ба + Бт) +

-^--1—^— = 0, ди1 = -а,и1-т--\-Ьи, в С} и ,

дг дг дг

(15)

где

Чт ~ ~аи1 ^ш) ° ~ + Ьш(5+),

Чш — ~аш о ~

Определим функции

(Ба(г, г) - Ьа(г) при (г, г) € д-, и(г,г) = <

[о при (г, г) € д+ и Г,

{Бы(г,г) + (г) при (г, г) € д-,

Бш при (г, г) € д+ и Г.

Функции и и ад те терпят разрыва при переходе через Г, и для них, в предположении, что известны Ьа(€) и Ьт(г), будет приведена сквозная постановка задачи (14), (15).

Пусть п(г, г) = (г, г) € д ; ^и (г, г) € д+} - характеристическая

функция области д- и пусть уи(г, г) - произвольная, достаточно гладкая в д

функция, обращающаяся в нуль на границе дд. Обозначим через п = (п , п4)

д- Г

= — п4/п2, традиционным методом доказываем эквивалентность системы (14) интегральному тождеству

J + Ь0 ^ + Шо-тр^ ЛгсИ = - J ,-.,,¡-.,¡1 \/ии.

Я Я

Поскольку в подобласти д+ справедливо равепство и = 0 , то слагаемое ди д«и д-уи

ца0——-— . входящее в выражение — в этом интегральном тождестве, можно дг дг дг

ди д«и

заменить на а0——-—. Таким образом, задача (14) при известных го, Ьш, а и в дг дг

предположении требуемой гладкости всех функций эквивалентна интегральному тождеству

д(ад — п)

дг дг

— цЪо^г-) = — [ г1—-1>исЬси Уг>и. (16) дг ) у от

Аналогично, система (15) эквивалентна интегральному тождеству

J ^т(-ги + Ьш + Я

( _ д-ш + дгю\ дУы , _ + д'Оы

+ ( »7<»ш-- + [Х-П1)а1— \-qZ- (ЯК + (! -тЮ — +

_ д(и + Ьр ??) + д(и + Ьр ??) \ душ \

—~д1— + —о;—(Ь)

где - произвольная, достаточно гладкая функция, обращающаяся в нуль

на границе д( области Q.

Из (16), (17) следует, что в области ( справедливы следующие уравнения, которые понимаются в смысле обобщенных функций:

д(т(и + Ьа п)) д ди д ( д(т - Ьшп)\ , д да

Що-т-- + т~(Фо) = (18)

dt dz ° dz dz \ 1 ° dz 1 ^ dzy 1 °J iï dt-'

d[m{W-t LW 11)] ~ l (Varw + (1 - „) + | blb- + (1 - „) Ъ+) -

" -qZ (')«» + i1 " V)K)-Q;-=°- (19)

Здесь аю = ao (и + Lo n,w - Lw n), 60 = 60 (u + Lo n), aw = aw (w - Lw ц,и + Lo n) и bw = bw (w - Lw n) ■

В дальнейшем считаем, что на границе z = 0 насыщенность So равна максимальной So = 1 Swr , поэтому справедлива следующая физически мотивированная гипотеза:

S0(z,t) > Lo(t) при всех (z,t) G Q-.

Поскольку согласно этой гипотезе функция и положительна в Q-, то n = H (и), где H (и) = и < 0; 1 при и > 0} - функция Хевисайда. Далее под

H (и) понимаем многозначный график функции Хевисайда: H (и) = {0 пр и и < < 0; [0,1] при и = 0; 1 при и > 0}. Заменим характеристическую функцию ц области Q- в уравнениях (18), (19) функцией, удовлетворяющей включению

n(z,t) G H(u(z,t)) при (z,t) G Q. (20)

a

д a a* - a ^

-7Г7 = V- при (z,t) G Q (21)

dt t

и начальные н краевые условия

Uo (z, 0)= Lo (z, 0)=0, w (z, 0) = Sw (z, 0)= Sw (z), (22)

w (0,t) = Sw (0), и (0,t) = 1 - Sw (0) - Lo(t), w (L,t) = 1; и (L, t) = 0, (23)

a (z, 0)=0. (24)

Назовем обобщенным решением задачи (14), (15) при известных Lw(t) и Lo(t) тройку функций (u,w,ri), удовлетворяющую системе соотношений (18)—(24). Мы

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

Для определения функций Ьт(£) и Ьа(£) (то есть скачков насыщенностей на фронте) используем условия (9) и (10):

—(тг^ + —-Рс + ^ - = рс ' (25)

\ 1 5 дт / ^тд \ 1 / \ 1 /

57 = (26)

В этих уравнениях под 57, понимаются значения соответствующих функций в точке 5 где

£"(£) = а^ир{г : п(г,£) = 1}, а под 5+ - значение в точке £+(£), где

(¿) = а^т1 {г : п(г,£)=0}.

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

57 = «(Г(^) + Ь0(«), = ™(Г (*),*)+ Ьш(*), = (£),£).

3. Конечномерная аппроксимация и алгоритм решения

Уравнения (18) (21) аппроксимируем конечно-разностной схемой на равномерной сетке. Пусть Н > 0 - шаг сетки по пространственной переменной г, Д£ > 0 — шаг по времени. Далее для сеточных функций, аппроксимирующих функции от непрерывных аргументов г и будем использовать те же обозначения. Пусть разностные производные по г и по £ определены равенствами

_ - - дг) _ ф + /?) - ф) _ ф) - ф - /?.)

31 = Ш ' = Тг ' 85 = Тг '

Тогда неявная сеточная аппроксимация уравнений (18) (20) имеет вид

(т(и + Ьа п))( - («аИ*)г - (п«а(^ - Ь^п)*)г + (пЬа)г = -п«ь (27)

(т(ад - Ьт - ((п«- + (1 - П) «+) )г + (п Ь- + (1 - П) Ь+)г -

- ((п «7 + (1 - п) ) (и + Ьа п)*)г = 0, (28)

п е Н(и). (29)

Коэффициенты в сеточных уравнениях определены в «полуцелых точках». Например, в выражении (ааи2)г под аа понимается

а0 = а0 ^(г + (г + /г)), ^(ф) + гф + Н)), ^(гф) + гф + /г))^ .

В коэффициентах Ь^ и Ь^ производная а* аппроксимируется разностным отношением аг.

Сеточная аппроксимация (21) имеет вид

*

« - «

0( = ?7-. (30)

Скалярные уравнения (25) и (26) используются для определения Ьо и Ь» на расчетном слое. Алгоритм решения задачи описан ниже.

Итерационный алгоритм решения системы сеточных уравнений на слое £ состоит в следующем.

Задаем начальные приближения для Ьо(£) и Ь»(£), например, Ь0(£) = Ьо(£—Д£) и Ь»(£) = Ь»(£ — Д£). Для к > 1 находим пару сеточных функций (ик,пк) из уравнений

п* — пк-1 1

+ У^Ь + »7* ' г' = + Ь0 77))|(4_до + («ГЧ^Ь, (31)

пк е н(ик), (32)

затем функции и пк - из уравнений

¿(т*-1^* " ^ V")) = ¿М™ - »?))|(*-до +

+ (пк п-'к-Чк-1 )- — (пк п-'к-1 (ик + ¿£-1 пк )г )- +

+ ((1 — пк)п+' *-Ч*-1 )- — ((1 — пк)п+' к-V + Ьк-1 ^)*)- +

+ (пк&-'к-1 )- + ((1 — пк№к-1 )-, (33) 1 1 п* —

¿»'=¿»(,-^+,'4^ (34)

Отметим, что и пк находятся из уравнений (33), (34) по явным формулам. Для определения (ик,пк) из уравнений (31), (32) требуется для каждой точки сетки (последовательно от г = К до г = Ь — К) решить систему из скалярного уравнения и скалярного включения

Аик + Впк = Р пк е н(ик)

с известными А, В и Р, зависящими от предыдущих итераций. Решение этой системы находится в явном виде.

По найденной сеточной функции пк определяем точки сетки £-'к (£) = = а^тах{г : пк(г, £) = 1} и £+'к(£) = а^тт{г : пк(г,£) = 0} и находим новое приближение

5-,к = ™к (г,к (*),*) + ь»-1.

Далее находим новые приближения 51-'к и к из системы уравнений (25) и (26). Для их вычисления используем явные формулы

5-'к = ),

_ 5-, к + 5-, к

1 —

Л 5-, к — 5

"»о 0 ,\ , о—.к л+,к\ \ , е

Н--1 _ у-) + (4 "С '

г = \ksKt)

Рис. 1. Положение свободной границы х = £(4)

Теперь равенства Ька = , Ь^ = Б:+'к — дают новые приближения для

скачков Ьа(€) и (Ь) на расчетном слое. Критерием остановки итерационного метода выбрана малость разности текущих приближений для скачков |Ь^ — Ь^—11

и ь — ьк—1|.

4. Вычислительные эксперименты

Были проведены вычислительные эксперименты по описанной схеме для трех разных наборов исходных данных, различающихся между собой видом функции ф(Б—) и максимальным значением набухания а* . Значения остальных параметров были следующими:

Параметры в уравнениях Вап-Гепухтепа и в формуле для капиллярного давления п = 3 б= 0.5 7 = 0.33

Поверхностные патяжепия (Н/кг) (7 од = 0.037 аюд = 0.07 а.юо = 0.022

Остаточные насыщенности Я от =0.1 = 0.1 V = о.1

Вязкость воды и нефти (кг/м3) = 0.001 р,а = 0.002

Глубина колонки (м) Ь = 4

Начальное значение пористости то = 0.5

Плотность воды и нефти (кг/м3) рт = 1000 Ро = 800

Скорость набухания т = 10ь

Константа в формуле для абсолютной проницаемости К = кг11

SJz),

t=0.004, 0.502, 1

5

- - a* = 0.01, A = 0.5 - a* = 0.03, A = 0.5

— a* = 0.01, A = 0.3

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Рис. 2. Насыщенность углеводорода

Рис. 3. Насыщенность воды

Z

Sw(z), t=0.004, 0.502, 1

Зависимость ф между следами иасыщеииостей на фронте во всех расчетах задавалась формулой S— = AS— (1 — S-), где A = const - некоторый параметр.

Вычислительные эксперименты были проведены для различных значений макси-а* A

A. а* = 0.01, A = 0.5;

B. а* = 0.03, A = 0.5;

C. а* = 0.01, A = 0.3.

На рис. 1 3 приведены графики численных решений: пунктирная линия соответствует первому эксперименту, сплошная утолщенная линия второму эксперименту. а сплошная тонкая линия третьему эксперименту.

На рис. 1 представлен график функции z = £(t), которая характеризует положение фронта в зависимости от времени. На рис. 2. 3 показаны значения абсолютных насыщенностей воды и углеводорода в разные моменты времени.

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

Работа выполнена при финансовой поддержке РФФИ (проект Х- 04-01-00484) и МНТЦ (проект Л* 2419).

Summary

R.F. Katlyrov, A.V. Lapin. Mathematical model and numerical solution of filtration problem for two immiscible fluids.

In this paper, a one-dimensional nonlinear time-dependent, problem is studied, the problem being a mathematical model of a flow of two immiscible fluids in porous medium. Fluid flow is described by two nonlinear parabolic degenerated differential equations, which act in subdomains separated by a priori unknown moving boundary. The main feature of this problem is that, both functions and fluxes have jumps through unknown boundary. On this boundary some additional conditions are imposed to connect, jumps of functions with jumps of fluxes.

The studied problem is formulated in a fixed domain. Scalar relations between values of solution functions on a moving boundary give additional conditions for finding jumps of these functions.

The implicit, finite-differences approximation of the problem is constructed and algorithm for its implementation is suggested. The numerical experiments demonstrate the adequat.eness of constructed model to physical process and good efficiency of proposed algorithm.

Литература

1. Современные проблемы гидромеханики и гидрогеологии. Сб. докл. копф. СПб., 2002. 572 с.

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

3. Бреус И.П., Игнатьев Ю.А., Ефстифеева Е.В., Неклюдов С.А., Бреус В.А. Миграция углеводородов в увлажненном до полевой влагоемкости выщелоченном черноземе // Докл. РАСХН. 2002. Л» 4. С. 34 37.

4. Finder G.F., Abrióla L.M. On the simulation of nonaqueous phase organic compounds in the subsurface // Water Res. Res. 1986. V. 22, No 9. P. 109 111.

5. Guarnaccia J., Finder G., Fishman M. NAPL: Simulator Documentation. U.S. Environmental Protection Agency. EPA/600/SR-97/102. Washington, DC, 1997. 9 p.

6. Кадыров Р.Ф., По'тагиев К.А., Якимов Н.Д. Транспорт жидкого углеводородного загрязнителя в увлажненной суглинистой почве: математическая модель и численная реализация решения задачи // Аппот. докл. IX Всеросс. съезда по теорет. и

прикл. механике, Нижний Новгород, 22-28 авг. 2006 г. Нижний Новгород: Изд-во Нижегородск. ун-та, 2006. С. 98.

7. Van Genuchten M.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils // Soil Sci. Soc. Am. J. 1980. V. 44. P. 892 898.

Поступила в редакцию 12.06.06

Лапин Александр Васильевич доктор физико-математических паук, профессор кафедры математической статистики Казанского государственного университета. E-mail: alapinQksu.ru

Кадыров Рафаэль Фаридович научный сотрудник НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета. E-mail: rkadyrovQksu.ru