1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 62, № 4, 2022

Журнал вычислительной математики и математической физики, 2022, T. 62, № 4, стр. 705-720

Дискретная модель уравнения Больцмана для девяти скоростей: решение в виде суммы вайлда и приложения к моделированию несжимаемых течений

О. В. Ильин *

ВЦ ФИЦ ИУ РАН
119133 Москва, ул. Вавилова, 40, Россия

* E-mail: oilyin@gmail.com

Поступила в редакцию 27.07.2021
После доработки 05.09.2021
Принята к публикации 16.12.2021

Полный текст (PDF)

Аннотация

Рассмотрена дискретная кинетическая модель уравнения Больцмана на плоскости с девятью скоростями. В пределе малых длин свободного пробега и малых потоковых скоростей модель описывает течения вязкой несжимаемой жидкости или газа. Полная дискретизация модели по пространственным и временным переменным, которая в частности необходима для численного решения, осуществляется с помощью усеченной суммы Вайлда. Показано, что построенная схема имеет второй порядок точности. В качестве приложения предлагаемого метода построены численные решения эталонных задач: вихри Тэйлора–Грина, течение в полости с подвижной границей. Проведено сравнение результатов моделирования с решениями на основе классической решеточной модели уравнения Больцмана с девятью скоростями. Библ. 43. Фиг. 6.

Ключевые слова: решеточные уравнения Больцмана, уравнение Бхатнагара–Гросса–Крука (БГК), уравнения вязкой жидкости.

1. ВВЕДЕНИЕ

Метод решеточных уравнений Больцмнана (РУБ) в настоящее время популярен для моделирования как академических, так и прикладных задач в области механики сплошной среды и разреженных течений [1]. Данный подход основан на кинетической модели БГК [2], которая дискретизирована специальным образом. В общем виде решеточная модель уравнения Больцмана записывается следующим образом:

${{D}_{i}}{{f}_{i}} = \frac{1}{\tau }(f_{i}^{{eq}} - {{f}_{i}}),\quad i = 1,\; \ldots ,\;N,$
где ${{f}_{i}}(t,{\mathbf{r}})$ – функция распределения (ФР), отвечающая частицам, движущимся со скоростью ${{{\mathbf{c}}}_{i}}$, $t$, ${\mathbf{r}}$ – временные и пространственные переменные, $\tau $ – характерное время релаксации системы, $f_{i}^{{eq}}$ – локальное распределение Максвелла, $N$ есть число дискретных скоростей, ${{D}_{i}} = \frac{\partial }{{\partial t}} + {{{\mathbf{c}}}_{i}} \cdot \frac{\partial }{{\partial {\mathbf{r}}}}$. Для наиболее популярных РУБ локальное состояние равновесия имеет вид полиномиальной функции от потоковой скорости, получаемой из распределения Максвелла разложением в ряд Тейлора по потоковой скорости с конечным числом членов. Отметим, что для данной модели взаимодействия частиц, т.е. столкновения, описываются феноменологическим образом, правая часть (интеграл столкновений) отражает стремление ФР к локальному равновесию со скоростью, пропорциональной отклонению от этого локального равновесия.

Другим способом дискретизации в пространстве скоростей является метод дискретных скоростей для уравнения Больцмана [3]–[5]. В данном подходе столкновения частиц рассматриваются явным образом, в общем виде модель записывается в виде:

(1)
${{D}_{i}}{{f}_{i}} = \sum\limits_{jkl = 1}^N {A_{{kl}}^{{ij}}} ({{f}_{k}}{{f}_{l}} - {{f}_{i}}{{f}_{j}}) \equiv {{I}_{i}}\left[ {{{f}_{1}},\; \ldots ,\;{{f}_{N}}} \right],\quad i = 1,\; \ldots ,\;N,$
где $A_{{kl}}^{{ij}} = A_{{ij}}^{{kl}} \geqslant 0$ есть величины, пропорциональные вероятности того, что частицы со скоростями ${{{\mathbf{c}}}_{i}}$, ${{{\mathbf{c}}}_{j}}$ после столкновения имеют скорости ${{{\mathbf{c}}}_{k}}$, ${{{\mathbf{c}}}_{l}}$. Важным свойством дискретных кинетических моделей уравнения Больцмана является то, что для всех уравнений вида (1) выполняется $H$-теорема, т.е.
$\frac{{dH}}{{dt}} \leqslant 0,\quad H(t) = \int d {\mathbf{r}}\sum\limits_i^N \,{{f}_{i}}\lg ({{f}_{i}}),$
здесь следует отметить, что для стандартных решеточных уравнений Больцмана $H$-теорема, вообще говоря, неверна [6]–[8], чтобы обеспечить рост энтропии (убывание $H$ функции) в системе, описываемой РУБ, требуется выбирать состояние локального равновесия таким образом, чтобы максвелловское распределение зависело от потоковой скорости существенно неполиномиальным образом [9]–[12].

На первый взгляд метод дискретных скоростей для уравнения Больцмана более перспективен, чем метод РУБ. Действительно, столкновения задаются явным образом аналогично уравнению Больцмана, $H$-теорема верна по построению. Однако в настоящее время метод дискретных скоростей для уравнения Больцмана менее популярен. Одна из главных проблем распространенных дискретных моделей уравнения Больцмана заключается в том, что они не обладают достаточно изотропным набором дискретных скоростей, а данное свойство необходимо для того, чтобы в пределе малых длин свободного пробега правильно описывалось поведение вязкой жидкости или газа [13]. Например, уравнение Броудэлла на плоскости, которое подробно исследовалось во многих работах [3]–[19], имеет всего четыре возможные скорости движения частиц в сплошносредном пределе не может правильно воспроизвести даже уравнения Эйлера. Также следует отметить еще одну особенность дискретных моделей уравнения Больцмана: гипотеза молекулярного хаоса не выполняется для таких моделей, например, для уравнения Бродуэлла на плоскости две сталкивающиеся частицы имеют существенно ненулевую вероятность столкнуться еще раз после одного дополнительного (промежуточного) столкновения, т.е. данные частицы оказываются скоррелированными перед столкновением [20]. С теоретической точки зрения данный эффект нежелателен, но его влияние в пределе сплошной среды на прикладные задачи, рассматриваемые в работе, незначительно. Дискретные модели с достаточным числом скоростей для описания уравнений гидродинамики были рассмотрены относительно недавно [21], [22]. Кроме того, была построена модель для девяти скоростей, правильно описывающая в пределе малых длин свободного пробега уравнения Навье–Стокса для случая медленных несжимаемых течений [23]. Приближенные решения данной модели рассмотрены в настоящей работе.

Для уравнения Больцммана в случае максвелловских псевдомолекул и для дискретных моделей уравнения Больцмана возможно построить решение в виде рядов Вайлда [24], [25], представляющие собой бесконечные суммы со свертками интеграла столкновений разных порядков. Данное свойство позволяет построить приближенные решения этих уравнений в виде конечных сумм Вайлда [26], [27]. В настоящей работе рассмотрена аппроксимация дискретной модели уравнения Больцмана с девятью скоростями с помощью варианта конечной суммы Вайлда. Для данной аппроксимации получены выражения вязкости в сплошносредном пределе и численно решены модельные задачи: двумерные вихри Тэйлора–Грина, течение в квадратной полости с твердыми границами, причем верхняя граница является подвижной. Проведено сравнение результатов моделирования с результатами расчетов на основе стандартной РУБ с девятью скоростями.

2. ДИСКРЕТНАЯ КИНЕТИЧЕСКАЯ МОДЕЛЬ НА ПЛОСКОСТИ С ДЕВЯТЬЮ СКОРОСТЯМИ

В настоящей работе удобно ввести следующие обозначения. Скалярное произведение двух векторов ${\mathbf{X}} \in {{\mathbb{R}}^{D}}$ и ${\mathbf{Y}} \in {{\mathbb{R}}^{D}}$ (здесь $D$ – размерность пространства) с координатами ${{X}_{\alpha }}$, ${{Y}_{\beta }}$; $\alpha ,\beta = 1,\; \ldots ,\;D$ обозначается как ${\mathbf{X}} \cdot {\mathbf{Y}} = \sum\nolimits_\alpha \,{{X}_{\alpha }}{{Y}_{\alpha }}$, кроме того, ${\mathbf{XY}}$ есть тензор второго ранга с компонентами ${{X}_{\alpha }}{{Y}_{\beta }}$. Если ${\mathbf{T}}$, ${\mathbf{S}}$ – тензоры второго ранга, то обозначает оператор свертки, т.е. ${\mathbf{T}}:{\mathbf{S}} = \sum\nolimits_{\alpha \beta } \,{{T}_{{\alpha \beta }}}{{S}_{{\alpha \beta }}}$. Также, если ${\mathbf{X}}$ вектор, а ${\mathbf{T}}$ – тензор второго ранга, то произведение ${\mathbf{X}} \cdot {\mathbf{T}}$ есть вектор ${\mathbf{Z}}$ с координатами ${{Z}_{\alpha }} = \sum\nolimits_\beta \,{{X}_{\beta }}{{S}_{{\alpha \beta }}}$.

Состояние локального равновесия для дискретных кинетических моделей (1) имеет вид (см. [29]):

(2)
$f_{i}^{{eq}} = \exp (a + {\mathbf{b}} \cdot {{{\mathbf{c}}}_{i}} + dc_{i}^{2}),\quad i = 1,\; \ldots ,\;N,$
где $a$, ${\mathbf{b}}$, $d$ – параметры, вычисляемые из законов сохранения
$\sum\limits_i {f_{i}^{{eq}}} = \rho ,\quad \sum\limits_i {f_{i}^{{eq}}} {{{\mathbf{c}}}_{i}} = \rho {\mathbf{u}},\quad \sum\limits_i {f_{i}^{{eq}}} c_{i}^{2} = \rho \left( {{{u}^{2}} + D\theta } \right),$
здесь $\rho $, ${\mathbf{u}}$, $\theta $ – плотность, потоковая скорость, температура, $N$ – число дискретных скоростей. В случае ${\mathbf{u}} = 0,$ а также постоянных плотности и температуры имеем состояние абсолютного равновесия, которое обозначим через ${{w}_{i}}$, т.е.
${{w}_{i}} = \exp \left( {{{a}_{0}} + {{c}_{0}}c_{i}^{2}} \right),$
где ${{a}_{0}}$, ${{d}_{0}}$ – числовые параметры, по аналогии с теорией РУБ будем называть величины ${{w}_{i}}$, соответствующие дискретной скорости ${{{\mathbf{c}}}_{i}}$, весами. Из законов сохранения получаем следующие соотношения для весов:
(3)
$\sum\limits_i {{{w}_{i}}} = 1,\quad \sum\limits_i {{{w}_{i}}} {{{\mathbf{c}}}_{i}} = 0,\quad \sum\limits_i {{{w}_{i}}} {{{\mathbf{c}}}_{i}}{{{\mathbf{c}}}_{i}} = {{\theta }_{0}}\delta ,$
где ${{{\mathbf{c}}}_{i}}{{{\mathbf{c}}}_{i}}$ есть тензор с элементами ${{c}_{{i,\alpha }}}{{c}_{{i,\beta }}}$, $\alpha ,\beta = 1,\; \ldots ,\;D$, где $D$ – число пространственных переменных.

Рассматриваемая дискретная модель уравнения Больцмана для двух пространственных переменных состоит из скоростей трех типов [23] (см. фиг. 1): нулевая скорость ${{{\mathbf{c}}}_{0}} = (0,\;0)$ с весом ${{w}_{0}} = 16{\text{/}}36$; четыре скорости, параллельные осям $x$, $y$ , т.е. ${{{\mathbf{c}}}_{{ \pm 1}}} = ( \pm 1,\;0)c$, ${{{\mathbf{c}}}_{{ \pm 2}}} = (0,\; \pm 1)c$ с весом ${{w}_{0}} = 4{\text{/}}36$; четыре скорости, направленные по диагоналям ${{{\mathbf{c}}}_{{ \pm 3}}} = ( \pm 1,\; \pm 1)c$, ${{{\mathbf{c}}}_{{ \pm 4}}} = ( \pm 1,\; \mp 1)c$ с весами ${{w}_{0}} = 1{\text{/}}36$, где $c$ – положительная величина. Также данная модель описывает среду с температурой ${{\theta }_{0}} = \sum\nolimits_i {{{w}_{i}}} c_{i}^{2} = {{c}^{2}}{\text{/}}3$. Очевидно, что набор дискретных скоростей для данной модели совпадает с дискретными скоростями решеточной модели Больцмана $D2Q9$ (см. [1]).

Фиг. 1.

Дискретные скорости для рассматриваемой модели.

Девятискоростная дискретная модель уравнения Больцмана модель записывается следующим образом [23]:

(4)
$\frac{{\partial {{f}_{1}}}}{{\partial t}} + c\frac{{\partial {{f}_{1}}}}{{\partial x}} = \alpha {{J}_{0}} + \beta ({{J}_{1}} + {{J}_{2}}) + \lambda ({{J}_{6}} + {{J}_{7}}),$
(5)
$\frac{{\partial {{f}_{{ - 1}}}}}{{\partial t}} - c\frac{{\partial {{f}_{{ - 1}}}}}{{\partial x}} = \alpha {{J}_{0}} + \beta ({{J}_{3}} + {{J}_{4}}) - \lambda ({{J}_{6}} + {{J}_{7}}),$
(6)
$\frac{{\partial {{f}_{2}}}}{{\partial t}} + c\frac{{\partial {{f}_{2}}}}{{\partial y}} = - \alpha {{J}_{0}} + \beta ({{J}_{1}} + {{J}_{3}}) + \lambda ({{J}_{8}} + {{J}_{9}}),$
(7)
$\frac{{\partial {{f}_{{ - 2}}}}}{{\partial t}} - c\frac{{\partial {{f}_{{ - 2}}}}}{{\partial y}} = - \alpha {{J}_{0}} + \beta ({{J}_{2}} + {{J}_{4}}) - \lambda ({{J}_{8}} + {{J}_{9}}),$
(8)
$\frac{{\partial {{f}_{3}}}}{{\partial t}} + c\frac{{\partial {{f}_{3}}}}{{\partial x}} + c\frac{{\partial {{f}_{3}}}}{{\partial y}} = \gamma {{J}_{5}} - \beta {{J}_{1}} - \lambda ({{J}_{6}} + {{J}_{9}}),$
(9)
$\frac{{\partial {{f}_{{ - 3}}}}}{{\partial t}} - c\frac{{\partial {{f}_{{ - 3}}}}}{{\partial x}} - c\frac{{\partial {{f}_{{ - 3}}}}}{{\partial y}} = \gamma {{J}_{5}} - \beta {{J}_{4}} + \lambda ({{J}_{7}} + {{J}_{8}}),$
(10)
$\frac{{\partial {{f}_{4}}}}{{\partial t}} + c\frac{{\partial {{f}_{4}}}}{{\partial x}} - c\frac{{\partial {{f}_{4}}}}{{\partial y}} = - \gamma {{J}_{5}} - \beta {{J}_{2}} + \lambda ( - {{J}_{7}} + {{J}_{9}}),$
(11)
$\frac{{\partial {{f}_{{ - 4}}}}}{{\partial t}} - c\frac{{\partial {{f}_{{ - 4}}}}}{{\partial x}} + c\frac{{\partial {{f}_{{ - 4}}}}}{{\partial y}} = - \gamma {{J}_{5}} - \beta {{J}_{3}} + \lambda ({{J}_{6}} - {{J}_{8}}),$
(12)
$\frac{{\partial {{f}_{0}}}}{{\partial t}} = - \beta ({{J}_{1}} + {{J}_{2}} + {{J}_{3}} + {{J}_{4}}),$
где $\alpha $, $\beta $, $\gamma $, $\lambda $ неотрицательны. Кроме того,
(13)
${{J}_{0}} = {{f}_{{ - 2}}}{{f}_{2}} - {{f}_{{ - 1}}}{{f}_{1}},$
отвечает за взаимодействия частиц $1$ и $ - 1$, которые после столкновения переходят в состояния $2$ и $ - 2$, схематично это взаимодействие обозначим как $(1,\; - {\kern 1pt} 1) \to (2,\; - {\kern 1pt} 2)$.

Также имеются столкновения, которые связывают частицы всех трех типов (т.е. с разными величинами кинетической энергии), имеется четыре реакции: $(1,\;2) \to (0,\;3)$, $(1,\; - {\kern 1pt} 2) \to (0,\;4)$, $( - {\kern 1pt} 1,\; - {\kern 1pt} 2) \to (0,\; - {\kern 1pt} 3)$, $( - 1,\;2) \to (0,\; - {\kern 1pt} 4)$, т.е.

(14)
$\begin{gathered} {{J}_{1}} = {{f}_{0}}{{f}_{3}} - {{f}_{1}}{{f}_{2}},\quad {{J}_{2}} = {{f}_{0}}{{f}_{4}} - {{f}_{1}}{{f}_{{ - 2}}}, \\ {{J}_{3}} = {{f}_{0}}{{f}_{{ - 4}}} - {{f}_{{ - 1}}}{{f}_{2}},\quad {{J}_{4}} = {{f}_{0}}{{f}_{{ - 3}}} - {{f}_{{ - 1}}}{{f}_{{ - 2}}}. \\ \end{gathered} $

Cтолкновения между частицами, направленными по диагоналям, т.е. $(3,\; - {\kern 1pt} 3) \to (4,\; - {\kern 1pt} 4)$, они задаются слагаемым вида

(15)
${{J}_{5}} = {{f}_{{ - 4}}}{{f}_{4}} - {{f}_{{ - 3}}}{{f}_{3}}.$

Cтолкновения между частицами, направленными по диагоналям, и частицами, движущимися параллельно осям, т.е. $( - 4,\;1) \to ( - 1,\;3)$, $( - 3,\;1) \to ( - 1,\;4)$, $( - 3,\;2) \to ( - 4,\; - {\kern 1pt} 2)$, $(4,\;2) \to ( - 2,\;3)$, получаем слагаемые вида

(16)
$\begin{gathered} {{J}_{6}} = {{f}_{3}}{{f}_{{ - 1}}} - {{f}_{{ - 4}}}{{f}_{1}},\quad {{J}_{7}} = {{f}_{{ - 1}}}{{f}_{4}} - {{f}_{{ - 3}}}{{f}_{1}}, \\ {{J}_{8}} = {{f}_{{ - 4}}}{{f}_{{ - 2}}} - {{f}_{{ - 3}}}{{f}_{2}},\quad {{J}_{9}} = {{f}_{3}}{{f}_{{ - 2}}} - {{f}_{4}}{{f}_{2}}. \\ \end{gathered} $

Модель (4)–(16) обладает следующим важным свойством. Разложение Чепмена-Энскога [1] для данной модели дает уравнения вязкой жидкости, которые отличаются от уравнений Навье–Стокса для несжимаемой среды только членами вида $O\left( {{{{\text{M}}}^{3}}} \right)$, где ${\text{M}}$ есть число Маха, определяемое как $u{\text{/}}\sqrt {{{\theta }_{0}}} $, если выполняется равенство [23]

(17)
$4\alpha = \gamma + 4\beta + 4\lambda ,$
причем вязкость $\nu $ равна

(18)
$\nu = \frac{3}{{4\alpha }}.$

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

3. АППРОКСИМАЦИЯ РЕШЕНИЙ УСЕЧЕННОЙ СУММОЙ ВАЙЛДА

Следуя предыдущим результатам [26], [27], рассмотрим обыкновенное дифференциальное уравнение вида

(19)
$\frac{{df}}{{dt}} = \frac{1}{\epsilon }(P(f,f) - \lambda f),\quad {{\left. f \right|}_{{t = 0}}} = {{f}_{0}},$
где $P$ – билинейный оператор, $\epsilon $, $\lambda $ – положительные константы. Уравнение (19) имеет решение, которое можно записать в следующем виде (см. [26], [27]):
(20)
$f = \exp \left( { - \frac{{\lambda t}}{\epsilon }} \right)\sum\limits_{n = 0}^\infty {} (1 - \exp ( - \lambda t{\text{/}}\epsilon )){{f}^{{(n)}}},\quad {{f}^{{(n + 1)}}} \equiv \frac{1}{{(n + 1)\lambda }}\sum\limits_{m = 0}^n P \left( {{{f}^{{(m)}}},{{f}^{{(n - m)}}}} \right),$
решение в форме (20) называется суммой Вайлда (см. [24]).

Теперь снова рассмотрим дискретную кинетическую модель (1), но в пространственно однородном случае перепишем ее следующим образом:

(21)
$\frac{{d{{f}_{i}}}}{{dt}} = {{R}_{i}}(f,f) - A\rho {{f}_{i}},$
где
(22)
${{R}_{i}}(f,f) = \sum\limits_{jkl = 1}^N {A_{{kl}}^{{ij}}} {{f}_{k}}{{f}_{l}} + \sum\limits_{j = 1}^N {\left( {A - \sum\limits_{kl = 1}^N {A_{{kl}}^{{ij}}} } \right)} {{f}_{i}}{{f}_{j}},$
также константа $A$ выбрана так
(23)
$A = \mathop {\max }\limits_{ij} \sum\limits_{kl} {A_{{kl}}^{{ij}}} ,$
отметим, что выбор $A$ в виде (23) гарантирует, что оператор ${{R}_{i}}$ в (21), (22) неотрицателен. Допустим, что нам известно решение системы (21), (22) в момент времени $t$. Решение в виде суммы Вайлда для (21), (22) в момент $t + \Delta t$ (где $\Delta t$ мало) по аналогии с [26], [27] аппроксимируем конечной суммой Вайлда
(24)
$\begin{gathered} {{f}_{i}}(t + \Delta t) = \exp \left( { - A\rho (t)\Delta t} \right)\sum\limits_{n = 0}^M {{{{(1 - \exp \left( { - A\rho (t)\Delta t} \right))}}^{n}}} f_{i}^{{(n)}}(t) + \\ + \;{{(1 - \exp \left( { - A\rho (t)\Delta t} \right))}^{{M + 1}}}f_{i}^{{(\infty )}},\quad i = 1,\; \ldots ,\;N, \\ \end{gathered} $
(25)
$f_{i}^{{(n)}}(t) = \frac{1}{{A\rho (t)(n + 1)}}\sum\limits_{m = 0}^n {{{R}_{i}}} \left( {{{f}^{{(m)}}}(t),{{f}^{{(n - m)}}}(t)} \right),\quad i = 1,\; \ldots ,\;N,$
где $M \geqslant 1$, а также $f_{i}^{{(\infty )}} = \mathop {\lim }\nolimits_{n \to \infty } f_{i}^{{(n)}}$. Отметим, что (24) выводится из бесконечной суммы Вайлда, если $f_{i}^{{(M + k)}} = f_{i}^{{(\infty )}}$, $k = 1,\;2,\; \ldots $ . Возникает вопрос, как вычислить величину $f_{i}^{{(\infty )}}$, с физической точки зрения данная ФР соответствует частице, испытавшей бесконечное число столкновений, следовательно, $f_{i}^{{(\infty )}}$ близко к локальному максвелловскому распределению, т.е. можно положить
$f_{i}^{{(\infty )}} = f_{i}^{{(eq)}},\quad i = 1,\; \ldots ,\;N.$
Можно показать, что схема (24), (25) в пространственно однородном случае сходится к точному решению дискретной кинетической модели, при этом ошибка монотонно уменьшается с ростом $M$ [26], [27]. Для доказательства используется свойство неотрицательности оператора ${{R}_{i}}$, которое гарантируется равенством (23).

В текущей работе основной интерес представляет вопрос описания гидродинамики с помощью дискретной кинетической модели, а не рост порядка аппроксимации суммой Вайлда кинетической модели, таким образом, величину $A$ будем выбирать не с помощью условия (23), а из физических соображений. В дальнейшем с помощью разложения Чепмена–Энскога будет показано, что данный подход позволяет получить уравнения вязкой несжимаемой жидкости. C физической точки зрения $\rho (t)A$ характеризует величину, обратно пропорциональную времени релаксации системы $\tau $, которое, в свою очередь, пропорционально времени свободного пробега, последнее может быть выражено из вязкости рассматриваемой жидкости или газа. Имеем

(26)
$\frac{1}{{A\rho (t)}} = \tau = \frac{\nu }{{{{\theta }_{0}}}},$
так как моделируется несжимаемый газ, то $\rho (t) = {{\rho }_{0}}$, т.е. постоянная величина, соответственно в левой и правой части (26) стоят выражения, не зависящие от времени.

Тогда, учитывая все изложенное выше, решение системы (1) аппроксимируем суммой Вайлда для $M = 1$, при этом адвекцию учтем аналогично РУБ (см. [1]), в итоге в пространственно неоднородном случае имеем следующую схему:

(27)
$\begin{gathered} {{f}_{i}}(t + \Delta t,{\mathbf{r}} + {{{\mathbf{c}}}_{i}}\Delta t) = \exp \left( { - \frac{{\Delta t}}{\tau }} \right){{f}_{i}}(t,{\mathbf{r}}) + \\ + \;\exp \left( { - \frac{{\Delta t}}{\tau }} \right)\left( {1 - \exp \left( { - \frac{{\Delta t}}{\tau }} \right)} \right)f_{i}^{{(1)}}(t,{\mathbf{r}}) + {{\left( {1 - \exp \left( { - \frac{{\Delta t}}{\tau }} \right)} \right)}^{2}}f_{i}^{{eq}}(t,{\mathbf{r}}), \\ \end{gathered} $
(28)
$f_{i}^{{(1)}}(t,{\mathbf{r}}) = \tau {{R}_{i}}(f,f)(t,{\mathbf{r}}) = \tau \sum\limits_{jkl = 1}^N {A_{{kl}}^{{ij}}} ({{f}_{k}}{{f}_{l}} - {{f}_{i}}{{f}_{j}})(t,{\mathbf{r}}) + {{f}_{i}}(t,{\mathbf{r}}),$
где время релаксации $\tau $ выражается через вязкость $\nu $ и температуру ${{\theta }_{0}}$

(29)
$\tau = \frac{\nu }{{{{\theta }_{0}}}}.$

4. ГИДРОДИНАМИЧЕСКИЙ ПРЕДЕЛ ДЛЯ СЛУЧАЯ МЕДЛЕННЫХ ТЕЧЕНИЙ И ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ

Уравнения гидродинамики, которые получаются из системы (27)–(29), в случае, когда столкновения и дискретные скорости задаются девятискоростной моделью (4)–(12), могут быть получены с помощью разложения Чепмена–Энскога. Детали разложения представлены в Приложении ниже. Схема (27)–(29) аппроксимирует уравнения Навье–Стокса (для несжимаемых течений) со вторым порядком точности по пространственным и временной переменным, причем вязкость выражается формулой

(30)
$\nu = \frac{{{{\theta }_{0}}\Delta t}}{{{{{\left( {1 - {{e}^{{ - \frac{{\Delta t}}{\tau }}}}} \right)}}^{2}} + \frac{4}{9}\tau {{\rho }_{0}}{{e}^{{ - \frac{{\Delta t}}{\tau }}}}\left( {1 - {{e}^{{ - \frac{{\Delta t}}{\tau }}}}} \right)\alpha }} - \frac{{{{\theta }_{0}}\Delta t}}{2},\quad \tau = \nu {\text{/}}{{\theta }_{0}},$
также дополнительно требуется выполнение условия

$4\alpha = \gamma + 4\beta + 4\lambda .$

Для решения тестовых задач используется схема (27)–(29) , при этом необходимо определить параметры $\alpha $, $\beta $, $\gamma $, $\lambda $. При решении модельных задач обычно известна вязкость или число Рейнольдса (из которого можно выразить вязкость). Тогда, используя равенство (30), находим $\alpha $:

$\alpha = \frac{1}{{\frac{4}{9}\tau {{\rho }_{0}}{{e}^{{ - \frac{{\Delta t}}{\tau }}}}(1 - {{e}^{{ - \frac{{\Delta t}}{\tau }}}})}}\left( {\frac{{{{\theta }_{0}}\Delta t}}{{\nu + \frac{{{{\theta }_{0}}\Delta t}}{2}}} - {{{(1 - {{e}^{{ - \frac{{\Delta t}}{\tau }}}})}}^{2}}} \right),\quad \tau = \nu {\text{/}}{{\theta }_{0}},$
здесь следует отметить, что для малых вязкостей экспоненциальные члены, стоящие в знаменателе, приводят к очень большим величинам $\alpha $, поэтому на практике удобно использовать аппроксимацию ${{e}^{{ - \frac{{\Delta t}}{\tau }}}} \approx {1 \mathord{\left/ {\vphantom {1 {\left( {1 + \frac{{\Delta t}}{\tau }} \right)}}} \right. \kern-0em} {\left( {1 + \frac{{\Delta t}}{\tau }} \right)}}$, в этом случае
(31)
$\alpha = \frac{{{{{\left( {1 + \frac{{\Delta t}}{\tau }} \right)}}^{2}}}}{{\frac{4}{9}{{\rho }_{0}}\Delta t}}\left( {\frac{{{{\theta }_{0}}\Delta t}}{{\nu + \frac{{{{\theta }_{0}}\Delta t}}{2}}} - {{{\left( {1 - \frac{1}{{1 + \frac{{\Delta t}}{\tau }}}} \right)}}^{2}}} \right),\quad \tau = \nu {\text{/}}{{\theta }_{0}},$
параметры $\alpha $, $\beta $, $\gamma $, $\lambda $ определяются из (17), в текущей работе положим $\lambda = 0$, т.е. столкновения между частицами, направленными по диагоналям, и частицами, движущимися параллельно осям, исключаются, это несколько упрощает систему. Кроме того, положим $\beta = 0.1\alpha $ и тогда $\gamma = 4\alpha - 4\beta = (18{\text{/}}5)\alpha $.

Схема решается в два шага. На первом шаге вычисляется столкновительная часть (правая часть) в (27)–(29), затем проводится шаг адвекции – ФР, полученные на первом шаге, разносятся по узлам ${\mathbf{x}} + {{{\mathbf{c}}}_{{\mathbf{i}}}}\Delta t$ в соответствии со скоростями ${{{\mathbf{c}}}_{i}}$. Данный подход соверешенно аналогичен применяемому в теории РУБ [1].

Также из (27)–(29) видно, что для решения необходимо знать локально равновесную функцию распределения (2), коэффициенты которой выражаются из законов сохранения. Уравнения для коэффициентов существенно нелинейные, но в текущей работе рассматриваются медленные изотермические сплошносредные течения, поэтому локально равновесное распределение можно разложить в ряд Тейлора, считая, что отклонения коэффициентов $a$, $b$, $d$ от их значений в случае абсолютно максвелловского распределения малы [23]. В случае медленных изотермических течений локальная максвелловская ФР имеет тот же вид, как и для стандартной РУБ $D2Q9$ (см. [1], [23])

$f_{i}^{{eq}} = \rho {{w}_{i}}\left( {1 + \frac{{{{{\mathbf{c}}}_{i}}{\mathbf{u}}}}{{{{\theta }_{0}}}} + \frac{{{\mathbf{uu}}}}{{2\theta _{0}^{2}}}:({{{\mathbf{c}}}_{i}}{{{\mathbf{c}}}_{i}} - {\mathbf{I}}{{\theta }_{0}})} \right),$
очевидно, что локальная максвелловская ФР в данной форме может быть легко вычислена на каждом временном шаге в любой пространственной точке, так как она выражается через потоковую скорость и плотность.

Отметим, что для всех тестовых задач брались следующие значения температуры и скорости частиц: ${{\theta }_{0}} = 1{\text{/}}3$, $c = 1$, кроме того, временной шаг $\Delta t = 1$.

4.1. Вихри Тэйлора–Грина

Рассмотрим область в виде квадрата размером $L \times L$, заполненную несжимаемой жидкостью или газом, в которой начальные скорости заданы следующим образом:

${{u}_{x}}(x,y,t = 0) = - {{U}_{0}}\cos (kx)\sin (ky),\quad {{u}_{y}}(x,y,t = 0) = {{U}_{0}}\sin (kx)\cos (ky),$
где ${{U}_{0}} = 0.01$ есть максимальная амплитуда скорости, вязкость $\nu = 0.01$, временной шаг $\Delta t = 1$, также $k = \frac{{2\pi }}{L}$. На границах ставятся периодические граничные условия. Аналитическое решение задачи следующее:
${{u}_{x}}(x,y,t) = - {{U}_{0}}\cos (kx)\sin (ky){{e}^{{ - 2\nu {{k}^{2}}t}}},\quad {{u}_{y}}(x,y,t) = {{U}_{0}}\sin (kx)\cos (ky){{e}^{{ - 2\nu {{k}^{2}}t}}},$
из данного решения видно, что вихревая структура сохраняется во все моменты времени, но уменьшается амплитуда скоростей. Для численного решения будем измерять ошибку ${\text{error}}$ относительно аналитического решения по формуле
(32)
${\text{error}} = \frac{{\sqrt {\sum\limits_{x,y} {{{{\left( {u_{x}^{{(n)}}(T,x,y) - {{u}_{x}}(T,x,y)} \right)}}^{2}}} + {{{\left( {u_{y}^{{(n)}}(T,x,y) - {{u}_{y}}(T,x,y)} \right)}}^{2}}} }}{{\sqrt {\sum\limits_{x,y} {{{u}_{x}}} {{{(T,x,y)}}^{2}} + {{u}_{y}}{{{(T,x,y)}}^{2}}} }},$
где суммирование ведется по всем пространственным узлам схемы $x$, $y$, также $u_{x}^{{(n)}}$, $u_{y}^{{(n)}}$ есть численные решения, $T = 1{\text{/}}\left( {2{{k}^{2}}\nu } \right)$. Число пространственных узлов $N$ в каждом направлении $x$ или $y$ выбиралось равным $K = 15,\;25,\;49,\;73$. Линии тока для разрешения $25 \times 25$ в момент времени $1{\text{/}}\left( {2{{k}^{2}}\nu } \right)$ представлены на фиг. 2а. Пространственные переменные $x$, $y$ нормализованы на длину $L$. Логарифмы ошибок схемы (32) относительно логарифма шага пространственной решетки $\lg (h) = \lg (L{\text{/}}K)$, $K = 15,\;25,\;49,\;73$ представлены на фиг. 2б, наклон линии линейной регрессии $a$, где $\lg ({\text{error}}) = a\lg (h) + b$, построенной на логарифмах ошибок, равен $a = 1.88$ , т.е. близок $2$, что подтверждает второй порядок сходимости схемы.

Фиг. 2.

Вихри Тэйлора–Грина. Линии тока для разрешения $25 \times 25$ в момент времени $1{\text{/}}(2{{k}^{2}}\nu )$, см. (а). Пространственные переменные $x$, $y$ нормализованы на длину $L$. Логарифмы ошибок схемы относительно логарифма шага пространственной решетки $\lg (h) = \lg (L{\text{/}}K)$, $K = 15,\;25,\;49,\;73$, см. (б), наклон линии линейной регрессии $a$, где $\lg ({\text{error}}) = a\lg (h) + b$, построенной на данных логарифма ошибок, равен $a = 1.88$.

4.2. Течение в полости с подвижной границей

В данной задаче снова рассматривается квадратная область, но границы области представляют собой твердые непроницаемые стенки. При этом верхняя стенка является подвижной, ее скорость, направленная параллельно стенке, равна ${{U}_{0}} = 0.03$. Число Рейнольдса определяется следующим образом: ${\text{Re}} = {{U}_{0}}H{\text{/}}\nu $, где $L = K \times c$ – длина или ширина области, $K$ есть число узлов в направлении оси $x$ или $y$ (всего в рассматриваемом квадрате $K \times K$ узлов), $\nu $ – вязкость. На границах ставятся стандартные условия отскока, применяемые в теории РУБ [1], для данных условий считается, что твердая граница находится на расстоянии половины длины шага решетки относительно узлов, находящихся рядом с границей. Для таких приграничных узлов имеем

${{f}_{{\overline i }}} = f_{i}^{ + } - \frac{{2{{w}_{i}}{{\rho }_{w}}}}{{{{\theta }_{0}}}}{{{\mathbf{c}}}_{i}}{{{\mathbf{U}}}_{w}},$
где $\overline i $ – индексы, соответствующие скоростям ${{{\mathbf{c}}}_{{\overline i }}}$, которые направлены от стенки, ФР для этих скоростей не могут быть найдены из дискретной кинетической модели. Также индексы $i$ соответствуют скоростям ${{{\mathbf{c}}}_{i}} = - {{{\mathbf{c}}}_{{\overline i }}}$, а ФР $f_{i}^{ + }$ в приграничных узлах, соответствующие скоростям ${{{\mathbf{c}}}_{i}}$, получаются после применения шага столкновений в этих узлах (но без адвекции); кроме того, ${{U}_{w}}$ – скорость вещества на стенке, она равна $0$ или ${{U}_{0}}$; ${{\rho }_{w}}$ – плотность вещества на стенке, в данной задаче считаем ${{\rho }_{w}} = 1$.

Для данной задачи неизвестны аналитические решения, вопрос скорости сходимости здесь не исследуется, но имеются результаты численного моделирования [30], которые часто используют для оценки точности и устойчивости РУБ для разных чисел Рейнольдса и разрешения сеток [28], [31], [32]. Начальное распределение задается в форме абсолютного максвелловского распределения, в ходе решения схемы (27)–(29) профили скорости могут сходиться, не сходиться к какому-то пределу, или вообще разрушаться. Будем считать, что численное решение сходится к некоторому пределу, если существует такой момент времени $t$, что имеем

$\frac{{\sqrt {\sum\limits_{x,y} {{{{\left( {u_{x}^{{(n)}}(t + 1000,x,y) - u_{x}^{{(n)}}(t,x,y)} \right)}}^{2}}} + {{{\left( {u_{y}^{{(n)}}(t + 1000,x,y) - u_{y}^{{(n)}}(t,x,y)} \right)}}^{2}}} }}{{\sqrt {\sum\limits_{x,y} \,u_{x}^{{(n)}}{{{(t,x,y)}}^{2}} + u_{y}^{{(n)}}{{{(t,x,y)}}^{2}}} }}{{ < 10}^{{ - 6}}},$
где суммирование идет по всем узлам $x$, $y$, также $u_{x}^{{(n)}}$, $u_{y}^{{(n)}}$ суть численные решения. Сходимость зависит от числа Рейнольдса, разрешения сетки. В настоящей работе мы рассмотрим случай числа Рейнольдса (${\text{Re}}$), равного $1000$. Особенностью настоящей схемы является то, что для данной задачи она позволяет получать устойчивые решения при грубом разрешении сетки. Для ${\text{Re}} = {\text{10}}00$ рассматриваемая схема позволяет получить сходящиеся устойчивые решения при очень маленьком разрешении, например $K \times K$, равном $23 \times 23$. Результаты моделирования для данного разрешения представлены на фиг. 3 и фиг. 4. На фиг. 3 представлены линии тока, решение устойчиво, качественно получается верная картина, следует отметить, что дополнительные вихри в углах не воспроизводятся, это связано с тем, что разрешение сетки очень грубое. На фиг. 4 представлено сравнение скоростей ${{u}_{x}}(L{\text{/}}2,y)$ и ${{u}_{y}}(x,L{\text{/}}2)$, обозначенные как ${{U}_{x}}$, ${{U}_{y}}$, с численным эталонным решением [30].

Фиг. 3.

Течение в полости с подвижной верхней границей. Показаны линии тока для разрешения $23 \times 23$ и числа Рейнольдса ${\text{Re}} = 1000$.

Фиг. 4.

Течение в полости с подвижной верхней границей. Показаны скорости ${{u}_{x}}(L{\text{/}}2,y)$ и ${{u}_{y}}(x,L{\text{/}}2)$ для разрешения $23 \times 23$, ${\text{Re}} = 1000$. Квадратами отмечено численное решение Ghia, Ghia, Shin [30].

Отдельно возникает вопрос минимального разрешения сетки, при котором сходятся решения стандартной решеточной девятискоростной модели Больцмана $D2Q9$. Для повышения устойчивости модели $D2Q9$ дополнительно проведем стандартную процедуру регуляризации [33]. Для ${\text{Re}} = {\text{10}}00$ регуляризованная модель $D2Q9$ при разрешении $K \times K$, где $K < 100$ неустойчива. Решения становятся устойчивыми примерно при $K > 140$. Этот результат согласуется с предыдущими исследованиями [28]. Можно заключить, что для данной задачи схема (27)–(29) для дискретной кинетической модели (4)–(12) позволяет получать сходящиеся решения для значительно меньших разрешений решетки, чем в случае стандартных РУБ.

При увеличении разрешения численные решения для схемы (27)–(29) приближаются к эталонному решению [30]. Результаты моделирования при $K = 127$ показаны на фиг. 5 и фиг. 6. Очевидно, дополнительные вихри в нижних углах воспроизводятся схемой (фиг. 5), как и должно быть, профили скоростей ${{u}_{x}}(L{\text{/}}2,y)$ и ${{u}_{y}}(x,L{\text{/}}2)$, обозначенные как ${{U}_{x}}$, ${{U}_{y}}$ на фиг. 6, практически совпадают с эталонным решением [30] (фиг. 6).

Фиг. 5.

Течение в полости с подвижной верхней границей. Показаны линии тока для разрешения $127 \times 127$ и числа Рейнольдса ${\text{Re}} = 1000$.

Фиг. 6.

Течение в полости с подвижной верхней границей. Показаны скорости ${{u}_{x}}(L{\text{/}}2,y)$ и ${{u}_{y}}(x,L{\text{/}}2)$ для разрешения $127 \times 127$, ${\text{Re}} = 1000$. Квадратами отмечено численное решение Ghia, Ghia, Shin [30].

5. ЗАКЛЮЧЕНИЕ

В настоящей работе построена аппроксимация для дискретной кинетической модели уравнения Больцмана для девяти скоростей на основе суммы Вайлда, которая в сплошносредном пределе описывает вязкую несжимаемую жидкость. Данная аппроксимация представляет собой разностную схему второго порядка точности по пространственным координатам и времени. Рассмотрены две тестовые задачи: вихри Тэйлора–Грина и течение в полости с подвижной верхней границей. Показано, что схема на основе суммы Вайлда может описывать течение в полости с подвижной верхней границей при очень грубом разрешении сетки и не малых числах Рейнольдса. При этом регуляризованное решеточное уравнение Больцмана $D2Q9$ для тех же условий неустойчиво.

Следует отметить, что для дискретной кинетической модели (4)–(12) выполняется $H$-теорема, но дальнейшая дискретизация этой системы по пространству и времени может нарушать это свойство. Выполнение $H$-теоремы, конечно, очень желательно, так как оно означает не только линейную, но и нелинейную устойчивость, что важно для приложений. Для четырехскоростной модели Броудэлла схема, описывающая столкновения в виде сумм Вайлда, обладает $H$-теоремой [27]. Таким образом, исследование схем на основе сумм Вайлда перспективно, так как здесь имеется большая вариативность – можно построить очень большое количество разностных схем с желаемыми свойствами. Для схемы в виде суммы Вайлда (27)–(29) для модели (4)–(12) доказать $H$-теорему не получится, хотя моделирование течения в полости с подвижной верхней границей показывает лучшую устойчивость, чем моделирование с использованием регуляризованной модели $D2Q9$.

Это не означает, что сумма Вайлда (27)–(29) устойчивеe РУБ $D2Q9$ в любых задачах. Действительно, линейная устойчивость РУБ исследовалась многими авторами [34]–[43]. Важно, что устойчивость РУБ существенно зависит от исследуемой задачи. Например, для возмущений в виде бегущих волн неустойчивость нерегуляризованных РУБ может наблюдаться в узком интервале волновых векторов, причем возмущение должно быть направлено под определенными углами относительно скорости потока. Кроме того, ряд регуляризованных РУБ, которые в модельных задачах показывают хорошую устойчивость, неожиданно оказываются неустойчивы для длинноволновых возмущений [42]. Предварительный анализ показывает, что у схемы (27)–(29) для уравнений (4)–(12) также существуют возмущения, ведущие к неустойчивости, что, впрочем, не означает, что неустойчивости будут наблюдаться для основных модельных течений. Можно заключить, что анализ линейной устойчивости и поиск схемы с выполняющейся $H$-теоремой представляют интересную задачу для будущих исследований.

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

  1. Krüger T., Kusumaatmaja H., Kuzmin A., Shardt O., Silva G., Viggen E. The Lattice Boltzmann Method. Principles and Practice. Swizerland: Springer, 2017.

  2. Kogan M. Rarefied Gas Dynamics. New York: Plenum Press, 1969.

  3. Broadwell J. Shock structure in a simple discrete velocity gas // Phys. Fluids. 1964. V. 7. P. 1243–1247.

  4. Godunov S., Sultangazin U. On discrete models of the kinetic Boltzmann equation // Russian Math. Surveys. 1971. V. 26. P. 1–56.

  5. Gatignol R. The hydrodynamical description for a discrete velocity model of gas // Complex Systems. 1987. V. 1. P. 709–725.

  6. Wagner A. An H-theorem for the lattice Boltzmann approach to hydrodynamics // Europhys. Lett. 1998. V. 44. P. 144–149.

  7. Yong W.-A., Luo L.-S. Nonexistence of H theorems for the athermal lattice Boltzmann models with polynomial equilibria // Phys. Rev. E. 2003. 051105.

  8. Yong W.-A., Luo L.-S. Nonexistence of H Theorem for some Lattice Boltzmann models // J. Stat. Phys. 2005. V. 121. P. 91–103.

  9. Karlin I., Ferrante A., Öttinger H. Perfect entropy functions of the Lattice Boltzmann method // Europhys. Lett. 1999. V. 47. P. 182–188.

  10. Ansumali S., Karlin I., Öttinger H. Minimal entropic kinetic models for hydrodynamics // Europhys. Lett. 2003. V. 63. P. 798–804.

  11. Karlin I., Ansumali S., Frouzakis C., Chikatamarla S. Elements of the Lattice Boltzmann Method I: Linear Advection Equation // Commun. Comput. Phys. 2006. V. 1. P. 616–655.

  12. Karlin I., Chikatamarla S., Ansumali S. Elements of the lattice Boltzmann method II: kinetics and hydrodynamics in one dimension // Commun. Comput. Phys. 2007. V. 2. P. 196–238.

  13. Chen H., Goldhirsch I., Orszag S. Discrete rotational symmetry, moment isotropy, and higher order lattice Boltzmann models // J. Sci. Comput. 2008. V. 34. P. 87–112.

  14. Bobylev A., Spiga G. On a class of exact two-dimensional stationary solutions for the Broadwell model of the Boltzmann equation // J. Phys. A: Math. Gen. 1994. V. 27. P. 7451–7459.

  15. Bobylev A. Exact solutions of discrete kinetic models and stationary problems for the plane Broadwell model // Math. Methods Appl. Sci. 1996. V. 19. P. 825–845.

  16. Bobylev A., Toscani G. Two dimensional half-space problems for the Broadwell discrete velocity model // Contin. Mech. Termodyn. 1996. V. 8. P. 257–274.

  17. Bobylev A., Caraffini G., Spiga G. Non-stationary two-dimensional potential flows by the Broadwell model equations // Eur. J. Mech. B, Fluids. 2000. V. 19. P. 303–315.

  18. Ilyin O. The analytical solutions of 2D stationary Broadwell kinetic model // J. Stat. Phys. 2012. V. 146. P. 67–72.

  19. Ильин О. Симметрии, функция тока и точные решения для двумерной стационарной кинетической модели Бродуэлла // Теор. Мат. Физ. 2014. Т. 179. С. 350–359.

  20. Uchiyama K. On the Boltzmann-Grad limit for the Broadwell model of the Boltzmann equation // J. Stat. Phys. 1988. V. 52. P. 331–355.

  21. Babovsky H. “Small” kinetic models for transitional flow simulations // AIP Conf. Proc. 2012. V. 1501. P. 272–278.

  22. Babovsky H. Discrete kinetic models in the fluid dynamic limit // Comput. and Math. with Appl. 2014. V. 67. P. 256–271.

  23. Ilyin O. Discrete Velocity Boltzmann Model for Quasi- Incompressible Hydrodynamics // Math. 2021. V. 9. P. 993.

  24. Wild E. On Boltzmann’s Equation in the Kinetic Theory of Gases // Proc. Camb. Philos. Soc. 1951. V. 47. P. 602–609.

  25. McKean H. An exponential formula for solving Boltzmann’s equation for a Maxwellian gas // J. of Combinat. Theory. 1967. V. 2. P. 358–382.

  26. Gabetta E., Pareschi L., Toscani G. Wild’s sums and numerical approximation of nonlinear kinetic equations // Transp. Theory and Stat. Phys. 1996. V. 25. P. 515–530.

  27. Gabetta E., Pareschi L., Toscani G. Relaxation Schemes for Nonlinear Kinetic Equations // SIAM J. Numer. Anal. 1997. V. 34. P. 2168–2194.

  28. Latt J. Hydrodynamic Limit of Lattice Boltzmann Equations. Dissertation. Geneva: University of Geneva, 2007.

  29. Chauvat P., Gatignol R. Euler and Navier-Stokes description for a class of discrete models of gases with different moduli // Transp. Theory Stat. Phys. 1992. V. 21. P. 417–435.

  30. Ghia U., Ghia K., Shin C. High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method // J. Comp. Phys. 1982. V. 48. P. 387–411.

  31. Montessori A., Falcucci G., Prestininzi P., La Rocca M., Succi S. Regularized lattice Bhatnagar–Gross–Krook model for two- and three-dimensional cavity flow simulations // Phys. Rev. E. 2014. V. 89. P. 053317.

  32. Wu J.-S., Shao Y.-L. Simulation of lid-driven cavity flows by parallel lattice Boltzmann method using multi-relaxation-time scheme // Int. J. Numer. Meth. Fluids. 2004. V. 46. P. 921–937.

  33. Latt J., Chopard B. Lattice Boltzmann method with regularized pre-collision distribution functions // Math. Comput. in Simul. 2006. V. 72. P. 165–168.

  34. Sterling J., Chen S. Stability Analysis of Lattice Boltzmann Methods // J. Comp. Phys. 1996. V. 123. P. 196–206.

  35. Lallemand P., Luo L.-S. Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability // Phys. Rev. E. V. 61. P. 6546–6562.

  36. Siebert D., Hegele Jr. L., Philippi P. Lattice Boltzmann equation linear stability analysis: Thermal and athermal models // Phys. Rev. E. 2008. V. 77. P. 026707.

  37. Servan-Camas B., Tsai F. Non-negativity and stability analyses of lattice Boltzmann method for advection–diffusion equation // J. Comp. Phys. 2009. V. 228. P. 236–256.

  38. Marié S., Ricot D., Sagaut P. Comparison between lattice Boltzmann method and Navier–Stokes high order schemes for computational aeroacoustics // J. Comp. Phys. 2009. V. 228. P. 1056–1070.

  39. Hosseini S.A., Darabiha N., Thèvenin D., Eshghinejadfard A. Stability limits of the single relaxation-time advection–diffusion lattice Boltzmann scheme // Int. J. Mod. Phys. C. 2017. V. 28. 1750141.

  40. Wissocq G., Sagaut P., Boussuge J.-F. An extended spectral analysis of the lattice Boltzmann method: modal interactions and stability issues // J. Comp. Phys. 2019. V. 380. P. 311–333.

  41. Masset P.-A., Wissocq G. Linear hydrodynamics and stability of the discrete velocity Boltzmann equations // J. Fluid Mech. A. 2020. V. 897. P. 29.

  42. Wissocq G., Coreixas C., Boussuge J.-F. Linear stability and isotropy properties of athermal regularized lattice Boltzmann methods // Phys. Rev. E. 2020. V. 102. 053305.

  43. Wissocq G., Sagaut P. Hydrodynamic limits and numerical errors of isothermal lattice Boltzmann schemes // arXiv:2104.14217v1 2021.

Дополнительные материалы отсутствуют.

Инструменты

предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал