CC BY
14
УДК 517.53/.55
Численное исследование задачи фильтрации жидкости в тонком пороупругом слое*
М.А. Токарева, А.Н. Сибин
Алтайский государственный университет (Барнаул, Россия)
Numerical Study of a Problem of Fluid Filtration in a Thin Poroelastic Layer
M.A. Tokareva,A.N. Sibin
Altai State University (Barnaul, Russia)
Излагаются результаты численного исследования математической модели фильтрации вязкой жидкости в пороупругой среде, обладающей вязкоупругими свойствами. Основное внимание уделяется выводу модели, постановке задачи, разработке численного алгоритма решения поставленной задачи, а также анализу предварительных численных результатов. Рассматриваемую в работе модель можно использовать при изучении процессов, происходящих в ледяном покрове. В данном подходе лед представляет собой двухфазную среду, состоящую из твердой фазы (лед) и жидкой (вода). В качестве твердой фазы рассматривается твердый деформируемый ледяной скелет, обладающий вязкоупруги-ми свойствами. Таким образом, в рассматриваемой модели ледяному покрову приписываются свойства неньютоновской жидкости. Данная модель не учитывает фазовые переходы и изменение температуры. В процессе обезразмерива-ния исходной системы уравнений вводится малый параметр (по времени). После предельного перехода по параметру (рассмотрен случай медленных процессов) система уравнений характеризует твердый скелет как среду, обладающую больше упругими свойствами, чем вязкими. Для полученной системы уравнений проведены тестовые численные расчеты, определены поле скоростей, пористость и критическое напряжение.
Ключевые слова: многофазная фильтрация, по-роупругость, ледовый покров, конечно-разностная схема, вязкоупругость.
The paper presents results of a numerical study of a mathematical model of viscous fluid filtration in a poroelastic medium with viscoelastic properties. The focus of this research is on model development, problem formulation, and elaboration of a numerical algorithm to solve this problem, as well as a preliminary analysis of numerical study results. The proposed model can be used in a study of processes that occur in the ice cover. This approach treats the ice as a biphasic medium consisting of a liquid (water) phase and a solid (ice) phase being a solid elastic ice skeleton with viscoelastic properties. Thus, the ice cover has properties of a non-Newtonian fluid in this model, and phase transitions and temperature changes are out of concern. A small time parameter is introduced for the process of nondimensionalizing of the original equation system. After passing to the limit (for slow processes), the equation system describes the solid skeleton as a medium with elastic properties greater than viscous. Test numerical calculations are performed, and the field of velocities, porosity, and critical stress values are obtained.
Key words: multiphase flow, poroelasticity, sheet ice, finite-difference scheme, viscoelasticity.
DOI 10.14258/izvasu(2017)1-25
Уравнения модели. Для каждой составляющей двухфазной среды (ледяного скелета в и содержащейся в нем жидкости /) вводятся понятия объемов твердого скелета У и пор Ур. Удельный объем пор (пористость)
* Работа выполнена при финансовой поддержке гранта РФФИ №16-08-00291.
Ф
V, где Vt
Ур + У - общий объем. Законы сохранения масс для жидкости и твердой фазы в отсутствие фазовых переходов выглядят следующим образом [1]:
d(Pf ф) dt
+ V • (Pf №f) = 0,
(1)
9(1 ~ Ф)Ре dt
+ V • ((1 - Ф)Рзйз
0,
(2)
где t — время, pf — плотность жидкости, рь — плотность твердого скелета, Vf,УУЬ — скорости жидкой и твердой фаз соответственно, V = (дХ, д, тЬ), (х, у, г) — переменные Эйлера.
Уравнение сохранения импульса для жидкости берется в форме закона Дарси [2, 3]:
Ф(Vf - Vs) = -К№)(Vpf + Pfg),
(3)
где К — коэффициент фильтрации, pf — давление жидкости, д — плотность массовых сил.
В зависимости от механических свойств льда используются различные подходы, определяющие его свойства. Реологическое соотношение, определяющее вязко-упругие свойства пористой деформируемой среды, имеет вид [3—9]
divvs
— а1 (ф)Ре — а-2(ф) <dp, (4)
где ре = (1 — ф)(рв — Pf) — эффективное давление, ах(ф) — коэффициент, определяющий объемную вязкость среды, а2(ф) — коэффициент, определяющий объемную сжимаемость среды, ^ = д + ■ V — материальная производная.
Для замыкания используется уравнение равновесия всей системы в целом, в предположение, что вязкость жидкости мала [9]:
рд + ¿гу((1 — ф)у(дХ + (дХ)*)) — = (5)
где р = (1 — ф)р.э + фрf — средняя плотность среды, а = аз(1 — ф) + o^fф — общее напряжение, определенное как средневзвешенное, р0 = фpf + (1 — ф)ра — общее давление, V — динамическая вязкость твердой среды.
Таким образом, в области О = (х, г) = [0, Ь] х [0,Н] рассматривается система уравнений (1)—(5) составного типа. В случае преобладания упругих свойств скелета локальная разрешимость по времени начально-краевой задачи для данной модели доказана в [10], свойство конечной скорости распространения возмущений установлено в [11]. В работе [12] рассмотрены различные режимы движения в зависимости от поведения возникающего в задаче малого параметра. Подобные модели были рассмотрены в работах [13, 14, 15, 16].
Введение малого параметра. Предельный переход. Проведем обезразмеривание уравнений (1)—(5). Пусть Х,г,~Е — безразмерные переменные, определенные равенствами
V
z k H
V = H, t = £kТоt, £ = V < 1,
Положим:
Pf (t,x,z) = apf (t,x,g) = ap f (£kTot, )>
Ps(t,x,z) = ap s(t,x, t) = aps(£kTot,X,HH), vi(t, x, z) = ßiVis(t, x, z) = ßiVis(£kTot, X, H), i = !, 2, vf (t, x, z) = ß'vf (g, x, z) = ß'gf (£kTot, ),i = 1, 2,
Ptot(t, x, z) = ap tot(t , x, z),
a__ a__ a
Pfg = HPfv, psg = Hpsg, P9 = h vv,
K (Ф) = koK (ф), ах(ф) = 01а1(ф), а2(ф) = a2g2 (ф). Здесь [ß>] = [м/с], [a] = [Па], [ko] = [П^],
[a1] = [пЬ], [a2] = \1/Па].
Для получения безразмерной формы уравнений следует положить
ß1 = £k toL, ß2 = £k toH.
После предельного перехода (£ ^ 0 ) при k = —2 система (1)—(5) примет вид
+ div((1 — ф)уа)=0,
dt
дф dt
+ div^v f) = 0,
vi = vj,
(6) (7)
,dp f
ToV-Фv — v2) = —К (ф)(д- — Pf v) (8)
dvS dvS 2 d(ptot — p f) , ,
dP + dV = —ааа2(Ф)-dt—, (9)
dv ((1—Ф)%) ^
\ dttl \ d
(10)
2dt ((1 — Ф)#) + dl ((1 — Ф)di) =0. (11)
Тензор напряжений пороупругой среды в этом случае имеет вид
0
a =, (1 — Ф) 9-dZ
\(1 — Ф) a-dZ 2(1 — Ф) dz
Максимальное нормальное напряжение дается формулой
dv2 / dv2 dv1
= (1 - ds W< ds)2 + (-ss)2
Модельная задача. Как частный случай рассмотрим задачу для уравнений (6), (10), (11), дополненных следующими начально-краевыми условиями:
dv2
(! - Ф)Iz |z=o= 0, vS \z=h = L = const,
z=o
где [L] = [H] = [м], [т0] = [1/с]; к - произвольное вещественное число.
=H = B = const, (1 — Ф) Ф |t=o = Ф0(Х,z), d§f \Z=H =0,
dvS dz
0,
Pf \t=0= Pf (x,z),Ps \t=0= p°o(X,z),Pf \ z=0= P0(x,t)
X
1
v
s \z
Тогда, пользуясь краевыми условиями, находим
ф\г=0 = ф0(х,г).
= ь,
а уравнение для в перепишется в виде
дв ^дв ьдв 0
дЬ дх дг
Его решение есть
ф = ф0(х - ВЬ,г - ЬЬ).
Тогда
Введем равномерную сетку с узлами х^ = ¿Нъ { = 0...^; г^ = ]Н2, ] = 0...Н2; Ьп = пт, п = 0...Т; Н\ = 1/М1 и Н2 = 1/^2 - шаги по пространственным координатам, т — шаг по времени. В данной работе используются обозначения, принятые в [18].
2 аК(ф0(х - ВЬ, -ЬЬ)) т = 9 г0Ь2ф-" + Ь,
Решения для р£ ,р8 выглядят следующим обра-
зом:
х
РГ = ро + Р£9(г - У
К(ф0(х - ВЬ,Н - ЬЬ)) К(ф0(х - ВЬ, г - ЬЬ))
с1т),
р8 = р£ - р£(х - ВЬ, г - ЬЬ) + р0(х - ВЬ, г - ЬЬ).
Алгоритм численного решения задачи.
Пусть ф = 1 - ф, и = и V = тогда уравнения (6)-(11) примут вид
дф дфи дф 0 дЬ дх дг
д
ди
дг \ф дг
2 д ( д дг\фдг) дх \ф дг
ди
0.
(12)
(13)
(14)
Заметим, что уравнения (12), (13) и (14) образуют замкнутую систему относительно искомой функции ф и компонент вектора скорости скелета = (и,^).
Основную трудность для численного решения системы (12)-(14) представляют сильные разрывы решения нелинейного гиперболического уравнения (12). Эффективный прием расчета задач с разрывными решениями изложен в работе [17 с. 359]. Вместо гиперболического уравнения (12), решение которого рвется, численно решается параболическое уравнение с псевдовязкостью еАф
дф _ л ,ди дv. . .
д + ^ -УФ - ^ = -Ф( д* + Тг], (15)
где коэффициент е мал.
На границах области Г = (х, г) = [0,1] х [0,1] и в начальный момент времени заданы условия
дф дх
= 0, и(х, 0, Ь) = и0(х, Ь), и(х, 1,Ь) = и1(х, Ь),
Рис. 1. Пористость
Рис. 2. Максимальное нормальное напряжение
v(x, 0,Ь) = v0(x,Ь),v(x, 1,Ь) = v1(x,Ь).
1
1
vs = v
£
г
Алгоритм вычисления следующий: сначала, используя начальное значение пористости, методом прогонки по г находятся поле скоростей из уравнений (13) и (14). Затем на известном поле скоростей методом переменных направлений определяется значение пористости из уравнения (15). Значения на п+1 слое находятся методом прогонки по х, а значения на п + 1 слое прогонкой по г.
В численных расчетах использовался следующий набор модельных параметров:
и0 = 2, и1 = 1, v0 = 0, v1 = 0, ф0 = 0.7.
На рисунке 1 представлено распределение пористости по оси х при г = 0.5. Пористая среда сжимается (пористость стремится к нулю), пока
напряжения не равны нулю. Максимальное нормальное напряжение уменьшается с течением времени (см. рис. 2).
Другие начально-краевые задачи уравнений динамики ледового покрова были численно исследованы в [19, 20].
Заключение. В работе рассмотрена математическая модель фильтрации жидкости в поро-упругом тонком слое. Предложен алгоритм численного решения задачи. Проведены тестовые численные расчеты.
Авторы статьи признательны С.С. Кузикову и А.А. Папину за обсуждение задачи и конструктивные замечания.
Библиографический список
1. Антонцев С.Н., Кажихов А.В., Монахов В.Н. Краевые задачи механики неоднородных жидкостей. — Новосибирск, 1983.
2. Bear J. Dynamics of Fluids in Porous Media. - Elsevier, New York, 1972.
3. Connolly J.A.D., Podladchikov Y.Y. Compaction-driven fluid flow in viscoelastic rock, Geodin. Acta. - 1998. - № 11.
4. Domenico P.A. and Schwartz F.W., Phisical and Chemical Hydrogeology, Jhon Wiley. - New York, 1990.
5. Fowler A.C., Yang X. Pressure solution and viscous compaction in sedimentary basins // J. Geophys. Res. - 104, 12,989-12,997, 1999.
6. Birchwood R.A., Turcotte D.L. A unified approach to geopressuring, low-permeability zone formation, and secondary porosity generation in sedimentary basins // J. Geophys. Res. 99, 20,05120,058, - 1994.
7. Fowler A.C. A compaction model for melt transport in the Earth's asthenosphere, part 1, the basic model, in Magma Transport and Storage, edited by M.P. Ryan, Jhon Wiley. - New York, 1990.
8. Mc.Kzenzie D.P. The generation and compaction of partial melts // J. Petrol., 25, 713-765, 1984.
9. Morency C., 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 Research. - 2007. - Vol. 112.
10. Tokareva M.A. Solvability of initial boundary value problem for the equations of filtration in poroelastic media // Journal of Physics: Conference Series. - 2016. - Т. 722, № 1.
11. Tokareva M.A. Localization of solutions of the equations of filtration in poroelastic medium // Журнал Сибирского федерального ун-та. Серия: Математика и физика. - 2015. - Т. 8, № 4.
12. Токарева М.А. Двумерная задача фильтрации в тонком пороупругом слое // Известия Алтайского гос. ун-та. - 2013. - № 1-1 (77).
13. Папин А.А., Подладчиков Ю.Ю. Изотермическое движение двух несмешивающихся жидкостей в пороупругой среде // Известия Алтайского гос. ун-та. - 2015. - Вып. 1/2. DOI:10.14258/izvasu(2015)1.2-24
14. ПапинА.А., Сибин А.Н. Автомодельное решение задачи поршневого вытеснения жидкостей в пороупругой среде // Известия Алтайского гос. ун-та. - 2016. - № 1 (89). DOI: 10.14258/izvasu(2016)1-27
15. Папин А.А., Сибин А.Н. О разрешимости первой краевой задачи для одномерных уравнений внутренней эрозии // Известия Алтайского гос. ун-та. - 2015. - Вып. 1/2 (85). DOI:10.14258/izvasu(2015)1.2-25
16. Папин А.А. Существование решения «в целом» уравнений одномерного неизотермического движения двухфазной смеси. I. результаты о разрешимости // Сибирский журнал индустриальной математики. - 2006. - Т. IX, № 3.
17. Калиткин Н. Н. Численные методы. — М., 1978.
18. Самарский А.А. Введение в теорию разностных схем. - М. 1971.
19. Shishmarev K., Khabakhpasheva T., Korobkin A. The response of ice cover to a load moving along a frozen channel. Applied Ocean Research 59 (2016) 313-326.
20. Шишмарев К.А. Постановка задачи о вязкоупругих колебаниях ледовой пластины в канале в результате движения нагрузки. Известия Алтайского гос. ун-та. - 2015. - Вып. 1/2 (85). D0I:10.14258/izvasu(2015)1.2-35
