
\documentclass[10pt]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}

\def\udcs{519.63, 533.5, 533.2} %Здесь автор определяет УДК своей работы
\def\mscs{76P05, 65N08} %Здесь автор определяет классификаторы AMS своей работы
\setcounter{page}{144}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}

\def\logo{{\bf\huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}

\def\semrtop
     {
  \vbox{
     \noindent\logo
     \hspace{80mm}\raisebox{1ex}{ISSN 1813-3304 }

     \vspace{5mm}

     \begin{center}
     {\huge СИБИРСКИЕ \ ЭЛЕКТРОННЫЕ} \\[2mm]
     {\huge МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ} \\[2mm]
     {\large Siberian Electronic Mathematical Reports} \\[1mm]
     {\LARGE\tt{http://semr.math.nsc.ru}}\\[0.5mm]
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
     \vspace{-3mm}
     \noindent
     \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa} \\
     \hline\hline
     \end{tabular}

     \vspace{1mm}
     {\flushleft\it Том XX, стр. XXX--XXX (20XX) \hspace{65mm}{\rm\small УДК \udcs}}
     \newline
        {\rm\small DOI~10.33048/semi.2019.16.xxx}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}{\rm\small MSC\ \ \mscs }
  }%\vbox
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% окружение для вставки рисунка
%  1 аргумент - имя файла без расширения,
%               это же имя является меткой рисунка
%  2 аргумент - ширина рисунка - доли страницы
%  3 аргумент - текст подрисуночной подписи
%
% максимальная ширина рисунка
\newlength{\figwidth}
\figwidth=\textwidth
%
\newenvironment{i-fig}[3]{
\begin{figure}[!ht]
\centering\includegraphics[width={#2\figwidth},keepaspectratio=true]{#1.eps}\vspace{-1.5em}\\
\noindent
\parbox[t]{\figwidth}{\renewcommand{\baselinestretch}{0.65}
\caption{#3}\label{#1}}
\end{figure}\par}\par{}

% вектор
\def \mb#1{{\boldsymbol{#1}}}

\begin{document}
\renewcommand{\refname}{References}


\thispagestyle{empty}

\title[Дискретизация уравнения Больцмана]{Дискретизация уравнения Больцмана при помощи метода конечных объемов и явно-неявных схем}
\author{{К.Н. ВОЛКОВ, В.Н. ЕМЕЛЬЯНОВ, А.В. ПУСТОВАЛОВ}}%
\address{Konstantin Volkov, Vladislav Emelyanov, Alexander Pustovalov % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Baltic State Technical University,
\newline\hphantom{iii} 1-ya Krasnoarmeyskaya ul., 1,
\newline\hphantom{iii} 190005, St Petersburg, Russia}%
\email{dsci@mail.ru}%

\thanks{\sc Volkov, K.N., Emelyanov V.N., Pustovalov, A.V,
Discretization of Boltzmann equation with finite volume method and explicit-implicit schemes}
\thanks{\copyright \ 2022 Волков К.Н., Емельянов В.Н., Пустовалов А.В.}
%\thanks{\rm Работа поддержана РФФИ (грант 03-01-11111)}
\thanks{\it Поступила 10 декабря 2022 г., опубликована XX декабря 20XX г.}%


\semrtop \vspace{1cm}
\maketitle {\small
\begin{quote}
\noindent{\sc Abstract. } The features of discretization of the Boltzmann equation using the finite volume method are
considered. Finite-difference schemes for calculation of fluxes and finite-difference schemes for discretization in time are
discussed. A TVD-type scheme is used for flux discretization, and an explicit-implicit scheme is applied to time
discretization. The results of numerical simulation of rarefied gas flow in a shock tube for various Knudsen numbers are
presented. For small Knudsen numbers, the solution of the Boltzmann equation is compared with the solution obtained from the
Euler equation.\medskip

\noindent{\bf Keywords:} finite volume method, Boltzmann equation, rarefied gas, shock tube.
 \end{quote}
}

\section{Введение}

Решение интегро-дифференциального уравнения Больцмана является основой моделирования течений разреженного газа. Функция
распределения, для которой записывается уравнение Больцмана, зависит от семи переменных~--- времени, пространственных
координат, компонент скорости и времени, а интеграл столкновений является пятикратным \cite{Bhatnagar54}. Трудности
численного решения кинетического уравнения состоят в его высокой размерности, присутствии многомерного интеграла
столкновений и сложности построения неявных методов интегрирования по времени.

Число операций, необходимых для расчета интеграла столкновений (\ref{ha01a}) прямым методом, имеет оценку $O(N^7)$
\cite{Bobylev95}, где $N$~--- число узлов сетки (дискретных скоростей) в скоростном пространстве в каждом координатном
направлении (при этом подход имеет первый порядок точности по скорости). При использовании модели дискретных скоростей и
модели твердых сфер число операций оценивается как $O(\overline N^3 N^3\log N)$, где $\overline N\ll N$ \cite{Mouhot13}. Для
расчета интеграла столкновений широкое применение находит метод дискретного преобразования Фурье и спектральные методы
\cite{Hu20, Wu13}. Вычислительные затраты, связанные с использованием дискретного преобразования Фурье, оцениваются как
$O(N^6\log N)$ \cite{Boelens20}. Применение спектрального метода для столкновительного ядра специального вида приводит к
следующей оценке числа операций $O(N^6)$ \cite{Pareschi96}. Такая же оценка имеет место при использовании разрывного метода
Галеркина \cite{Alekseenko19}. Подход к расчету столкновительного оператора для максвелловских молекул, развитый в работе
\cite{Bobylev97}, приводит к оценке числа операций $O(N^4)$. Однако подход имеет порядок точности $O(N^{-1/2})$.
Использование метода быстрого преобразования Фурье (Fast Fourier Transformation, FFT) в продольном направлении,
преобразования Хенлеля (Hankel transformation) в поперечном направлении и свойства цилиндрической симметрии позволяет
уменьшить число операций до $O(N^2\log N)$ \cite{Watch09}. Применение представления Карлемана (Carleman representation) и
преобразований Радона (Radon transformation) позволяет повысить порядок точности до $O(N^{-2})$ при числе операций порядка
$O(N^6\log N)$ \cite{Bobylev99}. В подходе \cite{Ibragimov02} число операций, необходимых для расчета интеграла столкновений
для модели твердых сфер, составляет $O(N^6)$, а точность метода имеет порядок $O(N^{-2})$.

Для ускорения расчетов находят применение технологии машинного обучения \cite{Holloway15}. Их применение позволяет уменьшить
число операций до $O(N^3)$ (оценка получена для пространственно-однородного уравнения Больцмана).

Наиболее развитыми численными методами решения уравнения Больцмана являются методы прямого статистического моделирования
Монте-Карло (Direct Simulation Monte Carlo, DSMC). Решение нестационарных задач динамики разреженного газа является
достаточно затратным с вычислительной точки зрения. В связи с этим, получили распространение методы прямого решения
уравнения Больцмана, основанные на конечно-объемном подходе (Finite Volume Method, FVM) к его дискретизации \cite{Dimarco13,
Dimarco18}.

Жесткость оператора столкновений вида при больш\'их частотах столкновений накладывает определенные требования на выбор
разностных схем для дискретизации по времени. Использование явных разностных схем приводит к достаточно малым шагам по
времени. Для дискретизации по времени находят применение неявная схема Эйлера (backward Euler scheme) \cite{Coron91} и
явно-неявные (implicit--explicit, IMEX) схемы Рунге-Кутты (Runge--Kutta scheme) \cite{Pieraccini07}, в которых для
дискретизации конвективных слагаемых используется явная схема, а для расчета столкновительного оператора~--- неявная схема
(IMEX-RK scheme). Для получения консервативной устойчивой схемы высокого порядка точности применяются дискретизации TVD и
WENO.

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


\section{Уравнение Больцмана и его свойства}

Уравнение Больцмана, описывающее эволюцию функции распределения молекул $f(t,\mb x,\mb v)$, имеющих скорость $\mb v$ в точке
$\mb x$ в момент времени $t$, имеет вид \cite{Bhatnagar54}
\begin{align}
\frac{\partial f}{\partial t}+\mb v\cdot\nabla f=Q(f,f), \label{ha01}
\end{align}
где $t$~--- время, $\mb x=(x,y,z)$~--- положение частицы в физическом пространстве, $\mb v=(v_x,v_y,v_z)$~--- вектор
скорости частицы, $Q(f,f)$~--- интеграл столкновений. В случае упругих столкновений с конечным радиусом взаимодействия
столкновительный оператор в уравнении Больцмана (\ref{ha01}) расщепляется на две части, соответствующие прямым и обратным
столкновениям
\begin{align}
Q(f,f)=-\nu(f)f+N(f,f), \label{ha01a}
\end{align}
где $\nu(f)$~--- частота столкновений. При этом слагаемое $\nu(f)f$ представляет собой интеграл прямых столкновений, а $N(f,
f)$~--- интеграл обратных столкновений. Для решения уравнения Больцмана (\ref{ha01}) с интегралом столкновений (\ref{ha01a})
применяется схема расщепления по физическим процессам. В течение малого шага по времени сначала происходит
бесстолкновительный перенос частиц, за которым следует пространственно-однородный этап столкновений между частицами.
Гидродинамическая (макроскопическая) консервативность, связанная с выполнением дискретных аналогов соотношений сохранения
для моментных уравнений, обеспечивается процедурой коррекции для функции распределения после этапа релаксации.

Решение уравнения Больцмана осложняется не только высокой размерностью фазового пространства, но и нелинейностью и
нелокальным характером столкновительного оператора. Для упрощения рассматривается уравнение BGK (Bhatnagar--Gross--Krook) со
столкновительным оператором вида \cite{Bhatnagar54}
\begin{align}
Q_\text{BGK}(f,f)=\lambda\left[M_f(t,\mb x,\mb v)-f(t,\mb x,\mb v)\right], \label{ha02}
\end{align}
где $\lambda=1/\tau$~--- частота столкновений между молекулами, $\tau$~--- характерное время. Существование неотрицательного
решения уравнения BGK в неограниченном пространстве и ограниченной области доказывается в работе \cite{Perthame89}, а в
области с периодическими граничными условиями~--- в работе \cite{Mischler96}.

Равновесная (максвелловская) функция распределения в пространстве размерности $d$ представляется в виде
\begin{align}
M_f(t,\mb x,\mb v)=n\left(\frac{m}{2\pi\theta}\right)^{d/2}\exp\left(-\frac{m\left|\mb v-\mb u\right|^2}{2\theta}\right),
\label{ha03}
\end{align}
где $m$~--- масса частицы, $n$~--- концентрация частиц, $\theta$~--- температура. В частности, функция распределения по
скоростям для частиц, входящих в расчетную область, задается в виде максвелловской функции распределения (\ref{ha03}).

Вектор консервативных переменных $\mb W=\{\rho,\rho\mb v,\rho e\}'$ определяется при помощи моментов функции распределения
\begin{align*}
\mb W=\int\mb\psi(\mb v)f\,d\mb v,
\end{align*}
где $\mb\psi=\{1,\mb v, \left|\mb v\right|^2\}'$~--- вектор столкновительных инвариантов. Функция $\psi(\mb v)$ является
столкновительным инвариантом (collision invariant), если выполняется условие
\begin{align*}
\int\limits_{R^d}\psi(\mb v)\left(M_f-f\right)d\mb v=0.
\end{align*}

Макроскопические свойства среды определяются первыми тремя моментами функции распределения. Вычисляя моменты различного
порядка функции распределения, получим макроскопические распределения концентрации частиц, скорости, энергии и температуры
\begin{align*}
&n(t,\mb x)=\int\limits_{R^d}f\,d\mb v;\\
&\mb u(t,\mb x)=\frac 1n\int\limits_{R^d}\mb v f \,d\mb v;\\
&E(t,\mb x)=\int\limits_{R^d}\frac 12\left|\mb v\right|^2 f\,d\mb v;\\
&\frac d2 n(t,\mb x)\theta(t,\mb x)=\int\limits_{R^d}\frac m2\left|\mb v-\mb u\right|^2 f\,d\mb v.
\end{align*}
Тензор вязких напряжений и вектор потока тепла определяются при помощи соотношений
\begin{align*}
\mb\tau=-\int\mb c\mb c\left(f-M\right)d\mb v,\quad \mb q=\frac 12\int\mb c c^2 f\,d\mb v,
\end{align*}
где $\mb c=\mb v-\mb u$. Уравнения, описывающие течение сплошной среды (в пределе при $\lambda\to\infty$ или $\tau\to 0$),
получаются путем умножения уравнения BGK на столкновительные инварианты и последующего интегрирования полученных уравнений.

Оператор столкновений в виде (\ref{ha02}) удовлетворяет консервативным и энтропийным свойствам исходного оператора
столкновений. Определим H-функционал (H-functional)
\begin{align*}
H[f](t,\mb x)=\int\limits_{R^d}f(t,\mb x,\mb v)\ln\left[f(t,\mb x,\mb v)\right]\,d\mb v.
\end{align*}
Функционал $H[f](t,\mb x)$ представляет собой энтропию (в термодинамике энтропия определяется с противоположным знаком).
Функция $\ln(M_f)$ является столкновительным инвариантом. При этом имеет место свойство диссипации
\begin{align}
\frac{d}{dt}\int\limits_{R^d}H[f]d\mb x\leqslant 0.\label{Hcond}
\end{align}
При реализации численных методов необходимо проверять выполнение условия диссипации (\ref{Hcond}) на дискретном уровне. В
общем случае проверка энтропийного условия является достаточно затратным с вычислительной точки зрения.

Произведя замену переменных $\mb x'=\varepsilon\mb x$ и $t'=\varepsilon t$, нетрудно показать, что при $\varepsilon\to 0$
решение уравнения BGK стремится к решению уравнений Эйлера. При этом сходимость уравнения BGK к решению уравнений
Навье--Стокса не имеет места \cite{Nassios13}. Оператор столкновений в виде BGK не позволяет получить корректное значение
числа Прандтля, поскольку частота столкновений не зависит от скорости. Имеются другие приближения уравнения Больцмана,
свободные от этого недостатка. К ним, в частности, относятся уравнение Шахова (Shakhov equation) \cite{Shakhov68}, уравнение
ES-BGK \cite{Holway66}, а также модель, предложенная в работе \cite{Mieussens04} (Mieussens--Struchtrup model).

Для получения корректного значения числа Прандтля в гидродинамическом пределе частота столкновений $\lambda(\mb x,t)>0$
принимается пропорциональной концентрации частиц и температуре \cite{Mieussens00b}. Время релаксации находится из
соотношения $\tau=\mu/p$. Для вязкости принимается степенная зависимость от температуры $\mu=\mu_0(T/T_0)^\omega$, где
$\mu_0$~--- динамическая вязкость при характерной температуре $T_0$. Показатель степени $\omega$ зависит от модели
взаимодействия молекул и типа газа (для модели твердых сфер $\omega=0.5$).

В работах \cite{Mieussens00b, Mieussens00a} построены разностные схемы для уравнения BGK, позволяющие обеспечить выполнение
условий консервативности массы, количества движения и энергии на дискретном уровне. Реализация такого подхода требует
решения системы из 5 нелинейных уравнений для модели BGK и 10 уравнений для модели ES-BGK в каждом узле сетки. Полностью
неявные схемы для диcкретизации кинетического уравнения общего вида, описывающего поведение разреженного газа и
бесстолкновительной плазмы (уравнение Больцмана--Власова), предлагаются в работе \cite{Garrett18}.


\section{Конечно-объемная дискретизация}

Полагая $d=1$, рассмотрим двумерное фазовое пространство ($1x\times 1v$). Для дискретизации уравнения BGK используется метод
конечных объемов в прямоугольной области $\Omega\times V$.

\subsection{Расчет потоков}

Дискретизация проводится в усеченном пространстве скоростей $V=[-v_\text{max},+v_\text{max}]$, где $v_\text{max}>0$, на
пространственном интервале $\Omega=[-L,+L]$, где $L>0$. В общем случае допускается выбор несимметричного интервала изменения
скорости. Строятся равномерные координатные и скоростные сетки с шагами
\begin{align*}
h_x=\frac{2L}{N_x},\quad h_v=\frac{2v_\text{max}}{N_v},
\end{align*}
где $N_x$ и $N_v$~--- число узлов в координатной и скоростной области. Координаты центров ячеек в пространстве координат и
скоростей находятся из соотношений
\begin{align*}
x_l=-L+\left(l-\frac12\right)h_x,\quad v_l=-v_\text{max}+\left(l-\frac12\right)h_v,
\end{align*}
где $l=1/2,1,\ldots$\,.

Расчетная область $\Omega\times V\subset R^2$ разбивается на прямоугольные ячейки в пространстве скоростей и координат
\begin{align*}
&\Omega\times V=\bigcup_{i,j=1}^{N_x,N_v}C_{i,j}=\bigcup_{i,j=1}^{N_x,N_v}C^x_{i}\times
C^v_{j}=\\
&=\bigcup_{i,j=1}^{N_x,N_v}\left[x_{i-1/2},x_{i+1/2}\right]\times\left[v_{j-1/2},v_{j+1/2}\right].
\end{align*}
Точка $(x_i,v_j)$ располагается в центре ячейки $C_{i,j}=C^x_{i}\times C^v_{j}$.

Среднее по контрольному объему (ячейке) значение функции распределения определяется следующим образом (теорема о среднем)
\begin{align*}
\overline f_{i,j}(t)=\overline f(t,x_i,v_j)=\frac{1}{|C^x_{i}\times C^v_{j}|}\int\limits_{C^x_{i}\times
C^v_{j}}f(t,x,v)\,dx\,dv.
\end{align*}

Среднее по контрольному объему значение равновесной функции распределения находится из соотношения
\begin{align*}
\overline M_{i,j}(t)=\frac{1}{|C_{i,j}|}\int\limits_{C_{i,j}}M_f(t,x,v)\,dx\,dv.
\end{align*}

Уравнение Больцмана, записанное с учетом аппроксимации BGK (\ref{ha01}) с интегралом столкновений (\ref{ha02}), в
интегральной форме имеет вид
\begin{align*}
\partial_t\overline f_{i,j}(t)+\frac{1}{|C_{i,j}|}\int\limits_{C_{i,j}}\partial_x(vf)dx\,dv=\frac{\lambda}{|C_{i,j}|}\int\limits_{C_{i,j}}M_f(t,x,v)-\lambda\overline f_{i,j}(t).
\end{align*}
Введем поток функции распределения
\begin{align*}
F(t,x,v)=vf(t,x,v).
\end{align*}
Используя теорему о дивергенции, получим
\begin{align*}
&\partial_t\overline f_{i,j}(t)+\lambda\overline f_{i,j}(t)-\lambda\overline M_{i,j}(t)=\\[0.2cm]
&=-\frac{1}{h_x h_v}\int\limits_{C_j^v}\left[F(t,x_{i+1/2},v)-F(t,x_{i-1/2},v)\right]dv.
\end{align*}
Для дискретизации потокового интеграла используется соотношение
\begin{align*}
&\frac{1}{h_x
h_v}\int\limits_{C_j^v}\left[F(t,x_{i+1/2},v)-F(t,x_{i-1/2},v)\right]dv=\\
&=\frac{F(t,x_{i+1/2},v)-F(t,x_{i-1/2},v)}{h_x}.
\end{align*}

Кусочно-линейная реконструкции функции в ячейке $C_{i,j}$ имеет вид
\begin{align*}
p_{i,j}(t,x,v)=\overline f_{i,j}(t)+\sigma_{i,j}^x(x-x_i)+\sigma_{i,j}^v(v-v_j).
\end{align*}
Для ограничения потока используется ограничитель minmod
\begin{align*}
&\sigma_{i,j}^x=\text{minmod}\left\{\left(\frac{\overline f_{i+1,j}-\overline f_{i-1,j}}{2h_x}\right),\theta^x\left(\frac{\overline f_{i,j}-\overline f_{i-1,j}}{h_x}\right),\theta^x\left(\frac{\overline f_{i+1,j}-\overline f_{i,j}}{h_x}\right)\right\};\\
%
&\sigma_{i,j}^v=\text{minmod}\left\{\left(\frac{\overline f_{i,j+1}-\overline
f_{i,j-1}}{2h_v}\right),\theta^v\left(\frac{\overline f_{i,j}-\overline f_{i,j-1}}{h_v}\right),\theta^v\left(\frac{\overline
f_{i,j+1}-\overline f_{i,j}}{h_v}\right)\right\}.
\end{align*}
Здесь $\theta^x,\theta^v\in[1,2]$. Функция $\text{minmod}$ имеет вид
\begin{align*}
\text{minmod}\{x,y,z\}=\left\{\begin{array}{ll}
\min(x,y,z), &\text{если}\,\,x,y,z>0,\\
\max(x,y,z), &\text{если}\,\,x,y,z<0,\\
0,                 &\text{иначе}.
\end{array}\right.
\end{align*}
В расчетах полагается, что $\theta^x=2$ и $\theta^v=2$.

Функция распределения слева (индекс $-$) и справа (индекс $+$) от грани контрольного объема находится из соотношений
\begin{align*}
&f^-_{i+1/2,j}(t)=p_{i,j}(t,x_{i+1/2},v_j),\quad f^+_{i+1/2,j}(t)=p_{i+1,j}(t,x_{i+1/2},v_j);\\[0.2cm]
&f^-_{i,j+1/2}(t)=p_{i,j}(t,x_i,v_{j+1/2}),\quad f^+_{i,j+1/2}(t)=p_{i,j+1}(t,x_i,v_{j+1/2}).
\end{align*}
Численный поток через грань контрольного объема находится из соотношения
\begin{align*}
F(t,x_{i+1/2},v_j)=\widetilde F\left[f^+_{i+1/2,j}(t),f^-_{i+1/2,j}(t)\right]=\widehat F_{i+1/2,j}(t).
\end{align*}
Для нахождения численных потоков используются противопоточные разности (upwind difference)
\begin{align*}
\widehat F_{i+1/2,j}(t)=\widetilde F\left[f^+_{i+1/2,j}(t),f^-_{i+1/2,j}(t)\right]=\left\{\begin{array}{ll}
v_j f^-_{i+1/2,j}(t), &\text{если}\,\,v_j\geqslant 0,\\[0.2cm]
v_j f^+_{i+1/2,j}(t), &\text{если}\,\,v_j<0.\\
\end{array}\right.
\end{align*}

Разностная схема для уравнения BGK принимает вид
\begin{align*}
\partial_t\overline f_{i,j}(t)+\lambda\overline f_{i,j}(t)-\lambda\overline M_{i,j}(t)=-\frac{\widehat F_{i+1/2,j}(t)- \widehat F_{i-1/2,j}(t)}{h_x}.
\end{align*}
Равновесная функция распределения в узлах сетки находится из соотношения
\begin{align*}
\overline M_{i,j}(t)=\widehat
M_{i,j}(t)=\frac{n_i(t)}{[2\pi\theta_i(t)]^{1/2}}\exp\left\{-\frac{[v_j-u_i(t)]^2}{2\theta_i(t)}\right\}.
\end{align*}
Макроскопические свойства среды находятся из соотношений
\begin{align*}
&n_i(t)=h_v\sum\limits_{j=1}^{N^v}\overline f_{i,j}(t);\\
&u_i(t)=\frac{h_v}{n_i(t)}\sum\limits_{j=1}^{N^v}v_j\overline f_{i,j}(t);\\
&\theta_i(t)=\frac{h_v}{n_i(t)}\sum\limits_{j=1}^{N^v}|v_j|^2\overline f_{i,j}(t)-|u_i(t)|^2.
\end{align*}
Замена интегралов квадратурными формулами приводит к потере столкновительных инвариантов на дискретном уровне.

\subsection{Дискретизация по времени}

Традиционные явно-неявные разностные схемы высокого порядка точности имеют ограничения на максимально возможных шаг по
времени (условие CFL), что делает разностные схемы такого типа достаточно дорогими с вычислительной точки зрения. На
практике шаг интегрирования по времени ддля схем высокого порядка точности оказывается сравнимым с шагом по времени,
используемым в схеме Эйлера \cite{Gottlieb01, Michel22}.

Безусловно устойчивые неявные схемы TVD типа ($L_\infty$-stable RK scheme) имеют первый порядок точности \cite{Gottlieb01}.
Для построения устойчивых IMEX-RK схем, удовлетворяющих условию TVD, используется подход, предложенный в работе
\cite{Dimarco18}, в котором схема Эйлера первого порядка комбинируется с потенциально осциллирующими схемами IMEX-RK более
высокого порядка точности. Схемы такого типа находят применение для дискретизации уравнений гиперболического типа, в которых
жесткость определяется потоковыми слагаемыми, а не источниковым членом \cite{Michel22}.

Для дискретизации по времени используется явно-неявная (IMEX) схема Рунге--Кутты (RK), удовлетворяющая условию TVD
\cite{Pieraccini07, Michel22}. Уравнение BGK записывается в виде
\begin{align*}
\frac{\partial}{\partial t}f(t,x,v)=T[f](t,x,v)+Q[f](t,x,v),
\end{align*}
Оператор переноса $T$ и столкновительный оператор $Q$ имеют вид
\begin{align*}
&T[f](t,x,v)=-v\partial_x f(t,x,v);\\[0.2cm]
&Q[f](t,x,v)=\lambda\left[M_f(t,x,v)-f(t,x,v)\right].
\end{align*}

В общем виде $s$-этапная схема Рунге-Кутта описывается набором следующих соотношений
\begin{align*}
&y^{(r)}=y^k+\Delta t\sum\limits_{l=1}^{r-1}\widetilde a_{rl}T[y^{(l)}]+\Delta t\sum\limits_{l=1}^{r}a_{rl}Q[y^{(l)}];\\
%
&y^{k+1}=y^k+\Delta t\sum\limits_{r=1}^{s}\widetilde b_r T[y^{(r)}]+\Delta t\sum\limits_{r=1}^{s}b_r Q[y^{(r)}].
\end{align*}
Здесь $r\in\{1,\ldots,s\}$. Неявные $s$-этапные методы требуют вычисления вычислений искомой функции на шаге не менее, чем
$2s$. Привлекательность неявных методов состоит в том, что при всех $s$ существуют методы, которые имеют порядок $p=2s$.
Такой порядок достигается путем специального выбора коэффициентов.

Для явной схемы (Explicit Runge--Kutta, ERK) лишь находящиеся ниже главной диагонали элементы матрицы Бутчера отличны от
нуля. В них переход на новый слой происходит по явным формулам. Явные схемы просты в реализации, а трудоемкость вычислений
мала. В диагонально-неявных схемах (Diagonal Implicit Runge--Kutta, DIRK) на каждой стадии для нахождения очередного
необходимо решить систему алгебраических уравнений при помощи итерационных методов. В полностью неявных схемах (Fully
Implicit Runge--Kutta, FIRK) вычисление всех стадий происходит не последовательно, а одновременно. Для нахождения решения
требуется решить систему алгебраических уравнений. Схемы ERK достаточно хорошо изучены. При этом нельзя построить схему
порядка точности $p>s$, но можно построить схему с $p=s$. При построении неявных схем не удается обеспечить одновременно
высокий порядок аппроксимации и хорошие свойства устойчивости.

Матрица коэффициентов $\{a_{rl}\}$ является матрицей Бутчера. Векторы коэффициентов $b_r$ и $c_r$ рассматриваются как
столбцы, дополняющие квадратную матрицу Бутчера до прямоугольной. Коэффициенты $\widetilde a_{rl}$, $\widetilde b_r$,
$a_{rl}$, $b_r$ приводятся в соответствующих таблицах (Butcher tableaux), имеющих вид (явная схема~--- слева, неявная
схема~--- справа)
\begin{center}
\begin{tabular}{r|l}
$\widetilde{\mb c}$ & $\widetilde A$\\
\hline
                    & $\widetilde{\mb b}^T$\vphantom{\Huge M}\\
\end{tabular}
%
\qquad
\begin{tabular}{r|l}
$\mb c$ & $A$\\
\hline
        & $\mb b^T$\vphantom{\Huge M}\\
\end{tabular}
\end{center}
Матрица $\widetilde A$ (явная часть) является строго нижне-треугольной матрицей (нули на главной диагонали) размерности
$s\times s$, а матрица $A$~--- ниже-треугольной матрицей размерности $s\times s$. Для коэффициентов $\widetilde c$ и $c$
выполняются условия
\begin{align*}
\widetilde c_i=\sum\limits_{j=1}^{i-1}\widetilde a_{ij},\quad c_i=\sum\limits_{j=1}^{i}a_{ij}.
\end{align*}
Условия, которым удовлетворяют другие коэффициенты схем различного порядка точности приводятся в работе \cite{Michel22}.

В развернутом виде имеем (явная схема~--- слева, неявная схема~--- справа)
\begin{center}
\begin{tabular}{c|cccc}
0              & 0                 & 0          & \ldots & 0\\
$\widetilde c_2$ & $\widetilde a_{21}$ & 0          & \ldots & 0\\
\vdots         & \vdots            & \vdots     & \vdots & \vdots\\
$\widetilde c_s$ & $\widetilde a_{s1}$ & \ldots     & $\widetilde a_{s,s-1}$ & 0\\
 \hline
               & $\widetilde b_1$  & \ldots & $\widetilde b_{s-1}$ & $\widetilde b_s$\vphantom{\Huge M}\\
\end{tabular}
%
\qquad
\begin{tabular}{c|cccc}
$c_1$  & $a_{11}$ & 0        & \ldots & 0\\
$c_2$  & $a_{21}$ & $a_{22}$ & \ldots & 0\\
\vdots & \vdots   & \vdots   & \vdots & \vdots\\
$c_s$  & $a_{s1}$ & $a_{s2}$ & \ldots & $a_{ss}$\\
 \hline
       & $b_1$    & $b_2$    & \ldots & $b_{s}$\vphantom{\Huge M}\\
\end{tabular}
\end{center}
Для явно-неявной схемы Эйлера $\widetilde a=0$, $\widetilde b=1$, $\widetilde c=0$ (явная схема, forward Euler) и $a=1$,
$b=1$, $c=1$ (неявная схема, backward Euler). В частности, имеют место следующие таблицы:

--- схема IMEX-SSP2(2,2,2)
\begin{center}
\begin{tabular}{c|cc}
0 & 0 & 0\\
1 & 1 & 0\\
\hline
  & 1/2  & 1/2\\
\end{tabular}
\qquad
\begin{tabular}{c|cc}
$\gamma$ & $\gamma$ & 0\\
$1-\gamma$ & $1-2\gamma$ & $\gamma$\\
\hline
  & 1/2  & 1/2\\
\end{tabular}
\end{center}

--- схема IMEX-SSP2(3,3,2)
\begin{center}
\begin{tabular}{c|ccc}
0 & 0 & 0 & 0\\
1 & 1 & 0 & 0\\
1/2 & 1/4 & 1/4 & 0\\
\hline
  & 1/6 & 1/6 & 2/3\\
\end{tabular}
\qquad
\begin{tabular}{c|ccc}
$\gamma$ & $\gamma$ & 0 & 0\\
$1-\gamma$ & $1-2\gamma$ & $\gamma$ & 0\\
1/2 & $1/2-\gamma$ & $\gamma$ & 0\\
\hline
  & 1/6  & 1/6 & 2/3\\
\end{tabular}
\end{center}
Здесь $\gamma=1/-1/\sqrt{2}$.

Расчеты проводятся при помощи соотношения
\begin{align*}
y^{(r)}=\frac{1}{1+\Delta t a_{rr}\lambda}\widehat y^{(r-1)}+ \frac{\Delta t a_{rr}\lambda}{1+\Delta t
a_{rr}\lambda}M_{\widehat y^{(r-1)}},
\end{align*}
где
\begin{align*}
\widehat y^{(r-1)}=y^k+\Delta t\sum\limits_{l=1}^{r-1}\widetilde a_{rl}T[y^{(l)}]+\Delta
t\sum\limits_{l=1}^{r-1}a_{rl}Q[y^{(l)}].
\end{align*}
Для максвелловской функции распределения $M_{y^{(r)}}=M_{\widehat y^{(r-1)}}$.

\subsection{Полностью дискретная схема}

Оператор переноса и столкновительный оператор представляются в следующем дискретном виде
\begin{align*}
&T^k_{i,j}=-\frac{\widehat F^k_{i+1/2,j}-\widehat F^k_{i-1/2,j}}{h_x};\\[0.2cm]
&Q^k_{i,j}=\lambda\left(M^k_{i,j}-f^k_{i,j}\right).
\end{align*}

Макроскопические свойства среды находятся из соотношений
\begin{align*}
&n^k_i=h_v\sum\limits_{j=1}^{N^v}f^k_{i,j};\\
&u^k_i=\frac{h_v}{n^k_i}\sum\limits_{j=1}^{N^v}v_j f^k_{i,j};\\
&\theta^k_i=\frac{h_v}{n^k_i}\sum\limits_{j=1}^{N^v}|v_j|^2 f^k_{i,j}-|u^k_i|^2.
\end{align*}

Равновесная функция распределения в узлах сетки находится из соотношения
\begin{align*}
M^k_{i,j}=\frac{n^k_i}{(2\pi\theta^k_i)^{1/2}}\exp\left[-\frac{(v_j-u^k_i)^2}{2\theta^k_i}\right].
\end{align*}

Разностная схема Рунге--Кутта для уравнения BGK имеет вид
\begin{align*}
&\widehat f^{(r-1)}_{i,j}=f^k_{i,j}+\Delta t\sum\limits_{l=1}^{r-1}\widetilde a_{rl}T^{(l)}_{i,j}+\Delta t\sum\limits_{l=1}^{r-1}a_{rl}Q^{(l)}_{i,j};\\
%
&f^{(r)}_{i,j}=\frac{1}{1+\Delta t a_{rr}\lambda}\widehat f^{(r-1)}_{i,j}+\frac{\Delta t a_{rr}\lambda}{1+\Delta t a_{rr}\lambda}\widehat M^{(r-1)}_{i,j};\\
%
&f^{k+1}_{i,j}=f^k_{i,j}+\Delta t\sum\limits_{r=1}^{s}\widetilde b_r T^{(r)}_{i,j}+\Delta t\sum\limits_{r=1}^{s}b_r
Q^{(r)}_{i,j}.
\end{align*}
При этом $\widehat M^{(r-1)}_{i,j}$ представляет собой равновесную функцию распределения, рассчитанную в соответствии с
дискретной аппроксимацией $\widehat f^{(r-1)}_{i,j}$.


\section{Тестовые примеры (релаксация)}

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

\subsection{Пример 1}

Решение уравнения BGK ищется в пространственно-однородном случае в скоростной области $v\in[-6,+6]$ на интервале времени
$(0,T]$, где $T=10$. Сетка состоит из $N_v=32$ узлов. Шаг по времени составляет $\Delta t=0.02$. Полагается, что
$\lambda=1$.

Начальное условие для функции распределения в момент времени $t=0$ имеет вид
\begin{align*}
f(0,v)=\frac{v^2}{\pi}\exp\left(-v^2\right).
\end{align*}

Точное решение уравнения BGK имеет вид
\begin{align*}
f(t,v)=\frac{1}{2\pi S^2}\exp\left(-\frac{v^2}{2S}\right)\left(2S-1+\frac{1-S}{2S}v^2\right),
\end{align*}
где $S=1-0.5\exp(-t/8)$.

Функцию распределения в различные моменты времени показывает рис.~\ref{has00}. Сплошные линии соответствуют численному
решению задачи, а значки~--- точному решению. В конечный момент времени функция распределения является максвелловской.

\begin{i-fig}{has00}{0.5}{
Функция распределения в моменты времени $t=1$ (1, $\circ$), 2 (2, $\bullet$), 10 (3, $\square$)
}\end{i-fig}

\subsection{Пример 2}

Решение уравнения BGK ищется в пространственной области $x\in[-5,+5]$, скоростной области $v\in[-5,+5]$ на интервале времени
$(0,T]$, где $T=1$. Сетка состоит из $N_x=128$ и $N_v=130$ узлов. Полагается, что $\lambda=1$. Конвективным переносом
пренебрегается.

Начальное условие для функции распределения в момент времени $t=0$ имеет вид (bump function)
\begin{align*}
f(0,x,v)=f_0(x,v)=b(v)g(x),
\end{align*}
где $g(x)=\exp(-|x|)$, а функция $b(v)$ определяется на бесконечном интервале, но принимает ненулевые значения только в
области $v\in(-2,+2)$. Начальное условие принимает вид
\begin{align*}
f(0,x,v)=\left\{\begin{array}{ll}
\exp(-|x|)\dfrac{5(v^2+4)h(v)}{(v^2-4)[1+h(v)]^2} &\text{при}\,\,v\in(-2,2),\\
0 &\text{при}\,\,v\notin(-2,2).
\end{array}\right.
\end{align*}
Здесь $h(v)=\exp[5v/(v^2-4)]$. Начальную функции распределения показывает рис.~\ref{has01}.

\begin{i-fig}{has01}{0.5}{
Функция распределения в начальный момент времени }\end{i-fig}

Точное решение уравнения BGK имеет вид
\begin{align*}
f(t,x,v)=f(0,x,v)\exp\left(-\lambda t\right)+M_f(0,x,v)\left[1-\exp\left(-\lambda t\right)\right].
\end{align*}

Вычисляя моменты, получим, что $n=g(x)$ (концентрация частиц). Функция $b(v)$ является четной, а функция $v b(v)$~---
нечетной. Интеграл от нечетной функции по симметричному интервалу обращается в ноль, поэтому $u=0$ (скорость потока).
Распределение энергии (температуры) представляется в приближенном виде $E=0.3713 g(x)$ (интегралы не вычисляются в конечном
виде). Функции распределения, соответствующие точному и численному решениям задачи, показывает рис.~\ref{has02} в конечный
момент времени. С течением времени происходит размытие начального распределения (конвективный перенос отсутствует).

\begin{i-fig}{has02}{0.9}{
Функции распределения, соответствующие точному (а) и численному (б) решениям задачи }\end{i-fig}

Относительная погрешность находится как разность между точным и численным решением (рис.~\ref{has03}). Сравнение численного
и точного решений показывает, что максимальная относительная ошибка не превосходит $2.5\cdot 10^{-5}$.

\begin{i-fig}{has03}{0.5}{
Относительная погрешность расчета функции распределения в конечный момент времени }\end{i-fig}

\subsection{Пример 3}

Решение уравнения BGK ищется в пространственной области $x\in[-6,+6]$, скоростной области $v\in[-6,+6]$ на интервале времени
$(0,T]$, где $T=1$. Сетка состоит из $N_x=64$ и $N_v=256$ узлов. Полагается, что $\lambda=1$. Конвективным переносом
пренебрегается.

Начальное условие для функции распределения в момент времени $t=0$ представляется в виде суммы четырех гауссовых функций,
центры которых расположены в точках $(x,v)=(\pm 2,\pm 2)$. Начальную функции распределения показывает рис.~\ref{has04}.

\begin{i-fig}{has04}{0.5}{
Функция распределения в начальный момент времени }\end{i-fig}

Функции распределения, соответствующие точному и численному решениям задачи, показывает рис.~\ref{has05} в конечный момент
времени. С течением времени происходит размытие начального распределения (конвективный перенос отсутствует).

\begin{i-fig}{has05}{0.9}{
Функции распределения, соответствующие точному (а) и численному (б) решениям задачи }\end{i-fig}

Относительная погрешность находится как разность между точным и численным решением (рис.~\ref{has06}). Сравнение численного
и точного решений показывает, что максимальная относительная ошибка не превосходит $2\cdot 10^{-3}$.

\begin{i-fig}{has06}{0.5}{
Относительная погрешность расчета функции распределения в конечный момент времени }\end{i-fig}


\section{Течение в ударной трубе}

Решение уравнения BGK ищется в пространственно-временной области $\Omega\times V=[-0.5,+0.5]\times[-10,+10]$ на интервале
времени $(0,T]$, где $T=0.2$. Сетка состоит из $N_x=256$ и $N_v=258$ узлов. Полагается, что $\lambda=10^4$ ($\tau=10^{-4}$).
Конвективным переносом пренебрегается.

Для задачи Сода (течение в ударной трубе) начальные распределения плотности, скорости и температуры задаются в виде:

--- при $x\in(-0.5,0)$
\begin{align*}
\left(\begin{array}{c}
n_L\\
u_L\\
\theta_L
\end{array}\right)=
\left(\begin{array}{c}
1.0\\
0.0\\
1.0
\end{array}\right);
\end{align*}

--- при $x\in(0,+0.5)$
\begin{align*}
\left(\begin{array}{c}
n_R\\
u_R\\
\theta_R
\end{array}\right)=
\left(\begin{array}{c}
0.125\\
0.0\\
0.8
\end{array}\right).
\end{align*}
Давление и внутренняя энергия находятся из соотношений: $p=n\theta$, $e=0.5\theta$.

Начальное распределение функции распределения задается в виде:

--- при $x\in [-0.5,0]$, $v\in[-10,+10]$
\begin{align*}
f_0(x,v)=\frac{n_L}{(2\pi\theta_L)^{1/2}}\exp\left(-\frac{|v-u_L|^2}{2\theta_L}\right);
\end{align*}

--- при $x\in (0,+0.5]$, $v\in[-10,+10]$
\begin{align*}
f_0(x,v)=\frac{n_R}{(2\pi\theta_R)^{1/2}}\exp\left(-\frac{|v-u_R|^2}{2\theta_R}\right).
\end{align*}

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

Решение уравнения BGK при $\lambda\to\infty$ стремится к решению уравнений Эйлера (функция распределения становится
максвелловской, а ее моменты дают распределения плотности, скорости и температуры). Моменты функции распределения находятся
путем интегрирования по усеченной скоростной области. Линии уровня начальной ($t=0$) и конечной ($t=0.2$) функции
распределения показывает рис.~\ref{has07}.

\begin{i-fig}{has07}{0.9}{
Функция распределения в начальный $t=0$ (а) и конечный $t=0.2$ (б) моменты времени }\end{i-fig}

Сравнение распределений максроскопических характеристик потока, полученных из решения уравнения BGK и уравнений Эйлера,
показывает рис.~\ref{has08}. Наблюдается достаточно хорошее согласование решений, полученных при помощи различных моделей.

\begin{i-fig}{has08}{0.9}{
Распределения плотности (а), скорости (б) и внутренней энергии (в) в конечный момент времени. Линии 1 соответствуют решению
уравнения BGK, а линии 2~--- решению уравнений Эйлера }\end{i-fig}

Другие расчеты проводятся при различных числах Кнудсена. При использовании модели твердых сфер показатель степени в
зависимости динамической вязкости от температуры равняется $\omega=1/2$, а средняя длина свободного пробега молекул
находится из соотношения
\begin{align*}
\lambda_0=\frac{16}{5}\frac{\mu_0}{p_0}\left(\frac{R\theta_0}{2\pi}\right)^{1/2},
\end{align*}
где $\mu$~--- динамическая вязкость, $R$~--- газовая постоянная. Индекс $0$ относится к характерным параметрам. В качестве
характерных параметров задачи принимаются параметры, соответствующие левой части расчетной области. Учитывая, что длина
расчетной области $L=1$, получим $\text{Kn}=\lambda_0$. Изменение динамической вязкости от $10^{-5}$ до $10$ соответствует
изменению числа Кнудсена в нервале от $1.28\cdot 10^{-5}$ до $12.8$. При этом число Прандтля принимает значение
$\text{Pr}=2/3$.

Линии уровня функции распределения в конечный момент времени показывает рис.~\ref{has09}.

\begin{i-fig}{has09}{0.9}{
Функция распределения в конечный момент времени при $\text{Kn}=0.1$ (а) и $\text{Kn}=10$ (б) }\end{i-fig}

Распределения максроскопических характеристик потока, полученных из решения уравнения BGK при различных числах Кнудсена,
показывает рис.~\ref{has10}. Результаты численных расчетов для различных чисел Кнудсена достаточно хорошо согласуются с
данными работ \cite{Guo15, Zhu17, Guo21}.

\begin{i-fig}{has10}{0.9}{
Распределения плотности (а), скорости (б) и внутренней энергии (в) в конечный момент времени при $\text{Kn}=0.1$ (линии 1) и
$\text{Kn}=10$ (линии 2)
}\end{i-fig}


\section{Заключение}

Разработан метод конечных объемов для решения кинетического уравнения Больцмана с интегралом столкновений в виде BGK. Для
дискретизации по времени используется явно-неявная схема Рунге--Кутты (простран\-ствен\-но-однородный столкновительный
этап). Численные потоки находятся при помощи схемы TVD типа с ограничителем minmod. Следует отметить, что данный метод не
обеспечивает кинетическую (микроскопическую) консервативность, связанную с выполнением дискретных аналогов условий
сохранения для интеграла столкновений.

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

Исследование выполнено в рамках научной программы Национального центра физики и математики (проект <<Математическое
моделирование на суперЭВМ экса-и зеттафлопсной производительности>>).


\bigskip
\begin{thebibliography}{99} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).
\bibitem{Bhatnagar54}
P.L. Bhatnagar , E.P. Gross, M. Krook, {\it Model for collision processes in gases. I. Small amplitude processes in charged
and neutral one-component systems}, Physical Review, \textbf{94}:3 (1954), 511--525.

\bibitem{Bobylev95}
A.V. Bobylev, A. Palczewski, J. Schneider, {\it On approximation of the Boltzmann equation by discrete velocity models},
Comptes Rendus de l'Academie des Sciences, \textbf{320} (1995), 639--644.

\bibitem{Mouhot13}
C. Mouhot, L. Pareschi, T. Rey, {\it Convolutive decomposition and fast summation methods for discrete-velocity
approximations of the boltzmann equation}, Mathematical Modelling and Numerical Analysis, \textbf{47}:5 (2013), 1515--1531.

\bibitem{Hu20}
J. Hu, K. Qi, {\it A fast Fourier spectral method for the homogeneous Boltzmann equation with non-cutoff collision kernels},
Journal of Computational Physics, \textbf{423} (2020), 109806.

\bibitem{Wu13}
L. Wu, C. White, T.J. Scanlon, J.M. Reese, Y. Zhang, {\it Deterministic numerical solutions of the Boltzmann equation using
the fast spectral method}, Journal of Computational Physics, \textbf{250} (2013), 27--52.

\bibitem{Boelens20}
A.M.P. Boelens, D. Venturi, D.M. Tartakovsky, {\it Tensor methods for the Boltzmann-BGK equation}, Journal of Computational
Physics, \textbf{421} (2020), 109744.

\bibitem{Pareschi96}
L. Pareschi,, B. Perthame, {\it A Fourier spectral method for homogeneous Boltzmann equation}, Transport Theory and
Statistical Physics, \textbf{25} (1996), 369--382.

\bibitem{Alekseenko19}
A. Alekseenko, J. Limbacher, {\it Evaluating high order discontinuous Galerkin discretization of the Boltzmann collision
integral in $O(N^2)$ operations using the discrete Fourier transform}, Kinetic and Related Models, \textbf{12} (2019), 703.

\bibitem{Bobylev97}
A.V. Bobylev, S. Rjasanow, {\it Difference scheme for the Boltzmann equation based on fast Fourier transformation}, European
Journal of Mechanics. B/Fluids, \textbf{16} (1997), 293--306.

\bibitem{Watch09}
C. Watchararuangwita, Y.N. Grigoriev, S.V. Meleshko, {\it A deterministic spectral method for solving the Boltzmann equation
for one-dimensional flows}, Science Asia, \textbf{35} (2009), 70--79.

\bibitem{Bobylev99}
A.V. Bobylev, S. Rjasanow, {\it Fast deterministic method of solving the Boltzmann equation for hard spheres}, European
Journal of Mechanics. B/Fluids, \textbf{18} (1999), 869--887.

\bibitem{Ibragimov02}
I. Ibragimov, S. Rjasanow, {\it Numerical solution of the Boltzmann equation on the uniform grid}, Computing, \textbf{69}
(2002), 163--186.

\bibitem{Holloway15}
I. Holloway, A. Wood, A. Alekseenko, {\it Acceleration of Boltzmann collision integral calculation using machine learning},
Mathematics, \textbf{9} (2021), 1384.

\bibitem{Dimarco13}
G. Dimarco, L. Pareschi, {\it Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations},
SIAM Journal of Numerical Analysis, \textbf{51}:2 (2013), 1064--1087.

\bibitem{Dimarco18}
G. Dimarco, R. Loubere, J. Narski, T. Rey, {\it An efficient numerical method for solving the Boltzmann equation in
multidimensions}, Journal of Computational Physics, \textbf{353} (2018), 46--81.

\bibitem{Coron91}
F. Coron, B. Perthame, {\it Numerical passage from kinetic to fluid equations}, SIAM Journal of Numerical Analysis,
\textbf{28} (1991), 26--42.

\bibitem{Pieraccini07}
S. Pieraccini, G. Puppo, {\it Implicit-explicit schemes for BGK kinetic equations}, Journal of Scientific Computing,
\textbf{32} (2007), 1--28.

\bibitem{Perthame89}
B. Perthame, {\it Global existence to the BGK model of Boltzmann equation}, Journal of Differential Equations, \textbf{82}
(1989), 191--205.

\bibitem{Mischler96}
S. Mischler, {\it Uniqueness for the BGK-equation in $R^N$ and rate of convergence for a semi-discrete scheme}, Differential
and Integral Equations, \textbf{9}:5 (1996), 1119--1138.

\bibitem{Nassios13}
J. Nassios, J.E. Sader, {\it High frequency oscillatory flows in a slightly rarefied gas according to the Boltzmann–BGK
equation}, Journal of Fluid Mechanics, \textbf{729} (2013), 1--46.

\bibitem{Shakhov68}
E.M. Shakhov, {\it On the generalization of the Krook kinetic equation}, Fluid Dynamics, \textbf{5} (1968), 142--145.

\bibitem{Holway66}
H. Holway, {\it New statistical models for kinetic theory: methods of construction}, Physics of Fluids, \textbf{9}:9 (1966),
1658--1673.

\bibitem{Mieussens04}
L. Mieussens, H. Struchtrup, {\it Numerical comparison of Bhatnagar--Gross--Krook models with proper Prandtl number},
Physics of Fluids, \textbf{16} (2004), 2797--2813.

\bibitem{Mieussens00b}
L. Mieussens, {\it Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric
geometries}, Journal of Computational Physics, \textbf{162}:2 (2000), 429--466.

\bibitem{Mieussens00a}
L. Mieussens, {\it Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics}, Mathematical
Models and Methods in Applied Sciences, \textbf{10} (2000), 1121--1149.

\bibitem{Garrett18}
C.K. Garrett, C.D. Hauck, {\it A fast solver for implicit integration of the Vlasov--Poisson system in the Eulerian
framework}, SIAM Journal of Scientific Computing, \textbf{40}:2 (2018), B483--B506.

\bibitem{Gottlieb01}
S. Gottlieb, C.-W. Shu, E. Tadmor, {\it Strong stability-preserving high-order time discretization methods}, SIAM Review,
\textbf{43}:1 (2011), 89--112.

\bibitem{Michel22}
V. Michel-Dansac, A. Thomann, {\it TVD-MOOD schemes based on implicit explicit time integration}, Applied Mathematics and
Computation, \textbf{433} (2022), 127397.

\bibitem{Guo15}
Z.L. Guo, R.J. Wang, K. Xu, {\it Designing a consistent implementation of the discrete unified gas-kinetic scheme for the
simulation of three-dimensional compressible natural convection}, Physical Review E, \textbf{91} (2015), 033313.

\bibitem{Zhu17}
L. Zhu, S. Chen, Z. Guo, {\it dugksFoam: an open source OpenFOAM solver for the Boltzmann model equation}, Computer Physics
Communications, \textbf{213} (2017), 155--164.

\bibitem{Guo21}
Z. Guo, K. Xu, {\it Progress of discrete unified gas-kinetic scheme for multiscale flows}, Advances in Aerodynamics,
\textbf{3}:6 (2021), 1--42.
\end{thebibliography}

\end{document}
