doctoral candidate from 01.01.2019 until now
Modelling porous media finds its application in many industries and human life, namely from the oil industry and construction to washing dishes. The article considers the processes of filtrating single-phase liquid through a porous structure, the description of which is based on lattice Boltzmann method in two-dimensional space. The problem of taking into account the structure of a simulated porous medium when applying the algorithm of lattice Boltzmann equations has many solutions. Nevertheless, at the microlevel, large errors in approximating the pore boundaries can occur, due to the fact that the lattice sites are assumed to be absolutely liquid or absolutely impenetrable. The paper proposes a model in which the lattice nodes are considered to be partially permeable. Partial permeability in this case refers to having a proportion of solids in the lattice site. As a result, it is possible to derive an equation that obeys the laws of mass and momentum conservation, the validity of which has been proven empirically. Due to the equation symmetry, it can be easily used not only on a plane, but also in three-dimensional space.The new lattice Boltzmann method can be employed, for example, in modelling water filtration through sandy soils, which is used in calculating subsidence of the foundations
area, density, experiment, Boltzmann equation, computational experiment
Введение
Метод решетчатых уравнений Больцмана (LBM) широко используется при моделировании процессов течения жидкости и воздушных потоков, а гибкость задания граничных условий позволяет учитывать сложную структуру среды.
При исследовании сред со сложной геометрической структурой современное применение метода охватывает такие проблемы как моделирование реактивных транспортных процессов в нефтегазовых резервуарах [1]; расчет потока жидкости и теплопередачи в средах со сложной геометрией, например, металлической пены [2]; моделирование взаимодействия волн и пористых структур [3] и др.
Наряду с методом решетчатых уравнений Больцмана для моделирования фильтрации потока жидкости через пористую среду также используются уравнения Навье-Стокса (NSE) и модели поровой сети (PNM). Уравнения Навье-Стокса решаются на основе метода конечных элементов, конечных объемов или конечно-разностных схем.
В [4] предлагается алгоритм, использующий модифицированные уравнения Навье-Стокса для моделирования фильтрации жидкости в анизотропной пористой среде, в [5] проводится исследование зависимостей результатов, полученных решением уравнений Навье-Стокса и закона Дарси для пористой среды.
Модели поровой сети [6] представляют реальное пространство множеством сферических пор, соединенных цилиндрическими капиллярами. Гидродинамическая модель основывается на законе сохранения массы и на известных характеристиках течения Пуазёйля в цилиндрических капиллярах. К PNM, например, относятся алгоритм преобразования медиальных осей, максимального шара и сегментации по водоразделам.
Каждая из перечисленных групп моделей имеет свои преимущества и недостатки. Уравнение Навье-Стокса всегда дает решение для потока жидкости, находящегося в стационарном режиме, в то же время в LBM поток переходит в стационарный режим движения жидкости после большого числа итераций (10 000…20 000, [6]), поэтому LBM может работать медленнее, чем NSE.
С другой стороны, несомненным достоинством LBM является легкость создания параллельного кода, что объясняется тем, что в LBM используются явные схемы расчетов в связи с чем возможно распараллеливание программы в достаточно высокой степени.
Также отметим, что в сравнении с PNM LBM зачастую не способен проанализировать всю структуру порового объекта, а алгоритмы PNM не имеют подобных ограничений.
Методология
Общие сведения о методе решетчатых уравнений Больцмана. LBM представляет собой дискретный аналог уравнения Больцмана [7]. В соответствии с LBM молекулы могут находиться только в узлах решетки, причем число направлений из одного узла в другой ограничено. Обычно решетка разбивается на клетки, имеющие форму квадратов или кубов.
Будем использовать модель D2Q9, которая имеет 2 измерения и 9 возможных направлений скорости. На рис. 1 изображено расположение решетки в декартовой системе координат со скоростями ci, где i = 0,1…8 индекс направления и c0 = 0 соответствует неподвижной молекуле.
Рис. 1. Модель D2Q9
Fig. 1. Model D2Q9
Функции распределения f(x, u, t) описывают состояние системы каждого узла решетки. Функция распределения определяет [6], какая часть частиц расположена в окрестности точки x(x, y) в момент времени t с координатами от x до x+Δx, от y до y+Δy и диапазоном скоростей от до
.
Все методы LBM основаны на принципе разделения этапов столкновения молекул (collision step) и потоковой передаче (streaming step).
На первом этапе [8] происходит релаксация функции (x, t) в направлении равновесной функции распределения
(x, u, t).
Во время потоковой передачи каждое распределение частиц перемещается в соответствии со скоростями ci в соседние узлы.
Объединив этапы столкновения и распространения молекул получим формулу (1):
где – оператор столкновений.
Равновесная функция распределения [7] определяется по формуле (2):
где веса wi = 4/9 для частиц находящихся в покое, 1/9 при i = 1, 2, 3, 4 и 1/36 при i = 5, 6, 7, 8; cs – скорость звука, принимается равной 1/ ; ρ(x) – плотность жидкости в точке x; u – скорость потока жидкости в точке x.
Макроскопические параметры скорости и плотности жидкости могут быть получены вычислением нулевого и первого момента функции распределения, то есть:
Схемы столкновений с двумя параметрами релаксации (TRT). Самой широко применяемой схемой столкновений является оператор BGK в совокупности со стандартными граничными условиями упругого столкновения. Из-за того, что оператор использует лишь один параметр релаксации, [9] появляется вычислительная нестабильность и наблюдается наличие зависимости расположения границ от вязкости жидкости.
С другой стороны, в [9] было показано, что необходимо как минимум два параметра релаксации, чтобы предотвратить наличие нелинейных ошибок, обусловленных вязкостью, при применении граничных условий типа упругого столкновения.
В связи с этим, в основе разрабатываемой модели будет лежать оператор столкновений с двумя параметрами релаксации (TRT-оператор). TRT-оператор имеет [10] параметр ω+, который задает вязкость жидкости, и свободный параметр ω-. Стабильность и точность оператора обусловлена так называемым «магическим параметром» :
Существуют известные значения , которые обеспечивают определенную степень точности, например,
= 1/12 отвечает третьему порядку точности, а
=1/6 – четвертому, а
= 1/4 зачастую стабилизирует вычисления.
Модель TRT-оператора основана [10] на декомпозиции ожидаемого фазового пространства молекул на симметричные и ассиметричные части. В методе решетчатых уравнений Больцмана для каждой скорости ci существует и таким образом множество скоростей ci обладает симметрией. По аналогии можно переписать формулы для расчета плотностей распределений:
Тогда стандартная TRT модель принимает вид:
И кинематическая вязкость жидкости:
Граничные условия. Будем считать, что исследуемая область жидкости имеет прямоугольную форму. Тогда на левой и правой границе применим периодические граничные условия, в верхней части зададим вектор скорости, направленный в область жидкости (Zou and He boundary conditions, например, описаны в [11]), в нижней условия нулевого градиента по функции распределения.
Внутри исследуемой области в качестве граничных условий используется видоизмененный half-way bounce-back (описание классического алгоритма приведено, например, в [10]).
Введение внешних сил в уравнение Больцмана. Учтем воздействие внешних сил на молекулы жидкости, добавив элемент в исходное уравнение (1) [12]:
Чтобы определить значение рассмотрим схему определения внешних сил, основанную на LGA ( Lattice gas automata – метод, предшествующий LBE):
которая соблюдает законы сохранения массы и импульса:
Заметим, что в рамках представленной схемы, расчетные формулы равновесной функции распределения и скорости течения жидкости остаются неизменными.
Для удобства определим [13] переменную , задающую изменение функции распределения вследствие влияния внешней силы и рассчитываемую по формуле:
где – ускорение (обусловливающее фиктивную силу тяжести), имитирующая действие градиента давления в направлении оси OY, а действующая сила представлена выражением
.
Внесение изменений в этап распространения частиц. Этап потоковой передачи зачастую включает в себя граничные условия вида упругого столкновения (bounce-back).
Будем применять half-way bounce-back, учитывая частичную непроницаемость узлов.
Введем матрицу , содержащую информацию о проницаемости каждого узла решетки, в таком случае
– где x указывает позицию узла;
если узел абсолютно жидкий и
если жидкость не может оказаться в точке x (например, в случае если узел полностью состоит из частицы грунта).
Объединим формулы упругого столкновения с непроницаемым узлом и распространения частиц в соседние узлы.
Взаимодействие узлов охарактеризуем с позиции функции вероятности, имеющей направление «в» узел, т.е. с позиции узла с координатами , а не x.
Пусть в момент времени t частицы перемещаются из узла в
по направлению k (рис. 2).
Рис. 2. Момент времени t
Fig. 2. Moment of time t
Тогда во время в узел
вернется
часть частиц, двигающихся в направлении
, с другой стороны в момент
в
из
переместится
частиц (рис. 3, 4).
Рис. 3. Момент времени t
Fig. 3. Moment of time t
Рис. 4. Момент времени t+Δ t
Fig. 4. Moment of time t+Δt
Таким образом, формула видоизменяется следующим образом:
Сравнение с аналогичными уравнениями. Пористую структуру среды можно описать двумя способами:
1. Ввести промежуточный шаг porous medium step [7]:
Обозначим , рассчитанной для этапа столкновения молекул
Тогда porous medium step принимает вид:
или сгруппируем
где макропараметр n – коэффициент пористости среды.
2. Считать каждый узел решетки частично непроницаемым, как например, в граничных условиях вида partially saturated bounce-back [8].
Выражение (11) схоже по структуре с (12), однако коэффициенты учитываются для каждого узла, что характерно для partially saturated bounce-back. В то же время коэффициенты
не зависят от параметра релаксации и в уравнении не используется дополнительный оператор столкновения. Кроме того, в (11) применяется TRT-оператор, в то время как в partially saturated bounce-back BGK-оператор.
Проверим, выполняются ли законы сохранения массы и импульса для выведенного уравнения.
Чтобы удостовериться в их справедливости проверим эквивалентность формул для нахождения плотности жидкости и скорости потока [12, 13]:
Тогда достаточно убедиться, что:
Выполнение законов будем проверять эмпирическим [14] путем на 8000 итераций:
Найдем нормы Фробениуса для матриц и
на каждой из итераций (рис. 5 – 7).
Из рисунков следует, что нормы матриц не превышают значений 8e-15 для любой итерации, т.е. соблюдаются законы сохранения массы и импульса.
Рис. 5. График зависимости от номера итерации
Fig. 5. Graph of dependence from the iteration number
Рис. 6. График зависимости от номера итерации
Fig. 6. Graph of dependence from the iteration number
Рис. 7. График зависимости от номера итерации
Fig. 7. Graph of dependence from the iteration number
Покажем, что
(рис. 8 – 10).
Рис. 8. График зависимости от номера итерации n
Fig. 8. Graph of the dependence of on the iteration number n
Рис. 9. График зависимости от номера итерации n
Fig. 8. Graph of the dependence of on the iteration number n
Рис. 10. График зависимости от номера итерации n
Fig.10. Graph of the dependence of on the iteration number n
Из рис. 10 следует, что , значит с некоторого номера N
,
.
Тогда , начиная с n = N преобразуется в
и при N
, следовательно,
, т.е. закон сохранения массы выполняется
начиная с номера N, из рис. 10 возьмем N = 8000. Законы сохранения импульса для
доказываются аналогично.
Заключение
Предложенный подход к описанию пористой структуры был проанализирован на предмет существования аналогичных методов. Были проведены теоретические и вычислительные эксперименты, в ходе которых удалось доказать справедливость законов сохранения массы и импульса.
Уравнение (11) можно применить при расчете поля скорости течения жидкости в пористых средах при низкой пористости среды и в случае малых размеров непроницаемых частиц, т.е. соизмеримых с 1 l.u. Если частицы значительно больше 1 l.u., то уравнение (11) фактически преобразуется в классическое с применением граничных условий типа half-way bounce-back, кроме того с точки зрения уменьшения вычислительных затрат целесообразно ввести промежуточный шаг porous medium step и перейти к уравнению с макропараметром n.
1. Jiang M., Xu Z.G., Zhou Z.P. Pore-Scale Investigation on Reactive Flow in Porous Media Considering Dissolution and Precipitation by LBM. Journal of Petroleum Science and Engineering. 2021;204:14579.
2. Hamidi E. Lattice Boltzmann Method Simula-tion of Flow and Forced Convective Heat Transfer on 3D Micro X-ray Tomography of Metal Foam Heat Sink. International Journal of Thermal Sciences. 2022;172(A):107240.
3. Xing E.A Three-Dimensional Model of Wave Interactions With Permeable Structures Using the Lattice Boltzmann Method. Applied Mathematical Modelling. 022;104:67-95.
4. Zhilenkov AA. High Productivity Numerical Computations for Gas Dynamics Modelling Based on DFT and Approximation. In: Proceedings of the 2018 IEEE Conference of Russian Young Researchers in Electrical and Electronic Engineering, ElConRus; 2018 Jan 29 - Feb 01; Saint Petersburg, Moscow: Institute of Electrical and electronics Engineers Inc.: 2018. p. 400-403. doi:https://doi.org/10.1109/EIConRus.2018.8317117
5. Narsilio G.A., Upscaling of Navier-Stokes Equations in Porous Media: Theoretical, Numerical and Experimental Approach. Computers and Geotechnics. 009;36(7):1200-1206.
6. Zhilenkov A.A. Developing a Method for Solv-ing Heat Conduction Equations With Non-Uniform Discretization for Modelling Processes in Gas-Phase Epitaxy Reactors. Control Systems and Information Technologies. 2017;3(69):11-15.
7. Sukop M.C. Lattice Boltzmann Modelling: An Introduction for Geoscientists and Engineers. Berlin Heidelberg: Springer; 2006.
8. Najuch T. Analysis of Two Partially-Saturated-Cell Methods for Lattice Boltzmann Simulation of Granular Suspension Rheology. Computers & Fluids. 2019;189:1-12.
9. Gharibi F. Darcy and Inertial Fluid Flow Simulations in Porous Media Using the Non-Orthogonal Central Moments Lattice Boltzmann Method. Journal of Petroleum Science and Engineering. 2020;194.
10. Kruger T. The Lattice Boltzmann Method, Principles and Practice. Switzerland: Springer; 2017.
11. Zhilenkov A.A. Hybrid Solution of the Navier-Stokes Equations in Spaces of Analytic Functions Using Bilinear Forms and Green’s Functions. Control Systems and Information Technologies. 2018;1(71):4-7.
12. Guo Zh. Lattice Boltzmann Method and Its Applications in Engineering. World Scientific Publishing Company; 2013.
13. Sauro S. The Lattice Boltzmann Equation. United Kingdom: Oxford University Press; 2018.
14. Zhilenkov A.A., Cherny S.G. Extracting Information From Big Data Using Neural Network Architectures as Networks of Information Granule Associations. Proceedings of the Institute for Systems Analysis of the Russian Academy of Sciences. 2022;72(3):81-90. doi:https://doi.org/10.14357/20790279220308