Возникновение концентрационной конвекции в длинной прямоугольной области пористой среды Текст научной статьи по специальности «Физика»

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

2023

• ФИЗИКА •

Вып. 4

УДК 532.546; 532.5.013.4 PACS 47.56.+r, 47.55.pd

Возникновение концентрационной

конвекции в длинной прямоугольной

пористой среды

М. Р. Хабин121, Б. С. Марышев12 *

1 Институт механики сплошных сред УрО РАН, Пермь, Россия

2 Пермский государственный национальный исследовательский университет, Пермь, Россия

tmikhail.khabin@psu.ru

ibmaryshev@mail.ru

Решается задача конвективной устойчивости течения в вытянутой прямоугольной области пористой среды. Транспорт растворенного вещества в жидкости описывается с помощью MIM (mobile-immobile media) подхода. Данный подход заключается в разделении концентрации растворенного вещества на подвижную и неподвижную компоненты. Влияние неоднородности плотности на течение смеси в пористой среде учитывается в уравнении для скорости фильтрации в приближении Дарси-Буссинеска. Задача решалась численно с использованием метода конечных разностей. Получены поля давления и концентрации. Проанализированы влияния параметров задачи на профили возмущений. Решена модифицированная задача устойчивости. Необходимость модификации возникает, поскольку требуется находить колебательные режимы конвекции регулярным образом. Для этого производился расчет спектра показателей Ляпунова с ортогонализацией Грама-Шмидта. Найдены нейтральные кривые, позволяющие обнаружить порог возникновения конвекции. Получены области параметров, в которых реализуются колебательные возмущения. Проанализировано влияние параметров сорбции на возникновение колебаний.

Ключевые слова: пористая среда; концентрационная конвекция; устойчивость течения

Поступила в редакцию 13.09.2023; после рецензии 09.10.2023; принята к опубликованию 23.10.2023

The onset of concentration convection in a long rectangular domain of a porous medium

M. R. Khabin121, B. S. Maryshev12*

1 Perm State University, Perm, Russia

2 Institute of Continuous Media Mechanics UB RAS, Perm, Russia tmikhail.khabin@psu.ru

ibmaryshev@mail.ru

We solve a problem of the convective flow stability in an elongated rectangular region of a porous medium. The transport of a mixture is described using the MIM (mobile-immobile media) approach. This approach consists in dividing the concentration of the solute into mobile and immobile components. The effect of density inhomogeneity on the flow of a solution in a porous medium is taken into account by the filtration equation within the Darcy-Boussinesq approximation. The problem is solved numerically with the use of the finite difference method. The fields of pressure and concentration are obtained. The influence of the parameters of the problem on the disturbance profiles is analyzed. A modified problem of stability is solved. The modification is necessary for finding the oscillatory convective mode. To obtain the solution, the spectrum of Lyapunov exponents is calculated with Gram-Schmidt orthogonalization. Neutral curves are found that make it possible to detect

© Хабин М. Р., Марышев Б. С., 2023

распространяется на условиях лицензии

Creative Commons Attribution 4.0 International (CC BY 4.0).

the threshold for the onset of convection. The ranges of parameters in which oscillatory perturbations are realized are obtained. The effect of sorption parameters on the occurrence of oscillations is analyzed.

Keywords: porous media; concentration convection; stability of flow

Received 13 September 2023; revised 09 October 2023; accepted 23 October 2023

doi: 10.17072/1994-3598-2023-4-10-21

1. Введение 1.1. Актуальность работы

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

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

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

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

Исходя из этого, изучение адсорбционно-десорбционных свойств необходимо при описании большого числа технологических процессов, связанных с фильтрацией, а также сохранения окружающей среды. Стандартная схема эксперимента по изучению адсорбционно-десорбционных свойств пористой среды - прокачка заданного объема смеси через рабочую область при заданном перепаде давления между ее границами. В ряде таких экспериментов [2] наблюдалось влияние гравитационной неустойчивости на ход экспериментов. Исследование концентрационной конвекции, возникающей вследствие этой неустойчивости, необходимо для корректного учёта её влияния при подготовке и проведении экспериментов по изучению процессов сорбции в пористых средах.

1.2. MIM подход и модели

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

Наиболее популярным подходом к описанию процессов сорбции в пористых средах является MIM (mobile-immobile media) подход [3], заключающийся в разделении примеси на две фазы: свободную и связанную (адсорбированную). Свободная - может перемещается вместе с потоком несущей жидкости, а связанная взаимодействует (оседает) с твердым скелетом пористой среды. Связанную примесь называют иммобилизованной,

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

- В стандартной модели MIM [4] поток примеси в иммобильную фазу линейно растет с увеличением концентрации примеси в мобильной фазе и линейно уменьшается пропорционально возрастанию концентрации в иммобильной фазе, такую модель также называют линейной;

- Фрактальная (дробная) модель MIM [5] использует предположение, что частицы примеси в пористой среде могут оставаться связанными на случайные сколь угодно длинные промежутки времени, распределение такой случайной величины не имеет среднего значения (что является признаком субдиффузионного переноса). В таком случае кинетическое уравнение выражается линейной зависимостью между потоком частиц в иммобильной фазе и дробной производной Капуто от концентрации частиц в мобильной фазе по времени. Данная модель имеет экспериментальное подтверждение, в основном для малых концентраций примеси [6];

- Нелинейная MIM модель с насыщением [7] описывает транспорт в пористых средах в присутствии значительного количества примеси. Модель предполагает наличие некоторого предельного значения концентрации иммобилизованной примеси, при котором процесс адсорбции прекращается. В [7] данный эффект описывается в рамках кинетической модели второго порядка, в которой скорость адсорбции прямо пропорциональна разнице между предельной концентрацией насыщения и концентрацией частиц в неподвижной фазе.

В настоящей работе рассматривается двумерная задача устойчивости течения смеси через прямоугольную область пористой среды. Предполагается, что вертикальный размер области меньше горизонтального в 10 раз. Между вертикальными границами задаётся перепад давления и концентрации, на горизонтальных границах - условие отсутствия потока несущей жидкости и примеси. Задача решается в приближении Дарси-Буссинеска [8], а транспорт примеси описывается в рамках линейной MIM модели [4]. Задача устойчивости исследуется с помощью расчета спектра показателей Ляпунова [9] с ортогонализацией Грама-Шмидта [10].

Статья содержит восемь параграфов. В первом параграфе описаны мотивация решения задачи и

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

2. Постановка задачи

Рис. 1. Принципиальная схема двумерного транспорта примеси в рабочей области

Рассматривается двумерная задача (рис. 1) об однородном горизонтальном просачивании смеси через замкнутую область пористой среды в поле силы тяжести. В прямоугольной области длиной L и высотой H, заполненной пористой средой, на вертикальных стенках задается постоянный перепад концентрации AC = C+ - C- и давления Ap = P+ - P-.

Моделирование транспорта примеси в пористой среде производится в рамках MIM подхода [3], предполагающего разделение концентрации примеси на мобильную и немобильную фазы. Это разделение позволяет учесть тот факт, что примесь может оседать на стенках пористой среды. Таким образом, уравнение транспорта примеси описывается с учетом несжимаемости несущей жидкости следующей системой уравнений, вывод которой приводится в [11]:

д(фе + q)

dt

div u = 0.

-u -Vc + $DAc,

(2.1)

где е,д - объемные концентрации мобильной и немобильной компонент, соответственно, и - скорость фильтрации, ф - пористость среды, В -эффективный коэффициент диффузии.

Для связи скорости фильтрации с приложенным внешним давлением в присутствии неоднородности плотности, вызванной неоднородностью концентрации, воспользуемся законом Рэлея- Дар-си [8]:

к

u =--(Vp + Pßcg).

1

(2.2)

где и - скорость фильтрации, к - коэффициент проницаемости, п - динамическая вязкость, р -

плотность, ß - коэффициент концентрационного расширения, c - концентрация, g - ускорение свободного падения.

Для описания процессов сорбции, т.е. перехода примеси между мобильной и иммобильной фазами, в рамках MIM подхода вводится кинетическое уравнение. Для линейной MIM [4] модели кинетическое уравнение имеет вид

d( c + q )

|q = а(c -Kdq),

(2.3)

(2.4)

dt

= -u • Vc + Ac,

dq = ac - bq, dt

u = - PeVp + Rpck, div u = 0;

(3.1)

где а - коэффициент переноса примеси, К -коэффициент распределения примеси, с,д - концентрации мобильной и немобильной компонент, соответственно.

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

д(рс + а )

-^ = и - Ус + рDДc,

где Ре=к(р+-р-)/цБф - число Пекле, Rp=pgкP(C+-C-УцDф - число Рэлея-Дарси, а=аЬ2/фБ - коэффициент адсорбции, Ь= аКй Ь2Ю - коэффициент десорбции. Число Пекле характеризует отношение сил давления к вязким силам, а число Рэлея-Дарси - отношение сил плавучести к силам вязкости. Коэффициенты сорбции отвечают за интенсивность перехода примеси между фазами.

Исключим скорость фильтрации из уравнений (3.1), используя условие несжимаемости:

д(с + а) _ д2с д2с

dt

dqq = а(c - Kdq),

и = -~(Vp + pßcg), div u = 0.

Граничные условия для данной задачи: постоянный перепад давления и концентрации между вертикальными границами, горизонтальные границы непроницаемы, что математически выражается в следующем виде:

х _ 0: р _ Р+, с _ С+,

х _ L : р _ Р_, с _ С _, дс

У _ 0: и _ 0, — _ 0, (2.5)

' у ду

дс

У _ Н : и _ 0, — _ 0.

у ду

3. Постановка двумерной задачи о переносе примеси

Удобно переписать уравнения (2.4) в безразмерной форме. Для этого в качестве единиц измерения выберем следующие величины:

- концентрация мобильной компоненты

[с] = С+ - С-;

- концентрация немобильной компоненты

Ы = Ф(С+ - С-);

- время [г] = L2/D;

- длина [1] = L;

- скорость фильтрации [и] = Т ;

- давление - [р] = р+ - р-.

Тогда система (2.4) перепишется в безразмерном виде:

dt

2 +—7 + dx2 dy2

„ I dc dp dc dp | „ dc

+ Pe I--- +--- I- Rpc—,

\ dx dx dy dy) dy

^ = ac - bq, dt

d2 p d2 p Rp dc

—TT + —---— = 0,

dx2 dy2 Pe dy

Граничные условия перепишутся в виде:

Rp |

(3.2)

dp dy

Pe iy=°.h;

y=0,h

p (0, y, t ) = 1, p (1, y, t) = 0,

dc = 0,

(3.3)

дУ у_0,И

с (0, у, г)_ 1, с (1, у, г)_ 0.

В отсутствие конвекции (Яр = 0) задача допускает следующее основное решение:

Ре-х Ре

е - е

c0 (x ) =

С Pe^x Pe Л

q (x )=a {^) •

p0 (x) = 1 - x. и0 = {Pe, 0};

(3.4)

4. Метод решения задачи

Задача (3.3)-(3.4) решается числено с использованием метода конечных разностей на базе схемы второго порядка точности по пространству и первого - по времени. Решение задачи находится в области х е[0, L], у е [0, Н] на сетке с квадратной ячейкой с шагом И = L/N, где N - число узлов по

горизонтали. Шаг по времени т=Т/К, следовательно, обозначим:

q(x,y,t) = q(ih,jh,kт) = qkijt

р(х,у, ^ = р(Ш,]И, кт) =/„■,

Xi=i■h , i=0,N, yj=j■h , ]=0,М, при этом ,=10М, где М - число узлов по вертикали. Сетка по временной координате к=0,К.

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

Рок, = 1,

р,, ■ = о,

Рим _ Р

Рко = Ра ■

hRp

Ре

^¡М-1'

hRp

--<

Ре

(4.1)

с0, ■ _1

гк

-0, ]

СN, ] = о,

к _ к

С/,1 _ С1,0, ск _ ск

С1,М _ С1,М

р0] _ 1 - i • h,

с0] _

q0 ] _

Ре^А Ре

е - е

1 - еР

(4.2)

, Ре • /• И Ре

а | е - е

1 - еР

Л ск+1/2 + В ск+1/2 + С Ск+1/2 _ Р

1+1,] т -"1°I,] т I-1,] 11

г Лт

А _--7--

1 2И2 2

В _ 1 + 4 - аЬт

+

ат

(

Р _

И2 2 (2 + Ьт) 2

С1 _-т + -1 2И2 2

^ „ Л

т+гв - крс]

2И2 2 4И

V у

< ]+1 +

(4.3)

Н' - ]

+ъ]1

2И2 2 4И

V у

ск]-1+

Ьт к +(2+Ьг)^

Л Ск+1 + В Ск+1 + С Ск+1 _ Р

Л2 _-т-Гв+^,

2 2И2 2 4И

В2 _1 + 72-

И2 2 (1 + Ьт) 2

С, _-т+Гв -

2И2 2

4И

2И2 2

р _ | т+— |ск,+1/2 + [ 1 -Г Iск+1/2 +

-¿+1, ]

+ 1 -Гг - —I с^1]2 + ,ЬГ , дк.;

-1, ] кк

2И2 2 У '-1,] 2 (1 + Ьт)

Л _ Ре (РЦ] - Р^]) В _ Ре (Р+и - р£и).

4И

атСк+1 + ак к+1 _ ] ]

qi,j

4И2

(1 + Ьт)

_ 4 (Р/+1, ] + р1-1] + Р ]+1 + Р ]-1 )-

(4.4)

м 4

1 Rp И

4Ре21 и+1 С'-]-1''

(4.5)

тах

- Р!

<

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

где е - задаваемый порядок точности, определяющий окончание расчёта давления.

Уравнения системы (4.3) решаются методом прогонки, на первом полушаге по времени находится ск+1/2^ , а на втором - ск+1^ . С помощью соотношения (4.4) вычисляется qk+1ij . Давление на

к+ 1

следующем шаге по времени ру определяется с помощью (4.5) по схеме Ричардсона [12], где 5 -номер итерации в схеме.

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

I _Л^(х,>-Др,Ре)|(1Б .

(4.6)

Численный расчет проводился для сетки размером 100^10. Шаг сетки И=0.01, шаг по времени составлял т=0.00001 (К=10000). Параметры сорбции для представленных расчетов а=1, Ь=2. Результаты расчета приведены для t=1.

5. Конвективная устойчивость течения в рамках модифицированной задачи

5.1. Постановка модифицированной задачи

Задача устойчивости решается в модифицированном виде. Модификация заключатся в том, что

2

И

к

к

¿л

5

Ь

2

Б

+

в отличие от задачи (3.2)-(3.4) рассматриваемое течение имеет вертикальную компоненту. Это означает изменение граничных условий (2.5) - горизонтальные границы становятся проницаемыми:

« {Ре,0} ^ и {РеДрс0}. Задача может быть записана в форме (3.2):

д( с+а)

(5.1)

дt

д с д c

-7 +-7 +

дх2 ду2

„ [ дс др дс др | „ дс

+ Pel--- +--- I-Rpc—,

^ дх дх дУ дУ) ду

^ = ас -bq,

at

д2 р д2 р Rp дс „

—тт + —т----— = 0,

дх2 дУ2 Pe ду

с измененными граничными условиями:

р| , = 1 - х,

(5.2)

р (0, у, t) = 1, р (1, у, t) = 0,

дс

дУ

= 0,

(5.3)

у_0,И

с (0, у, г)_ 1, с (1, у, г)_ 0.

Основное решение частично совпадает с решением (3.4) за исключением выражения для скорости, и может быть записано в форме:

Ре-х Ре

е - е

с ( х, у, t ) = с, ( х) + с ( х, у, t),

q ( х, У, t) = q0 ( х)+q ( х, у, t), р ( х У,t ) = р ( х) + р ( х У,t) •

(5.5)

Подставим (5.5) в (5.2), учитывая решение в основном состоянии (5.4), линеаризуем уравнения (5.2):

д(д+qq) Л~ P дс,

—-- = Ас - Pe---+

дt

дх

+Pe ^ -Rpc *

дх дх ду

= ас - b

дq ~dt

.„ Rp дс Ар = ^—, Pe ду

(5.6)

I у =0,h

= 0, р| = 0,

' 1х=0,1 '

дс?

ду

= 0, с = 0.

' 1х=0,1

у =0,h

Задача (5.6) решается методом Галеркина [13] (в модификации Галеркина-Канторовича) для следующих базисных функций:

: (t, х) = с0 (t, х) + с1 (t, х) cos —у +

+ с2 (t, х) cos

2—у;

q (t, х) = q0 (t, х) + ql (t, х) cos —у +

z ч 2—у

+ q2 (t, х) cos—; р (t, х) = р1 (t, х) sin у + р2 (t, х) sin .

(5.7)

с0 ( x, У )_ —

а ( еРе х — еРе ^ а0(х,у)_ ь ^"Г—^ ], (5.4)

р0 (x, у )_1 - х.

«0 _(Pe,Rpco).

5.2. Задача для возмущений основного состояния

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

В результате получается следующая система уравнений для коэффициентов:

, г 2И

И/0 +—й _ 0, -

И г 4И

~2 /1 + _

И 2И - /2 - —Й = 0, 2 3-

дс да д2 с

п , 1п п ,

f =_n +__in -

n дt дt дх'2

( — n Y „ дсп

+ сп ["Г J + Pe ~дх'

g = Pe2 .Pe-х дрп _ (1 - ePe) дх

п—Rp ( ePel - ePe

h I 1 - eP

^-(асп -bqn) = 0, n = 0,1,2, д2 рп ( —n V Rp —n

дх2 ( h

р+^Тсп=а n=1,2. (5.8)

h

Система (5.8) решается численно методом конечных разностей [11] с точностью второго порядка по пространству и первого порядка по времени.

5.3. Расчет спектра показателей Ляпунова с ортогонализацией Грама-Шмидта

Задача устойчивости состоит в определении значений параметров системы уравнений (5.8), при которых возмущения основного состояния остаются нейтральными. Для этого производится расчет спектра показателей Ляпунова по следующей процедуре [12]: рассмотрим вектор '(ов0, 0*1, ов2), где состоящий из значений концентраций получаемых в ходе решения задачи конечно-разностным методом. Пусть задача задаёт некоторое преобразование таких векторов со временем, обозначим факт этого преобразования как К('), очевидно, результатом будет некий вектор той же размерности, что и Выберем некоторую систему взаимоортогональных векторов (в данном случае мы ограничимся тремя, поскольку выбранный базис Галеркина содержит три функции), каждый из которых будем использовать как начальное условие для задачи. Тогда, численно решая задачу для каждого из выбранных векторов в течение времени t, получим:

К

К

лвХ

о{ = 0, о2 = 0

=

(в в в \ С0 , 01 , 02 ) ,

1 00 = 0,0 = sinо2 = 0

= ' 0

: (с0 ,01 ,02 ),

(5.9)

31 0 = 0, о* = 0, о2 = sin ^^

=

( в в в \ С0, 01, 02 ),

К [' 01 (о0,01в, 02)] = '11 (С0в, о, о2), К02 (о0, о, о2)] = ',2 (С0в, о, о2),

К03 (ов, о;, о2)] = (о0, о;, о2).

(5.10)

Заметим, что поскольку вектора '01, '02, 'оз ортогональны, то любой вектор '(ов0, ов1, ов2) может быть разложен в их базисе, тогда:

' (0)= 01 + В'02 + Cw03,

' (Т) = А'п + В'п + С'13.

(5.11)

Предположим, что для решения задачи имеет место отображение

(t = Т ) = и' ^ = 0)

(5.12)

где / - мультипликатор решения. Согласно [20] 1/1=1 отвечает нейтральным возмущениям (при большем значении возмущения нарастают, при меньшем затухают). Умножив скалярно выражение (5.12) на соответствующие вектора '01, '02, '03, получим систему алгебраических уравнений:

(('11- '01 )~м) А + ('12- '01 ) В +

(5.13)

+ ('13-'01 ) С = 0

('ц-'02 ) А + (('12-'02 )-м) В +

+ ('13-'02 ) С = 0 ('11-' 03 ) А +('12-'03 ) В + + (('13-'03 )-и) С =

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

О^3 + 02/и2 + 03/и + О4 = 0.

(5.14)

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

Для применения метода [12] необходимо иметь набор ортогональных начальных условий. Поэтому полученные вектора ортогонализируют методом Грама-Шмидта [10]. Применим к ортого-нализованным начальным векторам преобразование ещё раз, только в качестве временного интервала выберем Т:

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

(Rp,Pe, а, Ь ) = 1

(5.15)

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

= агс^ (1т и ^ и)

Т

(5.16)

в

I о0 = Sin

к

У1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

х

-4.0 -3.5 -2.9 -2.3 -1.8 -1.3 -0.7 -0.1 0.4 1.0 1.5 2.1 2.6 3.2

с* ( х, у )-10

0.1 у 0.0

-г

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

х

-4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4

| с* (х,у)-10

0-Н

у 0.05-1

(б)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

х

-2.2 -2 -1.8-1.6-1.4-1.2 -1 -0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

с* ( х, у )-10

(в)

Рис. 2. Пространственное распределение возмущений концентрации с*(х,у,г)=с(х,у,г)-с0(х,у) для параметров: а) Ре=1, Яр=0.05 б) Ре=1, Яр=30 в) Ре=10, Яр=30 в момент времени г=1.

0.1

у °.°5-

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

х

-1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2

(а)

р* ( х, у )-10

0.1

у 0.05

0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

х

-1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2

р* (х у)

0.

у 0.0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

х

р* (х, у )-10 1

(в)

Рис. 3. Пространственное распределение возмущений давления р*(х,у,г)=р(х,у,г)-р0(х,у) для параметров: а) Ре=1, Яр=0.05 б) Ре=1, Яр=30 в) Ре=10, Яр=30 в момент времени г=1

6. Результаты

эти параметры не варьируются (рис. 2-4). Главным образом структура течения определяется соотношением чисел Рэлея-Дарси и Пекле. Цель численного расчета - изучить структуру Значения чисел Рэлея-Дарси и Пекле были вы-течения. Исследование показало, что от парамет- браны для наиболее наглядной демонстрации ме-ров сорбции а и Ь зависит прежде всего интенсив- ханизмов возникновения конвективного течения. ность течения, поэтому в представленных данных Они не относятся к известным экспериментальным

6

3

10-2^

10-3 -Л

3

ад

10-4

10-5

10"2-э

10-

9

ад

10-4-

10-

10-6

(а) (б)

Рис. 4. а) Зависимость интегральной характеристики возмущений полей функции тока от параметра Яр. остальные параметры Ре=5, а=1, Ь=2 б) Зависимость интегральной характеристики возмущений полей функции тока от параметра Ре, остальные параметры Яр=5, а=1, Ь=2. Все расчеты производились для г=1

200 ______ 1.6 -1

160

120

а К

80

40

1.2-

3 0.8-

0.4-

Т

I

I

П

100

200 300 Ре

(а)

400

500

100

200 300

Ре

(б)

400 500

Рис. 5. а) Зависимость критического числа Рэлея-Дарси от числа Пекле. б) Зависимость частоты нейтральных возмущений от числа Пекле. Для параметров: а=35, Ь=7

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

жести (рис. 2, а, б). Увеличение числа Рэлея-Дарси приводит к росту отклонений концентрации

с*(х,у,г) = с(х,у,г) - с0(х,у)

и давления

р*(х,у,г) = р(х,у,г)-ро(х,у) от основного состояния. Из (рис. 2, в) и (рис. 3, в) видно, что увеличение числа Пекле делает более интенсивными возмущения вблизи выхода из области (правая граница).

На графике рис. 4 представлена зависимость интегральной характеристики (4.5), описывающей поведение возмущений функции тока от параметров задачи. Рис. 4, а демонстрирует стабилизирую-

0

0

0

0

о.

1.6 -I

1.2-

3 0.8-

0.4-

—I-1-1-1-1-1-1-1—

290 295 300 305

Ре

-I-1

310

285

(а)

—г—1—I—1—I—1—Г 290 295 300 305 Ре

(б)

Рис. 6. а) Зависимость критического числа Рэлея-Дарси от числа Пекле. б) Зависимость частоты нейтральных возмущений от числа Пекле. Для параметров: а=35, Ь = 7

200

200

200 220 240

260 280 Ре

300 320

320

(а)

(б)

Рис. 7. а) Зависимость критического числа Рэлея-Дарси от числа Пекле (а>Ь). б) Зависимость критического числа Рэлея-Дарси от числа Пекле (а<Ь) (правый график). Для параметров: (а=35, Ь = 7), (а=7, Ь=35). Синий цвет на нейтральной кривой соответствует колебательным возмущениям, черный цвет соответствует монотонным возмущениям

щее влияние роста числа Пекле: увеличение скорости потока уменьшает интенсивность возникшего конвективного течения. Рис 4, б показывает, что рост значений числа Рэлея-Дарси приводит к интенсификации возникшего конвективного течения.

На рис. 5 приведены нейтральная кривая и зависимость частоты нейтральных возмущений (5.16) от числа Пекле для значений параметров, соответствующих эксперименту [2].

Из рис. 5 видно, что при малых значениях числа Пекле возникает неустойчивость основного состояния монотонным образом при Rp~40, что соответствует работе [14]. В этом случае течение способствует возникновению конвективного дви-

жения. При дальнейшем увеличении числа Пекле происходит рост критического числа Рэлея-Дарси, т.е. основное течение становится более устойчивым. Для значений Ре > 200 наблюдаются области параметров, для которых возникает колебательная неустойчивость. Зависимость частоты таких колебаний от числа представлена на рис. 4, б. С целью более подробной иллюстрации смены режимов течения с монотонного на колебательный, на рис. 6 те же зависимости построены в диапазоне 285 < Ре < 310. Изломы нейтральной кривой на рис. 5, а при значениях числа Ре > 205 не являются ошибкой численного метода. Существует множество решений уравнения (5.15), на рис. 6, а по-

196.2-

196-

195.8-

195.6-

195.2-

195-

194.8-

0

194.6

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

На рис. 7 представлены те же зависимости, что и выше, но для различных соотношений параметров сорбции.

Из рис. 7 видно, что колебания возникают только в случае, когда параметр адсорбции а превышает параметр десорбции b (рис 7, а). Интервалы чисел Пекле в которых наблюдаются колебания выделены пунктирными линиями на рис. 7 (левая граница - черная линия, правая граница - красная линия). В обратном случае реализуются только переходы между различными решениями соответствующими монотонным возмущениям (рис. 7, б). Кроме того, видно, что с увеличением числа Пекле в диапазоне от 205 < Pe < 300 области колебательных возмущений становятся шире.

7. Заключение

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

Работы выполнены в рамках бюджетной темы № 121112200078-7

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

1. Schincariol R. A., Schwartz F. W. An experimental investigation of variable density flow and mixing in homogeneous and heterogeneous media // Water Resources Research. 1990. Vol. 26. N. 10. P. 2317-2329. DOI: 10.1029/WR026I010P02317

2. Ma^shev B. S., Khabin M. R., Evgrafova A. V. Identification of transport parameters for the solute filtration through porous media with clogging // Journal of Porous Media. 2023. Vol. 26. N. 6. P. 31-35. doi: 10.1615/JPorMedia.2022044645

3. Deans H. A. A mathematical model for dispersion in the direction of flow in porous media // Society of Petroleum Engineers Journal. 1963. Vol. 3. N. 01. P. 49-52. DOI: 10.2118/493-PA

4. van Genuchten M. T., Wierenga P. J. Mass transfer studies in sorbing porous media I. Analytical solutions 1 // Soil Science Society of America Journal. 1976. Vol. 40. N. 4. P. 473-480. DOI: 10.2136/sssaj1976.03615995004000040011x

5. Schumer R., Benson D. A., Meerschaert M. M., Baeumer B. Fractal mobile/immobile solute transport // Water Resources Research. 2003. Vol. 39. N. 10. DOI : 10.1029/2003WR002141

6. Gouze P., Le Borgne T., Lefrovost R., Lods G., Poidras T., Pezard P. Non-Fickian dispersion in porous media: 1. Multiscale measurements using single-well injection withdrawal tracer tests // Water Resources Research. 2008. Vol. 44. N. 6, W06426.

7. Selim H. M. Prediction of contaminant retention and transport in soils using kinetic multireaction models // Environmental Health Perspectives. 1989. Vol. 83. P. 69-75. DOI: 10.2136/sssaj2006.0422

8. Nield D., Bejan A. Convection in Porous Media. Cham, Switzerland: Springer, 2017. 988 p.

9. Кузнецов С. П. Динамический хаос. М.: Физма-тлит, 2001. 296 с.

10. Корн Г. А., Корн Т. М. Справочник по математике. СПб.: Лань, 2003. 832 с.

11.Хабин М. Р., Марышев Б. С. Идентификация параметров транспорта примеси через пористую среду с учётом закупорки // Вестник Пермского университета. Физика. 2021. №. 3. С. 44-55. DOI: 10.17072/1994-3598-2021-3-4455.

12. Самарский А. А. Введение в численные методы. СПб.: Лань, 2009. 288 c.

13. Михлин С. Г. Вариационные методы в математической физике. М.: Наука, 1970. C. 512

14. Horton C. W., Rogers F. T. Convection currents in a porous medium // Journal of Applied Physics. 1945. Vol. 16. N. 6. P. 367-370

References

1. Schincariol R. A., Schwartz F. W. An experimental investigation of variable density flow and mixing in homogeneous and heterogeneous media. Water Resources Research, 1990, vol. 26, no. 10. pp. 2317-2329.

2. Maryshev B. S., Khabin M. R., Evgrafova A. V. Identification of transport parameters for the solute filtration through porous media with clogging.

Journal of Porous Media, 2023, vol. 26, no. 6, pp. 31-35.

3. Deans H. A. A mathematical model for dispersion in the direction of flow in porous media Society of Petroleum Engineers Journal, 1963, vol. 3, no. 1. pp. 49-52. DOI : 10.2118/493-PA

4. Van Genuchten M. T., Wierenga P. J. Mass transfer studies in sorbing porous media I. Analytical solutions 1 Soil Science Society of America Journal, 1976, vol. 40, no. 4, pp. 473-480. DOI: 10.2136/sssaj1976.03615995004000040011x

5. Schumer R., Benson D. A., Meerschaert M. M., Baeumer B. Fractal mobile/immobile solute transport. Water Resources Research, 2003, vol. 39, no. 10, WR002141.

6. Gouze P., Le Borgne T., Leprovost R., Lods G., Poidras T., Pezard P. Non-Fickian dispersion in porous media: 1. Multiscale measurements using single-well injection withdrawal tracer tests. Water Resources Research. 2008, vol. 44, no. 6, W06426. DOI: 10.1029/2007WR006278

7. Selim H. M. Prediction of contaminant retention and transport in soils using kinetic multireaction

models Environmental Health Perspectives, 1989, vol. 83, pp. 69-75. DOI: 10.2136/sssaj2006.0422.

8. Nield D., Bejan A. Convection in Porous Media. Cham, Switzerland: Springer, 2017. 988 p.

9. Kuznetsov S. P. Dinamicheskii khaos (Dynamic Chaos). Moscow: Fizmatlit, 2001. 296 p. (In Russian)

10. Korn G. A., Korn T. M. Mathematical Handbook for Scientists and Engineers. New York: McGraw-Hill, 1968, 1130 p.

11. Khabin M. R., Maryshev B. S. Identification of the parameters of transport through porous media with clogging. Bulletin of the Perm University. Physics, 2021, no. 3, pp. 44-55 (In Russian).

12. Samarsky A. A. Vvedenye v chislennye metody (Introduction to numerical methods), Saint-Petersburg, Lan', 2009. 288 p. (In Russian).

13. Mikhlin S. G. Variatsionnye metody v ma-tematicheskoi fizike (Variational Methods in Mathematical Physics). Moscow: Nauka, 1970. 512 p. (In Russian).

14. Horton C. W., Rogers F. T. Convection currents in a porous medium. Journal of Applied Physics, 1945, vol. 16, no. 6, pp. 367-370.

Просьба ссылаться на эту статью в русскоязычных источниках следующим образом: Хабин М. Р., Марышев Б. С. Возникновение концентрационной конвекции в длинной прямоугольной области пористо среды // Вестник Пермского университета. Физика. 2023. № 4. С. 10-21. doi: 10.17072/1994-3598-2023-4-10-21

Please cite this article in English as:

Khabin M. R., Maryshev B. S. The onset of concentration convection in a long rectangular domain of a porous medium. Bulletin of Perm University. Physics, 2023, no. 4, pp. 10-23. doi: 10.17072/1994-3598-2023-4-10-21

Сведения об авторах

1. Михаил Романович Хабин, аспирант, инженер-исследователь, Институт механики сплошных сред УрО РАН, ул. Академика Королева, 1, Пермь, 614013; инженер, Пермский государственный национальный исследовательский университет, ул. Букирева, 15, Пермь, 614068.

2. Борис Сергеевич Марышев, канд. физ-мат. наук, н.с., Институт механики сплошных сред УрО РАН, ул. Академика Королева, 1, Пермь, 614013; доцент, Пермский государственный национальный исследовательский университет, ул. Букирева, 15, Пермь, 614068.

Author information

1. Mikhail R. Khabin, Postgraduate, Research Engineer, Institute of Continuous Media Mechanics UB RAS; 1, Akademika Koroleva st., Perm, 614013, Russia; Engineer, Perm State University; 15, Bukireva st., Perm, 614068, Russia.

2. Boris S. Maryshev, Candidate of Physical and Mathematical Sciences, Researcher, Institute of Continuous Media Mechanics UB RAS; 1, Akademika Koroleva st., Perm, 614013, Russia; Associate Professor, Perm State University; 15, Bukireva st., Perm, 614068, Russia.