Abstract
The low-variance direct simulation Monte Carlo (LVDSMC) is a powerful method to simulate low-speed rarefied gas flows. However, in the near-continuum flow regime, due to limitations on the time step and spatial cell size, it takes plenty of time to find the steady-state solution. Here we remove these deficiencies by coupling the LVDSMC with the general synthetic iterative scheme (GSIS) which permits the simulation at the hydrodynamic scale rather than the much smaller kinetic scale. As a proof of concept, we propose the stochastic-deterministic coupling method based on the Bhatnagar-Gross-Krook kinetic model. First, macroscopic synthetic equations are derived exactly from the kinetic equation, which not only contain the Navier-Stokes-Fourier constitutive relation, but also encompass the higher-order terms describing the rarefaction effects. Then, the high-order terms are extracted from LVDSMC and fed into synthetic equations to predict the macroscopic properties which are closer to the steady-state solution than LVDSMC. Finally, the state of simulation particles in LVDSMC is updated to reflect the change of macroscopic properties. As a result, the convergence to steady state is greatly accelerated, and the restrictions on cell size and the time step are removed. We conduct the Fourier stability analysis and simulate several canonical rarefied gas flows to demonstrate the advantages of LVDSMC-GSIS: when the Knudsen number is lower than 0.1, it can use the grid size about 10 times larger than that in traditional DSMC, and it can reduce the computational cost by two orders of magnitude in the flow regime.
Similar content being viewed by others
1 Introduction
In many modern engineering problems, e.g., the space re-entry capsule, microelectromechanical systems and shale gas extraction, multiscale gas flows that span a wide range of Knudsen numbers are frequently encountered, and accurate/efficient simulation methods are urgently needed. At the macroscopic level, the Navier-Stokes-Fourier (NSF) equations, which are the pillars in traditional computational fluid dynamics, provide a mathematical model incorporating the linear constitutive relations, such as Newton’s law of stress and Fourier’s law of heat conduction. However, they are only valid in flows where the characteristic length is much larger than the mean free path of gas molecules, i.e., when the Knudsen number is small. When the Knudsen number is large, the Boltzmann equation has to be adopted, which provides a universal approach from continuum flow to free molecular flow [1]. That is, while the NSF equations are only valid in the continuum flow regime, the Boltzmann equation can be applied in the whole flow regime, from the continuum to free molecular flow regimes, for the dilute gas.
One of the most widely used methods to model the rarefied gas flow is the direct simulation Monte Carlo (DSMC) method. Although DSMC does not solve the Boltzmann equation directly [2], it has been proven that its solution converges to that of the Boltzmann equation, when the number of simulation particles tends to infinity [3, 4]. DSMC prevails due to the following advantages: it is unconditionally stable; it automatically captures the discontinuity of velocity distribution function; it is convenient to add complicated physical and chemical processes without increasing the computational cost significantly [5]. However, there are some difficulties when using the DSMC method in real applications. First, since DSMC is a stochastic method, the macroscopic information sampled from the simulation particles is inevitably subject to fluctuations. Several remedies are proposed to reduce the fluctuation. For instances, the moment-guided Monte Carlo method [6] that matches the kinetic solutions to deterministic solutions of moment equations; the variance reduction technique [7,8,9] divides the particle distribution function into the equilibrium distribution and the deviation from the equilibrium distribution, and it only simulates the deviational part from the equilibrium. When the magnitude of the deviation distribution is small compared to the equilibrium distribution, smaller statistical error is achieved when compared to the traditional DSMC. Second, due to the splitting of advection and collision, the spatial cell size and time step should be smaller than the mean free path and mean collision time of gas molecules, respectively [5], which lead to slow convergence and high computational cost for near-continuum flows. To reduce the computational cost, the asymptotic preserving DSMC method [10, 11] and the NS-DSMC hybrid method [12,13,14] are developed, which partly solve these problems.
On the contrary, the deterministic methods to solve the Boltzmann equation have made remarkable achievements over the past few decades. Based on the discrete velocity method (DVM), the unified gas-kinetic scheme (UGKS) has been developed [15,16,17,18]. And with a similar process, the discrete UGKS is constructed by applying a simpler numerical characteristic solution of the kinetic equation [19]. The convection and collision are simultaneously solved, thus the restrictions on cell size and time step are removed. When the discrete scale is smaller than the kinetic scale, it has the same mechanism as traditional DVM. When the discrete scale is much larger than the kinetic scale, the scheme is asymptotically preserving the NSF equations, making it efficient in dealing with multiscale problems [20]. However, for near-continuum flows, UGKS still needs many iterations to obtain the steady-state solution. To overcome this challenge, the general synthetic iterative scheme (GSIS) that alternately solves the steady-state macroscopic synthetic equation and mesoscopic gas kinetic equation [21, 22] is proposed, which can not only find the steady-state solutions within dozens of iterations at any Knudsen number, but also use a larger cell size than UGKS [23].
Although the existing mature deterministic methods can be used to simulate many multiscale problems accurately and efficiently, they are still inferior to the DSMC method in dealing with hypersonic flows with complicated chemical and physical processes, where the number of discrete velocities is huge. Therefore, it is desirable to have a scheme mingling the advantages of both deterministic and stochastic methods. In recent years, the stochastic particle methods based on Bhatnagar-Gross-Krook [24,25,26] model and Fokker-Planck model [27,28,29,30] have been proposed. By simplifying the collision process, the computational efficiency in the continuum regimes has been much improved compared to the original DSMC. The former preserves the NSF limit based on the BGK model [31, 32], which is demonstrated to have high-order accuracy in space and time in the continuum regime, while in the latter method, a time integration scheme is applied, which is demonstrated to be more efficient than DSMC. The unified gas kinetic wave-particle method uses the wave-particle description to recover the non-equilibrium gas distribution function [33,34,35], where the particles without collisions are selectively sampled, while the unsampled particles are calculated by the deterministic method in the next evolution process. In the continuum regime, only a few non-collision particles are sampled and the evolution is dominated by the deterministic method. Thus, this method is efficient for high-speed multiscale problems.
Since in rarefied gas flows the steady-state solution is frequently needed, it is not necessary to follow the time evolution of the kinetic equation. If there exists a scheme that directly pulls/guides the solution to the final steady state, then the computational efficiency will be further improved on top of the above-mentioned methods. The deterministic solver GSIS serves this purpose. Inspired by the success of GSIS in multiscale problems, based on the Boltzmann equation and simplified kinetic models, a similar technique is expected to accelerate the slow convergence of DSMC in the simulation of near-continuum flows. As a preliminary work, we couple the GSIS with the low-variance (LV) DSMC to achieve asymptotic preservation and fast convergence in all flow regimes.
A brief outline of the remaining paper is sketched below. The linearized Boltzmann equation with BGK operator is introduced in Section 2; the synthetic equations are derived and the GSIS-LVDSMC coupling method is proposed in Section 3; the Fourier stability analysis is illustrated in Section 4; the accuracy and efficiency of the proposed method are assessed in canonical rarefied gas flows in Sections 5 and 6. Finally, conclusions as well as future perspectives are given in Section 7.
2 The linearized BGK equation
The BGK equation is widely used because of its simplicity:
where \(f(t, \varvec{x}, \varvec{c})\) is the velocity distribution function, with t being the time, \(\varvec{x}\) the spatial coordinates, and \(\varvec{c}\) the molecular velocity; \(\varvec{a}\) is the external acceleration, p is the gas pressure, \(\mu\) is the shear viscosity of the gas, and
is the local equilibrium distribution function, with n, \(\varvec{u}\), T, and R being the number density, flow velocity, temperature, and gas constant, respectively.
Introducing the dimensionless variables \(\hat{f}=c_0^{3}f/n_0\), \(\hat{t}=c_0 t/L\), \(\hat{\varvec{x}}=\varvec{x}/L\), \(\hat{\varvec{c}}=\varvec{c}/c_0\), and \(\hat{\varvec{a}}=\varvec{a}L/c_0^2\), where \(n_0\) is the reference number density, \(c_0=\sqrt{2RT_0}\) is the most probable speed at the reference temperature \(T_0\), L is the characteristic flow length, the BGK equation is normalized to the following form:
where \(\delta _{rp}\) is the rarefaction parameter defined as (it is related to the Knudsen number \({Kn}\) as \(\delta _{rp}=\frac{\sqrt{\pi }}{2 {Kn}}\)):
In the linearized kinetic model, the velocity distribution function can be written as the combination of the equilibrium distribution function \(\hat{f}_0=\pi ^{-3/2}\exp \left( -\hat{\varvec{c}}^2\right)\) at the reference state and the perturbed distribution function \(\hat{h}\):
where the small constant \(\alpha\) is related to the amplitude of perturbation, with \(\alpha \hat{h}\) satisfying \(|\alpha \hat{h}/\hat{f}_0|\ll 1\). However, the distribution function \(\hat{h}\) is not necessarily smaller than the fixed equilibrium distribution function \(\hat{f}_0\). With this in mind, Eq. (3) can be linearized as [22]:
where \(\rho\), \(\hat{\varvec{u}}\) and \(\tau\) are the perturbed dimensionless number density, flow velocity and temperature, respectively. Note that the normalized acceleration is also small so that the system permits linearization; the corresponding term is usually treated as the source term S. For instances, in the Poiseuille flow and thermal transpiration, due to the small pressure gradient and temperature gradient in the \(x_3\) direction, the “equivalent” source terms are
The dimensionless macroscopic quantities (which are further normalized by the constant \(\alpha\)) are defined as the moments of the perturbed velocity distribution function:
where \(\hat{\sigma }_{ij}\) and \(\hat{\varvec{q}}\) are the dimensionless deviatoric stress and heat flux, respectively; \(\delta _{ij}\) is the Kronecker delta function, and the subscripts \(i,j=1,2,3\) indicate directions in the Cartesian coordinate system.
3 The coupling of LVDSMC and GSIS
We adopt the LVDSMC method developed in Ref. [36], which solves the linearized BGK equation. Although the velocity distribution function f is non-negative, the perturbed distribution function \(\hat{h}\) can be positive or negative. Therefore, in LVDSMC, by introducing the positive and negative deviational particles, the perturbed distribution function over a single computational cell can be expressed as [36]:
where \(N_{\text {eff}}\) is the number of gas molecules represented by one simulated particle, \(s_p=\pm 1\) indicates the signs of positive and negative particles, and the subscript p indicates the information of a single simulation particle. The above equation shows that the p-th simulation particle in the cell has a position \(\hat{\varvec{x}}_{p}\) and velocity \(\hat{\varvec{c}}_{p}\). Combined with Eq. (8), the local macroscopic properties can be sampled over each computational cell (with volume \(V_{\text {cell}}\)) as:
Different from the conventional DSMC method, the number of deviational particles in LVDSMC changes with time. The adding and deleting processes of deviational particles are required in three subroutines: advection, applying body force, and collision [36]. Apart from the simple advection of deviational particles, as is performed in the conventional DSMC, additional particles have to be generated from the cell interfaces at every time step, which are sampled from a distribution determined by the difference in the macroscopic properties of adjacent cells. To enforce the boundary conditions during the advection, diffusely reflected particles are deleted at the wall and additional particles are redrawn from the distribution function evaluated at the boundary. When applying the body force, such as the pressure and thermal force, the generation of particles is based on the distribution function given by the linearized body force term. In the BGK version of LVDSMC, the collision process is implemented by deleting deviational particles with a probability determined by the mean collision time, and updating the local macroscopic properties at the same time.
LVDSMC is a stochastic method; therefore, macroscopic properties outputted from LVDSMC are the cumulative average of each time step after the flow field reaches the steady state (see the flowchart in Fig. 1 without the blue box). It is efficient when the Knudsen number is large. However, due to the restriction on cell size and time step, the computational cost is large when the Knudsen number is small, albeit the variance reduction is achieved by simulating only the deviation from equilibrium: First, the total evolution (iteration) steps required to find the steady-state solution at least scale as \(1/ {Kn}^2\). Second, the number of collisions in each cell and time step increases with the decrement of the Knudsen number, thus leading to not only more computational time on collision subroutine, but also more deviational particles.
In order to increase the efficiency of LVDSMC, the evolution steps should be reduced significantly. Since we are interested in the steady-state solution, a scheme that guides the LVDSMC evolution directly to the steady state, without considering the intermediate evolution, is highly desired. Furthermore, to remove the restriction on cell size, the scheme should also asymptotically preserve the NSF limit, so that the hydrodynamic scale (i.e., characteristic flow length), which is much larger than the kinetic scale (i.e., mean free path), can be used. The recently developed GSIS perfectly meets both requirements [21, 23], but has only been successfully applied to DVM-GSIS, which is a deterministic-deterministic coupling; this kind of coupling requires many discrete velocities and hence much computational memory and time in hypersonic flow simulations. In order to reduce the computational memory and time, here we explore the possibility of a stochastic-deterministic coupling scheme called GSIS-LVDSMC.
According to Ref. [21], the essential idea of GSIS is that macroscopic synthetic equations are exactly derived from the kinetic equation, which not only contains the NSF constitutive relations, but also encapsulates high-order terms to capture the rarefaction effects. On respectively multiplying Eq. (6) with 1, \(2\hat{c}_i\), and \(\hat{c}^2-\frac{3}{2}\), and integrating with respect to the molecular velocity space, we have
Since the stress and heat flux are not closed, we again consider their governing equations by multiplying Eq. (6) with \(2 \left( \hat{c}_i\hat{c}_j-\frac{\delta _{ij}}{3}\hat{c}^2\right)\) and \(\hat{c}_i\left( \hat{c}^2-\frac{5}{2}\right)\), respectively, and integrate with respect to the molecular velocity space, resulting
and
Note that in Eqs. (12) and (13), the high-order terms \(\text {HoT}_{\sigma _{ij}}\) and \(\text {HoT}_{q_{i}}\) are computed directly from the perturbed velocity distribution functions, thus no approximations are introduced. Inevitably, the high-order terms evaluated from the stochastic method at each time step are subject to significant fluctuations, which may lead to numerical instability. In order to reduce the fluctuations, the time-averaged values of the high-order terms are adopted. Also note that the velocity gradients in Eq. (12) and temperature gradients in Eq. (13) can not be canceled, since the one in the high-order term is statistically sampled from LVDSMC before solving the synthetic equations, while the one in Newton/Fourier’s law will be solved from the synthetic equations, in order to guide the evolution of velocity and temperature in LVDSMC.
The boundary conditions associated with the synthetic equations are extracted from LVDSMC. For example, the macroscopic properties of the cell adjacent to the wall surface can be sampled after the subroutine of advection and body force and then served as the boundary conditions. Note that the location of the sampled macroscopic properties is usually not exactly the boundary, but half of the cell size away from the boundary.
There are two tasks in performing the collision subroutine in the ordinary LVDSMC: changing the Maxwell-Boltzmann properties and deleting deviational particles. In GSIS-LVDSMC, the local Maxwell-Boltzmann distributions are updated by obtaining the macroscopic flow properties from the synthetic equations, while the process of deleting deviational particles is exactly the same as that in LVDSMC.
Since the steady state synthetic equations with Newton’s law of stress and Fourier’s law of heat conduction lead to diffusion-type of equations for the flow velocity and temperature, the information of updated deviational particles will reach the steady state much faster than the evolution of pure mesoscopic equations in the near-continuum flow regime, so that the fast convergence can be achieved by coupling the LVDSMC and GSIS. The flowchart of GSIS-LVDSMC algorithm is visualized in Fig. 1, and summarized as follows:
-
1.
System initialization with local Maxwellian distributions (deviational particles number \(= 0\));
-
2.
Half advection and linearized body force step as they are in the original LVDSMC [7];
-
3.
Coupling of GSIS for fast convergence:
-
i.
Calculate the time-averaged high-order terms \(\text {HoT}_{\sigma _{ij}}\) and \(\text {HoT}_{q_{i}}\) based on Eqs. (12) and (13) in all cells, and sample the macroscopic properties that need to be used as boundary conditions in the synthetic equations;
-
ii.
Solve the synthetic macroscopic equations (11) with constitutive relations (12) and (13), and boundary conditions sampled from deviational particles;
-
i.
-
4.
Full collision step:
-
i.
Update the local macroscopic properties \(\rho\), \(\hat{\varvec{u}}\) and \(\tau\) based on the solutions of the synthetic equations;
-
ii.
Delete deviational particles based on the probability determined by mean collision time as that in the original LVDSMC.
-
i.
-
5.
Another half advection and linearized body force step (repeat step 2);
-
6.
Sample and accumulate the macroscopic properties based on Eq. (10), and calculate the relative error of time-averaged values;
-
7.
Repeat steps 2-6 until the convergence criteria are met.
4 Fourier stability analysis
Although LVDSMC is a time-explicit scheme solving unsteady problems, GSIS-LVDSMC proposed in this work is dedicated to finding steady state solutions with fast convergence. Here we conduct the Fourier stability analysis to demonstrate that the GSIS-LVDSMC, as a combination of a time-dependent mesoscopic solver and a time-independent macroscopic solver, is always stable and has a convergence rate faster than LVDSMC.
4.1 Convergence rate of LVDSMC
The Fourier stability analysis is usually applied to investigate the efficiency of different numerical schemes of deterministic methods. Nevertheless, LVDSMC solves the distribution function of the linearized kinetic equations, and shares the same stability feature as the deterministic method where the temporal accuracy is of second-order. In order to analyze the convergence rate of macroscopic properties in LVDSMC, we first consider the temporal discretization of Eq. (6):
where the subscript k represents the macroscopic properties obtained at the time \(\hat{t}^k\), \(\varDelta \hat{t}=\hat{t}^{k+1}-\hat{t}^k\) is the time step. Note that the coefficient 1/2 at the right-hand-side of Eq. (14) represents the Crank-Nicolson scheme, and it gives the second-order accuracy in time.
Then, the Fourier stability analysis is applied to find the convergence rate of the above iteration [23, 37]. The error function between distribution functions at two consecutive iterations is defined as:
According to Eq. (8), the error functions for macroscopic properties \(M=[\rho ,\,\varvec{u},\,\tau ]\) between two consecutive iteration steps are expressed as:
where
On substituting Eqs. (15) and (16) into Eqs. (6) and (14), we have:
The error decay rate e is determined by performing the Fourier stability analysis and seeking the eigenfunctions \(y(\hat{\varvec{c}})\) and \(\alpha _M = [\alpha _\rho ,\,\varvec{\alpha }_{\varvec{u}},\,\alpha _\tau ]\) of the following forms:
where i is the imaginary unit and \(\varvec{\theta }=[\theta _1,\,\theta _2,\,\theta _3]\) is the wave vector of perturbation. Note that the factor \(e^k\), rather than \(e^{k+1}\), emerges in Eq. (19), since the error function \(Y^{k+1}(\hat{\varvec{x}},\hat{\varvec{c}})\) is obtained by macroscopic properties \({\it{\Phi}} ^{k}(\hat{\varvec{x}})\) in the k-th step. From Eqs. (19) and (20), it can be seen that the iteration scheme is stable when \(|e|\le 1\), and the fast (slow) convergence occurs when |e| approaches zero (one).
Then, by substituting the expressions of \(Y^{k+1}\) and \({\it{\Phi}} ^{k+1}\) into Eqs. (16) and (18), we have:
where
Now we obtain five linear algebraic equations for five unknown elements in \(\alpha _M\), which can be expressed as:
where the superscript \(\top\) represents the transpose operator; the matrix \(C_5\) is determined as follows:
The error decay rate e can be further obtained by numerically finding the eigenvalues of the matrix \(C_5\) and taking the maximum absolute value of e, as shown in Fig. 2. It can be seen from Eq. (23) that the function \(y_0(\hat{\varvec{c}})\) approaches zero when the Knudsen number is large, and hence e goes to zero, which implies that LVDSMC has a fast convergence rate in solving highly rarefied gas flows. On the other hand, the function \(y_0(\hat{\varvec{c}})\) approaches \(\hat{f}_0\) when the Knudsen number is small, and the matrix \(C_5\) has the maximum eigenvalue going to 1. Thus, the convergence rate of LVDSMC in the near-continuum regime is extremely slow.
4.2 Convergence rate of GSIS-LVDSMC
In GSIS-LVDSMC, to boost the convergence of iteration, additional synthetic equations (11) are solved with the constitutive relation of stress (12) and heat flux (13). The stress and heat flux at the \((k+1)\)-th iteration step can be constructed in the following manner:
where the high-order terms \(\text {HoT}_{\sigma _{ij}}^{k+1/2}\) and \(\text {HoT}_{q_{i}}^{k+1/2}\) are obtained directly from the distribution function \(\hat{h}^{k+1/2}\), which have to be sampled from the particle information in LVDSMC.
In order to find the convergence rate of GSIS, the error function is defined as:
Clearly, the expression of \(y(\hat{\varvec{c}})\) is still given by Eq. (22), since the distribution function is governed by the same kinetic equation. However, the error functions \({\it{\Phi}} _M\) are not calculated from Y but from the synthetic equations. To be specific, the macroscopic properties \(\rho\), \(\varvec{u}\), \(\tau\) in Eq. (11) is replaced by \({\it{\Phi}} _{\rho }^{k+1}\), \({\it{\Phi}} _{\varvec{u}}^{k+1}\), \({\it{\Phi}} _{\tau }^{k+1}\), respectively, and \(\hat{h}^{k+1/2}\) is replaced by \(Y^{k+1/2}\) (note that the body force term vanishes because of the subtraction). On substituting Eqs. (28) and (29) into the governing equations for \({\it{\Phi}} _{M}^{k+1}\), we obtain five linear algebraic equations,
where the source terms \(S_j\) derived from the synthetic equations are:
Note that \({\it{\Theta}} = \hat{c}_1\theta _1+\hat{c}_2\theta _2+\hat{c}_3\theta _3\) and the form of \(y_0(\hat{\varvec{c}})\) is the same as Eq. (23).
Then, the error decay rate of GSIS can be obtained by solving the linear algebraic equations (30) and (31), which can be rewritten in the matrix form as:
where \(L_5\) and \(R_5\) are both \(5\times 5\) matrices. Introducing the matrix \(G_5=L_5^{-1}R_5\), the maximum eigenvalue of matrix \(G_5\) can be numerically obtained, and the error decay rates are shown in Fig. 2. It can be clearly observed that when \(1/\delta _{rp}\rightarrow 0\), the error decay rate |e| approaches zero, demonstrating that the GSIS-LVDSMC has high efficiency in the continuum regime. Meanwhile, the error decay rate of GSIS-LVDSMC is always much smaller than 1, showing the stability of the method, and even lower than that of LVDSMC.
5 Numerical tests for the Poiseuille flow and thermal transpiration
We assess the accuracy and efficiency of GSIS-LVDSMC in the Poiseuille flow and thermal transpiration, from the continuum to free-molecular flow regimes. The mechanism of convergence-boosting is discussed in terms of the convergence rate, size cell, and time step.
5.1 Algorithm validation and efficiency assessment
It is noted from the source terms (7) of the Poiseuille flow and thermal transpiration that the velocity distribution function \(\hat{h}\) in the steady state is an odd function of \(\hat{c}_3\). Thus, the perturbed density and temperature become zero. Therefore, only the flow velocity has to be solved from the synthetic equation to boost convergence. To this end, the moment equation in Eq. (11) is reduced to
and the constitutive relation (12) of the shear stress remains unchanged in the presence of the source term (7):
where the superscript * indicates the macroscopic properties obtained from local distribution function \(\hat{h}\).
Substituting Eq. (34) into Eq. (33), the diffusion-type synthetic equation for the flow velocity can be expressed as:
where the high-order terms are defined and statistically sampled from LVDSMC as follows:
Meanwhile, the flow velocity \(\hat{u}_3\) in the cells adjacent to walls are sampled from the simulation particles, which are used as the boundary conditions in solving the macroscopic equation (35).
5.1.1 1D cases
In 1D simulations of the Poiseuille flow and thermal transpiration, the half spatial region \(0\le \hat{x}_1\le 0.5\) is uniformly divided into N cells. The symmetric condition is applied at \(\hat{x}_1=0.5\); the wall at \(\hat{x}_1=0\) is fully diffuse, i.e., in this problem the velocity distribution function for particles entering into the simulation domain is zero.
The synthetic equation (35) can be transformed into the following simple form:
which is numerically solved by the central difference scheme with the same spatial discretization as in the LVDSMC:
In the first cell, \({\it{\Phi}} _1\) is the combination of the local velocity sampled at each time step and the time-averaged high-order terms, which is sampled from LVDSMC and used as the boundary condition in solving Eq. (38). Due to symmetry, values in the N-th cell are the same as those in the \(N+1\) virtual cell: \({\it{\Phi}} _{N}={\it{\Phi}} _{N+1}\).
The relative error in flow velocity between two successive steps, k and \(k+1\), is recorded during the simulation:
In the transition state, the velocity \(\hat{u}_3\) is the time-averaged value sampled from the beginning of a simulation, and the system is regarded as reaching a steady state when \(\varepsilon <10^{-3}\). After that, the time averaging of all macroscopic values is restarted, and the simulation is terminated when \(\varepsilon <10^{-6}\). We set a maximum number of steps \(N_{\text {step,max}}\) to stop the simulations, even if the convergence criterion is not met, meaning that the simulation cannot be finished within an acceptable computational time.
Since the high-order terms (36) are essential to accurately capture the rarefaction effects, we compare their values obtained from GSIS-LVDSMC to the reference values from GSIS-DVM [22]. Figure 3 shows good agreement in \(F_{2,0,1}\): only tiny deviations occur when \({Kn}=0.1\), due to the relative significant fluctuations of the particle method in the GSIS-LVDSMC algorithm.
Velocity profiles of the Poiseuille flow and thermal transpiration obtained from GSIS-LVDSMC and GSIS-DVM are compared in Fig. 4. In the Poiseuille flow, excellent agreement of velocity profiles is seen for all Knudsen numbers. In the thermal transpiration, although the high accuracy of GSIS-LVDSMC is demonstrated in the free-molecular and transition flow regimes, it predicts higher values of creep velocity than those from GSIS-DVM when \({Kn}=0.01\). This is because the flow velocity approaches zero when the Knudsen number decreases, and large fluctuation leads to inaccuracy of the macroscopic properties and high-order terms sampled from LVDSMC, when the sample size is not large enough.
Next, we calculate the mass flow rate of the Poiseuille flow to quantify the accuracy and efficiency of GSIS-LVDSMC, where the dimensionless mass flow rate \(\dot{m}_p\) is defined as
Table 1 compares the essential simulation parameters and computational efficiency between GSIS-LVDSMC and LVDSMC. It is shown that GSIS-LVDSMC and LVDSMC have the same efficiency when the system is in high nonequilibrium, e.g., when \({Kn}=1\) and 10, and these two methods require the same number of times steps to converge in both transition and steady states. Besides, since the macroscopic equations have to be solved at each iteration step in GSIS-LVDSMC, it takes a little more CPU time than LVDSMC. However, as the Knudsen number decreases, the computational cost increases dramatically in LVDSMC. Remarkably, the GSIS-LVDSMC improves the efficiency by nearly 50 times when \({Kn}=0.01\); meanwhile, the accuracy is also improved, leading to the relative error lower than 0.02%, compared to 6.22% in LVDSMC.
5.1.2 2D cases
Consider the 2D Poiseuille flow in an infinite long channel with a square cross section. Due to symmetry, only the lower left quarter of the cross-section is simulated (\(0\le \hat{x}_1\le 0.5\), \(0\le \hat{x}_2\le 0.5\)), which is uniformly discretized into \(N_{x_1}\times N_{x_2}\) distributed spatial cells. The bottom (\(0\le \hat{x}_1\le 0.5\), \(\hat{x}_2=0\)) and left (\(\hat{x}_1=0\), \(0\le \hat{x}_2\le 0.5\)) boundaries are fully diffuse, while the top (\(0\le \hat{x}_1\le 0.5\), \(\hat{x}_2=0.5\)) and right (\(\hat{x}_1=0.5\), \(0\le \hat{x}_2\le 0.5\)) boundaries satisfy the symmetric conditions.
Based on the central difference scheme, Eq. (35) can be solved as:
where \(\hat{U}\) represents \(\hat{u}_3\) for simplicity, M, R, N represent the high-order terms \(\frac{1}{4}F_{2,0,1}\), \(\frac{2}{4}F_{1,1,1}\), \(\frac{1}{4}F_{0,2,1}\) in Eq. (36), respectively; i, j are the indexes for the cells in the \(x_1\) and \(x_2\) directions, respectively, and \(\alpha ={\it{\Delta}} x_2/{\it{\Delta}} x_1\) is the cell aspect ratio. The velocities sampled from LVDSMC and the symmetric conditions provide the boundary conditions:
-
When \(i=1\) or \(j=1\), the flow velocities \(\hat{U}_{i,j}\) sampled at each time step and the time-averaged high-order terms M, R, N are applied to solve the macroscopic equation (41).
-
When \(i=N_{x_1}\), the symmetric conditions read: \(\hat{U}_{i+1,j}=\hat{U}_{i,j}\), \(N_{i+1,j}=N_{i,j}\), \(R_{i+1,j}=-R_{i,j}\); When \(j=N_{x_2}\), the symmetric conditions read \(\hat{U}_{i,j+1}=\hat{U}_{i,j}\), \(M_{i,j+1}=M_{i,j}\), \(R_{i,j+1}=-R_{i,j}\).
The relative error in flow velocity between two successive time steps defined in Eq. (39) is extended in 2D cases as:
where the velocity \(\hat{U}\) is the time-averaged value, and the simulation is switched to steady state and terminated when \(\varepsilon <10^{-3}\) and \(\varepsilon <10^{-6}\), respectively.
Figure 5 shows the velocity contours from GSIS-LVDSMC and GSIS-DVM when \({Kn}=0.01,0.1,1\), and 10. Good agreement is achieved for all Knudsen numbers. The system parameters, computational accuracy and efficiency are compared in Table 2. Similar to the 1D case, the efficiency of GSIS-LVDSMC is approximately the same as that of LVDSMC when the Knudsen number is large. When \({Kn}=0.01\), results predicted by LVDSMC are not converged after \(N_{\text {step,max}}=5\times 10^5\) steps, which has a relative difference of around 42% to the reference solution. On the contrary, GSIS-LVDSMC obtains accurate results within \(4\times 10^4\) time steps. Therefore, the GSIS-LVDSMC is at least 50 times faster than LVDSMC.
5.2 Convergence rate
We take the 1D Poiseuille flow to further discuss the mechanisms of convergence-boosting in GSIS-LVDSMC. As shown in Table 1, the numbers of iteration steps required in GSIS-LVDSMC and LVDSMC are approximately the same when Kn is large, implying the same computational efficiency in simulating the high nonequilibrium gas flow. However, in the slip flow regime (\({Kn}=0.01\)), GSIS-LVDSMC is about 50 times more efficient than LVDSMC.
Figure 6 compares the convergence history of GSIS-LVDSMC and LVDSMC when \({Kn}=0.01\). Both the flow velocity and the mass flow rate are time-averaged values sampled from the beginning of simulations. Therefore, the number of time steps for the system to reach the steady state can be determined. The GSIS-LVDSMC takes around \(10^4\) iteration steps to make \(\varepsilon\) below \(10^{-3}\) and \(\dot{m}_p\) converged. However, \(3\times 10^5\) time steps are required in LVDSMC to reach the same criteria, yet the oscillation in mass flow rate is still much stronger. This can be understood as follows. In such a slip flow regime, frequent intermolecular collisions slow down the evolution of the gas flow from its initial state, i.e., slow down the information exchange across the computational domain. Also, LVDSMC introduces significant fluctuations in the flow velocities that lead to the oscillation in mass flow rate, as shown in Fig. 6(b). On the contrary, due to the coupling of macroscopic synthetic equations (which is of diffusion-type), which exchange the gas information across the whole computational domain, correct the gas properties, and drive the molecular distribution to the final solution quickly, the GSIS-LVDSMC greatly reduces the simulation step.
Figure 7 compares the decay of relative error \(\varepsilon\) sampled from the two algorithms when the steady state is reached. When \({Kn}=0.1\), the GSIS-LVDSMC and LVDSMC need about \(2.5\times 10^5\) and \(3\times 10^5\) time steps, respectively, to make \(\varepsilon <10^{-6}\). When \({Kn}=0.01\), the time steps become \(7\times 10^4\) and \(2\times 10^5\), respectively. Theoretically, there are two main factors that could influence the number of time steps: the particle number in each cell, and the correlation between two successive sampling steps. It can be seen from Table 1 that, the particle numbers are almost the same in GSIS-LVDSMC and LVDSMC. Therefore, we draw the conclusion that the coupling of synthetic equations and adjusting the simulation particles accordingly reduce the correlation of successive time steps and hence the fluctuations. Clearly, GSIS-LVDSMC not only accelerates the transition state of the particle method, but also achieves higher accuracy with fewer sample steps in the steady state.
5.3 Cell size and time step
In LVDSMC, due to the splitting of advection and collision, the cell size and time step should be smaller than the kinetic scales (i.e., the mean free path and mean collision time of gas molecules, respectively). Consequently, the cell number and the computational cost increase dramatically when the system is in the near continuum regime. On the contrary, in GSIS-LVDSMC, these restrictions are removed by coupling the synthetic equations, in which the NSF equations are predominant in the continuum limit so that the hydrodynamic scale (which is much larger than the kinetic scale) can be used. Detailed evidences are given below.
Figure 8(a) shows the mass flow rate obtained from the two methods when the cell number is changed from 10 to 300, which correspond to the cell size 10 and 1/3 times of the molecular mean free path, respectively. It can be seen that the mass flow rates from LVDSMC with cell size larger than the mean free path are wrong, which is due to the large numerical dissipation. When the cell size decreases to 1/3 of the mean free path, the result has a relative difference of 5% to the reference solution. On the other hand, results from GSIS-LVDSMC converge much faster and monotonically when the cell size is decreased. The relative errors in mass flow rate are less than 1% as long as the cell size is smaller than 3 times of the mean free path in this problem. In GSIS-LVDSMC, when the cell size is larger than 3 times of the mean free path, the deviation is mainly caused by the inaccurate boundary conditions sampled from LVDSMC. If finer grids near the boundary are applied, we can obtain accurate results by using fewer cells in the bulk region.
Figure 8(b) shows the mass flow rate of 1D Poiseuille flow with the time step \(N{\it{\Delta}} \hat{t}\) varying from 0.01 to 5, where N is the number of grid cells (\(N=100\) and \(N=300\) are used in GSIS-LVDSMC and LVDSMC, respectively). Therefore, \(N{\it{\Delta}} \hat{t}=1\) means that a particle with the most probable speed \(c_0\) travels a distance of one simulation cell during the time step \({\it{\Delta}} \hat{t}\). In GSIS-LVDSMC, results are not sensitive to the time step, and the relative error in mass flow rate is around 1% even when \({\it{\Delta}} \hat{t}\) is 5 times of the mean collision time. However, as it is commonly acknowledged, the time step in LVDSMC should be smaller than 1/3 of the mean collision time to guarantee a reliable result, which corresponds to \(N{\it{\Delta}} \hat{t}=1\) here (the grid size is 1/3 of the mean free path in this case). As shown in Fig. 8(b), larger \({\it{\Delta}} \hat{t}\) in LVDSMC leads to wrong mass flow rate. However, smaller \({\it{\Delta}} \hat{t}\) requires so tremendous computational cost that the results of LVDSMC are not converged within the maximum simulation time step \(N_{\text {step,max}}\).
6 Numerical tests for planar Fourier flow
Consider the planar Fourier flow of gas between two infinite parallel plates located at \(\hat{x}_1=0\) and \(\hat{x}_1=1\) with perturbed temperature \(\tau =-0.5\) and 0.5 respectively. Assuming the symmetric condition \(\hat{h}(-\hat{c}_1,\hat{c}_2,\hat{c}_3)=-\hat{h}(\hat{c}_1,\hat{c}_2,\hat{c}_3)\) can be applied at \(\hat{x}_1=0.5\). Based on the synthetic equations (11), (12) and (13), we have \(\hat{\varvec{u}}=0\), \(\hat{\sigma }_{ij}=0~(i\ne j)\), \(\hat{q}_2=\hat{q}_3=0\), \({\partial \hat{q}_1}/{\partial \hat{x}_1}=0\), and the constitutive relations are simplified to:
where the high-order terms \(F_{q_1}\) and \(F_{\sigma _{11}}\) are statistically sampled from LVDSMC as follows:
The problem can be solved by either following the coupling algorithm strictly, or using the simplified procedure as follows. Since the heat flux \(\hat{q}_1\) is a constant through the simulation domain, it is sampled according to Eqs. (13) and (9) as:
where \(\langle \cdot \rangle\) represents the ensemble average over the entire simulation domain. Therefore, both \(F_{q_{1}}\) and \(\hat{q}_1\) are sampled from the distribution of particles before solving macroscopic equations, which is indicated by step \((k+1/2)\) below. And then, the perturbed temperature at the \((k+1)\)-th step can be determined from Eq. (43) satisfying:
Meanwhile, the stress \(\hat{\sigma }_{11}^{(k+1/2)}\) is obtained by calculating the spatial derivative of term \(F_{\sigma _{11}}\) from Eq. (43). Then, the perturbed density \(\rho\) at the \((k+1)\)-th step can be solved based on Eqs. (11) and (8):
where the constant can be evaluated at \(\hat{x}_1=0.5\) using the symmetric condition, which is found to be zero.
The relative error between successive time steps k and \(k+1\) is defined as:
The simulation is switched to steady state and terminated when \(\varepsilon <10^{-3}\) and \(\varepsilon <10^{-6}\) are satisfied, respectively.
6.1 Transition flow regime: \({Kn}=0.1, 1, 10\)
When \({Kn}=0.1, 1, 10\), the system parameters are set to be the same for both GSIS-LVDSMC and LVDSMC, where 25 uniform cells are used in the computational domain (\(0\le \hat{x}_1 \le 0.5\)) and \(N{\it{\Delta}} \hat{t}=1\). Results obtained by GSIS-DVM are regarded as the reference solutions to assess the accuracy of GSIS-LVDSMC and LVDSMC. First, variations of terms \(F_{\sigma _{11}}\) and \(F_{q_1}\) are shown in Fig. 9. Based on Eq. (43), high-order terms are the spatial derivatives of \(F_{\sigma _{11}}\) and \(F_{q_1}\); thus the nonequilibrium effects mainly occur through the constitutive relation of heat flux in the region several mean free path away from the walls. When the Knudsen number decreases, high-order terms begin to vanish in the central region of the Fourier flow. Figure 9(c) and (d) show the perturbed converged density and temperature, where the accuracy of GSIS-LVDSMC is clearly demonstrated. Additionally, when compared to the reference solution obtained by GSIS-DVM, the heat flux predicted by GSIS-LVDSMC has only \(0.92\%\), \(0.27\%\) and \(0.087\%\) relative difference, when \({Kn}=0.1,\,1\) and 10, respectively, which shows nearly the same accuracy compared with the original LVDSMC algorithm.
The simulation parameters, especially the number of iteration steps and computational cost, are summarized in Table 3. Similar to the situations in the Poiseuille flow and thermal transpiration, the efficiency of GSIS-LVDSMC is close to that of LVDSMC when the Knudsen number is large, although the CPU time cost in GSIS-LVDSMC is relatively large due to the solving of synthetic equations. However, when \({Kn}=0.1\), the number of time steps (\(3\times 10^4\)) in LVDSMC significantly increases, but the GSIS-LVDSMC takes only 2000 iteration steps to reach the steady state in the transition regime. Eventually, the CPU time can be reduced by 6 times when GSIS-LVDSMC is applied. The computational cost is reduced by about 1800 times when \(Kn=0.01\), which will be even bigger when the Knudsen number is further reduced.
6.2 Slip flow regime: \({Kn}=0.01\)
When \({Kn}=0.01\), due to the restriction on cell size and time step, 300 spatial cells and the time step \(N{\it{\Delta}} \hat{t}=1\) are used in LVDSMC. In the GSIS-LVDSMC algorithm, 50 cells and \(N{\it{\Delta}} \hat{t}=0.5\) are applied instead. Figure 10 illustrates the difference in the decay of relative error \(\varepsilon\). The discrepancy in convergence history between the two methods appears when \({Kn}=0.1\), where GSIS-LVDSMC converges faster than LVDSMC by around \(10^6\) iteration steps. When \({Kn}=0.01\), LVDSMC needs about 10 times of the number of time steps than GSIS-LVDSMC to reach the same relative error.
Besides, the fluctuation of LVDSMC is much more significant, although the number of simulation particles in each cell is 7 times than that in GSIS-LVDSMC. Figure 11 shows the perturbed density and temperature at different time steps when \({Kn}=0.01\). GSIS-LVDSMC provides accurate solutions, with little fluctuations, after \(10^5\) time steps, which have less than 3% difference in heat flux compared to the reference solution. However, both the density and temperature from LVDSMC have significant fluctuations and discrepancies to the reference solutions, even \(N_{\text {step,max}}=10^7\) and 40 hours CPU time have been spent. In the slip flow regime, LVDSMC takes large number of time steps to make the influence from the walls pass through the entire bulk region, while the synthetic equations in GSIS-LVDSMC, because of its diffusion-type, help to pass the perturbation through the whole simulation domain immediately, thus boosting the convergence. Furthermore, the correction of macroscopic properties of the flow field also reduces fluctuations of the stochastic method, leading to a smaller sample size required in doing time-averaged sampling.
7 Conclusions
In summary, we have applied the GSIS to improve the computational efficiency of the LVDSMC method. The stability and fast convergence in all flow regimes of the proposed GSIS-LVDSMC method have been rigorously demonstrated by Fourier stability analysis. The accuracy of the coupling method is validated in the Poiseuille flow, thermal transpiration, and Fourier flow. The fast convergence of GSIS-LVDSMC is achieved by solving the synthetic macroscopic equations at each time step, which not only explicitly contain the exact constitutive relations for shear stress and heat flux extracted from LVDSMC, but also asymptotically preserve the NSF limit. The former feature guarantees the algorithm’s accuracy in all flow regimes, since the rarefaction effects are captured by high-order terms of the constitutive relations extracted from LVDSMC, the latter removes the constraints on spatial cell size and time step (if the non-uniform grids with finer grids near the boundary are applied, the expected results can be obtained by applying fewer computational cells in the bulk region). With the coupling of synthetic equations, the velocity distribution of particles is adjusted according to the solutions and thus approaches the steady state quickly. Therefore, the number of time steps required in the transition state is significantly reduced when the Knudsen number is small. Meanwhile, the sampling fluctuations are found to be much smaller in GSIS-LVDSMC when the Knudsen number is lower than 0.1, thus its efficiency is further improved.
Theoretically speaking, with small modifications, the essential idea of GSIS-LVDSMC coupling algorithm should be able to improve the computational efficiency of the conventional DSMC method in the near-continuum regime, which will have a much wider application scope than the low-variance version. We expect the essential idea in this scheme can be extended to efficiently and accurately solve the DSMC for multiscale hypersonic flows with chemical reactions, which has strong applications in space exploration.
Availability of data and materials
All data generated or analyzed during this study are included in this published article.
References
Cercignani C (2000) Rarefied gas dynamics: from basic concepts to actual calculations. Cambridge University Press, Cambridge
Bird GA (1970) Direct simulation and the Boltzmann equation. Phys Fluids 13(11):2676–2681
Li ZH, Fang M, Jiang XY et al (2013) Convergence proof of the DSMC method and the gas-kinetic unified algorithm for the Boltzmann equation. Sci China Phys Mech Astron 56(2):404–417
Wagner W (1992) A convergence proof for Bird’s direct simulation Monte Carlo method for the Boltzmann equation. J Stat Phys 66(3–4):1011–1044
Bird GA (1994) Molecular gas dynamics and the direct simulation of gas flows. Oxford University Press, Oxford
Degond P, Dimarco G, Pareschi L (2011) The moment-guided Monte Carlo method. Int J Numer Methods Fluids 67(2):189–213
Baker LL, Hadjiconstantinou NG (2005) Variance reduction for Monte Carlo solutions of the Boltzmann equation. Phys Fluids 17(5):051703
Homolle TMM, Hadjiconstantinou NG (2007) A low-variance deviational simulation Monte Carlo for the Boltzmann equation. J Comput Phys 226(2):2341–2358
Radtke GA, Hadjiconstantinou NG, Wagner W (2011) Low-noise Monte Carlo simulation of the variable hard sphere gas. Phys Fluids 23(3):030606
Pareschi L, Russo G (2000) Asymptotic preserving Monte Carlo methods for the Boltzmann equation. Transp Theory Stat Phys 29(3–5):415–430
Pareschi L, Russo G (2001) Time relaxed Monte Carlo methods for the Boltzmann equation. SIAM J Sci Comput 23(4):1253–1273
Patronis A, Lockerby DA, Borg MK et al (2013) Hybrid continuum-molecular modelling of multiscale internal gas flows. J Comput Phys 255:558–571
Stephani KA, Goldstein DB, Varghese PL (2013) A non-equilibrium surface reservoir approach for hybrid DSMC/Navier–Stokes particle generation. J Comput Phys 232(1):468–481
Zhang J, John B, Pfeiffer M et al (2019) Particle-based hybrid and multiscale methods for nonequilibrium gas flows. Adv Aerodyn 1(1):12
Huang JC, Xu K, Yu P (2012) A unified gas-kinetic scheme for continuum and rarefied flows II: Multi-dimensional cases. Commun Comput Phys 12(3):662–690
Huang JC, Xu K, Yu P (2013) A unified gas-kinetic scheme for continuum and rarefied flows III: Microflow simulations. Commun Comput Phys 14(5):1147–1173
Liu C, Xu K (2020) A unified gas-kinetic scheme for micro flow simulation based on linearized kinetic equation. Adv Aerodyn 2(1):21
Xu K, Huang JC (2010) A unified gas-kinetic scheme for continuum and rarefied flows. J Comput Phys 229(20):7747–7764
Guo Z, Xu K (2021) Progress of discrete unified gas-kinetic scheme for multiscale flows. Adv Aerodyn 3(1):6
Zhu Y, Zhong C, Xu K (2016) Implicit unified gas-kinetic scheme for steady state solutions in all flow regimes. J Comput Phys 315:16–38
Su W, Zhu L, Wang P et al (2020) Can we find steady-state solutions to multiscale rarefied gas flows within dozens of iterations? J Comput Phys 407:109245
Wu L, Zhang J, Liu H et al (2017) A fast iterative scheme for the linearized Boltzmann equation. J Comput Phys 338:431–451
Su W, Zhu L, Wu L (2020) Fast convergence and asymptotic preserving of the general synthetic iterative scheme. SIAM J Sci Comput 42(6):B1517–B1540
Fei F, Jenny P (2021) A hybrid particle approach based on the unified stochastic particle Bhatnagar-Gross-Krook and DSMC methods. J Comput Phys 424:109858
Fei F, Ma Y, Wu J et al (2021) An efficient algorithm of the unified stochastic particle Bhatnagar-Gross-Krook method for the simulation of multi-scale gas flows. Adv Aerodyn 3(1):18
Fei F, Zhang J, Li J et al (2020) A unified stochastic particle Bhatnagar-Gross-Krook method for multiscale gas flows. J Comput Phys 400:108972
Gorji MH, Jenny P (2014) An efficient particle Fokker-Planck algorithm for rarefied gas flows. J Comput Phys 262:325–343
Gorji MH, Torrilhon M, Jenny P (2011) Fokker-Planck model for computational studies of monatomic rarefied gas flows. J Fluid Mech 680:574–601
Jenny P, Torrilhon M, Heinz S (2010) A solution algorithm for the fluid dynamic equations based on a stochastic model for molecular motion. J Comput Phys 229(4):1077–1098
Pfeiffer M, Gorji MH (2017) Adaptive particle–cell algorithm for Fokker–Planck based rarefied gas flow simulations. Comput Phys Commun 213:1–8
Bhatnagar PL, Gross EP, Krook M (1954) A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys Rev 94(3):511–525
Tumuklu O, Li Z, Levin DA (2016) Particle ellipsoidal statistical bhatnagar-gross-krook approach for simulation of hypersonic shocks. AIAA J 54(12):3701–3716
Li W, Liu C, Zhu Y et al (2020) Unified gas-kinetic wave-particle methods III: Multiscale photon transport. J Comput Phys 408:109280
Liu C, Zhu Y, Xu K (2020) Unified gas-kinetic wave-particle methods I: Continuum and rarefied gas flow. J Comput Phys 401:108977
Zhu Y, Liu C, Zhong C et al (2019) Unified gas-kinetic wave-particle methods. II. Multiscale simulation on unstructured mesh. Phys Fluids 31(6):067105
Radtke GA, Hadjiconstantinou NG (2009) Variance-reduced particle simulation of the Boltzmann transport equation in the relaxation-time approximation. Phys Rev E Stat Nonlin Soft Matter Phys 79(5 Pt 2):056711
Zhu L, Pi X, Su W et al (2021) General synthetic iterative scheme for nonlinear gas kinetic simulation of multi-scale rarefied gas flows. J Comput Phys 430:110091
Acknowledgements
Not applicable.
Funding
This work is supported by the National Natural Science Foundation of China under the grant No. 12172162.
Author information
Authors and Affiliations
Contributions
QL and LW conceptualized the idea, while LL implemented and wrote the code. All authors contributed to the manuscript writing, draft or revision. All authors read and approved the final manuscript.
Corresponding authors
Ethics declarations
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Luo, L., Li, Q. & Wu, L. Boosting the convergence of low-variance DSMC by GSIS. Adv. Aerodyn. 5, 10 (2023). https://doi.org/10.1186/s42774-023-00138-0
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s42774-023-00138-0