Журнал вычислительной математики и математической физики, 2019, T. 59, № 11, стр. 1961-1972

ВВЕДЕНИЕ

Движение водных растворов в подземных пластах является важным фактором, влияющим на добычу полезных ископаемых, эффективное природопользование, эксплуатацию объектов альтернативной энергетики. Задачи фильтрации возникают при вытеснении нефти водой, глубинном захоронении углекислоты, извлечении горячих флюидов, транспортирующих на поверхность Земли геотермальную энергию. Во многих случаях численное моделирование является наиболее эффективным методом исследования, позволяющим изучать гидродинамические процессы и оптимизировать связанные с ними технологии. Для моделирования фильтрации многофазных многокомпонентных флюидов с учетом различных дополнительных факторов (фазовые переходы, энерговыделение в химических реакциях, анизотропия твердой матрицы и т.д.) применяются различные вычислительные коды, например, TOUGH2 [1], [2] и MUFITS [3], охватывающие широкий спектр явлений и открытые для некоммерческого коллективного пользования. Применяются также коды, разработанные авторами для решения конкретных задач, например, представленные в [4], [5].

Можно выделить класс задач, в которых основной движущей силой является сила тяжести – она инициирует во флюиде переменной плотности развитие естественной конвекции. Если же плотность флюида зависит, главным образом, от количества растворенной примеси, то возникает концентрационно-конвективное течение. Рассмотрение динамики подземных вод сводится к задаче о концентрационной конвекции в пористой среде, если образуется более плотный (неустойчивый) слой жидкости при растворении углекислоты в воде или нефти, при формировании концентрированных растворов солей, при испарении грунтовых вод [6]–[11]. Многие авторы используют собственные вычислительные коды, основанные на конечно-разностном методе [8]–[11], методе контрольного объема [7], спектральном методе – в [6] при решении уравнения Пуассона для функции давления выполняется преобразование Фурье. Преимуществами таких кодов являются простота использования, гибкость, возможность учесть основные факторы, пренебрегая второстепенными, а следовательно, и экономичность.

Одна из основных моделей, описывающих концентрационно-конвективные течения в пористых средах, включает уравнения неразрывности, движения (в форме уравнения Дарси) и переноса примеси [12], [13]. Особое внимание уделяется аппроксимации и методу численного интегрирования уравнения переноса примеси – нестационарного уравнения конвекции-диффузии. Простейшие разностные схемы, аппроксимирующие конвективный член, имеют существенные недостатки при использовании не слишком подробных пространственных сеток. Центрально-разностная схема, обеспечивающая второй порядок точности, может генерировать нефизические осцилляции переменной, в то время как схема против потока, имеющая первый порядок точности, может вносить существенную схемную вязкость и сглаживать решение в области резкого изменения переменной [14]. Чтобы обеспечить высокую точность численного моделирования в сочетании с экономичностью, эффективные разностные схемы для уравнения конвекции-диффузии активно разрабатываются и в настоящее время [15]–[17].

В данной работе описывается вычислительный код, основанный на конечно-разностном методе, который создан автором для моделирования концентрационно-конвективных течений в пористых средах. Фильтрационное течение и массоперенос описываются в рамках гидродинамической модели, конвективный член в уравнении конвекции-диффузии аппроксимируется с помощью схемы Quadratic Upstream Interpolation for Convective Kinematics (QUICK). На тестовой задаче о движении концентрационной ступени выполняется сравнение схемы QUICK с центрально-разностной схемой и схемой против потока. Возможности численного кода демонстрируются при решении задачи о начале и развитии концентрационной конвекции от источника примеси на границе в полубесконечной однородной области и в двухслойной области.

1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Движение жидкости с растворенной примесью внутри пористого образования в поле силы тяжести можно описать в рамках математической модели, которая включает уравнения неразрывности, Дарси и переноса примеси. Модель дополняется линейным уравнением состояния, связывающим плотность раствора с плотностью растворенной примеси:

Здесь ${\mathbf{U}}$, $P$, $\rho $, ${{\rho }_{c}}$ – скорость фильтрации, давление, плотности раствора и растворенной примеси; $k$, $\phi $, $\mu $, $D$ – проницаемость и пористость твердой фазы, вязкость и коэффициент дифузии раствора (в общем случае считаются переменными); $g$, ${\mathbf{e}}$, ${{\rho }_{0}}$ – модуль ускорения силы тяжести, единичный вектор, сонаправленный с вектором силы тяжести, плотность чистой жидкости; константа $\alpha $ определяет изменение плотности раствора при добавлении примеси.

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

Далее вместо плотности $\rho $ и давления $P$ в качестве переменных будем рассматривать отклонения этих величин от значений в чистой неподвижной жидкости: $\rho - {{\rho }_{0}}$ и $P - {{P}_{0}}$ соответственно. Для этого из уравнения (1.2) вычтем (1.4). Перейдем к безразмерным переменным, используя в качестве масштабов некоторую характерную длину $H$, диффузионную скорость ${{D}_{0}}{\text{/}}H$, время ${{H}^{2}}{\text{/}}{{D}_{0}}$, плотность ${{\rho }_{{{\text{sat}}}}} - {{\rho }_{0}}$, давление $({{\rho }_{{{\text{sat}}}}} - {{\rho }_{0}})gH$ и некоторые референсные значения проницаемости ${{k}_{0}}$, вязкости ${{\mu }_{0}}$ и коэффициента диффузии ${{D}_{0}}$. Индекс ${\text{sat}}$ соответствует насыщенному раствору. Уравнения (1.1)–(1.3) преобразуются в следующую безразмерную систему:

Здесь ${\mathbf{u}}$, $\Pi = (P - {{P}_{0}}){\text{/}}(({{\rho }_{{{\text{sat}}}}} - {{\rho }_{0}})gH)$, $S = (\rho - {{\rho }_{0}}){\text{/}}({{\rho }_{{{\text{sat}}}}} - {{\rho }_{0}})$ – безразмерные скорость фильтрации, давление, плотность; ${{C}_{k}} = k{\text{/}}{{k}_{0}}$, ${{C}_{\mu }} = \mu {\text{/}}{{\mu }_{0}}$, ${{C}_{D}} = D{\text{/}}{{D}_{0}}$ – переменные, задающие физические свойства твердой и жидкой сред. Скорость фильтрации ${\mathbf{u}}$ связана со скоростью движения в порах ${\mathbf{v}}$ соотношением ${\mathbf{u}} = \phi {\mathbf{v}}$. Наряду с безразмерной плотностью $S$ для представления результатов используется концентрация $c = {{\rho }_{c}}{\text{/}}\rho $ – физически более понятная величина, которая выражается через $S$ следующим образом: Константа ${{c}_{{{\text{sat}}}}}$ – это концентрация насыщенного раствора.

Система уравнений (1.5)–(1.7) содержит единственный критерий подобия – число Рэлея–Дарси ${\text{Rd}}$:

2. ЧИСЛЕННЫЙ МЕТОД

Уравнения (1.5)–(1.7) записываются в плоской декартовой системе координат. Интегрирование производится с использованием конечно-разностного метода на разнесенной неравномерной пространственной сетке: границы основной сетки совпадают с физическими границами области. Скалярные величины определяются в центре ячейки основной сетки, компоненты скорости – на ее гранях. Применяется двухслойная по времени аппроксимация.

На первом этапе совместно интегрируются уравнения неразрывности и Дарси – рассчитываются компоненты скорости и давления по алгоритму типа SIMPLE [18]. При корректировке для приращения давления получается уравнение Пуассона, которое в дискретном виде дает трехдиагональную по каждому направлению систему линейных алгебраических уравнений. Решение линейных уравнений осуществляется методом прогонки [14].

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

На третьем этапе выполняется совместная корректировка скорости, давления и плотности итерационным методом Якоби. Итерации выполняются до достижения условия $\left\| {\delta {{S}^{k}}} \right\| \leqslant \varepsilon $, где норма $\left\| {\delta {{S}^{k}}} \right\|$ представляет собой максимальное значение $({{S}^{k}} - {{S}^{{k - 1}}}){\text{/}}{{S}^{{k - 1}}}$, рассчитанное по всем узлам сетки; $k$ – номер итерации. В стандартном варианте $\varepsilon = {{10}^{{ - 6}}}$ и необходимое число итераций составляет несколько единиц.

Исходный вычислительный код, детально описанный в [19], включал центрально-разностную аппроксимацию конвективного члена и успешно применялся в течение ряда лет [9], [11], [20–22]. Чтобы центрально-разностная схема была монотонной, должно выполняться условие ${\text{Pe}} < 2$, где ${\text{Pe}} = V\delta x{\text{/}}D$ – сеточное число Пекле, построенное по скорости движения в порах $V$ и шагу сетки $\delta x$ [15]. Здесь $V$, $\delta x$, $D$ – размерные величины. Такое условие накладывает ограничение на шаг сетки, которая должна быть достаточно подробной в случае интенсивных течений. Для преодоления указанного ограничения в настоящей работе центрально-разностная схема заменена на схему QUICK.

Схема QUICK представляет собой трехточечную квадратичную схему, впервые предложенную в [23]; она сочетает схему против потока с градиентной поправкой и центрально-разностную схему. Фиг. 1 поясняет применение QUICK для одномерного потока; для простоты рассматривается равномерная сетка. Определим значение переменной ${{S}_{r}}$ на правой грани ячейки. Выбор узлов сетки зависит от направления движения. Если поток через грань положительный: ${{u}_{{i + 1/2}}} > 0$, то интерполяция производится по ${{S}_{{i - 1}}}$, ${{S}_{i}}$, ${{S}_{{i + 1}}}$, если ${{u}_{{i + 1/2}}} < 0$ – по ${{S}_{i}}$, ${{S}_{{i + 1}}}$, ${{S}_{{i + 2}}}$. На параболе, проведенной через точки ${{S}_{{i - 1}}}$, ${{S}_{i}}$, ${{S}_{{i + 1}}}$, переменная ${{S}_{r}}$ принимает значение

Правую часть уравнения (2.1) можно преобразовать в сумму двух слагаемых с весами где первое слагаемое ${{S}_{{r1}}}$ получено интерполяцией, включающей схему против потока, а второе ${{S}_{{r2}}}$ – интерполяцией с использованием центрально-разностной схемы: Запись в виде (2.2) наглядно демонстрирует смешанный тип схемы.

На левой грани ячейки значение переменной ${{S}_{l}}$ определяется аналогично. Конвективный член $\partial uS{\text{/}}\partial x$ аппроксимируется следующим образом:

Схема QUICK имеет третий порядок точности на равномерной сетке, она условно-устойчивая [24] и успешно работает при значениях сеточного числа Пекле ${\text{Pe}}$ порядка десятка. При больших ${\text{Pe}}$ схема становится немонотонной. Однако в задачах о естественной конвекции выбор пространственного шага диктуется не только ограничением на ${\text{Pe}}$, но и характерным масштабом образующейся конвективной структуры. При выборе сетки так, чтобы разрешить пространственную структуру, значение ${\text{Pe}}$, как правило, оказывается умеренным, поскольку скорости естественно-конвективных движений невысокие. Схема QUICK достаточно эффективна для моделирования конвективных течений; она применялась, в частности в [25], [26].

3. ТЕСТ: ДВИЖЕНИЕ КОНЦЕНТРАЦИОННОЙ СТУПЕНИ

В качестве теста рассматривается задача о движении концентрационной ступени при отсутствии силы тяжести. На плоскости задана прямоугольная расчетная область $\Omega $:

Скорость фильтрации имеет компоненты ${\mathbf{u}} = ({{u}_{x}},{{u}_{y}})$, скорость движения жидкости в порах ${\mathbf{v}} = ({{v}_{x}},{{v}_{y}})$. В начальный момент область заполнена чистой водой, которая совершает фильтрационное движение с постоянной горизонтальной скоростью $u_{x}^{{{\text{in}}}} > 0$. Начальные значения отмечены верхним индексом ${\text{in}}$. Давление линейно падает слева направо, что следует из уравнения Дарси (1.6). Начальные условия следующие: На левой границе поддерживается постоянная концентрация примеси, так что примесь начинает поступать в область вместе с потоком воды, образуя концентрационную ступень, которая отделяет чистую воду от раствора. Концентрационная ступень перемещается с потоком слева направо и сглаживается за счет диффузии. Правая и горизонтальные границы непроницаемы для примеси. Через горизонтальные границы вода не протекает, скорости потока на левой и правой границах совпадают. Граничные условия имеют вид

Физические свойства твердой и жидкой фаз считаются постоянными: ${{C}_{k}} = {{C}_{\mu }} = {{C}_{D}} = 1$. Задаются концентрация на входе и пористость: $c* = 0.1$, $\phi = 0.2$. Плотность раствора на входе $S{\text{*}}$ можно рассчитать по $c{\text{*}}$, выразив $S(c)$ из (1.8). Получается $S* = 0.230$ при характерных для воды константах $\alpha = 0.815$ и ${{c}_{{{\text{sat}}}}} = 0.342$. Взято число Рэлея–Дарси ${\text{Rd}} = {{10}^{{ - 5}}}$, что соответствует крайне слабой силе тяжести (на несколько порядков меньше силы тяжести Земли), влиянием которой на движение ступени можно пренебречь.

Безразмерная высота расчетной области составляет ${{h}_{y}} = 0.2$. Область разбивается на одинаковые квадратные ячейки со стороной $\Delta x = \Delta y = 0.02$, т.е. покрывается равномерной сеткой 50 × 10. Вершины сдвинутой сетки лежат в центрах ячеек. Выбраны два значения скорости фильтрации: $u_{x}^{{{\text{in}}}}$ = 20, 200, которым соответствуют значения скорости движения $v_{x}^{{{\text{in}}}} = {{10}^{2}}$, ${{10}^{3}}$. Сеточное число Пекле ${\text{Pe}}$ рассчитывается по безразмерным параметрам как ${\text{Pe}} = v_{x}^{{{\text{in}}}}\Delta x$ и принимает значения ${\text{Pe}} = 2$, $20$. Берется шаг интегрирования по времени $\tau = 2 \times {{10}^{{ - 6}}}$, $2 \times {{10}^{{ - 7}}}$ в первом и втором случаях. Решается начально-краевая задача (1.5)–(1.7), (3.2)–(3.5) в области (3.1).

На фиг. 2 приведены распределения концентрации $c$ в моменты времени $t = 4 \times {{10}^{{ - 4}}}$ (a), $t = 4 \times {{10}^{{ - 3}}}$ (б), когда концентрационная ступень (область резкого изменения $c$) достигла координаты $x = 0.4$. Задача близка к одномерной – решение не зависит от координаты $y$. Для сравнения показаны кривые $c(x)$, полученные при аппроксимации конвективного члена по центрально-разностной или противопоточной схемам. Схема QUICK переходит в центрально-разностную, если при определении значения ${{S}_{r}}$ на грани ячейки положить $\beta = 0$ в (2.2). Чтобы получить схему против потока, нужно взять $\beta = 1$. Также дана аналитическая зависимость $c(x)$, найденная из решения уравнения диффузионного массопереноса $\partial S{\text{/}}\partial t = \Delta S$. К последнему сводится уравнение конвекции–диффузии (1.7) в неподвижной системе с постоянными свойствами. При начальных условиях $S = S*,$ $0 \leqslant x \leqslant 0.4$ и $S = 0,$ $0.4 < x \leqslant 1.0$ решение одномерной задачи имеет вид [27]

Функция $S(x,t)$, взятая при двух выбранных $t$, преобразована в $c(x)$ по (1.8).

Расчеты показывают, что при сеточном числе Пекле ${\text{Pe}} = 2$ (фиг. 2а) все три схемы дают достаточно точные решения. При ${\text{Pe}}$ (фиг. 2б) схема QUICK моделирует движение концентрационной ступени с высокой точностью – форма профиля концентрации практически совпадает с аналитическим профилем. Центрально-разностная схема приводит к осцилляциям за ступенью, т.к. возмущение распространяется вверх по потоку из-за невыполнения свойства транспортивности; заметно некоторое сглаживание профиля вниз по потоку. Противопоточная схема дает осцилляции вниз по потоку и сглаживание профиля за ступенью. Таким образом, тест наглядно демонстрирует преимущества схемы QUICK.

4. НАЧАЛО И РАЗВИТИЕ КОНВЕКЦИИ В ПОЛУБЕСКОНЕЧНОЙ ОБЛАСТИ

Численный код применен для решения задачи о начале и развитии концентрационной конвекции в поле силы тяжести. Имеется полубесконечная однородная пористая среда, ограниченная сверху горизонтальной границей. Исходно твердая матрица заполнена чистой неподвижной водой. На границе поддерживается постоянная концентрация примеси, равная концентрации насыщения ${{c}_{{{\text{sat}}}}}$ (это равнозначно $S = 1$), что обозначает источник массы. Это может быть, например, пласт кристаллической соли. Примесь растворяется в воде и распространяется вниз за счет диффузии. Под границей образуется слой раствора, плотность которого больше, чем плотность чистой воды. Под действием силы тяжести слой стремится опуститься, что приводит к развитию естественной концентрационной конвекции.

Задается прямоугольная расчетная область $\Omega $, верхняя сторона которой лежит на границе пористой среды. Геометрические размеры нормируются на высоту области $H$.

Считается, что все границы, кроме верхней, примесь не пропускают. Постановка задачи показана на фиг. 3. Начальные условия имеют вид Граничные условия следующие:

Свойства скелета и воды полагаются постоянными: ${{C}_{k}} = {{C}_{\mu }} = {{C}_{D}} = 1$. Используются параметры: ${{\rho }_{0}} = 843$ кг/м$^{3}$, ${{\rho }_{{{\text{sat}}}}} = 1168$ кг/м$^{3}$, ${{\mu }_{0}} = 1.18 \times {{10}^{{ - 4}}}$ Па с, ${{D}_{0}} = 9 \times {{10}^{{ - 10}}}$ м$^{2}$/с, $g = 9.81$ м/с$^{2}$, $\phi = 0.1$, ${{h}_{x}} = 2$. В реальных условиях проницаемость горных пород варьируется на несколько порядков и характеризуется значениями ${{k}_{0}} = {{10}^{{ - 17}}}{\kern 1pt} - {\kern 1pt} {{10}^{{ - 13}}}$ м$^{2}$. При расчете числа Рэлея–Дарси ${\text{Rd}}$ (1.9) следует учесть, что выбор высоты области $H$ произволен, т.к. система не имеет характерного геометрического масштаба, который определялся бы условиями задачи. Кроме того, параметры $H$ и ${{k}_{0}}$ входят в ${\text{Rd}}$ как произведение, поэтому значение ${\text{Rd}}$ рассчитывается по комбинации $H{{k}_{0}}$. Для определенности возьмем $H = 20$ м и ${{k}_{0}} = {{10}^{{ - 15}}}$ м$^{2}$и получим ${\text{Rd}} = 600$. Равное значение ${\text{Rd}}$ и такое же решение в безразмерных переменных получаются, например, в парах $H = 2$ м и ${{k}_{0}} = {{10}^{{ - 14}}}$ м$^{2}$, $H = 50$ м и ${{k}_{0}} = 4 \times {{10}^{{ - 16}}}$ м$^{2}$ и т.д.

Исследования показывают, что распространение примеси внутри области проходит в три стадии: 1) диффузионный режим – формируется диффузионный пограничный слой под границей области, профиль примеси устойчив, конвекции нет; 2) линейное развитие конвекции (движение и массоперенос можно описать уравнениями, линеаризованными относительно малых возмущений) – профиль примеси теряет устойчивость, в пограничном слое возникают периодически растущие возмущения переменных и начинается конвекция; 3) конвекция переходит в нелинейный режим. Задача на этапах 1) и 2) допускает аналитическое решение. В частности, в [6] проведен линейный анализ устойчивости и найдено критическое значение длины волны возмущения ${{\lambda }_{c}} = 2\pi \phi {{\mu }_{0}}{{D}_{0}}{\text{/}}(0.07\Delta \rho gH{{k}_{0}})$; при использовании числа Рэлея–Дарси (1.9) получается

В [6] при сравнении аналитического и численного решений продемонстрировано, что начальный этап развития конвекции описывается линейным приближением с хорошей точностью. На этапе 3) решение может быть получено только численно.

Численное моделирование проводилось на пространственной сетке 1600 × 800, равномерной по $x$ и сгущающейся по $y$ около границы-источника. Шаг сетки по $x$ составлял $\Delta x = 1.25 \times {{10}^{{ - 3}}}$, по $y$ менялся от $\Delta {{y}_{{{\text{min}}}}} = 2.61 \times {{10}^{{ - 4}}}$ до $\Delta {{y}_{{{\text{max}}}}} = 1.31 \times {{10}^{{ - 3}}}$, шаг интегрирования по времени – $\tau = {{10}^{{ - 8}}}$.

На фиг. 4 дано поле концентрации примеси в различные моменты времени; белым цветом показана чистая вода, самым темным – насыщенный раствор. Картина на фиг. 4а соответствует диффузионному режиму – примесь скапливается в узком невозмущенном пограничном слое. Время начала конвекции определяется по скорости движения жидкости ${\mathbf{v}}$. Находится скорость, максимальная по области (она обозначена как ${{\left| {\mathbf{v}} \right|}_{{{\text{max}}}}}$), и фиксируется момент начала роста ${{\left| {\mathbf{v}} \right|}_{{{\text{max}}}}}$. Получено, что конвекция начинает развиваться с момента $t \approx 8.86 \times {{10}^{{ - 6}}}$. Возмущения пограничного слоя на начальной стадии являются периодическими с длиной волны ${{\lambda }_{c}}$. Чтобы найти ${{\lambda }_{c}}$ в численном моделировании, анализируется вертикальная компонента скорости ${{v}_{y}}(x)$ в горизонтальном сечении $y = - 0.02$ непосредственно под границей-источником. На фиг. 5 представлена часть кривой ${{v}_{y}}(x)$ при $x \in [0,0.04]$ после возникновения движения. Первые растущие возмущения ${{v}_{y}}(x)$ наблюдаются около боковых границ, поэтому ${{\lambda }_{c}}$ находится по максимуму и минимуму (между которыми содержится полволны) около левой границы, они отмечены маркерами. Получается ${{\lambda }_{c}} = 0.018$, что близко к теоретической оценке ${{\lambda }_{c}} = 0.015$, следующей из (4.1). Длина волны в моделировании вмещает ${{\lambda }_{c}}{\text{/}}\Delta x = 0.018{\text{/}}0.00125 \approx 14$ узлов сетки. Со временем длина волны возмущения увеличивается.

Возмущения на начальной (линейной) стадии конвекции могут расти в произвольных местах пограничного слоя. Картина на фиг. 4б близка к линейной стадии – в различных местах формируются периодически расположенные солевые пальцы. Далее солевые пальцы удлиняются, стремясь вниз, они искривляются, объединяются и становятся толще, превращаясь в солевые языки, скорость движения которых, толщина и глубина проникновения произвольны. Картина течения становится стохастической – конвекция переходит в нелинейную стадию, показанную на фиг. 4в–д. Интенсивность течения в области в целом можно охарактеризовать максимальной скоростью ${{\left| {\mathbf{v}} \right|}_{{{\text{max}}}}}$, график которой дан на фиг. 6 (сплошная линия). Видно, что в развитом течении достигаются значения ${{\left| {\mathbf{v}} \right|}_{{{\text{max}}}}} \approx 3 \times {{10}^{3}}$, т.е. конвективный перенос на три порядка превосходит диффузионный перенос примеси. Сеточное число Пекле ${\text{Pe}}$ не превышает величину ${\text{Pe}} = {{\left| {\mathbf{v}} \right|}_{{{\text{max}}}}}\Delta {{y}_{{{\text{max}}}}} \approx 4.0$.

Решена задача о начале и развитии конвекции в двухслойной области. В отличие от рассмотренной выше задачи твердая порода состоит из двух частей, верхняя характеризуется прежними значениями пористости $\phi $ и проницаемости ${{k}_{0}}$, а в нижней пористость в два раза больше. Величина проницаемости, которая связана с $\phi $ зависимостью $k \sim {{\phi }^{3}}{\text{/}}{{(1 - \phi )}^{2}}$ [12], в нижней части оказывается равной $k = 10{{k}_{0}}$. Нижняя высокопроницаемая подобласть начинается на уровне $y = - 0.3$. Таким образом, $\phi = 0.1$, ${{C}_{k}} = 1$ при $ - 0.3 < y \leqslant 0$ и $\phi = 0.2$, ${{C}_{k}} = 10$ при $y \leqslant - 0.3$; остальные параметры не меняются.

В этой задаче развитие конвекции вначале происходит так же, как и в однородной области, рассмотренной выше. До тех пор, пока конвективное движение не достигнет нижней подобласти, решения обеих задач будут качественно одинаковыми. На фиг. 7 дано поле концентрации примеси, когда солевые языки уже дошли до высокопроницаемой зоны и наблюдается нелинейный режим движения. Моменты времени на фиг. 7а, б совпадают с моментами на фиг. 4г, д соответственно. Можно заметить, что языки, попадая в высокопроницаемую подобласть, становятся узкими, поскольку начинают двигаться в несколько раз быстрее. Это подтверждается данными на фиг. 6, показывающими возрастание скорости ${{\left| {\mathbf{v}} \right|}_{{{\text{max}}}}}$ в 2–4 раза, начиная с момента $t \approx 6.50 \times {{10}^{{ - 4}}}$ (пунктирная линия). Кроме того, кривая скорости после этого момента становится пилообразной – это указывает на то, что язык, вошедший в зону повышенной проницаемости, сначала сильно ускоряется, затем, достигнув максимальной скорости, начинает замедляться. При этом кривая ${{\left| {\mathbf{v}} \right|}_{{{\text{max}}}}}$ идет сначала вверх, затем вниз до тех пор, пока другой язык не войдет в высокопроницаемую зону. Толщина солевых языков в самом узком месте составляет величину порядка $0.05$ и включает около 40 узлов пространственной сетки. Как видно на фиг. 6, скорость движения не превосходит значения ${{\left| {\mathbf{v}} \right|}_{{{\text{max}}}}} \approx 1.4 \times {{10}^{4}}$, следовательно, сеточное число Пекле ограничивается величиной ${\text{Pe}} \approx 17.5$.

ЗАКЛЮЧЕНИЕ

Разработан вычислительный код для моделирования концентрационно-конвективных течений в пористых средах на основе гидродинамической модели, включающей уравнения неразрывности, Дарси и переноса примеси с переменными свойствами твердой и жидкой фаз. Дискретизация уравнений осуществляется конечно-разностным методом на разнесенной неравномерной сетке. Уравнения неразрывности и Дарси интегрируются совместно с помощью алгоритма типа SIMPLE. Затем решается уравнение переноса примеси (нестационарное уравнение конвекции-диффузии), в котором конвективный член аппроксимируется по схеме QUICK. Тестовые расчеты, проведенные на задаче о движении концентрационной ступени, показали преимущества схемы QUICK по сравнению с центрально-разностной и противопоточной аппроксимациями – получено, что область резкого изменения переменной воспроизводится с хорошей точностью при более высоких значениях сеточного числа Пекле, а значит, можно использовать более грубые сетки.

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

Автор благодарит О.А. Бессонова и Г.Г. Цыпкина за полезные обсуждения.