1 Introduction↩︎

Rarefied gas dynamics constitutes a pivotal branch of fluid mechanics, essential for understanding multiscale flow phenomena in high-altitude aerothermodynamics, micro-electromechanical systems (MEMS), and specialized gas transport processes. The flow regimes are characterized by the Knudsen number (\({\rm Kn}\)), defined as the ratio of the molecular mean free path to a characteristic macroscopic length scale. Although the Euler or Navier-Stokes equations suffice for the continuum regime (\({\rm Kn}\ll 1\)), they become physically inadequate as the rarefaction effects increase. Consequently, the Boltzmann equation, which governs the evolution of the single-particle probability distribution function in phase space, serves as the fundamental framework that is valid across all flow regimes. However, the high dimensionality of phase space and the complex quadratic collision term pose substantial challenges for efficient numerical simulations.

Traditional numerical methods for solving the Boltzmann equation are generally categorized into stochastic methods and deterministic methods. Direct Monte Carlo simulation (DSMC) [1]–[3] is one of the most popular stochastic methods that is effective in simulating high-speed rarefied gas flows. However, it suffers from computational inefficiency and inherent statistical noise in low-speed regimes. The Unified Gas-Kinetic Wave-Particle (UGKWP) method couples deterministic waves for continuum flows with stochastic particles for kinetic non-equilibrium, covering the full Knudsen number regime. Unlike DSMC, which is costly and time-step limited in near-continuum flows, UGKWP adaptively weights waves and particles by local Knudsen number, recovering a fluid solver in continuum and a particle method in rarefied regimes. This hybrid reduces noise and cost while preserving kinetic accuracy for multiscale gas dynamics [4].

Conversely, deterministic methods are often preferred for continuous or low-speed flows where the statistical fluctuations must be minimized. The discrete velocity method (DVM) [5]–[7] is a classical deterministic method, which discretizes the distribution function at a series of microscopic velocity points. The unified gas-kinetic scheme (UGKS) can eliminate restrictions on cell size and time step by simultaneously handling free flow and collision of gas molecules [4], [8]. However, since information exchange depends on the evolution of the velocity distribution function, UGKS still requires a large number of iterations to obtain the steady-state solution of the near-continuum flow [9], [10]. Furthermore, spectral methods, including the Fourier spectral method [11]–[13], the Hermite spectral method [14], [15], the Burnett spectral method [16], [17], and the mapped Chebyshev spectral method [18], [19] utilize global basis functions to approximate the distribution function, leading to a higher order of convergence. Alternatively, the moment method [20] originally proposed by Grad simplifies the kinetic description by approximating the distribution function through a truncated series of macroscopic moments.

For the steady-state problems, which are of great engineering interest, the efficiency of the numerical solver is critical. The conventional source iteration (SI) is efficient in highly rarefied regimes, but exhibits prohibitively slow convergence as the flow approaches the continuum limit (\({\rm Kn}\to 0\)) [21]. In this diffusive limit, the “converged" solution is often contaminated by numerical errors arising from velocity discretization and collision operator approximations. To overcome the slow convergence in the near-continuum regime, Synthetic Iterative Schemes (SIS) have been developed to accelerate the solution of the Boltzmann equation [22], [23]. The core mechanism of SIS is to couple the microscopic kinetic transport with macroscopic fluid equations, where the macroscopic fluid equations are utilized to guide the evolution of the distribution function. While early SIS versions were often limited to specific flow scenarios or linearized equations, the subsequently proposed General Synthetic Iterative Scheme (GSIS) extended this capability to general rarefied gas flows [24]–[26]. By rigorously deriving the macroscopic equations from the kinetic level, GSIS ensures asymptotic preserving properties and achieves fast convergence across the entire range of Knudsen numbers, typically converging within dozens of iterations.

Alternatively, the Newton method addresses the collision nonlinearity by solving a sequence of Newton equations, each of which is essentially a linearized Boltzmann equation with a residual source term [27]. The competitiveness of this method has been substantially enhanced by the recent developments in the fast Fourier spectral method (FFT), where the complexity of the linearized collision operator is reduced to \(\mathcal{O}(N^4 \log N)\), with \(N\) denoting the number of modes in each velocity direction [27]. The computational cost is significantly lower than that of the quadratic collision operator with complexity of \(\mathcal{O}(MN^4 \log N)\) [13], where \(M\) is the number of quadrature points on the unit sphere. This new efficient FFT to solve the linear collision operator makes the Newton iteration a highly viable and efficient choice.

However, efficiently evaluating the collision operator does not equate to efficiently solving the linear system. Although the linearized collision operator can be evaluated efficiently, the corresponding source iteration still converges slowly for small \({\rm Kn}\), where macroscopic information is propagated inefficiently only through the update of the kinetic equation. This observation motivates a macroscopic synthetic preconditioning strategy, in which a macroscopic system is utilized to accelerate the slow macroscopic components of the Newton correction. However, the macroscopic equations here cannot be directly borrowed from the existing GSIS. The key difference lies in the state around which the equation is linearized. Precisely, nonlinear GSIS is built for the original nonlinear equation [26], and linearized GSIS is derived around a global equilibrium [25], which is spatially uniform. For the Newton correction, it is linearized around a local Maxwellian that varies in space and is updated at each outer step. As a result, a different form of the macroscopic synthetic system should be derived specifically.

In this work, a Newton method equipped with a macroscopic synthetic system (Newton-MS) is developed for the steady-state Boltzmann equation with the quadratic collision operator. The main idea is retaining the modified Newton iteration for the nonlinear collision term, while constructing a macroscopic synthetic system for the inner correction equation to accelerate the convergence of the slow macroscopic components. As discussed above, the corresponding macroscopic equations differ from existing GSIS and are derived for the Newton correction around a spatially varying local Maxwellian. These continuous equations, however, contain several background-dependent terms induced by the spatial variation of the local Maxwellian, which make direct discretization cumbersome. Rather than discretizing them directly, the matrix form of the macroscopic system at the fully discrete level is constructed, starting from the discrete Newton correction equation. This yields a compact macroscopic linear system that can be assembled once per Newton step and reused in the inner iterations, without explicitly deducing the complicated continuous macroscopic synthetic systems.

The spatial discretization is achieved by the discontinuous Galerkin (DG) method, ensuring high-order accuracy and geometric flexibility for multidimensional problems. Extensive numerical experiments, ranging from the 1D Fourier and Couette flows to 2D cavity flows, are presented to validate the efficiency and accuracy of the proposed Newton-MS method.

The rest of this paper is organized as follows. Sec. 2 briefly reviews the Boltzmann equation and the fast Fourier spectral method for the collision operator. The construction of the Newton-MS, especially the deduction of the continuous macroscopic synthetic system, is presented in Sec. 3, while Sec. 4 provides the fully discrete form of Newton-MS. Several numerical examples are provided in Sec. 5 with some concluding remarks and the appendix listed in Sec. 6 and App. 7.

2 Preliminaries↩︎

In this section, several basic properties of the steady-state Boltzmann equation are introduced, and the fast Fourier spectral method, which is utilized to approximate the quadratic collision model, is proposed in Sec. 2.2.

2.1 The steady-state Boltzmann equation↩︎

The dynamics of a rarefied gas are characterized by the evolution of the particle distribution function \(f(t, \boldsymbol{x}, {\boldsymbol{v}})\) in the phase space, where \(t \geqslant 0\) is time, \(\boldsymbol{x}\in \Omega \subset \mathbb{R}^d,d=1,2,3\) is position, and \({\boldsymbol{v}}\in \mathbb{R}^3\) is the microscopic velocity of gas molecules. Here, the steady-state Boltzmann equation is considered as \[\label{eq:Boltz} {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} f(\boldsymbol{x},{\boldsymbol{v}}) = \frac{1}{{\rm Kn}}\mathcal{Q}[f, f](\boldsymbol{x},{\boldsymbol{v}}),\tag{1}\] where \({\rm Kn}\) is the Knudsen number and \(\mathcal{Q}[f, f]\) is the quadratic collision operator, which describes binary collisions between particles and has the following integral form \[\label{eq:quad95col} \mathcal{Q}[f, f](\boldsymbol{x},{\boldsymbol{v}}) = \int_{\mathbb{R}^3} \int_{\mathbb{S}^2} \mathcal{B}({\boldsymbol{v}}- {\boldsymbol{v}}_*, \boldsymbol{\sigma}) \left[ f(\boldsymbol{x},{\boldsymbol{v}}') f(\boldsymbol{x},{\boldsymbol{v}}'_*) - f(\boldsymbol{x},{\boldsymbol{v}}) f(\boldsymbol{x},{\boldsymbol{v}}_*) \right] \; \mathrm{d}\boldsymbol{\sigma} \; \mathrm{d}{\boldsymbol{v}}_*,\tag{2}\] where \(({\boldsymbol{v}}, {\boldsymbol{v}}_*)\) and \(({\boldsymbol{v}}', {\boldsymbol{v}}'_*)\) denote the pre- and post-collision velocity pairs, respectively. Under the constraints of momentum and energy conservation, the post-collision velocities are parameterized by the unit vector \(\boldsymbol{\sigma} \in \mathbb{S}^2\) \[{\boldsymbol{v}}' = \frac{{\boldsymbol{v}}+ {\boldsymbol{v}}_*}{2} + \frac{|{\boldsymbol{v}}- {\boldsymbol{v}}_*|}{2}\boldsymbol{\sigma}, \quad {\boldsymbol{v}}'_* = \frac{{\boldsymbol{v}}+ {\boldsymbol{v}}_*}{2} - \frac{|{\boldsymbol{v}}- {\boldsymbol{v}}_*|}{2}\boldsymbol{\sigma}.\] The collision kernel \(\mathcal{B}\) reflects the intermolecular potential between gas molecules. A widely adopted simplification is the variable hard-sphere (VHS) model [1], which assumes \(\mathcal{B}\) is independent of the scattering angle \(\boldsymbol{\sigma}\) \[\mathcal{B} = C |{\boldsymbol{v}}- {\boldsymbol{v}}_*|^\alpha, \quad \alpha = 2(1 - \omega),\] where \(\omega\) represents the viscosity index of the gas and the constant \(C\) is determined by the molecular diameter. Due to the complex form of the quadratic collision, several simplified collision models are proposed, such as the linearized collision model \[\label{eq:linear95col} \mathcal{L}[g] = \mathcal{Q}[\mathcal{M},g]+\mathcal{Q}[g,\mathcal{M}],\tag{3}\] where \(g\) denotes the distribution acted on by the linearized collision operator. The Maxwellian \(\mathcal{M} = \mathcal{M}_{[\rho,\boldsymbol{u},T]}({\boldsymbol{v}})\) is given by \[\label{eq:Maxwell} \mathcal{M}_{[\rho,\boldsymbol{u},T]}({\boldsymbol{v}}) = \frac{\rho}{\sqrt{2\pi T}^{3}}\exp\left(-\frac{|{\boldsymbol{v}}-\boldsymbol{u}|^2}{2T}\right).\tag{4}\] Here, \(\rho\) is the density, \(\boldsymbol{u}=(u_1,u_2,u_{3})^T\) is the macroscopic velocity, and \(T\) is the temperature. When the reference Maxwellian is associated with a distribution function \(f\), the corresponding macroscopic quantities are given by \[\label{eq:mac95var} \rho = \int_{\mathbb{R}^{3}} f(\boldsymbol{x},{\boldsymbol{v}})\; \mathrm{d}{\boldsymbol{v}},\qquad \boldsymbol{u} =\frac{1}{\rho} \int_{\mathbb{R}^{3}} {\boldsymbol{v}}f(\boldsymbol{x},{\boldsymbol{v}}) \; \mathrm{d}{\boldsymbol{v}},\qquad T= \frac{1}{3\rho}\int_{\mathbb{R}^{3}} |{\boldsymbol{v}}-\boldsymbol{u}|^2 f(\boldsymbol{x},{\boldsymbol{v}}) \; \mathrm{d}{\boldsymbol{v}}.\tag{5}\]

2.2 Fast Fourier spectral approximation of the collision operator↩︎

Due to the complex form of the quadratic collision model \(\mathcal{Q}[f, f]\), numerical methods for it have been extensively studied in recent years. The fast Fourier spectral method (FFS) [12], [13] is adopted here, due to its high efficiency and easy implementation. As stated in [13], the time complexity for FFS to approximate the quadratic collision operator is \(\mathcal{O}(MN^4\log N)\). For VHS models, the linearized collision operator \(\mathcal{L}[g]\) can be approximated with time complexity \(\mathcal{O}(N^4\log N)\), according to the recent research work [27]. In this section, these algorithms will be briefly reviewed.

For FFS, it truncates the microscopic velocity space to a finite domain \(\mathcal{D}_v = [-L, L]^3\). Let \(N\) denote the one-dimensional spectral truncation parameter, so that \(2N\) velocity grid points are used in each direction. The distribution function \(f(\boldsymbol{x},{\boldsymbol{v}})\) is approximated by \[\label{eq:dis95f} f({\boldsymbol{v}})=\frac{1}{(2L)^3}\sum_{\boldsymbol{k}=-N}^{N}\hat{f}_{\boldsymbol{k}}\,e^{i\pi \boldsymbol{k} \cdot {\boldsymbol{v}}/L}, \qquad \hat{f}_{\boldsymbol{k}}=\frac{1}{c_{\boldsymbol{k}}}\left(\frac{L}{N}\right)^3\sum_{\boldsymbol{\ell}=-N}^{N-1} f\!\left(\frac{\boldsymbol{\ell} L}{N}\right)e^{-i\pi \boldsymbol{k}\cdot \boldsymbol{\ell}/N},\tag{6}\] where \(\boldsymbol{k} = (k_1, k_2, k_3)\) and \[\sum_{\boldsymbol{k}=-N}^{N} = \sum_{k_1=-N}^{N} \sum_{k_2=-N}^{N} \sum_{k_3=-N}^{N}, \quad c_{\boldsymbol{k}} = c_{k_1} c_{k_2} c_{k_3}, \quad c_j = \begin{cases} 2, & \text{if } j = \pm N, \\ 1, & \text{otherwise.} \end{cases}\] Moreover, the discretization of the quadratic collision model \(\mathcal{Q}[f, f]\) requires a truncation of the collision kernel by assuming \(\mathcal{B}({\boldsymbol{v}}- {\boldsymbol{v}}_*, \boldsymbol{\sigma}) = 0\) when \(|{\boldsymbol{v}}- {\boldsymbol{v}}_*| > R\), where \(R > 0\) is a problem-dependent parameter [13]. Then the Fourier coefficients of \(\mathcal{Q}\) are calculated as \[\label{eq:coe95Q} \widehat{\mathcal{Q}}_{\boldsymbol{k}} = \sum_{j=1}^J \sum_{m=1}^M C_{j,m,\boldsymbol{k}} \sum_{\boldsymbol{l} = -N}^{N-1} (\alpha^{(1)}_{j,m,\boldsymbol{k}-\boldsymbol{l}} \beta^{(1)}_{j,m,\boldsymbol{l}} + \alpha^{(2)}_{j,m,\boldsymbol{k}-\boldsymbol{l}} \beta^{(2)}_{j,m,\boldsymbol{l}}) \hat{f}_{\boldsymbol{k}-\boldsymbol{l}} \hat{f}_{\boldsymbol{l}} - \sum_{\boldsymbol{l} = -N}^{N-1} \gamma_{\boldsymbol{l}} \hat{f}_{\boldsymbol{k}-\boldsymbol{l}} \hat{f}_{\boldsymbol{l}},\tag{7}\] where all coefficients \(C_{j,m,\boldsymbol{k}}\), \(\alpha^{(1)}_{j,m,\boldsymbol{k}}\), \(\alpha^{(2)}_{j,m,\boldsymbol{k}}\), \(\beta^{(1)}_{j,m,\boldsymbol{k}}\), \(\beta^{(2)}_{j,m,\boldsymbol{k}}\) and \(\gamma_{\boldsymbol{l}}\) can be precomputed, and we refer to [13] for the detailed calculation of these coefficients. Thus, if the fast Fourier transform (FFT) is adopted to compute the convolution, the final computational cost of 7 is \(\mathcal{O}(JM N^3 \log N)\), where \(J\) is the number of quadrature points for a one-dimensional integral over \([0,R]\). Generally speaking, \(R\) is chosen to be proportional to \(L\), and it is expected that \(J = \mathcal{O}(N)\). Thus the computational cost can also be written as \(\mathcal{O}(M N^4 \log N)\) [13].

For the VHS models, the Fourier coefficients of the linearized collision operator \(\mathcal{L}[g]\) can be computed more efficiently. The algorithm in [27] provides the following form of \(\widehat{\mathcal{L}}_{\boldsymbol{k}}\): \[\label{eq:coe95L} \widehat{\mathcal{L}}_{\boldsymbol{k}} = \sum_{j=1}^J \alpha_{j,\boldsymbol{k}} \sum_{\boldsymbol{l} = -N}^{N-1} \mathcal{M}_{\boldsymbol{l}}^h \psi_{j,\boldsymbol{l}} e^{-i\pi \boldsymbol{k} \cdot \boldsymbol{l}/N} - \sum_{\boldsymbol{l} = -N}^{N-1} \gamma_{\boldsymbol{l}} (\hat{f}_{\boldsymbol{k}-\boldsymbol{l}} \hat{\mathcal{M}}_{\boldsymbol{l}} + \hat{\mathcal{M}}_{\boldsymbol{k}-\boldsymbol{l}} \hat{f}_{\boldsymbol{l}}),\tag{8}\] where \(\hat{\mathcal{M}}_{\boldsymbol{l}}\) is the Fourier coefficient of the local Maxwellian 4 , and \[\begin{gather} \psi_{j,\boldsymbol{l}} = \left(\frac{\pi}{N}\right)^3 \sum_{\boldsymbol{m}} \beta_{j,\boldsymbol{m}} r_{\boldsymbol{l}-\boldsymbol{m}} \mathcal{M}^c_{\boldsymbol{m}}, \qquad r_{\boldsymbol{l}} = g\left( \frac{\boldsymbol{l}L}{N} \right) \Bigg/\mathcal{M}\left( \frac{\boldsymbol{l}L}{N} \right), \\ \mathcal{M}_{\boldsymbol{m}}^h = \frac{\rho}{(\pi T)^{3/2}} \exp \left( -\frac{1}{T} \left| \frac{\boldsymbol{m} L}{N} - \boldsymbol{u} \right|^2 \right), \qquad \mathcal{M}_{\boldsymbol{m}}^c = \frac{\rho}{(\pi T)^{3/2}} e^{-\frac{|\boldsymbol{m}|^2 L^2}{N^2 T}}. \end{gather}\] The coefficients \(\alpha_{j,\boldsymbol{k}}\) and \(\beta_{j,\boldsymbol{m}}\) can again be precomputed. With FFT, this leads to a total computational cost \(\mathcal{O}(N^4 \log N)\), removing the dependence on \(M\), which is the degree of freedom for the integration on the unit sphere compared with the quadratic collision operator.

With the lower computational cost of the linearized collision operator, one can achieve a cheaper model by replacing \(\mathcal{Q}[f,f]\) with \(\mathcal{L}[g]\), while still maintaining much better accuracy than BGK-type models. Moreover, as briefly tested in [27], the lower computational cost of \(\mathcal{L}[g]\) can also help develop efficient algorithms for the steady-state Boltzmann equation 1 , whose efficiency will be further improved in the following sections, especially in the case of small Knudsen number.

3 Construction of Newton-macroscopic synthetic system↩︎

Solving the steady-state Boltzmann equation efficiently is hindered by two primary challenges: the prohibitive computational cost of the high-dimensional collision integral, and the slow convergence caused by the stiffness of the transport term in the near-continuum regime. To address these issues, a Newton iteration accelerated by a macroscopic synthetic system in the inner iteration (Newton-MS) is constructed. Precisely, the nonlinear Boltzmann equation is solved by an outer modified Newton iteration, and each Newton correction equation is solved by an inner source iteration preconditioned by a macroscopic synthetic system. The linearized collision operator is used in the Newton correction equation to reduce the cost of collision evaluation, while the macroscopic correction accelerates the propagation of large-scale flow information when \({\rm Kn}\) is small.

For the modified Newton’s method to solve the steady-state Boltzmann equation, it aims to seek the numerical solution to 9 \[\label{eq:ss95bolt} \mathcal{R}(f) := {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} f - \frac{1}{{\rm Kn}}\mathcal{Q}[f, f] = 0.\tag{9}\] A similar Newton method proposed in [27] is utilized here, where the distribution function is updated iteratively by \[\label{eq:Newton95update} f^{(n+1)} = f^{(n)} - g^{(n)},\tag{10}\] where the correction \(g^{(n)}\) is solved from the linearized equation \[\label{eq:Newton95eq} {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g^{(n)} - \frac{1}{{\rm Kn}}\Big(\mathcal{Q}[f^{(n)},g^{(n)}]+\mathcal{Q}[g^{(n)},f^{(n)}]\Big) = r^{(n)},\tag{11}\] with \[\label{eq:r} r^{(n)} = \mathcal{R}(f^{(n)})\tag{12}\] the residual at the \(n\)-th Newton iterate. This Newton iteration 10 to update the distribution function is called the outer Newton iteration. The computational cost can be extremely high when numerically solving 11 with an iterative method, regardless of which method is used, due to the high cost of computing the binary collision operator. To solve this, the quadratic collision term \(\mathcal{Q}[f^{(n)},g^{(n)}]+\mathcal{Q}[g^{(n)},f^{(n)}]\) is approximated with a linearized collision operator around the local Maxwellian, with a similar method utilized in [27]. This strategy is conceptually similar to the quasi-Newton method by replacing the Jacobian matrix with a computationally efficient approximation. Therefore, the governing equation of the correction \(g^{(n)}\) is reduced to \[\label{eq:Newton95eq2} {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g^{(n)} -\frac{1}{{\rm Kn}}\mathcal{L}^{(n)}[g^{(n)}] = r^{(n)},\tag{13}\] where the linearized operator \(\mathcal{L}^{(n)}[\cdot]\) is defined around the local Maxwellian as \[\label{eq:linear95coll} \mathcal{L}^{(n)}[\cdot] = \mathcal{Q}[\mathcal{M}^{(n)}, \cdot]+\mathcal{Q}[\cdot,\mathcal{M}^{(n)}].\tag{14}\] Here, \(\mathcal{M}^{(n)}\) is the Maxwellian corresponding to \(f^{(n)}\). It has been demonstrated in [27] that this modified Newton’s method shows fast convergence rates and can greatly reduce the computational cost at the same time.

Within each Newton’s iteration, to solve the linear equation 13 , the iterative method is also utilized, which is called the inner iteration. For example, the classical source iteration (SI) was employed in [27], which updates the numerical solution by \[\label{eq:SI} {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g^{(n,l+1)} + \frac{{\nu}}{{\rm Kn}} g^{(n,l+1)} = \frac{1}{{\rm Kn}}\mathcal{L}^{(n)}[g^{(n,l)}] + \frac{{\nu}}{{\rm Kn}} g^{(n,l)} + r^{(n)},\tag{15}\] where \(n\) stands for the index of the outer Newton iteration, and \(l\) stands for the index of the inner source iteration. Here, the penalization terms involving \(\nu\) are introduced to obtain a stable source iteration, where \(\nu\) is usually chosen as the local collision frequency depending on the density and temperature. The resulting outer-inner solver is referred to as the modified Newton method combined with source iteration (Newton–SI). However, it is well known that SI leads to slow convergence when the Knudsen number is small [10], [25], [28], [29].

Improving the computational efficiency of SI is full of challenges, and several numerical methods have been proposed to accelerate it, especially in the near-continuum regime. For example, a macroscopic synthetic preconditioner is constructed by solving a macroscopic system, and a similar macroscopic coupling principle is also proposed in the general synthetic iterative scheme (GSIS) [24], [25].

The general idea of GSIS is to accelerate the convergence of macroscopic variables by solving an NSF-like system, where the stress tensor and the heat flux include higher-order contributions from the microscopic solution in the previous iteration. Precisely, the synthetic equations are the full nonlinear NSF system with high-order corrections in nonlinear GSIS [26], while the linearized counterpart around a global equilibrium is adopted in [25]. The main contribution of this work is to construct an efficient macroscopic synthetic preconditioner to solve 15 by coupling the macroscopic governing equations.

3.1 Macroscopic moment system for the Newton correction equation↩︎

To improve the convergence of SI in the near-continuum regime, we follow the idea of adopting the macroscopic system to accelerate convergence, but the specific form of the macroscopic equations is different from GSIS. Compared with GSIS, the correction equation 13 in this work is linearized around the local Maxwellian \(\mathcal{M}^{(n)}(\boldsymbol{x}, {\boldsymbol{v}})\), where the macroscopic variables, such as the density \(\rho^{(n)}\) and macroscopic velocity \(\boldsymbol{u}^{(n)}\), vary in space and will only be refreshed at the outer Newton step when the distribution function is updated by 10 . This may bring a different macroscopic synthetic preconditioner, which will be derived through the Chapman-Enskog expansion of the correction \(g^{(n)}\) around \(\mathcal{M}^{(n)}\).

Before presenting the macroscopic synthetic preconditioner, the local equilibrium \(\mathcal{E}^{(n)}\) of the correction \(g^{(n)}\) is introduced. Defining the peculiar velocity \(\boldsymbol{c}^{(n)} = {\boldsymbol{v}}- \boldsymbol{u}^{(n)}(\boldsymbol{x})\), the local Maxwellian \(\mathcal{M}^{(n)}\) of the distribution function \(f^{(n)}\) in 14 can be rewritten as \[\label{eq:local95Max} \mathcal{M}^{(n)}(\boldsymbol{x},{\boldsymbol{v}}) = \frac{\rho^{(n)}(\boldsymbol{x})}{(2\pi T^{(n)}(\boldsymbol{x}))^{3/2}} \exp \left( -\frac{|\boldsymbol{c}^{(n)}|^2}{2T^{(n)}(\boldsymbol{x})} \right),\tag{16}\] where \(\rho^{(n)}(\boldsymbol{x})\), \(\boldsymbol{u}^{(n)}(\boldsymbol{x})\), and \(T^{(n)}(\boldsymbol{x})\) denote the macroscopic moments of \(f^{(n)}\) through 5 . During the inner iteration in one Newton step, these macroscopic variables are known and fixed. Then the kernel space of the linearized collision operator 14 holds the form as \[\label{eq:ker95linear} \operatorname{ker} \mathcal{L}^{(n)} = \{\varphi(\boldsymbol{c}^{(n)})^T \alpha(\boldsymbol{x}) \mathcal{M}^{(n)}(\boldsymbol{x},{\boldsymbol{v}}) \mid \alpha: \Omega \rightarrow \mathbb{R}^5\},\tag{17}\] where \[\label{eq:phi} \varphi(\boldsymbol{c}^{(n)}) = \left( 1, \boldsymbol{c}^{(n)}, \displaystyle \frac{|\boldsymbol{c}^{(n)}|^2 - 3T^{(n)}}{2} \right)^T.\tag{18}\] Then, the local equilibrium \(\mathcal{E}^{(n)}(\boldsymbol{x},{\boldsymbol{v}})\) of the correction \(g^{(n)}\) is defined as the projection of \(g^{(n)}\) onto \(\operatorname{ker} \mathcal{L}^{(n)}\) \[\label{eq:linear95Max} \mathcal{E}^{(n)}(\boldsymbol{x},{\boldsymbol{v}}) \triangleq \varphi(\boldsymbol{c}^{(n)})^T {m}^{(n)}(\boldsymbol{x}) \mathcal{M}^{(n)}(\boldsymbol{x},{\boldsymbol{v}}),\tag{19}\] where \(\boldsymbol{m}^{(n)}\) is the vector of local macroscopic variables defined as \[\label{eq:local95mac} {m}^{(n)}(\boldsymbol{x}) = \begin{pmatrix} \displaystyle \frac{\delta \rho^{(n)}(\boldsymbol{x})}{\rho^{(n)}(\boldsymbol{x})}, & \displaystyle \frac{\delta \boldsymbol{u}^{(n)}(\boldsymbol{x})}{T^{(n)}(\boldsymbol{x})}, & \displaystyle \frac{\delta T^{(n)}(\boldsymbol{x})}{(T^{(n)}(\boldsymbol{x}))^2} \end{pmatrix}^T,\tag{20}\] with \(\delta \rho\), \(\delta \boldsymbol{u}\), and \(\delta T\) being the perturbations of macroscopic variables related to \({g}^{(n)}\) \[\label{eq:macro95g} {m}^{(n)}(\boldsymbol{x}) = \int_{\mathbb{R}^3} \overline{\varphi}(\boldsymbol{c}^{(n)}) g^{(n)} \; \mathrm{d}{\boldsymbol{v}}, \qquad \overline{\varphi}(\boldsymbol{c}^{(n)}) = \left(\frac{1}{{\rho}^{(n)}}, \frac{\boldsymbol{c}^{(n)}}{{\rho}^{(n)}{T}^{(n)}}, \frac{|{c}^{(n)}|^2 - 3{T}^{(n)}}{3{\rho}^{(n)}({T}^{(n)})^2}\right)^T.\tag{21}\] It is easy to verify that \[\label{eq:conver95g} \int_{\mathbb{R}^3} \phi({\boldsymbol{v}}) g^{(n)} \; \mathrm{d}{\boldsymbol{v}}= \int_{\mathbb{R}^3} \phi({\boldsymbol{v}}) \mathcal{E}^{(n)} \; \mathrm{d}{\boldsymbol{v}}, \qquad \phi({\boldsymbol{v}}) = (1, {\boldsymbol{v}}, |{\boldsymbol{v}}|^2)^T.\tag{22}\] Defining \(\mathcal{P}\) as the orthogonal projection operator from any distribution function to \(\operatorname{ker} \mathcal{L}\), it holds that \[\label{eq:operator95P} \int_{\mathbb{R}^3} \phi({\boldsymbol{v}}) (I - \mathcal{P}) f \; \mathrm{d}{\boldsymbol{v}}= 0, \qquad \mathcal{P}g^{(n)} = \mathcal{E}^{(n)}.\tag{23}\] The detailed derivation of the local equilibrium of \(g^{(n)}\) can be found in [30].

In the following, the macroscopic system associated with the Newton correction will be derived with a Chapman-Enskog type expansion of the correction \(g^{(n)}\) in terms of the Knudsen number, and the outer-iteration superscript “\((n)\)” is omitted for simplicity. Thus, the linearized equation 13 is reduced to \[\label{eq:Kn95scaled} {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g - \frac{1}{{\rm Kn}} \mathcal{L}[g] = r.\tag{24}\] Moreover, with the local Maxwelilan 19 , the Chapman-Enskog ansatz of the correction \(g\) has the form \[\label{eq:ce95expan} g = \mathcal{E} + {\rm Kn}g_1 + {\rm Kn}^2 g_2,\tag{25}\] where \(g_1\) represents the first-order Navier-Stokes correction, and \(g_2\) denotes the high-order non-equilibrium remainder. With the definition of \(\mathcal{P}\), it holds that \[\label{eq:g195g2} \mathcal{P} g_1 = 0, \qquad \mathcal{P} g_2 = 0.\tag{26}\] Substituting 25 into 24 , with \(\mathcal{L}[\mathcal{E}] = 0\), we can obtain \[\label{eq:CE} {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} \mathcal{E} + {\rm Kn}\, {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g_1 + {\rm Kn}^2 {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g_2 - \mathcal{L}[g_1] - {\rm Kn}\, \mathcal{L}[g_2] = r.\tag{27}\] To obtain \(g_1\), by comparing the terms in order \(\mathcal{O}(1)\), we obtain \[\label{eq:order951} \mathcal{L}[g_1] = {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} \mathcal{E} - r.\tag{28}\] Modding out the kernel of \(\mathcal{L}\) by taking the operator \((I-\mathcal{P})\), then \(g_1\) is derived as \[\label{eq:ce95expan2} g_1 = \mathcal{L}^{-1}(I - \mathcal{P})\!\left({\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} \mathcal{E} - r\right),\tag{29}\] where \(\mathcal{L}^{-1}\) is the pseudoinverse of \(\mathcal{L}\). From 23 , it holds for \(g_1\) and \(g_2\) that \[\label{eq:con95g1} \int_{\mathbb{R}^3} \phi({\boldsymbol{v}}) g_1 \; \mathrm{d}{\boldsymbol{v}}= \int_{\mathbb{R}^3} \phi({\boldsymbol{v}}) g_2 \; \mathrm{d}{\boldsymbol{v}}= 0.\tag{30}\] Moreover, since \(\mathcal{L}^{-1} (I-\mathcal{P}) = \mathcal{L}^{-1}\), the projection \((I-\mathcal{P})\) will be omitted below.

The macroscopic equations for \(\delta \rho\), \(\delta \boldsymbol{u}\) and \(\delta T\) are obtained by taking moments of 27 with respect to the collision invariants \(\phi({\boldsymbol{v}})\). With the conservation property of the linearized collision operator and the first-order correction 29 , the resulting macroscopic equations can be written as \[\tag{31} \begin{align} \tag{32} \mathcal{S}_1 &= \int_{\mathbb{R}^3} r \; \mathrm{d}{\boldsymbol{v}}, \\ \tag{33} \mathcal{S}_2 + {\rm Kn}\, \nabla_{\boldsymbol{x}} \cdot \boldsymbol{\sigma} &= \int_{\mathbb{R}^3} {\boldsymbol{v}}r \; \mathrm{d}{\boldsymbol{v}}- {\rm Kn}^2 \nabla_{\boldsymbol{x}} \cdot \int_{\mathbb{R}^3} {\boldsymbol{v}}\otimes {\boldsymbol{v}}g_2 \; \mathrm{d}{\boldsymbol{v}}, \\ \tag{34} \mathcal{S}_3 + 2{\rm Kn}\, \nabla_{\boldsymbol{x}} \cdot \left( \boldsymbol{\sigma} \boldsymbol{u}+ \boldsymbol{q}\right) &= \int_{\mathbb{R}^3} |{\boldsymbol{v}}|^2\, r \; \mathrm{d}{\boldsymbol{v}}- {\rm Kn}^2 \nabla_{\boldsymbol{x}} \cdot \int_{\mathbb{R}^3} |{\boldsymbol{v}}|^2 {\boldsymbol{v}}\, g_2\; \mathrm{d}{\boldsymbol{v}}, \end{align}\] with \[\label{eq:eq95S} \begin{align} \mathcal{S}_1= &\nabla_{\boldsymbol{x}} \cdot \left( \boldsymbol{u}\, \delta\rho + \rho\, \delta\boldsymbol{u}\right), \\ \mathcal{S}_2 = &\nabla_{\boldsymbol{x}} \cdot \Big[ \delta\rho\, (\boldsymbol{u}\otimes \boldsymbol{u}+ T \mathbf{I}) + \rho\, (\boldsymbol{u}\otimes \delta \boldsymbol{u}+ \delta \boldsymbol{u}\otimes \boldsymbol{u}) + \rho\, \delta T\, \mathbf{I}\Big], \\ \mathcal{S}_3=& \nabla_{\boldsymbol{x}} \cdot \Big[ \delta\rho \boldsymbol{u}(|\boldsymbol{u}|^2 + 5T) + \rho\, \delta\boldsymbol{u}\, (|\boldsymbol{u}|^2 + 5T) + 2\rho \boldsymbol{u}(\boldsymbol{u}\cdot \delta\boldsymbol{u}) + 5\rho\, \boldsymbol{u}\, \delta T \Big]. \end{align}\tag{35}\] Here, the variables \(\boldsymbol{\sigma}\) and \(\boldsymbol{q}\) are defined as \[\label{eq:sigma95q} \boldsymbol{\sigma} = \int_{\mathbb{R}^3} \left( \boldsymbol{c} \otimes \boldsymbol{c} - \frac{|\boldsymbol{c}|^2}{3}\mathbf{I} \right) g_1 \; \mathrm{d}{\boldsymbol{v}}, \qquad \boldsymbol{q}= \int_{\mathbb{R}^3} \frac{1}{2} |\boldsymbol{c}|^2 \boldsymbol{c}\, g_1 \; \mathrm{d}{\boldsymbol{v}},\tag{36}\] with \(\otimes\) denotes the tensor product of two vectors. The left-hand side of 31 has a similar form to the standard linearized Navier-Stokes-Fourier equations (NSF), which are approximated in the reference state \((\rho, \boldsymbol{u}, T)\). But the variables \(\boldsymbol{\sigma}\) and \(\boldsymbol{q}\) defined in 36 depend not simply on \(\delta \boldsymbol{u}\) and \(\delta T\) as the standard NSF, but also on the macroscopic variables \(\rho, \boldsymbol{u}\) and \(\theta\) due to the non-homogeneous Maxwellian \(\mathcal{M}^{(n)}\) in the outer Newton iteration. To obtain \(\boldsymbol{\sigma}\) and \(\boldsymbol{q}\), they are split into two parts \[\label{eq:split95s95q} \boldsymbol{\sigma} = \boldsymbol{\sigma}^{(\nabla)} - \boldsymbol{\sigma}^{(r)}, \qquad \boldsymbol{q} = \boldsymbol{q}^{(\nabla)} - \boldsymbol{q}^{(r)},\tag{37}\] where \[\begin{align} \tag{38} \boldsymbol{\sigma}^{(\nabla)} &= \int_{\mathbb{R}^3} \left( \boldsymbol{c}\otimes \boldsymbol{c} - \frac{|\boldsymbol{c}|^2}{3} \mathbf{I} \right) \mathcal{L}^{-1} ({\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} \mathcal{E}) \; \mathrm{d}{\boldsymbol{v}}, \qquad \boldsymbol{q}^{(\nabla)} = \int_{\mathbb{R}^3} \frac{|\boldsymbol{c}|^2}{2} \boldsymbol{c}\, \mathcal{L}^{-1} ({\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} \mathcal{E}) \; \mathrm{d}{\boldsymbol{v}}, \\ \tag{39} \boldsymbol{\sigma}^{(r)} &= \int_{\mathbb{R}^3} \left( \boldsymbol{c}\otimes \boldsymbol{c} - \frac{|\boldsymbol{c}|^2}{3} \mathbf{I} \right) \mathcal{L}^{-1} r \; \mathrm{d}{\boldsymbol{v}}, \qquad \boldsymbol{q}^{(r)} = \int_{\mathbb{R}^3} \frac{|\boldsymbol{c}|^2}{2} \boldsymbol{c}\, \mathcal{L}^{-1} r \; \mathrm{d}{\boldsymbol{v}}. \end{align}\] Here, \(\boldsymbol{\sigma}^{(\nabla)}\) and \(\boldsymbol{q}^{(\nabla)}\) are calculated as \[\begin{align} \label{eq:sigma95q95full} \boldsymbol{\sigma}^{(\nabla)} &= \frac{C_{2,0}}{\nu}\,\mathbf{S} + \frac{C_{2,1}}{\nu}\,\mathbf{R}, \qquad \boldsymbol{q}^{(\nabla)} = -\frac{1}{2\nu} \left({C_{1,1}}\,\,\boldsymbol{P} + {C_{1,2}}\,\,\boldsymbol{Q} \right), \end{align}\tag{40}\] where \[\tag{41} \begin{align} \tag{42} \mathbf{S} &= \frac{1}{T}(\nabla_{\boldsymbol{x}}\delta \boldsymbol{u})_{\mathrm{stf}} + \left(\frac{\delta \boldsymbol{u}}{T}\otimes\frac{\nabla_{\boldsymbol{x}}\rho}{\rho}\right)_{\mathrm{stf}} + \left(\frac{\delta\rho}{\rho T} - \frac{5\delta T}{2T^2} \right)\left(\nabla_{\boldsymbol{x}}\boldsymbol{u}\right)_{\mathrm{stf}} \notag \\ &\quad - \frac{5}{2T^2}\left(\delta\boldsymbol{u}\otimes\nabla_{\boldsymbol{x}}T\right)_{\mathrm{stf}} + \left(\frac{\delta\boldsymbol{u}}{T}\otimes(\boldsymbol{u}\cdot\nabla_{\boldsymbol{x}})\boldsymbol{u}\right)_{\mathrm{stf}}, \\ \tag{43} \mathbf{R} &= \frac{\rho\delta T}{3T}\,\left(\nabla_{\boldsymbol{x}}\boldsymbol{u}\right)_{\mathrm{stf}} + \frac{\rho}{2T}\,\left(\delta\boldsymbol{u}\otimes\nabla_{\boldsymbol{x}}T\right)_{\mathrm{stf}}, \end{align}\] and \[\tag{44} \begin{align} \tag{45} \boldsymbol{P} &= \rho T\,\nabla_{\boldsymbol{x}}\delta T + T\delta T\,\nabla_{\boldsymbol{x}}\rho + T\delta\rho\,\nabla_{\boldsymbol{x}} T - 5\rho\delta T\,\nabla_{\boldsymbol{x}} T \notag \\ & \quad + \frac{2\rho T}{5}\bigl[(\delta\boldsymbol{u}\cdot\nabla_{\boldsymbol{x}})\boldsymbol{u} + (\nabla_{\boldsymbol{x}}\boldsymbol{u})^{\top}\delta\boldsymbol{u} + (\nabla_{\boldsymbol{x}}\cdot\boldsymbol{u})\,\delta\boldsymbol{u}\bigr] \\ \notag &\quad + \rho\delta T\,(\boldsymbol{u}\cdot\nabla_{\boldsymbol{x}})\boldsymbol{u} + \rho(\boldsymbol{u}\cdot\nabla_{\boldsymbol{x}} T)\,\delta\boldsymbol{u}\\ \tag{46} \boldsymbol{Q} &= \frac{\rho\delta T}{2}\,\nabla_{\boldsymbol{x}} T. \end{align}\] Here, \(\left(\mathbf{X}\right)_{\mathrm{stf}}\) denotes the symmetric and trace-free part of the matrix \(\mathbf{X}\), defined by \((\mathbf{X} + \mathbf{X}^T)/2 - \operatorname{tr}(\mathbf{X})\mathbf{I}/3\). The detailed deduction is listed in App. 7. For the spatially homogeneous background where \(\nabla_{\boldsymbol{x}}\rho\), \(\nabla_{\boldsymbol{x}}\boldsymbol{u}\) and \(\nabla_{\boldsymbol{x}}T\) all equal zero, 41 and 44 reduce to the standard constitutive laws \[\label{eq:sigma95q95frozen} \boldsymbol{\sigma}^{(\nabla)} = \frac{C_{2,0}}{\nu T} \left(\nabla_{\boldsymbol{x}} \delta\boldsymbol{u}\right)_{\mathrm{stf}}, \qquad \boldsymbol{q}^{(\nabla)} = -\frac{C_{1,1}}{2\nu}\rho T\, \nabla_{\boldsymbol{x}} \delta T.\tag{47}\] Substituting 37 in 31 , with 40 , the macroscopic system 31 is reduced into \[\tag{48} \begin{align} \tag{49} \mathcal{S}_1 &= \int_{\mathbb{R}^3} r \; \mathrm{d}{\boldsymbol{v}}, \\ \tag{50} \mathcal{S}_2 + {\rm Kn}\nabla_{\boldsymbol{x}}\cdot \boldsymbol{\sigma}^{(\nabla)} &= {\rm Kn}\nabla_{\boldsymbol{x}}\cdot \boldsymbol{\sigma}^{(r)} + \int_{\mathbb{R}^3} {\boldsymbol{v}}r \; \mathrm{d}{\boldsymbol{v}}- \mathcal{G}_{2,1},\\ \tag{51} \mathcal{S}_3 + 2{\rm Kn}\, \nabla_{\boldsymbol{x}} \cdot \left( \boldsymbol{\sigma}^{(\nabla)} \boldsymbol{u}+ \boldsymbol{q}^{(\nabla)}\right) &= 2{\rm Kn}\, \nabla_{\boldsymbol{x}} \cdot \left( \boldsymbol{\sigma}^{(r)} \boldsymbol{u}+ \boldsymbol{q}^{(r)}\right) + \int_{\mathbb{R}^3} |{\boldsymbol{v}}|^2\, r \; \mathrm{d}{\boldsymbol{v}}- \mathcal{G}_{2,2}, \end{align}\] with \[\label{eq:G2} \mathcal{G}_{2,1} = {\rm Kn}^2 \nabla_{\boldsymbol{x}} \cdot \int_{\mathbb{R}^3} {\boldsymbol{v}}\otimes {\boldsymbol{v}}g_2 \; \mathrm{d}{\boldsymbol{v}},\qquad \mathcal{G}_{2,2} = {\rm Kn}^2 \nabla_{\boldsymbol{x}} \cdot \int_{\mathbb{R}^3} |{\boldsymbol{v}}|^2 {\boldsymbol{v}}\, g_2\; \mathrm{d}{\boldsymbol{v}}\tag{52}\] where the left parts are the linear system with the unknown variables \(\delta\rho\), \(\delta \boldsymbol{u}\) and \(\delta T\), who depend on \(\rho\), \(\boldsymbol{u}\) and \(T\), which keep constant during the inner iteration. The right parts are unknown terms left to be closed, which will be introduced in detail in the following section. Then, the macroscopic system 31 will be completed, which form a linear system of \(\delta \rho\), \(\delta \boldsymbol{u}\), and \(\delta T\), while most of the coefficients are spatially dependent.

3.1.0.1 BGK-type approximation of the coefficients \(C_{i,j}\)

For the linearized collision operator, to obtain the coefficients \(C_{i,j}\) in 40 , the operator \(\mathcal{L}^{-1}\) should be calculated exactly, which is quite difficult [30]. In the numerical implementation, the BGK model is utilized as a penalty term to obtain \(C_{i,j}\). Precisely, 24 is reformulated by inserting a simpler BGK collision operator: \[\label{eq:BGK95penalized} {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g - \frac{1}{{\rm Kn}}\mathcal{L}_{\rm BGK}[g] = r + \frac{1}{{\rm Kn}}\Big(\mathcal{L}[g] - \mathcal{L}_{\rm BGK}[g]\Big), \qquad \mathcal{L}_{\rm BGK} = \nu(\mathcal{E} - g).\tag{53}\] Thus, define \[\label{eq:tilde95r} \tilde{r} = r + {\rm Kn}^{-1} \big(\mathcal{L}[g] - \mathcal{L}_{\rm BGK}[g]\big).\tag{54}\] It is easy to verify that utilizing \(\tilde{r}\) instead of \(r\) in 53 , the deduction of the macroscopic synthetic system is the same as 48 , due to \[\int \phi({\boldsymbol{v}}) r \; \mathrm{d}{\boldsymbol{v}}= \int \phi({\boldsymbol{v}}) \tilde{r} \; \mathrm{d}{\boldsymbol{v}}.\] Moreover, the pseudoinverse of the BGK operator has the form \[\label{eq:BGK95inv} \mathcal{L}_{\rm BGK}^{-1} =-\frac{1}{\nu} (I-\mathcal{P}),\tag{55}\] thus, the coefficients \(C_{i,j}\) can be obtained explicitly as \[C_{2,0} = -\rho T^2,\qquad C_{2,1} = -7 T,\qquad C_{1,1} = 5,\qquad C_{1,2} = 70.\] The detailed deduction is presented in App. 7.2.

3.2 Closure and synthetic iteration↩︎

To obtain the right-side terms of 48 , we recall the numerical scheme of the inner iteration 15 . Omitting the superscript \((n)\) in the outer iteration, 15 is reduced into \[\label{eq:SI95update} {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g^{(l+1)} + \frac{{\nu}}{{\rm Kn}} g^{(l+1)} = \frac{1}{{\rm Kn}}\mathcal{L}^{(n)}[g^{(l)}] + \frac{{\nu}}{{\rm Kn}} g^{(l)} + r,\tag{56}\] with \(r\) defined in 12 , which keeps constant in the inner iteration. Moreover, the background macroscopic variables \(\rho,\boldsymbol{u}\), and \(T\) also remain constant. Let \[\label{eq:g195l} g_1^{(\nabla,l)}= \mathcal{L}^{-1}\left({\boldsymbol{v}}\cdot\nabla_{\boldsymbol{x}}\mathcal{E}^{(l)}\right), \qquad g_1^{(r)} = g_1^{(r,l)} = \mathcal{L}^{-1}r, \qquad g_1^{(l)} = g_1^{(\nabla,l)} - g_1^{(r)},\tag{57}\] where \(\mathcal{E}^{(l)}\) is the corresponding local equilibrium of \(g^{(l)}\). Here, the superscript \((l)\) is omitted in \(g_1^{(r, l)}\) due to the constancy of the residual \(r\) 12 in each inner iteration. With the Chapman-Enskog decomposition 25 and 29 , it holds at the \(l\)-th inner iteration that \[\label{eq:g2} {\rm Kn}^2 g_2^{(l)} = g^{(l)} - \mathcal{E}^{(l)} - {\rm Kn}\,g_1^{(\nabla,l)} + {\rm Kn}g_1^{(r,l)}.\tag{58}\] Moreover, due to 30 , we can derive that \[\label{eq:G2951} \int_{\mathbb{R}^3} {\boldsymbol{v}}\otimes {\boldsymbol{v}}g_2^{(l)} \; \mathrm{d}{\boldsymbol{v}}= \int \left( \boldsymbol{c} \otimes \boldsymbol{c} - \frac{|\boldsymbol{c}|^2}{3}\mathbf{I} \right) g_2^{(l)} \; \mathrm{d}{\boldsymbol{v}}, \qquad \int_{\mathbb{R}^3} |{\boldsymbol{v}}|^2{\boldsymbol{v}}g_2^{(l)} \; \mathrm{d}{\boldsymbol{v}}= \int_{\mathbb{R}^3} |\boldsymbol{c}|^2 \boldsymbol{c} g_2^{(l)} \; \mathrm{d}{\boldsymbol{v}}.\tag{59}\]

Since the terms 52 are small variables, they are closed with variables obtained at the \(l\)-th inner iteration to close 48 . Precisely, substituting 58 into 52 , with 59 , 39 , then 48 is reduced into \[\tag{60} \begin{align} \tag{61} \mathcal{S}_1 &= \int_{\mathbb{R}^3} r \; \mathrm{d}{\boldsymbol{v}}, \\ \tag{62} \mathcal{S}_2 + {\rm Kn}\nabla_{\boldsymbol{x}}\cdot \boldsymbol{\sigma}^{(\nabla)} &= \int_{\mathbb{R}^3} {\boldsymbol{v}}r \; \mathrm{d}{\boldsymbol{v}}- {\rm Kn}\nabla_{\boldsymbol{x}} \cdot \left(\boldsymbol{\sigma}^{(\mathrm{neq},l)} - \boldsymbol{\sigma}^{(\nabla,l)}\right),\\ \tag{63} \mathcal{S}_3 + 2{\rm Kn}\, \nabla_{\boldsymbol{x}} \cdot \left( \boldsymbol{\sigma}^{(\nabla)} \boldsymbol{u}+ \boldsymbol{q}^{(\nabla)}\right) &= \int_{\mathbb{R}^3} |{\boldsymbol{v}}|^2\, r \; \mathrm{d}{\boldsymbol{v}}- 2{\rm Kn}\nabla_{\boldsymbol{x}} \cdot \left[ \left( \boldsymbol{\sigma}^{(\mathrm{neq},l)} \boldsymbol{u}+ \boldsymbol{q}^{(\mathrm{neq},l)} \right) - \left( \boldsymbol{\sigma}^{(\nabla,l)} \boldsymbol{u}+ \boldsymbol{q}^{(\nabla,l)} \right) \right]. \end{align}\] Here, the terms 39 are canceled and there is no need to compute \(\mathcal{L}^{-1}r\). \(\boldsymbol{\sigma}^{(\mathrm{neq},l)}\) and \(\boldsymbol{q}^{(\mathrm{neq},l)}\) are non-equilibrium moments defined as \[\label{eq:sigma95q95neq} \boldsymbol{\sigma}^{(\mathrm{neq},l)} = \frac{1}{{\rm Kn}} \int_{\mathbb{R}^3} \left( \boldsymbol{c}\otimes\boldsymbol{c}-\frac{|\boldsymbol{c}|^2}{3}\mathbf{I} \right) \left( g^{(l)}- \mathcal{E}^{(l)}\right)\,\; \mathrm{d}{\boldsymbol{v}},\qquad \boldsymbol{q}^{(\mathrm{neq},l)}= \frac{1}{{\rm Kn}} \int_{\mathbb{R}^3} \frac{1}{2}|\boldsymbol{c}|^2\boldsymbol{c}\,\left( g^{(l)}- \mathcal{E}^{(l)}\right)\,\; \mathrm{d}{\boldsymbol{v}}.\tag{64}\] They are calculated directly by 64 in the \(l\)-th inner iteration since \(g^{(l)}\) and \(\mathcal{E}^{(l)}\) are already known for the moment. \(\boldsymbol{\sigma}^{(\nabla, l)}\) and \(\boldsymbol{q}^{(\nabla,l)}\) are defined in 38 with \(\mathcal{E}\) replaced by \(\mathcal{E}^{(l)}\), and are calculated by 40 . Besides, the terms \(\int_{\mathbb{R}^3} r \; \mathrm{d}{\boldsymbol{v}}, \int_{\mathbb{R}^3} {\boldsymbol{v}}r \; \mathrm{d}{\boldsymbol{v}}\) and \(\int_{\mathbb{R}^3} |{\boldsymbol{v}}|^2 r \; \mathrm{d}{\boldsymbol{v}}\) keep constant in the inner iteration, and only one calculation is needed.

For now, the macroscopic moment system is completely derived. Next, we will introduce how it is adopted to accelerate inner iteration 56 . At each outer Newton iteration, the correction is initialized as \(g^{(0)} = 0\). Then, at the \(l\)-th inner iteration, the following three steps are applied

• Obtain the macroscopic variables \(\boldsymbol{m}^{(l)}=(\delta\rho^{(l)}/\rho, \delta\boldsymbol{u}^{(l)}/T,\delta T^{(l)}/T^2)\) related to \(g^{(l)}\) by 21 , the non-equilibrium moments \(\boldsymbol{\sigma}^{(\mathrm{neq},l)}, \boldsymbol{q}^{(\mathrm{neq},l)}\) by 64 , and \(\boldsymbol{\sigma}^{(\nabla,l)}, \boldsymbol{q}^{(\nabla,l)}\) by 40 .

• Obtain the intermediate macroscopic variables \[\label{eq:mid95macro} {m}^{(l+1,\ast)} = \left(\frac{\delta \rho^{(l+1, \ast)}}{ \rho}, \frac{\delta \boldsymbol{u}^{(l+1,\ast)}}{T}, \frac{\delta T^{(l+1,\ast)}}{ T^2}\right)^T\tag{65}\] by solving the closed macroscopic system 60 , once \(\boldsymbol{\sigma}^{(\mathrm{neq},l)}\), \(\boldsymbol{\sigma}^{(\nabla,l)}\) and \(\boldsymbol{q}^{(\mathrm{neq},l)}\), \(\boldsymbol{q}^{(\nabla,l)}\) are known. Then, reconstruct the local equilibrium part \(\mathcal{E}^{(l+1,\ast)}\) from \({m}^{(l+1,\ast)}\).

• Obtain the Newton correction \(g^{(l+1)}\) by the modified source iteration \[\label{eq:acc95SI} {\boldsymbol{v}}\cdot\nabla_{\boldsymbol{x}}g^{(l+1)} + \frac{{\nu}}{{\rm Kn}} \,g^{(l+1)} = \frac{1}{{\rm Kn}}\mathcal{L}\!\left[g^{(l)}\right] + \frac{{\nu}}{{\rm Kn}}\,g^{(l)} + \alpha^{(l)} \frac{{\nu}}{{\rm Kn}}\!\left(\mathcal{E}^{(l+1,\ast)} - \mathcal{E}^{(l)}\right) + r,\tag{66}\] where \(\alpha^{(l)} \in [0,1]\) is a relaxation parameter.

  Remark 1. With the macroscopic correction \[\label{eq:macro95corr} \alpha \frac{{\nu}}{{\rm Kn}}(\mathcal{E}^{(l+1,\ast)} - \mathcal{E}^{(l)}),\qquad{(1)}\] the information of the updated macroscopic variables is injected into the microscopic inner iteration, and therefore, the convergence of the macroscopic component in the correction \(g^{(l)}\) is accelerated, especially when \({\rm Kn}\) is small. A similar acceleration method can also be found in [24].

To distinguish from Newton-SI, we call this Newton method accelerated by the macroscopic synthetic system in the inner iteration as Newton-MS. Though the convergence of the inner iteration 66 can be accelerated with the macroscopic synthetic system 60 , directly solving 60 is still expensive. In the following section, the discrete form of 60 is directly deduced in the framework of the discontinuous Galerkin (DG) method, which is reduced into a linear system of \(g^{(n)}\) and can be directly solved.

4 Discrete form of Newton-MS↩︎

In this section, the fully discrete form of Newton-MS will be introduced, and the whole process can be treated as deducing the discrete form of the macroscopic synthetic system 60 in the framework of DG. We begin from the discretization of the spatial space where the DG method is adopted. The physical domain \(\Omega \subset \mathbb{R}^d\) (\(d=1,2\)) is partitioned into \(N_{\mathrm{el}}\) elements, and on each element, the numerical solution is approximated by polynomials of degree \(N_p\). For the microscopic velocity space, as in Sec. 2.2, the microscopic velocity space is first truncated to a finite domain as \([-L, L]^3\), and then the uniform mesh with number \(2N\) is utilized to discretize this microscopic velocity space. Thus, the discrete velocity set is \(\{\boldsymbol{v}_k\}_{k=1}^{(2N)^3}\), where \(\boldsymbol{v}_k=(v_{1k},v_{2k},v_{3k})^T\), and the corresponding quadrature weight is \(\Delta v=(L/N)^3\). For simplicity, we introduce the total degrees of freedom \(N_{t}\), those in the spatial space \(N_{s}\) and those in the microscopic space \(N_{v}\) as \[\label{eq:dof} N_t = N_s \times N_v, \qquad N_{s} = N_{\mathrm{el}}(N_p+1)^d, \qquad N_{v} = (2N)^3.\tag{67}\] With this discretization, let \(\boldsymbol{f}^{(n)}, \boldsymbol{g}^{(n)} \in \mathbb{R}^{N_t}\) be the global numerical solution and the global discrete correction, and their entries are ordered first by the spatial index. In particular, \[\boldsymbol{f}^{(n)} =\left((\boldsymbol{f}_1^{(n)})^T, \cdots, (\boldsymbol{f}_{N_{s}}^{(n)})^T \right)^T,\quad \boldsymbol{g}^{(n)} =\left((\boldsymbol{g}_1^{(n)})^T, \cdots, (\boldsymbol{g}_{N_{s}}^{(n)})^T \right)^T,\qquad \boldsymbol{f}_i^{(n)}, \boldsymbol{g}_i^{(n)}\in\mathbb{R}^{N_v},\] where \(\boldsymbol{f}_i^{(n)}, \boldsymbol{g}_i^{(n)}\) represent the distribution function \(\boldsymbol{f}\), and the correction \(\boldsymbol{g}\) at the \(i\)-th mesh, respectively. Thus, the update in each Newton step 10 becomes \[\label{eq:dis95Newton95update} \boldsymbol{f}^{(n+1)} = \boldsymbol{f}^{(n)} - \boldsymbol{g}^{(n)}.\tag{68}\] Let \(\mathbf{T}\) denote the DG discretization of the transport operator \({\boldsymbol{v}}\cdot\nabla_{\boldsymbol{x}}\). With the assumption of the upwind numerical flux with a homogeneous inflow condition, the operator \(\mathbf{T}\) is reduced to a matrix as \(\mathbf{T}\in\mathbb{R}^{N_{t}\times N_{t}}\). Moreover, it is block-tridiagonal in 1D and block-pentadiagonal in 2D on a structured mesh, with each block of size \((N_v (N_p+1)^d) \times (N_v (N_p+1)^d)\), due to the coupling with adjacent spatial elements through the upwind flux. In this case, the physical boundary condition is incorporated explicitly through a boundary source term. Thus, the discrete form 13 has the form as \[\label{eq:fully95discrete} \mathbf{T}\boldsymbol{g}^{(n)} - \frac{1}{{\rm Kn}}\mathbf{L}^{(n)}\boldsymbol{g}^{(n)} = \boldsymbol{r}^{(n)}+\boldsymbol{b}_{\rm bc}^{(n)},\tag{69}\] where \(\mathbf{L}^{(n)} \in \mathbb{R}^{N_t \times N_t}\) is the discrete form of the linearized collision operator 14 , \(\boldsymbol{r}^{(n)}\) is the discrete outer residual, and \(\boldsymbol{b}_{\rm bc}^{(n)}\) is the boundary source term, depending on \(\boldsymbol{g}^{(n)}\). Since the linear collision operator \(\mathcal{L}^{(n)}\) is a local operator, \(\mathbf{L}^{(n)}\) is block-diagonal \[\label{eq:diag95L} \mathbf{L}^{(n)} = \operatorname{diag}\left\{ \mathbf{L}_1^{(n)},\, \mathbf{L}_2^{(n)},\, \ldots,\, \mathbf{L}_{N_{s}}^{(n)} \right\},\tag{70}\] with each block \(\mathbf{L}_i^{(n)}\in\mathbb{R}^{N_v\times N_v}\). In the following, the outer iteration superscript “\((n)\)” will be omitted for brevity.

To deduce the discrete form of the macroscopic synthetic system 60 , the discrete form of the deduction process as in Sec. 3 is proposed. Precisely, with the Chapman-Enskog ansatz 25 , the discrete correction \(\boldsymbol{g}\) is expanded as \[\label{eq:dis95ce} \boldsymbol{g} = \boldsymbol{\mathcal{E}} + {\rm Kn}\tilde{\boldsymbol{g}}_1 + {\rm Kn}^2\tilde{\boldsymbol{g}}_2, \qquad \tilde{\boldsymbol{g}}_1, \tilde{\boldsymbol{g}}_2 \in \mathbb{R}^{N_t},\tag{71}\] where \(\boldsymbol{\mathcal{E}}\) is the discrete local equilibrium of \(\boldsymbol{g}\), \(\tilde{\boldsymbol{g}}_1\) is the discrete first-order Navier–Stokes correction, and \(\tilde{\boldsymbol{g}}_2\) represents the discrete high-order non-equilibrium remainder. Substituting 71 into the discrete Newton correction equation 69 , and using \(\mathbf{L} \boldsymbol{\mathcal{E}}=0\), we obtain \[\label{eq:dis95CE95eq} \mathbf{T} \boldsymbol{\mathcal{E}} + {\rm Kn}\, \mathbf{T} \tilde{\boldsymbol{g}}_1 + {\rm Kn}^2 \mathbf{T} \tilde{\boldsymbol{g}}_2 - \mathbf{L}\tilde{\boldsymbol{g}}_1 - {\rm Kn}\, \mathbf{L}\tilde{\boldsymbol{g}}_2 = \boldsymbol{r}+\boldsymbol{b}_{\rm bc}.\tag{72}\] To obtain the discrete macroscopic equation related to 31 , the discrete collision invariant operator \(\Phi \in \mathbb{R}^{5N_{s}\times N_{t}}\) is introduced as \[\label{eq:dis95V} \Phi = \operatorname{diag}\{\bar{\Phi}, \ldots, \bar{\Phi}\}, \qquad \bar{\Phi} = \Delta v \left( \phi({\boldsymbol{v}}_1), \cdots, \phi({\boldsymbol{v}}_{N_v}) \right) \in \mathbb{R}^{5\times N_v},\tag{73}\] where \(\phi({\boldsymbol{v}})\) is the collision invariants 22 . Multiplying 72 by \(\Phi\) from the left, and with the discrete conservation property \[\label{eq:dis95conver} \Phi\mathbf{L}=0,\tag{74}\] we can derive the discrete form of 31 as \[\label{eq:dis95mac95eq1} \Phi \mathbf{T}\left( \boldsymbol{\mathcal{E}}+ {\rm Kn}\tilde{\boldsymbol{g}}_1 + {\rm Kn}^2\tilde{\boldsymbol{g}}_2\right) = \Phi\boldsymbol{r} + \Phi\boldsymbol{b}_{\rm bc}.\tag{75}\] Then, a similar closure is displayed to obtain the closed discrete macroscopic synthetic system. We first introduce two block-diagonal operators \(\mathbf{S}\) and \(\Gamma\) as \[\label{eq:S95G} \mathbf{S} = \operatorname{diag}\{\mathbf{S}_1,\ldots,\mathbf{S}_{N_s}\}, \qquad \Gamma = \mathrm{diag}\{\Gamma_1,\ldots,\Gamma_{N_s}\},\tag{76}\] with \(\mathbf{S}_i\) and \(\Gamma_i\) as \[\begin{align} \label{eq:SG} \mathbf{S}_i &= \Delta v\left(s_i({\boldsymbol{v}}_1),\cdots, s_i({\boldsymbol{v}}_{N_v})\right)\in \mathbb{R}^{5\times N_v}, \qquad s_i({\boldsymbol{v}}_k) = \overline{\varphi}(\boldsymbol{c}_{ki}), \\ \Gamma_i &= \left(\gamma_i({\boldsymbol{v}}_1), \cdots, \gamma_i({\boldsymbol{v}}_{N_v})\right)^T \in \mathbb{R}^{N_v\times 5}, \qquad \gamma_i({\boldsymbol{v}}_k) = \mathcal{M}_i({\boldsymbol{v}}_k)\varphi(\boldsymbol{c}_{ki}). \end{align}\tag{77}\] Here, \(\boldsymbol{c}_{ki} = {\boldsymbol{v}}_k - \boldsymbol{u}_i\) is the peculiar velocity, \(\overline{\varphi}\) and \(\varphi\) are defined in 21 , and 18 , respectively. \(\mathcal{M}_i\) is the Maxwellian 4 , whose corresponding macroscopic variables are \((\rho_i, \boldsymbol{u}_i, T_i)\) at the \(i\)-th mesh. It is easy to derive that \[\label{eq:oper95S} m_i = \mathbf{S}_i \boldsymbol{g}_i, \qquad \boldsymbol{\mathcal{E}}_i = \Gamma_i m_i, \qquad i = 1, \cdots, N_s,\tag{78}\] It means that the operator \(\mathbf{S}\) maps the discrete correction \(\boldsymbol{g}\) to the local macroscopic variables \(m_i\) 20 , and \(\Gamma\) reconstructs local equilibrium from the local macroscopic variables. Moreover, we can deduce that \(\Gamma\mathbf{S}\) is equal to the discrete form of the orthogonal projection operator \(\mathcal{P}\) 23 . With the definition 57 , it holds that \[\label{eq:dis95g1} \tilde{\boldsymbol{g}}_1^{(l)} = \tilde{\boldsymbol{g}}_1^{(\nabla, l)} - \tilde{\boldsymbol{g}}_1^{(r)}, \qquad \tilde{\boldsymbol{g}}_1^{(\nabla, l)}=\mathbf{L}^{-1} \mathbf{T}\Gamma\boldsymbol{m}^{(l)},\tag{79}\] where \(\boldsymbol{m} = (m_1^T, \cdots, m_{N_s}^T)^T\in \mathbb{R}^{5N_s}\), and \(\mathbf{L}^{-1}\) is the discrete form of the operator \(\mathcal{L}^{-1}\). In this case, with the discrete Chapman-Enskog expansion 58 , 71 and 78 , the high order \(\tilde{\boldsymbol{g}}_2^{(l)}\) is obtained as \[\label{eq:dis95g2} {\rm Kn}^2 \tilde{\boldsymbol{g}}_2^{(l)} = \boldsymbol{g}^{(l)} - \Gamma \boldsymbol{m}^{(l)} - {\rm Kn}\tilde{\boldsymbol{g}}_1^{(\nabla,l)}+ {\rm Kn}\tilde{\boldsymbol{g}}_1^{(r)} \triangleq \hat{\boldsymbol{g}}^{(l)} +{\rm Kn}\tilde{\boldsymbol{g}}_1^{(r)}.\tag{80}\] Substituting 78 , 79 and 80 into the discrete moment equation 75 , the residual-driven term \(\boldsymbol{g}_1^{(r)}\) is canceled as in the continuous case 60 . Then, the discrete macroscopic system for the local macroscopic variables \(\boldsymbol{m}^{(l+1,\ast)}\) in the inner iteration is driven as \[\label{eq:dis95mac95eq2} \Phi \mathbf{T}\left(\Gamma \boldsymbol{m}^{(l+1,\ast)} + {\rm Kn}\,\tilde{\boldsymbol{g}}_1^{(\nabla, l+1,\ast)} + \hat{\boldsymbol{g}}^{(l)}\right)=\Phi\boldsymbol{r}+ \Phi \boldsymbol{b}_{\rm bc}^{(l)}.\tag{81}\] The macroscopic system 81 is the discrete form of the continuous closed system 60 and can be rewritten as \[\label{eq:dis95mac95eq3} \left(\Psi_{\mathrm{E}}+{\rm Kn}\,\Psi_{\mathrm{NS}}\right)\boldsymbol{m}^{(l+1,\ast)}=\boldsymbol{b}^{(l)},\tag{82}\] where \[\label{eq:matrix95K} \Psi_{\mathrm{E}}=\Phi\mathbf{T}\Gamma,\qquad\Psi_{\mathrm{NS}}=\Phi\mathbf{T} \mathbf{L}^{-1} \mathbf{T}\Gamma, \qquad \boldsymbol{b}^{(l)}=-\Phi\mathbf{T}\hat{\boldsymbol{g}}^{(l)}+\Phi\boldsymbol{r}+\Phi\boldsymbol{b}_{\rm bc}^{(l)}.\tag{83}\] Here, \(\Psi_{\mathrm{E}}\) represents the contribution of the equilibrium as 35 , while \(\Psi_{\mathrm{NS}}\) corresponds to the discrete first-order Chapman-Enskog correction, or the terms \((\cdot)^{(\nabla)}\) on the left side of 60 .

Note that \(\Psi_{\rm E}\) and \(\Psi_{\rm NS}\) depend only on the outer iteration variables and remain constant in the inner iterations and their reconstruction is based on applications of the sparse transport operator \(\mathbf{T}\) on the block-diagonal matrix \(\Gamma\), and the applications of the block-diagonal matrix \(\Phi\), \(\mathbf{S}\) and \(\Gamma\).

Finally, the discrete form of the augmented source iteration 66 is \[\label{eq:dis95inner95si} \left(\mathbf{T} + \frac{\boldsymbol{\nu}}{{\rm Kn}}\right)\boldsymbol{g}^{(l+1)} = \frac{1}{{\rm Kn}}\mathbf{L}\boldsymbol{g}^{(l)} + \frac{\boldsymbol{\nu}}{{\rm Kn}}\boldsymbol{g}^{(l)} + \alpha \frac{\boldsymbol{\nu}}{{\rm Kn}}\Gamma\bigl(\boldsymbol{m}^{(l+1,\ast)}-\boldsymbol{m}^{(l)}\bigr) + \boldsymbol{r} + \boldsymbol{b}_{\rm bc}^{(l)}.\tag{84}\]

For the computational complexity, the one-time setup within each Newton step consists of assembling the synthetic macroscopic operators \(\Psi_{\mathrm{E}}\) and \(\Psi_{\mathrm{NS}}\), whose construction is approximated as \(\mathcal{O}(N_s N_p^dN^3)\). For one inner iteration, the dominant cost comes from the implement of the linearized collision operator \(\mathbf{L}\boldsymbol{g}^{(l)}\), whose complexity is \(\mathcal{O}(N_sN^4\log N)\) with the FFT-based algorithm in Sec. 2.2. The computational cost brought by the transport operator \(\mathbf{T}\) is approximated as \(\mathcal{O}(N_sN_p^{d}N^3)\), while the numerical cost of solving the macroscopic synthetic system is \(\mathcal{O}(N_s)\). Therefore, compared with Newton-SI, Newton-MS has essentially the same leading-order cost in the microscopic iteration. For the completeness of Newton-MS, the algorithm is listed in Alg. 1.

5 Numerical results↩︎

In this section, several numerical experiments are presented to validate the accuracy and efficiency of the proposed Newton-MS scheme. The main comparison is made with the standard Newton source iteration (Newton-SI), focusing on convergence behavior and computational cost over a range of Knudsen numbers, especially for the small Knudsen number. The collision operator is evaluated by the fast Fourier spectral method described in Sec. 2.2, and the transport term is discretized by the discontinuous Galerkin (DG) method described in Sec. 4.

The stop criteria for the outer Newton iteration 10 is that the residual of the nonlinear steady Boltzmann equation satisfies \[\label{eq:residual} R_{\mathrm{out}}(n) = \left(\int_{\Omega}\int_{[-L,L]^3}\left|r^{(n)} \right|^2 \; \mathrm{d}{\boldsymbol{v}}\; \mathrm{d}\boldsymbol{x}\right)^{1/2} < \epsilon_{\rm out}.\tag{85}\] Moreover, the stop criteria for the inner iteration 56 are \[\label{eq:resi95inner} R_{\mathrm{in}}(n,l)=\left(\int_{\Omega}\int_{[-L,L]^3}\left| \mathcal{R}_{\mathrm{in}}^{(n,l)}\right|^2\; \mathrm{d}{\boldsymbol{v}}\; \mathrm{d}\boldsymbol{x}\right)^{1/2}< \epsilon_{\rm in,1}, \quad\text{or}\quad \frac{R_{\mathrm{in}}(n,l)}{R_{\mathrm{out}}(n)}<\epsilon_{\rm in, 2},\tag{86}\] where the residual of the inner iteration \(\mathcal{R}_{\mathrm{in}}^{(n,l)}\) is defined as \[\label{eq:inner95residual95def} \mathcal{R}_{\mathrm{in}}^{(n,l)} = {\boldsymbol{v}}\cdot \nabla_{\boldsymbol{x}} g^{(n,l)} - \frac{1}{{\rm Kn}}\mathcal{L}^{(n)}[g^{(n,l)}] - r^{(n)}.\tag{87}\] Here, the relative error in the stop criterion 86 prevents unnecessary over-iteration of the linearized correction equation for \(g\) when a highly accurate inner solution is not required for the outer Newton iteration. In the simulation, the related parameters \(\epsilon_{\rm out}\), \(\epsilon_{\rm in, 1}\) and \(\epsilon_{\rm in, 2}\) are problem dependent, and unless otherwise specified, they are set as 88 in this work. \[\label{eq:ep} \epsilon_{\rm out} = 10^{-5}, \qquad \epsilon_{\rm in, 1} = 10^{-6}, \qquad \epsilon_{\rm in, 2} = 10^{-2}.\tag{88}\]

Moreover, the macroscopic correction \[\label{eq:macro95cor} \alpha^{(l)} \frac{\nu}{{\rm Kn}}\left(\mathcal{E}^{(l+1, \ast)} - \mathcal{E}^{(l)}\right)\tag{89}\] in 66 mainly speeds up the convergence of macroscopic variables, including density, macroscopic velocity, and temperature. Therefore, the inner residual drops quickly in the first few iterations, but it may stop decreasing once these macroscopic variables have converged, and the residual mainly comes from the non-equilibrium parts. In this case, the macroscopic correction 89 is no longer needed to achieve convergence of the correction \(g\). Thus, the relaxation parameter \(\alpha^{(l)}\) in 89 is chosen as \[\alpha^{(n,l)} = \begin{cases} \alpha_0, & l < l_{\mathrm{sw}}^{(n)},\\ 0, & l \geqslant l_{\mathrm{sw}}^{(n)}, \end{cases} \label{eq:relax95para}\tag{90}\] where \(l_{\mathrm{sw}}^{(n)}\) is the switching index at the \(n\)-th Newton step, which is determined by the following residual-decay criterion \[l_{\mathrm{sw}}^{(n)} = \min\left\{ l\geqslant p+1: \frac{R_{\mathrm{in}}(n,l)}{R_{\mathrm{in}}(n,l-j)} > \eta_{\mathrm{sw}}, \quad j=1,\ldots,p \right\}. \label{eq:switch95criterion}\tag{91}\] Here, \(p\) chosen as \(p = 3\) is the monitoring window length and \(\eta_{\mathrm{sw}}=0.9\) is a threshold to measure the rate of convergence. In the following tests, \(\alpha_0=0.4\) and \(\alpha_0 = 1\) are utilized for the one-dimensional Fourier and Couette flow problems, and for the two-dimensional cavity-flow problems, respectively.

For the truncation region of the microscopic velocity space, \(L\) is set as \[\label{eq:L} L=\frac{3+\sqrt{2}}{2}R,\qquad R=4.\tag{92}\] Besides, when the fast Fourier spectral method is utilized to solve the linearized collision term, the total mass conservation can not be conserved [27]. Therefore, a post-processing is added to keep mass conservation. Precisely, after each outer Newton iteration, the distribution function is rescaled by \[\label{eq:mass95correction} f^{(n+1)} \leftarrow f^{(n+1)}\cdot \frac{m_{\text{tot}}^{(0)}}{m_{\text{tot}}^{(n+1)}}, \qquad m_{\text{tot}}^{(n+1)} = \int_\Omega\int_{[-L,L]^3} f^{(n+1)}\,\mathrm{d}{\boldsymbol{v}}\,\mathrm{d}\boldsymbol{x},\tag{93}\] where \(m_{\text{tot}}^{(0)}\) is the initial total mass.

In the simulation, the efficiency comparison between Newton-MS, Newton-SI and GSIS is displayed. For all three methods, they are all implemented with the same spatial and velocity discretization. For Newton-MS and Newton-SI, the outer iteration is the same nonlinear Newton iteration 10 , but the macroscopic moment equations are utilized to accelerate the convergence 66 of the inner iteration in Newton-MS, while 15 is directly adopted in Newton-SI. For GSIS, only the outer iteration with macroscopic synthetic acceleration is adopted, and no inner iteration is needed.

For all three methods, to obtain the nonlinear outer residual \(r^{(n)}\) 12 , although the quadratic collision term only needs to be calculated once in each outer iteration, this is still quite expensive, and the CPU time to obtain \(r^{(n)}\) for all \((n)\) is labeled \(T_{\rm out}\). For Newton-MS and Newton-SI, to obtain the linear inner residual \(\mathcal{R}_{\rm in}^{(n,l)}\) 87 , only the linear collision term is calculated, but it needs to be calculated once for each inner iteration, which will also be quite expensive if the number of inner iterations is large. Here, the CPU time to obtain \(\mathcal{R}_{\rm in}^{(n,l)}\) for all \((n,l)\) is labeled \(T_{\rm in}\). Moreover, the CPU time for solving the macroscopic synthetic equations for all \((l,n)\) is labeled \(T_{\rm m}\) with \(T_{\rm tol}\) the total CPU time.

5.1 1D Fourier flow problem↩︎

We first consider the classical one-dimensional Fourier flow problem. The scenario consists of two parallel, stationary plates at \(x=0\) and \(x=1\), with wall temperatures \(T_L=1\) and \(T_R=1.2\), respectively. The discretization of the spatial domain is \(N_{\rm el} = 40\) with the degree of polynomial \(N_p = 2\). For the discretization of the microscopic velocity space, the number of the Fourier modes is set as \(N = 24\), corresponding to the degree of freedom \(2N=48\) in each direction of the microscopic velocity space. The fully diffuse reflection boundary condition [27] is imposed on both walls. The similar Fourier flow problem is also studied in [31].

The numerical solution of the density \(\rho\), temperature \(T\) at the steady state for \({\rm Kn}= 1, 0.1, 0.01\) and \(0.001\) by Newton-MS is shown in Fig. 2, where the reference solution for \({\rm Kn}= 0.001\) is obtained by directly solving the Navier-Stokes equation with temperature slip boundary condition imposed and those for other \({\rm Kn}\) are obtained by Newton-SI. Fig. 2 indicates that for both \(\rho\) and \(T\), the numerical solution matches well with the reference solution for all Knudsen numbers.

The efficiency comparison of Newton-MS, Newton-SI, and GSIS is summarized in Tab. 1. It indicates that for Newton-MS, the inner iteration number \(N_{\rm in}\) remains small with decreasing \({\rm Kn}\), while for Newton-SI, \(N_{\rm in}\) is increasing rapidly, and it even fails to converge when \({\rm Kn}= 0.001\). The behavior of \(T_{\rm in}\) is similar, which is at the same order for \({\rm Kn}= 1\), but that of Newton-SI is \(11\) times that of Newton-MS for \({\rm Kn}= 0.01\). Moreover, the total computational time \(T_{\rm tol}\) of Newton-MS is also greatly reduced to more than \(7\) times compared to Newton-SI for \({\rm Kn}= 0.01\). We can expect that efficiency can be further improved for smaller Knudsen numbers. Compared to GSIS, the outer iteration number is much smaller for Newton-MS. Thus, the total CPU time for Newton-MS is only half that of GSIS for \({\rm Kn}= 1\) and \(0.001\), and about three-quarters for \({\rm Kn}= 0.1\) and \(0.01\). This means that though the inner iteration is added, the total computational cost can still be reduced for Newton-MS compared to GSIS for all \({\rm Kn}\). Besides, Tab. 1 also shows that for all Knudsen numbers, the computational cost for solving the macroscopic synthetic system is quite small, all less than \(5\%\) of the total computational cost.

The convergence histories of the outer Newton and inner iterations for \({\rm Kn}= 0.01\) are plotted in Fig. 3. It shows that the evolution of the outer Newton iteration residual \(r^{(n)}\) for Newton-MS and Newton-SI is almost the same as in Fig. 3 (a), indicating that the macroscopic synthetic system does not deteriorate the convergence of the outer Newton iteration. Fig. 3 (b) to 3 (d) present the evolution of the inner iteration residual \(\mathcal{R}_{\rm in}^{(n,l)}\) for each outer Newton iteration of both methods. It shows that compared to Newton-SI, the inner iteration residual of Newton-MS decreases much more quickly for each outer Newton iteration, which is also consistent with the results in Tab. 1.

5.2 1D Couette flow problem↩︎

In this section, the planner Couette flow problem is studied to validate the performance of Newton-MS in the shear-dominated problems. The scenario consists of two parallel plates at \(x=0\) and \(x=1\), which move in opposite tangential directions with wall velocities \(u_w=\pm 0.5\) and fixed temperature \(T_w=1.0\). Fully diffuse reflection boundary conditions are imposed on both walls. The similar example is also tested in [31]. The discretization of the spatial and microscopic velocity space is the same as that in Sec. 5.1.

The steady state numerical solution of the density \(\rho\), the tangential velocity \(u_y\) and the temperature for \({\rm Kn}= 1, 0.1, 0.01\) and \(0.001\) is shown in Fig. 4, where the reference solution is also plotted. Here, the reference solution of \({\rm Kn}= 0.001\) is obtained by the compressible Navier-Stokes equations with slip boundary conditions, while the others are obtained by Newton-SI. Fig. 4 shows that for all Knudsen numbers, the numerical solution agrees well with the reference solution.

The comparison of the computational efficiency between Newton-MS, Newton-SI, and GSIS is summarized in Tab. 2, where the similar behavior is found as in Sec. 5.1. For \({\rm Kn}=1\), Newton-MS has a comparable cost to Newton-SI, since the source iteration is already efficient in the rarefied regime and the additional macroscopic correction brings only limited benefit. As \({\rm Kn}\) decreases, however, the advantage of Newton-MS becomes increasingly evident. Newton-SI suffers from the slow convergence of the inner source iteration, whereas Newton-MS keeps the number of inner iterations small. Compared with GSIS, Newton-MS also requires fewer outer iterations and achieves a lower total CPU time, especially in the near-continuum regime. The computational cost \(T_{\rm m}\) to solve the macroscopic synthetic system all remains a small fraction of the total cost for all Knudsen numbers.

The convergence histories of the outer Newton and inner iterations for \({\rm Kn}= 0.01\) are plotted in Fig. 5, the behavior of which is similar to that in Sec. 5.1. Precisely, the evolution of the outer Newton iteration residual \(r^{(n)}\) for Newton-MS and Newton-SI is almost the same as in Fig. 5 (a), indicating that the macroscopic synthetic system does not deteriorate the convergence of the outer Newton iteration. Fig. 5 (b) to 5 (c) present the evolution of the inner iteration residual \(\mathcal{R}_{\rm in}^{(n,l)}\) for each outer Newton iteration of both methods. It shows that compared to Newton-SI, the inner iteration residual of Newton-MS decreases much more quickly for each outer Newton iteration, which is also consistent with the results in Tab. 2. The inner iteration residual reaches the tolerance much faster for Newton-MS in each outer Newton iteration, compared to Newton-SI, indicating the effect of the acceleration for the macroscopic synthetic system when applied to the shear-dominated problems.

5.3 2D lid-driven cavity flow problem↩︎

In this section, the two-dimensional lid-driven cavity flow problem is studied to demonstrate the capability of Newton-MS for the high dimensional problems. The scenario is a unit square domain \((x, y) \in [0,1]^2\). The top wall at \(y = 1\) is moving with a constant tangential velocity \(\boldsymbol{u}_w = (0.5, 0, 0)\), while the other three walls are stationary. All walls are maintained at the fixed temperature \(T_w=1.0\), and fully diffuse reflection boundary conditions are imposed. Similar problems are also studied in [31], [32]. The mesh size in the spatial space is \(N_{\rm el} = 40 \times 40\) with the DG polynomial degree set to \(N_p = 2\). For the microscopic velocity space, the number of Fourier modes is chosen as \(N = 16\).

The steady state numerical solution of the density \(\rho\), temperature \(T\), the macroscopic velocity \(u_x\) and \(u_y\) for \({\rm Kn}= 0.01\) is shown in Fig. 6, where the reference solution obtained by Newton-SI is also plotted. It shows that the numerical solution matches well with the reference solution. The computational efficiency of Newton-MS and Newton-SI is summarized in Tab. 3. It shows that though the number of outer Newton iteration is the same for both methods, the inner iteration number of Newton-MS is much smaller for each outer Newton iteration. Therefore, the total computational time of Newton-MS is less than half of Newton-SI. Moreover, the additional cost of solving the macroscopic synthetic system is quite small, less than \(5\%\) of the total cost.

The convergence histories of the outer Newton and inner iterations for Kn = 0.01 are plotted in Fig. 7. The evolution of the outer Newton residuals for Newton-MS and Newton-SI is similar, as shown in Fig. 7 (a). The evolution of \(\mathcal{R}_{\rm in}^{(n,l)}\) for each Newton iteration are presented in Fig. 7 (b) to 7 (d), respectively. They illustrate that the decay of \(\mathcal{R}_{\rm in}^{(n,l)}\) is much faster for Newton-MS compared to Newton-SI, which is also consistent with the results in Tab. 3.

5.4 2D thermal cavity flow problem↩︎

In this section, the 2D thermal cavity flow problem is studied. The scenario is also a unit square domain \((x,y)\in[0,1]^2\) as in Sec. 5.3, but all four walls are stationary. The temperature of the top wall is set as \(T_{\mathrm{top}}=1.2\), while all others are fixed as \(T = 1\). The fully diffuse reflection boundary conditions are imposed on all walls. Different from the lid-driven cavity flow in Sec. 5.3, there is no external momentum added, and the flow is driven entirely by the difference of the wall temperature. The similar example is also studied in [31].

The same discretization of the spatial space and microscopic velocity as in Sec. 5.3 is adopted. The steady state numerical solution of the density \(\rho\), and the temperature \(T\) for \({\rm Kn}= 0.01\) is plotted in Fig. 8, where it matches well with the reference solution obtained by Newton-SI. The computational efficiency is summarized in Tab. 4, where the behavior of both methods is similar as that in Sec. 5.3. The total computational time of Newton-MS is reduced to one third of that by Newton-SI, which is also mainly due to the reduced number of the inner iterations. Moreover, the additional cost of solving the macroscopic synthetic system still remains quite small, compared to the total computational time saved.

The convergence histories of the outer Newton and inner iterations for \({\rm Kn}= 0.01\) are plotted in Fig. 9. The outer Newton residuals decay at almost the same rates for Newton-SI and Newton-MS, as shown in Fig. 9 (a), indicating that the synthetic acceleration does not affect the outer nonlinear convergence. For the inner iteration, the evolution of \(\mathcal{R}_{\rm in}^{(n,l)}\) is the similar to that in Sec. 5.3, which all shows a much more rapid decay rate for Newton-MS compared to that of Newton-SI.

6 Conclusion↩︎

In this work, a modified Newton’s method accelerated by a macroscopic synthetic system (Newton-MS) is proposed for the steady-state Boltzmann equation. For Newton-MS, the outer iteration is the normal Newton method, while the linearized collision operator is utilized instead of the quadratic Boltzmann collision operator in the inner iteration to reduce computational cost. For the inner iteration, a macroscopic synthetic system based on a Chapman-Enskog closure is first derived, and then the information of the updated macroscopic variables obtained by solving this macroscopic synthetic system is added in the general source iteration to accelerate inner iteration convergence. The discrete matrix form of the macroscopic synthetic system is derived in the framework of the discontinuous Galerkin method for the practical implementation, where the computational cost can be greatly reduced compared to directly discretizing the continuous macroscopic synthetic system. Numerical results indicate that Newton-MS maintains consistent convergence properties from the free-molecular regime to the continuum limit and offers a practical and effective alternative to standard Newton-SI schemes for rarefied gas simulations.

Acknowledgments↩︎

We thank Prof. Lei Wu from Southern University of Science and Technology for providing the GSIS code. This work of Yanli Wang is partially supported by the Science Challenge program (NO. TZ2025016). This work of Pei Zhang was supported by the China Scholarship Council (NO. 202404890001). The work of Zhenning Cai was supported by the Academic Research Fund of the Ministry of Education of Singapore under grant A-8002392-00-00.

7 Derivation of \(\sigma^{(\nabla)}\) and \(q^{(\nabla)}\)↩︎

In this appendix, the detailed deduction of the stress tensor \(\sigma^{(\nabla)}\) and heat flux \(q^{(\nabla)}\) in 38 is presented. The derivation is based on the expansion of \({\boldsymbol{v}}\cdot\nabla_{\boldsymbol{x}}\mathcal{E}\) on the microscopic velocity space, and the rotational invariance of the linear collision operator \(\mathcal{L}\).

7.1 Deducing \(\sigma^{(\nabla)}\) and \(q^{(\nabla)}\)↩︎

We first introduce some abbreviations for convenience as \[\tag{94} \begin{gather} \tag{95} \partial_i = \frac{\partial {}}{\partial {x_i}}, \qquad A_i = \frac{\partial_i\rho}{\rho}, \qquad U_{ij} = \partial_i \boldsymbol{u}_{j}, \qquad G_i = \partial_i T, \qquad D_{ij} = \partial_i\delta \boldsymbol{u}_j, \\ \mathcal{A}_{ij}(\boldsymbol{c}) = c_i c_j - \frac{|\boldsymbol{c}|^2}{3}\delta_{ij}, \qquad \mathcal{B}_i(\boldsymbol{c}) = \frac{1}{2}|\boldsymbol{c}|^2 c_i, \qquad \alpha = \frac{\delta\rho}{\rho}, \qquad \beta_i = \frac{\delta u_i}{T}, \qquad \vartheta = \frac{\delta T}{2T^2}. \tag{96} \end{gather}\] Then the local equilibrium 19 can be rewritten as \[\label{eq:H95def} \mathcal{E}(\boldsymbol{x},{\boldsymbol{v}}) = \mathcal{M}\,\chi, \qquad \chi = \alpha + \beta_i c_i + \vartheta (|\boldsymbol{c}|^2 - 3T).\tag{97}\] Here, the Einstein summation convention is utilized. Thus, it holds that \[{\boldsymbol{v}}\cdot\nabla_{\boldsymbol{x}}\mathcal{E} = \mathcal{M}\,{\boldsymbol{v}}\cdot\nabla_{\boldsymbol{x}}\chi + \chi\,{\boldsymbol{v}}\cdot\nabla_{\boldsymbol{x}}\mathcal{M} \triangleq \mathcal{M}\bigl(\mathcal{T}_\rho + \mathcal{T}_u + \mathcal{T}_T\bigr), \label{eq:vgradH95decomp}\tag{98}\] with \[\tag{99} \begin{align} \tag{100} \mathcal{T}_\rho &= \beta_k A_\ell\,\mathcal{A}_{k\ell}(\boldsymbol{c}) + \vartheta A_k\,|\boldsymbol{c}|^2 c_k + \mathcal{R}_{\rho},\\ \tag{101} \mathcal{T}_u &= S^{(u)}_{k\ell}\,\mathcal{A}_{k\ell}(\boldsymbol{c}) + R^{(u)}_{k\ell}\,|\boldsymbol{c}|^2\mathcal{A}_{k\ell}(\boldsymbol{c}) + C^{(u,3)}_k\,|\boldsymbol{c}|^2 c_k + \mathcal{R}_{u},\\ \tag{102} \mathcal{T}_T &= S^{(T)}_{k\ell}\,\mathcal{A}_{k\ell}(\boldsymbol{c}) + R^{(T)}_{k\ell}\,|\boldsymbol{c}|^2\mathcal{A}_{k\ell}(\boldsymbol{c}) + C^{(T,3)}_k\,|\boldsymbol{c}|^2 c_k + C^{(T,5)}_k\,|\boldsymbol{c}|^4 c_k + \mathcal{R}_{T}, \end{align}\] and \[\begin{align} &S^{(u)}_{k\ell} = \left(\frac{\alpha}{T} - 3\vartheta\right)U_{k\ell} + \frac{D_{k\ell}}{T} - \frac{\delta u_\ell\,G_k}{T^2} + \frac{\beta_k\,u_{m}U_{m\ell}}{T}, \qquad R^{(u)}_{k\ell} = \frac{\vartheta}{T}\,U_{k\ell}, \tag{103}\\ & C^{(u,3)}_k = \frac{\vartheta}{T}\,u_{\ell}U_{\ell k} + \frac{1}{5T}\bigl( b_\ell U_{\ell k} + b_\ell U_{k\ell} + \beta_k U_{\ell\ell} \bigr), \tag{104}\\ &S^{(T)}_{k\ell} = -\frac{3b_k G_\ell}{2T} - \frac{\delta T\,U_{k\ell}}{T^2}, \qquad R^{(T)}_{k\ell} = \frac{\beta_k G_\ell}{2T^2}, \tag{105}\\ &C^{(T,3)}_k = \frac{\alpha\,G_k}{2T^2} - \frac{3\vartheta\,G_k}{T} + \frac{\partial_k\delta T}{2T^2} - \frac{\delta T\,G_k}{T^3} + \frac{\beta_k\,u_{\ell}G_\ell}{2T^2}, \qquad C^{(T,5)}_k = \frac{\vartheta\,G_k}{2T^2}. \tag{106} \end{align}\] Here, only the terms that will contribute when calculating \(\sigma^{(\nabla)}\) and \(q^{(\nabla)}\) are listed, namely \[\bigl\{\, \mathcal{A}_{k\ell}(\boldsymbol{c}),\quad |\boldsymbol{c}|^2\,\mathcal{A}_{k\ell}(\boldsymbol{c}),\quad |\boldsymbol{c}|^2\,c_k,\quad |\boldsymbol{c}|^4\,c_k \,\bigr\}, \label{eq:relevant95modes}\tag{107}\] while all the others are summarized in \(\mathcal{R}_{s},s = \rho, u, T\). Moreover, due to the isotropy of the linearized collision operator \(\mathcal{L}\) [30], we can deduce that the integrals have the form below \[\tag{108} \begin{align} &\int_{\mathbb{R}^3}\mathcal{A}_{ij}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{A}_{k\ell}(\boldsymbol{c})\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} =: \frac{C_{2,0}}{\nu} \!\left(\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk} -\tfrac{2}{3}\delta_{ij}\delta_{k\ell}\right), \tag{109} \\ &\int_{\mathbb{R}^3}\mathcal{A}_{ij}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[|\boldsymbol{c}|^2\mathcal{A}_{k\ell}(\boldsymbol{c})\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} =: \frac{C_{2,1}}{\nu}\,\rho T^2 \!\left(\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk} -\tfrac{2}{3}\delta_{ij}\delta_{k\ell}\right), \tag{110} \\ & \int_{\mathbb{R}^3}\mathcal{B}_i(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[|\boldsymbol{c}|^2 c_k\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} =: -\frac{C_{1,1}}{\nu}\,\rho T^3\,\delta_{ik}, \tag{111} \\ & \int_{\mathbb{R}^3}\mathcal{B}_i(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[|\boldsymbol{c}|^4 c_k\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} =: -\frac{C_{1,2}}{\nu}\,\rho T^4\,\delta_{ik}, \tag{112} \\ &\int_{\mathbb{R}^3}\mathcal{A}_{ij}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[|\boldsymbol{c}|^2 c_k\,\mathcal{M}\bigr] = 0, \qquad \int_{\mathbb{R}^3}\mathcal{A}_{ij}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[|\boldsymbol{c}|^4 c_k\,\mathcal{M}\bigr] = 0, \\ & \int_{\mathbb{R}^3}\mathcal{B}_i(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{A}_{k\ell}(\boldsymbol{c})\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}}= 0, \qquad \int_{\mathbb{R}^3}\mathcal{B}_i(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[|\boldsymbol{c}|^2\mathcal{A}_{k\ell}(\boldsymbol{c})\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}}= 0, \end{align}\] where \(\nu\) is the local collision frequency and \(C_{2,0}\), \(C_{2,1}\), \(C_{1,1}\), \(C_{1,2}\) are scalar coefficients determined by the linearized collision operator. Therefore, with the definition \[\sigma^{(\nabla)}_{ij} = \int_{\mathbb{R}^3} \mathcal{A}_{ij}(\boldsymbol{c})\, g_1^{(\nabla)} \,\mathrm{d}{\boldsymbol{v}}, \qquad q^{(\nabla)}_i = \int_{\mathbb{R}^3} \mathcal{B}_i(\boldsymbol{c})\, g_1^{(\nabla)} \,\mathrm{d}{\boldsymbol{v}}, \label{eq:constitutive95def}\tag{113}\] it holds that \[\tag{114} \begin{align} &\sigma^{(\rho)}_{ij}\triangleq \int_{\mathbb{R}^3}\mathcal{A}_{ij}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{T}_{\rho}\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} = \beta_k A_\ell \int_{\mathbb{R}^3}\mathcal{A}_{ij}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{A}_{k\ell}(\boldsymbol{c})\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} = \frac{C_{2,0}}{\nu} \bigl(\left(\beta \otimes A\right)_{\mathrm{stf}}\bigr)_{ij}, \tag{115} \\ & q^{(\rho)}_i \triangleq \int_{\mathbb{R}^3}\mathcal{B}_{i}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{T}_{\rho}\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} = \vartheta A_k \int_{\mathbb{R}^3}\mathcal{B}_i(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[|\boldsymbol{c}|^2 c_k\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} = -\frac{C_{1,1}}{\nu}\rho T^3\,\vartheta A_i. \tag{116} \end{align}\] Similarly, we can obtain that \[\tag{117} \begin{align} & \sigma^{(u)}_{ij} \triangleq \int_{\mathbb{R}^3}\mathcal{A}_{ij}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{T}_{u}\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} = \frac{C_{2,0}}{\nu}\left(\left(S^{(u)}\right)_{\mathrm{stf}}\right)_{ij} + \frac{C_{2,1}}{\nu}\rho T^2\left(\left(R^{(u)}\right)_{\mathrm{stf}}\right)_{ij}, \tag{118}\\ & q^{(u)}_i \triangleq \int_{\mathbb{R}^3}\mathcal{B}_{i}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{T}_{u}\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}}= -\frac{C_{1,1}}{\nu}\rho T^3\,C^{(u,3)}_i, \tag{119} \end{align}\] and \[\tag{120} \begin{align} & \sigma^{(T)}_{ij} \triangleq \int_{\mathbb{R}^3}\mathcal{A}_{ij}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{T}_{T}\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}}= \frac{C_{2,0}}{\nu}\left(\left(S^{(T)}\right)_{\mathrm{stf}}\right)_{ij} + \frac{C_{2,1}}{\nu}\rho T^2\left(\left(R^{(T)}\right)_{\mathrm{stf}}\right)_{ij}, \tag{121} \\ & q^{(T)}_i \triangleq \int_{\mathbb{R}^3}\mathcal{B}_{i}(\boldsymbol{c})\, \mathcal{L}^{-1}\bigl[\mathcal{T}_{T}\,\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} = -\frac{C_{1,1}}{\nu}\rho T^3\,C^{(T,3)}_i -\frac{C_{1,2}}{\nu}\rho T^4\,C^{(T,5)}_i. \tag{122} \end{align}\] Collecting 114 , 117 , and 120 , we finally achieve that \[\begin{align} \sigma^{(\nabla)} = \frac{C_{2,0}}{\nu}\,\mathbf{S} + \frac{C_{2,1}}{\nu}\,\mathbf{R}, \qquad q^{(\nabla)} = -\frac{1}{2\nu}\bigl(C_{1,1}\,\boldsymbol{P} + C_{1,2}\,\boldsymbol{Q}\bigr), \end{align}\] where \(\mathbf{S}\), \(\mathbf{R}\), \(\boldsymbol{P}\), \(\boldsymbol{Q}\) are defined in 41 and 44 .

7.2 Calculating the coefficients \(C_{i,j}\)↩︎

We first present the integrals of the Maxwellian as below, which are utilized to calculate \(C_{i,j}\). \[\begin{gather} \label{eq:gauss} \int_{\mathbb{R}^3} c_ic_j \mathcal{M}\,\mathrm{d}{\boldsymbol{v}} = \rho T\delta_{ij}, \qquad \int_{\mathbb{R}^3} c_ic_jc_kc_\ell \mathcal{M}\,\mathrm{d}{\boldsymbol{v}} = \rho T^2(\delta_{ij}\delta_{k\ell} +\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk}), \\ \int_{\mathbb{R}^3} |\boldsymbol{c}|^2c_ic_j \mathcal{M}\,\mathrm{d}{\boldsymbol{v}} = 5\rho T^2\delta_{ij},\quad \int_{\mathbb{R}^3} |\boldsymbol{c}|^4c_ic_j \mathcal{M}\,\mathrm{d}{\boldsymbol{v}} = 35\rho T^3\delta_{ij}, \quad \int_{\mathbb{R}^3} |\boldsymbol{c}|^6c_ic_j \mathcal{M}\,\mathrm{d}{\boldsymbol{v}} = 315\rho T^4\delta_{ij}. \end{gather}\tag{123}\] To obtain the coefficients \(C_{i,j}\) in 108 , the technique of utilizing the BGK operator \(\mathcal{L}_{\rm BGK}\) instead of the linearized collision operator \(\mathcal{L}\) is applied. The pseudoinverse of \(\mathcal{L}_{\rm BGK}\) has the form as \[\mathcal{L}_{\mathrm{BGK}}^{-1}[g] = -\frac{1}{\nu}(I-\mathcal{P})g. \label{eq:BGK95inv951}\tag{124}\] Then, it is straightforward to verify \[\label{eq:orth95A} \int_{\mathbb{R}^3} \phi({\boldsymbol{v}}) \mathcal{A}_{k\ell}\mathcal{M} \; \mathrm{d}{\boldsymbol{v}}= 0, \qquad \phi({\boldsymbol{v}}) = (1, {\boldsymbol{v}}, |{\boldsymbol{v}}|^2)^T,\tag{125}\] so that \[\label{eq:c2095L} (I-\mathcal{P})[\mathcal{A}_{k\ell}\mathcal{M}] = \mathcal{A}_{k\ell}\mathcal{M},\qquad \mathcal{L}_{\mathrm{BGK}}^{-1}[\mathcal{A}_{k\ell}\mathcal{M}] = -\frac{1}{\nu}\mathcal{A}_{k\ell}\mathcal{M}.\tag{126}\] Substituting 126 into 109 , with 123 it yields \[\begin{align} \int_{\mathbb{R}^3}\mathcal{A}_{ij}\,\mathcal{L}_{\mathrm{BGK}}^{-1}[\mathcal{A}_{k\ell}\mathcal{M}] \,\mathrm{d}{\boldsymbol{v}} &= -\frac{1}{\nu}\int\mathcal{A}_{ij}\mathcal{A}_{k\ell} \mathcal{M}\,\mathrm{d}{\boldsymbol{v}} = -\frac{\rho T^2}{\nu} \!\left(\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk} -\frac{2}{3}\delta_{ij}\delta_{k\ell}\right). \label{eq:C20} \end{align}\tag{127}\] Comparing 127 with 109 , it gives \[\label{eq:C9520} C_{2,0}=-\rho T^2.\tag{128}\] Similarly, it holds for \(|\boldsymbol{c}|^2\mathcal{A}_{k\ell}\mathcal{M}\) that \[\label{eq:orth95bA} \int_{\mathbb{R}^3} \phi({\boldsymbol{v}}) |\boldsymbol{c}|^2\mathcal{A}_{k\ell}\mathcal{M} \; \mathrm{d}{\boldsymbol{v}}= 0,\tag{129}\] with 124 and 123 , we can derive that \[\begin{align} \int_{\mathbb{R}^3} \mathcal{A}_{ij}\, \mathcal{L}_{\mathrm{BGK}}^{-1}\bigl[|\boldsymbol{c}|^2\mathcal{A}_{k\ell}\mathcal{M}\bigr] \,\mathrm{d}{\boldsymbol{v}} = -\frac{7\rho T^3}{\nu} \left(\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk} -\frac{2}{3}\delta_{ij}\delta_{k\ell}\right). \label{eq:C21} \end{align}\tag{130}\] Comparing 130 with 110 , it gives \[\label{eq:C9521} C_{2,1} = -7T.\tag{131}\] Moreover, for \(|\boldsymbol{c}|^2c_k\mathcal{M}\), and \(|\boldsymbol{c}|^4c_k\mathcal{M}\), we have that \[\label{eq:bcM} \int_{\mathbb{R}^3} c_i|\boldsymbol{c}|^2c_k\mathcal{M}\,\mathrm{d}{\boldsymbol{v}}= 5\rho T^2\delta_{ik}, \qquad \int_{\mathbb{R}^3} c_i\cdot|\boldsymbol{c}|^4c_k\mathcal{M}\,\mathrm{d}{\boldsymbol{v}} = 35\rho T^3\delta_{ik}.\tag{132}\] Then, it holds that \[\mathcal{P}[|\boldsymbol{c}|^2c_k\mathcal{M}] = 5T c_k\mathcal{M}, \qquad \mathcal{P}[|\boldsymbol{c}|^4c_k\mathcal{M}] = 35T^2\,c_k\mathcal{M}. \label{eq:C1195L}\tag{133}\] Using 124 ,123 , 132 and 133 , we obtain \[\int_{\mathbb{R}^3} \mathcal{B}_i\, \mathcal{L}_{\mathrm{BGK}}^{-1}[|\boldsymbol{c}|^2c_k\mathcal{M}] \,\mathrm{d}{\boldsymbol{v}} = -\frac{5\rho T^3}{\nu}\delta_{ik}, \qquad \int_{\mathbb{R}^3} \mathcal{B}_i\,\mathcal{L}_{\mathrm{BGK}}^{-1}[|\boldsymbol{c}|^4c_k\mathcal{M}] \,\mathrm{d}{\boldsymbol{v}} = -\frac{70\rho T^4}{\nu}\,\delta_{ik}. \label{eq:C11}\tag{134}\] With 134 and 108 , it gives \[\label{eq:C1195C12} C_{1,1} = 5, \qquad C_{1,2} = 70.\tag{135}\]

In summary, with the BGK approximation, the coefficients \(C_{i,j}\) take explicit values as \[C_{2,0} = -\rho T^2,\quad C_{2,1} = -7T,\quad C_{1,1} = 5,\quad C_{1,2} = 70. \label{eq:coeff95values}\tag{136}\]

References↩︎

---

1. Beijing Computational Science Research Center, Beijing, China, email: zhangpei@csrc.ac.cn.↩︎

2. Department of Mathematics, National University of Singapore, Singapore, 119076, email: matcz@nus.edu.sg.↩︎

3. Beijing Computational Science Research Center, Beijing, China, email: ylwang@csrc.ac.cn.↩︎