keywords. Stiff kinetic equations, BGK/ES-BGK models, IMEX Runge–Kutta methods, Compressible Euler equations, Compressible Navier–Stokes equations

1 Introduction↩︎

One of the most well-known kinetic models for rarefied gas dynamics is the Boltzmann transport equation (BTE). Its dimensionless form is written as \[\label{Beq} \partial_t f + v\cdot \nabla_x f = \frac{1}{\varepsilon}Q(f,f),\tag{1}\] where \(f(t,x,v)\) is the distribution function which depends on time \(t > 0\), on the position of particles \(x \in \mathbb{R}^{d_x}\) and on their velocity \(v \in \mathbb{R}^{d_v}\). Here \(\varepsilon\) is the so-called Knudsen number, defined as the ratio of the mean free path of molecules and the characteristic length scale of the physical problem. The dimension of space and velocity domains are denoted by \(d_x\) and \(d_v\), respectively.

The Boltzmann collision operator \(Q(f,f)\) is a non-linear operator that describes the binary collisions between molecules. It acts only on the velocity dependence of the distribution function \(f\), and has the following fundamental properties

• of conserving mass momentum and energy: \[\label{rel1} \left \langle Q(f,f)\phi(v) \right \rangle = \boldsymbol{0} \in \mathbb{R}^{d_v+2},\tag{2}\] for \(\phi(v) = \left(1, v, \frac{1}{2}|v|^2\right)^\top\), and \(\left \langle g \right \rangle := \int_{\mathbb{R}^{d_v}}g(v) dv\),

• to satisfy \(H\)-theorem: \[\label{rel2} \int_{\mathbb{R}^{d_v}}Q(f,f) \ln f dv \le 0,\tag{3}\]

• to vanish, i.e., \(Q(f,f)=0\) when \(f\) is the local Maxwellian: \[\label{Max} \mathcal{M}[f](t,x,v) = \frac{\rho(t,x)}{(2 \pi T(t,x))^{d_v/2}}\exp\left(-\frac{|v-u(t,x)|^2}{2T(t,x)}\right)\tag{4}\] where \(\rho\), \(u\), and \(E\) are density, mean velocity and energy associated to \(f\): \[\left \langle f \phi(v) \right \rangle = \left( \rho, \rho u, E\right)^\top = U\in \mathbb{R}^{d_v+2},\] and temperature \(T\) is given by \(\displaystyle \frac{d_v \rho T}{2}=E- \frac{1}{2}\rho |u|^2\). We introduce the vector \(U\) for later use.

When the Knudsen number is small, it is well known that BTE is closely related to fluid models such as compressible Euler or compressible Navier-Stokes (CNS) equations. The form of these fluid models associated to BTE are traditionally derived using the perturbation techniques like the Hilbert or Chapman-Enskog expansions [3]–[5]. Thus, BTE can be used for various Knudsen number, i.e., from rarefied to continuum gas dynamics.

In spite of its good predictability and close relationship with fluid models, computation of the Boltzmann collision operator is very expensive due to its high dimensionality. Furthermore, the problem becomes more severe as the Knudsen number gets closer to zero (fluid regime). In this case, solving BTE by a standard explicit numerical scheme requires the use of a time step of the order of \(\varepsilon\), which leads to very expensive numerical computations. Even if one adopts an implicit or semi-implicit time discretization for the collision part, it is still numerically challenging to construct an efficient implicit solver due to the complicate structure of the Boltzmann collision operator.

To circumvent the issue on computational cost of the Boltzmann collision operator \(Q(f,f)\) in (1 ), simpler kinetic models have been proposed to mimic the main properties of the full integral operator \(Q(f,f)\). One such model is the BGK model [6]: \[\label{BGKo} \frac{\partial f}{\partial t} + v \cdot \nabla_x f = {\frac{\tau}{\varepsilon}}\left(\mathcal{M}[f] -f \right),\tag{5}\] where \(\tau\) is the collision frequency that depends on \(\rho\) and \(T\). The BGK model replaces the Boltzmann collision operator with a simple relaxation toward the local Maxwellian. Note that the BGK model still maintains the conservation of mass, momentum, and energy, as well as the entropy inequality [4]. In addition, the BGK model describes the correct fluid limit as \(\varepsilon \to 0\), i.e., at the leading order term \(\varepsilon=0\), it yields the compressible Euler equations (see [3], [4]). Unfortunately, at the first-order correction in \(\varepsilon\), the transport coefficients obtained at the Navier–Stokes level are not satisfactory. In particular, the Prandtl number defined by \(\displaystyle \textrm{Pr} = \frac{\gamma}{\gamma -1}\frac{\mu}{\kappa},\) which relates the viscosity \(\mu\) to the heat conductivity \(\kappa\) of gases is fixed by \(1\). Here, the polytropic constant \(\gamma\) for monatomic molecule with transional motions is given by \(\displaystyle\gamma=\frac{d_v+2}{d_v}\). Note that for most realistic gases, we have \(Pr < 1\). In addition, in the hard-sphere model for monoatomic gases (\(\gamma = 5/3\)) in Boltzmann equation, its Prandtl number is very close to \(2/3\).

Many variants of the BGK model have been proposed in order to give the correct transport coefficeints at the Navier–Stokes level. In [7], Holway proposed a model that possesses several desirable properties. It not only satisfies the correct conservation laws and the entropy condition, but also yields the Navier-Stokes equations with a Prandtl number less than one through the Chapman-Enskog expansion. This model is known as the ellipsoidal statistical model (ES-BGK) and reads: \[\begin{align} \label{ES-BGK} \begin{aligned} &\frac{\partial f}{\partial t} + v \cdot \nabla_x f = \frac{\tau}{\varepsilon}(\mathcal{G}[f] - f), \end{aligned} \end{align}\tag{6}\] where the collision frequency \(\tau =\frac{\rho T}{\mu(1-\nu)}\) depends on the free parameter \(-\frac{1}{2} \le \nu < 1\). The anisotropic Gaussian distribution \(\mathcal{G}[f]\) is defined by \[\label{Gauss} \mathcal{G}[f](t,x,v) = \frac{\rho(t,x)}{\sqrt{\textrm{det}(2\pi \mathcal{T}(t,x))}}\exp\left(-\frac{(v-u(t,x))^\top\mathcal{T}(t,x)^{-1}(v-u(t,x))}{2}\right).\tag{7}\] The temperature tensor \(\mathcal{T}(t,x)\) and stress tensor \(\Theta(t,x)\) are defined as \[\label{TTT} \mathcal{T}(t,x) = (1-\nu)T(t,x) Id + \nu \Theta(t,x),\tag{8}\] \[\label{rhoTheta} \Theta(t,x) = \frac{1}{\rho}\int_{R^{d_v}}(v-u) \otimes (v-u) f(t,x,v)dv = \frac{1}{\rho}\left \langle (v-u) \otimes (v-u) f(v) \right \rangle,\tag{9}\] respectively, and \(Id\) is the \(d_v\times d_v\) identity matrix.

As emphasized before, ES-BGK model also has a close relationship with fluid models. The Chapman-Enskog expansion applied to ES-BGK model gives the compressible Euler equations [8] at the leading order term: \[\label{Eq:Euler} \partial_t \left( \begin{array}{l} \rho\\ \rho u \\ E \end{array} \right) + \nabla_x \cdot\left( \begin{array}{l} \rho u \\ \rho u \otimes u + p Id\\ (E + p ) u \end{array}\right) = 0,\tag{10}\] where \(p=\rho T\). While at the first-order correction in \(\varepsilon\), we obtain the CNS equations [4] (for details see Appendix 6): \[\begin{align} \label{NSeq} \partial_t \left( \begin{array}{l} \rho\\ \rho u \\ E \end{array} \right) + \nabla_x\left( \begin{array}{l} \rho u \\ \rho u \otimes u + p Id\\ (E + p ) u \end{array}\right)=\left( \begin{array}{l} 0 \\ \varepsilon \nabla_x\cdot (\mu\sigma(u))\\ \varepsilon \nabla_x\cdot( \mu\sigma(u)u + q) \end{array}\right) \end{align}\tag{11}\] where the stress tensor \(\sigma(u)\) and heat flux \(q\) are given by \[\sigma(u) = \nabla_x u + (\nabla_x u)^\top - \frac{2}{d_v} \nabla_x \cdot uId,\quad q = \kappa \nabla_x T,\] with viscosity \(\mu = \frac{p}{(1-\nu)\tau}\) and thermal conductivity \(\kappa = \frac{d_v + 2}{2}\frac{p}{\tau}.\) Thus, the Prandtl number is given by \[\frac{2}{3} \leq Pr = \frac{d_v+2}{2}\frac{\mu}{\kappa} = \frac{1}{1 - \nu} < \infty,\] and this implies that desired such number can be recovered by choosing \(\nu\) appropriately. Note that the equation for the total energy can be replaced by the equation for the temperature \(T\), [9] \[\label{Temperature} \frac{d_v}{2}\rho \left( \partial_t T + u \cdot \nabla_x T \right) + \rho T \nabla_x \cdot u = \mathcal{O}(\varepsilon).\tag{12}\]

Considering the relationship with fluid models, when developing numerical methods for ES-BGK model, such methods should have a correct asymptotic behavior, i.e., for small parameter \(\varepsilon\), the schemes should degenerate into a good approximation of the fluid asymptotic (compressible Euler or CNS equations).

In the fluid dynamic limit, i.e., in the case of the limiting compressible Euler equations, numerical methods that work effective while keeping the mesh size and time step fixed as the Knudsen number approaches zero, are referred to as asymptotic preserving (AP) [10]. Additionally, a scheme that not only preserves the correct asymptotic behavior but also maintains high accuracy in the fluid dynamic regime is called asymptotically accurate (AA) [11]. In other words, a scheme that satisfies both the AP and AA properties exhibits robustness and high accuracy in the fluid dynamic regime, without resolving the small Knudsen number.

On the construction of AP and AA methods for kinetic models, IMEX time discretization [12]–[17] have been successfully applied and studied.

We focus our literature reviews to the works on the NS limit based on IMEX time discretization applied to BGK or ES-BGK model. There is, of course, a considerable amount of literature on the use of IMEX-RK methods for BGK-type equations, (see, e.g.,[8], [11], [18]). For instance, in [19], [20], the authors considered a micro-macro decomposition of the BGK equation and then applied IMEX-RK schemes to the resulting coupled system. In [21], authors introduced a first order IMEX method for ES-BGK model and showed that it is consistent to first order time discretization of CNS equations. High order IMEX-RK methods of type CK (or type II) [1] and IMEX multistep methods [18] are also similarly applied to ES-BGK model, and the methods are shown to be able to capture the NS limit under suitable conditions. However, reference [1] focuses solely on a specific subclass of IMEX schemes (namely, type CK or type II methods, [11]).

This paper provides a broader and more unified analysis of the various IMEX-RK schemes, specifically type I and type II, extending and improving upon existing results in the literature for the ES-BGK model. In particular, we study their asymptotic behavior and generalize the results presented in [1], showing that both types of schemes are capable of capturing the compressible Navier-Stokes (CNS) limit without resolving the small parameter \(\varepsilon\).

Another novelty of the present work is that we numerically demonstrate the uniform accuracy of IMEX-RK methods applied to ES-BGK model for a wide range of Knudsen numbers presented in [2]. These IMEX-RK schemes were originally developed and analyzed for a specific class of hyperbolic relaxation systems. In the context of the ES-BGK model, we derive additional order conditions (see Theorem 11) that ensure uniform accuracy for schemes of type II. We show that, under suitable assumptions on the coefficients of such schemes, the considered schemes satisfy these additional conditions and maintain the order of accuracy throughout the entire range of Knudsen numbers. These conditions are sufficient but not necessary. However, the price to pay is that the IMEX-RK schemes proposed in [2] require more stages than those commonly used.

The outline of this paper is the following. In Section 2, we explain IMEX-RK methods for ES-BGK model. In Section 3, we perform asymptotic analysis at the Navier-Stokes level. Then in Section 4, we give several numerical tests in order to validate our theoretical findings. Finally, conclusion will be provided.

2 IMEX-RK schemes for ES-BGK equation↩︎

For the time discretization of 6 we consider an implicit-explicit (IMEX) Runge-Kutta (RK) scheme because the convection term in 6 is not stiff and the collision term is stiff when Knudsen number is small.

An IMEX-RK scheme usually can be represented by the double Butcher tables: \[\label{BT} \begin{array}{c|c} \tilde{c} & \tilde{A}\\ \hline \\[-0.3cm] & \tilde{b}^\top \end{array}, \quad \begin{array}{c|c} c & {A}\\ \hline \\[-0.3cm] & {b}^\top \end{array}.\tag{13}\] Here \(\tilde{A} = (\tilde{a}_{ij})\) with \(\tilde{a}_{ij} = 0\) for \(j \ge i\) and \(A = ({a}_{ij})\) with \({a}_{ij} = 0\) for \(j > i\) are \(s \times s\) matrices, which are associated to the explicit and implicit time discretizations, respectively. The coefficients vectors \({\tilde{b}} = (\tilde{b}_1, ..., \tilde{b}_s)^\top\) and \({{b}} = ({b}_1, ..., {b}_s)^\top\) represent the weights, and the vectors \({\tilde{c}} = (\tilde{c}_1, ..., \tilde{c}_s)^\top\) and \({{c}} = ({c}_1, ..., {c}_s)^\top\) are the nodes defined as \[\label{c95relation} \tilde{c}_i = \sum_{j = 1}^{i-1} \tilde{a}_{ij}, \quad {c}_i = \sum_{j = 1}^{i} {a}_{ij}, \quad i = 1,...,s.\tag{14}\]

Now we give some preliminary definitions. Based on the structure of matrix \(A\) in the implicit table, the IMEX schemes can be classified into the following types [11], [15], [16]:

We note that IMEX-RK schemes of type II are very attractive because they allow some simplifying assumptions that make order conditions easier to treat, therefore permitting the construction of higher order IMEX-RK schemes [15], [17]. On the other hand, schemes of type I are more amenable to a theoretical analysis since the matrix \(A\) of the implicit scheme is invertible. Later, we start our analysis with the latter scheme, type I.

In the following, we will also make use of the following representation of the matrix \(\tilde{A}\) in the explicit part \[\tilde{A}= \left( \begin{array}{ll} 0 & 0 \\ \tilde{a}& \hat{\tilde{A}}\end{array} \right), \label{CK2bis}\tag{15}\] where \(\tilde{a}=(\tilde{a}_{21},\ldots,\tilde{a}_{s 1})^\top\in\mathbb{R}^{s-1}\) and \(\hat{\tilde{A}}\in\mathbb{R}^{(s-1)\times (s-1)}\). This representation of matrix \(\tilde{A}\) is useful for the analysis of IMEX-RK methods of type II.

Finally, we give the following definition [22], [23]:

The first condition in Definition 3 guarantees that an \(A\)-stable implicit tableau is \(L\)-stable³. The AP property of IMEX-RK schemes are strongly related to the \(L\)-stability of the implicit part of the scheme. Finally, if the IMEX-RK scheme is GSA we have \(f^{n+1} = f^{(s)}\), i.e., the numerical solution coincides with the last stage of the method.

Below, we review IMEX-RK methods of type I and type II applied to ES-BGK model 6 .

• IMEX-RK method of type I. Applying an IMEX-RK method of type I to 6 , we get in vector form \[\begin{align} \label{ES95BGK} \begin{aligned} &{\boldsymbol{F}}= f^n {\boldsymbol{e}} - \Delta t \tilde{A} L({\boldsymbol{F}}) + \frac{\Delta t}{\varepsilon} A {\boldsymbol{\bar{\tau}}}(\mathcal{G}[{\boldsymbol{F}}] - {\boldsymbol{F}}),\\ &f^{n+1} = f^{n} - \Delta t \tilde{b}^\top L({\boldsymbol{F}}) + \frac{\Delta t}{\varepsilon} b^\top {\boldsymbol{\bar{\tau}}}(\mathcal{G}[{\boldsymbol{F}}] - {\boldsymbol{F}}), \end{aligned} \end{align}\tag{16}\] where \({\boldsymbol{F}}= (f^{(1)}, f^{(2)},..., f^{(s)})^\top\), \(L({\boldsymbol{F}}) = (L(f^{(1)}),...,L(f^{(s)}))^\top\), being \(L(f^{(k)}) = v \cdot \nabla_x f^{(k)}\) for all \(k =1 ,..., s\), \(\mathcal{G}[{\boldsymbol{F}}] := (\mathcal{G}(f^{(1)}),...,\mathcal{G}(f^{(s)}))^\top\), \(\bar{{\boldsymbol{\tau}}} := diag(\tau^{(1)},..., \tau^{(s)})\) is a diagonal matrix, and \({\boldsymbol{e}}=(1,1,...,1)^\top\) is a unit vector of length \(s\).

  From the first equation in (16 ), we get \[\label{eps} \Delta t {\bar \tau} (\mathcal{G}[{\boldsymbol{F}}] - {\boldsymbol{F}}) = \varepsilon A^{-1}\left( {\boldsymbol{F}}- f^n{\boldsymbol{e}} + \Delta t \tilde{A} L({\boldsymbol{F}})\right)\tag{17}\] and inserting in the numerical solution, we obtain \[\label{GSAeq} f^{n+1} = (1- b^\top A^{-1}{\boldsymbol{e}})f^n + b^\top A^{-1}{\boldsymbol{F}}-\Delta t(\tilde{b}^\top- b^\top A^{-1}\tilde{A}) L({\boldsymbol{F}}).\tag{18}\] Therefore, for IMEX-RK schemes of type I the numerical solution is independent on \(\varepsilon\) and we are able to pass to the limit \(\varepsilon \to 0\) in 16 .

• IMEX-RK method of type II. Regarding the IMEX-RK method of type II applied to 6 , with the notations in (?? ) and (15 ), the scheme reads in vector form \[\begin{align} \label{ES95BGK295CK} \begin{aligned} & F^{(1)} = f^n,\\ &\hat{{\boldsymbol{F}}} = f^n \hat{{\boldsymbol{e}}} - \Delta t \tilde{a} L({f^n}) - \Delta t \hat{\tilde{A}} L({\hat{\boldsymbol{F}}}) + \Delta t a \frac{\tau^n}{\varepsilon}(\mathcal{G}[f^n] - f^n) +\Delta t \hat{A} \frac{\hat{{\boldsymbol{\bar{\tau}}}}}{\varepsilon}(\mathcal{G}[\hat{{\boldsymbol{F}}}] - \hat{{\boldsymbol{F}}}),\\ &f^{n+1} = f^{n} -\Delta t \tilde{b}_1 L({f^n}) - \Delta t \hat{\tilde{b}}^\top L(\hat{{\boldsymbol{F}}}) + \Delta t b_1 \frac{\tau^n}{\varepsilon}(\mathcal{G}[f^n] - f^n) + \Delta t \hat{b}^\top \frac{\hat{{\boldsymbol{\bar{\tau}}}}}{\varepsilon}(\mathcal{G}[{\boldsymbol{F}}] - {\boldsymbol{F}}), \end{aligned} \end{align}\tag{19}\] with \(\hat{{\boldsymbol{e}}}=(1,1,...,1)^\top\) a unit vector of length \(s-1\). From the second equation in (19 ) we get \[\Delta t \hat{\bar{\tau}}(\mathcal{G}[\hat{{\boldsymbol{F}}}] - \hat{{\boldsymbol{F}}})=\varepsilon \hat{A}^{-1}\left[\hat{\boldsymbol{F}}-f^{n}\hat{\boldsymbol{e}}+\Delta t \tilde{a} L(f^{n})+\Delta t \hat{\tilde{A}}L(\hat{\boldsymbol{F}})\right]- \Delta t \hat{A}^{-1}a {\tau^n}(\mathcal{G}[f^n] - f^n). \label{eq:ll1}\tag{20}\] and substituting into the numerical solution in (19 ), we obtain \[\begin{align} \label{numsolII} \begin{aligned} f^{n+1} &= (1-\hat{b}^\top\hat{A}^{-1}\hat{{\boldsymbol{e}}})f^{n} -\Delta t \left(\tilde{b}_1 - \hat{b}^\top\hat{A}^{-1}\tilde{a}\right) L({f^n}) - \Delta t \left(\hat{\tilde{b}}^\top - \hat{b}^\top\hat{A}^{-1}\hat{\tilde{A}}\right) L(\hat{{\boldsymbol{F}}}) \\ &+ \Delta t \left(b_1 - \hat{b}^\top\hat{A}^{-1}a\right)\frac{\tau^n}{\varepsilon}(\mathcal{G}[f^n] - f^n) +\hat{b}^\top\hat{A}^{-1} {\hat{\boldsymbol{F}}}. \end{aligned} \end{align}\tag{21}\] Unfortunately, here the numerical solution depends on \(\varepsilon\). Now, in order to be able to pass to the limit \(\varepsilon \to 0\), we can either require that the initial conditions are well-prepared⁴ [11], or impose that the additional condition \(b_1-\hat{b}^\top\hat{A}^{-1}a=0\) has to be satisfied, which allows us to take the limit \(\varepsilon \to 0\) in 19 . Note that the condition \(b_1-\hat{b}^\top\hat{A}^{-1}a=0\), is automatically satisfied if the scheme is SA, i.e., \(\hat{b}^\top \hat{A}^{-1}={\boldsymbol{\hat{e}}}^\top_{s-1}\), with \(\hat{{\boldsymbol{e}}}^\top_{s-1}=(0,\cdots,0,1)\) vector of length \(s-1\). Alternatively, this condition also holds for IMEX schemes of type ARS, where having \(a = 0\) and \(b_1 = 0\) condition \(b_1-\hat{b}^\top\hat{A}^{-1}a=0\) is also satisfied.

In both types, the methods require implicit solver for computing relaxation terms. However, the use of collision invariants allows us to treat the implicit terms explicitly. The details will be provided in Section 3.

3 Asymptotic Properties of the IMEX RK Schemes↩︎

In this section, we discuss in detail the asymptotic properties of the IMEX-RK schemes of type I and II with respect to the Euler and Navier-Stokes limits. We begin with a brief overview and proof of the results related to the Euler limit, highlighting the preservation of leading-order asymptotic. More details can be found in [1], [8], [11], [21]. Next, we will analyze the asymptotic properties of IMEX-RK schemes of type I and II in the Navier-Stokes limit.

3.1 Preserving the Euler Limit.↩︎

We start with an arbitrary IMEX-RK scheme written for each stage \(k\): \[\label{IMEXcomp} f^{(k)} = f^{n} - \Delta t\sum_{\ell=1}^{k-1}\tilde{a}_{k\ell} v \cdot \nabla_x f^{(\ell)} + \frac{\Delta t}{\varepsilon}\sum_{\ell=1}^{k}a_{k\ell}\tau^{(\ell)}\left(\mathcal{G}[f^{(\ell)}] - f^{(\ell)}\right), \quad k = 1,...,s\tag{22}\] with numerical solution \[\label{solnum} f^{n+1} = f^{n} - \Delta t\sum_{k=1}^{s}\tilde{b}_{k} v \cdot \nabla_x f^{(k)} + \frac{\Delta t}{\varepsilon}\sum_{k=1}^{s}{b}_{k}\tau^{(k)}\left(\mathcal{G}[f^{(k)}] - f^{(k)}\right).\tag{23}\] At every \(k\)th stage of (22 ), the computation of \(f^{(k)}\) requires implicit treatment of \(\tau^{(k)}\) and \(\mathcal{G}[f^{(k)}]\). This difficulty can be circumvented by approximating \(\tau^{(k)}\) and \(\mathcal{G}[f^{(k)}]\) explicitly. For this, we rewrite (22 ) as \[\label{fK} f^{(k)} = \frac{\varepsilon}{\varepsilon + \tau^{(k)}\Delta t a_{kk}} f_*^{(k)} + \frac{\Delta t \tau^{(k)} a_{kk}}{\varepsilon + \tau^{(k)}\Delta t a_{kk}}\mathcal{G}[f^{(k)}]\tag{24}\] with \[f_*^{(k)} = f^{n} - \Delta t\sum_{\ell=1}^{k-1}\tilde{a}_{k\ell} v \cdot \nabla_x f^{(\ell)} + \frac{\Delta t}{\varepsilon}\sum_{\ell=1}^{k-1}a_{k\ell} \tau^{(\ell)}\left(\mathcal{G}[f^{(\ell)}] - f^{(\ell)}\right).\] Taking the moments \(\langle \cdot, \phi \rangle := \int_{\mathbb{R}^{d_v}} \cdot\, \phi(v) dv\) with \(\phi(v) := \left( 1, v, |v|^2/2\right)^\top\) on both sides of (24 ), and using the conservation properties (82 ), the implicit part is canceled, i.e. \[\biggl\langle(\mathcal{G}[f^{(k)}] -f^{(k)})\phi \biggr\rangle = 0.\] Therefore, one obtains the macroscopic quantities \(U:= (\rho, \rho u, E)\) at every stage \(k\): \[\begin{align} \label{LimitEq951} \begin{aligned} U^{(k)} &= \int_{\mathbb{R}^{d_v}} \left(f^{n} - \Delta t\sum_{\ell=1}^{k-1}\tilde{a}_{k\ell} v \cdot \nabla_x f^{(\ell)} \right) \phi(v) dv = \biggl\langle f^n \phi \biggr\rangle - \Delta t\sum_{\ell=1}^{k-1}\tilde{a}_{k\ell} \biggl\langle v \cdot \nabla_x f^{(\ell)} \phi \biggr\rangle, \end{aligned} \end{align}\tag{25}\] and numerical solution \[\begin{align} \label{LimitEq952} \begin{aligned} U^{n+1} &= \biggl\langle f^n \phi \biggr\rangle - \Delta t\sum_{k=1}^{s}\tilde{b}_{k} \biggl\langle v \cdot \nabla_x f^{(k)} \phi \biggr\rangle. \end{aligned} \end{align}\tag{26}\] Note that using (25 ) we can obtain \(\rho^{(k)}\), \(u^{(k)}\) and \(T^{(k)}\) at every stage \(k\), which enable us to compute explicitly \(\tau^{(k)}\). To evaluate \(\mathcal{G}[f^{(k)}]\), however, we need to determine the stress tensor \(\Theta^{(k)}\) in (9 ). To achieve this, we first define the tensor \(\Sigma\) as follows \[\label{SigmaAAA} \Sigma = \int_{\mathbb{R}^{d_v}} {v} \otimes { v} f d { v} = \rho(\Theta + u \otimes u ),\tag{27}\] and consider the stage values: \[\label{SigmaAAA95bis} \Sigma^{(k)} = \int_{\mathbb{R}^{d_v}} {v} \otimes { v} f^{(k)} d { v} = \rho^{(k)}(\Theta^{(k)} + u^{(k)} \otimes u^{(k)} ).\tag{28}\] Note that \(\Theta^{(k)}\) can be obtained by computing \(\Sigma^{(k)}\). Next, we use (8 ) and 27 to get \[\rho \mathcal{T} = \rho[(1-\nu)T Id + \nu \Theta] = \rho(1-\nu)T Id + \nu \Sigma - \nu \rho u \otimes u,\] and combine this with \[\int_{\mathbb{R}^{d_v}} v \otimes v \, \mathcal{G}[f](v) dv = \rho(\mathcal{T} + u \otimes u),\] to obtain \[\begin{align} \label{identity} \int_{\mathbb{R}^{d_v}} v \otimes v \,\left( \mathcal{G}[f] -f\right)dv = (1-\nu)\left(\rho\left({T}Id + u \otimes u\right)\;- \Sigma \right). \end{align}\tag{29}\] Now, we multiply the scheme (22 ) by \(v \otimes v\) and integrate it over \(v\), and use the relation 29 to derive the equation for \(\Sigma^{(k)}\): \[\begin{align} \label{sigma95comp} \Sigma^{(k)} = \frac{\varepsilon }{\varepsilon + (1-\nu)\tau^{(k)}\Delta t a_{kk}}\Sigma^* + \frac{\Delta t (1-\nu)\tau^{(k)} a_{kk}}{\varepsilon + (1-\nu)\tau^{(k)}\Delta t a_{kk}} \rho^{(k)}\left({T}^{(k)} Id + u^{(k)} \otimes u^{(k)}\right), \end{align}\tag{30}\] where \[\Sigma^* = \Sigma^{n} - \Delta t \sum_{\ell = 1}^{k-1} \tilde{a}_{k, \ell} \nabla_x \biggl\langle v \otimes v \, v f^{(\ell)} \biggr\rangle + \frac{\Delta t}{\varepsilon}\sum_{\ell = 1}^{k-1} a_{k\ell}(1-\nu)\tau^{(\ell)}[ \rho^{(\ell)}(T^{(\ell)}Id + u^{(\ell)}\otimes u^{(\ell)}) - \Sigma^{(\ell)}].\] Then, the IMEX-RK scheme reads \[\begin{align} \label{SCHEME95imp} \begin{aligned} \Sigma^{(k)} & = \frac{\varepsilon }{\varepsilon + (1-\nu)\tau^{(k)}\Delta t a_{kk}}\Sigma^* + \frac{\Delta t (1-\nu)\tau^{(k)} a_{kk}}{\varepsilon + (1-\nu)\tau^{(k)}\Delta t a_{kk}} \rho^{(k)}\left({T}^{(k)} Id + u^{(k)} \otimes u^{(k)})\right),\\ f^{(k)} & = \frac{\varepsilon}{\varepsilon + \tau^{(k)}\Delta t a_{kk}} f_*^{(k)} + \frac{\Delta t \tau^{(k)} a_{kk}}{\varepsilon + \tau^{(k)}\Delta t a_{kk}}\mathcal{G}[f^{(k)}], \end{aligned} \end{align}\tag{31}\] and \[\begin{align} \label{SCHEME95imp95bis} \begin{aligned} \Sigma^{n+1} &= \Sigma^{n} - \Delta t\sum_{k=1}^{s}\tilde{b}_{k} \nabla_x \biggl\langle v \otimes v \, v f^{(k)} \biggr\rangle + \frac{\Delta t}{\varepsilon}\sum_{k=1}^{s}{b}_{k}\biggl\langle v \otimes v \, \left(\tau^{(k)}\left(\mathcal{G}[f^{(k)}] - f^{(k)}\right)\right)\biggr\rangle,\\ f^{n+1} &= f^{n} - \Delta t\sum_{k=1}^{s}\tilde{b}_{k} v \cdot \nabla_x f^{(k)} + \frac{\Delta t}{\varepsilon}\sum_{k=1}^{s}{b}_{k}\tau^{(k)}\left(\mathcal{G}[f^{(k)}] - f^{(k)}\right). \end{aligned} \end{align}\tag{32}\] Regarding the leading order limit, i.e., Euler equation, the following results for IMEX-RK schemes of type I and II have been proved in [1], [11].

Usually, to guarantee that an IMEX-RK scheme of type I or II remains robust and stable for any value of \(\varepsilon\), that is, performs correctly across both the rarefied and fluid regimes, assuming GSA condition may achieve this task, particularly in the fluid regime, without additional order conditions, [2], [11]. However, as noted in Remark 8, we can consider Type ARS schemes that are not GSA, but characterized by some assumptions on the coefficients of the scheme in order to guarantee both AP and AA properties.

We conclude this section by considering an explicit RK scheme characterized by \((\tilde{A}, \tilde{b})\) applied to the limit compressible Euler equations (10 ) with \(T = p/\rho\), which takes the form: \[\begin{align} \label{ExS2} \begin{aligned} \rho^{(k)} &= \rho^{n} - \Delta t \sum_{\ell = 1}^{k-1} \tilde{a}_{k\ell}\textrm{div}_x (\rho^{(\ell)}u^{(\ell)}),\\ (\rho u)^{(k)} &= (\rho u)^{n} - \Delta t \sum_{\ell = 1}^{k-1} \tilde{a}_{k\ell}\textrm{div}_x ((\rho u)^{(\ell)}\otimes u^{(\ell)} + \rho^{(\ell)} T^{(\ell)}I) ,\\ E^{(k)} &= E^{n} - \Delta t \sum_{\ell = 1}^{k-1}\tilde{a}_{k\ell} \textrm{div}_x ([E^{(\ell)} + \rho^{(\ell)}T^{(\ell)}] u^{(\ell)}), \end{aligned} \end{align}\tag{33}\] for \(k=1,...,s\), with numerical solution \[\begin{align} \label{ExS3} \begin{aligned} \rho^{ n+1} &= \rho^{n} - \Delta t \sum_{k = 1}^{s} \tilde{b}_k \nabla_x (\rho^{(k)}u^{(k)}) ,\\ (\rho u)^{n+1} &= (\rho u)^{n} - \Delta t \sum_{k = 1}^{s} \tilde{b}_k \nabla_x ((\rho u)^{(k)}\otimes u^{(k)} + \rho^{(k)} T^{(k)}I) ,\\ E^{n+1 } &= E^{n} - \Delta t \sum_{k = 1}^{s} \tilde{b}_k \nabla_x ([E^{(k)} + \rho^{(k)}T^{(k)}] u^{(k)}). \end{aligned} \end{align}\tag{34}\] Note that in (33 )-(34 ), to the leading-order, we can exchange the energy equation with the temperature one (12 ), so we get \[\begin{align} \label{Temp} T^{(k)} = T^n - \Delta t\sum_{\ell = 1}^{k-1} \tilde{a}_{k\ell}\left( u^{(\ell)}\cdot \nabla_x T^{(\ell)} + \frac{2}{d_v}T^{(\ell)}\textrm{div}_x u^{(\ell)} \right), \end{align}\tag{35}\] \[\begin{align} \label{Temp2} T^{n+1} = T^n -\Delta t \sum_{k = 1}^{s} \tilde{b}_{k}\left( u^{(k)}\cdot \nabla_x T^{(k)} + \frac{2}{d_v} T^{(k)}\textrm{div}_x u^{(k)}\right). \end{align}\tag{36}\]

3.2 Preserving the Navier-Stokes Limit.↩︎

In this section, we analyze the asymptotic behavior of IMEX-RK schemes of type I and II for ES-BGK equations (6 ), and prove that for small value of \(\varepsilon\), these schemes asymptotically capture the NS limit without resolving \(\varepsilon\), providing a numerical scheme for the corresponding CNS equations (11 ). At the Navier-Stokes level, (for small but non-negligible Knudsen numbers, \(\varepsilon \ll 1\)), the GSA condition is not crucial to ensure consistency with CNS equations for both types of IMEX-RK schemes. Therefore, for the subsequent analysis, we do not impose the GSA condition on the scheme.

IMEX-RK scheme of type I. We first analyze the IMEX-RK method of type I.

Note that the pair \((\tilde{A}, \tilde{b})\) is the explicit table of the IMEX-RK method of type I in 13 .

Theorem 9 implies that, to accurately capture the Navier-Stokes equation, it is necessary that the local truncation error satisfies \(L.T.E = o(\varepsilon)\), which holds true if \(\displaystyle \Delta t^p + \frac{\varepsilon^2}{\Delta t}= {o}(\varepsilon)\).

In the stiff case, i.e., \(\Delta t \gg \varepsilon\), the phenomenon of order reduction may occur, particularly in the worst case when \(\Delta t \approx \varepsilon\) (intermediate regime of \(\varepsilon\)). In this scenario, the classical order of the method can drop off, leading to a loss of accuracy in IMEX-RK schemes within this regime.

In 57 , in the worse case scenario \(\Delta t \approx \varepsilon\) the maximum error is \(\mathcal{O}(\Delta t)\) for \(p >1\). For instance, this result states that the a third-order schemes \(p = 3\) may exhibit only first-order accuracy in the worst case. However, this is a lower bound, and in practice, better behavior is often observed (see numerical results).

This phenomenon was investigated in detail in [2], [11], where the authors considered a prototype hyperbolic relaxation system and performed an asymptotic expansion up to \(\mathcal{O}(\varepsilon)\) for IMEX-RK methods of type I. They showed that, under additional order conditions, these schemes effectively reduce to explicit Runge-Kutta methods at the \(\mathcal{O}(\varepsilon)\) level. In particular, these methods were carefully designed to accurately match the \(\mathcal{O}(\varepsilon)\) terms up to a desired order \(p\), ensuring that effectively behave as explicit Runge-Kutta schemes of order \(p\) for the convection-diffusion equation when \(\varepsilon\) is small but not negligible. However, constructing high-order IMEX-RK schemes of type I requires additional stages to satisfy extra-order conditions. In this work, we do not consider high-order IMEX-RK schemes of type I that satisfy these additional conditions introduced in [2], as their construction is particularly challenging. Instead, we focus on type II schemes, which are easier to construct. In Section 4, we compare these methods with other IMEX-RK schemes of type I or II existing in literature, and we will show that the later schemes that do not satisfy the extra order conditions may suffer from order reduction in the intermediate regime of \(\varepsilon\).

Finally, considering the results in [2], [11], we state the following.

In this case, viscosity and thermal conductivity \({\bar \mu}\) and \({\bar \kappa}\) are defined by \[\label{wdjnhvbx} {\bar \mu} = \frac{1}{(1-\nu)}\bar{p}, \quad {\bar\kappa} =\frac{d_v + 2}{2} \bar{p},\tag{58}\] with \(\bar{p}=(p^{(1)}, p^{(2)},...,p^{(s)})^\top\).

Similarly, analogous considerations can be made for type II by introducing a new definition of the matrix \(B\) and weights \(w\) (for details, see paper [2]). However, in practice, \(\tau\) is closely related to the viscosity and heat conductivity of gases, and hence the form of \(\tau\) should be set carefully (see for example [3], [8], [19], [21]).

IMEX-RK scheme of type II. In this section we analyze the asymptotic behavior of type II IMEX-RK schemes. As type I scheme, at this stage, we do not need to assume that the scheme satisfies the GSA condition.

We conclude this section with the following observations. If conditions ?? are not satisfied, order reduction is likely to occur in these intermediate regimes of \(\varepsilon\). Therefore, we adopt IMEX-RK schemes of type II that satisfy conditions ?? ensuring the validity of the estimates in Theorem 11 and preventing order reduction in the intermediate regime.
In this paper, we examine two IMEX-RK schemes of type ARS, introduced in [2] and called IMEX-II-GSA3 and IMEX-II-ISA3. These two schemes are specifically designed in [2] to improve the resolution of hyperbolic relaxation systems in the intermediate regime \(\varepsilon \approx \Delta t\).
Before presenting these two schemes, we note that from the second condition in ?? , under the assumption of \(\tilde{b} = b\) we get \[\label{bb95b} \hat{\tilde{b}}^\top \hat{A}^{-1}({a} - \tilde{a}) = \hat{{b}}^\top \hat{A}^{-1}({a} - \tilde{a}) = \hat{\boldsymbol{e}}_s^\top (a - \tilde{a})\tag{75}\] which equals zero if \(\tilde{a}_{s1} = a_{s1}\).
In Section 4, we will show that IMEX-II-GSA3 scheme, which is GSA and \(\tilde{b} \neq b\) leads to a mild order reduction for small values of \(\varepsilon\). On the contrary, IMEX-II-ISA3 with \(\tilde{b} = b\), does not exhibit order reduction, as it satisfies 75 , (as confirmed by the numerical results in the next section).

4 Numerical Results↩︎

In this section, we test the performance of different types of IMEX RK schemes using two distinct models: the 1 + 1 BGK model and the 1 + 2 ES-BGK model. We use the fifth-order finite difference WENO method [25] to approximate spatial derivatives. The discretization of the velocity domain is achieved through uniform grid points within a sufficiently large interval or domain. The choice of free parameter \(\nu\) and collision frequency \(\tau\) will be explained in each problem.

Throughout this section, we consider IMEX-RK schemes of type I and II. In particular, we take into account the GSA IMEX-II-GSA3 scheme and IMEX-II-ISA3 with \(\tilde{b} = b\). For comparison, we consider the IMEX-RK(4,3,3) scheme of type I not GSA in [26], and the scheme GSA ARS(4,4,3) of type II in [12].

In the following, these schemes are represented, as usual, by the double Butcher tableau.

• Third-order IMEX-RK(4,3,3) scheme of type I in [26]. (Top: explicit method. Bottom: implicit method).

  \[\begin{align} \begin{aligned} \begin{array}{c|cccc} 0& 0 & 0 & 0 & 0\\ \gamma & \gamma & 0 & 0 & 0\\ c_3 & 1.243893189483362 & -0.525959928729133 & 0& 0\\ 1 & 0.630412558152867 & 0.786580740199155 & -0.416993298352022 & 0\\ \hline & 0 & 1.208496649176010 & -0.644363170684468 & \gamma \end{array} \\[4mm] \nonumber \begin{array}{c|cccc} \gamma & \gamma & 0 & 0& 0\\ \gamma & 0 & \gamma & 0& 0\\ c_3 & 0 & 0.282066739245771 & \gamma & 0 \\ 1 & 0 & 1.208496649176010 & -0.644363170684468 & \gamma\\ \hline & 0 & 1.208496649176010 & -0.644363170684468 & \gamma \end{array} \end{aligned} \end{align}\]

  \(\gamma = 0.435866521508459, \, c_3 = 0.7179332607542294\).

• Third-order GSA ARS(4,4,3) in [12]:

  \[\label{A3:ARS2-IMEX3} \begin{array}{c|ccccc} 0 & 0 & 0 & 0 &0 &0 \\ 1/2& 1/2& 0 & 0 & 0 &0\\ 2/3 & 11/18 & 1/18 & 0 & 0 &0 \\ 1/2 & 5/6 &-5/6 & 1/2 & 0 &0\\ 1 & 1/4 &7/4 & 3/4 & -7/4 &0\\ \hline \\ & 1/4 &7/4 & 3/4 & -7/4 &0 \end{array} \qquad \begin{array}{c|ccccc} 0 & 0 & 0 & 0 &0 &0 \\ 1/2& 0&1/2 & 0 & 0 &0\\ 2/3 & 0 & 1/6 & 1/2 & 0 &0 \\ 1/2 & 0 &-1/2 & 1/2 & 1/2 &0\\ 1 & 0 &3/2 & -3/2 & 1/2 &1/2\\ \hline \\ &0& 3/2 &-3/2 & 1/2 & 1/2 \end{array}\tag{76}\]

• Third-order GSA IMEX-II-GSA3 of type ARS, [2]. (Top: explicit method. Bottom: implicit method): \[\begin{align} \label{GSA:CK-IMEX3} \begin{aligned} \begin{array}{c|ccccccc} 0 & 0 & 0 & 0 &0 &0 &0 &0 \\ 43/100& 43/100& 0 & 0 & 0 &0&0 &0 \\ 336/929 & 0 & 336/929 & 0 & 0 &0 &0 &0 \\ -29/42 & 0 & -29/42 & 0 & 0 &0&0 &0 \\ 581/527 & 0 &-1213/770 & 2491/956 & 267/3701 &0&0 &0 \\ 2/3 & 0 & -197/1238 & 499/743 & 0 &581/3768 &0 &0 \\ 1 & 0 & 263/620 & 134/16589 & 1040/22119 &0 &4777/9174 &0 \\ \hline \\ & 0 & 263/620 & 134/16589 & 1040/22119 &0 &4777/9174 &0 \\ \end{array}\\ \begin{array}{c|ccccccc} 0 & 0 & 0 & 0 &0 &0 &0 &0 \\ 43/100& 0 & 43/100 & 0 & 0 &0&0 &0 \\ 336/929 & 0 & -168/2459 & 43/100 & 0 &0 &0 &0 \\ -29/42 & 0 & -2353/2100 & 0 & 43/100 &0&0 &0 \\ 581/527 & 0 &889/1322 & 0 & 0 &43/100 &0 &0 \\ 2/3 & 0 & 247/2416 & 0 &408/3035 &0 & 43/100&0 \\ 1 & 0 & 872/1201 & 0 & 139/4081 & -50/237 &434/20817 &43/100 \\ \hline \\ & 0 & 872/1201 & 0 & 139/4081 & -50/237 &434/20817 &43/100 \\ \end{array} \end{aligned} \end{align}\tag{77}\]

• Third-order IMEX-II-ISA3, [2]. (Top: explicit method. Bottom: implicit method):

  \[\begin{align} \label{A3:CK-IMEX3} \begin{aligned} \begin{array}{c|ccccccc} 0 & 0 & 0 & 0 &0 &0 &0 &0 \\ 1/5& 1/5& 0 & 0 & 0 &0&0 &0 \\ 1/3 & 0 & 1/3 & 0 & 0 &0 &0 &0 \\ 2/3 & 0 & 557/867 & 7/289 & 0 &0&0 &0 \\ 3/4 & 0 &16/289 & 803/1156 & 0 &0&0 &0 \\ 1 & 0 & 13348/3993 & -9355/3993 & 0 & 0 &0 &0 \\ 1 & 0 & 75/154 & 0 & -3/14 &8/11&0 &0 \\ \hline \\ & 0 & -155/112 & 251/80 &-547/280 &2/3 &1/3& 1/5 \\ \end{array}\\ \begin{array}{c|ccccccc} 0 & 0 & 0 & 0 &0 &0 &0 &0 \\ 1/5& 0 & 1/5 & 0 & 0 &0&0 &0 \\ 1/3 & 0 & 2/15 & 1/5 & 0 &0 &0 &0 \\ 2/3 & 0 & 7/15 & 0 & 1/5 &0&0 &0 \\ 3/4 & 0 &1137/1004 & -731/1255 & 0 &1/5 &0 &0 \\ 1 & 0 & 447/565 & 0 & -636/613 &519/496 & 1/5&0 \\ 1 & 0 & -155/112 & 251/80 &-547/280 &2/3 &1/3& 1/5 \\ \hline \\ & 0 & -155/112 & 251/80 &-547/280 &2/3 &1/3& 1/5 \\ \end{array} \end{aligned} \end{align}\tag{78}\]

4.1 Accuracy test.↩︎

Test 1. In this test we investigate numerically the convergence rate solving the BGK model for a smooth solution. The initial data is taken as \[\begin{align} \label{init32well-prepared321D} \begin{aligned} \rho_0(x)&=1+0.2\sin(\pi x),\quad u_0(x)=\left(1,0\right),\quad T_0(x)=\frac{1}{\rho_0(x)}\cr f_0(x,v)&=\mathcal{M}(x,v) - \frac{\varepsilon}{\tau}\left(I-\Pi_\mathcal{M}\right)(v_1\partial_x \mathcal{M}) \end{aligned} \end{align}\tag{79}\] where \[\mathcal{M}(x,v)=\frac{\rho_0(x)}{2\pi T_0(x)}\exp\left(-\frac{|v-u_0(x)|^2}{2T_0(x)}\right)\] For an explicit definition of the term \(\left(I-\Pi_\mathcal{M}\right)(v_1\partial_x \mathcal{M})\), please refer to Remark 14 in the appendix.
We use \(32\) uniform points for velocity in the interval \(v \in [-10, 10]\), and \(N_x\) uniform points for space in the interval \(x \in [0, 2]\) with periodic boundary conditions. We set \(\Delta t = \text{CFL}\Delta x /10\) with \(\text{CFL} = 0.9\) and \(\tau=1\).

Since the exact solution is not available for this test, as well as for the subsequent test, the \(L^1\) errors are computed as the difference of the numerical solutions on two consecutive meshes in space, i.e., we use the numerical solution on a finer mesh with mesh size \(\Delta x/2\) as the reference solution to compute the error for solutions on the mesh size of \(\Delta x\).

Table [1d32accuracy] shows results for third order accurate IMEX-RK schemes at \(t = 0.25\). First we can see order reduction in the intermediate regime for ARS(4,4,3) and IMEX-RK(4,3,3) scheme. Although we observe order reductions for IMEX-RK at small Knudsen numbers, such as \(10^{-4}\), it is a typical issue of IMEX-RK methods. In general, it is highly nontrivial to avoid order reduction of IMEX-RK schemes in the intermediate regime (for example, see papers [15]–[17]).

For comparison, we test the two third-order schemes introduced in [2] IMEX-II-GSA3 (77 ) and IMEX-II-ISA3 (78 ) for this example. The third order accuracy is achieved for IMEX-II-ISA3 for all values of \(\varepsilon\) even in the intermediate regime, i.e., the scheme maintains uniform accuracy, on the other hand IMEX-II-GSA3 shows a slight order reduction for the value \(\varepsilon = 10^{-4}\). Note that, as discussed in Section 2, since the IMEX-II-ISA3 scheme has same weights, it performs better than the other methods for intermediate values of \(\varepsilon\).

Test 2. In this test we investigate numerically the convergence rate solving the ES-BGK model for a smooth solution considering \(d_v = 2\) with \(v = (v_1, v_2)\), and the initial well-prepared initial data (79 ), where \[\label{init32well-prepared322D} \mathcal{M}(x,v)=\frac{\rho_0(x)}{2\pi T_0(x)}\exp\left(-\frac{(v_1-u_0(x))^2 + v_2^2}{2T_0(x)}\right)\tag{80}\] For the velocity discretization, we use \(32 \times 32\) uniform points in the domain \(v \in [-10, 10] \times [-10,10]\). Note that, for the IMEX-RK(4,3,3) scheme in order to achieve the expected order of convergence, we required more points in velocity. In this test we use \(N_v = 80\) uniform points for the velocity. We use \(N_x\) uniform points for space in the interval \(x \in [0,2]\) with periodic boundary conditions. The time step is taken as \(\Delta t = \text{CFL} \Delta x/10\) with \(\text{CFL} = 0.9\), \(\tau=1\) and \(\nu = -1/2\).

Table [2d95bis32accuracy] shows the result for the third order accurate IMEX-RK schemes at \(t = 0.25\). As the one dimensional case we can see order reduction in the intermediate regime (\(\varepsilon=10^{-4}\)) for the ARS(4,4,3) scheme. The third order scheme IMEX-II-ISA3 achieves the third order accuracy for all the values of \(\varepsilon\) even in the intermediate regime, whereas IMEX-II-GSA3 and IMEX-RK(4,3,3) schemes show a slight order reduction for the value \(\varepsilon = 10^{-4}\).

Test 3. Riemann Problem. In this problem, we consider a Riemann problem in 1D space and 2D velocity domain. The same test has been adopted in [21] to show the consistency between the ES-BGK model and BTE or NSE. As initial macroscopic states, we use \[\left(\rho_0, u_{x0}, u_{y0}, T_0 \right) = \Biggl\{ \begin{array}{cc} (1, M\sqrt{2}, 0,1), & -1 \le x \le 0.5,\\[2mm] (\frac{1}{8},0,0,\frac{1}{4}), & \textrm{otherwise}, \end{array}\] where \(M = 2.5\) is the Mach number. Here we impose the free-flow boundary condition in space \(x \in [-1, 2]\). We truncate the velocity domain with \(v_{max} = 15\). We compare numerical solutions at time \(t=0.15,\,0.2,\,0.4\) with the time step that corresponds to \(CFL = 0.5\). Here we set \(\nu = -1\) and \(\tau = 0.9\pi \rho/2\) following [21], which gives the following viscosity and heat conductivity: \[\mu = \frac{1}{0.9 \pi}T, \quad \kappa = \frac{1}{0.9\pi}T.\] We consider these transport coefficients when computing the reference solutons to NSE. We remark that this choice matches the viscosity of ES-BGK model and the one derived from BTE for Maxwellian molecules [21].

In Figs. 1-2, we present numerical results for BTE, NSE and ES-BGK model for different Knudsen numbers \(\varepsilon = 0.5, 0.1\) respectively. For BTE, we use an explicit fourth-order semi-Lagrangian scheme [26], while for NSE, we use the classical RK4 and WENO23. For ES-BGK, we use the combination of IMEX-II-ISA3 for the time discrsetization and WENO23 for the space one. For \(\varepsilon=0.5\), the numerical solution of the ES-BGK model is very close to that of BTE, while solutions to NSE are deviated from the other solutions. On the other hand, for \(\varepsilon=0.1\), the discrepancy between three solutions become smaller, which confirms the consistency of the three models.

Test 4. The Lax Shock Tube Problem. Finally we test IMEX-RK schemes solving the ES-BGK model for the Lax shock tube problem, [1]. We set the collision frequency \(\tau\) for ES-BGK model as \[\tau = \frac{2}{3}\rho \sqrt{T},\] which yields the viscosity and heat conductivity: \[\mu = \sqrt{T}, \quad \kappa = \frac{15}{4}\sqrt{T},\] where \(\nu = -1/2\), and this implies that the Prandtl number is \(2/3\). We consider the initial macroscopic variables \[\left(\begin{array}{c} \rho \\ u \\ p \end{array} \right) = \Biggl\{ \begin{array}{cc} (0.445, 0.698, 3.528)^\top, & -0.5 \le x \le 0,\\[2mm] (0.5,0,0.571)^\top, & 0 < x \le 0.5, \end{array}\] on the free-flow condition \(x \in [-5,5]\). We take well-prepared initial conditions \[f_0(x,v) =\mathcal{M}(x,v) - \frac{\varepsilon}{\tau}\left(I-\Pi_\mathcal{M}\right)(v_1\partial_x \mathcal{M}),\] where \(I\) is the identity operator and \(\Pi_\mathcal{M}\) is the projection operator defined as ?? . We first use \(80 \times 80\) uniform points for the velocity domain \([-20, 20] \times [-20, 20]\), and \(N_x = 200\) uniform points for the spatial discretization. The CFL number is taken as \(0.2\) and final time \(1.3\). We consider IMEX-II-ISA3 scheme coupled with WENO35 for solving the ES-BGK model.

In Fig. 3, we compare the numerical solution of the ES-BGK model with those of the compressible Navier-Stokes equations and the compressible Euler equations. The reference solution of the Euler equations and of the Navier-Stokes equations was generated by using spectral method for the spacial derivatives and fourth-order RK method for the time on a grid of 3200 for \(\varepsilon = 0, 10^{-2}, 10^{-4}\) respectively. As we can see from the Fig. 3 the numerical solution of the ES-BGK model is very close to that of the NS equations. Furthermore for a better visualization, we zoomed the two solutions on a portion of the domain.

In Fig. 4 we reported the moments \[\begin{align} \label{195A} \begin{aligned} -\mu \sigma(u) &= \int_{\mathbb{R}^{d_v}} f_1 (v-u) \otimes (v-u)dv = \left\langle f_1(v -u)\otimes(v-u)\right\rangle,\\ -\kappa \nabla T &= \int_{\mathbb{R}^{d_v}} f_1 \frac{1}{2}(v-u) |v-u|^2 dv = \left\langle f_1(v -u)|v-u|^2\right\rangle, \end{aligned} \end{align}\tag{81}\] i.e., the shear stress and heat flux, computed (1) from the function \(f_1 = \frac{\mathcal{G}-f}{\varepsilon}\) (2) from the macroscopic quantities and 3) from the NS equation for the case of \(\varepsilon = 10^{-2}\) and \(\varepsilon = 10^{-4}\).

5 Conclusions↩︎

In this work, we investigated the asymptotic behavior at the Navier-Stokes level, reproducing theoretical results, Theorems 9-11. We showed that existing IMEX-RK schemes including type I and II can accurately capture the NS limit without resolving the small scales determined by the Knudsen number. To address the issues on order reduction of these schemes at the Navier-Stokes level, we propose IMEX-RK schemes of type II developed in [2]. In particular, the IMEX-II-ISA3 scheme of type ARS satisfying additional conditions ?? guarantees uniform accuracy avoiding order reduction. Finally, we provided numerical examples that validate the theoretical results.

Acknowledgements↩︎

Sebastiano Boscarino is supported for this work by 1) the Spoke 1 "FutureHPC \(\&\) BigData" of the Italian Research Center on High-Performance Computing, Big Data and Quantum Computing (ICSC) funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di "campioni nazionali di R\(\&\)S (M4C2-19 )", (CN00000013); by 2) the Italian Ministry of Instruction, University and Research (MIUR) to support this research with funds coming from PRIN Project 2022 (2022KA3JBA), entitled "Advanced numerical methods for time dependent parametric partial differential equations and applications"; 3) from Italian Ministerial grant PRIN 2022 PNRR "FIN4GEO: Forward and Inverse Numerical Modeling of hydrothermal systems in volcanic regions with application to geothermal energy exploitation.", (No. P2022BNB97). S. Boscarino is a member of the INdAM Research group GNCS. S. Y. Cho was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. RS-2022-00166144), and Learning & Academic research institution for Master’s\(\cdot\)PhD students, and Postdocs (LAMP) Program of the National Research Foundation of Korea (NRF) grant funded by the Ministry of Education (No. RS-2023-00301974).

6 Appendix↩︎

6.1 Basic properties of Gaussian distribution function↩︎

By definition, \(\mathcal{G}[f]\) satisfies \[\begin{align} \label{G9595U} \begin{aligned} &\biggl\langle f \biggr\rangle = \int_{\mathbb{R}^{d_v}} f(v) dv = \int_{\mathbb{R}^{d_v}} \mathcal{G}[f](v) dv = \biggl\langle \mathcal{G}[f] \biggr\rangle = \rho,\\ &\biggl\langle v f \biggr\rangle = \int_{\mathbb{R}^{d_v}} vf(v) dv = \int_{\mathbb{R}^{d_v}} v\mathcal{G}[f](v) dv = \biggl\langle v \mathcal{G}[f](v)\biggr\rangle = \rho u,\\ &\biggl\langle \frac{|v|^2}{2} f\biggr\rangle = \int_{\mathbb{R}^{d_v}}\frac{|v^2|}{2} f(v) dv = \int_{\mathbb{R}^{d_v}} \frac{|v^2|}{2}\mathcal{G}[f](v) dv= \biggl\langle \frac{|v|^2}{2} \mathcal{G}[f](v)\biggr\rangle = E. \end{aligned} \end{align}\tag{82}\] Note that from \[\begin{align} \label{G9595Um} \begin{aligned} &\int_{\mathbb{R}^{d_v}} (v-u) \otimes (v-u) f(v) dv = \rho \Theta,\\ & \int_{\mathbb{R}^{d_v}} (v-u) \otimes (v-u) \mathcal{G}[f](v) dv = \rho \mathcal{T}, \end{aligned} \end{align}\tag{83}\] we further have \[\begin{align} \label{Prop1} \begin{aligned} f &= \mathcal{G}[f] \Longleftrightarrow f = \mathcal{M}[f];\\[3mm] f &= \mathcal{G}[f] + \mathcal{O}(\varepsilon) \quad \textrm{implies} \quad \mathcal{G}[f] = \mathcal{M}[f] + \mathcal{O}(\varepsilon). \end{aligned} \end{align}\tag{84}\] The proof is straightforward. In fact, on both side of \(f = \mathcal{G}[f]\) using (83 ) we get \(\Theta = \mathcal{T}\) and from (8 ) we get \(\mathcal{T} = TId\) hence \(\mathcal{G}[f]\) is just the isotropic Maxwellian \(\mathcal{M}[f]\). Similarly if on both side of \(f = \mathcal{M}[f]\) takes (83 ). For the second point in (84 ) the proof is similar and we omit the details.

6.2 Chapmann-Enskog expansion of ES-BGK model↩︎

Here we perform the first order Chapman-Enskog expansion of ES-BGK model with respect to \(\varepsilon\): \[\label{firstEx} f = \mathcal{M}(f) + \varepsilon f_1,\tag{85}\] which implies that \(f_1\) satisfies the so-called compatibility relations: \[\label{Comp} \langle \phi f_1 \rangle=0, \quad \phi = (1,v,|v|^2).\tag{86}\] Now we take the expansion in terms of \(\varepsilon\) for the stress tensor \[\label{ST} \Theta = T Id + \varepsilon \Theta_1,\tag{87}\] and the heat flux \[\label{HF} \mathbb{Q} := \left\langle \frac{|v-u|^2}{2}(v-u)f \right\rangle = 0 + \varepsilon \mathbb{Q}_1\tag{88}\] with \[\begin{align} \label{lrzitkae} \Theta_1 = \frac{1}{\rho }\int_{\mathbb{R}^{d_v}} f_1 (v-u) \otimes (v-u)dv = \frac{1}{\rho } \left\langle f_1(v -u)\otimes(v-u)\right\rangle, \end{align}\tag{89}\] and \[\label{395A} \mathbb{Q}_1 = \int_{\mathbb{R}^{d_v}} f_1 \frac{1}{2}(v-u) |v-u|^2 dv = \left\langle f_1\frac{1}{2}(v -u)|v-u|^2\right\rangle.\tag{90}\] Inserting (85 ) and these latter expansions for \(\Theta\) and \(\mathbb{Q}\) into the conservation laws \[\langle \phi f \rangle + \langle \phi v \cdot \nabla_x f \rangle = 0,\] or \[\langle \phi \mathcal{M}(f) \rangle + \langle \phi v \cdot \nabla_x \mathcal{M}(f) \rangle = - \varepsilon \nabla_x\langle \phi v f_1 \rangle,\] it gives the equations: \[\label{DiscCL95bis} \partial_t U + \nabla_x F(U) = -\varepsilon \nabla_x H({\boldsymbol{U}})\tag{91}\] with \[U = \langle \phi f \rangle = \left( \begin{array}{c} \rho\\ \rho u\\ E \end{array} \right), \quad F(U) = \left( \begin{array}{c} \rho u\\ \rho u \otimes u + \rho T Id \\ (E + \rho T)u \end{array} \right), \quad H({\boldsymbol{U}}) = \left( \begin{array}{c} 0\\ \rho \Theta_1\\ \mathbb{Q}_1 + \rho \Theta_1 u \end{array} \right).\] To evaluate the quantities \(\rho \Theta_1\) and \(\mathbb{Q}_1\), we seek the form of \(f_1\). For this aim, we consider the expansion of the anisotropic Gaussian \(\mathcal{G}(f)\) with respect to \(\varepsilon\), i.e. \[\label{G95ex} \mathcal{G}(f) = \mathcal{M}(f) + \varepsilon g.\tag{92}\] By (8 ) and considering the expansion (87 ), we get \(\mathcal{T} = T Id + \nu \varepsilon \Theta_1\), and using the fact that \(tr(\Theta_1) = 0\), this gives \(det(\mathcal{T}) = T^{d_v} + \mathcal{O}(\varepsilon^2)\) and \[[\mathcal{T}]^{-1} = \frac{1}{T}\left(Id - \frac{\nu \varepsilon}{T}\Theta_1\right) + \mathcal{O}(\varepsilon^2),\] and therefore by (7 ) we obtain \[\label{gVect} g := \frac{\mathcal{G}(f) - \mathcal{M}(f)}{\varepsilon} = \nu \frac{ \mathcal{M}(f)}{2 T^2}(v-u)^\top\Theta_1(v-u).\tag{93}\] where we neglect the term \(\mathcal{O}(\varepsilon^2)\). Now the next step is to insert expansions (85 ) and (92 ) into the equation \[\partial_t f + v \cdot \nabla_x f = \frac{\tau}{\varepsilon}(\mathcal{G}[f]-f).\] Then, we get \[\label{1} \tau (f_1 - g) = - \left( \partial_t \mathcal{M} + v\cdot \nabla_x \mathcal{M} +\varepsilon\left(\partial_t {f}_1 + v\cdot \nabla_x {f}_1 \right)\right).\tag{94}\]

Before proceeding, we review a Hilbert space and an associated operator considered in [19] and some Lemmas.

By Remark 14, we can rewrite (94 ) as \[\label{f11} f_1= g-\frac{1}{\tau}\left(\mathcal{M}\left(A(V):\frac{\sigma(u)}{2} + 2B(V)\cdot \nabla_x \sqrt{T}\right)\right) + \mathcal{O}(\varepsilon).\tag{95}\]

Next, we multiply both sides of (95 ) by \(v \phi\) and take integration and use Lemmas 15-17 to obtain \[\begin{align} \label{f1195bis} \begin{aligned} \left\langle v\phi f_1\right\rangle&= \left\langle v\phi g\right\rangle- \frac{1}{\tau}\left\langle v\phi \mathcal{M}\left(A(V):\sigma(u) + 2B(V)\cdot \nabla_x \sqrt{T}\right)\right\rangle + \mathcal{O}(\varepsilon),\\ &H(U) = \mathcal{G}(U) - \frac{1}{\tau}\mathcal{F}(U) + \mathcal{O}(\varepsilon), \end{aligned} \end{align}\tag{96}\] that is, \[\label{f1195tris} \left( \begin{array}{c} 0\\ \rho \Theta_1\\ \mathbb{Q}_1 + \rho \Theta_1 u \end{array} \right)= \left( \begin{array}{c} 0\\ \nu \rho \Theta_1\\ \nu \rho \Theta_1 u \end{array} \right) - \frac{1}{\tau} \left( \begin{array}{c} 0\\ \rho T\sigma(u)\\ \rho T\sigma(u)u + \frac{d_v +2}{2}\rho T \nabla_x T \end{array} \right) + \mathcal{O}(\varepsilon).\tag{97}\] Then, we get \[\rho \Theta_1 =- \mu\sigma(u) + \mathcal{O}(\varepsilon),\] with viscosity \(\mu = \frac{\rho T}{(1-\nu)\tau}\) and \[\mathbb{Q}_1 =- \kappa \nabla_x T + \mathcal{O}(\varepsilon),\] with the thermal conductivity \[\kappa = \frac{d_v + 2}{2}\frac{p}{\tau}\] where \(p = \rho T\). Finally, from (91 ) we derive the CNS equations \[\label{Mom295bis} \partial_t U + \nabla_x \cdot F(U) = \varepsilon \nabla_x \cdot S(U),\tag{98}\] with \[S(U)=\left( \begin{array}{c} 0\\ \mu \sigma(u)\\ \mu\sigma(u)u -\kappa \nabla_x T \end{array} \right),\] \[\sigma(u) = \nabla_x u + (\nabla_x u)^ \top - \frac{2}{d}\nabla_x \cdot u I.\]

References↩︎