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

УДК 519.218.82

Численное решение одномерной задачи фильтрации несжимаемой жидкости в вязкой пористой среде*

Р.А. Вирц1, А.А. Папин1, В.А. Вайгант2

1 Алтайский государственный университет (Барнаул, Россия)

2 Боннский университет (Бонн, Германия)

Numerical Solution of the One-Dimensional Problem of Filtration of an Incompressible Fluid in a Viscous Porous Medium

R.A. Virts1, A.A. Papin1, W.A. Weigant2

1 Altai State University (Barnaul, Russia)

2 Bonn University (Bonn, Germany)

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

Ключевые слова: пористость, фильтрация, поро-

упругость, закон Дарси.

БОТ 10.14258/izvasu(2018)4-11

The process of fluid filtration in a deformable porous medium is described by a system of equations consisting of the equations of mass conservation for the liquid and solid phases, the Darcy law, the Maxwell-type rheological relation, and the law of conservation of the balance of forces. It is assumed that the poroelastic medium has predominantly viscous properties and the phase densities are constant. In the case of one spatial variable, the transition to Lagrange variables makes it possible to reduce the initial system of determining equations to one equation for the required porosity. The paper is aimed to study the initial — boundary value problem numerically. In clause 1, the task is set and a brief review of the literature on related work. In paragraph 2, the system of equations is transformed, resulting in a nonlinear third-order equation for porosity. In paragraph 3 we propose an algorithm for the numerical solution of a one-dimensional initial-boundary problem. For numerical realization, we use a homogeneous difference scheme for a second-order equation with variable coefficients and a Runge — Kutta scheme of second-order approximation. The obtained solution satisfies the physical principle of maximum. In paragraph 4, we consider a more general case of reducing the original system to a single equation.

Key words: porosity, filtration, poroelasticity, Darcy law.

1. Постановка задачи. В работе изучается следующая квазилинейная система уравнений составного типа:

* Работа выполнена при финансовой поддержке гранта РФФИ № 16-08-00291.

d^Pf dt

dps(l - Ф) dt

+ div(^vf pf ) = 0,

(1)

+ div((1 - ф)щр^)=0, (2)

- Яз) = -k(ф)(yРf + Pfg),

Ф0(9))3(9,Ь), где 3(9,Ь) = §(9,Ь) (3) рехода. Вместо (1)-(6) имеем

якобиан пе-

divvз = -а\(ф)ре - а2(Ф)(+ • Уре),

Ре = Ргаг - Pf, У Ргаг - Ршд = 0,

(4)

(5)

Ргаг = ФРf + (1-Ф)Рз, Ргаг = ФРf + (1-Ф)рз. (6)

Система (4)-(6) решается в области (х,Ь) е Qт = о х (0,Т), о = (0,1), при краевых и начальных условиях

\х=0,х = 1= vf \х=0,х = 1 = ° ф

г=0=

Ф0 (х).

дх

дТ = 'из(х,Т Т),

\т

г= х.

дФ (1 - Ф ) д

дЬ + 1 - Ф0 дх

(фvf)

(1 - Ф) д

1 - Ф0 дх

Ф,

д(1 - Ф) (1 - Ф)2 дУз

дЬ

+

1 - Ф0 дх

Ф('з - vf) = дЬ - ^9),

(1 - ф) д'з 9 а 9 дре

т-ф0дхз = -а1(ф)ре - а2(Ф)Ж,

(1 - Ф) дргаг

1 - Ф0 дх

= - 9гаг9.

Данная начально-краевая задача описывает одномерное нестационарное движение магмы в горной породе. Для описания процесса используются законы сохранения масс для каждой из фаз, закон Дарси для жидкой фазы, учитывающий движение твердого скелета, реологическое соотношение типа Максвелла, закон сохранения импульса системы в целом [1-4]. Здесь Рf, рз,vf з — соответственно истинные плотности и скорости жидкой и твердой фаз, Ф — пористость, pf ,рз — соответственно давления жидкой и твердой фаз, ре — эффективное давление, ргаг — общее давление, гаг — плотность двухфазной среды, д — плотность массовых сил; к(Ф) — коэффициент фильтрации, а1 (Ф) — коэффициент объемной вязкости. Плотности жидкой и твердой фаз считаются постоянными. Задача записана в эйлеровых координатах (х,Ь) е Qт.

Близкие по структуре системы рассматривались в работах [5-12]. В этих работах на основе ряда упрощающих предположений исходные системы сводились к одному уравнению высокого порядка. В [5] установлена локальная разрешимость задачи Коши пространствах С.Л. Соболева. В [6, 7] исследованы решения типа «простой волны». Вопросы обоснования моделей фильтрации в деформируемой пористой среде в модельных случаях исследованы в [13-16].

2. Одномерная задача. Рассмотрим одномерный случай и перейдем в этой системе уравнений к переменным Лагранжа [15]. Пусть х = х(т, х, Ь) — решение задачи Коши

Поскольку

„ дФ д - 9д'з

дх дх дх

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

1 дФ 1 д

+

1 - Ф дЬ 1 - Ф0 дх

1 д' з

(Ф(vf - 'з 3)) + т-ф0 Ф1х =0

Используя уравнение неразрывности для твердой фазы, получим

д 9 \ 1 д . 9 ,А А .. п Ц--Ф) +Т-Ф0дх(Ф(vf- 9)) = °.

Переходя от (х , Ь) к массовым лагранжевым переменным (у,Ь) по правилу (1 - Ф°(х))dX =

Х

dy, у(х) = 1(1 - ф0(n))dn е [0,1] и формально

0

заменяя у на х, приходим к следующей системе

+ - Ф)дх = (7)

й (1-Ф) + 1{Ф("> - П = 0- (8)

дРf

Ф' - '■) = -к(Ф)((1 - Ф)дх - Рf д), (9)

,, ,лдРгаг

(1 - Ф)^о1Г = -Ргагg,

дх

(10)

Положим 9 = х(0, х, Ь) и возьмем за новые переменные х и Ь. Тогда 1 - Ф(х,Ь) = (1 -

(1 - Ф) = -ач(Ф)Ре - а2(Ф)дЬ. (11)

'

■

Второе уравнение системы с учетом закона Дарси принимает вид

т^- дх(к(ф)((1 -ф)д£-= (12)

Далее рассмотрим случай, когда a2 (ф) (7) и (11) следует уравнение

0. Из

3. Численное исследование. В области Qт = [0,1] х [0,1] построим сетку с шагами Ь = 1/М,г = 1/М. Значение сеточной функции ф(х,£] в узлах сетки будем обозначать фп, i =

0,1,...,М; п = 0,1,...,М. Уравнение (16) представим в виде

дФ=-ф'ч — ф),

1 дф 1 - ф dt

= -а1(ф)(Рш - Pf ).

Это уравнение можно представить в виде

Ptot - Pf =

dG(ф)

дГ''

(13)

где функция G(ф) определяется равенством

dG _ 1

dф а1(ф)(1 — ф)

Уравнение (12) с учетом (10) и (13) перепишется в виде

д ( ф

dt \ 1 - ф

= IH - ф) д-Ш ) -

- — (к(ф)д(рм + Pf )). (14)

Перейдем к безразмерным переменным V = Т ,хХ = р, следовательно, область изменения хХ и V есть единичный отрезок [0,1], и так как а1 (ф) = фт/п,к(ф) = kфnто в уравнении (14) возникнут безразмерные параметры А = кп/^2 и е = gTLps/п и его можно представить в виде (штрихи опускаются)

1

дф

— ( ф ) = \— ( (ф) — dt (1 - ф' дх\а'[ф) дх \фт(1 - ф) dt

— Ь(ф) (15)'

где а(ф) = фп(1 — ф), Ь(ф) = ефпр/ра + (1 — ф) +

Для уравнения (15) рассмотрим следующую краевую задачу. Пусть в исходной области в переменных Эйлера на границах берутся условия Vs |ж=0,ж=1= Vf |ж=о,ж=1 = 0, тогда, учитывая (9), приходим к задаче для отыскания функций г и ф:

1

дф

фт(1 - ф) dt '

(16)

гЛ(ф) - дХ(а(ф) g - Ь(ф))=0, (17)

dz

(а(ф) дХ - Ь(ф)) Ix=o,x=i = 0, ф |t=o= ф°(х), (18)

где З.(ф) = фт/\(1 - ф).

аппроксимируем его двухэтапным методом Рун-ге - Кутта [17, с. 219].

Уравнение (17) системы аппроксимируем разностной схемой для однородного уравнения второго порядка с переменными коэффициентами [18, с. 153]. В уравнениях (7), (8) и (9) также перейдем к безразмерным переменным V = Т ,х' = р = и = и ,Р/ = р ,и = Тр, тогда уравнения для нахождения скоростей фаз и давления жидкой фазы с учетом (16) примут вид (штрихи опускаем):

dv

dvs дх

1 дф

гфт 1 - ф '

^ _ ф^V " ^) " *ф>

dpf _ £- + \^1-n(vs — vf) дх 1 — ф

где AT _ Tkp, £i _ - — безразмерные параметры.

Рассмотрим следующие значения параметров [19, с. 563]: м _ 1 Па • с, L _ 103 м, g _ 10 м/с2, T _ 105 с, pf _ 2, 5 • 103 кг/м3, ps _3 • 103 кг/м3, n _ 3,m _ 1. Вязкость горной среды ц принимает значения от 1013 Па • с до 1020 Па • с, проницаемость k может меняться от 10-19 м2 до 10-7 м2. Меняя значения этих двух величин, получим три варианта для параметров A и £. Численный расчет для пористости представлен на рис. 1. Когда g _ 0, решение имеет вид ф _ const и не зависит от выбора параметров A и £. Скорости фаз представлены на рис. 2. Давления pf ,ps и ptot построены с точностью до констант (см. рис. 3). Безразмерные параметры принимают значения A1 _ 104, £1 _ 2.5 • 102. График функции p'e = ape представлен на рис. 4, где a _ Tp.

Для получения точного решения рассмотрим тестовый пример. Пусть s _ -—ф. Далее полагаем

^^ФУ, следовательно,

а1(ф) = (1 - ф), тогда Щ

G _ -т^ф _ s. Уравнение (14) запишем в виде

ds д д2 s

- _ -тт^т-) ^)))

и рассмотрим случай, когда величины ф _ k^(s))(1 — ф(^)),ш _ k^(s))g(ptot + Pf) можно

z

Рис. 1. Изменение пористости при t = 1

Рис. 4. Эффективное давление при х = 1

Рис. 2. Скорости жидкой и твердой фазы при х = 1

считать постоянными. Тогда приходим к следующей задаче

дs д д'2s дЬ дх дхдЬ '

(19)

д21;

(^ьдт - Х=0д=о, s Ь=о= s0. (20)

Точное решение задачи (19)—(20) имеет вид

—х 2х 1

п Ьше ^ (е ^ — е ^)

s = s0 + у '

(е ^ +

(21)

Задача (19)-(20) решалась аналогично задаче (16)-(18). Численное решение с хорошей погрешностью сходится к точному решению на различных сетках.

4. Учет упругих свойств среды. Рассмотрим другой случай сведения системы типа (1)-(6) к одному уравнению высокого порядка. Далее полагаем а2(ф) = 0.

Проинтегрировав уравнение (12) от 0 до Ь, получим

ф дВ 1 — ф дх' где функция В имеет вид

Гх ^0

В

Ф0(.0

1 — Ф0(£)

t Цф)((1 — ф) дХ — рг дЦт

Тогда

Вг = Кф)((1 — ф) дх — Pf д), (22)

и с учетом (4), (10) уравнение (22) представим в виде

Вг = —Кф)((1 — ф) дре + д(ро + Pf)).

Выразив дг, получим

Рис. 3. Давления жидкой и твердой фаз и общее давление при х = 1

дре

Вг , д(ргог + Pf)

+

дх Щ)(1 — ф) 1 — ф

0

Уравнение (11) с учетом (7) представим в виде

= -о.1(ф)ре - a,2(ф)д-. (24)

Для простоты изложения ограничимся случаем о1(ф) = а1ф,о2(ф) = а2ф, (а1,а2) = const. Тогда уравнение (24) примет вид

Gif) дре

(25)

где

дG2(ф)

1

дф

в = Oi (ф)/02 (ф) .

d2G2W = вВг

dxdt к(ф)(1 - ф)

+ g(ptot + Pf) +

+

Bt

т\к(ф)(1 - ф)

+

1-ф

g(Ptot + Pf) 1 - ф

дЬ (1 - ф)а2(ф) дЬ

Продифференцировав уравнение (25) по х и воспользовавшись (23), получим уравнение для функции В(х,Ь) вида

где ф = Вх/(1 — Вх). Вопросам разрешимости начально-краевых задач для уравнения вида (26) посвящена книга [20].

5. Заключение. В работе проведено численное исследование одномерной нелинейной задачи фильтрации несжимаемой жидкости в вязкой пористой среде.

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

1. Bear J. Dynamics of Fluids in Porous Media // Elseiver. - New York, 1972.

2. Connoly J.A.D., Podladchikov Y.Y. Compaction-driven fluid flow in viscoelastic rock // Geodin. - Acta. - 11 (1998).

3. Morency S., Huismans R.S., Beaumont C, Fullsack P. A numerical model for coupled fluid flow and matrix deformation with applications to disequilibrium compaction and delta stability // Journal of Geophysical Redearch. — 112 (2007). -B10407.

4. Нигматулин Р.И. Динамика многофазных сред. — М., 1987. — Ч. 1.

5. Simpson M., Spiegelman M., Weinstein C.I. Degenerate dispersive equations arising in the stady of magma dynamics // Nonlinearty. — 20 (2007).

6. Abourabia A.M., Hassan K.M., Morad A.M. Analytical solutions of the magma equations for molten rocks in a granular matrix // Chaos Solutions Fract. — 42 (2009).

7. Geng Y., Zhang L. Bifurcations of traveling wave solutions for the magma equations // Applied Mathematics and computation. — 217 (2010).

8. Гоман В.А., Папин А.А., Шишмарев К.А. Численное решение двумерной задачи движения воды и воздуха в тающем снеге // Известия Алтайского государственного университета. — 2014. — № 1-2.

9. Шишмарев К.А. Тепломассоперенос в тающем снеге // Труды молодых ученых Алтайского государственного университета. — 2011. — № 8.

10. Токарева М.А., ВирцР.А. Аналитическое и численное исследование задачи фильтрации в пороупругой среде : cборник трудов Всероссийской конференции по математике. — 2016.

11. Байкин А. Н. Динамика трещины гидроразрыва пласта в неоднородной пороупругой среде: дис. ... физ.-мат. наук. — Новосибирск, 2016.

12. Dushin V.R., Nikitin V.F., Legros J.C., Silnikov M.V. Mathematical modeling of flows in porous media // WSEAS Transactions on Fluid Mechanics. — 2014. — T. 9.

13. Tokareva M.A. Solvability of initial boundary value problen for the equations of filtration poroelastic media // Journal of Physics: Conference Series. — 2016. — T. 722. — № 1.

14. Papin A.A., Tokareva M.A. Correctness of the initial-boundary problem of the compressible fluid filtration in a viscous porous medium // Journal of Physics: Conference Series. — 2017. — T. 894. — № 1.

15. Papin A.A., Tokareva M.A. On Local solvability of the system of the equation of one dimensional motion of magma // Журнал Сибирского федерального университета. Серия: Математика и физика. — 2017. — T. 10. — № 3.

16. Токарева М.А. Конечное время стабилизации уравнений фильтрации жидкости в по-роупругой среде // Известия Алтайского государственного университета. — 2015. — Т. 2. — № 1.

17. Самарский А.А., Гулин А.А. Численные методы. — М., 1989.

18. Самарский А.А. Теория разностных схем. — М., 1977.

19. Fowler A. Mathematical Geoscience. — Springer-Verlag London Limited, 2011.

20. Ларькин Н.А., Новиков В.А., Янен-ко Н.Н. Нелинейные уравнения переменного типа. — Новосибирск, 1983.