Further acceleration of multiscale simulation of rarefied gas flow via a generalized boundary treatment

To simultaneously account for real gas effects while avoiding complex kinetic models, we consider the kinetic model of polyatomic gas with rotational degrees of freedom $d_{r}$, excluding vibrational and electronic degrees of freedom. The gas system can be described by the distribution function $f^{T}(t,\bm{x},\bm{\xi},I_{r})$, where the variable $\mathbf{x}$ denotes spatial position, $t$ represents time, $\bm{\xi}$ represents the molecular velocity, and $I_{r}$ is the rotational energy. Two velocity distribution functions (VDFs), $f_{0}(\mathbf{x},\bm{\xi},t)=\int f^{T}dI_{r}$ and $f_{1}(\mathbf{x},\bm{\xi},t)=\int I_{r}f^{T}dI_{r}$, can be introduced to reduce the rotational energy variable, and the final dimensionless kinetic model in the absence of external force can be expressed as [38, 39]

The function $g$ denotes the reference distribution function, where the subscripts $t$ and $r$ stand for the relaxation processes associated with translational and rotational motions, respectively:

with the equilibrium distribution functions

It should be noted that, from Eq. (1) we known that the VDFs $f_{0}$ and $f_{1}$ evolve towards the reference VDFs in Eq. (4). Meanwhile, during the evolution of $f_{0}$ and $f_{1}$, the heat fluxes in the reference VDFs gradually go to zero in isolated system. Therefore, in the final steady state, the VDFs $f_{0}$ and $f_{1}$ will be described by the corresponding equilibrium distribution functions in Eq. (5), which is consistent with the Boltzmann H-theorem.

Macroscopic quantities including the density $\rho$, velocity $\mathbf{u}$, translational temperature $T_{t}$, rotational temperature $T_{r}$, as well as the stress $\bm{\sigma}$, translational heat flux $\mathbf{q}_{t}$ and rotational heat flux $\mathbf{q}_{r}$ are expressed through the moment integration of the distribution functions as

Furthermore, relationships pertaining to the total temperature $T$, translational pressure $p_{t}$, and total pressure $p$ are delineated as follows:

It is noted that the derivations mentioned here are based on dimensionless kinetic equations. The dimensionless treatment of macroscopic variables and distribution functions is presented in A for reference.

We briefly introduce the CIS to solve the kinetic equation by the implicit finite-volume method. Considering that the evolution equations for translational and rotational distribution functions in Eq. (1) exhibit identical structural formulations, here we employ the notations $f$ and $g$ to represent a unified form of translational/rotational distribution functions. The discrete form of the mesoscopic kinetic equation within the control volume cell and a time step $\Delta t_{f}=t^{n+1}-t^{n}$ from the $n-$th step to the $(n+1)-$th step is shown as follows

where the subscript $i$ denotes the average variables at the cell center, while the subscript $ij$ represents the cell interface located between cells $i$ and $j$. $S_{ij}$ and $\Omega_{i}$ respectively denote the area of the cell interface and the volume of the cell. ${\xi}_{n}=\bm{{\xi}}_{n}\cdot\mathbf{n}_{ij}$ represents the molecular velocity component along the interface normal vector $\mathbf{n}_{ij}$.

To facilitate matrix-free implicit iterations, an incremental form of the distribution function is introduced:

which allows the rewriting of Eq. (9) as

The CIS is found to be efficient when the Knudsen number is large, but is highly dissipative and inefficient when the Knudsen number is small [14]. To expedite the convergence of the kinetic equations, the VDF is firstly evaluated at the intermediate step, i.e., $n+1$ in Eqs. (9) and (10) is replaced by $n+1/2$. Subsequently, macroscopic synthetic equations are introduced in GSIS which will guide the evolution of the VDF towards steady-state solutions.

The proper design of macroscopic synthetic equation is crucial in GSIS, which the guidelines are: (i) the synthetic equation should be derived exactly from the kinetic equation, (ii) the synthetic equation should facilitate efficient exchange of flow information in the whole computational domain, and (iii) the synthetic equation should be as simple as possible, so that it is compatible with the numerical techniques in the computational fluid dynamics. In this case, the second version of GSIS (which contains the evolution equations for the density, velocity, and temperature only) is used [27, 24], rather than the first version (which not only contains the evolution equations for the density, velocity, and temperature, but also the evolution equations for the stress and heat flux) [20].

Exploiting the mesoscopic-macroscopic relationships in Eq. (6) and applying them to Eq. (1), the governing equations for macroscopic quantities can be derived as

When the Knudsen number is small, the stress and heat flux in Eq. (13) can be approximated via the first-order Chapman-Enskog expansion of the kinetic equation as [41]:

where $\kappa_{t}$ and $\kappa_{r}$ are the translational and rotational thermal conductivities:

The macroscopic synthetic equation is constructed as follows:

where

It should be emphasized that, in the numerical simulation, the variables at the left-hand-side of Eq. (16) are evaluated at the $(n+1)$-th time step, while the one at the right-hand-side are evaluated at the $(n+1/2)$-th time step. Only when the solution finally converges (the solution at the $(n+1/2)$-th time step is the same as that at the $(n+1)$-th time step) can the synthetic equation turns back to Eq. (13) which is exactly derived from the kinetic equation. In the intermediate iteration step, Eq. (16) is not the same as the exact one (13), but such treatment plays an important role in boosting the convergence to the steady state.

Compared to the original macroscopic equation (13), the synthetic equation introduces high-order corrections in the stress and heat flux $\Delta\bm{\sigma},\Delta\mathbf{q}_{t},\Delta\mathbf{q}_{r}$ on the right-hand-side of equations, ensuring the accuracy of macroscopic equation iterations in the presence of strong rarefaction effects. On the other hand, the explicit inclusion of the Newton law of viscosity and the Fourier law of thermal conductivity in the left-hand-side of Eq. (16) could lead to fast exchange of the flow field, since in the steady state the diffusion-type equations for the velocity and temperature can be obtained, where the information propagation speed is infinite.

Under the finite-volume discretization, the macroscopic synthetic equation within the control volume cell $i$ and a inner time step $\Delta t_{W}=t^{m+1}-t^{m}$ can be formulated as follows:

Similarly, an incremental form of macroscopic quantities is introduced as $\Delta\bm{W}_{i}^{m}=\bm{W}_{i}^{m+1}-\bm{W}_{i}^{m}$, and Eq. (18) can be reformulated as follows:

It is worth noting that the high-order corrections in Eq. (17) are computed based on the VDF and macroscopic quantities at the time step $n+1/2$, and they remain constant during inner iterations. Therefore, the high-order residual terms $\mathbf{R}_{i}^{\mathrm{HoT}}$ in Eq. (21) should be treated as source terms in the computation.

The Rusanov scheme and the central-difference scheme are employed to compute the convective flux $\mathbf{F}_{c,ij}$ and viscous flux $\mathbf{F}_{v,ij}$ in the residual $\mathbf{R}_{i}^{m}$ of Eq. (21), respectively. Meanwhile, the incremental form of the macroscopic flux $\Delta\mathbf{F}_{ij}$ can be approximated by the Eulerian form of the convective flux. Equations (11) and (21) can be iteratively solved using the classical Lower Upper Symmetric Gauss–Seidel (LU-SGS) technology. Given the widespread application of LU-SGS, specific details are omitted here for the sake of brevity. Relevant literature specifying the settings for implicitly solving these equations using LU-SGS are provided in the reference [27, 24].

The framework of the GSIS algorithm, aimed at accelerating convergence by alternately solving mesoscopic kinetic equation and macroscopic synthetic equation within one iteration step, is summarized below:

• 1)

  Based on the macroscopic quantities $\mathbf{W}_{i}^{n}$ and mesoscopic CFL number $\vartheta_{f}$, the time step $\Delta t_{f}$ for the mesoscopic kinetic equations can be set according to Eq. (12).

• 2)

  The implicit solution of the mesoscopic kinetic equations evolves the VDF to $f_{i}^{n+1/2}$ according to Eq. (11). The kinetic gas-surface boundary condition is used.

• 3)

  The corresponding macroscopic quantities $\mathbf{W}_{i}^{n+1/2}$, as well as the stress $\bm{\sigma}^{n+1/2}$ and the heat flux $\mathbf{q}_{t}^{n+1/2},\mathbf{q}_{r}^{n+1/2}$, are obtained through moment integration of the distribution function as presented in Eq. (6).

• 4)

  Calculate the derivatives of the macroscopic quantities by the least square method. Then, the high-order correction including $\Delta\bm{\sigma}^{n+1/2},\Delta\mathbf{q}^{n+1/2},\Delta\mathbf{q}^{n+1/2}$ can be obtained from Eq. (17).

• 5)

  The LU-SGS technology is employed for the implicit iteration of the macroscopic synthetic equations (21). In contrast to the kinetic equations, which are computed only once in each step, the synthetic equations involve hundreds of inner iterations to expedite the flow information exchange. Note that the traditional no-velocity-slip and no-temperature boundary conditions do not work any more in rarefied gas flows [29], therefore the boundary condition should be carefully designed; this will be elaborated in section 3.

• 6)

  Based on the macroscopic quantities $\mathbf{W}_{i}^{k}$ from the inner iterations ($k$ is the total steps of inner iterations), the next-step VDF $f_{i}^{n+1}$ are computed as

  		$\displaystyle f_{0,i}^{n+1}=f_{0,i}^{n+1/2}+\left[f_{0}^{\mathrm{eq}}\left(\mathbf{W}_{i}^{k}\right)-f_{0}^{\mathrm{eq}}\left(\mathbf{W}_{i}^{n+1/2}\right)\right],$		(22)
  		$\displaystyle f_{1,i}^{n+1}=f_{1,i}^{n+1/2}+\left[f_{1}^{\mathrm{eq}}\left(\mathbf{W}_{i}^{k}\right)-f_{1}^{\mathrm{eq}}\left(\mathbf{W}_{i}^{n+1/2}\right)\right].$

  Then the macroscopic quantities $\mathbf{W}_{i}^{n+1}$ are consistent with those obtained through moment integration of $f_{i}^{n+1}$ in Eq. (6).

• 7)

  Repeat the steps (1-7) until the following convergence criterion is satisfied:

  $$\left.\frac{\sqrt{\sum_{i}\left(\psi_{i}^{n+1}-\psi_{i}^{n}\right)^{2}\Omega_{i}}}{\sqrt{\sum_{i}\left(\psi_{i}^{n}\right)^{2}\Omega_{i}}}\right|_{\max}<10^{-6},\quad\text{ for }\quad\psi\in\left\{\rho,\mathbf{u},T,T_{r}\right\}.$$

  (23)

It is noteworthy that, in the simulation of hypersonic and complex three-dimensional flows, it is a common practice to perform thousands steps of macroscopic iterations by first-order Euler equations to initialize the flow field [42, 43] and 10 steps of CIS iterations to initialize the VDF [44] prior to GSIS iterations.

The fast convergence of GSIS lies in the use of macroscopic synthetic equation (16) to guide the evolution of gas systems towards the final steady state [20, 26, 22]. Since the deterministic mesoscopic solver is much more time-consuming to solve than the macroscopic synthetic equation, in the outer iteration, Eq. (11) is only solved once at each iteration, while in the inner iteration, Eq. (16) is solved multiple times or till convergence. Since the macroscopic synthetic equation can be solved quickly and for many steps, the flow field has the time to exchange information globally, and hence faster convergence to the steady state is achieved in GSIS than solving the kinetic equation solely.

However, the macroscopic synthetic equation needs boundary conditions, but in general rarefied gas flows the no-velocity-slip and no-temperature-jump boundary conditions are not valid anymore. But, since the kinetic solver is already slow to converge when the Knudsen number is small, directly feeding its boundary information into macroscopic synthetic equation is not very efficient [27]. In order to further accelerate the GSIS, we need to figure out a way to provide the boundary condition to the macroscopic synthetic equation at each inner iteration [24], and this boundary condition should be compatible with that in the kinetic equation.

In the moment methods [35, 36, 37], the boundary conditions are introduced in three steps. First, the VDF is re-constructed from conserved macroscopic quantities and high-order moments in the moment equation. Second, the reflected VDF from the wall is calculated by applying the kinetic gas-surface boundary condition. Third, the macroscopic boundary conditions are explicitly derived for the wall VDF. The similar idea can be borrowed to GSIS. To be specific, the truncated distribution function applicable to polyatomic gases with rotational degrees of freedom is derived to approximate the unknown distribution function at the wall for an explicit integration. Then, the numerical flux can be directly reconstructed recursively based on the conserved macroscopic quantities, as well as the stress and heat flux at the interface. The stress and heat flux embedded in the truncated distribution function can be closed based on explicit constitutive relations of the Newton and Fourier. Additionally, the concept of high-order correction to the constitutive equations, which is employed in GSIS, is also incorporated into this macroscopic boundary condition, ensuring stability and accuracy in hypersonic flows and strong rarefied flows.

Within the finite-volume framework, the boundary condition can be treated as an special cell interface. Therefore, the implementation of boundary conditions requires the VDF and macroscopic numerical flux at the interface based on the gas-surface kinetic boundary condition. Taking the classical Maxwell boundary condition as an example, assuming the left side of the interface $ij$ is the wall and the right side is the interior cell $j$, the distribution function at the interface $ij$ can be determined as follows

In the first GSIS for nonlinear flows [27], the numerical fluxes of the macroscopic equations at the interface $\mathbf{F}_{ij}$ are calculated by the VDFs at the interface $ij$ as

While this boundary treatment ensures that the macroscopic equations attain consistent boundary conditions with the mesoscopic kinetic equations, it fails to maintain updates during macroscopic iterations, which does not realize the best possible performance of GSIS. In the second GSIS for nonlinear flows [24], we introduce a correction for macroscopic boundaries based on the increment of macroscopic quantities at the cell center. That is, in the inner iteration, after a single LU-SGS iteration, we obtain the global macro quantity increment and use a method similar to the Roe scheme to modify the interface flux using the physical quantity increment of the boundary element:

This approach involves updating the boundary flux based on the updates in the interior flow, ensuring that the macroscopic wall flux evolves in tandem with the interior field to boost convergence. Nevertheless, the macroscopic boundary, unable to evolve independently, would result in a delay in the evolution of the macroscopic equations at the boundary during the inner iterations. This delay hinders the acceleration of convergence for the kinetic equations at the boundary.

Hence, with the hinder-sight, we propose a GBT for the macroscopic synthetic equation. The fundamental concept lies in approximating the incident distribution function at the interface using the truncated distribution function, thereby achieving an explicit reconstruction of the macroscopic numerical flux at the interface. In order to approximate the distribution function and close the higher-order terms in the macroscopic equations, Grad proposed a series of distribution functions based on Hermite orthogonal polynomials. Following this expansion, we present the truncated distribution function $f_{ij}^{T}(t,\bm{x}_{ij},\bm{\xi},I_{r})$ applicable to polyatomic gas with rotational energy $I_{r}$ as

With the approximation of the incident distribution function $f_{k,ij}^{R}$ by the truncated distribution function $f_{ij}^{T}$, the macroscopic numerical flux at the interface can be explicitly reconstructed. Taking the case of complete diffuse reflection ($\alpha_{M}$ = 1) in Eq. (24) as an example, the macroscopic numerical flux at the interface can be expressed as follows