Введение

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

Кинетическая модель идеального бесстолкновительного газа — динамическая система, состоящая из статистически большого числа частиц, которые не взаимодействуют друг с другом, но взаимодействуют с границами. Основным вычислительным методом такой модели является метод Монте-Карло (ММК). Для проведения вычислительных экспериментов необходимо разработать комплекс проблемно-ориентированных программ, позволяющий рассчитывать огромное число траекторий движения частиц, взаимодействующих с границами. В процессе моделирования такой задачи получается огромный объем числовых данных. Визуализация таких данных в виде графиков необходимых величин и изображений процессов течения газа в моделируемой области позволяет представить в наглядной форме и качественно оценить полученные результаты.

1. Описание математической модели фильтрации идеального газа

Математическая модель фильтрации идеального газа в пористой среде описывается следующей системой уравнений [5-9]:

(1)

где  — проницаемость среды;  — динамическая вязкость среды;  — средняя скорость молекул;  — средняя длина свободного пробега молекул.

Далее рассматривается одномерный случай при условии, что течение стационарно

(2)

После интегрирования (1) уравнения непрерывности приобретает следующий вид:

(3)

При подстановке 1-го и последнего уравнений из системы (1) в (3) с последующим интегрированием получается формула фильтрации идеального газа в пористой среде (рис. 1):

(4)

а) конфигурация 1 (K1)

б) конфигурация 2 (K2)

в) конфигурация 3 (K3)

Рис. 1. Вычислительная схема: красным цветом изображена пористая среда; зеленым цветом — границы элементарных областей; желтым цветом — частицы

где  и  – номера элементарных областей,  – расстояние между элементарными областями (длина фильтрующей части пористой среды).

Элементарная область — выделенная часть некоторой моделируемой области. Каждой такой подобласти в определенный момент времени может принадлежать некоторое количество частиц. С помощью локальных параметров этих частиц рассчитываются значения статистических оценок макроскопических величин.

2. Постановка вычислительной задачи

Пусть моделируемая область, имеющая форму прямоугольного параллелепипеда, не ограничена в положительном направлении оси X. Область условно разбита на 3 подобласти.

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

Таблица 1. Параметры конфигураций пористой среды

С другой стороны, часть подобласти, в которой движутся частицы, содержит 6 элементарных областей, в которых вычисляются статистические оценки макроскопических параметров. Границы элементарных областей выделены зеленым цветом. Нумерация элементарных областей осуществляется слева направо вдоль оси Х, начиная с нуля.

Число сфер для -й конфигурации пористой среды определяется по следующей формуле:

(5)

3. Описание комплекса программ

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

(6)

(7)

где  — пространственные координаты частицы;  — скорость частицы после -го взаимодействия с границей;  — временной шаг;  — число столкновений частицы с границами;  — функция взаимодействия частицы с границей;  — нормаль к поверхности -й границы в точке столкновения в момент времени  на -м шаге.

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

(8)

где ; .

В случае, когда граница — неподвижная сфера, формула, по которой вычисляется момент времени , имеет вид [10]:

(9)

где ; ; .

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

Рис. 2. Иллюстрация, поясняющая корректировку с помощью поправочного коэффициента : жирная линия — граница; линия — траектория движения частицы с учетом взаимодействия с границей и поправочным коэффициентом ; точечный пунктир — траектория движения частицы с учетом взаимодействия с границей без дополнительного поправочного коэффициента

В случае, когда частица отражается от границы по закону зеркального отражения, скорость частицы после -го столкновения с -й движущейся границей находится [11]

(10)

На рис. 3 представлена схема расчета траектории движения частиц с учетом их возможного взаимодействия с границами.

Рис. 3. Схема расчета траекторий движения частиц с учетом их возможного взаимодействия с границами

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

Известно, что операции деления и извлечения корня медленнее операций сложения (вычитания) и умножения. В циклах «медленные» операции, стоящие перед дополнительными условиями, если возможно, переносятся внутрь условия, изменяя также и условие дополнительными операциями сложения и умножения, что приводит к увеличению скорости расчетов. Для быстрого извлечения обратного квадратного корня, который используется при нормализации векторов, применяется решение, предложенное в [13].

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

Рис. 4. Ячейка кубической сетки пористой среды

На рис. 5 представлена общая схема работы комплекса программ [14].

Рис. 5. Общая схема работы комплекса программ

Под установлением начального состояния модели подразумевается инициализация параметров частиц с помощью генератора псевдослучайных чисел, границ и элементарных областей. Для каждой элементарной области накапливаются суммы параметров частиц (общая масса, сумма компонентов скоростей и квадратов компонентов скоростей), принадлежащих определенной области в момент времени . С помощью этих сумм вычисляются статистические оценки макроскопических параметров по формулам, представленным в [5]. Сохранение промежуточных данных включает сохранение фазового пространства и сумм параметров частиц, которые используются на этапе обработки данных для визуализации процессов и вычисления статистических оценок макроскопических параметров соответственно. Также последний этап (обработка данных) может включать построение графиков изменения статистических оценок макроскопических параметров с течением времени и вычисления максимальных и относительных погрешностей, полученных в результате сравнения с точным решением.

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

Известно, что подавляющем большинстве случаев для управления современными супервычислительными системами используется ОС Linux. Разработанный комплекс программ ориентирован для семейства таких операционных систем. Для согласования всех программ, каждая из которых имеет свое специальное назначение, используются shell-скрипты. Алгоритмы, рассчитывающие траектории движения частиц с учетом взаимодействия с границами и вычисляющие статистические оценки макроскопических параметров, были реализованы на языке Си, позволяющим относительно эффективно (по сравнению с большинством других языков, исключениями составляют языки Fortran и С++) задействовать вычислительные ресурсы современных классических процессоров. Язык Си позволяет абстрагироваться от аппаратной архитектуры процессоров, но учитывать их особенности с помощью директив предпроцессора и использования специальных опций (ключей) на этапе компиляции. В качестве генератора псевдослучайных чисел был использован вихрь Мерсена, реализация которого существует в библиотеке MKL компании Intel. Визуализация течения газа выполнена с помощью графической библиотеки OpenGL. Визуализация поля статистических оценок макроскопических параметров строится с помощью пакета прикладных программ TecPlot, который требует файлы с входными данными в специальном формате. Для построения графиков использованы скрипты, написанные на Python. Этот язык поддерживает динамическую типизацию и содержит огромное количество вспомогательных библиотек (например, matplotlib).

4. Моделирование и анализ результатов

Были проведены серии вычислительных экспериментов с различным количеством частиц для трех конфигураций пористой среды. На рис. 6 представлена визуализация течения газа в пористой среде при .

а) ;

б) ;

в) ;

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

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

Поскольку коэффициент  содержит некоторые величины, которые можно вычислить (11), тогда для каждой конфигурации достаточно найти следующее отношение .

(11)

Вычисление коэффициента  делится на следующие этапы:

1. Для каждой пары скоростей  и  с одинаковыми номерами  и  вычисляется значение коэффициента  по следующей формуле:

(12)

Формула (12) основана на МНК с временным интервалом , когда течение газа выходит на квазистационарный режим, при . Значения  представлены в таблице 2 в зависимости от  и  при различном количестве частиц.

Таблица 2

Значение коэффициента  в зависимости от количества частиц , конфигурации пористой среды и номеров элементарных областей ,

2. Полученные значения  усредняются для каждой конфигурации пористой среды. На рис. 7-10 представлены некоторые результаты, где значения коэффициентов  для каждой конфигурации пористой среды взяты из таблицы 2 при .

а) ;                                              б) ;

в) ;                                               г) ;

д) ;                                              е) ;

Рис. 7. Графики зависимости скорости течения  и скорости фильтрации  для  при различном количестве частиц  от времени

а) ;                                               б) ;

в) ;                                               г) ;

д) ;                                               е) ;

Рис. 8. Графики зависимости скорости течения  и скорости фильтрации , статистической оценки давления и температуры для  при  от времени

а) ;                                              б) ;

в) ;                                              г) ;

д) ;                                              е) ;

Рис. 9. Графики зависимости скорости течения  и скорости фильтрации , статистической оценки давления и температуры для  при  от времени

а) ;                                              б) ;

в) ;                                              г) ;

д) ;                                              е) ;

Рис. 10. Графики зависимости скорости течения  и скорости фильтрации , статистической оценки давления и температуры для  при  от времени

Поле статистических оценок макроскопических параметров строится разбиением области на подобласти, где каждая подобласть — элементарная область. На рис. 11 представлены распределения температуры и давления для трех конфигураций пористой среды, где размер сетки (число элементарных областей по каждому измерению для моделируемой части области) равен 128х32х32 при .

а)                                       б)

в)                                     г)

д)                                      е)

Рис. 11. Области распределения статистических оценок давления и температуры для различных конфигураций пористой среды при : левый столбец — поле давления; правый столбец — поле температуры

Оценка производительности разработанного комплекса программ производилась на процессоре Intel Xeon E5-2690 V2 с пиковой производительностью 240 GFlops DP для 10 ядер.

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

Таблица 3

Оценка производительности на Intel Xeon E5-2690 V2 на одно ядро в зависимости от конфигурации пористой среды и шага по времени

Заключение

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

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

1. Сопоставление графиков скорости фильтрации и скорости течения газа позволило качественно подтвердить, что существует зависимость этих величин при различных конфигурациях пористой среды.

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

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

Работа выполнена при поддержке РФФИ, гранты 16-29-15105, 18-01-00343.

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

1. Галкин В. А., Быковских Д. А., Гавриленко Т. В., Стулов П. А. Фильтрационная модель движения идеального газа в пористой среде // Вестник кибернетики. Электр. Журн. 2016. Т. 4 (24). С. 50-57.

2. Быковских Д. А., Галкин В. А. О вычислительном тесте для одной модели бесстолкновительного идеального газа // Вестник кибернетики. Электр. Журн. 2017. 3 (27). С. 100-120.

3. Cercignani C. Mathematical methods in kinetic theory. Milano: Politecnico di Milano, 1969.

4. Галкин В. А. Анализ математических моделей: системы законов сохранения, уравнения Больцмана и Смолуховского. М.: БИНОМ. Лаборатория знаний, 2011.

5. Леонтьев Н. Е. Основы теории фильтрации. М.: Изд-во ЦПИ при мех.-мат. фак-ту МГУ, 2009.

6. Маскет М. Течение однородных жидкостей в пористой среде. М.: НИЦ «Регулярная и хаотическая динамика», 2004.

7. Wu Y.-S., Pruess K., Persoff P. Steady and Transient Analytical Solutions for Gas Flow in Porous Media with Klinkenberg Effects. California: Lawrence Berkeley National Laboratory, 1996.

8. Трапезникова М. А., Белоцерковская М. С., Четверушкин Б. Н. Аналог кинетически-согласованных схем для моделирования задачи фильтрации // Матем. моделирование. 2002. Т.14, No 10. С. 69–76.

9. Schmidt B. E. Compressible flow through porous media with application to injection: internal report for Caltech hypersonics group FM. Pasadena: California Institute of Technology, 2014.

10. Циссарж В. В., Марусик Р. И. Математические методы компьютерной графики. К.: ФАКТ, 2004.

11. Гальперин Г. А., Чернов Н. И. Биллиарды и хаос. М.: Знание, 1991.

12. Bik A. J. The Software Vectorization Handbook: Applying Intel Multimedia Extensions for Maximum Performance. Intel Press, 2004.

13. Robertson M. A Brief History of InvSqrt. Department of Computer Science & Applied Statistics, 2012.

14. Bird G. A. Molecular gas dynamics and the direct simulation of gas flows. N.Y.: Oxford University Press, 1994.

---

An ideal gas flow modeling in porous medium by Monte Carlo method

Authors: D.A. Bykovskih^1,A, V. A.  Galkin^2,B

^A Surgut State University

^B System Research Institute, Russian Academy of Sciences

¹ ORCID: 0000-0002-5796-3786, dmitriy.bykovskih@gmail.com

² ORCID: 0000-0002-9721-4026, val-gal@yandex.ru

 

Abstract

The article is devoted to research of collisionless gas flow in porous medium. Darcy’s law description of an ideal gas in porous medium was presented. The computing experiments with different particle numbers and the porous medium configurations were conducted. The relation between the flow speed (a statistical assessment of the hydrodynamic particle velocity) and the filtration rate of an ideal collisionless gas has been identified at the steady-state process.

The aim of the created software is computing experiments on modern high performance systems and visualization of the particle distribution evolution. The description and the performance evaluation of the software are included in the article.

 

Keywords: ideal gas filtration, Darcy’s law, Monte Carlo method.

 

References

               

1. Galkin V. A., Bykovskikh D. A., Gavrilenko T. V., Stulov P. A. Filtratsionnaya model dvizheniya idealnogo gaza v poristoy srede // Vestnik kibernetiki. Elektr. Zhurn. 2016. 4 (24). pp. 50–57.

2. Bykovskikh D. A., Galkin V. A. O vychislitelnom teste dlya odnoy modeli besstolknovitelnogo idealnogo gaza // Vestnik kibernetiki. Elektr. Zhurn. 2017. 3 (27). pp. 100-120.

3. Cercignani C. Mathematical methods in kinetic theory. Milano: Politecnico di Milano, 1969.

4. Galkin V. A. Analiz matematicheskikh modeley: sistemy zakonov sokhraneniya, uravneniya Boltsmana i Smolukhovskogo . M.: BINOM. Laboratoriya znaniy, 2011.

5. Leontev N. E. Osnovy teorii filtratsii. M.: Izd-vo TSPI pri mekh.-mat. fak-tu MGU, 2009.

6. Masket M. Techenie odnorodnykh zhidkostey v poristoy srede. M.: NITS «Regulyarnaya i khaoticheskaya dinamika», 2004.

7. Wu Y.-S., Pruess K., Persoff P. Steady and Transient Analytical Solutions for Gas Flow in Porous Media with Klinkenberg Effects. California: Lawrence Berkeley National Laboratory, 1996.

8. Trapeznikova M. A., Belotserkovska M. S., Chetverushkin B. N. Analog of kinetically consistent schemes for simulation of a filtration problem // Matem. Mod. 2002. Vol. 14, No 10, pp.  69-76.

9. Schmidt B. E. Compressible flow through porous media with application to injection: internal report for Caltech hypersonics group FM. Pasadena: California Institute of Technology, 2014.

10. Tsissarzh V. V., Marusik R. I. Matematicheskie metody kompyuternoy grafiki. K.: FAKT, 2004.

11. Galperin G. A., Chernov N. I. Billiardy i khaos. M.: Znanie, 1991.

12. Bik A. J. The Software Vectorization Handbook: Applying Intel Multimedia Extensions for Maximum Performance. Intel Press, 2004.

13. Robertson M. A Brief History of InvSqrt. Department of Computer Science & Applied Statistics, 2012.

14. Bird G. A. Molecular gas dynamics and the direct simulation of gas flows. N.Y.: Oxford University Press, 1994.