Оценка максимального перемещения конечного объема легкой жидкости в насыщенной пористой среде

Полный текст

Климатические изменения стимулируют исследование процессов, сопровождающих размещение углекислого газа (CO₂) в проницаемых геологических пластах [1]. Закачка в недра парниковых газов, в частности CO₂, снижает воздействие человека на окружающую среду и, в перспективе, способствует замедлению эффекта глобального потепления. Водонасыщенные пласты являются перспективными объектами для геологического хранения CO₂, так как они повсеместно распространены в осадочных бассейнах. Выбор наиболее перспективного объекта сопряжен с оценкой различных показателей хранения, например, емкости — объема газа, который может быть размещен в пласте, или приемистости — максимального темпа нагнетания газа. Также важна оценка безопасности размещения CO₂, т.е. целостности и герметичности объекта хранения [2, 3].

При закачке CO₂ в пласт важно учитывать эффект плавучести газа, представляющего собой более легкую фазу. Под действием силы Архимеда CO₂ стремится перетечь в зоны меньших глубин. Газ всплывает в водонасыщенной пористой среде. Если пласт имеет структурный наклон, то сила Архимеда приводит к течению газа вдоль пласта в направлении уменьшения глубины. В результате CO₂ может переместиться на значительное расстояние от скважины. Это напрямую влияет на безопасность хранения, так как в направлении течения могут быть геологические разломы и заброшенные скважины, представляющие собой возможные пути утечки CO₂. В этой связи актуальна оценка максимального расстояния, на которое CO₂ переместится от скважины [4–6].

В работе [6] проведено обширное параметрическое исследование максимального перемещения CO₂ в наклонном пласте. Методом Монте-Карло рассчитывались различные сценарии закачки CO₂ с учетом многих осложняющих факторов, как, например, гистерезис относительной фазовой проницаемости [7]. Показано, что максимальное перемещение CO₂ от скважины определяется критерием подобия, характеризующим силу Архимеда, и разностью критических насыщенностей газа в процессах пропитки и дренажа. Эта разность фактически характеризует эффект гистерезиса.

Настоящая работа направлена на подробное исследование влияния гистерезиса на максимальное перемещение CO₂. Цель работы — определить критерии подобия, имеющие наибольшее влияние, и разработать методы оценки максимального расстояния, на которое газ распространится от скважины. Для того чтобы упростить исследование и сконцентрироваться на влиянии гистерезиса, рассмотрена упрощенная одномерная постановка задачи с осью x, направленной вдоль наклонного пласта (рис. 1). Предполагается, что газ и вода являются несмешивающимися несжимаемыми жидкостями, далее обозначаемыми символами g и l. Закачка CO₂ приводит к образованию локализованной зоны, насыщенной газом. При этом объем размещенного в пласте газа конечен. Считая, что такое состояние существует при t = 0, исследуется течение CO₂ вдоль пласта (т.е. в направлении оси x) и оценивается его максимальное перемещение x_max.

 

Рис. 1. Схема распределения насыщенности легкой жидкости g при t = 0 и t → ∞. Сила тяжести действует в противоположном оси x направлении.

 

Сделанное предположение о несжимаемости жидкостей не является ограничительным для приложений, связанных с геологическим хранением CO₂. Во-первых, коэффициент сжимаемости воды мал, а, во-вторых, в пластовых условиях (при давлении более 100 бар) CO₂ — сверхкритическая жидкость, имеющая плотность 500–800 кг/м³ [1]. Коэффициент сжимаемости этой жидкости значительно меньше, чем у совершенного газа. Предполагается, что газ остается в пласте на глубинах более 1000 м, где он находится в сверхкритическом состоянии, поэтому сделанное предположение будем считать не слишком грубым.

Предположение о несмешивающейся фильтрации конечно же можно назвать сильным, учитывая, что 10–20% закачанного CO₂ растворяется в воде [3]. Однако подробное исследование влияния растворения на максимальное перемещение, проведенное нами ранее в [5, 6], не является целью настоящей теоретической работы. Еще одним упрощающим допущением рассматриваемой модельной постановки задачи является предположение о равенстве фазовых давлений в жидкости и газе. Фактически рассматривается фильтрация в пределе бесконечно больших чисел Бонда. Оценка влияния капиллярного давления может являться предметом отдельных исследований. Предварительные результаты показывают, что при уменьшении числа Бонда максимальное перемещение газа уменьшается.

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим фильтрацию в бесконечной вертикальной колонне, наполненной пористой средой. Предполагаем, что пористость и проницаемость имеют однородное распределение и не изменяются со временем. Введем координатную ось x, направленную вдоль колонны вверх против направления действия силы тяжести (рис. 1). В начальный момент времени, t = 0, пористая среда при 0 ≤ x ≤ L насыщена легкой жидкостью g (газом). При x < 0 и x > L среда насыщена тяжелой жидкостью l (водой). Полагаем, что на значительном удалении от зоны 0 ≤ x ≤ L, т.е. при |x / L| >> 1, жидкость l покоится. Это может достигаться установкой при |x / L| >> 1 непроницаемых стенок, препятствующих течению жидкости вдоль оси x.

Предполагается, что обе жидкости l и g — несжимаемые несмешивающиеся среды, характеризующиеся постоянными значениями плотности ρ_l и ρ_g, причем ρ_l > ρ_g. Также жидкости характеризуются постоянными динамическими вязкостями μ_l и μ_g.

Так как ρ_l > ρ_g, то на жидкость g действует сила Архимеда. Следовательно, при t > 0 она начнет перемещаться вверх. Это движение будет сопровождаться капиллярным захватом жидкости g в виде неподвижных пузырьков при их относительно малом критическом объемном содержании в поровом пространстве. Таким образом, подвижный объем фазы g будет уменьшаться со временем. В результате при t → ∞ вся жидкость g будет обездвижена, так как ее насыщенность при любом x будет меньше или равна критическому значению. При этом g достигнет некоторой максимальной координаты x, которую дальше обозначим x_max. Цель настоящей работы заключается в получении простых оценок для x_max при различных параметрах жидкостей и пористой среды.

2. ОСНОВНЫЕ УРАВНЕНИЯ

Задачу рассматриваем в одномерном приближении, предполагая, что при каждом x параметры течения имеют одинаковые значения в поперечном сечении колонны. Также предположим, что давления в обеих жидкостях одинаковые. Тогда уравнения, описывающие двухфазную фильтрацию жидкостей, принимают вид [8-10]:

$\frac{\partial }{\partial t}\left(\varphi {\rho }_i{s}_i\right)+\frac{\partial }{\partial x}\left({\rho }_i{u}_i\right)=\mathrm{0,} i=g, l$; (1)

$u_i=-k\frac{k_{ri}}{{\mu }_i}\left(\frac{\partial p}{\partial x}+{\rho }_{i}g\right)$; (2)

$s_g+s_l=1$, (3)

где φ — пористость, ρ — плотность, s — насыщенность, u — скорость фильтрации, k — абсолютная проницаемость, k_r — относительная фазовая проницаемость, μ — динамическая вязкость, p — поровое давление, g — проекция ускорения свободного падения на направление оси x, а индексами g и l обозначены параметры фаз g и l. Уравнения (1) — законы сохранения массы, записанные для каждой жидкости, (2) — закон фильтрации Дарси, a (3) — замыкающее соотношение для насыщенностей.

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

$\begin{array}{l}{s}_g=s_{g,\max }\left(0\le x\le L\right),\\ s_g=0\left(x<0\cup x>L\right),\end{array}$ (4)

где s_g,max — максимальная насыщенность фазы g.

Условие отсутствия потока при больших x записывается в виде:

$\left|x/L\right|\gg 1:u_i=0$. (5)

Характерный вид кривых относительной фазовой проницаемости приводится на рис. 2. Предполагается, что кривая k_rl (s_g) не зависит от направления изменения насыщенности s_g. Фазовая проницаемость k_rl равна 1 при s_g = 0. Она убывает с ростом s_g, обращаясь в ноль при s_g = s_g,max = 1 − s_lc, где s_lc — критическая насыщенность жидкости l [9]. Отметим, что согласно (4), при t = 0 насыщенность s_g не превышает максимального значения s_g,max.

 

Рис. 2. Схема кривых относительной фазовой проницаемости для условий дренажа и пропитки.

 

Фазовая проницаемость k_rg зависит от направления изменения насыщенности s_g, т.е. она проявляет гистерезисное поведение. Предположим, что при заданном x s_g = 0 и далее s_g начинает расти, т.е. полагая, что g — газ, происходит процесс дренажа. В этом случае $k_{rg}\left(s_g\right)=k_{rg}^{\left(d\right)}\left(s_g\right)$ (рис. 2). Если s_g меньше критического значения $s_{gc}^{\left(d\right)}$, то $k_{rg}^{\left(d\right)}=0$, т.е. жидкость g не движется. При $s_g>s_{gc}^{\left(d\right)}$, $k_{rg}^{\left(d\right)}$ растет с s_g, достигая максимального значения k_rg,max при s_g = s_g,max. Если далее процесс дренажа сменяется на процесс пропитки, т.е. s_g уменьшается, то $k_{rg}=k_{rg}^{\left(i\right)}$. Кривая пропитки $k_{rg}^{\left(i\right)}\left(s_g\right)$ лежит ниже кривой дренажа $k_{rg}^{\left(d\right)}\left(s_g\right)$. Фазовая проницаемость $k_{rg}^{\left(i\right)}$ обращается в ноль при $s_g=s_{gc}^{\left(i\right)}\ge s_{gc}^{\left(d\right)}$. Таким образом, если s_g достигла значения s_g,max, то она уже не может опуститься ниже значения $s_{gc}^{\left(i\right)}$.

Фактически кривые $k_{rg}^{\left(d\right)}\left(s_g\right)$ и $k_{rg}^{\left(i\right)}\left(s_g\right)$ соответствуют максимальной петле гистерезиса, в которой s_g достигает значения s_g,max. Предположим, что в процессе дренажа s_g достигает меньшего значения, которому на рис. 2 соответствует точка I при s_g = s_g,hy. Если далее дренаж сменяется пропиткой, то  соответствует кривая DI, расположенная между $k_{rg}^{\left(d\right)}\left(s_g\right)$ и $k_{rg}^{\left(i\right)}\left(s_g\right)$. Уменьшению s_g соответствует движение от I к D, где точка D определяет критическую насыщенность s_gc, а точка I — исторически максимальную насыщенность s_g,hy. Если пропитка сменяется дренажем, то k_rg возрастает вдоль отрезка DI пока s_g не достигнет значения s_g,hy, а k_rg (s_g) не совпадет с $k_{rg}^{\left(d\right)}\left(s_g\right)$. При дальнейшем росте s_g выполняется равенство $k_{rg}=k_{rg}^{\left(d\right)}$. Если в чередующихся процессах пропитки и дренажа достигается следующее максимальное значение, которому, например, соответствует точка I′, то процессу пропитки соответствует кривая D′I′.

Таким образом, гистерезис k_rg определяется кривыми $k_{rg}^{\left(d\right)}\left(s_g\right)$ и $k_{rg}^{\left(i\right)}\left(s_g\right)$ и исторически максимальным значением s_g,hy, достигнутым к данному моменту времени. Для математического описания гистерезиса k_rg используется модель [7], которая основывается на следующих соотношениях:

$s_{gc}=s_{gc}^{\left(d\right)}+\frac{s_{g,hy}-s_{gc}^{\left(d\right)}}{1+C(s_{g,hy}-s_{gc}^{\left(d\right)})}, C=\frac{1}{s_{gc}^{\left(i\right)}-s_{gc}^{\left(d\right)}}-\frac{1}{s_{g,\max }-s_{gc}^{\left(d\right)}}$; (6)

$k_{rg}\left(s_g\right)=k_{rg}^{\left(i\right)}\left(s_g^{*}\right)\frac{k_{rg}^{\left(d\right)}\left(s_{g,hy}\right)}{k_{rg,\max }}\left(s_g<s_{g,hy}\right), s_g^{*}=s_{gc}^{\left(i\right)}+(s_{g,\max }-s_{gc}^{\left(i\right)})\frac{s_g-s_{gc}}{s_{g,hy}-s_{gc}}$, (7)

где уравнения (6) определяют критическую насыщенность s_g, т.е. абсциссу точки D на рис. 2, а уравнения (7) — кривую DI.

В настоящей работе положим, что

$$\begin{gathered} k_{rg}^{\left(j\right)}=k_{rg,\max }{\left(\frac{s_g-s_{gc}^{\left(j\right)}}{s_{g,\max }-s_{gc}^{\left(j\right)}}\right)}^{n_g}\left(s_g\ge s_{gc}^{\left(j\right)}\right) \\ {k}_{rl}={\left(\frac{s_{g,\max }-s_g}{s_{g,\max }}\right)}^{n_l} \end{gathered}$$

где n_g и n_l — константы [11], а .

3. УРАВНЕНИЯ В БЕЗРАЗМЕРНОМ ВИДЕ

Выберем L в качестве характерного масштаба длины. Так как фильтрация жидкости g происходит под действием силы Архимеда, то в качестве масштаба давления естественно взять (ρ_l − ρ_g) gL. Тогда u_scl = (k / μ_l) (ρ_l − ρ_g) g есть характерная скорость фильтрации (2). Масштаб u_scl характеризует интенсивность течения, вызванного силой Архимеда. Так как истинная скорость течения жидкости характеризуется u_scl / φ (1 − s_lc), то в качестве масштаба времени естественно выбрать φ (1 − s_lc) L / u_scl.

В соответствии с выбранными масштабами введем безразмерные переменные

$\begin{array}{l}x=\widehat{x}L, t=\widehat{t}\varphi (1-s_{lc})L/u_{scl}, u=\widehat{u}{u}_{scl},\\ p+{\rho }_{l}gx=\widehat{p}\left({\rho }_l-{\rho }_g\right)gL, s_i=(1-s_{lc}){\widehat{s}}_i, k_{rg}=k_{rg,\max }{\widehat{k}}_{rg,}\end{array}$ (8)

где «крышкой» обозначены безразмерные и масштабированные величины. Выбор 1 − s_lc в качестве масштаба насыщенности обусловлен тем, что доля s_lc порового пространства всегда остается насыщенной неподвижной связной фазой l. Эта доля фактически не влияет на течение и, таким образом, может быть исключена из рассмотрения с помощью введения масштабного коэффициента.

Подставляя (8) в уравнения (1) — (5) и для упрощения обозначений опуская у безразмерных переменных знак «крышки», получим

$\frac{\partial s_i}{\partial t}+\frac{\partial u_i}{\partial x}=\mathrm{0,} i=g, l$; (9)

$u_g=-M{k}_{rg}\left(\frac{\partial p}{\partial x}-1\right), u_l=-k_{rl}\frac{\partial p}{\partial x}$; (10)

$t=0:s_g=1\left(0\le x\le 1\right), s_g=0\left(x<0\cup x>1\right)$; (11)

$\left|x\right|\gg 1:s_g=\mathrm{0,} \frac{\partial p}{\partial x}=0$; (12)

$M=k_{rg,\max }\frac{{\mu }_l}{{\mu }_g}$, (13)

где М — отношение подвижностей жидкостей.

Складывая уравнения (9) и учитывая, что сумма насыщенностей равна единице, получим

$\frac{\partial }{\partial x}(u_g+u_l)=0$. (14)

Подставляя (10) в (14) и интегрируя получившееся соотношения с учетом граничных условий (12), получим

$\frac{\partial p}{\partial x}=\frac{M{k}_{rg}}{M{k}_{rg}+k_{rl}}$. (15)

Подставляя (15), (10) в (9) для i = g, получим одно уравнение на насыщенность s_g

$\frac{\partial s_g}{\partial t}+\frac{\partial G}{\partial x}=0$, (16)

$G(s_g,s_{g,hy})=\frac{M{k}_{rg}\left(s_g,s_{g,hy}\right)k_{rl}\left(s_g\right)}{M{k}_{rg}\left(s_g,s_{g,hy}\right)+k_{rl}\left(s_g\right)}$, (17)

$s_{g,hy}(t,x)=\underset{t^{\text{'}}\le t}{\max }{s}_g(t^{\text{'}},x)$. (18)

Таким образом, исследование фильтрации сводится к решению одного уравнения гиперболического типа (16), в котором функция потока G (17) зависит не только от значения s_g в данный момент времени t, но и от исторически максимального значения насыщенности s_g,hy (18). При s_g = s_g,hy имеем кривую G (s_g, s_g) для процесса дренажа, а при s_g,hy = 1 кривая G (s_g, 1) описывает процесс пропитки (рис. 3).

 

Рис. 3. Типичный вид функции G (s_g, s_g,hy) при s_g = s_g,hy (дренаж) и s_g,hy = 1 (пропитка). Используются масштабированные в соответствии с (8) значения s_g.

 

Решение уравнения (16) при начальных условиях (11) может быть получено методом характеристик, скорость которых пропорциональна производной ∂G / ∂s_g при s_g,hy = const в случае пропитки или s_g,hy = s_g в случае дренажа [8, 9]. Зависимость G от s_g,hy существенно осложняет процесс аналитического решения задачи. В этой связи дальше построим приближенное решение задачи и оценим его точность с помощью сравнения с численным решением уравнения (16) при начальных условиях (11).

4. ПРИМЕР ТЕЧЕНИЯ

Рассмотрим результаты прямого численного моделирования течения при M = 3.7, n_l = n_g = 2, $s_{gc}^{\left(d\right)}=0.125$ и $s_{gc}^{\left(i\right)}=0.375$. В соответствии с (11), при t = 0 вся жидкость g находится при 0 ≤ x ≤ 1, причем ее насыщенность равна 1. Таким образом, безразмерный объем жидкости

$V=\underset{-\infty }{\overset{+\infty }{\int }}{s}_{g}dx=1$ (19)

равен единице при t = 0. Так как газ моделируется несжимаемой средой, интеграл (19) не изменяется со временем. Следовательно площадь под графиком функции s_g (x) при любом t равна 1.

При 0 < x << 1 изменения s_g со временем сконцентрированы при x = 1. Эти изменения происходят из-за распада разрыва, который при t = 0 находится при x = 1. В положительном направлении оси x распространяется сильный разрыв S, на котором s_g скачком возрастает от 0 до некоторого значения $s_g^{\text{'}}$ (рис. 4, t = 2.2). Причем, при выбранных параметрах течения $s_g=s_g^{\text{'}}=const$ во всей области, простирающейся от x = 1 до S. В отрицательном направлении оси x от x = 1 распространяется последовательность волн, состоящая из сильного разрыва S_r и присоединенной к нему волны Римана. Таким образом, сила Архимеда приводит к распаду начального разрыва при x = 1 за счет течения легкой жидкости g вверх, а тяжелой жидкости l вниз.

 

Рис. 4. Рассчитанные распределения s_g (x) в последовательные моменты времени. При t ≥ 100 параметры течения практически не изменяются со временем. Жирная кривая показывает распределение s_g,hy при t >> 1.

 

Начальный разрыв при x = 0 устойчивый, так как на нем s_g возрастает с x, т.е. тяжелая жидкость расположена под легкой.

В соответствии с общими принципами графического решения задачи о распаде произвольного разрыва в среде, описываемой уравнением (16), необходимо на плоскости {s_g, G} построить огибающую к функции G [9]. В рассматриваемом случае нахождение решения осложняется тем, что G зависит не только от s_g, но и от s_g,hy. Опуская детали, опишем решение задачи при x ∼ 1, которое является автомодельным относительно переменной ξ = (x − 1) / t. Состоянию сверху и снизу от начального разрыва на плоскости  соответствуют точки U и L при s_g = 0 и s_g = 1 (рис. 3).

При малых t и 0 ≤ x ≤ 1 происходит пропитка, а s_g,hy = 1. Следовательно, течение в этой области описывается кривой G (s_g, 1). Разрыву S_r соответствует касательная к G (s_g, 1), проведенная из точки L (рис. 3). Точка касания R соответствует состоянию за разрывом S_r. Отрезок RA кривой G (s_g, 1) соответствует присоединенной к S_r волне Римана, которая простирается от S_r до x = 1. В точке A при  характеристическая скорость ∂G / ∂s_g равна нулю.

При малых t и x > 1 происходит процесс дренажа, а s_g,hy = s_g. Следовательно, течение в этой зоне описывается кривой G (s_g, s_g). Разрыву при x = 1 и t > 0 соответствует переход с линии G (s_g, 1) из зоны пропитки при x ≤ 1 на линию G (s_g, s_g) в зоне дренажа при x > 1. Причем, так как скорость этого разрыва равна нулю, то на плоскости {s_g, G} ему соответствует горизонтальный отрезок AI (рис. 3). Разрыву S, распространяющемуся в область x > 1, соответствует отрезок UI. Таким образом, решение задачи о распаде произвольного разрыва описывается кривой LRAIU, огибающей LRA и UI.

Описанная волновая конфигурация существует только при малых t, до тех пор пока разрыв S_r не достигнет границы x = 0. После взаимодействия с этой границей от x = 0 вверх распространяются отраженные волны, которые догоняют разрыв S (рис. 4, t = 4.2 и 5.7). В зависимости от вида функции G (s_g, s_g,hy), в частности от расположения ее точек перегиба, это могут быть различные последовательности сильных разрывов и волн Римана. Их взаимодействие с S приводит к снижению s_g за разрывом S. Таким образом, s_g,hy тоже опускается ниже s_g у верхнего края зоны, насыщенной жидкостью g.

С ростом t насыщенность s_g при каждом x за передним фронтом газа S снижается, приближаясь к критическому значению s_gc, а скорость фильтрации стремится к нулю. Скорость разрыва S также снижается по мере того как S достигает максимальной координаты x_max. Фактически при t ≥ 100 распределение s_g можно уже считать не зависящим от t (рис. 4, t = 100). В значительном диапазоне x, исключая малую зону за фронтом S, это распределение кусочно-постоянное. При 0 ≤ x ≤ 1 $s_g=s_{gc}^{\left(i\right)}$, а при x > 1 $s_g=s_{gc}^{\text{'}}$.

Заметим, что последовательности волн, в том числе и в задаче о распаде произвольного разрыва при x ∼ 1 и t << 1, могут отличаться от описанных. Например, при $s_{gc}^{\left(i\right)}\approx s_{gc}^{\left(d\right)}$ в область x > 1 будет распространяться не уединенный разрыв S, а разрыв с присоединенной волной Римана. Однако, такие изменения — детали, не имеющие значительного влияния на проводимые ниже рассуждения.

5. СЛУЧАЙ ${\mathit{s}}_{\mathbf{g}\mathbf{c}}^{\mathbf{\left(}\mathbf{i}\mathbf{\right)}}\mathbf{=}{\mathit{s}}_{\mathbf{g}\mathbf{c}}^{\mathbf{\left(}\mathbf{d}\mathbf{\right)}}$

При отсутствии гистерезиса фазовой проницаемости можно получить точную формулу для x_max. В этом случае во всей области, насыщенной жидкостью g, при t → ∞ имеем $s_g=s_{gc}^{\left(d\right)}$. Учитывая интеграл (19), получим, что $1=s_{gc}^{\left(d\right)}{x}_{\max }$. Следовательно

$x_{\max }=\frac{1}{s_{gc}^{\left(d\right)}}$. (20)

Согласно (20), уменьшение критической насыщенности приводит к увеличению x_max. Действительно, если жидкость g обездвиживается при меньшей насыщенности $s_{gc}^{\left(d\right)}$, то g очевидно может переместиться выше по колонне.

6. ОЦЕНКА ДЛЯ x_max

В случае общего положения $s_{gc}^{\left(i\right)}\ne s_{gc}^{\left(d\right)}$ получить простое уравнение для x_max не представляется возможным. Выведем приближенное соотношение для x_max, предположив кусочно-постоянное распределение s_g (x) при t → ∞ (рис. 4). Пренебрежем неоднородным распределением s_g в малой окрестности за разрывом S, положив $s_g=s_{gc}^{\text{'}}$ всюду при x > 1. Тогда площадь под графиком  можно представить в виде:

$s_{gc}^{\left(i\right)}+s_{gc}^{\text{'}}(x_{\max }^{*}-1)=1$, (21)

где $x_{\max }^{*}$ — оценка длины x_max, полученная в сформулированном приближении. Выражая $x_{\max }^{*}$ из уравнения (21), получим:

$x_{\max }^{*}=1+\frac{1-s_{gc}^{\text{'}}}{s_{gc}^{\left(i\right)}}$, (22)

где $s_{gc}^{\text{'}}$ — абсцисса точки D на рис. 3. Величину $s_{gc}^{\text{'}}$ можно легко найти, используя уравнение (6) и графическое построение на рис. 3. Очевидно, что $x_{\max }^{*}<x_{\max }$ (рис. 4). Таким образом, уравнение (22) дает заниженную оценку для x_max.

Для оценки точности оценки (22) сравним рассчитанные в соответствии с (22) $x_{\max }^{*}$ и значения x_max, полученные с помощью прямого численного моделирования фильтрации. Предварительный анализ результатов моделирования показал, что их удобно представить на плоскости $\{H,s_{gc}^{\left(d\right)}{x}_{\max }\}$, где

$H=\frac{s_{gc}^{\left(i\right)}-s_{gc}^{\left(d\right)}}{s_{gc}^{\left(d\right)}}$. (23)

На рис. 5, а представлены результаты расчетов при M = 3.7, n_l = n_g = 2 и различных значениях $s_{gc}^{\left(d\right)}$ и $s_{gc}^{\left(i\right)}$. Действительно, кривые соответствуют фиксированным значениям $s_{gc}^{\left(d\right)}$, а различным значениям критерия подобия H тогда соответствуют различные $s_{gc}^{\left(i\right)}$. Согласно (20), при H = 0, т.е. при $s_{gc}^{\left(i\right)}=s_{gc}^{\left(d\right)}$, имеем $s_{gc}^{\left(d\right)}{x}_{\max }=1$.

 

Рис. 5. Рассчитанные перемещения легкой жидкости g при различных критериях подобия. Точки соответствуют результатам численного моделирования, а кривые — приближенной оценке ${\mathit{x}}_{\mathbf{m}\mathbf{a}\mathbf{x}}^{\mathbf{*}}$.

 

Результаты расчетов демонстрируют, что при увеличении H перемещение x_max снижается. Действительно, большим H соответствуют большие критические насыщенности $s_{gc}^{\left(i\right)}$ и, следовательно, большие объемы фазы g, обездвиживаемой при $s_{gc}\le s_{gc}^{\left(i\right)}$. Таким образом, меньший объем g остается в подвижном состоянии, что в свою очередь снижает x_max. Этот результат сохраняется в широком диапазоне M, n_l и n_g (рис. 5, б–г).

Согласно рис. 5, $x_{\max }^{*}$ служит очень хорошей оценкой снизу для x_max. Действительно, результаты численного моделирования лежат выше соответствующих $x_{\max }^{*}$. Однако, ошибка $s_{gc}^{\left(d\right)}(x_{\max }-x_{\max }^{*})$ не превосходит 0.05 и она меньше в большинстве рассмотренных случаев. Таким образом, уравнение (22) позволяет рассчитать x_max с погрешностью не более 5%.

Проведенные расчеты демонстрируют, что  зависит от каждого из параметров $s_{gc}^{\left(d\right)}$, $s_{gc}^{\left(i\right)}$, M, n_l и n_g. Например, при увеличении M длина x_max увеличивается (рис. 5, б). Это связано с тем, что согласно (13), большим M соответствуют малые вязкости жидкости g. Можно ожидать, что менее вязкая жидкость g распространится дальше по вертикальное колонне. Малое влияние числа M на x_max связано с параметрами сильных разрывов и волн Римана, характеризующихся числом M. Так как процесс фильтрации зависит от максимальной насыщенности s_g,hy, которая в свою очередь зависит от M, то и при больших t расстояние x_max тоже зависит от M.

Несмотря на нелинейную зависимость x_max от $s_{gc}^{\left(d\right)}$, $s_{gc}^{\left(i\right)}$, M, n_l и n_g, при малых значения параметра H все рассчитанные кривые x_max (H) (и $x_{\max }^{*}\left(H\right)$) подобны друг другу (рис. 5). При H < 2, x_max сильно зависит только от H и слабо зависит от других критериев подобия. Действительно, увеличению H от 0 до 2 соответствует снижение $s_{gc}^{\left(d\right)}{x}_{\max }$ от 1 до 0.4, т.е. на 0.6 безразмерных единиц. При этом изменению других параметров (например, M) соответствуют гораздо меньшие изменения x_max. Таким образом, при H < 2 перемещение x_max фактически зависит только от критерия подобия H. При H ≥ 2 кривые x_max (H) на рис. 5 выполаживаются. При таких H влияние других критериев подобия может быть важным для точной оценки x_max. Тем не менее, один из основных результатов проведенных расчетов заключается в том, что для оценки x_max наиболее важным является критерий подобия H.

ЗАКЛЮЧЕНИЕ

Исследовано течение в насыщенной пористой среде конечного объема легкой жидкости, на которую действует сила Архимеда. В виде конечного соотношения предложена оценка максимального расстояния, на которое жидкость распространится вверх по пористой среде. Показано, что полученное соотношение для максимального расстояния дает заниженную оценку, которая тем не менее имеет погрешность не более 5%. Показано, что расстояние наиболее чувствительно к одному критерию подобия, пропорциональному разности критических насыщенностей в процессах пропитки и дренажа. Оценено влияние других критериев подобия, которое остается малым в широком диапазоне параметров течения.

Исследование выполнено за счет гранта Российского научного фонда № 19-71-10051, https://rscf.ru/project/19-71-10051/.

Об авторах

А. А. Афанасьев

МГУ им. М.В. Ломоносова

Автор, ответственный за переписку.
Email: afanasyev@imec.msu.ru
Россия, Москва

Е. А. Веденеева

МГУ им. М.В. Ломоносова

Email: el-vedeneeva@imec.msu.ru
Россия, Москва

И. Е. Михеев

МГУ им. М.В. Ломоносова

Email: afanasyev@imec.msu.ru
Россия, Москва

Список литературы

Дополнительные файлы