CC BY
5
УДК 519.6
DOI 10.33764/2618-981X-2022-2-2-32-37
Моделирование нестационарных течений многофазной вязкой жидкости в поровом пространстве
Е. А. Гондюл12*, В. В. Лисица3 1 Новосибирский национальный государственный университет, г. Новосибирск,
Российская Федерация 2 Институт математики СО РАН, г. Новосибирск, Российская Федерация 3 Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, г. Новосибирск, Российская Федерация *e-mail: e.gondyul@g.nsu.ru
Аннотация. В данной работе разработан и реализован численный алгоритм моделирования нестационарных течений вязкой многофазной изотермической жидкости методом конечных разностей с использованием проекционного метода для численного решения уравнения Навье-Стокса. Проекционный метод подразумевает расщепление исходной системы уравнений по физическим процессам, в котором раздельно учитываются конвективный перенос и действие градиента давления. В результате на каждом шаге по времени необходимо решать уравнение Пуассона для нахождения поля давления. Решение СЛАУ выполняется параллельным прямым решателем, основанный на LU разложении. Используется явная схема для решения уравнения Кана-Хилларда для обновления фазового поля, параметр которого учитывается при добавлении к уравнению Навье-Стокса поверхностных сил. Приведены результаты вычислительных экспериментов, показывающие качественно и количественно совпадение с экспериментальными и численными данными из литературы.
Ключевые слова: многофазная жидкость, численное моделирование
Modeling of unsteady flows of multiphase viscous fluid in a pore space
Е. А. Gondyul12*, V. V. Lisitsa3
1 Novosibirsk State University, Novosibirsk, Russia 2 Institute of Mathematics SB RAS, Novosibirsk, Russia 3Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, Russia
*e-mail: e.gondyul@g.nsu.ru
Abstract. The authors have developed and implemented a numerical algorithm to model unsteady flows of a viscous multiphase isothermal fluid by finite difference method using the projection method for the numerical solution of the Navier-Stokes equation. The projection method implies splitting the initial system of equations by physical processes, in which convective transport and the effect of the pressure gradient are separately taken into account. As a result, at each step, it is necessary to solve the Poisson equation to find the pressure field. The solution of SLAE is performed by a parallel direct solver based on LU decomposition. An explicit scheme is used to solve the Cahn-Hilliard equation to update the phase field, the parameter of which is taken into account when adding surface forces to the Navier-Stokes equation. Computational experiments showing qualitative and quantitative agreement with experimental and numerical data from the literature are presented.
Keywords: multiphase liquid, numerical simulation
Введение
Численное моделирование многофазных течений в областях с топологически сложной структурой позволяет дополнить и расширить лабораторные петро-физические исследования на керне. В нефтегазовой промышленности это требуется для получения корректных данных о параметрах нефтегазоносного пласта, повышения точности описания свойств системы "флюид-порода", оптимизации параметров вытеснения и выбора оптимального вытесняющего агента и т.д. Для описания течения многофазной жидкости в поровом пространстве используется модель фазового поля, параметр которого непрерывно меняется в исследуемой области на границах раздела фаз. Таким образом, граница заменяется тонкой, но ненулевой по толщине межфазной областью, в которой поверхностные силы меняются плавно. Модель фазового поля начала развиваться в работах [1,2]. В их работах был сделан вывод общего уравнения свободной энергии для неоднородной системы, который послужил отправной точкой для создания полной модели фазового поля. Классически уравнение Кана-Хилларда выводится путём минимизации свободной энергии Гинзбурга-Ландау. Популярность в настоящее время данной математической модели обусловлено термодинамически согласованной теорией и относительно лёгкой реализацией.
Наравне с другими методами, такими как level-set [3] и VOF (volume-of-fluid) [4], метод фазового поля позволяет явно не отслеживать резкую границу раздела фаз и не требует пересчитывать на каждом итерационном шаге кривизну для подсчёта силы поверхностного натяжения. При правильном выборе параметров системы метод достаточно точно описывает динамику границы и сохраняет объёмные фазы.
Постановка задачи
Рассмотрим открытую область П = П-l +П2, QER2 с границей 5П (рис.1), а также некоторый непустой промежуток времени [0, Т].
Задача заключается в нахождении распределения фазового параметра ф(х, t), который определяется следующим выражением:
Ф= \ ,-1 <ф<1, т1 +т2
где тг, т2 — массы первого и второго флюидов.
Решается система уравнений, описывающая изотермические течения вязкой многофазной жид- рис ^ Схематичное кости в области со сложной геометрией с учётом изображение ИССЛедуемой поверхностных сил, состоящая из уравнений На- области
вье-Стокса и Кана-Хилларда:
p(ut + u-Vu) = -Vp + V • [r)(Vu + VuT)] +SF + pg,
V-u = 0, 0t + u-V0 = V-(M($)Vii),
№
ß = F'(0) -еЧф,
где p{x, t), u(x, t), p(x, t), ц(х, t) — плотность, вектор скорости, давление и вязкость среды соответственно, SF(ф) — поверхностное натяжение, д — ускорение свободного падения, ф(х, t) — параметр фазового поля, М{ф), р{ф), е - подвижность, химический потенциал, свободная энергия Гельмгольца на единицу объема соответственно, константа, отвечающая за ширину границы раздела фаз.
В данной работе, подвижность М(ф)принимает константное значение, как в работе [5]. Свободная энергия Гельмгольца представлена на рис.2 и принимает следующее выражение:
¥{ф~) = 0.25(ф2 - 1)2.
Для моделирования влияния поверхностного натяжения на движение жидкости используется CSF (continuum surface force) [6]. Данный метод подразумевает, что вместо использования классического условия скачка на давление на границе раздела фаз, вводится сингулярная сила, непрерывно действующая в области перехода. Такая сила пропорциональная кривизне границе, поэтому в большинстве литературы [7] и нашей работе используется следующий вид поверхностного натяжения:
SF = -6V2 |V0|V0,
где о — коэффициент поверхностного натяжения.
Для решения уравнения Навье-Стокса используется проекционный метод второго порядка, изложенный в [9,10]. Применительно к решению задачи о течении несжимаемой вязкой жидкости проекционный метод может трактоваться как расщепление исходной системы уравнений по физическим процессами, в котором раздельно учитывается конвективный перенос и действие градиента давления. Одним из шагов алгоритма является решение уравнения Пуассона на давление, которое реализовывалось с помощью прямого решателя Intel MKL PARDISO. Уравнение Кана-Хилларда решалось с помощью явной конечно-разностной схемой второго порядка.
Численные эксперименты
Для верификации численной модели в случае эволюции фазового поля с различными углами смачивания без гравитационного и скоростного поля выполнялась серия вычислительных экспериментов, моделирующая формирование равновесной формы капли из начальной формы с углом смачивания в^ — 90° в равновесную форму с углом смачивания в, лежащей на твердой поверхности. Для
Рис. 2. Свободная энергия Гельмгольца
-lioo -0Л5 -0.50 -0^25 0.00 0.25 0.50 0.75 1.00 Параметр фазового поля ф
этого случая существует геометрически выведенная формула [11] расчёта радиуса кривизны, контакта и высоты капли (рис. 3):
Rf — R(
п
xf~nU 2в-зт2в' Tf — Rf sin в,
hf = Rf(l -cose),
где Rf — радиус кривизны, R0 — начальный радиус капли, Tf — радиус контакта, hf — высота капли, в — угол смачивания.
/Г / / < \ ■— ^^ / N / Rn Y\
т Jr \
R,
Рис.3. Схематическое изображение начальной и равновесной формы капли
с углом смачивания в
Моделирование выполняется на сетке размером Ых — Ыу — 200 с шагом по пространству кх — ку — 5 • 10"5 и шагом по времени Ы — 5.52 • 10"11. На рис. 4 показаны хорошо согласованные с теорией экспериментальные данные.
В нестационарном случае на рис. 5 показан качественно перенос фазового поля под действием силы тяжести в 6 миллиметровой трубке. Перепад давления составлял 1 атм. Параметры расчётной сетки: кх — ку — 1 • 10-4, Ых — 60, А1у = 100, 10"4.
2.5 -2.0 -1.5 -
1.0 -0.5 -0.0 -
Рис. 4 Зависимость параметров капли от угла смачивания
25 50 75 100 125 150 175
в°
Таблица 1
Параметры среды
Параметр Значение Единица измерения
Pi 1220 кг/м3
Р2 1.2 кг/м3
т 0.116 Н-м
т 2- 10"5 Н-м
а 0.063
М 1
е min (hx, hy)
О 20 40 60 0 20 40 60 о 20 40 60 о 20 40 60
Рис. 5 Отрыв капли
Расчёты проводились на CPU, время расчёта для случая, когда сетка имеет следующее разрешение: Nx — Ny — 200 и hx — hy — 1 • 10 "3 с шагом по времени At — 1 • 10-4 с одна итерация составляет ~44 с. Для повышения производительности в дальнейшей работе предполагается использовать GPU.
Заключение
В работе представлен численный алгоритм для моделирования течений вязкой многофазной жидкости. Используется проекционный метод, в котором необходимо решать уравнение Пуассона для коррекции давления. Алгоритм включает реализацию граничных условий с заданным углом смачивания. Проведена серия вычислительных экспериментов, имитирующих перенос двухфазной жидкости в поле скоростей и позволяющих сравнить результаты с экспериментальными [12] и численными данными, представленные в литературе [11,13].
Благодарности
Е.А. Гондюл разработала и реализовала алгоритм моделирования течений двухфазной жидкости в рамках подготовки дипломной работы в Новосибирском государственном университете. Е.А. Гондюл, работая в ИМ СО РАН, провела
численные эксперименты с использованием оборудования суперкомпьютерного центра «Политехнический» на базе СПбПУ при финансовой поддержке РНФ грант № 21-71-20003. Постановка задачи проведена В.В. Лисицей.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Cahn J.W., Hilliard J.E. Free energy of non-uniform system I: interfacial free energy // J. Chem. Phys- 1958.-28. - P.258-267.
2. Cahn J.W., Hilliard J.E. Free energy of a nonuniform system III: nucleation in a two-component incompressible fluid // J. Chem. Phys.— 1959.—31.— P.688-699.
3. Sussman M., Smereka P., Osher S. A level set approach for computing solutions to incompressible two-phase flow // J. Comput. Phys. - 1994. - 114. - P. 146-159.
4. Gueyffier D., Li J., Nasim A., Scardovelli R., Zaleski S. Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows // J. Comput. Phys. - 1999. - 152. - P. 423-456.
5. Borcia R., Bestehorn M. Phase-field for Marangoni convection in liquid-gas systems with a deformable interface // Phys. Rev. E. - 2003. - 67. - 066307.
6. Brackbill J.U., Kothe D.B., Zemach C. A continuum method for modeling surface tension // Journal of computational physics. - 1991. - 100. - P.335-354.
7. Kim J.S. A continuous surface tension force formulation for diffuse-interface models // Journal of computational physics. - 2005. - 204(2). - P.784-804.
8. Liu C., Shen J. A phase field model for the mixture of two incompressible fluids and its approximation by a Fourier-spectral method // Phys. D. - 2003. - 179. - P. 211-228.
9. Chorin J. Numerical solution of the Navier-Stokes equations // Mathematics of Computation.-1968.-22.-104. -P.745-762
10. Темам Р. Уравнение Навье-Стокса. Теория и численный анализ // — М.: Мир. - 1981.— с. 408.
11. Misrandi H., Rajkotwala A.H., Balrussen M.W., Peters E.A. Numerical simulation of bubble formation with a moving contact line using Local Front Reconstruction Method // Chemical Engineering Science. -2018. - 187. - P.415-431.
12. Kim J. Phase-field models for multi-component fluid flowa // Communications in Computational Physics. - 2012.- 12(3).
13. Muradoglu M., Tasoglu S. A front-tracking method for computational modeling of impact and spreading of viscous droplets on solid walls // Computers & Fluids. - 2010. - 39. - P.261-625.
© Е. А. Гондюл, В. В. Лисица, 2022
