MODELLING OF PROPAGATION OF STRONG PERTURBATIONS IN GAS OF ARBITRARY RAREFACTION

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

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

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

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

Состояние любого газа можно описать одночастичной функцией распределения

где 

J – интеграл столкновений – выражение, учитывающее влияние столкновений частиц или молекул.

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

В кинетической теории газов одной из популярных моделей столкновений для уравнения Больцмана является модель Бхатнагара-Гросса-Крука (БГК), которая заключается в том, что интеграл столкновений имеет следующий вид

где 

где 

Для модели БГК локально-равновесная функция распределения 

где n – числовая плотность или концентрация частиц, 1/м3;

m – масса молекулы газа, кг;

M – молярная масса газа, кг/моль;

R = 8,31 Дж/(моль*K) – универсальная газовая постоянная;

T – температура газа, К

u – вектор средней скорости газового потока.

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

где 

Введем масштабные параметры:

В формулах (1.6) характеристическая скорость 

где 

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

Также обезразмерим давление и выразим ее через числовую плотность и температуру, записав уравнение идеального газа:

После подстановки получим безразмерное уравнение Больцмана:

где 

Скорость частицы 

где 

Для простоты записи формул вместо

где 

Также обратим внимание на то, что уравнение (1.13) является неявным из-за включения члена столкновений 

Подставив преобразование (1.14) в уравнение (1.13), получим:

В уравнении (1.15) 

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

где 

Подставим (1.16) в (1.15), получим следующее выражение: 

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

1) При 

2) Подставив (1.18) в (1.17), получим так называемый шаг столкновений:

где 

Таким образом, чтобы найти функцию возмущения 

В приведенном выше уравнении (1.20) ключевым моментом является вычисление градиента функции возмущения 

где 

Как видно из уравнения (1.21), ключом к вычислению градиента 

Для вычисления локально-равновесной функции возмущения 

где 

Вычислив 

Теперь для удобства решения рассматриваемой задачи введем следующие обозначения:

Подставим (1.26) в уравнение (1.17), получим: 

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

где 

Уравнение (1.28) также можно записать в виде условия сходимости Куранта-Фридрихса-Леви, которое является необходимым условием сходимости при численном решении некоторых уравнений в частных производных. Оно возникает при численном анализе явных схем интегрирования по времени, когда они используются для решения рассматриваемой задачи. Как следствие, при численном решении нестационарных задач шаг по времени должен быть меньше определенного времени, в противном случае моделирование дает неправильные результаты из-за расходимости искомых величин

где 

Обезразмерим рассматриваемый интервал времени 

Было показано, что 

где 

Подставим (1.31) в (1.26), получим следующие выражения:

Теперь подставим (1.16) и (1.31) в (1.14) и получим окончательное выражение для поиска функции возмущения 

В итоге принципиальная схема метода дискретных скоростей показана на рисунке 1.

Рисунок 1 - Схема метода дискретных скоростей в двумерном случае

DOI:10.60797/IRJ.2024.143.150.1

Как видно из этого рисунка и уравнений (1.27) и (1.33), для поиска функции возмущения 

$Y\left({\mathbf{x}}_c,{\mathbf{\xi }}_{\alpha },t+\Delta t\right)$

 необходимо рассчитать только функцию возмущения 

$\bar{Y}\left({\mathbf{x}}_c-{\mathbf{\xi }}_{\alpha }\Delta t,{\mathbf{\xi }}_{\alpha },t\right)$

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

$Y_e\left({\mathbf{x}}_c,{\mathbf{\xi }}_{\alpha },t+\Delta t\right).$

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

${\bar{Y}}^{\ast }\left({\mathbf{x}}_c,{\mathbf{\xi }}_{\alpha },t+\Delta t\right)$

 с помощью потокового процесса (1.18) и переместить функции возмущения на следующий временной шаг с помощью шага столкновений (1.27). Таким образом, сделав обратное преобразование (1.33) и вычислив

$Y\left({\mathbf{x}}_c,{\mathbf{\xi }}_{\alpha },t+\Delta t\right)$

,  можно рассчитать искомые функции распределения 

$f\left({\mathbf{x}}_c,{\mathbf{\xi }}_{\alpha },t+\Delta t\right)$

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

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

Из формулировки функции распределения

характеризует среднее число частиц, находящихся в момент времени 

Из кинетической теории газов известно, что среднее значение любой функции

Тогда составляющие макроскопической средней скорости газового потока 

Аналогичное выражение получим для

Составляющей средней скорости газового потока по азимутальному направлению

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

Обезразмерим записанное выражение при помощи (1.6) и (1.7): 

Сократим, опустим индекс обезразмеривания и применим формулу (1.35). В итоге получим окончательное выражение для вычисления температуры газа:

Для вычисления функции возмущения 

Чтобы вычислить соответствующие частные производные, необходимо восстановить значения функции возмущения на границах ячейки. Для этого в данной работе используется метод MUSCL (Monotonic Upstream-centered Scheme for Conservations Laws) так называемая монотонная схема законов сохранения, ориентированная на восходящий поток

[10]

,

[11]

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

Рисунок 2 - Графическая интерпретация метода MUSCL

DOI:10.60797/IRJ.2024.143.150.2

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

Зная 

где 

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

Аналогичная формула для частной производной по радиусу получится, если молекулы движутся в расчетной области слева направо [10]. В этом случае мы находим правую границу ячейки 

Выражения для частных производных по высоте 

Перейдем к вычислению ограничителя Бергер 

Затем вычисляются индикатор местоположения 

Тогда ограничитель Бергер 

Таким образом, функция 

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

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

где 

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

где

Поскольку нижняя стенка непроницаемая, следовательно, плотность потока молекул газа через нее равна нулю. Запишем это, воспользовавшись выражением (1.36):

где 

Тогда выражение для числовой плотности 

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

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

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

Таким образом, задаются следующие кусочно-постоянные начальные условия

где 

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

Рисунок 3 - Начальное распределение температуры и концентрации по высоте цилиндра

DOI:10.60797/IRJ.2024.143.150.3

При моделировании, как отмечалось ранее, используется цилиндрическая система координат. При этом исследуемый интервал скоростей разделяется на 40 значений, сетка по радиусу – на 11 ячеек, по высоте цилиндра – на 101 ячейку, по азимуту – на 11 ячеек и по полярному углу в пространстве скоростей на 20 ячеек. Причем первая и последняя ячейки по высоте и последняя ячейка по радиусу представляют некоторую фиктивную расчетную область цилиндра, где задаются граничные условия.

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