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

2015 Математика и механика № 3(35)

УДК 532.546: 519.63 БОТ 10.17223/19988621/35/7

Х.М. Гамзаев

О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ДВИЖЕНИЯ ЖИДКОСТИ В ДВУХПЛАСТОВОЙ ВОДОНОСНОЙ СИСТЕМЕ

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

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

Исследование движения жидкостей в водоносных пластах имеет большое значение во многих областях: в орошении и осушении земель, накоплении подземных вод в природных емкостях, искусственном вытеснении ресурсов подземных вод, фильтрационных потерях из водохранилищ, сбросе и захоронении сточных вод, внедрении морской воды в пресноводные горизонты и т.д. [1-3]. Обычно под водоносным пластом понимается проницаемый или слабопроницаемый пласт горных пород, способный в естественных условиях пропускать значительное количество воды. Известно, что при моделировании движения жидкостей в безнапорных водоносных пластах применяются методы гидравлической теории [1]. Согласно этой теории гидравлический напор И в любом сечении пласта принимается за уровень жидкости в этом сечении, а горизонтальные составляющие скорости и и V выражаются через гидравлический напор по закону Дарси

к дИ к дИ

у дх ' у ду '

где х, у - горизонтальные оси; к - коэффициент фильтрации; у = pg - удельный вес жидкости, р - плотность жидкости, g - ускорение свободного падения.

Подставляя эти соотношения в уравнение неразрывности фильтрационного потока

дтрИ + дирИ + дvрИ = о д/ дх ду

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

дтрИ д дИ д дИ

—;— = — (кРИ —) + — (кРИ — ), (1

д/ дх дх ду ду

где / - время; т - коэффициент пористости.

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

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

Предположим, что рассматривается движение несжимаемой жидкости в неде-формируемом, горизонтально расположенном замкнутом цилиндрическом водоносном пласте протяженностью Я . В центре пласта расположена гидродинамически совершенная скважина радиусом гу/. Пренебрегая вертикальными составляющими скорости фильтрации, поток в пласте можно считать плоскорадиальным и осесимметричным. Тогда уравнение движения (1) для данного фильтрационного процесса в отсутствие инфильтрации на свободной поверхности примет вид

т — =1 — (гкк —), (г, Г) еП = [гК < г < Я, 0 < Г < Т} . (2)

дt г дг дг

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

w = -Х(к -Н), ё

где х - коэффициент фильтрации подстилающего прослоя, ё - его толщина.

Учитывая перетоки из нижнего пласта, уравнение (2) можно представить в виде

дк 1 д , дк х „ т— =--(гкк—)- — (к - Н),

д г дг дг ё

(г, 0 <еО. = [гч, < г < Я, 0 < t <Т} . (3)

Предположим, что в начальный момент времени t = 0 распределение гидравлического напора в пласте известно, т.е. для уравнения (3) имеем следующее начальное условие:

к(г,0) = ф(г). (4)

Предполагая, что изменение во времени расхода жидкости на скважине описывается функцией q(t), граничное условие для уравнения (3) при г = ^ можно

записать в виде

дк

2пгкк — = q(t), г = ^ , 0 < t < Т . (5)

дг

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

дк

— _ 0, г _ Я, 0 <Г < Т . (6)

дг

Необходимо отметить, что гидравлический напор в нижнем пласте Н (/) недоступен для непосредственного измерения. Следовательно, функция Н(/) неизвестна и также подлежит определению одновременно с решением к(г, ?) задачи (3)-(6). Очевидно, что для корректной постановки задачи необходимо задавать дополнительное условие. Предположим, что дополнительное условие для уравнения (3) задается на скважине

к(г, 0 = ДО, г = г, , 0 < / < Т . (7)

Таким образом, задача заключается в определении функций к(г, ?), Н(?), удовлетворяющих уравнению (3) и условиям (4)-(7). Задача (3)-(7) относится к классу обратных задач, связанных с восстановлением зависимости правых частей параболических уравнений от времени [4].

2. Метод решения

Уравнение (3) и условие (4) - (7) представим в безразмерной форме. Следует отметить, что представление уравнение в безразмерной форме позволяет выбрать диапазон изменения безразмерных переменных таким образом, чтобы улучшить обусловленность задачи. Введем следующие безразмерные переменные:

:_г, 1 й= к Н т= к -= 1

Я г к к к 1

t * t h — h h , h — H h*' k k — k*,

— — — h*' r w ii R|r T — T * t w=

f — J — — ^ — — w Т — ' — — R Х

h' h' R t* h*k*d

* * * mR2 * * *2 где h , k , t — * *, q — 2nk h - характерные размерные величины. k h

Опуская черточки над безразмерными переменными, задачи (3) - (7) запишем в виде

dh 1 д dh

— —--(rkh—) - wh + wH , (r, t) eO. — [rw < r <1,0 < t <T} ; (8)

dt r dr dr

h(r ,0) — —(r); (9)

dh

rkh— — q(t), r — rw , 0 < t < T ; (10)

dr

dh

— — 0, r — 1, 0 < t < T ; (11)

dr

h(r, t) — f (t), r — rw , 0 < t < T . (12)

Для численного решения задачи (8) - (12) используем подход, предложенный в [5, 6].

С этой целью введем равномерную разностную сетку

Jh%

= {(r,tj) :rt = Го + ;Дг, t}. = JAt, r0 = rw,i = 0,1,2,...N, j = 0,1,2,..M}

с шагами Дг = 1/ Ж по переменной г и Дt = Т /М по времени t. Пользуясь ин-тегро-интерполяционным методом, разностный аналог задачи (8)-(12) на сетке юкт представим в виде

hj+1 - h;j Дt

1

r Дг

. h+l - hj+1 r i k i h i--r i k i h i

; +— ;+— ;+— Дг ; — ;— ;— . 2 2 2 2 2 2

hj+i -

Дг

- whj+i +wHj+i,

; = i, 2, 3, ....N -1, j = 0,i, 2, 3,....,M -1

h =ф;, ; = 0, n.

- = q

hj+1 - hj+1 r k hj hi h0 = qj+1

Дг

j+1 - uj+1

h] 1 - h tl-Kj tl

■ = 0,

Дг

hj+1 = fj+1,

где

H

Ф =ф(г ), q1+1 = q(tj+i), f1+1 = f (tj+i), hj - h(r, tj ),

j+1 - H(tj+i), r± 1 = r ±ДТ, * 1 = «S^, hj 1 = ^ .

2

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

r; Дг

2

- ch+1 + Ъ] = -(-1-h + wr Дг2Hj+1), ; = 1, N -1, j = 0,M -1; (13)

Дt

h» =ф, ; = 0, N ; hj+1 = hj+1 + v ;

hj+1 = hj+1 •

"n "n-1 '

h0j+1 = fj+1, j = 0,1,2,..M-1,

r Дг 2 2

где a = r ik ih i; b; = r ik ihj i; c; = a + Ъг +J— + wrДг , v = -

(14)

(15)

(16) (17)

,j+i

; — ; — ;— 2 2 2

; +— ; +— ;+-

2 2 2

r0 k0 h0

Представим решение задачи (13)-(17) в виде

к/+1 =амк£ +Рг+1, I = 0,1,2,..., N -1, (18)

где аг+1, Рг-+1 - неизвестные пока коэффициенты. Запишем аналогичное выражение для :

кД1 =а к/+1 .

Подставляя выражения для к/+1, к/-! в уравнение (13), получим следующие нелинейные уравнения для определения коэффициентов аг+1, Рг+1:

г Аг 2

, аР + --к + Аг 2 Н}+1

п 1>1 . I I

а __ь_ в __А1_

"-¿+1 ' Рг+1 '

Сг - аг аг Сг - аг аг

I _ 1, 2, ..., N - 1.

Начальные значения коэффициентов а,, рг находим из требования эквивалентности соотношения (18) при г _ 0 , т.е. к^1 _ аД-^1 +р1, условию (15):

а1 _ 1, Р1 _ V.

Нелинейное уравнение для Рг+1 преобразуем к виду

вг+1 --а^- р, ы +^^ н]+1,

С - аг аг А (сг - аг аг ) Сг - аг аг

или р,+1 _ s1 рг + у + 2Н}+\ (19)

г Аг2 , ■ wг, Аг2

где - _-1-, Уг -Г2г _

Сг - аг аг А(сг - аг аг )

Введем новые переменные рг, 2г, г _ 1,2,..., N, удовлетворяющие уравнениям

Р г +1 _ ^ г + Уг , Ъ+1 _ ^г + Ъ ,

Р1 _в1, Ъ _ 0.

С учетом вновь введенных переменных, уравнение (19) можно представить в виде рекуррентного соотношения

Рг+1 _Рг+1 + ^+1Н+1, г _ 1,2,..., N -1. (20)

Теперь найдем зависимость между к+1 и к^1 в явном виде. Для этого соотношение (18) запишем при г _ 0 :

■ _а1к/+1 +Р1.

Подставив сюда выражение для к(+х, т.е. к(+1 _ а2к^+1 +Р2, будем иметь к0+1 _а1(а2 к{+1 +Р2) + Р1.

Далее, подставляя в последнее уравнение выражения для к^+1,к^+1,...,к^, получим формулу, в которой к10+1 выражается через к1^1:

N N г -1

к0+1 _ К+1 Паг +!Рг Паг +Р1. (21)

г _1 г_2 1_1

Теперь исключив из следующей системы уравнений:

кЖ-1 = аЖкЖ + Рж , и+1 = 1

"Ж "ж-1 '

получим соотношение, связывающее к^ и :

=. (22)

1 -аж

Подставляя соотношение (22) в уравнение (21), будем иметь

в N N I-1

к0+1 = ^ Паг + Паг +Р1.

1 аЖ г=1 I=2 1=1

Последнее уравнение с учетом рекуррентного соотношения (20) запишется в виде

в N = NN г-1 N I-1

ко;+1 = Паг + Н^1---^ +!вг Паг + Н+1 ^ Па,

1 аЖ 1=1 1 аЖ 1=1 1=2 I =1 1=2 1=1

Из полученного уравнения можно найти приближенное значение искомой функции Н ^) при t = tj+1 :

N N 1-1

Ч+1 -^ПаI -ЕРIПа,-в1

Н1 +1 =_1 аЖ I =1 I=2 1=1

Па +Е ^ Па,

1 -аЖ I=1 г=2 , =1

Определив Н1+1, по формуле (22) можно найти к^1 и далее по рекуррентной

формуле (18) последовательно определить к^,к^^,...,к11+1. При переходе на следующий временной слой описанная процедура вычислений снова повторяется.

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

3. Результаты численных расчетов

Для выяснения эффективности практического применения предложенного вычислительного алгоритма были проведены численные эксперименты для модельных задач. Схема численного эксперимента заключалась в следующем. Для заданных функций Нq(t),ф(г) решалась прямая задача (8)-(П). Найденная зависимость к(гМ1, t) = / ^) считалась точной и использовалась для численного решения обратной задачи по восстановлению Н ^).

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

/ ^) = / ^) + 5ст^),

где ст(/) - случайный процесс, моделируемый с помощью датчика случайных чисел, 5 - уровень погрешности.

Расчеты выполнялись на пространственно-временной разностной сетке с шагами к = 0.0499, т = 0.04. Результаты численного эксперимента, проведенного для

п

случая к = 2, ^ = 40, q = 0.02, Н(?) = 20 - ^ш—?, ф(г) = 15, с использованием

невозмущенных и возмущенных входных данных представлены в таблице; в ней ? - время, Н - точные значения функции Н ), Н и Н - вычисленные значения функции Н ) при невозмущенных и возмущенных данных. Для возмущения входных данных использовались погрешности уровня 5 = 0.5 и 5 = 1.2.

Результаты численного эксперимента

т Н Н Н

5 = 0.5 5 = 1.2

0.4 17.966 17.978 18.207 18.529

0.8 16.284 16.290 16.314 16.348

1.2 15.245 15.246 15.214 15.170

1.6 15.027 15.027 14.785 14.446

2.0 15.670 15.672 16.006 16.475

2.4 17.061 17.068 17.097 17.137

2.8 18.960 18.974 18.874 18.735

3.2 21.040 21.059 21.577 22.303

3.6 22.939 22.963 23.574 24.429

4.0 24.330 24.358 24.768 25.343

4.4 24.973 25.002 25.165 25.395

4.8 24.755 24784 24.923 25.117

5.2 23.716 23.742 24.100 24.600

5.6 22.034 22.057 22.369 22.805

6.0 20.00 20.018 20.393 20.919

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

Заключение

При исследовании процессов в многопластовых водоносных системах необходимо учитывать обменные процессы между пластами.

ЛИТЕРАТУРА

1. Полубаринова-Кочина П.Я. Теория движения грунтовых вод. М.: Наука, 1977.

2. Веригин Н.Н. и др. Гидродинамические и физико-химические свойства горных пород. М.: Недра, 1977.

3. Бондаренко Н.Ф. Физика движения подземных вод. М.: Недра, 1973.

4. Самарский А.А., Вабищевич П.Н. Численные методы решения обратных задач математической физики. М.: Издательство ЛКИ, 2009.

5. Гамзаев Х.М. Численный метод решения обратной задачи поршневого вытеснения нефти из пласта водой // Инженерно-физический журнал. 2012. Т. 85. № 5. С. 925-930.

6. Гамзаев Х.М. Численный метод решения обратной задачи упруговодонапорного режима разработки нефтяного пласта // Вычислительная механика сплошных сред. 2012. Т. 5. № 4. С. 392-396.

Статья поступила 16.12.2014 г.

Gamzaev Kh.M. ON NUMERICAL SIMULATION OF THE FLUID FLOW IN A DUAL-COMPLETION WATER-BEARING SYSTEM

DOI 10.17223/19988621/35/7

The motion of an incompressible fluid in a one-dimensional water-bearing stratum is considered in the presence of a pressure communication with the underlying stratum. The hydraulic pressure head in the lower stratum is considered to be unknown. The question of modeling the process is reduced to solving the inverse problem of reconstructing the right-hand side of differential equations of the fluid motion. A finite-difference analog of the inverse problem is constructed and a computational algorithm for solving the resulting system of finite-difference equations is proposed.

Keywords: water-bearing stratum, pressure communication, inverse problem of reconstructing the right-hand side, finite-difference method.

GAMZAEV KhanlarMekhvaly ogly (Doctor of Technical Sciences, Azerbaijan State Oil Academy, Baku, Azerbaijan) E-mail: xan.h@rambler.ru

REFERENCES

1. Polubarinova-Kochina P.Ya. Teoriya dvizheniya gruntovykh vod. Moskow, Nauka Publ., 1977. (in Russian)

2. Verigin N.N. et al. Gidrodinamicheskie i fiziko-khimicheskie svoystva gornykh porod. Moskow, Nedra Publ., 1977. (in Russian)

3. Bondarenko N.F. Fizika dvizheniya podzemnykh vod. Moskow, Nedra Publ., 1973. (in Russian)

4. Samarskiy A.A., Vabishchevich P.N. Chislennye metody resheniya obratnykh zadach matematicheskoy fiziki. Moskow, Izdatel'stvo LKI, 2009. (in Russian)

5. Gamzaev Kh.M. Chislennyy metod resheniya obratnoy zadachi porshnevogo vytesneniya nefti iz plasta vodoy. Inzhenerno-fizicheskiy zhurnal, 2012, vol. 85, no. 5, pp. 925-930. (in Russian)

6. Gamzaev Kh.M. Chislennyy metod resheniya obratnoy zadachi uprugovodonapornogo rez-hima razrabotki neftyanogo plasta. Vychislitel'naya mekhanika sploshnykh sred, 2012, vol. 5, no. 4, pp. 392-396. (in Russian)