A fast iterative scheme for the linearized Boltzmann equation

An iterative scheme can be used to find a steady-state solution to the Boltzmann equation, however, it is very slow to converge in the near-continuum flow regime. In this paper, a synthetic iterative scheme is developed to speed up the solution of the linearized Boltzmann equation. The velocity distribution function is first solved by the conventional iterative scheme, then it is corrected such that the macroscopic flow velocity is governed by a diffusion equation which is asymptotic-preserving in the Navier-Stokes limit. The efficiency of the new scheme is verified by calculating the eigenvalue of the iteration, as well as solving for Poiseuille and thermal transpiration flows. The synthetic iterative scheme is significantly faster than the conventional iterative scheme in both the transition and the near-continuum flow regimes. Moreover, due to the asymptotic-preserving properties, the SIS needs less spatial resolution in the near-continuum flow regimes, which makes it even faster than the conventional iterative scheme. Using this synthetic iterative scheme, and the fast spectral approximation of the linearized Boltzmann collision operator, Poiseuille and thermal transpiration flows between two parallel plates, through channels of circular/rectangular cross sections, and various porous media are calculated over the whole range of gas rarefaction. Finally, the flow of a Ne-Ar gas mixture is solved based on the linearized Boltzmann equation with the Lennard-Jones potential for the first time, and the difference between these results and those using hard-sphere intermolecular potential is discussed.

Iterative schemes to find steady-state solutions to the Boltzmann equation are efficient for highly rarefied gas flows, but can be very slow to converge in the near-continuum flow regime. In this paper, a synthetic iterative scheme is developed to speed up the solution of the linearized Boltzmann equation by penalizing the collision operator L into the form L = (L + Nδh) − Nδh, where δ is the gas rarefaction parameter, h is the velocity distribution function, and N is a tuning parameter controlling the convergence rate. The velocity distribution function is first solved by the conventional iterative scheme, then it is corrected such that the macroscopic flow velocity is governed by a diffusion-type equation that is asymptotic-preserving into the Navier-Stokes limit. The efficiency of this new scheme is assessed by calculating the eigenvalue of the iteration, as well as solving for Poiseuille and thermal transpiration flows. We find that the fastest convergence of our synthetic scheme for the linearized Boltzmann equation is achieved when Nδ is close to the average collision frequency. The synthetic iterative scheme is significantly faster than the conventional iterative scheme in both the transition and the near-continuum gas flow regimes. Moreover, due to its asymptotic-preserving properties, the synthetic iterative scheme does not need high spatial resolution in the near-continuum flow regime, which makes it even faster than the conventional iterative scheme. Using this synthetic scheme, with the fast spectral approximation of the linearized Boltzmann collision operator, Poiseuille and thermal transpiration flows between two parallel plates, through channels of circular/rectangular cross sections and various porous media are calculated over the whole range of gas rarefaction. Finally, the flow of a Ne-Ar gas mixture is solved based on the linearized Boltzmann equation with the Lennard-Jones intermolecular potential for the first time, and the difference between these results and those using the hard-sphere potential is discussed.

Introduction
The Boltzmann equation is fundamental to a broad range of applications from aerodynamics to microfluidics [1], and it is important to be able to solve it accurately and efficiently. Most often, the Boltzmann equation is solved by the stochastic Direct Simulation Monte Carlo (DSMC) technique, which uses a number of simulated particles to mimic the binary collisions and streaming of very large numbers of gas molecules [2]. In DSMC, the length of the spatial cell and the time step need to be smaller than the local molecular mean free path and the mean collision time, respectively, and for this reason the technique becomes very slow and costly for near-continuum flows. Although time-relaxed and asymptotic-preserving Monte Carlo methods allow larger time steps [3,4], the restriction on the size of the spatial cells has not yet been removed. The same problem, in fact, exists in deterministic numerical methods for the Boltzmann equation, where the streaming and the collisions are treated separately in the splitting scheme [5,6].
The unified gas-kinetic scheme (UGKS) provides an alternative approach. It was first developed for the Bhatnagar-Gross-Krook (BGK) kinetic model [7,8], then for the Shakhov model [9,10], and finally generalized to the Boltzmann equation [11]. It handles the streaming and binary collisions simultaneously so that, for time-dependent problems, the time step is only limited by the Courant-Friedrichs-Lewy condition. Also, the UGKS is asymptotic-preserving into the Navier-Stokes limit [12], so the length of the spatial cells can be significantly larger than the molecular mean free path. Moreover, the UGKS is a finite volume method, and the analytical integral solution of the BGK-type model enables accurate flux evaluation at the cell interface, so that the essential flow physics can be captured even with the coarse grids. These advanced properties make the UGKS (and its improved version: the discrete UGKS [13,14]) a multiscale method for efficient and accurate calculations of rarefied gas flows over a wide range of the gas rarefaction. Recently, an implicit UGKS has been proposed to eliminate the time step limitation and further improve the numerical efficiency [15].
To find a steady-state solution to the Boltzmann equation, an iterative scheme is often used. In the free-molecular flow regime, where binary collisions are negligible, an iterative scheme is efficient, because the gas molecules move in straight way (except the collision with wall surfaces) so that any disturbance at one point can be quickly felt by all other points. However, for near-continuum flows the iterations converge slowly and the results are very likely to be biased by accumulated rounding errors. Although the time and spatial steps can be large, the UGKS still needs a large number of iterations [15]. This is governed by the underlying physics: the exchange of information through streaming becomes very inefficient when binary collisions dominate. Therefore, it would be useful to develop an efficient numerical scheme to solve the Boltzmann equation that both has the asymptotic-preserving property in the Navier-Stokes limit (like the UGKS where the spatial resolution can be coarse) and converges to the steady-state rapidly.
Inspired by work on fast iterative methods for radiation transport processes [16], accelerated iterative schemes have been developed for the linearized BGK and Shahkov models [17] to overcome slow convergence in the near-continuum flow regime. The fast iterative scheme is called a "synthetic iterative scheme" (SIS) since kinetic model equations are solved in parallel with diffusion-type equations for macroscopic quantities such as the flow velocity and heat flux. The SIS has been successfully applied to Poiseuille flow in channels with two-dimensional cross sections of arbitrary shapes [18] using a BGK model for single-species gases, and flows of binary and ternary gas mixtures driven by local pressure, temperature, and concentration gradients [19][20][21]18,[22][23][24] using the McCormack model [25]. The fast convergence of the SIS is due to three factors: first, the macroscopic synthetic diffusion-type equations exchange the information very efficiently; second, the macroscopic flow quantities can be fed back into mesoscopic kinetic models; and third, macroscopic diffusion-type equations are solved more quickly than mesoscopic kinetic model equations.
In the present paper we develop a SIS to solve the linearized Boltzmann equation (LBE) for Poiseuille and thermal transpiration flows. Although for the single-species LBE these canonical flows have been extensively studied for hard-sphere [26,27], inverse power-law [28], and even Lennard-Jones [29,30] potentials, numerical results for near-continuum flows are scarce. Moreover, for gas mixtures these flows have only been solved based on the hard-sphere model in a one-dimensional geometry [31]. We will calculate these flows through two-dimensional cross sections of arbitrary shape, and investigate the influence of realistic intermolecular potentials for gas mixtures. The core methods we adopt are (i) the SIS originally developed for kinetic model equations, which is introduced in the previous paragraph, (ii) the penalization method [6,32,33], which makes the development of SIS for the LBE possible, and (iii) the fast spectral method, developed by Mouhot and Pareschi [34] and extended by us to gas mixtures and Lennard-Jones potentials [35,30], which enables the efficient and accurate computation of the linearized Boltzmann collision operator.
The rest of this paper is organized as follows. In Sec. 2, we briefly introduce the LBE for single-species gases and the conventional iterative scheme (CIS). Then, by analyzing the SIS for the BGK model, we develop a SIS for the LBE and test its performance by calculating both the eigenvalues of iterations and Poiseuille/thermal transpiration flows. We improve the efficiency of the proposed SIS by adjusting a parameter in the scheme, which can be determined prior to the numerical simulation. In Sec. 3, the SIS is used to solve rarefied gas flows in multiscale problems. In Sec. 4, the SIS in polar coordinates is proposed and numerical results of the LBE for Poiseuille flow through a tube are presented. In Sec. 5, the SIS is extended to the LBE for gas mixtures, and Poiseuille flow of a Ne-Ar mixture is solved for the first time based on the Lennard-Jones potential. In Sec. 6, we conclude with a summary of the new numerical method and future perspectives.

A synthetic scheme for the single-species LBE
Consider the steady flow of a single-species monatomic gas along a channel of arbitrary cross section in the x 1 -x 2 plane, subject to small pressure/temperature gradients in the x 3 direction. The velocity distribution function (VDF) can be is the equilibrium VDF and h(x 1 , x 2 , v) is the perturbed VDF satisfying |h/ f eq | 1. The LBE for h is: with the linearized Boltzmann collision operator [28]: is the molecular velocity vector normalized by the most probable speed v m = 2k B T 0 /m (k B is the Boltzmann constant, T 0 is the gas/wall temperature, and m is the gas molecular mass), B(|v − v * |, θ) is the collision kernel determined by the intermolecular potential [28,30], and is the equilibrium collision frequency. Finally, S is the source term: where X P and X T are the pressure and temperature gradients, respectively. For the LBE, since macroscopic quantities are proportional to X P and X T , we assume X P = X T = −1.
The macroscopic quantities of interest are the flow velocity normalized by the most probable speed: the shear stresses normalized by the equilibrium gas pressure p 0 : and the heat flux normalized by p 0 v m : The dimensionless mass flow rate M and heat flow rate Q are: where A is the area of the cross section.
The integro-differential system defined by Eqs. (2) and (3) is usually solved by the CIS. Given the value of h (k) at the k-th iteration step, the VDF at the next iteration step is calculated by solving the following equation: where derivatives with respect to spatial variables are usually approximated by a second-order upwind finite difference. The process is repeated until relative differences between successive estimates of macroscopic quantities are less than a convergence criterion . The number of iteration steps in CIS increases significantly when the ratio of the molecular mean free path to the characteristic flow length decreases, especially when the gas flow is in the near-continuum regime [17], see also the data in Table 1 below. It is our goal here to develop a fast iterative scheme to solve the LBE efficiently and accurately over the whole range of gas rarefaction.

SIS for the BGK equation
To begin with, we introduce the SIS for the BGK equation [17]. The linearized Boltzmann collision operator in Eq. (3) is replaced by that of the BGK kinetic model, yielding the following equation for the perturbed VDF h: is the rarefaction parameter, with μ being the gas shear viscosity.
Multiplying Eq. (11) by the Hermite polynomials and applying the recursion relation, a set of first-order partial differential equations can be obtained for various orders of moments [36]. Here, two relevant equations for the macroscopic flow velocity are listed: where (15) are the non-accelerated high-order moments, with H n (v) being the n-th order physicists' Hermite polynomial. The combination of Eqs. (13) and (14) leads to an equation of diffusion-type for the flow velocity: , thermal transpiration. (16b) Note that in obtaining the final equation we have used the relation ∂ P 13 /∂x 1 + ∂ P 23 /∂x 2 = 1 for Poiseuille flow and ∂ P 13 /∂x 1 + ∂ P 23 /∂x 2 = 0 for thermal transpiration flow. The SIS for the BGK equation then works as follows [17,36]: • When h (k) and U (k) 3 are known at the k-th iteration step, calculate the VDF h (k+1) by solving the following equation: • From h (k+1) , calculate the non-accelerated moments F 2,0,1 , F 1,1,1 , and F 0,2,1 .
• From h (k+1) , calculate the flow velocity U The above iterative procedure is continued until convergence. It should be emphasized that the relation ∂ P 13 /∂x 1 + ∂ P 23 /∂x 2 = 1 or 0 for Poiseuille flow or thermal transpiration flow, respectively, is crucial for the fast convergence of the SIS. This is because non-accelerated moments are negligible at large values of the rarefaction parameter δ, so the synthetic equation (16b) quickly adjusts the flow velocity to the solution of the Navier-Stokes equation, which is close to the solution of the linearized BGK equation. If the synthetic equation (16a) is used instead, with P 13 and P 23 calculated based on the VDF obtained at each iteration step, the slow convergence at large values of δ is not improved because it takes a lot of iterations to reach the condition ∂ P 13 /∂x 1 + ∂ P 23 /∂x 2 = 1 or 0; in the worst-case scenario, it may even lead to false convergence and incorrect solutions when the spatial resolution is not high enough.
However, Eq. (16b) guarantees the correctness of the solution at large values of δ, as it has the asymptotic-preserving property in the Navier-Stokes limit [12]. This point will be demonstrated in Sec. 2.6 below.

SIS for the LBE
The development of a SIS for the LBE is not straightforward, since the collision operator of the LBE is much more complicated than that of the linearized BGK model. Directly following the method in Sec. 2.1, the following diffusion-type equation for the flow velocity is obtained: which, like Eq. (16a), cannot improve the slow convergence at large values of δ.
To speed up the convergence, the relation ∂ P 13 /∂x 1 + ∂ P 23 /∂x 2 = 1 or 0 for Poiseuille flow or thermal transpiration flow, respectively, must be reflected in the diffusion-type equation. For instance, as in Eq. (16b), a term similar to −δ should appear on the right-hand-side of Eq. (19) for Poiseuille flow. To achieve this, we penalize the linearized Boltzmann collision operator by the linearized BGK operator [6], i.e., and let which is very close to Eq. (16b) for the linearized BGK equation. At large values of the rarefaction parameter δ, dv approach zero, and Eq. (22) possesses the asymptotic-preserving property in the Navier-Stokes limit [12]. Therefore, a SIS can be developed based on this equation. Note that for thermal transpiration flow, δ in Eq. (22) should be replaced by zero, as ∂ P 13 /∂x 1 + ∂ P 23 /∂x 2 = 0.
The SIS for the LBE then works as that for the BGK equations, with some changes: 3 are known at the k-th iteration step, we calculate We also calculate the VDF h (k+1/2) by solving the following equation: • From h (k+1/2) , we calculate the flow velocity U (k+1/2) , and the non-accelerated moments F 2,0,1 , F 1,1,1 , and F 0,2,1 .
• Near the boundary, we let U (k+1) , while we solve the diffusion-type equation (22) to obtain the flow velocity in the bulk.
• A correction of the VDF is introduced in accordance with the changed flow velocity: Eigenvalue ω versus the inverse rarefaction parameter 1/δ for different iterative schemes for the LBE with Maxwell and hard-sphere molecules (Note that Poiseuille and thermal transpiration flows have the same eigenvalue). Our method to calculate the eigenvalue for the SIS at large values of δ is not accurate, since the convergence is so fast (it converges after one iteration) that we have few data to calculate λ through numerical fitting. However, the trend that ω in the SIS decreases with 1/δ is clear when 1/δ → 0.

Numerical analysis of the convergence rate
Analytical solutions for the eigenvalue ω have previously been introduced in order to characterize the convergence rate of the iterative scheme for the linearized BGK equation [17,36]. However, this is more difficult for the LBE because of its intricate collision operator. Here we calculate the eigenvalue numerically in order to study the performance of both the SIS and the CIS. For simplicity, we consider a periodic system of length in the x 1 direction, while the system is homogeneous in the x 2 direction.
For the CIS described by Eq. (10), the VDF is solved in the following manner (hereafter in this section, h and U 3 should be viewed as their Fourier transforms in the x 1 direction): During iteration, the flow velocity U (k+1) 3 =´h (k+1) v 3 dv is recorded, and upon convergence the resultant series of the flow velocity is fitted by U 3 (k) = U 3∞ + Ce −λk . The eigenvalue ω is then calculated as ω = e −λ . It is obvious that the smaller ω is, the faster the convergence; the case of ω = 1 means no convergence.
For the SIS, the VDF is first updated according to Eq. (23): Then the flow velocity is calculated according to the diffusion-type equation (22) as where Finally, this flow velocity is used to correct the VDF according to Eq. (24):  For small values of the rarefaction parameter δ both schemes have the same convergence rate. However, for large values of δ, the CIS has extremely slow convergence (ω ≈ 1), while the SIS converges much faster. It is also interesting to note that the intermolecular potential greatly affects the convergence rate: at the same value of δ, the solution of the LBE for hard-sphere molecules converges faster than that for Maxwell molecules, 1 in both the SIS and the CIS.

Numerical results for spatially-inhomogeneous systems
We now present numerical simulations that demonstrate the efficiency and accuracy of the SIS for Poiseuille/thermal transpiration flows between infinite parallel plates and through a two-dimensional square channel. We first consider a gas flow between two infinite parallel plates located at x 1 = −1/2 and x 1 = 1/2 (note that x 1 has been normalized by the distance between two parallel palates ). Pressure and temperature gradients are applied in the x 3 direction only, so the flow is homogeneous in the x 2 direction and partial derivatives with respect to x 2 can be dropped. The discretization of the three-dimensional molecular velocity space, as well as the fast spectral method to solve the linearized Boltzmann collision operator, are given in Ref. [28]. We adopt the diffuse boundary condition for the gas-wall interaction. Due to symmetry, only half of the spatial region (−1/2 ≤ x 1 ≤ 0) is simulated, with a specular-reflection boundary condition at x 1 = 0, and the diffuse boundary condition h(v 2 > 0) = 0 at x 1 = −1/2. The spatial domain is divided into 100 nonuniform sections, with most of the discrete points placed near the wall: The size of the smallest section is 1.24 × 10 −6 , small enough to capture the Knudsen layer.
For the one-dimensional problem, the shear stress is P 13 = x 1 for Poiseuille flow and P 13 = 0 for thermal transpiration flow. The diffusion-type equation (22) is integrated to give the following first-order ordinary differential equation: which is solved by a second-order upwind finite difference (with a first-order scheme at the wall), with the boundary A comparison between the SIS and the CIS is tabulated in Table 1 for Poiseuille flow. The relative differences in mass/heat flow rates between the two schemes is less than 0.5%, which demonstrates the accuracy of the SIS. The superiority of the SIS over the CIS is immediately seen: for the CIS, the number of iteration steps increases rapidly with the rarefaction parameter, while for the SIS it only slightly increases with δ in the free-molecular and transition flow regimes and saturates in the near-continuum flow regime (δ ≥ 10). Since, compared to the fast spectral approximation to the Boltzmann collision operator the time for solving Eq. (28) is negligible, the CPU time saving is proportional to the time-step saving, and this is tremendous for the SIS. At δ = 10, the SIS is about 15 times faster than the CIS, while at δ = 50 it is about 220 times faster.
It is interesting to note that for both the SIS and CIS, solutions of the LBE for hard-sphere molecules converge about 1.5 times faster than for Maxwell molecules, a result which supports the convergence analysis in Sec. 2.3.
We also consider Poiseuille/thermal transpiration flows along a channel of square cross section. Due to symmetry, only one quarter of the spatial domain is simulated, which is divided into 50 × 50 non-uniform cells: in each direction, from the boundary to the center, the length of each cell side forms a geometric progression with a common ratio 1.05. The diffusion-type equation (22) is discretized by a five-point central difference, and solved by the successive-over-relaxation method [37]. Table 2 summarizes the numerical results from the SIS. The mass flow rate in thermal transpiration flow is not shown, as according to the Onsager-Casimir relation it is equal to the heat flow rate in Poiseuille flow. From this table it is seen that our SIS for the LBE works efficiently over the whole range of gas rarefaction. This efficient SIS can also be used to calculate the slip coefficients. In Ref. [38] it was stated in that the thermal slip coefficient is strongly affected by the intermolecular potential. Therefore, we calculate this coefficient based on the LBE for hard-sphere and Maxwellian molecules, and compare our results to that of the Shakhov kinetic model [39]. Since in the near-continuum regime, the dimensionless mass flow rate in a thermal transpiration flow can be expressed as [40] we calculate the thermal slip coefficient as σ T = 2δM T by choosing large values of δ. The numerical results are shown in Table 3. It is clear that there is a large difference (nearly 20%) between the LBE for the hard-sphere gas and the Shakhov model, but the difference between the LBE for the Maxwellian gas and the Shakhov model is small. This is probably because the collision frequencies in the LBE for Maxwellian molecules and the Shakhov model are constants, while that in the LBE for hard-sphere molecules is a function of the molecular velocity. This example shows that it is necessary to use the Boltzmann equation even in the near-continuum flow regime, and the SIS developed here is useful, for example, in assessing the accuracy of various kinetic gas-surface boundary conditions by comparing the numerical solutions for thermal transpiration flow with experimental data [41,42].

The most efficient scheme
In Sec. 2.4 we saw that the SIS for the LBE can be faster than the CIS by several orders of magnitude in the nearcontinuum regime. Now we look at the possibility of speeding up the convergence even more, without modifying the SIS too much. To this end, we penalize the linearized Boltzmann collision operator into the following form: where the constant N is a tuning parameter which affects the convergence rate.
With Eq. (20) replaced by Eq. (30), the diffusion-type equation (22) should be changed accordingly. Our numerical results for Poiseuille flow between two parallel plates show that, for a fixed δ, all synthetic schemes with different values of N converge to the same solution (with relative errors in flow rates less than 0.1%). However, the convergence rate varies with N. From the top row in Fig. 2 we see that in the free molecular regime all schemes have the same convergence rate, while in the transition regime the scheme with N < 1 (N > 1) converges slower (faster) than that with N = 1. The situation becomes complicated in the near-continuum regime: for hard-sphere molecules, the case with N = 1.5 converges fastest, followed by N = 1, 2, and 0.5. For Maxwell molecules, however, the cases with N = 1.5 and N = 2 have roughly the same fast convergence, followed by N = 1 and 0.5. Similar behaviors are observed for the thermal transpiration flow.
To further investigate the relationship between the convergence iteration step and N in the synthetic scheme, we fix δ = 100 and vary N. The numerical results in the bottom row of Fig. 2 show that the fastest convergence is achieved when N is approximately 2 and 1.5 for Maxwell and hard-sphere molecules, respectively. This may be interpreted in terms of the average collision frequency. In the LBE, the equilibrium collision frequency ν eq is in general a function of the molecular velocity. The average collision frequency, varies between different intermolecular potentials even when the shear viscosity is the same. We found that for Maxwell and hard-sphere molecules, the average collision frequencies are 2.22 and 1.25 times larger than the rarefaction parameter, respectively-which are very close to the two values for the fastest convergence as seen in the bottom row of Fig. 2. Therefore, to achieve the best performance of the synthetic scheme we suggest using N =ν δ , (32) which further reduces the iteration number by about 30% when compared to the case of N = 1.

Further benefit in using SIS
In addition to the significant speed-up of convergence, the SIS can also help to reduce the spatial resolution. It is wellknown that in order to solve the kinetic equations the size of the spatial cells in the traditional discrete velocity method or the DSMC method should be smaller than the molecular mean free path, so that numerical results are reliable because the artificial viscosity is much smaller than the physical viscosity. In the SIS, the macroscopic flow velocity is obtained by solving the synthetic diffusion-type equation (22) which is asymptotic-preserving into the Navier-Stokes regime, so the spatial resolution can be relatively coarser.
To demonstrate this, we run the test case in Sec. 2.4 again, but the half spatial domain is instead divided into 10 uniform cells; for δ = 200, this means that the spatial cell size is about 10 times larger than the molecular mean free path. Both the SIS and CIS on this coarse grid are compared to the reference solutions of Table 1. Since the mass flow rate in Poiseuille flow increases rapidly with δ, we study how the apparent gas permeability changes with the rarefaction parameter. Here, the apparent gas permeability, which is normalized by 2 , is defined as Note that according to the Navier-Stokes equation with the no-slip velocity boundary condition, the flow velocity satisfies  Thus we have U 3 = −(δ/2)(x 2 1 − 1/4), and κ = 1/12; this permeability is also known as the "intrinsic" or "liquid" permeability. The apparent gas permeability is always larger than the intrinsic permeability, and increases with 1/δ (or the Knudsen number). Figure 3 shows the apparent gas permeability obtained for these different spatial resolutions over a wide range of the rarefaction parameter. It is clear that the SIS, even with a coarse spatial resolution, can yield good results. This proves that the SIS is asymptotic-preserving into the Navier-Stokes regime. The CIS results, however, have larger errors at large values of δ. For instance, when δ = 150, the non-accelerated scheme underpredicts the apparent gas permeability by about 12.5%. This error continues to increases with δ: we have tested the BGK model for δ = 10 4 and found that the relative error is about 62.5% [43].

SIS in multiscale problems
We now investigate the performance of the SIS in more complex geometries, where the problems are multiscale in the sense that the rarefaction parameter varies by several orders of magnitude due to different characteristic flow length scales, e.g. flow in a fractal geometry.

Rarefied gas flow through the Sierpinski carpet
We first consider the gas flow through a two-dimensional cross section described by the Sierpinski carpet, which can be generated through recursion. Beginning with a square, the square is cut into 9 congruent subsquares in a 3 × 3 grid, and the central subsquare is removed. The same procedure is then applied recursively to the remaining 8 subsquares. Resulting geometries after several levels of recursion are presented in Fig. 4.
Due to the symmetry, the one quarter of the level 1, 2, and 3 Sierpinski carpets are divided into 60 × 60, 90 × 90, and 135 × 135 uniform cells, respectively. The molecular velocity space is represented by 32 × 32 × 12 discrete grids. We find that the iteration number, using N = 1.5 in Eq. (30) for a hard-sphere gas, is always less than 55 for each rarefaction parameter, when the relative error in the mass flow rate M P between two consecutive iterations is less than 10 −5 .   Figure 6 shows the mass and heat flow rates in the Poiseuille flow of the hard-sphere gas through the Sierpinski carpet.
The Knudsen minimum in the mass flow rate can be seen, however, the location of the minimum M P shifts towards larger values of δ as the recursion level of the Sierpinski carpet increases. This is because, in the calculation of δ according to Eq. (12), the characteristic flow length is chosen to be the side length of the largest square, which is larger than, say, the smallest side length of the solids near which the flow velocity is maximum. As the recursion level increases, the porosity (the void area fraction) of the Sierpinski carpet decreases, so the flow rates decrease. We also plot in Fig. 6 the thermomolecular pressure difference exponent, which is an important parameter determining the performance of a Knudsen pump. The exponent always increases with decreasing δ and approaches the value of 0.5 when δ → 0 if the diffuse gas-surface boundary condition is used [40]. This also indicates the correctness of our numerical simulations.

Rarefied gas flow through random structures
We also consider the gas flow through two-dimensional porous media, where the porosity is 0.6. The first porous medium is generated by adding circular solids of different radii randomly to a square. The radius ratio of the largest disc to the smallest is 10. The square is then divided into 200 × 200 uniform cells, and the discs are approximated by the "staircase", as visualized in Fig. 7(a). The second porous medium, shown in Fig. 7(b), also consisting of 200 × 200 uniform cells, is generated by the quartet structure generation set [44].
Following the numerical simulations, velocity contours are displayed in Fig. 8 for the free molecular and near-continuum flow regimes, while the mass and heat flow rates are shown in Fig. 9 from the free molecular to near-continuum flow regimes. It is interesting to note that, in the free molecular flow regime, the mass flow rates in the two random porous media are nearly the same. However, in the near-continuum regime, the mass flow rate of the porous medium consisting of random squares is about twice that in disc medium. This research may find applications in shale gas extraction.

A special case: SIS in polar coordinates
The SIS developed above for the LBE works well in Cartesian coordinates, for rarefied gas flows through general cross sections. For flows through circular cross sections, the use of polar coordinates can reduce the computational cost significantly, as previously demonstrated in the SIS for the McCormack kinetic model [22,24]. However, for the LBE, the velocity space cannot be represented in polar coordinates as that of the McCormack kinetic model, since the molecular velocity along the flow direction cannot be integrated due to the nonlinear structure of the Boltzmann equation. Therefore, the SIS for the  LBE will be developed in polar coordinates for spatial variables, while the three-dimensional molecular velocity space is represented by cylindrical coordinates.
We consider Poiseuille flow along a pipe as an example, where the axis of the pipe is along the x 3 direction, and its cross section is located in the x 1 -x 2 plane. The spatial coordinates are normalized by the radius of the tube. Introducing the transformation x 1 = r cos θ , x 2 = r sin θ , v 1 = v r cos θ , v 2 = v r sin θ , and defining the VDF h = h(r, θ, v r , v 3 ) in cylindrical (molecular velocity)-polar (space) coordinates v r ∈ [0, +∞), θ ∈ [0, 2π ], v z ∈ (−∞, ∞), and r ∈ [0, 1], the LBE can be written   Fig. 7 consisting of random discs (circles) and islands (squares).
by taking the velocity moment of the LBE (35). Multiplying Eq. (35) by 2v 3 and integrating over the molecular velocity space, we obtain the equation for the shear Multiplying Eq. (35) by v 3 v 1 , penalizing the linearized Boltzmann collision operator in the form of Eq. (30), and integrating over the molecular velocity space, we obtain ∂ ∂rˆv which is simplified, with the help of Eq. (37), into 1 r ∂ ∂r If we express ´2 v 3 v 2 1 hdv =´v 3 (2v 2 1 − 1)hdv + U 3 , the diffusion-type equation in the form of Eq. (36) can be derived. But for practical numerical calculations, the following first-order ordinary differential equation for the flow velocity is desirable: In the numerical simulation, the spatial coordinate r is discretized into 150 nonuniform points, with most of the points located near the pipe surface r = 1. Due to symmetry, the truncated velocity v r ∈ (0, 4) is discretized into 64 nonuniform points, with most of the points located near v r = 0, while θ ∈ [0, π] and v 3 ∈ (0, 6) are discretized into 40 and 12 uniform points, respectively. The linearized Boltzmann collision operator is approximated by the fast spectral method: first, the spectrum of the VDF is calculated by Fourier transform from the cylindrical molecular velocity space to the Cartesian frequency space. Second, the fast spectral method [28] is applied to find the spectrum of the linearized Boltzmann collision operator in the Cartesian coordinate. Finally, the inverse Fourier transform is used to find the collision operator in the cylindrical space. The SIS in polar coordinate is implemented in the following three steps: first, as usual, Eq. (35) is solved by the implicit iterative scheme, with the spatial derivatives being approximated by the second-order upwind finite difference. Then, Eq. (40), which is used to expedite convergence to the steady-state solution, is also solved using the second-order upwind finite difference, where the boundary condition of U 3 (r) at r = 1 is calculated from the VDF obtained in the first step. Finally, having obtained U 3 (r) from Eq. (40), a correction in the VDF is performed, see Eq. (24).
Our numerical results for mass and heat flow rates in Poiseuille flow of a hard-sphere gas through a pipe are summarized in Table 4. It can be seen that, as for rarefied gas flows between two parallel plates and in rectangular cross sections, the SIS in the polar coordinates is also very efficient.

SIS for the LBE for gas mixtures
In this section we develop a SIS for the LBE for binary gas mixtures. For simplicity, only Poiseuille flow is considered, but the method can be generalized to flows driven by temperature and concentration gradients.
Let f A and f B be, respectively, the VDFs of gas components A and B with molecular masses m A and m B , and molar fractions χ A and χ B = 1 − χ A . Introducing the equilibrium VDF (in which the velocity is normalized by the most probable and expressing the VDF in the form f α = f α,eq + h α , where h α are perturbed VDFs satisfying |h α / f α,eq | 1, the LBE for h α with the linearized Boltzmann collision operators L α = β=1,2 Q αβ ( f α,eq , h β ) + Q αβ (h α , f β,eq ), where the details of Q αβ can be found in Ref. [35]. The source term for Poiseuille flow is and again we take X P = −1.
When the perturbed VDFs are known, the flow velocity normalized by v mA is calculated as U α =´h α v 3 dv/χ α , shear stresses normalized by the total gas pressure p 0 are P α13 = 2m α´hα v 1 v 3 dv/m A and P α23 = 2m α´hα v 2 v 3 dv/m A , and the heat flux normalized by p 0 v mA is q α =´ m α |v| 2 /m A − 5/2 v 3 hdv. The dimensionless mass flow rate M and heat flow rate Q, normalized by the most probable speed of the gas mixture, are calculated as where m = χ A m A + (1 − χ A )m B is the average molecular mass of the mixture.

The synthetic scheme for a gas mixture
As emphasized above, the relation ∂ P 13 /∂x 1 + ∂ P 23 /∂x 2 = 1 is important in developing the SIS. For binary mixtures, this relation still holds, but now shear stresses are replaced by mixture shear stresses, i.e. P 13 = P A13 + P B13 and P 23 = P A23 + P B23 . This poses an additional difficulty.
Following the basic steps in developing the synthetic equation, we obtain the following two equations for the flow velocity of each component: To obtain diffusion-type equations which recover the Stokes equation (an important step, guaranteeing fast convergence) in the hydrodynamic regime, we rewrite the linearized collision operator as 2 : where N α are two constants, and let 2 while for component A, by using the relation ∂ P 13 /∂x 1 + ∂ P 23 /∂x 2 = 1 with P 13 = P A13 + P B13 and P 23 = P A23 + P B23 , We therefore propose the SIS for the LBE for binary gas mixtures: while the VDFs in Eq. (42) are first solved by the CIS [35], flow velocities are updated according to diffusion-type equations (47) and (48). Then the VDFs are corrected in a way similar to Eq. (24). Note that the fastest convergence is achieved when where ν α,eq is the equilibrium collision frequency of the α-component.  [25] for multiple gas mixture, by simply replacing L α in Eqs. (47) and (48) with that in the McCormack model; the resulting diffusion-type equations will be much simpler than those in Refs. [37,45].

Convergence analysis
To show the efficiency of the proposed SIS for binary gas mixtures, we calculate the eigenvalue of the iteration as a function of the inverse rarefaction parameter. The numerical procedure is essentially the same as that in Sec. 2.3 for single-species gases. Fig. 10 shows the eigenvalue of both the SIS and CIS for an equimolar Ne-Ar mixture, where Ne and Ar are treated as hard-sphere molecules with a molecular diameter ratio of 0.711. It is clear that the SIS is superior to the CIS at large values of the rarefaction parameter. Also, when compared to the SIS for a single-species hard-sphere gas, the synthetic scheme for a binary gas mixture has roughly the same maximum eigenvalue, i.e. ω 0.6. So it is expected that the synthetic scheme for a binary gas mixture will be as efficient as that for a single-species gas. It is also interesting to note that at small values of δ the eigenvalue of the SIS is slightly higher than that of the CIS; the reason for this is not clear. However, this does not affect the efficiency of the SIS because at small values of δ both the SIS and CIS converge rapidly, i.e. within a small number of iterations.

Numerical simulations of Poiseuille flow
Following the convergence analysis above, numerical simulations are now conducted for the Poiseuille flow of a Ne-Ar mixture between two parallel plates. First, the hard-sphere molecular potential is considered. Table 5 shows the mass and Table 5 Mass and heat flow rates in Poiseuille flow of Ne-Ar gas mixtures between two parallel plates, as well as the number of iterations (Itr) to reach the convergence criterion = 10 −10 in the SIS. The LBE with the hard-sphere model is used.  In the free molecular and transition regimes, the mass flow rate of the lighter species is always higher than that of the heavier species; however, in the near-continuum regime the mass flow rates are the same for both gas components. When compared to the results in Table 1 we see that the difference in the mass flow rates is within 1% in the near-continuum regime. A comparison of velocity and heat flux profiles in different flow regimes is presented in Fig. 11, which shows that the difference between the velocity profiles of each component is reduced as the rarefaction parameter increases.
As δ increases, the difference between the heat fluxes further normalized by the molar fraction of each component also decreases, but in the near-continuum limit the heat flux of Ne is roughly √ 2 times that of Ar, while the difference in the mass flow rate goes to zero.
We have also simulated for the first time the Ne-Ar mixture flow based on the LBE for Lennard-Jones potentials, where the fast spectral approximation of the Boltzmann collision operator is described in Ref. [30]. The influence of the molecular model on the velocity and the heat flux profiles is shown in Fig. 11: the hard-sphere model overestimates the velocity and the heat flux in the free-molecular regime, but in the transition and near-continuum regimes the differences between the two molecular models reduce as δ increases. This is in good agreement with observations of the single-species case [29,30]. Figure 12 shows the mass and heat flow rates as a function of the rarefaction parameter. The two flow rates obtained from the LBE with Lennard-Jones potentials are both smaller than that for the hard-sphere model when δ < 1 (e.g. when

Table 6
Mass and heat flow rates in Poiseuille flow of an equimolar Ne-Ar gas mixture along channels of rectangular cross section, of aspect ratios one and two, as well as the number of iterations (Itr) needed to reach the convergence criterion = 10 −7 in the SIS. N α = 1.5 is adopted in Eqs. (47)  Finally we calculate the Poiseuille flow of an equimolar Ne-Ar mixture along channels of rectangular cross section, based on the hard-sphere model. To the best of our knowledge, the LBE for a gas mixture has not previously been solved in a two-dimensional geometry, because of the numerical complexity; we tackle the problem here by using the SIS and the fast spectral method [35]. The discretization of the spatial domain for a square cross section is the same as that in Sec. 2.4, while for a rectangular cross section of aspect ratio 2 the spatial domain is discretized by 50 × 100 cells, and the characteristic length is chosen to be the shorter side. Table 6 summarizes the LBE solution for the mass and heat flow rates in two-dimensional Poiseuille flow over a wide range of the rarefaction parameter. The iterations needed to achieve the convergence criterion = 10 −7 are fewer than 40 when δ ≤ 100, demonstrating again the efficiency of the SIS.
The normalized mass flow rates of Ne and Ar through the rectangular cross section with aspect ratio 2 are always larger than those in the aspect ratio 1 case. However, although the heat flow rates through the rectangular cross section with an aspect ratio of 2 are larger than those with the aspect ratio 1 when δ < 20, they are roughly the same for δ ≥ 20.
Typical velocity and heat flux profiles of the Ne-Ar gas mixture through the square cross section in the free molecular, transition, and near-continuum regimes are shown in Fig. 13, and are also compared to those of the single-species gas in the same geometry. From this figure we see that the velocity and heat flux of Ne is always larger than that of Ar, while the corresponding results for a single-species gas lie between those of Ne and Ar.

Conclusions
We have proposed a synthetic iterative scheme to accelerate the convergence of the linearized Boltzmann equation for gas flows driven by pressure and temperature gradients in long channels. By penalizing the linearized Boltzmann collision operator L into the form L = (L + N L BG K ) − N L BG K or L = (L + Nδh) − Nδh, a diffusion-type equation has been derived for the macroscopic flow velocity. The velocity distribution function in the linearized Boltzmann equation was first solved by the conventional iterative scheme, where the linearized Boltzmann collision operator was approximated by the fast spectral method. Then the flow velocity was obtained by solving the diffusion-type equation, which was used to correct the velocity distribution function. In this way the slow convergence of the conventional iterative scheme in the near-continuum flow regime has been ameliorated: we found, through the numerical solution of Poiseuille and thermal transpiration flows, that the synthetic iterative scheme is faster than the conventional scheme by up to several orders of magnitude. More, it has been found that the synthetic iterative scheme is asymptotic-preserving into the Navier-Stokes level, so in the near-continuum regime the spatial resolution can be much larger than the molecular mean free path. This makes the synthetic iterative scheme much more faster and accurate than the conventional one, especially in multiscale problems with a wide variation of local flow length.
The tuning parameter N controls the convergence rate of the synthetic iterative scheme. In numerical investigations we found that for the linearized Boltzmann equation the fastest convergence is achieved when N roughly equals the ratio of the equilibrium collision frequency to the rarefaction parameter. Thus, N varies with the intermolecular potential: for a single-species gas, we found that the fastest convergence occurred with N ≈ 1.5 for the hard-sphere gas model and N ≈ 2 for the Maxwell gas model.
We also extended the synthetic iterative scheme to binary gas mixtures, and both the hard-sphere and Lennard-Jones intermolecular potentials have been considered. As an example, Poiseuille flow of a Ne-Ar mixture was simulated in order to test the computational performance as well as the influence of the intermolecular potential. The synthetic iterative scheme required only a limited number of iterations over the whole range of gas rarefaction. Based on this efficient scheme, the Poiseuille flow of a Ne-Ar mixture between two parallel plates was simulated for the first time using the realistic Lennard-Jones potential. We found that the hard-sphere gas model overestimates the mass and heat flow rates when δ < 1. Poiseuille flow of a Ne-Ar mixture through two-dimensional rectangular cross sections was also simulated using the linearized Boltzmann equation for the first time.
The developed method can also be extended to the efficient calculation of flows of multiple-species gas mixtures. In particular, our method can be applied to the McCormack model and we believe that the resulting diffusion-type equations for the flow velocity of each component will be much simpler than those in Refs. [18,23]. Our synthetic iterative scheme could also be applied straightforwardly to other canonical gas flows, such as Couette flow and the flow driven by a concentration gradient. However, it requires future work to investigate whether this method can be applied to other gas flow systems (such as the cavity flow) or not.