МОДЕЛИРОВАНИЕ ДВУХФАЗНЫХ ТЕЧЕНИЙ ЖИДКОСТЕЙ В ПОРИСТОЙ СРЕДЕ В РЕЖИМЕ ДОМИНИРОВАНИЯ КАПИЛЛЯРНЫХ СИЛ Текст научной статьи по специальности «Физика»

ГЕОРЕСУРСЫ / GEORESOURCES gr//V^ 2020. Т- 22- № Ь С- 4-12

оригинальная статья

Б01: https://doi.Org/10.18599/grs.2020.1.4-12 ~ УДК 532.546

Моделирование двухфазных течений жидкостей в пористой среде в режиме доминирования капиллярных сил

Т.Р. Закиров1*, М.Г. Храмченков1,2

1 Казанский (Приволжский) федеральный университет, Казань, Россия 2Федеральный научный центр Научно-исследовательский институт системных исследований РАН, Москва, Россия

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

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

Для цитирования: Закиров Т.Р., Храмченков М.Г. (2020). Моделирование двухфазных течений жидкостей в пористой среде в режиме доминирования капиллярных сил. Георесурсы, 22(1), с. 4-12. DOI: о^/10.18599^.2020.1.4-12

1. Введение

Понимание процессов, происходящих при многофазных течениях в пористых средах, играет огромную роль во многих отраслях науки и техники, среди которых гидрогеология, почвоведение и грунтоведение, динамика подземных вод, освоение месторождений нефти и газа. В настоящее время в связи развитием экспериментальной и вычислительной техники опубликовано много работ по данной тематике (Liu и др., 2016; Chen и др., 2018; Zakirov и др., 2018a; Aursj0 и др., 2011; Tsuji и др., 2016). Эти работы объединяет исследование одного явления -развитие неустойчивостей или пальцеобразование при двухфазном течении. В зависимости от механики развития пальцев различают три режима потоков: течения с вязкими пальцами, капиллярными пальцами и со стабильным фронтом вытеснения. Выявлено, что режим течения контролируется двумя безразмерными параметрами - числом капиллярности Ca, характеризующим соотношение вязких и капиллярных сил, и соотношением вязкостей несмачивающей и смачивающей фаз M.

Течения со стабильным фронтом формируются при высоких Ca (более 10-3) и M > 1. Данный режим течения характеризуется плоской границей раздела фаз в гомогенных моделях (Liu и др., 2016), однако, как показано в

* Ответственный автор: Тимур Рустамович Закиров

E-mail: tirzakirov@kpfu.ru

© 2020 Коллектив авторов

работе (Tsuji и др., 2016), геометрия фронта в природных образцах, ввиду их неоднородности, не является стабильной. Течения с формированием вязких пальцев возникают при M< 0.1 и высоких Ca (более 10-3). В данном режиме вязкие силы трения намного превышают капиллярное сопротивление. Совместное течение двух жидкостей характеризуется нестабильным фронтом вытеснения и формированием тонких доминантных каналов по направлению внешнего перепада давления (Zakirov и др., 2018b; Tsuji и др., 2016).

При низких Ca (обычно меньше 10-5) капиллярные силы являются доминирующими. Некоторые особенности совместного течения описаны в работах (Zakirov и др., 2018b; Tsuji и др., 2016; Yamabe и др., 2015), согласно которым для данного режима характерны нестабильность фронта вытеснения, а также рост капиллярных пальцев во всех направлениях, в том числе противоположных внешнему перепаду давления. Однако большинство современных работ, изучающих механику вытеснения в режиме доминирования капиллярных сил, направлены на изучение скачков Хайнса. Данное явление впервые было обнаружено в работе (Haines и др., 1930) и характеризуется эпизодами быстрого вытеснения при перемещении мениска из узкого горла поры в ее широкое тело. Скачки Хайнса являются ярким примером проявления капиллярных сил в пористой среде.

В работах (Moebius, Or, 2012, 2014) на примере гранулярных моделей обнаружены падения давления,

связанные с перемещением мениска из узкого горла поры в ее широкое тело, однако механика вытеснения после скачков Хайнса детально не изучалась. Скорость движения мениска при скачке, согласно работам (Moebius, Or, 2012; Zacharoudiou, Boek, 2016; Armstrong, Berg, 2013; Armstrong и др., 2015), намного превышает среднюю скорость совместного течения. Показано, что быстрые события вытеснения оказывают влияние на распределение жидкости в порах, соседствующих с той, где произошел скачок (Armstrong и др., 2015). Вследствие высокой скорости мениска, как показано в работе (Armstrong и др., 2015), процессы вытеснения в режиме доминирования капиллярных сил не контролируются числами Ca и M, ввиду значимости сил инерции. В работе (Zacharoudiou, Boek, 2016) предлагается использовать число Охнесорга, связывающего вязкие, инерционные и капиллярные силы.

В настоящей работе процессы двухфазного течения изучаются методами математического моделирования. В последние годы решеточные уравнения Больцмана (LBM - Lattice Boltzmann methods) являются широко используемым инструментом для изучения двухфазных течений в масштабе поровых каналов (Liu и др., 2016; Zakirov и др., 2018b; Aursj0 и др., 2011; Tsuji и др., 2016; Yamabe и др., 2015). Для описания межфазных явлений применяется наиболее современная версия модели градиента цветового поля. LBM с использованием MRT-оператора (MRT - Multi relaxation time) столкновений и модель градиента цветового поля для трехмерных решеток (D3Q15, D3Q19, и D3Q27) сформулированы в работе (Leclaire и др., 2017).

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

Также настоящее исследование направлено на установление закона статистического распределения падений давлений, возникающих при скачках Хайнса, в трехмерном природном пористом пространстве. Исследуется влияние соотношения капиллярных и вязких сил на статистические характеристики.

2. Материалы и методы

2.1. Математическая модель

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

распределения_Дг,иД). Данная функция показывает долю частиц в момент времени t, находящихся в окрестности точки r(x,y,z) с координатами от х до х+Ах, от y до y+Ay, от z до z+Az и со скоростями в диапазоне от u(ux u, и) до u(ux+Aux, uy+Auy, и+Аи) (Leclaire и др., 2017).

В качестве возможных для перемещения ансамбля частиц направлений используется трехмерная модель D3Q19. Поскольку подробное описание данной модели уже приведено в работах (Leclaire и др., 2017, Zakirov и др., 2018b), мы ограничимся кратким представлением основных формул и этапов.

Описание динамики ансамбля частиц каждого из флюидов проводится в несколько этапов. Первый этап является этапом переноса («streaming step»). На данном этапе частицы за время At смещаются в соседние узлы, причем направление скорости движения не меняется. На втором этапе рассматривается процесс столкновения частиц («collision step»), в результате которого функция распределения частиц стремится к равновесному состоянию. На третьем этапе описывается взаимодействие жидкостей друг с другом на границе раздела, а также с твердой фазой. Эволюция функции распределения каждой жидкости во времени и пространстве описывается при помощи уравнения (1):

fk (r + eг. At, t + At) = fk (r, t) + (Qk (r, t))1 +

(Q k (r, t ))2. (1)

Индекс k = 1,2 указывает на вид жидкости, т.е. на смачивающую и несмачивающую фазы, соответственно; индекс i = (1-19) характеризует направление возможного перемещения частиц в модели D3Q19; (Q к- оператор столкновения, (Q.k )2 - оператор, описывающий межфазное взаимодействие и эффекты смачивания.

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

pk (r, t) = £ f (r, t);

i=1 1 2 19

u(r, t)=-YLef (r, t).

P к=1 i=1

(2)

(3)

Ввиду более высокой точности решений по сравнению с SRT-моделью (SRT - Single relaxation time) (Zakirov и др., 2018a), в настоящей работе для вычисления оператора (Qk У используется MRT - модель:

{ф1 =-M-lS{m1-mf). (4)

19

В уравнении (4) т. = ^М^. • ff. Вид матрицы M и

м

формулы для расчета mi для модели D3Q19 приведены в работах (Zakirov и др., 2018a). Диагональные компоненты матрицы S также описаны в работе (Zakirov и др., 2018a).

Давление Pk в LBM, создаваемое каждым флюидом, связано с его плотностью следующим соотношением: Pk = pkc2 /3, где c = Al / At - сеточная скорость, Al -шаг сетки.

Параметр релаксации является определяющим в LBM. Кинематическая вязкость жидкости контролируется при

НАУЧНО-ТВШИЧЕСКИЙ ЖУРНАЛ

www.geors.ru ГЕОРЕСУРСЫ

помощи параметра релаксации тк и связана с ним следующим соотношением:

"Я^Ж- (5)

Для описания явлений, происходящих на границах раздела жидкостей и жидкостей с твердыми поверхностями, в данной работе используется метод градиента цветового поля ^есЫге и др., 2017). На первом этапе проводится расчет градиента цветового поля g, компоненты которого вычисляются по формуле (6) ^есЫге и др., 2017):

19

g(r, о = 2 (р2(г+0 -р\г+0). (6)

¿=1

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

На втором этапе проводится вычисление оператора

2

(А' Г :

2 S

(7)

В уравнении (7) | ^ - модуль градиента цветового поля; Лк - численный параметр, контролирующий межфазное натяжение, при этом Л1 = А2 = А; В. принимает следующие значения: В, = -2/9; В , = 1/27; В 10 = 1/54; wл = 1/3;

1 2—/ 8-19 1

^2-7 = 1/18; ^8-19 = 1/36 - весовые коэффициенты.

Форма оператора (П*)2 не обеспечивает несмешиваемости жидкостей. Чтобы этого избежать, применяется модификация функций /,к . Она проводится после решения уравнения (1) с учетом операторов У и )2 ^есЫге и др., 2017):

( f-У =-fi-ß^^- f? • cos(a ) ; p p

(f2)* =—f f? • C0S(a,) .

(8)

(9)

P P

В уравнениях (8) и (9) f. = f + f2, p = pl + p2 ; f,eq - равновесная функция распределения (Zakirov и др., 2018b), рассчитанная при плотности р и нулевой скорости в каждом узле:

(e/ • и),лЛеi • u)2

- 1.5—).(10)

/Г (р, и) = ^ р-(1 + + 4.5

с с с

Параметр в в уравнениях (8) и (9) управляет толщиной границы раздела жидкостей. Он может принимать значения от 0 до 1 и в данной работе равен 0.8; а. - угол между направлениями g и е.. Краевой угол смачивания моделируется при помощи подхода, описанного в работе (Akai и др., 2018).

2.2. Граничные и начальные условия

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

оси OX. Вычислительный эксперимент проводится при заданной постоянной скорости течения в каждой ячейке входного сечения области фильтрации; обе жидкости вытекают через условно правую грань образца, перпендикулярную оси OX, давление на которой считается постоянным. В LBM такие условия задаются при помощи соотношений Зоу и Хе (Zou, He, 1997). Остальные четыре внешние грани расчетной области считаются непроницаемыми. В качестве граничных условий на непроницаемых областях используются условия «bounce back» (Zakirov и др., 2018b), являющиеся аналогом условий непротекания и прилипания в классической постановке задачи.

Математическая модель была численно реализована на языке Intel Fortran в среде разработки Visual Studio 2017. Численные валидации алгоритма на примерах слоистого течения жидкостей с разными вязкостями в плоском канале, а также при моделировании дренирования в однородном плоском канале, приведены в работах (Zakirov и др., 2018a, 2018b).

2.3. Образец исследования. Цифровая модель

В качестве образца пористой среды использовался нефтенасыщенный песчаник Ашальчинского месторождения (Республика Татарстан). Цифровая модель пористой среды, в которой проводились вычислительные эксперименты, получена при рентгеновском микротомографическом сканировании на томографе GE Phoenix v|tome|x s 240. Для сканирования был изготовлен кубический керн размером 3.7 мм и отсканирован с разрешением 4 мкм. Цифровая модель размером 3503 ячеек показана на рис. 1. оттенки серого на рис. 1а характеризуют интенсивность ослабления рентгеновского излучения различными участками образца: светло-серые оттенки соответствуют гранулам песчаника, а черные - поровому пространству. Результаты вычислительного эксперимента, представленные в разделе 3.1., были проведены на фрагменте с размером 250x125x125 ячеек, обладающем пористостью 0.272 и компонентой абсолютной проницаемости kvv = 0.612 мкм2.

XX

Статистические исследования скачков Хайнса, представленные в разделе 3.2, проводились на четырех фрагментах с размерами 350x175x175 ячеек. Их характеристики представлены в таблице 1. Поровое пространство образца 1 показано на рис. 1b и 1c.

fi * » v Vi

Wi

I г « « p* ■

Рис. 1. Цифровая микротомографическая модель песчаника с разрешением 4 мкм: а - исходная модель с размером 3503 ячеек; Ь - поровое пространство фрагмента с размером 350^175^175 ячеек; с - срез, перпендикулярный направлению ОУ

2

^ЩЩ GEDRESDURCES www.geors.ru

Номер образца Пористость, отн. ед. Абсолютная проницаемость &хх, мкм2

1 0.270 0.645

2 0.263 0.456

3 0.277 0.621

4 0.276 0.541

Табл. 1. Пористости и компоненты кХХабсолютной проницаемости четырех фрагментов

3. результаты и обсуждение

3.1. связь между событиями вытеснения и численными характеристиками течения

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

Численный эксперимент проведен на образце размером 250x125x125 ячеек при скорости и0 = 0.2х10"3 м/с, поверхностном натяжении а = 20 мН/м и краевом угле смачивания в = 10°. На рис. 2 показана динамика перепада давления (далее перепад давления между входной и выходной гранями и падение давления при скачке Хайнса будут обозначены как Р и АР, соответственно). Замеры давления проводились каждые 500хА/, что составляет 0.08 мс. Синяя кривая на рис. 2, построенная при а = 0 мН/м, характеризует величину вязких сил трения.

Согласно поведению красной кривой на рис. 2, наблюдаются флуктуации динамики Р (также рис. 7). Для усиления связи между характеристиками потока и событиями вторжения в работе проводится измерение удельной поверхности объема порового пространства £у, заполненного нагнетаемой фазой. Эта величина определяется как соотношение между площадью поверхности и занимаемым объемом. Дальнейшие результаты в данном разделе представлены для промежутка П = [115 175] (область интереса на рис. 2 и 3). На рис. 3Ь показана динамика изменения Бу, вычисляемая как А£г (П +1) = Бг (П +1) - Бг (П ).

Пошаговое заполнение порового пространства показано на рис. 4-6. Рис. 4а иллюстрирует распределение жидкостей, соответствующее позиции «а» и событию «1» на рис. 3; желтый цвет указывает на область, заполненную

Рис. 2. Динамика перепада давления между входной и выходной границами образца

Рис. 3. Динамика характеристик течения в промежутке п' = [115 ^ 175]: а - флуктуации перепада давления; Ь - флуктуации изменения удельной поверхности порового пространства, заполненного нагнетаемой фазой. События и позиции «а» -«п» ассоциированы с рис. 4-6. Красные квадраты указывают на появление новых подвижных менисков.

в промежутке П = [115 120]. Согласно рис. 3, рост капиллярного давления сопровождается возрастанием удельной поверхности £ поскольку АБу > 0 . Такое согласование характеристик объясняется обратными зависимостями удельной поверхности и капиллярного давления от радиуса порового канала.

Далее на рис. 3 а после события «1» наблюдается резкое падение давления. Вытеснение между событиями «1» и «3» будет происходить в области, отмеченной на рис. 4а и 4Ь красными линиями. Течение в данной области проиллюстрировано в более крупном масштабе на рис. 5. Расположение межфазного мениска в узком горле поры, который переместится в широкое тело поры в следующий момент времени, показано на рис. 5а (позиция «а» на рис. 3). Событие вторжения в данной области обусловлено превышением перепада давления Р капиллярного сопротивления, создаваемого данным мениском. Рис. 5Ь и 5 с показывают вытеснение, соответствующее позициям «Ь» и «с» на рис. 3. Согласно иллюстрациям, наблюдается быстрое смещение мениска в широкое тело. Такой характер вытеснения происходит вследствие резкого увеличения размера мениска. Данное наблюдение согласуется с поведением А£у на рис. 3Ь и сопровождается резким падением капиллярного сопротивления (рис. 3а). Описанное явление носит название «скачок Хайнса».

Рис. 5d и 5е иллюстрируют распределение жидкостей в позициях и «е». Вытеснение продолжается в данном поровом канале и отсутствует в остальных областях

НАУЧНО-ТЕХНИЧЕСКИЙ ЖУРНАЛ

www.geors.ru ГЕОРЕСУРСЫ

Рис. 4. Распределения нагнетаемой жидкости между событиями 1 и 7 (рис. 3). Голубой, розовый и желтый - несмачиваю-щая фаза; темно-красный - смачивающая фаза. Области заполнения между текущим и следующим событиями окрашены в желтый и розовый цвета.

порового пространства. После окончания быстрого вдавливания из узкого горла в широкое тело поры, наблюдается рост давления после события «2» (рис. 3а). Согласно рис. 5d и 5e, причиной роста капиллярного давления является уменьшение радиуса мениска при сужении радиуса поры на пути его перемещения, что также подтверждается ростом удельной поверхности в интервале между позициями «с» и «e» (рис. 3b). Рис. 4b иллюстрирует распределение жидкостей на момент события «3».

Вследствие уменьшения радиуса мениска, вероятность появления новых подвижных интерфейсов возрастает. Вытеснение будет продолжено в поровом канале, где капиллярное давление, создаваемое мениском, меньше, чем сопротивление текущего подвижного интерфейса (границы раздела). Рис. 5f и 5g (позиции «f» and «g» на рис. 3) показывают окончание вытеснения в данном канале. Переключение (смена) подвижных интерфейсов происходит между событиями «e» и «f» и соответствует событию «3». На рис. 3 показано увеличение капиллярного сопротивления в данном временном интервале и наибольшая скорость роста SV (моменты переключения подвижных менисков отмечены на рис. 3 красными квадратами). Новая область вытеснения при событии «4», окрашенная в розовый цвет, показана на рис. 4с.

Последующее вытеснение между событиями «4» и «5», согласно рис. 4d, будет происходить в новой области, выделенной красными рамками. Течение в данной области проиллюстрировано в более крупном масштабе на рис. 6. В позиции «h» (рис. 3), показанной на рис. 6а, мениск снова расположен на границе раздела узкого горла и поры. При этом перепад давления, очевидно, превышает капиллярное давление мениска, поскольку дальнейшее вытеснение происходит в данном поровом канале. Скачок Хайнса наблюдается в позициях «i» и «j» и сопровождается резким уменьшением капиллярного давления и удельной поверхности (рис. 3). Новая область, заполненная между событиями «4» и «5», окрашена в желтый на рис. 4d.

Рис. 5. Динамика перемещения межфазного мениска между событиями «1» и «4» (рис. 3) в области, отмеченной на рис. 4а красными линиями: а - расположение мениска в месте соединения узкого горла и широкого тела поры, позиция «а»; Ь, с -явление быстрого вдавливания нагнетаемой фазы в тело поры (скачок Хайнса), позиции «Ь» и «с»; d, е - окончание явления быстрого вдавливания, продолжение вытеснения при сужающимся размере поры, позиции «й» и «е»; /, g - локализация мениска в новом узком горле и окончание вытеснения в данной области, позиции «/» и «g». Голубой - несмачивающая фаза; темно-красный - смачивающая фаза.

Событие «5» на рис. 3 соответствует завершению скачка Хайнса. Последующее вытеснение в данной поре, как показано на рис. бё и бе (позиции «к» и «1» на рис. 3), имеет равномерный характер и происходит при сужении канала. Появление новых подвижных интерфейсов (рис. 4е) наблюдается после позиции «1» и соответствует наибольшей скорости роста удельной поверхности (красный квадрат после позиции «1» на рис. 3).

После события 6 снова было обнаружено резкое падение капиллярного давления и удельной поверхности. Согласно уже установленным связям, наблюдается скачок Хайнса в новой области порового пространства (рис. 4^.

Необходимо подчеркнуть, что, в общем случае, переключение подвижных интерфейсов и / или продолжение вытеснения в данной поре после скачка Хайнса является результатом «соревнования» между капиллярными давлениями, создаваемыми текущим подвижным мениском и остальными неподвижными менисками. Таким образом, выявленная смена подвижных интерфейсов, наряду со скачками Хайнса, является еще одной отличительной особенностью двухфазных течений при доминировании капиллярных сил.

Основываясь на выявленных закономерностях, определен сценарий дренажного вытеснения, состоящий из последовательности следующих событий: 1 - равномерное вытеснение в узких поровых каналах при растущем капиллярном давлении и удельной поверхности; 2 - скачок Хайнса, сопровождающийся быстрым падением давления и удельной поверхности; 3 - окончание скачка Хайнса и продолжение равномерного вытеснения в данном поровом канале; возобновление роста капиллярного давления и удельной поверхности; 4 - возможное, в общем случае,

SCIENTIFIC AND TECHNICAL JOURNAL

GEDRESDURCES www.geors.ru

ПОЗ. «Ь» "...... - - - 1КУ5. «i»

Горло поры

Мениск в горле поры jj Тело поры Ь

v £

I_* % ^

Скачок Хайнса

1с

*

у-«-'

мох «1» Горло поры но. «т» -4 а

/ \ \ * , ' • Л- :

I ^ Ч, ! \»вая облас

Горло поры . ' лЛ ' течения

^^^^^^ й" I

Рис. 6. Динамика вытеснения между событиями «4» и «6» в области, отмеченной красными рамками на рис. 4d: a - расположение мениска в месте соединения узкого горла и широкого тела поры, позиция «h»; b, c - скачок Хайнса, позиции «i» и «j»; d-f - окончание скачка Хайнса, продолжение равномерного течения в данной поре при сужающемся радиусе, появление новых подвижных интерфейсов, позиции «k»-«m». Голубой -несмачивающая фаза; темно-красный - смачивающая фаза.

продолжение вытеснения в данной поре и / или появление новых подвижных интерфейсов; 5 - скачок Хайнса в новой области порового пространства. Последовательность событий 2-4 может считаться периодом дренажного течения, повторяющимся, согласно рис. 2 и 7, с различной длительностью и амплитудой флуктуаций давления.

3.2. Статистика скачков Хайнса

В данном разделе проводится исследование влияния скорости потока и межфазного натяжения на статистическое распределение падений давлений. Стоит отметить, что в отличие от двухмерных микромоделей, статистические исследования для трехмерных объектов затруднены и ограничены. Одна из причин заключается в высоком разрешении модели и, следовательно, ее малом объеме. Количество событий в 3D моделях, как правило, значительно меньше, чем в 2D образцах (более 500) (Furuberg и др., 1996; Soto и др., 2017). Кроме того, исследование затрудняется одновременными событиями вторжения, тогда как для статистической интерпретации необходимо наблюдение за единичными процессами. Для уменьшения частоты одновременных скачков Хайнса и получения достаточного количества этих событий вычислительные эксперименты проводятся на четырех небольших фрагментах пористых сред, свойства которых указаны в таблице 1.

3.2.1. Влияние скорости потока

Расчеты проведены в широком диапазоне скорости нагнетаемой жидкости: [0.05^1.5]*10-3 м/с: а = 20 мН/м во всех вариантах. Число Ca и средние оценки соотношения между капиллярными (Pcap) и вязкими силами

(Р™с) приведены в таблице 2. Все расчеты выполнялись до достижения 0.3 поровых объемов закачанной жидкости.

Динамика перепада давления для фрагмента №2, построенная для трех значений скорости, показана на рис. 7. Согласно столбцу 5 таблицы 2, количество обнаруженных падений давления (№р) сильно зависит от скорости потока. Эта особенность отчетливо видна на рис. 7. Также стоит отметить, что величина №гор не равна количеству единичных скачков Хайнса. Выявлено, что увеличение скорости потока приводит к росту количества одновременных скачков Хайнса. В результате большинство единичных скачков, особенно при высоких скоростях, не могут быть идентифицированы при единичном падении давления.

При визуализации течения обнаружено, что для и0 > 0.4х10-3 м/с три и более одновременных скачка происходят при единичном падении давления, тогда как для и0 < 0.2х10-3 м/с количество одновременных скачков не превышает двух. Данная особенность делает невозможным корректный статистический анализ при и0 > 0.4х10-3 м/с.

Для и0 = 0.2х10-3 м/с доля одновременных парных скачков при единичном падении давления составляет около 30 % и уменьшается до 5 % при и0 = 0.05*10-3 м/с. Влияние скорости потока на количество единичных скачков существует, но является незначительным в данном диапазоне. Это наблюдение также подтверждается вычисленным средним поровым объемом, заполненным нагнетаемой жидкостью в течение единичного падения давления (РУ)ишр в столбце 7, табл. 2). Для и0 < 0.2*10-3 м/с увеличение РУ-)итр в результате одновременных парных скачков составляет

Рис. 7. Динамика перепада давления при различных скоростях двухфазного течения

№ w0x10-3, (м/с) logwCa pcap/ pvisc N drop Среднее значение AP, (кПа) Pyjumpx10-3

1 0.05 -5.6 67 360 0.206 0.68

2 0.1 -5.3 33.5 307 0.173 0.81

3 0.2 -4.9 16.7 271 0.165 0.90

4 0.4 -4.7 8.3 193 0.152 1.34

5 0.8 -4.4 4.2 159 0.129 2.14

6 1.5 -4.1 2.0 119 0.111 2.07

Табл. 2. Численные параметры вычислительных экспериментов и результаты влияния скорости потока на статистику скачков Хайнса

НАУЧНО-ТЕХНИЧЕСКИЙ ЖУРНАЛ

www.geors.ru ГЕРРЕСУРСЫ

9

около 25 %, тогда как при высоких скоростях прирост pyjump является слишком большим.

Статистические распределения падений давления, нормированных к AP (столбец 6, табл. 2), для u0 < 0.2*10-3 м/с показаны на рис. 8а. Выявлено, что полученные распределения подчиняются логнормальному закону (сплошные линии на рис. 8а) с R2 = 0.93, 0.96 и 0.94 по мере возрастания u0. Можно заключить, что скорость потока оказывает незначительный эффект на распределения, но отрицательно влияет на среднее значение падения давления.

Согласно экспериментам, проводимым на 2D гомогенных микромоделях (Moebius, Or, 2012; Furuberg и др., 1996; Soto и др., 2017), кумулятивное распределение падений давления N подчиняется следующему закону:

1п(Ж) = fi-AP/ÁP . (11)

Кумулятивные распределения показаны на рис. 8b. Сплошные линии являются результатом аппроксимации

Рис. 8. Статистические распределения для различных скоростей потока: а - распределение относительной частоты; Ь -кумулятивные распределения N. Круглые символы показывают измеренные данные, сплошные линии - результаты аппроксимации.

при помощи соотношения (11) с R2 > 0.98. Таким образом, падения давления в трехмерных гетерогенных моделях, двухмерных гранулярных образцах (Moebius, Or, 2012; Furuberg и др., 1996) и почвах (Soto и др., 2017) подчиняются одному закону распределения. Близкие значения в равные -0.970, -0.983 и -0.969 по мере возрастания u0 также подтверждают слабое влияние расхода на статистические параметры.

3.2.2. Влияние межфазного натяжения

В данном разделе представлены результаты вычислительных экспериментов, проведенных при а = 13, 20 и 28 мН/м; скорость потока u0 = 0.1*10-3 м/с фиксирована. Все расчеты также выполнялись до достижения 0.3 по-ровых объемов закачанной жидкости.

Параметры и результаты расчетов представлены в таблице 3. Согласно приведенным данным, межфазное натяжение практически не влияет на количество падений давления (столбец 5, табл. 3), но прямо пропорционально воздействует на AP (столбец 6, табл. 3). Близкие значения PVjump (столбец 7, табл. 3) также характеризуют слабое влияние с на количество единичных скачков.

Статистические распределения представлены на рис. 9. Характер сплошных линий на рис. 9a, являющихся результатом аппроксимации по логнормальному закону с R2 > 0.93, характеризует слабое влияние межфазного натяжения на распределения.

Принимая во внимание, что в ~ -1 независимо от скорости течения и межфазного натяжения, можно заключить, что распределение AP / AP определяется только геометрией порового пространства. Например, в = -1.3 и -1.5 для двух образцов двухмерных моделей, рассмотренных в (Moebius, Or, 2012); в = -1.31 для двухмерной гранулярной модели, испытуемой в (Furuberg и др., 1996); в ~ -3 для почв (Soto и др., 2017).

Выводы

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

Статистические исследования падений давления при скачках Хайнса показали, что распределения подчиняются логнормальному закону. Показано, что параметры

№ а, мН/м logwCa p cap / p vise Ndrop Среднее значение AP, (kPa) PVjumpx10-3

1 13 -5.1 21.7 289 0.115 0.87

2 20 -5.3 33.5 307 0.173 0.81

3 28 -5.45 47 300 0.241 0.92

Табл. 3. Численные параметры вычислительных экспериментов и результаты влияния межфазного натяжения на статистику скачков Хайнса

Рис. 9. Статистические распределения для различных межфазных натяжений: а - распределение относительной частоты; Ь -кумулятивные распределения N. Круглые символы показывают измеренные данные, сплошные линии - результаты аппроксимации.

двухфазного течения не оказывают влияния на статистические распределения и воздействуют только на среднее значение падений давления.

Финансирование/Благодарности

Работа выполнена при поддержке гранта РФФИ №18-31-00134, а также за счет средств субсидии, выделенной в рамках государственной поддержки Казанского (Приволжского) федерального университета в целях повышения его конкурентоспособности среди ведущих мировых научно-образовательных центров.

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

литература

Akai T., Bijeljic B., Blunt M.J. (2018). Wetting boundary condition for the color-gradient lattice Boltzmann method: Validation with analytical and experimental data. Advances in Water Resources, 116, pp. 56-66. https://doi. org/10.1016/j.advwatres.2018.03.014

Armstrong R.T., Berg S. (2013). Interfacial velocities and capillary pressure gradients during Haines jumps. Physical Review E, 88, 043010. https://doi.org/10.1103/PhysRevE.88.043010

Armstrong R.T., Evseev N., Koroteev D., Berg S. (2015). Modeling the velocity field during Haines jumps in porous media. Advances in Water Resources, 77, pp. 57-68. https://doi.org/10.1016Zj.advwatres.2015.01.008 Aursj0 O., L0voll G., Knudsen H.A., Flekk0y E.G., Mal0y K.J. (2011). A Direct Comparison Between a Slow Pore Scale Drainage Experiment and a 2D Lattice Boltzmann Simulation. Transport in Porous Media, 86(1), pp. 125-134. https://doi.org/10.1007/s11242-010-9611-y

Chen Y.-F., Wu D.-S., Fang Sh., Hu R (2018). Experimental study on two-phase flow in rough fracture: Phase diagram and localized flow channel. International Journal of Heat and Mass Transfer. 122, pp. 1298-1307. https:// doi.org/10.1016/j.ijheatmasstransfer.2018.02.031

Furuberg L., Mal0y K.J., Feder J. (1996). Intermittent behavior in slow drainage. Physical Review E, 53, pp. 966-977. https://doi. org/10.1016/0378-4371(92)90542-X

Haines W.B. (1930). Studies in the physical properties of soil. Journal of Agricultural Science, 20, pp. 98-116. https://doi.org/10.1017/ S002185960008864X

Leclaire S., Parmigiani A., Malaspinas O., Chopard B., Latt J. (2017). Generalized three-dimensional lattice Boltzmann color-gradient method for immiscible two-phase pore-scale imbibition and drainage in porous media. Physical Review E, 95, 033306. https://doi.org/10.1103/PhysRevE.95.033306 Liu H., Kang Q., Leonardi C.R., Schmieschek S., Narvaez A., Jones B.D., Williams J.R., Valocchi A.J., Harting J. (2016). Multiphase lattice Boltzmann simulations for porous media applications. Computational Geosciences, 20(4). pp. 777-805. https://doi.org/10.1007/s10596-015-9542-3

Moebius F., Or D. (2012). Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries. Journal of Colloid and Interface Science, 377, pp. 406-415. https://doi.Org/10.1016/j.jcis.2012.03.070 Moebius F., Or D. (2014). Pore scale dynamics underlying the motion of drainage fronts in porous media. Water Resources Research, 50, pp. 84418457. https://doi.org/10.1002/2014WR015916

Soto D., Paradelo M., Corral Á., López Periago J.E. (2017) Pressure jumps during drainage in macroporous soils. Vadose Zone Journal, 16. https:// doi.org/10.2136/vzj2017.04.0088

Tsuji T., Jiang F., Christensen K.T. (2016). Characterization of immiscible fluid displacement processes with various capillary numbers and viscosity ratios in 3D natural sandstone. Advances in Water Recourses. 95, pp. 3-15. https://doi.org/10.1016Zj.advwatres.2016.03.005

Yamabe H., Tsuji T., Liang Y., Matsuoka T. (2015). Lattice Boltzmann Simulations of Supercritical CO2-Water Drainage Displacement in Porous Media: CO2 Saturation and Displacement Mechanism. Environmental Science and Technology, 49, pp. 537-543. https://doi.org/10.1021/es504510y

Zacharoudiou I., Boek E.S. (2016). Capillary filling and Haines jump dynamics using free energy Lattice Boltzmann simulations. Advances in Water Resources, 92, pp. 43-56. https://doi.org/10.1016/j.advwatres.2016.03.013 Zakirov T.R., Galeev A.A., Statsenko E.O., Khaidarova L.I. (2018a). Calculation of filtration characteristics of porous media by their digitized images. Journal of Engineering Physics and Thermophysics, 91(4), pp. 1069-1078. https://doi.org/10.1007/s10891-018-1833-9

Zakirov T.R., Galeev A.A., Khramchenkov M.G (2018b). Pore-scale Investigation of Two-Phase Flows in Three-Dimensional Digital Models of Natural Sandstones. Fluid Dynamics. 53(5), pp. 76-91. https://doi. org/10.1134/S0015462818050087

Zou Q., He X. (1997). On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Phys. Fluids, 9, pp. 1591-1598. https:// doi.org/10.1063/1.869307

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

Тимур Рустамович Закиров - кандидат физ.-мат. наук, доцент, Институт геологии и нефтегазовых технологий Казанский (Приволжский) Федеральный университет Россия, 420008, Казань, ул. Кремлёвская, д. 18 E-mail: tirzakirov@kpfu.ru

Максим Георгиевич Храмченков - доктор физ.-мат. наук, профессор, заведующий кафедрой математических методов в геологии, Казанский (Приволжский) федеральный университет; ведущий научный сотрудник, Федеральный научный центр Научно-исследовательский институт системных исследований РАН

Россия, 420008, Казань, ул. Кремлёвская, д. 18

Статья поступила в редакцию 16.09.2019;

Принята к публикации 25.10.2019; Опубликована 30.03.2020

НЮЧНО-7ЕХНИЧЕСКИЙ ЖУРНАЛ

www.geors.ru ГЕОРЕСУРСЫ

reOPECyPCbl / GEORESOURCES

2020. T. 22. № 1. C. 4-12

Pore-scale investigation of the displacement fluid mechanics during two-phase flows in natural porous media under the dominance of capillary forces

T.R Zakirov1*, M.G. Khramchenkov1,2

'Kazan Federal University, Kazan, Russian Federation

2Scientific Research Institute for System Analysis of the Russian Academy of Sciences, Moscow, Russian Federation *Corresponding authors: Timur R. Zakirov, e-mail: tirzakirov@kpfu.ru

Abstract. This paper presents the results of numerical simulations of two-phase flows in porous media under capillary forces dominance. For modeling of immiscible two-phase flow, the lattice Boltzmann equations with multi relaxation time operator were applied, and the interface phenomena was described with the color-gradient method. The objective of study is to establish direct links between quantitative characteristics of the flow and invasion events, using high temporal resolution when detecting simulation results. This is one of the few works where Haines jumps (rapid invasion event which occurs at meniscus displacing from narrow pore throat to its wide body) are considered in three-dimensional natural pore space, but the focus is also on the displacement mechanics after jumps. It was revealed the sequence of pore scale events which can be considered as a period of drainage process: rapid invasion event during Haines jump; finish of jump and continuation of uniform invasion in current pore; switching of mobile interfaces and displacement in new region. The detected interface change, along with Haines jump, is another distinctive feature of the capillary forces action. The change of the mobile interfaces is manifested in step-like behavior of the front movement. It was obtained that statistical distributions of pressure drops during Haines jumps obey lognormal law. When investigating the flow rate and surface tension effect on the pressure drop statistics it was revealed that these parameters practically don't affect on the statistical distribution and influence only on the magnitude of pressure drops and the number of individual Haines jumps.

Keywords: two-phase flow, Haines jumps, capillary forces, lattice Boltzmann equations, color-gradient method, X-ray computed tomography

Recommended citation: Zakirov T.R., Khramchenkov M.G. (2020). Pore-scale investigation of the displacement fluid mechanics during two-phase flows in natural porous media under the dominance of capillary forces. Georesursy = Georesources, 22(2), pp. 4-12. DOI: https:// doi.org/10.18599/grs.2020.1.4-12

References

Akai T., Bijeljic B., Blunt M.J. (2018). Wetting boundary condition for the color-gradient lattice Boltzmann method: Validation with analytical and experimental data. Advances in Water Resources, 116, pp. 56-66. https://doi. org/10.1016/j.advwatres.2018.03.014

Armstrong R.T., Berg S. (2013). Interfacial velocities and capillary pressure gradients during Haines jumps. Physical Review E, 88, 043010. https://doi.org/10.1103/PhysRevE.88.043010

Armstrong R.T., Evseev N., Koroteev D., Berg S. (2015). Modeling the velocity field during Haines jumps in porous media. Advances in Water Resources, 77, pp. 57-68. https://doi.org/10.1016Zj.advwatres.2015.01.008 Aursj0 O., L0voll G., Knudsen H.A., Flekk0y E.G., Makiy K.J. (2011). A Direct Comparison Between a Slow Pore Scale Drainage Experiment and a 2D Lattice Boltzmann Simulation. Transport in Porous Media, 86(1), pp. 125-134. https://doi.org/10.1007/s11242-010-9611-y

Chen Y.-F., Wu D.-S., Fang Sh., Hu R (2018). Experimental study on two-phase flow in rough fracture: Phase diagram and localized flow channel.

International Journal of Heat and Mass Transfer. 122, pp. 1298-1307. https:// doi.org/10.1016/j.ijheatmasstransfer.2018.02.031

Furuberg L., Mal0y K.J., Feder J. (1996). Intermittent behavior in slow drainage. Physical Review E, 53, pp. 966-977. https://doi. org/10.1016/0378-4371(92)90542-X

Haines W.B. (1930). Studies in the physical properties of soil. Journal of Agricultural Science, 20, pp. 98-116. https://doi.org/10.1017/ S002185960008864X

Leclaire S., Parmigiani A., Malaspinas O., Chopard B., Latt J. (2017). Generalized three-dimensional lattice Boltzmann color-gradient method for immiscible two-phase pore-scale imbibition and drainage in porous media. Physical Review E, 95, 033306. https://doi.org/10.1103/PhysRevE.95.033306 Liu H., Kang Q., Leonardi C.R., Schmieschek S., Narváez A., Jones B.D., Williams J.R., Valocchi A.J., Harting J. (2016). Multiphase lattice Boltzmann simulations for porous media applications. Computational Geosciences, 20(4). pp. 777-805. https://doi.org/10.1007/s10596-015-9542-3

Moebius F., Or D. (2012). Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries. Journal of Colloid and Interface Science, 377, pp. 406-415. https://doi.org/10.1016/jjcis.2012.03.070 Moebius F., Or D. (2014). Pore scale dynamics underlying the motion of drainage fronts in porous media. Water Resources Research, 50, pp. 84418457. https://doi.org/10.1002/2014WR015916

Soto D., Paradelo M., Corral Á., López Periago J.E. (2017) Pressure jumps during drainage in macroporous soils. Vadose Zone Journal, 16. https:// doi.org/10.2136/vzj2017.04.0088

Tsuji T., Jiang F., Christensen K.T. (2016). Characterization of immiscible fluid displacement processes with various capillary numbers and viscosity ratios in 3D natural sandstone. Advances in Water Recourses. 95, pp. 3-15. https://doi.org/10.1016Zj.advwatres.2016.03.005

Yamabe H., Tsuji T., Liang Y., Matsuoka T. (2015). Lattice Boltzmann Simulations of Supercritical CO2-Water Drainage Displacement in Porous Media: CO2 Saturation and Displacement Mechanism. Environmental Science and Technology, 49, pp. 537-543. https://doi.org/10.1021/es504510y

Zacharoudiou I., Boek E.S. (2016). Capillary filling and Haines jump dynamics using free energy Lattice Boltzmann simulations. Advances in Water Resources, 92, pp. 43-56. https://doi.org/10.1016/j.advwatres.2016.03.013 Zakirov T.R., Galeev A.A., Statsenko E.O., Khaidarova L.I. (2018a). Calculation of filtration characteristics of porous media by their digitized images. Journal of Engineering Physics and Thermophysics, 91(4), pp. 1069-1078. https://doi.org/10.1007/s10891-018-1833-9

Zakirov T.R., Galeev A.A., Khramchenkov M.G (2018b). Pore-scale Investigation of Two-Phase Flows in Three-Dimensional Digital Models of Natural Sandstones. Fluid Dynamics. 53(5), pp. 76-91. https://doi. org/10.1134/S0015462818050087

Zou Q., He X. (1997). On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Phys. Fluids, 9, pp. 1591-1598. https:// doi.org/10.1063/1.869307

About the Authors

Timur R. Zakirov - Cand. Sci. (Physics and Mathematics), Associate Professor, Institute of Geology and Oil and Gas Technologies, Kazan Federal University

18, Kremlevskaya st., Kazan, 420008, Russian Federation E-mail: tirzakirov@kpfu.ru

Maxim G. Khramchenkov - Dr. Sci. (Physics and Mathematics), Professor, Head of the Department of Mathematical Methods in Geology, Kazan Federal University; Leading Researcher, Scientific Research Institute for System Analysis of the Russian Academy of Sciences 18, Kremlevskaya st., Kazan, 420008, Russian Federation

Manuscript received 16August 2019;

Accepted 25 October 2019;

Published 30March 2020