Прикладная математика и механика, 2020, T. 84, № 1, стр. 44-63

МАКРОСКОПИЧЕСКАЯ МОДЕЛЬ ДВУХФАЗНОГО ТЕЧЕНИЯ СЖИМАЕМЫХ ЖИДКОСТЕЙ В СРЕДЕ С ДВОЙНОЙ ПОРИСТОСТЬЮ

М. Б. Панфилов 12*, Ж. Д. Байшемиров 34, А. С. Бердышев 34

1 Институт Эли Картан – Университет Лотарингии
Нанси, Франция

2 Институт д’Аламбера – Университет Сорбонна
Париж, Франция

3 Казахский национальный педагогический университет им. Абая
Алматы, Казахстан

4 Институт информационных и вычислительных технологий
Алматы, Казахстан

* E-mail: michel.panfilov@univ-lorraine.fr

Поступила в редакцию 27.06.2019
После доработки 20.09.2019
Принята к публикации 13.11.2019

Полный текст (PDF)

Аннотация

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

Ключевые слова: двойная пористость, осреднение, двухфазное течение, сжимаемая жидкость, нелокальность, память, запаздывание, неравновесность

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

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

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

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

Авторам известны две работы, посвященные двухфазному течению сжимаемых флюидов в среде с двойной пористостью [11, 12]. В первой из них [11] обе фазы представляли собой сжимаемый газ, а во второй одна из фаз полагалась несжимаемой [12]. Сделаны математические продвижения в области доказательства сходимости исходной системы к осредненной в условиях взаимодействия нелокальности, двойной нелинейности и двойной вырождаемости (параболических уравнений), однако полученные модели остаются не полностью осредненными.

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

Численное моделирование на уровне пор, основанное на уравнениях Навье–Стокса, использует метод диффузной поверхности Кана–Хильярда [15, 17, 18], позволяющий отслеживать поверхность раздела фаз и ее деформации. Численное моделирование, основанное на законе Дарси, сфокусировано в настоящее время, прежде всего на проблеме стыковки трехмерной задачи в блоке и двумерной задачи в плоской трещине [19, 20].

Во всех этих исследованиях получены сходные результаты относительно роли капиллярной пропитки на нефтеотдачу: чем больше капиллярное давление, тем интенсивнее пропитка и тем больше нефти вымывается из блоков. Однако есть и противоречивые результаты относительно влияния темпа вытеснения на нефтеотдачу: так, например, по результатам [14, 17, 18, 20] отдача падает с ростом темпа закачки, тогда как согласно [13, 16] – растет. Это указывает на то, что общая теория двухфазных течений в среде с сильными неоднородностями остается открытым вопросом.

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

– эффект запаздывания насыщенности, вызываемый несимметричным отжимом фаз из блоков за счет расширения фаз и компакции пор;

– эффект запаздывания насыщенности, вызываемый нелинейным наложением сжимаемости и капиллярности.

Ниже метод осреднения применяется для двухмасштабной формулировки задачи в вариационной форме, по технике близкий к конструктивной части метода двухмасштабной сходимости [21]. Предлагаемый подход позволяет провести все вычисления в наиболее компактной форме.

Следует подчеркнуть, что в рамках используемого подхода сжимаемыми являются прежде всего флюиды. Пористая среда может быть тоже сжимаемой, то есть изменять объем пор в зависимости от давления жидкости их заполняющей (такая деформация среды называется компакцией), но сжимаемость среды считается слабой. Полностью деформируемые среды здесь не рассматриваются. Подходы к осреднению деформируемых упругих насыщенных пористых сред также развиваются в последнее время [22, 23].

1. Формулировка задачи.

1.1. Уравнения двухфазного течения сжимаемых жидкостей. Сжимаемость означает, что плотности жидкостей и пористость/проницаемость пористой среды зависят от давления. Традиционно сжимаемость пористых сред считается слабой и описывается линейным законом, а сжимаемость жидкостей описывается более общим экспоненциальным законом [24]:

(1.1)
$\frac{{d{{\rho }_{w}}}}{{dp}} = {{\tilde {\beta }}_{w}}{{\rho }_{w}},\quad \frac{{d{{\rho }_{o}}}}{{dp}} = {{\tilde {\beta }}_{o}}{{\rho }_{o}},\quad \frac{{d\phi }}{{dp}} = {{\tilde {\beta }}_{\phi }}{{\phi }^{0}},$
где $p$ – давление, ${{\phi }^{0}}$ – характерная постоянная величина пористости; ${{\tilde {\beta }}_{w}}$, ${{\tilde {\beta }}_{o}}$, ${{\tilde {\beta }}_{\phi }}$ – изотермические коэффициенты сжимаемости, которые считаются константами. Их размерность Pa–1. Для воды, нефти и геологических пористых сред характерная величина коэффициента сжимаемаемости составляет 10–9–10–8 Pa–1. Положительная производная $d\phi {\text{/}}dp$ означает, что пористость уменьшается, если давление жидкости в порах падает. Последнее соответствует компакции пор под весом вышележащих пород.

Для однофазных сжимаемых флюидов в деформируемой среде используется классическое упрощение уравнений течения, которое заключается в следующем. Подставим уравнения сжимаемости (1.1) в уравнение сохранения массы флюида:

${{\partial }_{t}}\left( {\phi \rho } \right){\text{ }} + \nabla \cdot \left( {\rho {\mathbf{V}}} \right) = 0,$
где V – скорость фильтрации. Тогда получим:
$0 = (\phi \rho ){\kern 1pt} '{\kern 1pt} {{\partial }_{t}}p + \rho \nabla \cdot {\mathbf{V}} + \rho {\kern 1pt} '{\kern 1pt} {\mathbf{V}} \cdot \nabla p,$
где “штрих” означает производную по давлению. Последний член, ${\mathbf{V}} \cdot \nabla p{\text{ }}\sim {{(\nabla p)}^{{\text{2}}}}$, пренебрежимо мал, так как градиенты давления в подземных пластах малы. Тогда получаем классическое уравнение однофазной фильтрации сжимаемых жидкостей в деформируемых пластах:
${{\phi }^{0}}\beta {{\partial }_{t}}p + \nabla \cdot {\mathbf{V}} = 0,$
где $\beta = \frac{{\left( {\phi \rho } \right){\text{'}}}}{{{{\phi }^{0}}\rho }} = {{\tilde {\beta }}_{w}} + {{\tilde {\beta }}_{\phi }}$, а скорость фильтрации ${\mathbf{V}}$ считается пропорциональной параметру абсолютной проницаемости среды $K$.

Пористость ${{\phi }^{0}}$ и проницаемость $K$ в этих уравнениях считаются далее постоянными.

Применим такое же допущение и для двухфазных флюидов. Рассмотрим несмешивающиеся фазы. Будем называть их водой и нефтью. Фазы имеют разные давления: ${{p}_{w}}$ и ${{p}_{o}}$. Поэтому ${{\rho }_{w}} = {{\rho }_{w}}\left( {{{p}_{w}}} \right)$ и ${{\rho }_{o}} = {{\rho }_{o}}\left( {{{p}_{o}}} \right)$. Классические уравнения теории двухфазной фильтрации применимы для шнурковой структуры течений [25], когда каждая фаза занимает свои собственные поры. Поэтому в каждый момент времени $\phi = \phi \left( {{{p}_{w}}} \right)$ для пор, занятых водой, и $\phi = \phi \left( {{{p}_{o}}} \right)$ для пор, занятых нефтью. Исходные уравнения сохранения массы фаз имеют вид:

${{\partial }_{t}}\left( {\phi {{\rho }_{w}}s} \right) + \nabla \cdot \left( {{{\rho }_{w}}{{{\mathbf{V}}}_{w}}} \right) = 0,\quad {\text{и}}\quad {{\partial }_{t}}\left( {\phi {{\rho }_{o}}\left( {1 - s} \right)} \right) + \nabla \cdot \left( {{{\rho }_{o}}{{{\mathbf{V}}}_{o}}} \right) = 0,$
где $s$ – объемная доля воды (насыщенность), $\rho $ – плотность фазы, $\phi $ – пористость среды, $p$ – давление, ${\mathbf{V}}$ – вектор скорости фильтрации (скорость Дарси). Индексы $w$ и $o$ означают воду (water) и нефть (oil). Дифференцируя по частям, получим для воды:

$\phi {{{\rho }}_{w}}{{\partial }_{t}}s + \frac{{d\left( {\phi {{{\rho }}_{w}}} \right)}}{{d{{p}_{w}}}}s{{\partial }_{t}}{{p}_{w}} + {{{\rho }}_{w}}\nabla \cdot {{{\mathbf{V}}}_{w}} + \frac{{d{{{\rho }}_{w}}}}{{d{{p}_{w}}}}\nabla {{p}_{w}} \cdot {{{\mathbf{V}}}_{w}} = 0$

Пренебрегая членами порядка $\;{{\left( {\nabla p} \right)}^{2}}$, получим окончательно систему уравнений двухфазного течения сжимаемых жидкостей в слабосжимаемой пористой среде:

(1.2)
$\begin{gathered} \phi {{\partial }_{t}}s + \phi {{\beta }_{w}}s{{\partial }_{t}}{{p}_{w}} + \nabla \cdot {{{\mathbf{V}}}_{w}} = 0 \\ - \phi {{\partial }_{t}}s + \phi {{\beta }_{o}}\left( {1 - s} \right){{\partial }_{t}}{{p}_{o}} + \nabla \cdot {{{\mathbf{V}}}_{o}} = 0, \\ \end{gathered} $
где ${{\beta }_{w}} \equiv {{\tilde {\beta }}_{w}} + {{\tilde {\beta }}_{\phi }}$, ${{\beta }_{o}} \equiv {{\tilde {\beta }}_{o}} + {{\tilde {\beta }}_{\phi }}$, а величина ${{\phi }^{0}}$ обозначена просто как $\phi $.

Система (1.2) дополняется уравнениями сохранения импульса в форме закона Дарси для каждой фазы и уравнением капиллярного равновесия, связывающим давления в фазах:

(1.3)
${{{\mathbf{V}}}_{\alpha }} = - K{{\lambda }_{\alpha }}\nabla {{p}_{\alpha }},\quad \alpha = w,o;\quad {{\lambda }_{\alpha }}\left( {x,s} \right) \equiv \frac{{{{k}_{\alpha }}\left( {x,s} \right)}}{{{{\mu }_{\alpha }}}}$
(1.4)
${{p}_{o}} = {{p}_{w}} + {{p}_{c}}\left( s \right),$
где $K$ – абсолютная проницаемость среды, $\mu $ – динамическая вязкость фазы. Подчеркнем, что $K$ и $\phi $ считаются независящими от давления, что соответствует слабосжимаемой пористой среде.

Функции относительных фазовых проницаемостей ${{k}_{w}}\left( s \right)$, ${{k}_{o}}\left( s \right)$ и капиллярного давления ${{p}_{c}}\left( s \right)$ заданы. Они обладают следующими свойствами:

${{k}_{w}}(s)$ непрерывная и монотонно неубывающая функция водонасыщенности $s$, такая что ${{k}_{w}} \equiv 0$ при $s \in [0,{{s}_{ * }}]$, и ${{k}_{w}}{\text{(1)}} = {\text{1}}$;

${{k}_{o}}(s)$ непрерывная и монотонно невозрастающая функция такая, что ${{k}_{o}}\left( 0 \right) = {\text{1}}$ и ${{k}_{o}} \equiv 0$ при $s \in [s*,{\text{1}}]$;

${{p}_{c}}(s)$ непрерывная и монотонно невозрастающая функция такая, что ${{p}_{c}}\left( {\text{1}} \right) = 0$ и ${{p}_{c}} \to \infty $ при $s \to {{s}_{ * }}$. Функция ${{p}_{c}}$ не определена при $s \in [0,{{s}_{ * }}]$.

1.2. Задача осреднения. Пусть резервуар $\Omega \subset {{\Re }^{d}}$ ($d = 2,3$) – ограниченная, связная область с периодической микроструктурой с характерным масштабом ε, который является отношением периода неоднородности к длине всей области. Параметр $\varepsilon $ малый: $0 < \varepsilon \ll 1$ (рис. 1).

Рис. 1.

Структура среды.

Единичная ячейка этой микроструктуры выделена черным квадратом. Будучи растянутой в ${{{\varepsilon }}^{{ - 1}}}$ раз по каждой оси она образует куб $Y = {{\left( {0,1} \right)}^{d}}$ со стороной 1, где $d$ – размерность пространства. Этот куб состоит из двух подобластей: связной подобласти ${{Y}^{\mathcal{F}}}$ (“трещины”) и ${{Y}^{\mathcal{M}}}$ (“блок матрицы”). Обозначим через $\Gamma $ границу раздела между ${{Y}^{\mathcal{F}}}$ и ${{Y}^{\mathcal{M}}}$ в $Y$.

Четыре уравнения сохранения массы и импульса (1.2) и (1.3) можно свести к двум уравнениям, исключив скорости фильтрации:

(1.5)
$\phi {{\partial }_{t}}{{s}_{\alpha }} + \phi {{s}_{\alpha }}{{{\beta }}_{\alpha }}{{\partial }_{t}}{{p}_{\alpha }} = \nabla \cdot \left( {K{{{\lambda }}_{\alpha }}\nabla {{p}_{\alpha }}} \right),\quad \alpha = w,o,$
где $s \equiv {{s}_{w}}$, ${{s}_{o}} = 1 - {{s}_{w}}$.

Система (1.5) и соотношение (1.4) определяют три неизвестных функции: насыщенность воды $s$, давление в воде ${{p}_{w}}$ и давление в нефти ${{p}_{o}}$. На границе раздела сред $\Gamma $ нормальные потоки фаз и фазовые давления непрерывны.

Граничными условиями являются заданное давление воды и насыщенность. Начальные условия фиксируют заданное распределение насыщенности и давления в области:

$s\left( {x,0} \right) = s_{{}}^{0}\left( x \right),\quad {{p}_{w}}\left( {x,0} \right) = p_{w}^{0}\left( x \right)$

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

$\phi \left( x \right) = \left\{ \begin{gathered} {{\phi }^{\mathcal{F}}} \hfill \\ {{\phi }^{\mathcal{M}}} \hfill \\ \end{gathered} \right.,\quad {{\beta }_{\alpha }}\left( x \right) = \left\{ \begin{gathered} \beta _{\alpha }^{\mathcal{F}},\quad x \in {{\Omega }^{\mathcal{F}}} \hfill \\ \beta _{\alpha }^{\mathcal{M}},\quad x \in {{\Omega }^{\mathcal{M}}} \hfill \\ \end{gathered} \right.$
${{\lambda }_{\alpha }}\left( {s,x} \right) = \left\{ \begin{gathered} \lambda _{\alpha }^{\mathcal{F}}\left( s \right) \hfill \\ \lambda _{\alpha }^{\mathcal{M}}\left( s \right) \hfill \\ \end{gathered} \right.,\quad {{p}_{c}}\left( {s,x} \right) = \left\{ \begin{gathered} p_{c}^{\mathcal{F}}\left( s \right),\quad x \in {{\Omega }^{\mathcal{F}}} \hfill \\ p_{c}^{\mathcal{M}}\left( s \right),\quad x \in {{\Omega }^{\mathcal{M}}} \hfill \\ \end{gathered} \right.,\quad \alpha = w,o,$
где ${{\phi }^{\mathcal{F}}}$, ${{\phi }^{\mathcal{M}}}$, $\beta _{\alpha }^{\mathcal{F}}$, $\beta _{\alpha }^{\mathcal{M}}$, $\lambda _{\alpha }^{\mathcal{F}}\left( s \right)$, $\lambda _{\alpha }^{\mathcal{M}}\left( s \right)$, $p_{с}^{\mathcal{F}}\left( s \right)$, $p_{с}^{\mathcal{M}}\left( s \right)$ не зависят от $\varepsilon $. Параметры ${{\phi }^{\mathcal{F}}}$, ${{\phi }^{\mathcal{M}}}$, $\beta _{\alpha }^{\mathcal{F}}$, $\beta _{\alpha }^{\mathcal{M}}$ положительны. Функции $\lambda _{\alpha }^{\mathcal{F}}\left( s \right)$, $\lambda _{\alpha }^{\mathcal{M}}\left( s \right)$, $p_{с}^{\mathcal{F}}\left( s \right)$, $p_{с}^{\mathcal{M}}\left( s \right)$ неотрицательны. Полагаем, что $p_{c}^{\mathcal{M}}\left( s \right) \geqslant p_{c}^{\mathcal{F}}\left( s \right)$, так как капиллярное давление больше там, где проницаемость меньше.

В то же время, абсолютная проницаемость блоков и трещин может различаться очень значительно. Известно, что при вариации пористости в пределах 0.1–0.2 проницаемость меняется на несколько порядков. Поэтому принимаем основное условие двойной пористости:

(1.6)
$K\left( x \right) = \left\{ \begin{gathered} {{K}^{\mathcal{F}}},\quad x \in {{\Omega }^{\mathcal{F}}} \hfill \\ \delta {{K}^{\mathcal{M}}},\quad x \in {{\Omega }^{\mathcal{M}}} \hfill \\ \end{gathered} \right.$
где параметр $\delta $ мал, коэффициенты ${{K}^{\mathcal{M}}}$, ${{K}^{\mathcal{F}}}$ положительные и не зависят от $\varepsilon $.

Параметр $\delta > 0$ описывает степень контрастности между проницаемостями блоков и трещин. Большой контраст вызывает запаздывание или память. Степень запаздывания (или “длина памяти”) выражается другим параметром, $\omega $, который определяется следующим образом:

$\omega = \frac{{{{\varepsilon }^{2}}}}{\delta }$

Параметр $\omega $ входит в задачу в явном виде, если исходные уравнения (1.5) записаны не через макроскопические координаты ${{x}_{i}}$, а микроскопические (растянутые) координаты ${{y}_{i}} = {{x}_{i}}{\text{/}}\varepsilon $. Тогда в уравнениях (1.5) правая часть примет вид на блоке: $ \ldots = {{\omega }^{{ - 1}}}{{\nabla }_{y}} \cdot ({{K}^{\mathcal{M}}}{{{\lambda }}_{\alpha }}{{\nabla }_{y}}{{p}_{\alpha }})$.

По степени контрастности и длине памяти интерес представляют два типа сред:

– среды с долгой памятью: $\omega \sim 1$, или с сильным контрастом: $\delta \sim {{\varepsilon }^{2}}$ (так называемые ${{\varepsilon }^{2}}$-среды);

– среды с короткой памятью: $\varepsilon \ll \omega \ll 1$, или с умеренным контрастом: ${{\varepsilon }^{2}} \ll \delta \ll \varepsilon $.

Классификация сред по параметрам $\delta $ и $\omega $ представлена на рис. 2.

Рис. 2.

Классификация сред по параметрам контрастности свойств ($\delta $) и длины возникающей памяти ($\omega $).

Существуют также среды неконтрастные: $\delta \sim 1$, слабоконтрастные $\varepsilon \ll \delta \ll 1$ и среды с непроницаемыми блоками: $\delta \ll {{\varepsilon }^{2}}$. В первых и вторых обменный процесс между блоками и трещинами происходит без запаздывания, а в третьих полностью отсутствует из-за очень сильного запаздывания, поэтому ни те, ни другие, ни третьи для целей настоящего исследования интереса не представляют.

2. Метод осреднения с расщеплением нелокальности и нелинейности.

2.1. Расщепление нелокальности и нелинейности. Целью данной статьи является построение полностью осредненной модели. Полностью осредненной называется модель, макроскопические уравнения которой не содержат микроскопические переменные, а задача на ячейке, определяющая макроскопические коэффициенты, не зависит от макроскопических переменных и поэтому решается только один раз [4].

Чтобы обойти наложение нелинейности и нелокальности, описанное ранее, мы разнесем их на разные уровни асимптотического разложения по параметру длины памяти $\omega $, то есть нелинейность сохраняется в нулевом приближении, а нелокальность – в первом. Поскольку в первом приближении асимптотические разложения обычно приводят к линейным задачам, то нелокальность возможно будет описать явно и полностью. Ранее этот метод был применен для несжимаемых жидкостей [8].

Таким образом, параметр $\omega $ должен быть малым, то есть мы приходим к средам с умеренным контрастом между проницаемостями блоков и трещин. В общем случае мы вынуждены работать с двухпараметрическими разложениями: по периоду неоднородности $\varepsilon $ и по параметру длины памяти $\omega $, которые независимы. Учитывая упорядочение сред по параметрам $\delta $ и $\omega $, представленное на рис. 2, асимптотическое разложение любой функции $f$ будет выглядеть следующим образом:

(2.1)
${{f}^{{\varepsilon ,\omega }}}\left( {x,t} \right) = \sum\limits_{k = 0}^\infty {\sum\limits_{j = 0}^\infty {{{\omega }^{k}}{{\varepsilon }^{j}}} } f_{{kj}}^{\varepsilon }\left( {x,t} \right) = {{f}_{{00}}} + \omega f_{{01}}^{\varepsilon } + \varepsilon f_{{10}}^{\varepsilon } + \ldots ,$
где верхние индексы $\varepsilon $ и $\omega $ означают зависимость от соответствующих параметров. При этом учтено, что задача регулярна по $\omega $ (за исключением малых времен $t\sim \omega $ после возмущения) и сингулярно возмущена по $\varepsilon $, поэтому коэффициенты $f_{{kj}}^{\varepsilon }$ в общем случае зависят от $\varepsilon $ и не зависят от $\omega $.

Если $\varepsilon \to 0$, то неоднородности стягиваются в точку, то есть среда превращается в однородную осредненную. Если $\omega \to 0$, то в системе исчезает память. Поэтому нулевое приближение ${{f}_{{00}}}$ соответствует осредненному значению функции в системе без контраста свойств блоков и трещин.

Для построения осредненной модели достаточно нулевого и первого приближения: ${{f}_{{00}}}$, $f_{{01}}^{\varepsilon }$, $f_{{10}}^{\varepsilon }$. При этом член $\omega f_{{01}}^{\varepsilon }$ описывает вклад контрастности свойств, то есть эффекты памяти, а $\varepsilon f_{{10}}^{\varepsilon }$ ответствен за вклад пространственных осцилляций проницаемости.

Очевидно, что техника осреднения значительно упрощается, если выбрать $\omega $ так, чтобы например:

(2.2)
$\omega = \sqrt \varepsilon ,\quad {\text{или}}\quad \delta = {{\varepsilon }^{{3/2}}},$
что использовалось в [8]. Тогда полное асимптотическое разложение сводится к однопараметрическому разложению по $\sqrt \varepsilon $.

2.2. Двухмасштабная формулировка задачи. Прежде всего, введем расширение решения $\tilde {s}\left( {x,y,t} \right)$, $\tilde {p}\left( {x,y,t} \right)$, так что

${{p}_{w}}\left( {x,t} \right) = {{\left. {\tilde {p}\left( {x,y,t} \right)} \right|}_{{y = x{\text{/}}\varepsilon }}},\quad s\left( {x,t} \right) = {{\left. {\tilde {s}\left( {x,y,t} \right)} \right|}_{{y = x{\text{/}}\varepsilon }}}$

Тогда следующее справедливо для производных:

$\frac{{\partial p\left( {x,t} \right)}}{{\partial {{x}_{i}}}} \to {{\left. {\left( {\frac{{\partial{ \tilde {p}}\left( {x,y,t} \right)}}{{\partial {{x}_{i}}}} + \frac{1}{\varepsilon }\frac{{\partial{ \tilde {p}}\left( {x,y,t} \right)}}{{\partial {{y}_{i}}}}} \right)} \right|}_{{y = x{\text{/}}\varepsilon }}}$

Поскольку операции расширения и дифференцирования коммутируют, получаем, что аргумент y может рассматриваться как независимая переменная от $x$, а в конечных результатах надо положить его равным $x{\text{/}}\varepsilon $.

Двухмасштабная формулировка уравнений (1.5) имеет следующий вид:

(2.3)
$\begin{gathered} \phi \left( y \right){{\partial }_{t}}{{s}_{\alpha }} + \phi \left( y \right){{{\beta }}_{\alpha }}\left( y \right){{s}_{\alpha }}{{\partial }_{t}}{{p}_{\alpha }} = \left( {{{\partial }_{{xi}}} + \frac{1}{\varepsilon }{{\partial }_{{yi}}}} \right)\left( {K\left( y \right){{{\lambda }}_{\alpha }}\left( {{{\partial }_{{xi}}}{{p}_{\alpha }} + \frac{1}{\varepsilon }{{\partial }_{{yi}}}{{p}_{\alpha }}} \right)} \right) \\ x \in \Omega ,\quad y \in Y,\quad t \in \left( {0,T} \right) \\ {{p}_{o}} = {{p}_{w}} + {{p}_{с}},\quad {{s}_{o}} = 1 - {{s}_{w}} \\ \end{gathered} $

Здесь и далее по повторяющимся индексам $i$, $k$, $j$ идет суммирование от 1 до 3. С условиями:

(2.4)
$\begin{gathered} {{\left. {{{s}_{w}}} \right|}_{{t = 0}}} = {{s}^{0}}(x),\quad {{\left. {{{p}_{w}}} \right|}_{{t = 0}}} = p_{w}^{0}(x) \\ {{\left[ {K{{\lambda }_{\alpha }}\frac{{\partial {{p}_{\alpha }}}}{{\partial n}}} \right]}_{\Gamma }} = 0,\quad {{\left[ {{{p}_{\alpha }}} \right]}_{\Gamma }} = 0,\quad \alpha = w,o \\ {{s}_{w}},{{p}_{w}}\; - \;{\text{1 - периодичны по }}y \\ \end{gathered} $
и условиями Неймана или Дирихле на границе области $\partial \Omega $. Здесь [⋅] – обозначают скачок функции.

Задачу (2.3), (2.4) можно представить в вариационной формулировке:

(2.5)
$\begin{gathered} \int\limits_\Omega {\int\limits_Y {\phi w{{\partial }_{t}}{{s}_{\alpha }}dxdy} } + \int\limits_\Omega {\int\limits_Y {\phi {{{\beta }}_{\alpha }}{{s}_{\alpha }}w{{\partial }_{t}}{{p}_{\alpha }}dxdy} } = \\ = - \frac{1}{{{{\varepsilon }^{2}}}}\int\limits_\Omega {\int\limits_Y {K{{{\lambda }}_{\alpha }}\left( {{{\partial }_{{yi}}}{{p}_{\alpha }} + \varepsilon {{\partial }_{{xi}}}{{p}_{\alpha }}} \right)\left( {{{\partial }_{{yi}}}w + \varepsilon {{\partial }_{{xi}}}w} \right)dy} dx} = \\ = - \frac{1}{{{{\varepsilon }^{2}}}}\int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{F}}}} {{{K}^{\mathcal{F}}}{\lambda }_{\alpha }^{\mathcal{F}}\left( {{{\partial }_{{yi}}}{{p}_{\alpha }} + \varepsilon {{\partial }_{{xi}}}{{p}_{\alpha }}} \right)\left( {{{\partial }_{{yi}}}w + \varepsilon {{\partial }_{{xi}}}w} \right)dy} dx} - \\ - \;\frac{1}{{\sqrt \varepsilon }}\int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{M}}}} {{{K}^{\mathcal{M}}}{\lambda }_{\alpha }^{\mathcal{M}}\left( {{{\partial }_{{yi}}}{{p}_{\alpha }} + \varepsilon {{\partial }_{{xi}}}{{p}_{\alpha }}} \right)\left( {{{\partial }_{{yi}}}w + \varepsilon {{\partial }_{{xi}}}w} \right)dy} dx} \\ \end{gathered} $
для любой функции $w = w\left( {x,y} \right)$, действующей в пространстве $\Omega \times Y$, непрерывной и периодичной по $y$, а также такой, что ${{\left. w \right|}_{{\partial \Omega }}} = 0$.

Соотношение (2.5) получается умножением уравнения (2.3) на функцию $w$, с последующим интегрированием по частям, и использованием теоремы Гаусса–Остроградского. Возникающие при этом интегралы по границе области $\partial \Omega $ обнуляются из-за граничного условия на функцию $w$. Интегралы по границе периода $\partial Y$, в свою очередь, обращаются в нуль из-за периодичности всех функций по $y$.

2.3. Асимптотическое разложение. Осреднение строится методом двухмасштабных асимптотических разложений по параметру $\varepsilon $, а также регулярных асимптотических рядов по параметру нелокальности $\omega $, что, с учетом параметризации (2.2), означает появление членов, содержащих дробные степени $\sqrt \varepsilon $. Общая структура разложений имеет следующий вид, для $\alpha = w,o$:

${{p}_{\alpha }}\left( {x,y,t} \right) = \left\{ \begin{gathered} {{p}_{{\alpha 0}}}\left( {x,t} \right) + \sqrt \varepsilon p_{{\alpha ,1/2}}^{\mathcal{M}}\left( {x,y,t} \right) + \varepsilon p_{{\alpha 1}}^{\mathcal{M}}\left( {x,y,t} \right) + \ldots ,\quad y \in {{Y}^{\mathcal{M}}} \hfill \\ {{p}_{{\alpha 0}}}\left( {x,t} \right) + \sqrt \varepsilon p_{{\alpha ,1/2}}^{\mathcal{F}}\left( {x,t} \right) + \varepsilon p_{{\alpha 1}}^{\mathcal{F}}\left( {x,y,t} \right) + \ldots ,\quad y \in {{Y}^{\mathcal{F}}} \hfill \\ \end{gathered} \right.$
${{s}_{\alpha }}\left( {x,y,t} \right) = \left\{ \begin{gathered} {{s}_{{\alpha 0}}}\left( {x,t} \right) + \sqrt \varepsilon s_{{\alpha ,1/2}}^{\mathcal{M}}\left( {x,y,t} \right) + \varepsilon s_{{\alpha 1}}^{\mathcal{M}}\left( {x,y,t} \right) + \ldots ,\quad y \in {{Y}^{\mathcal{M}}} \hfill \\ {{s}_{{\alpha 0}}}\left( {x,t} \right) + \sqrt \varepsilon s_{{\alpha ,1/2}}^{\mathcal{F}}\left( {x,t} \right) + \varepsilon s_{{\alpha 1}}^{\mathcal{F}}\left( {x,y,t} \right) + \ldots ,\quad y \in {{Y}^{\mathcal{F}}} \hfill \\ \end{gathered} \right.$

Причем выполняются условия непрерывности фазовых давлений:

(2.6)
${{\left. {p_{{\alpha ,1/2}}^{\mathcal{M}}} \right|}_{{y \in \Gamma }}} = p_{{\alpha ,1/2}}^{\mathcal{F}}\left( {x,t} \right)$
(Независимость первых членов разложений от $y$ легко доказывается подстановкой разложения в (2.5)).

Тогда для нелинейных функций, входящих в (2.5), получим:

(2.7)
${{\lambda }_{\alpha }}\left( {x,y,t} \right) = \left\{ \begin{gathered} {{\lambda }_{{\alpha 0}}}\left( {x,t} \right) + \sqrt \varepsilon \lambda _{{\alpha ,1/2}}^{\mathcal{M}}\left( {x,y,t} \right) + \varepsilon \ldots ,\quad y \in {{Y}^{\mathcal{M}}} \hfill \\ {{\lambda }_{{\alpha 0}}}\left( {x,t} \right) + \sqrt \varepsilon \lambda _{{\alpha ,1/2}}^{\mathcal{F}}\left( {x,t} \right) + \varepsilon \ldots ,\quad y \in {{Y}^{\mathcal{F}}} \hfill \\ \end{gathered} \right.$

Интегральное тождество (2.5) примет вид:

(2.8)
$\begin{gathered} \int\limits_\Omega {\int\limits_Y {\phi w({{\partial }_{t}}{{s}_{{\alpha 0}}} + \sqrt \varepsilon {{\partial }_{t}}{{s}_{{\alpha 1/2}}} + \ldots )dxdy} } + \\ + \;\int\limits_\Omega {\int\limits_Y {\phi {{{\beta }}_{\alpha }}w({{s}_{{\alpha 0}}} + \sqrt \varepsilon {{s}_{{\alpha 1/2}}} + \ldots )({{\partial }_{t}}{{p}_{{\alpha 0}}} + \sqrt \varepsilon {{\partial }_{t}}{{p}_{{\alpha 1/2}}} + \ldots )dxdy} } = \\ = - \frac{1}{\varepsilon }\int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{F}}}} {{{K}^{\mathcal{F}}}({\lambda }_{{\alpha 0}}^{\mathcal{F}} + \sqrt \varepsilon {\lambda }_{{\alpha 1/2}}^{\mathcal{F}} + \ldots )({{\partial }_{{yi}}}{{p}_{{\alpha 1}}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}} + \sqrt \varepsilon ({{\partial }_{{yi}}}{{p}_{{\alpha 3/2}}} + {{\partial }_{{xi}}}{{p}_{{\alpha 1/2}}})) \times } } \\ \times \;({{\partial }_{{yi}}}w + \varepsilon {{\partial }_{{xi}}}w)dydx - \\ - \;\int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{M}}}} {{{K}^{\mathcal{M}}}{\lambda }_{{\alpha 0}}^{\mathcal{M}}({{\partial }_{{yi}}}{{p}_{{\alpha 1/2}}} + \sqrt \varepsilon ({{\partial }_{{yi}}}{{p}_{{\alpha 1}}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}}))({{\partial }_{{yi}}}w + \varepsilon {{\partial }_{{xi}}}w)dy} dx + \mathcal{O}\left( \varepsilon \right)} \\ \end{gathered} $

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

2.4. Результат осреднения – макроскопическая модель. Приведем сразу результат осреднения (полный вывод дается в разд. 3).

Определим осредненные фазовые давления и насыщенности в блоках и трещинах как:

(2.9)
$\begin{gathered} \mathcal{P}_{\alpha }^{\mathcal{F}} \equiv \frac{1}{{\left| {{{Y}^{\mathcal{F}}}} \right|}}\int\limits_{{{Y}^{\mathcal{F}}}} {({{p}_{{\alpha 0}}} + \sqrt \varepsilon p_{{\alpha 1/2}}^{\mathcal{F}})dy = {{p}_{{\alpha 0}}} + \sqrt \varepsilon p_{{\alpha 1/2}}^{\mathcal{F}}} \\ \mathcal{P}_{\alpha }^{\mathcal{M}} \equiv \frac{1}{{\left| {{{Y}^{\mathcal{M}}}} \right|}}\int\limits_{{{Y}^{\mathcal{M}}}} {({{p}_{{\alpha 0}}} + \sqrt \varepsilon p_{{\alpha 1/2}}^{\mathcal{M}})dy} = {{p}_{{\alpha 0}}} + \sqrt \varepsilon {{\left\langle {p_{{\alpha 1/2}}^{\mathcal{M}}} \right\rangle }_{\mathcal{M}}} \\ \end{gathered} $
(2.10)
$\begin{gathered} \mathcal{S}_{\alpha }^{\mathcal{F}} \equiv \frac{1}{{\left| {{{Y}^{\mathcal{F}}}} \right|}}\int\limits_{{{Y}^{\mathcal{F}}}} {({{s}_{{\alpha 0}}} + \sqrt \varepsilon s_{{\alpha 1/2}}^{\mathcal{F}})dy = {{s}_{{\alpha 0}}} + \sqrt \varepsilon s_{{\alpha 1/2}}^{\mathcal{F}}} \\ \mathcal{S}_{\alpha }^{\mathcal{M}} \equiv \frac{1}{{\left| {{{Y}^{\mathcal{M}}}} \right|}}\int\limits_{{{Y}^{\mathcal{M}}}} {({{s}_{{\alpha 0}}} + \sqrt \varepsilon s_{{\alpha 1/2}}^{\mathcal{M}})y} = {{s}_{{\alpha 0}}} + \sqrt \varepsilon {{\left\langle {s_{{\alpha 1/2}}^{\mathcal{M}}} \right\rangle }_{\mathcal{M}}}, \\ \end{gathered} $
где ${{\left\langle \cdot \right\rangle }_{\alpha }} \equiv \frac{1}{{\left| {{{Y}^{\alpha }}} \right|}}\int_{{{Y}^{\alpha }}}^{} {\left( \cdot \right)dy} $, $\alpha = \mathcal{M},\mathcal{F}$.

Между ними существуют следующие связи:

(2.11)
$\mathcal{S}_{o}^{} = 1 - \mathcal{S}_{w}^{\mathcal{F}},\quad \mathcal{S}_{o}^{\mathcal{M}} = 1 - \mathcal{S}_{w}^{\mathcal{M}},\quad \mathcal{P}_{o}^{\mathcal{F}} = \mathcal{P}_{w}^{\mathcal{F}} + p_{c}^{\mathcal{F}}\left( {\mathcal{S}_{w}^{\mathcal{F}}} \right),\quad \mathcal{P}_{o}^{\mathcal{M}} = \mathcal{P}_{w}^{\mathcal{M}} + p_{c}^{\mathcal{M}}\left( {\mathcal{S}_{w}^{\mathcal{M}}} \right)$

Заметим, что макроскопические капиллярные давления в трещинах и блоках, $\mathcal{P}_{c}^{\mathcal{F}}\left( {\mathcal{S}_{w}^{\mathcal{F}}} \right)$ и $\mathcal{P}_{c}^{\mathcal{M}}\left( {\mathcal{S}_{w}^{\mathcal{M}}} \right)$, не фигурируют в этой модели, так как они просто равны исходным функциям капиллярного давления, взятым от средних насыщенностей, с точностью $\mathcal{O}\left( \varepsilon \right)$:

(2.12)
$\mathcal{P}_{c}^{\mathcal{F}}\left( {\mathcal{S}_{w}^{\mathcal{F}}} \right) = p_{c}^{\mathcal{F}}\left( {\mathcal{S}_{w}^{\mathcal{F}}} \right),\quad \mathcal{P}_{c}^{\mathcal{M}}\left( {\mathcal{S}_{w}^{\mathcal{M}}} \right) = p_{c}^{\mathcal{M}}\left( {\mathcal{S}_{w}^{\mathcal{M}}} \right)$

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

(2.13)
$\begin{gathered} {{\phi }^{\mathcal{F}}}\left( {1 - \theta } \right){{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{F}} + {{\phi }^{\mathcal{F}}}{\beta }_{\alpha }^{\mathcal{F}}\left( {1 - \theta } \right)\mathcal{S}_{\alpha }^{\mathcal{F}}{{\partial }_{t}}\mathcal{P}_{\alpha }^{\mathcal{F}} = \\ = {{\partial }_{{xi}}}({{\mathbb{K}}_{{ik}}}\lambda _{\alpha }^{\mathcal{F}}(\mathcal{S}_{\alpha }^{\mathcal{F}}){{\partial }_{{xk}}}\mathcal{P}_{\alpha }^{\mathcal{F}}) + {{\xi }_{\alpha }}(\mathcal{P}_{\alpha }^{\mathcal{M}} - \mathcal{P}_{\alpha }^{\mathcal{F}}),\quad \alpha = w,o \\ \end{gathered} $
(2.14)
$\mathcal{P}_{w}^{\mathcal{M}} = \mathcal{P}_{w}^{\mathcal{F}} - \frac{{\tau _{w}^{{{\text{comp}}}}}}{{\beta _{w}^{\mathcal{M}}\mathcal{S}_{w}^{\mathcal{M}}}}[{{\partial }_{t}}\mathcal{S}_{w}^{\mathcal{M}} + \beta _{w}^{\mathcal{M}}\mathcal{S}_{w}^{\mathcal{M}}{{\partial }_{t}}\mathcal{P}_{w}^{\mathcal{M}}]$
(2.15)
$p_{с}^{\mathcal{M}}(\mathcal{S}_{w}^{\mathcal{M}}) = p_{с}^{\mathcal{F}}(\mathcal{S}_{w}^{\mathcal{F}}) + \frac{{(\tau _{{}}^{{{\text{cap}}}} + \tau _{{}}^{{{\text{cc}}}})}}{{{\beta }_{o}^{\mathcal{M}}}}{{\partial }_{t}}\mathcal{S}_{w}^{\mathcal{M}} + (\tau _{w}^{{{\text{comp}}}} - \tau _{o}^{{{\text{comp}}}}){{\partial }_{t}}\mathcal{P}_{w}^{\mathcal{M}},$
где $\theta \equiv \left| {{{Y}^{\mathcal{M}}}} \right|$ – объемная доля блока.

Эффективная проницаемость определена как:

(2.16)
${{\mathbb{K}}_{{ik}}} \equiv \int\limits_{{{Y}^{\mathcal{F}}}} {{{K}^{\mathcal{F}}}\left( {\frac{{\partial {{\psi }_{k}}}}{{\partial {{y}_{i}}}} + {{\delta }_{{ik}}}} \right)dy} ,$
а ячеечные функции ${{\psi }_{k}}(y)$ – решения первой задачи на ячейке (на трещине), для $k = 1,2,3$:

(2.17)
$\begin{gathered} \frac{\partial }{{\partial {{y}_{i}}}}\left( {{{K}^{\mathcal{F}}}\left( {\frac{{\partial {{\psi }_{k}}}}{{\partial {{y}_{i}}}} + {{\delta }_{{ik}}}} \right)} \right) = 0,\quad y \in {{Y}^{\mathcal{F}}} \\ {{\left. {{{K}^{\mathcal{F}}}\left( {\frac{{\partial {{\psi }_{k}}}}{{\partial {{y}_{i}}}} + {{\delta }_{{ik}}}} \right)n_{i}^{\Gamma }} \right|}_{{y \in \Gamma }}} = 0 \\ {{\psi }_{k}} - y{\text{ - периодична}} \\ {{\left\langle {{{\psi }_{k}}} \right\rangle }_{{{{Y}^{\mathcal{F}}}}}} = 0 \\ \end{gathered} $

Времена запаздывания между блоками и трещинами, вызванные сжимаемостью фаз, и времена запаздывания, вызванные капиллярностью и совместным действием капиллярности и сжимаемости, определены как:

(2.18)
$\tau _{\alpha }^{{{\text{comp}}}}(\mathcal{S}_{w}^{\mathcal{M}}) = \sqrt \varepsilon {{\left\langle \varphi \right\rangle }_{\mathcal{M}}}{{\phi }^{\mathcal{M}}}\frac{{{\beta }_{\alpha }^{\mathcal{M}}\mathcal{S}_{\alpha }^{\mathcal{M}}}}{{\lambda _{\alpha }^{\mathcal{M}}}}$
(2.19)
$\tau _{{}}^{{{\text{cap}}}}(\mathcal{S}_{w}^{\mathcal{M}}) = \sqrt \varepsilon {{\left\langle \varphi \right\rangle }_{\mathcal{M}}}{{\phi }^{\mathcal{M}}}{\beta }_{o}^{\mathcal{M}}\frac{{(\lambda _{w}^{\mathcal{M}} + \lambda _{o}^{\mathcal{M}})}}{{\lambda _{w}^{\mathcal{M}}\lambda _{o}^{\mathcal{M}}}}$
(2.20)
$\tau _{{}}^{{{\text{cс}}}}(\mathcal{S}_{w}^{\mathcal{M}}) = - {\beta }_{o}^{\mathcal{M}}\frac{{dp_{c}^{\mathcal{M}}}}{{d\mathcal{S}_{w}^{\mathcal{M}}}}\tau _{o}^{{{\text{comp}}}},$
а функция ${\varphi }(y)$ – решение второй ячеечной задачи (на блоке):

(2.21)
$\frac{\partial }{{\partial {{y}_{i}}}}\left( {{{K}^{\mathcal{M}}}\frac{{\partial \varphi }}{{\partial {{y}_{i}}}}} \right) = - 1,\quad y \in {{Y}^{\mathcal{M}}},\quad {{\left. \varphi \right|}_{{y \in \Gamma }}} = 0$

Параметры массообмена ${{\xi }_{\alpha }}$ определены следующим образом:

${{\xi }_{\alpha }} = \frac{{\theta \lambda _{\alpha }^{\mathcal{M}}}}{{\sqrt \varepsilon {{{\left\langle \varphi \right\rangle }}_{\mathcal{M}}}}}$

Модель (2.13)–(2.15) является формальным асимптотическим разложением исходной задачи, сохраняющим члены порядка $\mathcal{O}(\sqrt \varepsilon )$. Порядок отброшенных членов, таким образом, $\mathcal{O}\left( \varepsilon \right)$.

3. Вывод осредненной модели

3.1. Первый шаг осреднения – разложение на трещине. Пусть $w = w\left( {x,y} \right)$ в $\Omega \times Y$, ${{\left. w \right|}_{{\partial \Omega }}} = 0$, $w$ непрерывна и периодична по $y$. Тогда из интегрального тождества (2.8) для отрицательных степеней $\varepsilon $ получим:

(3.1)
$\begin{gathered} 0 = \frac{1}{\varepsilon }\int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{F}}}} {{{K}^{\mathcal{F}}}{\lambda }_{{\alpha 0}}^{{}}({{\partial }_{{yi}}}p_{{\alpha 1}}^{\mathcal{F}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}}){{\partial }_{{yi}}}wdydx} } \\ 0 = \frac{1}{{\sqrt \varepsilon }}\int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{F}}}} {{{K}^{\mathcal{F}}}[{\lambda }_{{\alpha 0}}^{{}}({{\partial }_{{yi}}}p_{{\alpha 3/2}}^{\mathcal{F}} + {{\partial }_{{xi}}}p_{{\alpha 1/2}}^{\mathcal{F}}) + {\lambda }_{{\alpha 1/2}}^{\mathcal{F}}({{\partial }_{{yi}}}p_{{\alpha 1}}^{\mathcal{F}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}})]{{\partial }_{{yi}}}wdydx} } \\ \end{gathered} $

Из системы (3.1) вытекает:

$0 = - \int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{F}}}} {w{{\partial }_{{yi}}}({{K}^{\mathcal{F}}}{\lambda }_{{\alpha 0}}^{{}}({{\partial }_{{yi}}}p_{{\alpha 1}}^{\mathcal{F}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}}))dxdy} } + \int\limits_\Omega {\int\limits_\Gamma {w{{K}^{\mathcal{F}}}{\lambda }_{{\alpha 0}}^{{}}({{\partial }_{{yi}}}p_{{\alpha 1}}^{\mathcal{F}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}})n_{i}^{\Gamma }dxdy} } $
$0 = \int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{F}}}} {w{{\partial }_{{yi}}}({{K}^{\mathcal{F}}}{\lambda }_{{\alpha 0}}^{{}}({{\partial }_{{yi}}}p_{{\alpha 3/2}}^{\mathcal{F}} + {{\partial }_{{xi}}}p_{{\alpha 1/2}}^{\mathcal{F}}) + } } {{K}^{\mathcal{F}}}{\lambda }_{{\alpha 1}}^{\mathcal{F}}({{\partial }_{{yi}}}p_{{\alpha 1}}^{\mathcal{F}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}}))dydx + $
$ + \;\int\limits_\Omega {\int\limits_\Gamma {w({{K}^{\mathcal{F}}}{\lambda }_{{\alpha 0}}^{{}}({{\partial }_{{yi}}}p_{{\alpha 3/2}}^{\mathcal{F}} + {{\partial }_{{xi}}}p_{{\alpha 1/2}}^{\mathcal{F}}) + } } {{K}^{\mathcal{F}}}{\lambda }_{{\alpha 1}}^{\mathcal{F}}({{\partial }_{{yi}}}p_{{\alpha 1}}^{\mathcal{F}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}}))n_{i}^{\Gamma }dxdy,$
где ${\mathbf{n}}_{{}}^{\Gamma } = \{ n_{1}^{\Gamma },n_{2}^{\Gamma },n_{3}^{\Gamma }\} $ – единичный нормальный вектор к границе Г, направленный от трещины к блоку.

С учетом разложений (2.7) и того факта, что ${\lambda }_{{\alpha 0}}^{{}}$ не зависит от y, это эквивалентно двум задачам в классической формулировке:

$\left\{ \begin{gathered} {{\partial }_{{yi}}}({{K}^{\mathcal{F}}}({{\partial }_{{yi}}}p_{{\alpha 1}}^{\mathcal{F}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}})) = 0,\quad y \in {{Y}^{\mathcal{F}}} \hfill \\ {{\left. {{{K}^{\mathcal{F}}}({{\partial }_{{yi}}}p_{{\alpha 1}}^{\mathcal{F}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}})n_{i}^{\Gamma }} \right|}_{{y \in \Gamma }}} = 0 \hfill \\ \end{gathered} \right.$
(3.2)
$\left\{ \begin{gathered} {{\partial }_{{yi}}}({{K}^{\mathcal{F}}}({{\partial }_{{yi}}}p_{{\alpha 3/2}}^{\mathcal{F}} + {{\partial }_{{xi}}}p_{{\alpha 1/2}}^{\mathcal{F}})) = 0,\quad y \in {{Y}^{\mathcal{F}}} \hfill \\ {{\left. {{{K}^{\mathcal{F}}}({{\partial }_{{yi}}}p_{{\alpha 3/2}}^{\mathcal{F}} + {{\partial }_{{xi}}}p_{{\alpha 1/2}}^{\mathcal{F}})n_{i}^{\Gamma }} \right|}_{{y \in \Gamma }}} = 0 \hfill \\ \end{gathered} \right.$

Это дает следующее представление для $p_{{\alpha 1}}^{\mathcal{F}}$ и $p_{{\alpha 3/2}}^{\mathcal{F}}$ через $p_{{\alpha 0}}^{{}}$ и $p_{{\alpha 1/2}}^{\mathcal{F}}$:

$p_{{\alpha 1}}^{} = {{\psi }_{k}}(y)\frac{{\partial {{p}_{{\alpha 0}}}}}{{\partial {{x}_{k}}}} + \bar {p}_{{\alpha 1}}^{\mathcal{F}}(x,t),\quad p_{{\alpha 3/2}}^{\mathcal{F}} = {{\psi }_{k}}(y)\frac{{\partial p_{{\alpha 1/2}}^{\mathcal{F}}}}{{\partial {{x}_{k}}}} + \bar {p}_{{\alpha 3/2}}^{\mathcal{F}}(x,t),$
где $\bar {p}_{{\alpha 1}}^{\mathcal{F}}(x,t)$ и $\bar {p}_{{\alpha 3/2}}^{\mathcal{F}}(x,t)$ – некоторые медленные функции, которые не войдут в осредненную модель. Для функций ${{\psi }_{k}}(y)$ получаем первую задачу на ячейке (2.17). Последнее условие в ней добавляется для единственности решения.

3.2. Второй шаг – осредненное уравнение на трещинах. Пусть $w = w\left( x \right)$ в $\Omega $ и ${{\left. w \right|}_{{\partial \Omega }}} = 0$. Тогда интегральное тождество (2.8) дает:

$\begin{gathered} \int\limits_\Omega {\int\limits_Y {w\phi ({{\partial }_{t}}{{s}_{{\alpha 0}}} + \sqrt \varepsilon {{\partial }_{t}}{{s}_{{\alpha 1/2}}})dxdy} } + \int\limits_\Omega {\int\limits_Y {w\phi {{{\beta }}_{\alpha }}({{s}_{{\alpha 0}}} + \sqrt \varepsilon {{s}_{{\alpha 1/2}}})({{\partial }_{t}}{{p}_{{\alpha 0}}} + \sqrt \varepsilon {{\partial }_{t}}{{p}_{{\alpha 1/2}}})dxdy} } = \\ = - \int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{F}}}} {{{K}^{\mathcal{F}}}({\lambda }_{{\alpha 0}}^{\mathcal{F}} + \sqrt \varepsilon {\lambda }_{{\alpha 1/2}}^{\mathcal{F}})({{\partial }_{{yi}}}{{p}_{{\alpha 1}}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}} + \sqrt \varepsilon ({{\partial }_{{yi}}}{{p}_{{\alpha 3/2}}} + {{\partial }_{{xi}}}{{p}_{{\alpha 1/2}}}))} } {{\partial }_{{xi}}}wdydx + \mathcal{O}\left( \varepsilon \right) \\ \end{gathered} $

Или с учетом задачи (3.2):

$\begin{gathered} \int\limits_\Omega {\int\limits_Y {w\phi {{\partial }_{t}}({{s}_{{\alpha 0}}} + \sqrt \varepsilon {{s}_{{\alpha 1/2}}})dxdy} } + \\ + \;\int\limits_\Omega {\int\limits_Y {w\phi {{{\beta }}_{\alpha }}({{s}_{{\alpha 0}}} + \sqrt \varepsilon {{s}_{{\alpha 1/2}}}){{\partial }_{t}}({{p}_{{\alpha 0}}} + \sqrt \varepsilon {{p}_{{\alpha 1/2}}})dxdy} } = \\ = \int\limits_\Omega {\int\limits_{{{Y}^{}}} {w{{\partial }_{{xi}}}({{K}^{}}({\lambda }_{{\alpha 0}}^{{}} + \sqrt \varepsilon {\lambda }_{{\alpha 1/2}}^{})({{\partial }_{{yi}}}{{\psi }_{k}} + {{\delta }_{{ik}}}){{\partial }_{{xk}}}({{p}_{{\alpha 0}}} + \sqrt \varepsilon p_{{\alpha 1/2}}^{\mathcal{F}}))dy} dx} + \mathcal{O}\left( \varepsilon \right) \\ \end{gathered} $

Вводя средние давления и насыщенности (2.9), получим:

$\begin{gathered} \int\limits_\Omega {w[{{\phi }^{\mathcal{F}}}\left( {1 - \theta } \right){{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{F}} + {{\phi }^{}}\theta {{\partial }_{t}}\mathcal{S}_{\alpha }^{}]dx} + \\ + \;\int\limits_\Omega {w[{{\phi }^{\mathcal{F}}}{\beta }_{\alpha }^{\mathcal{F}}\left( {1 - \theta } \right)\mathcal{S}_{\alpha }^{\mathcal{F}}{{\partial }_{t}}\mathcal{P}_{\alpha }^{\mathcal{F}} + {{\phi }^{}}{\beta }_{\alpha }^{}\theta \mathcal{S}_{\alpha }^{}{{\partial }_{t}}\mathcal{P}_{\alpha }^{}]dx} = \int\limits_\Omega {w{{\partial }_{{xi}}}({{\mathbb{K}}_{{ik}}}\lambda _{\alpha }^{\mathcal{F}}(\mathcal{S}_{\alpha }^{\mathcal{F}}){{\partial }_{{xk}}}\mathcal{P}_{\alpha }^{\mathcal{F}})dx} \\ \alpha = w,o, \\ \end{gathered} $
где эффективная проницаемость есть (2.16). Это дает первое осредненное уравнение системы (2.9) в почти законченном виде:

(3.3)
$\begin{gathered} {{\phi }^{\mathcal{F}}}\left( {1 - \theta } \right){{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{F}} + {{\phi }^{\mathcal{F}}}{\beta }_{\alpha }^{\mathcal{F}}\left( {1 - \theta } \right)\mathcal{S}_{\alpha }^{\mathcal{F}}{{\partial }_{t}}\mathcal{P}_{\alpha }^{\mathcal{F}} = \\ = {{\partial }_{{xi}}}({{\mathbb{K}}_{{ik}}}\lambda _{\alpha }^{\mathcal{F}}(\mathcal{S}_{\alpha }^{\mathcal{F}}){{\partial }_{{xk}}}\mathcal{P}_{\alpha }^{\mathcal{F}}) - {{\phi }^{\mathcal{M}}}\theta {{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{M}} - {{\phi }^{\mathcal{M}}}{\beta }_{\alpha }^{\mathcal{M}}\theta \mathcal{S}_{\alpha }^{\mathcal{M}}{{\partial }_{t}}\mathcal{P}_{\alpha }^{\mathcal{M}} \\ \end{gathered} $

3.3. Третий шаг – разложение на блоках. Пусть функция $w = w\left( {x,y} \right)$ – непрерывная функция обоих аргументов, действующая в пространстве $\Omega \times Y$, ${{\left. w \right|}_{{\partial \Omega }}} = 0$, такая, что выполняется условие $w \equiv 0$ в ${{Y}^{\mathcal{F}}}$. Тогда для членов нулевого порядка интегральное тождество (2.8) дает:

$\begin{gathered} \int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{M}}}} {w\phi ({{\partial }_{t}}{{s}_{{\alpha 0}}} + \sqrt \varepsilon {{\partial }_{t}}s_{{\alpha 1/2}}^{\mathcal{M}})dxdy} } + \\ + \;\int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{M}}}} {w\phi {{{\beta }}_{\alpha }}({{s}_{{\alpha 0}}} + \sqrt \varepsilon s_{{\alpha 1/2}}^{\mathcal{M}})({{\partial }_{t}}{{p}_{{\alpha 0}}} + \sqrt \varepsilon {{\partial }_{t}}p_{{\alpha 1/2}}^{\mathcal{M}})dxdy} } = \\ = - \int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{M}}}} {{{K}^{\mathcal{M}}}({\lambda }_{{\alpha 0}}^{{}} + \sqrt \varepsilon {\lambda }_{{\alpha 1/2}}^{\mathcal{M}})({{\partial }_{{yi}}}p_{{\alpha 1/2}}^{\mathcal{M}} + \sqrt \varepsilon ({{\partial }_{{yi}}}p_{{\alpha 1/2}}^{\mathcal{M}} + {{\partial }_{{xi}}}{{p}_{{\alpha 0}}})){{\partial }_{{yi}}}wdy} dx + \mathcal{O}\left( \varepsilon \right)} \\ \end{gathered} $

Интеграл по границе блока $\Gamma $ равен нулю, так как функции $w$ равны нулю на трещине и непрерывны.

Члены нулевого порядка дают:

$\int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{M}}}} {w\phi \left( {{{\partial }_{t}}{{s}_{{\alpha 0}}} + {{{\beta }}_{\alpha }}{{s}_{{\alpha 0}}}{{\partial }_{t}}{{p}_{{\alpha 0}}}} \right)dydx} } = \int\limits_\Omega {\int\limits_{{{Y}^{\mathcal{M}}}} {w{{\partial }_{{yi}}}({{K}^{\mathcal{M}}}{\lambda }_{{\alpha 0}}^{{}}{{\partial }_{{yi}}}p_{{\alpha 1/2}}^{\mathcal{M}})dy} dx} ,$
что определяет следующее выражение для $p_{{\alpha 1/2}}^{\mathcal{M}}$ с учетом (2.6):
(3.4)
$p_{{\alpha 1/2}}^{\mathcal{M}} = p_{{\alpha 1/2}}^{\mathcal{F}}\left( {x,t} \right) - \varphi \left( y \right)\frac{{{{\phi }^{\mathcal{M}}}}}{{{\lambda }_{{\alpha 0}}^{{}}}}[{{\partial }_{t}}{{s}_{{\alpha 0}}} + {\beta }_{\alpha }^{\mathcal{M}}{{s}_{{\alpha 0}}}{{\partial }_{t}}{{p}_{{\alpha 0}}}],$
где $\varphi \left( y \right)$ – решение второй задачи на ячейке (2.21). Граничное условие в (2.21) вытекает из условия непрерывности фазовых давлений (2.6).

3.4. Четвертый шаг – осредненное уравнение на блоках. Соотношение (3.4) дает возможность получить явное соотношение между функциями $\mathcal{P}_{\alpha }^{\mathcal{M}}$ и $\mathcal{P}_{\alpha }^{\mathcal{F}}$. В самом деле, беря среднее от (3.4) по ${{Y}^{\mathcal{M}}}$ (и учитывая, что в правой части (3.4) только функция $\varphi \left( y \right)$ зависит от $y$), умножая на $\sqrt \varepsilon $ и прибавляя ${{p}_{{\alpha 0}}}$, получим:

${{p}_{{\alpha 0}}} + \sqrt \varepsilon {{\left\langle {p_{{\alpha 1/2}}^{\mathcal{M}}} \right\rangle }_{\mathcal{M}}} = {{p}_{{\alpha 0}}} + \sqrt \varepsilon p_{{\alpha 1/2}}^{\mathcal{F}} - \sqrt \varepsilon {{\left\langle \varphi \right\rangle }_{\mathcal{M}}}\frac{{{{\phi }^{\mathcal{M}}}}}{{{\lambda }_{{\alpha 0}}^{{}}}}[{{\partial }_{t}}{{s}_{{\alpha 0}}} + {\beta }_{\alpha }^{\mathcal{M}}{{s}_{{\alpha 0}}}{{\partial }_{t}}{{p}_{{\alpha 0}}}]$

Учитывая определение средних давлений и насыщенностей (2.9), получим:

(3.5)
$\mathcal{P}_{\alpha }^{\mathcal{M}} = \mathcal{P}_{\alpha }^{\mathcal{F}} - \frac{{\tau _{\alpha }^{{{\text{comp}}}}}}{{\beta _{\alpha }^{\mathcal{M}}\mathcal{S}_{\alpha }^{\mathcal{M}}}}[{{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{M}} + \beta _{\alpha }^{\mathcal{M}}\mathcal{S}_{\alpha }^{\mathcal{M}}{{\partial }_{t}}\mathcal{P}_{\alpha }^{\mathcal{M}}],\quad \alpha = w,o,$
что совпадает с уравнением (2.14). Характерные времена запаздывания $\tau _{\alpha }^{{{\text{comp}}}}$ определены выражением (2.18).

В структуре времен запаздывания было учтено соотношение

${{\partial }_{t}}{{s}_{{\alpha 0}}} = {{\partial }_{t}}({{s}_{{\alpha 0}}} + \sqrt \varepsilon {{\left\langle {s_{{\alpha 1/2}}^{\mathcal{M}}} \right\rangle }_{\mathcal{M}}}) - \sqrt \varepsilon {{\partial }_{t}}{{\left\langle {s_{{\alpha 1/2}}^{\mathcal{M}}} \right\rangle }_{\mathcal{M}}} = {{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{M}} + \mathcal{O}(\sqrt \varepsilon )$
$\begin{gathered} {\lambda }_{{\alpha 0}}^{\mathcal{M}} \equiv {\lambda }_{\alpha }^{\mathcal{M}}\left( {{{s}_{{\alpha 0}}}} \right) = {\lambda }_{\alpha }^{\mathcal{M}}\left( {{{s}_{{\alpha 0}}}} \right) + \sqrt \varepsilon \frac{{d{\lambda }_{\alpha }^{\mathcal{M}}}}{{d{{s}_{{\alpha 0}}}}}{{\left\langle {s_{{\alpha 1/2}}^{\mathcal{M}}} \right\rangle }_{\mathcal{M}}} - \sqrt \varepsilon \frac{{d{\lambda }_{\alpha }^{\mathcal{M}}}}{{d{{s}_{{\alpha 0}}}}}{{\left\langle {s_{{\alpha 1/2}}^{}} \right\rangle }_{\mathcal{M}}}{\text{ = }} \\ = {\lambda }_{\alpha }^{\mathcal{M}}({{s}_{{\alpha 0}}} + \sqrt \varepsilon {{\left\langle {s_{{\alpha 1/2}}^{\mathcal{M}}} \right\rangle }_{\mathcal{M}}}) + \mathcal{O}(\sqrt \varepsilon ) = {\lambda }_{\alpha }^{\mathcal{M}}(\mathcal{S}_{\alpha }^{\mathcal{M}}) + \mathcal{O}(\sqrt \varepsilon ) \\ \end{gathered} $

В результате чего получаем:

$\frac{{\sqrt \varepsilon }}{{{\lambda }_{{\alpha 0}}^{\mathcal{M}}}}{{\partial }_{t}}{{s}_{{\alpha 0}}} = \frac{{\sqrt \varepsilon {{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{M}} + \mathcal{O}\left( \varepsilon \right)}}{{{\lambda }_{\alpha }^{\mathcal{M}}(\mathcal{S}_{\alpha }^{\mathcal{M}}) + \mathcal{O}(\sqrt \varepsilon )}} = \frac{{\sqrt \varepsilon {{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{M}}}}{{{\lambda }_{\alpha }^{\mathcal{M}}(\mathcal{S}_{\alpha }^{\mathcal{M}})}} + \mathcal{O}\left( \varepsilon \right)$
и т.п. для других членов.

Вместо двух уравнений (3.5) для фазовых давлений удобнее заменить одно из них уравнением для капиллярного давления. Вычитая одно уравнение (3.5) из другого, получаем уравнение связи средних капиллярных давлений в блоках и трещинах:

(3.6)
$\begin{gathered} \mathcal{P}_{с}^{\mathcal{M}} = \mathcal{P}_{с}^{} + \tau _{w}^{{{\text{comp}}}}{{\partial }_{t}}\mathcal{P}_{w}^{\mathcal{M}} - \tau _{o}^{{{\text{comp}}}}{{\partial }_{t}}(\mathcal{P}_{w}^{\mathcal{M}} + \mathcal{P}_{c}^{\mathcal{M}}) + \left( {\frac{{\tau _{w}^{{{\text{comp}}}}}}{{\beta _{w}^{\mathcal{M}}\mathcal{S}_{w}^{\mathcal{M}}}} + \frac{{\tau _{o}^{{{\text{comp}}}}}}{{\beta _{o}^{\mathcal{M}}\mathcal{S}_{o}^{\mathcal{M}}}}} \right){{\partial }_{t}}\mathcal{S}_{w}^{\mathcal{M}} = \\ = \mathcal{P}_{с}^{} + \left( {\frac{{\tau _{w}^{{{\text{comp}}}}}}{{\beta _{w}^{\mathcal{M}}\mathcal{S}_{w}^{\mathcal{M}}}} + \frac{{\tau _{o}^{{{\text{comp}}}}}}{{\beta _{o}^{\mathcal{M}}\mathcal{S}_{o}^{\mathcal{M}}}} - \tau _{o}^{{{\text{comp}}}}\frac{{d\mathcal{P}_{c}^{\mathcal{M}}}}{{d\mathcal{S}_{w}^{\mathcal{M}}}}} \right)\frac{{\partial{ \mathcal{S}}_{w}^{\mathcal{M}}}}{{\partial t}} + (\tau _{w}^{{{\text{comp}}}} - \tau _{o}^{{{\text{comp}}}})\frac{{\partial{ \mathcal{P}}_{w}^{\mathcal{M}}}}{{\partial t}}, \\ \end{gathered} $
что сводится к уравнению (2.15), если связать средние капиллярные давления $\mathcal{P}_{с}^{\mathcal{F}}$ и $\mathcal{P}_{с}^{\mathcal{M}}$ с исходными заданными кривыми капиллярного давления $p_{с}^{\mathcal{F}}\left( s \right)$ и $p_{с}^{\mathcal{M}}\left( s \right)$, то есть получить соотношения (2.12). Используя разложения выражений для $p_{c}^{\mathcal{F}}(\mathcal{S}_{w}^{\mathcal{F}})$ и$p_{c}^{\mathcal{M}}(\mathcal{S}_{w}^{\mathcal{M}})$ в ряд Тейлора, учитывая выражения (2.10), приходим к соотношениям

$p_{c}^{\mathcal{F}}(\mathcal{S}_{w}^{\mathcal{F}}) = p_{c}^{\mathcal{F}}\left( {{{s}_{{w0}}}} \right){\text{ + }}\sqrt \varepsilon \frac{{dp_{\alpha }^{\mathcal{F}}}}{{d{{s}_{{w0}}}}}s_{{w1/2}}^{\mathcal{F}} + \mathcal{O}\left( \varepsilon \right)$
$p_{c}^{\mathcal{M}}(\mathcal{S}_{w}^{\mathcal{M}}) = p_{c}^{\mathcal{M}}\left( {{{s}_{{w0}}}} \right){\text{ + }}\sqrt \varepsilon \frac{{dp_{\alpha }^{\mathcal{M}}}}{{d{{s}_{{w0}}}}}{{\left\langle {s_{{w1/2}}^{\mathcal{M}}} \right\rangle }_{}} + \mathcal{O}\left( \varepsilon \right)$

Поскольку все функции справа не зависят от $y$, получаем:

$\begin{gathered} \mathcal{P}_{c}^{\mathcal{F}}(\mathcal{S}_{w}^{\mathcal{F}}) \equiv {{\left\langle {p_{c}^{\mathcal{F}}(\mathcal{S}_{w}^{\mathcal{F}})} \right\rangle }_{}} = p_{c}^{\mathcal{F}}(\mathcal{S}_{w}^{\mathcal{F}}) + \mathcal{O}\left( \varepsilon \right) \\ \mathcal{P}_{c}^{\mathcal{M}}(\mathcal{S}_{w}^{\mathcal{M}}) \equiv {{\left\langle {p_{c}^{\mathcal{M}}(\mathcal{S}_{w}^{\mathcal{M}})} \right\rangle }_{\mathcal{M}}} = p_{c}^{\mathcal{M}}(\mathcal{S}_{w}^{\mathcal{M}}) + \mathcal{O}\left( \varepsilon \right), \\ \end{gathered} $
что подтверждает соотношения (2.12). Найденные соотношения дают возможность записать формулу (3.6) через исходные заданные функции капиллярного давления: (2.15).

Для выражений в правой части макроскопических уравнений (3.3) можно получить явное выражение через разницу давлений, используя уравнение (2.14):

${{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{M}} + \beta _{\alpha }^{\mathcal{M}}\mathcal{S}_{\alpha }^{\mathcal{M}}{{\partial }_{t}}\mathcal{P}_{\alpha }^{\mathcal{M}} = - \frac{{\beta _{\alpha }^{\mathcal{M}}\mathcal{S}_{\alpha }^{\mathcal{M}}}}{{\tau _{\alpha }^{{{\text{comp}}}}}}(\mathcal{P}_{\alpha }^{\mathcal{M}} - \mathcal{P}_{\alpha }^{\mathcal{F}}),$
что дает окончательную форму осредненных уравнений (2.13).

4. Физический анализ модели. Уравнения (2.13) описывают поведение насыщенности $\mathcal{S}_{{}}^{\mathcal{F}}$и давления $\mathcal{P}_{{}}^{\mathcal{F}}$в трещинах с обменными членами справа, пропорциональными разнице фазовых давлений. Соотношения (2.14), (2.15) представляют собой систему двух обыкновенных дифференциальных уравнений для насыщенности $\mathcal{S}_{{}}^{\mathcal{M}}$ и давления $\mathcal{P}_{{}}^{\mathcal{M}}$ в блоках (при известных насыщенности и давлении в трещинах).

4.1. Частный случай однофазной сжимаемой системы. В однофазном случае система (2.13)–(2.15) принимает вид хорошо известной модели двойной пористости [2, 3]:

(4.1)
${{\phi }^{\mathcal{F}}}{\beta }_{{}}^{\mathcal{F}}\left( {1 - \theta } \right){{\partial }_{t}}\mathcal{P}_{{}}^{\mathcal{F}} - {{\partial }_{{xi}}}({{\mathbb{K}}_{{ik}}}{{\partial }_{{xk}}}\mathcal{P}_{{}}^{\mathcal{F}}) = \xi (\mathcal{P}_{{}}^{\mathcal{M}} - \mathcal{P}_{{}}^{\mathcal{F}})$
(4.2)
$\mathcal{P}_{{}}^{\mathcal{M}} - \mathcal{P}_{{}}^{\mathcal{F}} = - \tau _{{}}^{{{\text{comp}}}}\frac{{\partial{ \mathcal{P}}_{{}}^{\mathcal{M}}}}{{\partial t}},$
где $\xi \equiv \frac{\theta }{{\sqrt \varepsilon {{{\left\langle \varphi \right\rangle }}_{\mathcal{M}}}}}$, $\tau _{{}}^{{{\text{comp}}}} = \sqrt \varepsilon {{\left\langle \varphi \right\rangle }_{\mathcal{M}}}{{\phi }^{\mathcal{M}}}{\beta }_{{}}^{\mathcal{M}}\mu $, $\mu $ – динамическая вязкость жидкости.

4.2. Частный случай несжимаемой двухфазной системы. В случае, когда фазы несжимаемы, система (2.13)–(2.15) принимает вид:

(4.3)
${{\phi }^{\mathcal{F}}}\left( {1 - \theta } \right){{\partial }_{t}}\mathcal{S}_{\alpha }^{\mathcal{F}} - {{\partial }_{{xi}}}({{\mathbb{K}}_{{ik}}}\lambda _{\alpha }^{\mathcal{F}}(\mathcal{S}_{\alpha }^{\mathcal{F}}){{\partial }_{{xk}}}\mathcal{P}_{\alpha }^{\mathcal{F}}) = {{\xi }_{\alpha }}(\mathcal{P}_{\alpha }^{} - \mathcal{P}_{\alpha }^{\mathcal{F}}),\quad \alpha = w,o$
(4.4)
$\mathcal{P}_{w}^{\mathcal{M}} = \mathcal{P}_{w}^{\mathcal{F}} - \frac{{\tau _{w}^{{{\text{comp}}}}}}{{\beta _{w}^{\mathcal{M}}\mathcal{S}_{w}^{\mathcal{M}}}}{{\partial }_{t}}\mathcal{S}_{w}^{\mathcal{M}}$
(4.5)
$p_{с}^{\mathcal{M}}(\mathcal{S}_{w}^{\mathcal{M}}) = p_{с}^{} + \frac{{\tau _{{}}^{{{\text{cap}}}}}}{{{\beta }_{o}^{\mathcal{M}}}}{{\partial }_{t}}\mathcal{S}_{w}^{\mathcal{M}},$
где отношения $\frac{{\tau _{w}^{{{\text{comp}}}}}}{{\beta _{w}^{\mathcal{M}}\mathcal{S}_{{}}^{\mathcal{M}}}}$ и $\frac{{\tau _{{}}^{{{\text{cap}}}}}}{{{\beta }_{o}^{\mathcal{M}}}}$ не зависят от коэффициентов сжимаемости. В самом деле, согласно соотношениям (2.18), (2.19) получаем:

$\frac{{\tau _{w}^{{{\text{comp}}}}}}{{\beta _{w}^{\mathcal{M}}\mathcal{S}_{{}}^{\mathcal{M}}}} = \sqrt \varepsilon {{\left\langle \varphi \right\rangle }_{\mathcal{M}}}\frac{{{{\phi }^{\mathcal{M}}}}}{{\lambda _{w}^{\mathcal{M}}}},\quad \frac{{\tau _{{}}^{{{\text{cap}}}}}}{{{\beta }_{o}^{\mathcal{M}}}} = \sqrt \varepsilon {{\left\langle \varphi \right\rangle }_{\mathcal{M}}}{{\phi }^{\mathcal{M}}}\frac{{(\lambda _{w}^{\mathcal{M}} + \lambda _{o}^{\mathcal{M}})}}{{\lambda _{w}^{\mathcal{M}}\lambda _{o}^{\mathcal{M}}}}$

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

Эта модель совпадает с моделью [3, 8] для несжимаемого течения в недеформируемой среде.

4.3. Капиллярное запаздывание и запаздывание, вызванное сжимаемостью. Сравнение со случаем несжимаемых флюидов (4.2) и однофазного течения (4.1) позволяет лучше понять суть полученной модели (2.13)–(2.15) и роль сжимаемости в двухфазной системе.

Неравновесное поведение системы (2.13)–(2.15) определяется подсистемой двух ОДУ (2.14), (2.15) для насыщенности $\mathcal{S}_{w}^{\mathcal{M}}$ и давления $\mathcal{P}_{w}^{\mathcal{M}}$ в блоках. Связь решения на блоках и трещинах через дифференциальные уравнения по времени указывает на неравновесный характер процесса на макроуровне. Что при этом происходит на микроуровне, какие конкретно процессы вызывают эту неравновесность? Прежде всего, выделим два основных механизма:

Неравновесное капиллярное перераспределение фаз на макромасштабе, вызванное тем, что средняя насыщенность в блоках меняется медленнее, чем в трещинах, что нарушает равенство средних капиллярных давлений (условие капиллярного равновесия). В итоге, возникает разница средних капиллярных давлений, которая зависит от скорости изменения насыщенности в блоке. В уравнении (2.15) этот процесс описывается первым слагаемым в правой части: $p_{с}^{\mathcal{M}} - p_{с}^{\mathcal{F}}$ = $\frac{{\tau _{{}}^{{{\text{cap}}}}}}{{{\beta }_{o}^{\mathcal{M}}}}\frac{{\partial{ \mathcal{S}}_{w}^{\mathcal{M}}}}{{\partial t}}$. Разница в капиллярных давлениях автоматически вызывает и разницу в фазовых давлениях, которая описывается первым слагаемым в правой части уравнения (2.14): $\mathcal{P}_{{}}^{\mathcal{M}} - \mathcal{P}_{{}}^{\mathcal{F}}$ = $ - \frac{{\tau _{w}^{{{\text{comp}}}}}}{{\beta _{w}^{\mathcal{M}}\mathcal{S}_{w}^{\mathcal{M}}}}\frac{{\partial{ \mathcal{S}}_{w}^{\mathcal{M}}}}{{\partial t}}$. В сильно контрастной среде такой процесс привел бы к долговременной памяти. В рассматриваемом же случае умеренного контраста времена запаздывания малы. Этот процесс не зависит от сжимаемости фаз и пород, и выражение для разницы давлений совпадает с формулами (4.4), (4.5), полученными для несжимаемых фаз.

Релаксация волн давления в сжимаемой среде, вызванная тем, что возмущение давления в низкопроницаемых блоках распространяется много медленнее, чем в трещинах. Возникает разница средних давлений, которая описывается вторым слагаемым в правой части уравнения (2.14): $\mathcal{P}_{{}}^{\mathcal{M}} - \mathcal{P}_{{}}^{\mathcal{F}}$ = $ - \tau _{w}^{{{\text{comp}}}}\frac{{\partial{ \mathcal{P}}_{{}}^{\mathcal{M}}}}{{\partial t}}$. Тогда процесс релаксации не зависит от двухфазности системы, а его уравнение совпадает с уравнением (4.1) для однофазной системы.

Механизмы неравновесности, связанные с перекрестными эффектами наложения сжимаемости и капиллярности, можно описать следующим образом.

Перераспределение фаз из-за асимметричного отжима (выдавливания). Сжимаемость приводит к расширению жидкостей и компакции пор под весом вышележащих пород (если давление падает), что ведет к «выдавливанию» обеих фаз из пор. В геонауках такой эффект называют «отжимом» жидкости из пор. Этот эффект похож на перистальтику, когда жидкость в канале выдавливается из него за счет резкого сжатия канала. Если такое выдавливание происходит симметрично для обеих фаз, то они отжимаются как единое целое, что не меняет их насыщенность. В случае же асимметричного отжима происходит перераспределение фаз в порах, что влияет на насыщенность. В блоках это движение, вызванное отжимом, запаздывает, что вносит дополнительную неравновесность в поведение насыщенности. Возникает дополнительная разница в капиллярных давлениях, которая описывается третьим членом в (2.15): $p_{с}^{\mathcal{M}} - p_{с}^{\mathcal{F}}$ = $(\tau _{w}^{{{\text{comp}}}}$$\tau _{o}^{{{\text{comp}}}})\frac{{\partial{ \mathcal{P}}_{w}^{\mathcal{M}}}}{{\partial t}}$. Как видно, этот эффект действительно отсутствует при симметричном отжиме фаз, то есть когда $\tau _{w}^{{{\text{comp}}}}$ = $\tau _{o}^{{{\text{comp}}}}$.

Перераспределение фаз из-за нелинейного отжима, вызванное наложением капиллярности и сжимаемости, при учете нелинейности капиллярных явлений. Этот эффект описывается вторым слагаемым в правой части уравнения (2.15): $p_{с}^{\mathcal{M}} - p_{с}^{\mathcal{F}}$ = = $\frac{{\tau _{{}}^{{{\text{cc}}}}}}{{{\beta }_{o}^{\mathcal{M}}}}\frac{{\partial{ \mathcal{S}}_{w}^{\mathcal{M}}}}{{\partial t}}$.

5. Подход к описанию сильно сжимаемых пористых сред. Полученные результаты относятся к слабосжимаемой пористой среде, так что коэффициенты сжимаемости пор входят в структуру исходных уравнений, но пористость и проницаемость в них считаются от давления зависящими очень слабо. Насколько можно снять это условие слабой сжимаемости? Допустим, что среда может сжиматься сильно, и, что при этом не меняются фазовые проницаемости и кривые капиллярного давления блоков и трещин. Тогда есть все основания полагать, что сильная сжимаемость среды не повлияет на структуру асимптотических разложений в первых его членах, а значит, полученные макроскопические уравнения (2.13)(2.15) формально не изменятся. Однако, тогда исходные пористости ${{\phi }^{\mathcal{F}}}$и ${{\phi }^{\mathcal{M}}}$, а также проницаемости ${{K}^{\mathcal{F}}}$ и ${{K}^{\mathcal{M}}}$ зависят от давления. В результате тензор эффективной проницаемости (2.16), (2.17) также зависит от давления через проницаемость трещин ${{K}^{\mathcal{F}}}$. Времена запаздывания (2.13) тоже будут зависеть от давления через пористость блоков ${{\phi }^{\mathcal{M}}}$ и проницаемость блоков ${{K}^{\mathcal{M}}}$, которая входит в ячеечную задачу (2.21). Более того, поскольку давления в воде и нефти разные, то проницаемости ${{K}^{\mathcal{F}}}$ и ${{K}^{\mathcal{M}}}$ и пористость ${{\phi }^{\mathcal{M}}}$ ведут себя по-разному в порах, занятых водой и нефтью. Это можно учесть, если считать их зависящими от давлений, осредненных по фазам:

(5.1)
${{\mathcal{P}}^{\mathcal{F}}} \equiv \mathcal{P}_{w}^{\mathcal{F}}\mathcal{S}_{w}^{\mathcal{F}} + \mathcal{P}_{o}^{\mathcal{F}}(1 - \mathcal{S}_{w}^{\mathcal{F}})\quad {\text{и}}\quad {{\mathcal{P}}^{\mathcal{M}}} \equiv \mathcal{P}_{w}^{\mathcal{M}}\mathcal{S}_{w}^{\mathcal{M}} + \mathcal{P}_{o}^{\mathcal{M}}(1 - \mathcal{S}_{w}^{\mathcal{M}})$

Таким образом, возникает ситуация, в которой ячеечные задачи (2.17) и (2.21) зависят от макроскопических переменных $\mathcal{P}_{w}^{\mathcal{M}}$, $\mathcal{P}_{w}^{\mathcal{F}}$, $\mathcal{P}_{o}^{\mathcal{F}}$, $\mathcal{P}_{o}^{\mathcal{F}}$, $\mathcal{S}_{w}^{\mathcal{M}}$, $\mathcal{S}_{w}^{\mathcal{F}}$, что является типичным случаем не полностью осредненной модели. И хотя ее можно использовать для качественного анализа процесса, но для численных расчетов она не проще, чем исходная неосредненная система уравнений.

Ситуация изменяется, если ячеечные задачи можно решить аналитически. Тогда макроскопические давления и насыщенности появятся в структуре эффективной проницаемости и времен релаксации в явном виде. Это можно сделать для простой геометрии блоков ${{Y}^{\mathcal{M}}}$. Например, для сферических блоков радиуса $R$, ячеечная задача (2.21) имеет точное аналитическое решение, а ячеечная задача (2.17) имеет приближенное решение, полученное методом потенциальных течений [3]:

$\varphi (y) = \frac{{{{R}^{2}} - {{r}^{2}}}}{{6{{K}^{\mathcal{M}}}}},\quad {{\psi }_{k}}(y) \approx \frac{{{{y}_{k}}\left[ {2 + {{{\left( {\frac{R}{r}} \right)}}^{3}}} \right]}}{{2 + \theta }} - {{y}_{k}}\quad \left( {k = 1,2,3} \right);\quad r \equiv \sqrt {y_{1}^{2} + y_{2}^{2} + y_{3}^{2}} ,$
где $R = l{{\left( {\frac{{3\theta }}{{4\pi }}} \right)}^{{1/3}}}$, а $l$ – размерная длина одного периода неоднородности среды. Тогда можно получить:

(5.2)
${{\left\langle \varphi \right\rangle }_{\mathcal{M}}} = \frac{{{{l}^{2}}}}{{15{{K}^{\mathcal{M}}}}}{{\left( {\frac{{3\theta }}{{4\pi }}} \right)}^{{2/3}}},\quad \mathbb{K} \approx 2K_{{}}^{\mathcal{F}}\frac{{1 - \theta }}{{2 + \theta }}$

Эти соотношения должны быть дополнены законом сжимаемости, записанным относительно проницаемостей ${{K}^{\mathcal{F}}}$ и ${{K}^{\mathcal{M}}}$, например, экспоненциальным законом типа (1.1):

(5.3)
${{K}^{\mathcal{M}}} = {{K}^{{\mathcal{M},0}}}{{e}^{{\tilde {\beta }_{K}^{\mathcal{M}}({{}^{\mathcal{M}}} - {{P}^{0}})}}},\quad {{K}^{\mathcal{F}}} = {{K}^{{\mathcal{F},0}}}{{e}^{{\tilde {\beta }_{K}^{\mathcal{F}}({{}^{\mathcal{F}}} - {{P}^{0}})}}},$
где $\tilde {\beta }_{K}^{\mathcal{F}}$ и $\tilde {\beta }_{K}^{\mathcal{F}}$ коэффициенты сжимаемости по проницаемости, давления ${{\mathcal{P}}^{\mathcal{F}}}$ и ${{\mathcal{P}}^{\mathcal{M}}}$ определены в (5.1), а верхним индексом “ноль” обозначена некая характерная величина.

Теперь модель (2.13)–(2.15) с соотношениями (5.2), (5.3) и (5.1) полностью осреднена, но существенно нелинейна по давлению.

В частности, для однофазного случая получаем замкнутую нелинейную систему течения (сильно-) сжимаемой жидкости в (сильно-) сжимаемой среде:

(5.4)
${\beta }_{{}}^{\mathcal{F}}\left( {1 - \theta } \right){{\phi }^{\mathcal{F}}}(\mathcal{P}_{{}}^{\mathcal{F}}){{\partial }_{t}}\mathcal{P}_{{}}^{\mathcal{F}} - {{\partial }_{{xi}}}(\mathbb{K}(\mathcal{P}_{{}}^{\mathcal{F}}){{\partial }_{{xi}}}\mathcal{P}_{{}}^{\mathcal{F}}) = \xi (\mathcal{P}_{{}}^{\mathcal{M}})(\mathcal{P}_{{}}^{\mathcal{M}} - \mathcal{P}_{{}}^{\mathcal{F}})$
$\mathcal{P}_{{}}^{\mathcal{M}} - \mathcal{P}_{{}}^{\mathcal{F}} = - \tau _{{}}^{{{\text{comp}}}}(\mathcal{P}_{{}}^{\mathcal{M}})\frac{{\partial{ \mathcal{P}}_{{}}^{\mathcal{M}}}}{{\partial t}}$
(5.5)
$\xi = \frac{{15{{\theta }^{{1/3}}}{{{\left( {\frac{{4\pi }}{3}} \right)}}^{{2/3}}}}}{{\sqrt \varepsilon {{l}^{2}}}}{{K}^{\mathcal{M}}}({{\mathcal{P}}^{\mathcal{M}}}),\quad \mathbb{K} \approx 2\frac{{1 - \theta }}{{2 + \theta }}K_{{}}^{\mathcal{F}}(\mathcal{P}_{{}}^{\mathcal{F}})$
$\tau _{{}}^{{{\text{comp}}}} = \frac{{\sqrt \varepsilon }}{{15}}{\beta }_{{}}^{\mathcal{M}}{{l}^{2}}\mu {{\left( {\frac{{3\theta }}{{4\pi }}} \right)}^{{2/3}}}\frac{{{{\phi }^{\mathcal{M}}}(\mathcal{P}_{{}}^{\mathcal{M}})}}{{K_{{}}^{\mathcal{M}}(\mathcal{P}_{{}}^{\mathcal{M}})}},$
где ${\beta }_{{}}^{\mathcal{F}}$, ${\beta }_{{}}^{\mathcal{M}}$, $\theta $, $\mu $, $l$, $\varepsilon $ – константы, а зависимости ${{\phi }^{\mathcal{M}}}(\mathcal{P}_{{}}^{\mathcal{M}})$, $\phi _{{}}^{\mathcal{F}}(\mathcal{P}_{{}}^{\mathcal{F}})$, $K_{{}}^{\mathcal{M}}(\mathcal{P}_{{}}^{\mathcal{M}})$, $K_{{}}^{\mathcal{F}}(\mathcal{P}_{{}}^{\mathcal{F}})$ описываются соотношениями типа (5.3).

Заключение. Построенная осредненная модель (2.13)–(2.15) представляет теоретический интерес. Во-первых, она полностью осреднена и в ней отсутствуют микроскопические переменные, а задачи ячейки (2.17), (2.21) не содержат макроскопических переменных, поэтому решаются только один раз.

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

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

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

Работа выполнена при поддержке комитета Науки министерства образования и науки Республики Казахстан: грант № AP05132680.

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

  1. Баренблатт Г.И., Желтов Ю.П., Кочина И.Н. Об основных представлениях теории фильтрации однородных жидкостей в трещиноватых породах // ПММ. 1960. Т. 24. № 5. С. 852–864.

  2. Arbogast T., Douglas J., Hornung U. Derivation of the double porosity model of single phase flow via homogenization theory // SIAM J. Math. Anal. 1990. V. 21. № 4. P. 823–836.

  3. Panfilov M. Macroscale Models of Flow through Highly Heterogeneous Porous Media. Dordrecht: Kluwer Acad. Publ., 2000. 363 p.

  4. Amaziane B., Jurak M., Pankratov L., Vrbaski A. Some remarks on the homogenization of immiscible incompressible two-phase flow in double porosity media // Discr. Cont. Dyn. Syst. Ser. B. 2018. V. 23. № 2. P. 629–665.

  5. Amaziane B., Milisic J.P., Panfilov M., Pankratov L. Generalized nonequilibrium capillary relations for two-phase flow through heterogeneous media // Phys. Rev. E. 2012. V. 85. P. 016304.

  6. Amaziane B., Pankratov L., Jurak M., Vrbaski A. A fully homogenized model for incompressible two-phase flow in double porosity media // Appl. Anal. 2015. April.

  7. Bourgeat A., Luckhaus S., Mikelic A. Convergence of the homogenization process for a double-porosity model of immiscible two-phase flow // SIAM J. Math. Anal. 1966. V. 27. № 6. P. 1520–1543.

  8. Bourgeat A., Panfilov M. Effective two-phase flow through highly heterogeneous porous media // Comput. Geosci. 1998. V. 2. P. 191–215.

  9. Yeh L.M. Homogenization of two-phase flow in fractured media // Math. Models Meth. Appl. Sci. 2006. V. 16. P. 1627–1651.

  10. Arbogast T. A simplified dual porosity model for two-phase flow // in: Computat. Meth. in Water Res. IX. V. 2: Math. Model. Water Res. Eds. T.F. Russell et al. Southampton: Comput. Mech. Publ., 1992.

  11. Ait Mahiout L., Amaziane B., Mokrane A., Pankratov L. Homogenization of immiscible compressible twophase flow in double porosity media // Electron. J. Diff. Equat. 2016. V. 52. P. 1–28.

  12. Amaziane B., Pankratov L. Homogenization of a model for water-gas flow through double-porosity media // Math. Meth. Appl. Sci. 2016. V. 39. P. 425–451.

  13. Jafari I., Masihi M., Nasiri Zarandi M. Experimental study on imbibition displacement mechanisms of two-phase fluid using micromodel: Fracture network, distribution of pore size, and matrix construction // Phys. Fluids. 2017. V. 29. № 11. P. 122004.

  14. Khoshkalam Y., Khosravi M., Rostami B. Visual investigation of viscous cross-flow during foam injection in a matrix-fracture system // Phys. Fluids. 2019. V. 31. P. 023102.

  15. Yao C.C., Yan P.Y. A diffuse interface approach to injection-driven flow of different miscibility in heterogeneous porous media // Phys. Fluids. 2015. V. 27. № 8. P. 083101.

  16. Li H., Guo H., Yang Z. et al. Evaluation of oil production potential in fractured porous media // Phys. Fluids. 2019. V. 31. P. 052104.

  17. Jafari I., Masihi M., Nasiri Zarandi M. Numerical simulation of counter-current spontaneous imbibition in water-wet fractured porous media: Influences of water injection velocity, fracture aperture, and 826 grains geometry // Phys. Fluids. 2017. V. 29. № 11. P. 113305.

  18. Panfilov M. Macroscale Models of Flow through Highly Heterogeneous Porous Media. Dordrecht: Kluwer Acad. Publ., 2000. 363 p.

  19. Hussein M. Multiphase flow simulations in heterogeneous fractured media through hybrid grid method // AIP Conf. Proc. 2013. № 1558. P. 2048.

  20. Saedi B., Ayatollahi S., Masihi M. Free fall and controlled gravity drainage processes in fractured porous media: Laboratory and modelling investigation // Can. J. Chem. Eng. 2015. V. 93. P. 2286.

  21. Allaire G. Homogenization and two-scale convergence // SIAM J. Math. Anal. 1992. V. 28. P. 1482–1518.

  22. Meirmanov A.M. Mathematical Models for Poroelastic Flows, Atlantis Studies in Differential Equations. V. 1. Beijing: Atlantis Press, 2014. 308 p.

  23. Мейрманов А.М. Повторное усреднение в задачах фильтрации подземных жидкостей // Научн. зап. Белгор. гос. ун-та. Сер. Математика–Физика. № 17 (136). Вып. 28. 2012. 16 с.

  24. Николаевский В.Н. Механика пористых и трещиноватых сред. М.: Недра, 1984. 233 с.

  25. Scheidegger A.E. The Physics of Flow through Porous Media. 3rd Edition. Toronto: Univ. Press, 1974. 372 p.

Дополнительные материалы отсутствуют.