Fokker–Planck–DSMC algorithm for simulations of rarefied gas flows

Journal of Computational Physics 287 (2015) 110–129

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Fokker–Planck–DSMC algorithm for simulations of rarefied gas flows M. Hossein Gorji ∗ , Patrick Jenny Institute of Fluid Dynamics, ETH Zentrum, Sonneggstrasse 3, 8092 Zürich, Switzerland

a r t i c l e

i n f o

Article history: Received 24 May 2014 Received in revised form 28 January 2015 Accepted 29 January 2015 Available online 4 February 2015 Keywords: DSMC Fokker–Planck equation Rarefied gas flows Kinetic models

a b s t r a c t A Fokker–Planck based particle Monte Carlo algorithm was devised recently for simulations of rarefied gas flows by the authors [1–3]. The main motivation behind the Fokker– Planck (FP) model is computational efficiency, which could be gained due to the fact that the resulting stochastic processes are continuous in velocity space. This property of the model leads to simulations where the computational cost becomes independent of the Knudsen number (Kn) [3]. However, the Fokker–Planck model which can be seen as a diffusion approximation of the Boltzmann equation, becomes less accurate as Kn increases. In this study we propose a hybrid Fokker–Planck–Direct Simulation Monte Carlo (FP–DSMC) solution method, which is applicable for the whole range of Kn. The objective of this algorithm is to retain the efficiency of the FP scheme at low Kn (Kn  1) and to employ conventional DSMC at high Kn (Kn  1). Since the computational particles employed by the FP model represent the same data as in DSMC, the coupling between the two methods is straightforward. The new ingredient is a switching criterion which would ideally result in a hybrid scheme with the efficiency of the FP method and the accuracy of DSMC for the whole Kn-range. Here, we adopt the number of collisions in a given computational cell and for a given time step size as a decision criterion in order to switch between the FP model and DSMC. For assessment of the hybrid algorithm, different test cases including flow impingement and flow expansion through a slit were studied. Both accuracy and efficiency of the model are shown to be excellent for the presented test cases. © 2015 Elsevier Inc. All rights reserved.

1. Introduction It is well established that the Navier–Stokes–Fourier description of gas flows breaks down as the Knudsen number becomes relatively large. Alternatively, given the dilute gas assumption, the Boltzmann equation is valid for the whole Knudsen number range. However, due to the complexity of the Boltzmann collision operator and due to the high dimensionality of the solution domain, different modeling attempts have been made in order to simplify the governing equation. The most widely used simulation method based on the Boltzmann equation is DSMC, which was originally proposed by Bird [4]. There, computational particles mimic the description by the Boltzmann equation in two steps, i.e. streaming and collisions. DSMC has been validated and enhanced for complex gas flows (including polyatomic effects, mixtures, chemical reactions, etc.) over the past decades [5]. Despite the huge success of DSMC, it is a rather expensive method for near continuum and low Mach flow scenarios. In the former case, which occurs at low Knudsen numbers, the cost is hight due

*

Corresponding author. E-mail address: [email protected] (M.H. Gorji).

http://dx.doi.org/10.1016/j.jcp.2015.01.041 0021-9991/© 2015 Elsevier Inc. All rights reserved.

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

111

to the spatio-temporal discretization in DSMC. In the latter case, the statistical noise becomes large, which requires many particles. In this paper, the focus lies on the issue related to low Kn. We show how an FP based particle algorithm can fill the gap between continuum flows and DSMC, thus extending the applicability of DSMC to problems involving not only large, but also very small Knudsen numbers. Note that certain modifications already exist for DSMC which address the near continuum problem [10,11]. Also in order to control the statistical error, the deviational approach (see e.g. [9]), the moment guided approach [12] and the information preserving approach [8] have been proposed. Here, we consider conventional DSMC as discussed in [5]. The Fokker–Planck approximation of the Boltzmann equation can be used for rarefied gas flows as shown in [1–3,24,25]. The objective behind this model is the numerical advantage, which can be gained due to the fact that the resulting diffusion model is computationally less challenging than the Boltzmann collision integral. A particle Monte Carlo scheme based on the Fokker–Planck equation was proposed in [1] and was later extended to the cubic model to honor the proper Prandtl number for monatomic gas [2,3]. The Fokker–Planck model then has been generalized in [24] and [25] for diatomic molecules and mixtures, respectively. It is worthwhile to mention that other research groups have also developed solution algorithms based on the Fokker–Planck model. For example, Bogomolov and co-workers developed a Fokker–Planck approximation based on hard sphere molecules [21,22]. Another development was a Fokker–Planck based information preserving algorithm in [28] for low Knudsen numbers. Different hybrid methods already exist to couple DSMC and continuum solvers for the Navier–Stokes–Fourier system (see e.g. [13–15]). In general, challenges have to be addressed during the development of such hybrid DSMC–CFD solution algorithms. First, since CFD is only applicable for the continuum regime, a careful positioning of the interface between the two methods should be made. Second, exchanging information between the two methods is not straightforward. However note that the continuum solvers are typically much more efficient than stochastic particle ones due to the fact that the results are free from the statistical noise and the computation only scales with the number of computational cells. A more natural coupling can be achieved, if DSMC is combined with another particle based kinetic algorithm; e.g. BGK based particle method (see e.g. [16]). In comparison to the BGK particle method, the algorithm presented here has some computational advantages. First, BGK based particle simulations with correct Prandtl number (i.e. ellipsoidal BGK [7] or S-BGK [6]) employ acceptance-rejection to sample particle velocities, which can become very expensive. In the FP algorithm, however, the required random numbers can be sampled from a normal distribution, which are generated very efficiently (e.g. using Ziggurat algorithm). Second, the BGK equation similar to the Boltzmann equation, describes a jump type process for the velocities. This leads to restrictions for grid spacing and time step size similar to DSMC. On the other hand, the FP approach can employ cell and time step sizes which are larger than the collisional scales [1–3]. Another possible hybrid approach is given by coupling DSMC with low-diffusion equilibrium or near-equilibrium particle methods [17–20]. Note that the former, i.e. equilibrium particle methods employ Maxwellian distribution for the velocity of particles and thus they are mostly suitable for inviscid gas flow simulations. The extension of the low-diffusion particle method to the near-equilibrium condition was provided in [19], where Stokes and Fourier laws were employed in order to model the effects of viscous stresses and heat fluxes on the bulk velocity and temperature of the gas. Further, by providing slip velocity and temperature jump for the boundary conditions, more accurate results compared to the equilibrium particle methods could be obtained [19]. The cubic FP model considered in this article however, provides an evolution equation for the distribution function. Thus no closure assumption for heat fluxes or viscous stresses is required here. The model is constructed such that it gives rise to correct viscosity and heat conductivity coefficients at the hydrodynamic limit, as well as accurate description of macroscopic flow properties for the non-equilibrium condition e.g. counter Fourier heat fluxes for the lid driven cavity flow [3]. Furthermore, similar to DSMC, no boundary condition for macroscopic flow properties are required by the FP model. Therefore in comparison to low diffusion particle methods, the FP model provides more physics in non-equilibrium flow conditions and thus more accurate hybrid simulations may become possible. Motivated by the fact that DSMC becomes less efficient at low Knudsen numbers, and that the cubic Fokker–Planck model is very accurate up to quite large Knudsen numbers [2,3], a hybrid Fokker–Planck–DSMC approach can be very attractive. Note that the cost of the FP solution is independent of the Knudsen number, since no collisions are computed. Moreover, time steps larger than the mean collision time and grid spacings larger than the mean free path are allowed [1,3]. At the same time, accurate descriptions of rarefied gas flow can be obtained for a considerable range of Knudsen numbers. For example counter Fourier heat fluxes in lid-driven cavity flows at Kn ≈ 0.05 [3] or accurate mass flow rates in Poiseuille flow up to Kn ≈ 5 [2] can be computed with the cubic FP model. Therefore, one can rely on the FP method for low Knudsen numbers and due to the similar description of particles in both DSMC and the FP models, a hybrid scheme can be constructed without any modification of either methods. The rigorous development of a hybrid Fokker–Planck DSMC method is the objective of this paper with the following structure. First, a short review of the cubic Fokker–Planck model for monatomic gas is presented. Second, it is explained how the model can be coupled to DSMC via a justifiable criterion. Section three presents some validations for homogeneous relaxation, one dimensional Couette flow, one dimensional shock tube, flow through a thin slit and plume–wall interaction. It is shown that the hybrid scheme can lead to very accurate result at a computational cost which is substantially lower than for full DSMC. The paper is concluded with final remarks and the future outlook.

112

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

2. Review of the Fokker–Planck kinetic model In this section, the cubic Fokker–Planck model as an approximation of the Boltzmann equation is briefly reviewed. The model was originally devised in [2]. Consider an ensemble of N gas molecules, each with a random velocity M i (t ) and a random position X i (t ) at time t, where i ∈ {1, 2, . . . , N }. The one point statistical description of the flow is fully determined by the probability density function (PDF) f ( V , x; t ), which is a function of the velocity sample space V , position sample space x and time t. The relation between molecular random variables and the PDF is

f ( V , x; t ) = lim

N →∞

N 1 

N

δ( M i (t ) − V )δ( X i (t ) − x),

(1)

i =1

where δ() is the Dirac delta. The PDF f ( V , x; t ) can be converted into the velocity PDF f V ( V ; x, t ) using the identity

f V ( V ; x, t ) =

f ( V , x; t ) f X ( x; t )

(2)

, 

where the marginal PDF f X (x; t ) = R3 f ( V , x; t )dV 1 dV 2 dV 3 is introduced. It can be shown that f X = ρ (x, t )/M, where is the gas density and M is the total mass in the sample space x (see [24] for derivation). The mass density function (or more precisely the mass weighted probability density function) is defined as

ρ

F ( V , x, t ) = M f ( V , x; t )

= ρ (x, t ) f V ( V ; x, t ). Under the dilute gas and molecular chaos assumptions, the Boltzmann collision operator



∂F ∂t

 = coll

1

 4π 

m

 F ( V )F ( V t ) − F ( V )F ( V t ) g I (, g )ddV t 1 dV t 2 dV t 3

(3)

R3 0

can be derived for the rate of change of F due to the binary collisions between gas molecules. Here, g is the magnitude of the relative velocity g = V − V t , I (, g ) is the differential cross section,  is the solid angle about g, and m is the molecular mass. Note that the velocity pairs ( V , V t ) are post collision velocities of the pair ( V , V t ). If M (t ) can be approximated by a continuous Markovian process, then the temporal evolution of F can be described by a Fokker–Planck equation [1–3,29], i.e.



∂F ∂t



≈− coll

∂ ∂Vi

( Ai F ) +

 ∂2  DijF . 2 ∂ Vi∂ V j 1

(4)

Here, A and D represent drift coefficient and the positive definite diffusion, which in general are functions of the MDF. Here and henceforth Einstein’s summation convention is adopted. Our motivation for developing a Fokker–Planck approximation was computational efficiency. Since here the underlining random process is continuous in velocity space, larger time steps and larger computational cells can be employed than in the case of solving Eq. (3). In the following, we review the closure of the drift and diffusion coefficients based on the cubic model of [2], which leads to correct viscosity and heat conductivity in the Navier–Stokes–Fourier limit. 2.1. Cubic model The closures for drift and diffusion coefficients can be found based on the comparison between the Fokker–Planck and the Boltzmann operators. Let’s consider the production terms resulting from the Boltzmann collision integral, i.e.



φ S Boltz dV 1 dV 2 dV 3 ,

(5)

R3

where φ ∈ { V i , V i V j  , V i V j V j } and S Boltz is the right hand side of Eq. (3). Note that the weights φ correspond to the mean velocity U , molecular stresses π and the heat flux vector q. Assuming Maxwell type molecules, the production terms yield



V i S Boltz dV 1 dV 2 dV 3 = 0

(6)

R3



V i V j  S Boltz dV 1 dV 2 dV 3 = −

p

μ

πi j

and

(7)

R3



V i V j V j S Boltz dV 1 dV 2 dV 3 = − R3

2p 3μ

qi ,

(8)

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

113

where p is the pressure and μ the viscosity. Now, by setting the production terms resulting from the Fokker–Planck operator equal to Eqs. (6)–(8) and fulfilling the conservation laws, the model coefficients can be derived. The final form of the closure reads

Ai = − and

Dij =

1 v i + c i j v j + γi



τ

2kT

v j v j −

3kT



  2qi + v i v j v j −

(9)

ρ

m

(10)

δi j ,

mτ

where c i j , γi and are the macroscopic coefficients and the fluctuating velocity vector is v = V − U with U being the bulk velocity. Here, the Fokker–Planck time constant is τ = μ/(2p ), T is the temperature, k is the Boltzmann constant and δi j the Kronecker delta. Note that the coefficients c i j , γi and depend on the moments of the distribution function (see Appendix A for corresponding relations). 2.2. Particle system The particle representation of the Fokker–Planck solution can be devised via Ito¯ processes. Let’s consider an ensemble of particles with positions X i (t ) and velocities M i (t ). Similar as in DSMC, each particle represents a certain number of real molecules and therefore, by introducing the positive statistical weights w i , one can find the associated mass density function M



 N 

F ( V , x, t ) = lim

N →∞

w

i =1

i

f ( V ,x;t )

N

i =1 δ( X

lim

N →∞

i

 (t ) − x)δ( M i (t ) − V ) w i .

N i i =1

(11)

w

˜ } which are inside an infinitesimal volume about position x, i.e. with X i (t ) ∈ Note that if only those particles i ∈ {1, . . . , N δ(x), are considered, then Eq. (11) simplifies to ρ (x,t )

f V ( V ;x,t )





N˜ ˜N i i  1 i =1 δ( M (t ) − V ) w F ( V , x, t ) = wi ,

˜ N δ(x) wi



i =1

=

1

δ(x)

˜ N 

(12)

i =1

δ( M i (t ) − V ) w i ,

(13)

i =1

which is convenient in particle simulations in order to calculate macroscopic moments. Further, let’s assume that F fulfills the Fokker–Planck equation

 ∂2  ∂ ∂F ∂F ∂F + Gi + Vi =− DijF , ( Ai F ) + ∂t ∂Vi ∂ xi ∂Vi ∂ Vi∂ V j

(14)

where G is an external force. According to [3], the corresponding random velocities and positions should follow the Ito¯ processes

dM i = A i dt + G i dt + dik dW k

and

(15)

d X i = M i dt ,

(16)

where dW k are independent Wiener increments with zero expectation and variance dt (see [30] for details on Wiener process). The diffusion coefficient is obtained from the relation dik dkj = D i j . 2.3. Time integration Solutions of the stochastic differential equations (15)–(16) can be computed by numerical integration. For the studies in this paper, the integration presented in [3] was employed. Here only its final form is presented. Let M ni be an approximation of M i (t n ). The velocity update for a time step size of t = t n+1 − t n then reads

M ni +1

=

1

α

where

α2 = 1 +





n −t /τ

Mi e

mτ  3kT

+ (1 − e

−t /τ

)τ

N ni

+



kT m



(1 − e −2t /τ )ξ 

i

+ U in + G i t ,

(17)



τ (1 − e−t /τ )2 N ni N ni  + 2 e−t /τ − e−2t /τ M i n N ni 

(18)

114

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

is the correction term, which ensures the energy conservation, and

N i = c i j M j + γi



M j M j −

3k m



T

  qi + M i M j M j −

ρ

(19)

is the nonlinear part of the drift term. M = M − U is the peculiar particle velocity, ξi are independent standard normal variables and . is the averaging operator. Similar to DSMC, the particle positions are updated as

X in+1 = X in + M ni t + G i

t 2 2

.

(20)

Note that more accurate position schemes can be constructed by integration of M (t ) similar to [1,2]. However, here due to the computational efficiency we use a simpler position updating scheme (20), which resembles DSMC. 2.4. Discussion Besides the computational advantages of the FP based solution algorithm at low Kn, the relation between the Fokker– Planck approximation and the Boltzmann equation, and the relation between the Fokker–Planck approximation and the Navier–Stokes–Fourier system were investigated. 1. The mathematical consistency between the Fokker–Planck and Boltzmann equations can be shown via the Kramers– Moyal expansion as carried out by Pawula [29] or by directly comparing the SDE (15) with respect to the Poisson random process based on the Boltzmann equation [21]. The first approach considers the fact that if the particle velocity evolution can be approximated by a continuous Markov process, then one can apply the Kramers–Moyal expansion for the temporal evolution of F . Finally, the theorem of Pawula leads to the Fokker–Planck equation. The second approach based on SDEs was followed by Skorokhod [23], and later by Bogomolov and co-workers [21]. One can show that by using the central limit theorem, if the density of the Poisson operator becomes large, i.e. if the Knudsen number becomes small enough, the Poisson process converges to the diffusion process. Thus the evolution of F can be described by a Fokker–Planck equation. The error introduced by the Fokker–Planck approximation of the Boltzmann collision operator depends on F , which is unknown a priori. Therefore, one has to follow an empirical approach in demonstrating the accuracy of the Fokker– Planck model for different flow scenarios (as was carried out in various studies by the authors [1–3,24,25]). Furthermore, note that the resulting error of this asymptotic approximation of velocity jumps by a set of continuous SDEs, scales inversely with the number of collisions that each particle experiences in a given time interval. This scaling behavior of the error bound is the direct result of the Berry–Esseen inequality for approximating a scaled sum of random variables with the Gaussian distribution (see [21,23,26]). Consequently the error of the Fokker–Planck approximation of the Boltzmann operator, scales approximately with 1/( N coll / N p ), where N coll is the number of total collisions for a given t and N p number of particles. This measure, i.e. N coll / N p will be used as the basis of the hybrid Fokker–Planck–DSMC approach presented in this article. 2. The transport properties in the Navier–Stokes–Fourier limit can be derived from the Chapmann–Enskog expansion [27]. They are related to the moments of the collision operator. By construction, one can set the Fokker–Planck moments equal to the Boltzmann moments. The important physical moments are viscous stresses and heat fluxes, which correspond to the second and third velocity moments. For the cubic model used here (derived in [2]), consistent relaxations of the viscous stresses and heat fluxes are honored. Therefore both viscosity and Prandtl number can be fixed at the same time. 2.5. Comparison with DSMC The conventional DSMC algorithm of Bird [5] is based on a particle representation of the Boltzmann equation, in which the particles evolve according to a two step scheme. First, streaming of the particles according to the transport term in the Boltzmann equation and second, collisions between particles according to the Boltzmann collision integral are computed. In the Fokker–Planck method the particles represent the same mathematical object, i.e. mass density function, their evolution is somewhat different however. The effect of collisions is treated by a continuous stochastic process, i.e. by Eq. (15). This replacement leads to certain computational advantages. Foremost, collisions between particles are removed and second, large time steps and grid spacings then can be employed. Obviously the latter advantage relies on an accurate time integration of Eq. (15). In practice, for low Knudsen numbers the computational cost of simulations based on the FP scheme is shown to be much lower than that of DSMC [3]. It has to be emphasized however, such a treatment of the collisions by a continuous stochastic process is only justified as long as the Knudsen number is not too large (i.e. up to order 1). 3. Hybrid solution algorithm Very often the Knudsen number varies significantly inside the computational domain. This, together with the fact that DSMC becomes very expensive for the continuum or near continuum flow scenarios, is a very good motivation for hybrid

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

115

methods, in which different models and numerical schemes are applied for different Knudsen regimes. By coupling DSMC with the Fokker–Planck scheme the whole Knudsen number range can be covered by an accurate solution algorithm. Since the computational particles represent the mass density function in both DSMC and the FP simulations, coupling can be achieved without introducing interfaces. In the hybrid algorithm either DSMC or Fokker–Planck velocity-position updates are performed depending on a criterion which is individually applied in each computational cell, such that both efficiency of the FP method and accuracy of DSMC are optimally exploited. It has been shown in various studies that the FP method can safely be used for Knudsen numbers up to order one [1–3]. At the same time, DSMC becomes less efficient for Kn < O (1). Next, the computational efficiency of DSMC is analyzed more carefully. The local Knudsen number can be defined based on the gradient of a relevant macroscopic quantity Q

K nl ≈ λ|∇ Q / Q |,

(21)

where λ is the local mean free path. Note that the required grid spacing to adequately resolve variations of Q is

≤

1

(22)

.

|∇ Q / Q |

In DSMC however, also the mean free path has to be resolved leading to the more constringent requirement

≤λ

(23)

If  is determined based on Eq. (22), then the number of collisions (assuming that a CFL condition is consulted to calculate the time step size) may become larger than the number of particles in the cell. In the proposed hybrid algorithm,  is selected based on condition (22) irrespective of condition (23), and the time step size will be evaluated based on the CFL condition

t ≤

 

(24)

|U | + kT /m

which means that the resolution of both space and time are chosen independent of λ. The cost of DSMC scales with the number of collisions N ccoll . For a given cell one obtains

N ccoll ∝

Nc K nc

(25)

,

where N c denotes the number of particles in that cell. Since the number of collisions in DSMC depends on the distribution function which is unknown in general, for the following analysis, we focus on hard sphere type molecules in equilibrium. In that case the total number of collisions for a given time step size t becomes

N ccoll



√

2

= 2 π d nN c t

kT m

(26)

,

where d is the molecular diameter according to the hard sphere model. Based on the previous CFL criterion, we have

t ≤

 λ   = . |U | + kT /m K nc (|U | + kT /m)

Therefore the inequality

N ccoll

√

2

≤ 2 π d nN 

(27)



kT /m

 |U | + kT /m

(28)

should hold, and thus

N ccoll ≤ 2

√

π d2nN .

(29)

On the other hand the equilibrium mean free path reads

√ λ = ( 2π d2n)−1

and thus

 N ccoll

≤ 

⇒ K nc ≤

(30)

2 Nc

(31)

π K nc 2 Nc

π N ccoll

≤

N N ccoll

.

(32)

From the above equation one can obtain a switching criterion between DSMC and FP based on N c / N ccoll . If the number of collisions in a cell becomes larger than the number of particles then FP will be performed and otherwise DSMC. Below we analyze this criterion for the hybrid approach in terms of efficiency and accuracy.

116

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

Table 1 Outline of the algorithm. 1234567-

Advect the particles according to Eq. (20) (similar to DSMC). Apply boundary conditions (similar to DSMC). Assign FP or DSMC to the computational cell based on condition (34). Calculate the macroscopic coefficients in each FP cell from Eq. (A.4). Update the velocities based on FP (Eq. (17)) in each FP cell. Update the velocities based on DSMC in each DSMC cell. Sample the results (similar to DSMC).

1. Accuracy: The accuracy of the scheme should be analyzed in both the DSMC and the FP regions. Note that in the DSMC region, where N ccoll ≤ N c , the local Knudsen number is larger than one. Therefore the macroscopic length scale is smaller than the mean free path and thus cell and time step sizes are adequate to resolve the collisional scales, which ensures accurate DSMC results. In the FP regions, on the other hand, K nc < 1, which ensures accurate FP results. 2. Efficiency: For DSMC, the computational cost depends on two factors: the number of collisions and the spatio-temporal discretization. However, since the resolution is chosen such that the macroscopic length and time scales are resolved, it is not required to consider the latter factor. Moreover, to ensure accuracy, DSMC is employed in regions where N ccoll ≤ N c and thus the cost scales with the number of particles. For FP regions, one should note that the cost also scales with the number of particles, but independent of K nc . Therefore provided there always exist at least few particles per computational cell, the overall cost of the solution algorithm approximately scales with the number of particles independent of the Knudsen number. Here an efficient discretization practice would be to employ a fine spatio-temporal resolution in the high Knudsen regime and coarser ones elsewhere. This can be achieved by utilizing local time stepping and multi-level grids. Note that the former involves correction of statistical particle weights [3]. As a result of this practice, the FP scheme will be chosen in regions where coarser computational cells and larger time steps are utilized, whereas DSMC will be employed at high Knudsen regimes where enough resolution for molecular scales is provided. For a general molecular model, one has to approximate the number of collisions for a given time step in a given computational cell. Here, the selected number of collision partners is estimated as

N csel =

1 N c N c  F N (σ T cr )max t 2

Vc

(33)

and is used as a parameter for the hybrid method. The number of molecules represented by each computational particle is F N = w /m, V c is the volume of the computational cell, N c  is the average number of particles in the cell, σ T is the total collision cross section and cr is the relative speed between the collision partners. Therefore, the choice between DSMC and FP is made as follows:



chosen algorithm ∈

DSMC if N csel ≤ N c , FP

if N c < N csel .

(34)

The algorithm for each time step is given in Table 1. 4. Results In this section, the performance of the hybrid method is assessed and demonstrated. In particular, the accuracy and computational efficiency are investigated. Therefore we present five test cases. First, the relaxation from a non-equilibrium initial condition is studied. Second, the method is evaluated for the hydrodynamic limit and a planar subsonic Couette flow is considered. For the third test case, the shock tube problem with Sod’s conditions is tackled by the proposed hybrid FP–DSMC method. Thereafter, flow through a slit is considered with an infinitely large pressure ratio and Kn = 0.01. In this case Kn varies substantially across the slit resulting in both DSMC and FP regions. The last test case deals with plume impingement in vacuum, where the Knudsen number varies from almost zero to infinity. In all presented test cases argon is considered as working fluid. The following data is based on [5]. The viscosity is μ0 = 2.17 × 10−5 kg m−1 s−1 at T 0 = 273 K, the viscosity exponent is ω = 0.81 and the molecular mass is m = 66.3 × 10−27 kg. The Prandtl number for the FP calculations was set to Pr = 2/3. The DSMC code used here is based on Bird’s no time counter (NTC) method for the variable hard sphere (VHS) model and was developed by the author [5]. All the computational times reported in this paper are from simulations performed on a single CPU-Intel ® @ 3.40 GHz. 4.1. Homogeneous relaxation Relaxation of the velocity MDF can be studied using the hybrid FP–DSMC method or DSMC. Here a meaningful comparison between the two schemes is possible, if time step sizes larger than the mean collision time are considered. Since

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

117

Fig. 1. Relaxation of the velocity PDF. (a) t = 0 s; (b) t = t; (c) t = 2t; (d) t = 3t. The DSMC results are shown by the red solid lines and the hybrid FP–DSMC results by the black dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

for smaller time step sizes the hybrid FP–DSMC method switches to DSMC (due to the fact that the number of collisions becomes smaller than the number of particles). Here we consider relaxation of 2 × 107 number of particles initially sampled from a non-equilibrium MDF



F (V , t0 ) =

ρV1 V ·V exp − 2π θ 3 2θ 2



,

(35)



where V 1 ≥ 0, θ = kT 0 /m and T 0 = 273 K. Further, a time step size of t = μ/ P was chosen, which results in around 5 × 107 number of collisions per time step. Note that the hybrid FP–DSMC scheme chooses the FP model for each time step (according to criterion (34)). The relaxations of MDFs are shown in Fig. 1, where overall good agreement between both methods is visible. Here one should note that the difference between the Fokker–Planck and DSMC vanishes for larger time steps, i.e. for higher N coll / N p per time interval. Another important observation here is that the FP model can capture non-equilibrium behavior of the velocity density in contrast to equilibrium particle methods (such as low diffusion particle models [17,18,20]). 4.2. Couette flow The problem considers flow inside a planar channel, driven by parallel motion of the walls. Assuming infinitely long walls, the flow field becomes one dimensional, i.e. only a function of the coordinate normal to the walls. The purpose of this problem here is twofold. First, we want to analyze the accuracy of the hybrid method in the hydrodynamic limit and second, the computational efficiency is estimated by calculating the viscosity and heat conductivity of the gas. Physical condition: Since we are interested in the hydrodynamic regime, Kn = λ/ L = 0.01 is considered, where L is the width of the channel and λ is calculated based on the VHS model [5]. The coordinate normal to the walls is denoted by x1 , while

118

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

Fig. 2. Subsonic Couette flow. Mean velocity U 2 [m s−1 ], temperature T [K], heat flux in normal direction q1 [W m−2 ] and the stress component π12 [N m−2 ] are shown in the upper left, upper right, bottom left and bottom right plots, respectively. Note that the DSMC results are shown by dashed red lines, whereas the hybrid FP–DSMC results are shown by solid black lines with symbols. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

x2 is the coordinate parallel to the walls. The walls are isothermal at temperature T w = 300 K and fully diffusive Maxwell boundary condition was employed. The velocities are ±150 m s−1 for the left and the right walls, respectively, and the initial density of the gas is ρ = 0.001 kg m−3 . Numerical parameters: Different numbers nc of equidistant computational cells were used for simulations, i.e. nc ∈ {20, 40, 60, 80}. To keep the statistical error small, in average 5000 particles were used per each computational cell. The time step size was calculated based on a CFL criterion, i.e.

t = 0.5



L

nc kT w /m

.

(36)

Around 10,000 time steps had to be performed until steady state was achieved in the case with nc = 80. From then on, time averaging was employed until the results became sufficiently smooth. Mean velocity, temperature, heat flux and molecular stress are shown in Fig. 2 for both DSMC and the hybrid method with nc = 80. Note that since here all cells are larger than the mean free path, the hybrid method detects no DSMC region and thus only the FP model is employed. The agreement is very good for all quantities. The computational time for DSMC was two times larger than for the hybrid simulations. A more detailed analysis is possible by viscosity and heat conductivity calculations based on simulation results. We employed

∂U2 and ∂ x1 ∂T q1 = −κ ∂ x1

π12 = μ

(37) (38)

to compute viscosity μ and heat conductivity κ , respectively. Note that attention has to be paid to the extraction of this data, e.g. the Knudsen layer should be discarded and also the coefficients should appropriately be scaled with local temperature (see [5] for details). The convergence study was conducted based on the errors of μ and κ with respect to the reference

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

119

Fig. 3. Efficiency analysis for Couette flow. The error curves for DSMC and the hybrid FP–DSMC method with respect to the normalized CPU time are marked by triangles and circles, respectively. The different CPU times correspond to simulations with nc ∈ {20, 40, 60, 80}. On the left, the errors in viscosity and on the right the errors in heat conductivity are depicted. Note that the CPU times are normalized by the CPU time of the hybrid FP–DSMC simulation corresponding to nc = 20.

values μ0 = 2.17 × 10−5 kg m−1 s−1 and Pr = 2/3. The errors of the coefficient with respect to computational cost are shown in Fig. 3. In general, the errors of both methods are in the same order for a given grid size, whereas the hybrid approach is more efficient (e.g. 6 times faster for nc = 20). 4.3. Sod tube For the third test case, we consider a one dimensional shock tube problem with the conditions given by Sod [31]. The gas experiences three phenomena through the shock tube, i.e. rarefaction, shock and contact discontinuity, which result in a challenging test case for assessment of numerical schemes and kinetic models. Here our objective is twofold. First, the accuracy of the hybrid FP–DSMC model will be evaluated in comparison to the pure DSMC and to the compressible Navier– Stokes. Second, the convergence behavior of the integration scheme employed by the hybrid FP–DSMC will be examined in comparison to DSMC. Note that the set of parameters for the Sod tube are chosen such that the Knudsen numbers of the left and right chambers become 0.01 and 0.08, respectively, similar to the test case 2 of [32]. The reason for this specific choice is the fact that at higher Knudsen numbers, the proposed hybrid scheme switches to DSMC almost in the whole domain. Physical parameters: The gas is initially at rest in a one dimensional tube of length L = 1 m with a discontinuity in density at the middle, i.e. at x = 0.5 m. The initial densities are ρl = 10−5 kg m−3 and ρr = 0.125 × 10−5 kg m−3 , where the left and right conditions are denoted by subscripts l and r, respectively. Similar to Sod’s test case, initial temperatures are T l = 273.008012 K and T r = 273.00641 K. The boundary conditions at the left and right ends of the tube are open boundaries with equilibrium conditions at temperatures T l and T r , respectively. Numerical parameters: The computational domain is discretized into 200 equidistant cells for both DSMC and the hybrid FP–DSMC simulations. Since the flow is unsteady and time averaging is not allowed, a large number of particles should be employed. Therefore, initially 2.5 × 106 and 20 × 106 computational particles are distributed in the right and left chambers, respectively. Similar to [32], the reference time step size is chosen t ref = 7.41 × 10−6 s. The simulation runs until t end = 6.8 × 10−4 s. First, the accuracy of the hybrid FP–DSMC model is examined. Fig. 4 shows the results for the mean velocity and temperature of the gas, using three different models, i.e. the hybrid FP–DSMC, DSMC and Navier–Stokes. Note that the Navier–Stokes results were taken from [32]. Excellent agreement between the hybrid FP–DSMC and DSMC can be observed, whereas Navier–Stokes provides good results for the mean velocity, but the resulting temperature is not accurate. Fig. 4d shows the chosen scheme for each computational cell in the hybrid simulation. Not surprisingly, in the high Knudsen regime at the right-end DSMC is employed, whereas the FP model was used elsewhere. For this set of simulations, the CPU time of the hybrid FP–DSMC simulation was 2 times lower than for pure DSMC. In order to better assess the time integration scheme used in the presented hybrid method, here we conduct the same simulation however with different time step sizes, i.e. t ∈ {2t ref , 6t ref , 10t ref }. For comparison, we consider the following error measure

 E T = max

ref

|T i − T i | ref

Ti

 ,

(39)

120

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

Fig. 4. Sod tube. (a) Streamwise velocity component U [m s−1 ]; (b) temperature T [K]; (c) normalized density ρ /ρl ; (d) corresponding FP/DSMC assignment in the hybrid simulation. The DSMC results are depicted by red dashed lines, the hybrid FP–DSMC results by solid black lines with symbols and the Navier–Stokes results of [32] by black dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) ref

where T i is the computed value of temperature at cell number i, and T i the corresponding value computed by DSMC with the time step size of t ref . Fig. 5 shows the convergence efficiency of both hybrid FP–DSMC and DSMC simulations in terms of E T . It is interesting to see that for the simulation with the largest time step size, i.e. t = 10t ref , the hybrid FP–DSMC scheme results in 4 times faster simulations and yet better accuracy compared to the corresponding DSMC simulation. 4.4. Flow through a slit Gas flow through a slit has many practical applications ranging from microfluidics to space systems [35,36]. The flow is generated due to the pressure gradient across the slit. If the length of the slit is small enough, non-continuum effects in the gas flow have to be taken into account. Physical parameters: The setup considered here consists of two reservoirs at pressures P 1 > 0 and P 2 = 0 connected through a slit of width δ as shown in Fig. 6. Here we make use of the symmetry and only simulate the upper half of the flow. The height of the reservoirs is L = 10δ , where δ = λ/Kn and Kn = 0.005. The mean free path λ is calculated based on the upstream condition, i.e. the number density n1 = 1020 m−3 and temperature T 1 = 273 K as



λ=

μ 2kT 1 /m P1

.

(40)

The resulting upstream pressure based on n1 and T 1 is approximately P 1 ≈ 0.38 Pa. For the right chamber vacuum is assumed. The two chambers are separated by an isothermal wall of 273 K, for which fully diffusive Maxwell boundary model was employed. The upstream reservoir is fed by the equilibrium molecular flux continuously. Numerical parameters: An equidistant Cartesian mesh was used with nc = 180 cells in the main flow direction X and 90 cells in the normal direction Y . Initially each cell was occupied by 100 particles. Due to the sharp acceleration and expansion in

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

121

Fig. 5. Efficiency analysis for Sod tube. The error curves of temperature for DSMC and the hybrid FP–DSMC method with respect to the normalized CPU time are marked by triangles and circles, respectively. The different CPU times correspond to simulations with t ∈ {2, 6, 10}t ref . Note that the CPU times are normalized by the CPU time of the hybrid FP–DSMC simulation corresponding to t = 10t ref .

Fig. 6. Test case of gas flow through a slit.

the vacuum chamber, the number of particles significantly drops in theright part of the domain. Therefore different time step sizes were used for the left and right chambers, i.e. t 1 = 0.5L /nc / 2kT 1 /m and t 2 = 0.1t 1 , respectively. Note that the particles weights were corrected accordingly [3,33], resulting in approximately 485,000 particles overall in steady state. Note that this local time stepping approach results in using the FP scheme in the left chamber, where the low Knudsen number is very low, and thus increases the efficiency of the hybrid FP–DSMC simulation. In order to reach steady state, 20,000 time steps were required before about 100,000 time steps were performed to get smooth averaged results. Color coded Mach contours are shown in Fig. 7 representing both DSMC and hybrid FP–DSMC results. As expected, the flow is subsonic on the upstream side of the slit, sonic in the slit and supersonic on the right side reaching Ma ≈ 5. The agreement between the results of both methods is very good. Fig. 7b depicts the computational cells composing DSMC and FP regions at a given instance in the hybrid simulation in red and blue, respectively. As expected it can be observed that on the upstream side, only the FP model was used, whereas the scheme switched to DSMC on the downstream side. In addition, Mach-, temperature- and velocity contour lines are shown in Fig. 8, where close again agreement between results of the two methods can be observed. In order to compare the results with literature data, the reduced mass flow rate, i.e.

W =

˙ M ˙0 M

(41)

, √ 

˙ is the mass flow rate and M ˙ 0 = δ P 1 /( π 2kT 1 /m) the mass flow rate in the free molecular was calculated, where M regime for the planar geometry. The corresponding reduced mass flow rates are shown in Table 2, where close agreements between our results and the results reported in [35] can be observed. The CPU time required for the hybrid method to reach was around 10 hours, while DSMC took around 70 hours for the same number of time steps on a single CPU. It is worthwhile to mention that the efficiency gain of the hybrid method compared to DSMC could become more pronounced, if e.g. the size of the problem is increased or the upstream Knudsen number is reduced.

122

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

Fig. 7. (a) Mach contours of gas flow through slit. The DSMC results are shown in the upper half plane and the hybrid FP–DSMC results are shown in the lower half plane. (b) FP/DSMC assignments for hybrid simulation of gas flow through slit. DSMC was assigned to computational cells depicted in red and the FP scheme to computational cells depicted in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Gas flow through slit. (a) Contours of streamwise velocity component U 1 [m s−1 ]; (b) contours of transverse velocity component U 2 [m s−1 ]; (c) contours of temperature; (d) Mach contours. The DSMC results are depicted by red solid lines and the hybrid FP–DSMC results are depicted by blue dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

123

Table 2 Reduced mass flow rate at Kn = 0.005 and P 2 = 0. Method

W

Hybrid FP–DSMC DSMC DSMC by Sharipov [35]

1.54 1.54 1.55

Fig. 9. Gas flow through slit. (a) Temperature contours; (b) Mach contours. The FP results are depicted by black solid lines and the hybrid FP–DSMC results are depicted by blue dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Sketch of plume impingement test case.

For completeness, we have also performed pure FP simulations while expecting inaccuracies on the vacuum side. Comparisons with respect to the hybrid FP–DSMC results are shown in Fig. 9 for temperature and Mach contours. Note that even in the high Knudsen regions there only exist moderate deviations. 4.4.1. Plume impingement As another representative example, where the gas flow experiences different Knudsen regimes, a plume impingement test case is studied here. The exhaust gas accelerates through a convergent–divergent nozzle and afterwards expand into the vacuum. Finally, the plume impinges on a solid wall downstream of the nozzle. Such flow scenarios can occur e.g. due to satellite attitude control [34]. Due to the extreme variation of Kn inside the flow domain, this problem is ideal to demonstrate the advantage of hybrid schemes; note that pure DSMC may not be feasible due to continuum condition at the inflow boundary. Problem setup: The geometrical setup is a two-dimensional convergent–divergent nozzle with parameters shown in Fig. 10. Argon enters the nozzle at a speed of U in = 5.3 m s−1 , which corresponds to Ma = 0.01. The number density is nin = 6.63 × 1023 m−3 . For the entrance section of the nozzle specular walls are considered while for the convergent–divergent

124

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

Fig. 11. Plume impingement. (a) Mach contours overlaid by velocity streamlines. The DSMC (with FP inflow) results are shown on the upper half plane and the hybrid FP–DSMC results are shown on the lower half. (b) FP and DSMC assignment color coded in blue and red, respectively. The assignment map for DSMC with FP inflow is shown on the upper half plane and the assignment map for hybrid FP–DSMC is shown on the lower half plane. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12. Plume impingement. (a) Contour lines of streamwise velocity component U 1 [m s−1 ]. (b) Velocity component U 2 [m s−1 ]. (c) Temperature contours T [K]. (d) Mach contours. The DSMC (with FP inflow) results are depicted by red solid lines and the hybrids FP–DSMC results by blue dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

125

Fig. 13. Plume impingement. (a)–(d) Grids D–A on the lower half-planes and corresponding FP/DSMC assignment maps on the upper half-planes. FP and DSMC regions are coded in blue and red, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

section diffusive wall boundary conditions with at T w = 587 K are applied. For the plate opposite the nozzle exit specular boundary conditions are assumed. At the outflow boundaries vacuum condition is employed. Numerical parameters: Simulations were performed for different numbers of computational cells, i.e. for nc = 3600 (grid A), nc = 8100 (grid B), nc = 14,400 (grid C) and nc = 22,500 (grid D) as shown in Fig. 13. Initially, each computational cell was populated by around 500 particles. The time step size was chosen based on an average CFL criterion resulting in t ∈ {9.68, 6.45, 4.84, 3.87} × 10−8 s for the grids A, B, C and D, respectively. For simplicity constant time step sizes were employed for each set of simulations. Approximately 10,000 time steps were performed to reach steady state and another 10,000 to perform time averaging in order to reduce the statistical noise. Like in the previous test case, symmetry was exploited to reduce the computational cost. With each grid two simulations were performed. One with the described hybrid method, while for the other DSMC and FP simulations were performed right and left of the nozzle bottle neck, respectively. As observed in Fig. 11b, DSMC with FP inflow simulations basically uses FP only in the hydrodynamic region inside the nozzle. Therefore one can expect that the accuracy of the method be as high as full DSMC. Note that full DSMC was not possible due to the high number of collisions in the converging part of the nozzle. For accuracy analysis grid D shown at Fig. 13a was chosen. The flow field is depicted in Fig. 11a, where hybrid and DSMC (with FP inflow) results are compared. The exhaust gas goes through a semi-circular bow shock, which interacts with the expansion waves at the end of the nozzle. Note that in the absence of the plate in front of the nozzle exit, the flow would accelerate to higher Mach numbers and no shock would exist at the nozzle exit. More detailed comparisons are presented in Fig. 12 for Mach-, temperature- and velocity contours. The agreement is excellent. To assess the efficiency of the hybrid scheme, simulations on grid A, grid B and grid C were performed. Fig. 13 shows the grids together with colormaps indicating FP and DSMC regions (in blue and red, respectively). It can be observed that the fraction of cells assigned to DSMC is smaller for coarser grids, which is based on condition (34). Therefore, the computational

126

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

Fig. 14. Plume impingement. (a), (b) and (c): Mach contours obtained with grids C, B and A, respectively. Red solid contours represent the Mach number resulting from DSMC with FP in the inflow region, blue dash-dotted contours result from the hybrid method and black solid lines (lower half-planes) result from DSMC with FP in the inflow region computed with the finest grid D (reference). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 3 Computational time comparison for plume impingement simulations. The CPU time per time step per particle is reported for DSMC and for the hybrid FP–DSMC method. Grid type

CPU time per particle per time step [μs] Hybrid FP–DSMC

DSMC

Grid Grid Grid Grid

0.60 1.10 1.20 1.25

43.7 30.1 22.8 18.60

A B C D

advantage of the hybrid method compared to DSMC becomes more pronounced as the grid gets coarser. To asses the quality of different hybrid results, Mach contours are shown in Fig. 14 for all grids in comparison with finest DSMC results, with FP in the inflow region. For this test case the hybrid method shows a slightly better convergence behavior than DSMC with FP in the inflow region. Table 3 shows the computational cost per time step and particle for both the hybrid method and pure DSMC. It can be observed that the hybrid method is much more efficient for all four grids and the comparison confirms that the difference becomes more dramatic for coarser grids. For example considering grid A, the hybrid is almost 70 times faster than pure DSMC. It is worthwhile to mention that the hybrid simulation with grid A took around 5 hours on the mentioned CPU. 5. Conclusion and outlook The Fokker–Planck approximation of the Boltzmann collision operator has been used in order to extend the applicability of DSMC for cases including near continuum and continuum regions. We have shown that the number of collisions between

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

127

gas molecules can serve as an indicator for switching between DSMC and the FP method. That is, if the selected number of collision pairs becomes larger than the number of particles in the computational cell, the FP method is employed and otherwise DSMC. Different test cases are considered, including planar Couette flow in the hydrodynamic limit, Sod tube, flow through a slit and two dimensional plume–wall interaction. Note that in the latter two, whole Knudsen number range has to be modeled. Comparisons with DSMC have shown that the efficiency of the presented hybrid method is excellent and that the computational cost is much lower. At the end it is necessary to mention certain developments that have to be carried out in the next development steps of the hybrid FP–DSMC solution algorithm: 1. Validation studies: More validation studies have to be performed. In particular, the performance of the hybrid FP–DSMC scheme should be assessed for three dimensional complex geometries. Note that the FP scheme, similar to DSMC, is intrinsically three dimensional and thus no special generalization is needed for the corresponding simulations. 2. Complex gas flows: Further, the diatomic FP model derived in [24] should be extended for reactive gas flows, where the effects of dissociation and ionization can be described. The accuracy and efficiency of the hybrid scheme then should be assessed for these scenarios. Acknowledgement Funding for this research was provided by the Swiss National Science Foundation under the grant number 153116. Appendix A. Macroscopic coefficients Here we review the system of equations, which are required to calculate the model coefficients c i j and γi . These coefficients appear in the drift term and are responsible for correct evolution of the second and third order moments by the FP operator with respect to the Boltzmann equation. The model coefficients have to fulfill the linear equations (2 )

(2 )

c il ulj + c jl uli + γi u j + γ j u i

(2 )

+ 2 u i j = 0

(A.1)

for stresses and

      5 (2 ) 2

3u (i 4) − ul(2) uli − u (i 2) u (2) + 2c jl uli j + c il ul(2) + 2γ j u (i j2) − u (2) u i j + γi u (4) − (u (2) ) = ui 3τ

(A.2)

for heat fluxes, where ( p)

ui

1 i 2 ...i n



1

=

F | v | p v i 1 . . . v in dV 1 dV 2 dV 3 .

ρ

(A.3)

R3

For convenience, the constitutive relations, i.e. Eqs. (A.1)–(A.2) can be written as the linear system



L M



c





=

γ

Z R



(A.4)

consisting of 9 equations (6 from (A.1) and 3 from (A.2)) and 9 unknowns (6 c i j and 3 one obtains

⎛

⎜ ⎜ ⎜ ⎜ L =⎜ ⎜ ⎜ ⎜ ⎝

( 0)

( 0)

2u 11

( 0)

2u 13

( 0)

u 11 + u 22

( 0)

u 23

u 12

( 0)

( 0)

u 13

⎛

( 0)

2u 12

( 0)

( 0)

u 23 ( 0)

2u 12

0

u 13

u 12

0

2u 13

0 (2 )

0

u2

( 0)

u 13

( 0)

u3

0

0

0

2u 2

( 0)

0

u3

0

0

u 12

( 0)

( 0)

2u 23 ( 0)

( 0)

u 23

u 22 + u 33

0

2u 23

( 0)

( 0)

( 0)

0 2u 22

( 0)

2u 1

u 13

( 0)

0

0

( 0)

( 0)

( 0)

(2 )

0

u 12

u 11 + u 33

0

( 0)

0

(2 ) (2 )

u 23

( 0)

( 0)

( 0)

u 33

⎞

u 1 + 2u 111 2u 112 2u 113 (2 ) ( 0) ( 0) ⎟ ⎜ u (2) + 4u (0) u 1 + 4u 122 4u 123 ⎟ ⎜ 2 112 ⎟ ⎜ ( 2 ) ( 0 ) ( 0 ) ( 2 ) ( 0 ) ⎜ u + 4u 4u 123 u 1 + 4u 133 ⎟ 3 113 ⎟ ⎜ ⎟ ⎜ ( 0) (2 ) ( 0) ( 0) 2u 122 u 2 + 2u 222 2u 223 ⎟ ⎜ ⎟ ⎜ ( 0) (2 ) ( 0) (2 ) ( 0) T ⎜ M =⎜ 4u 123 u 3 + u 223 u 2 + 4u 233 ⎟ ⎟ ⎟ ⎜ ( 0) ( 0) (2 ) ( 0) ⎜ 2u 133 2u 233 u 3 + 2u 333 ⎟ ⎟ ⎜ ⎜ 2u (2) − 2u (0) u (2) 2u (2) − 2u (0) u (2) 2u (2) − 2u (0) u (2) ⎟ ⎟ ⎜ 11 11 12 12 13 13 ⎟ ⎜ (2 ) ⎝ 2u − 2u (0) u (2) 2u (2) − 2u (0) u (2) 2u (2) − 2u (0) u (2) ⎠ 12 (2 )

12 ( 0)

2u 13 − 2u 13 u (2)

22 (2 )

22 ( 0)

2u 23 − 2u 23 u (2)

23 (2 )

23 ( 0)

2u 33 − 2u 33 u (2)

γi ). Based on the constitutive relations

0 (2 )

u1

(2 )

(2 )

0

⎞

⎟ ⎟ (2 ) ⎟ u1 ⎟ ⎟, ⎟ 0 ⎟ (2 ) ⎟ u2 ⎠ 0

(2 )

2u 3

(A.5)

128

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

⎛

0 0 0 0 0 0

0 0 0 0 0 0 0

0 0 0 0 0 0 0

⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ +⎜ ⎜ ⎜ ⎜ ( 4 )  ( 2 ) 2 ⎜u − u ⎜ 2  ⎜ 0 u (4 ) − u (2 ) ⎝ 0

⎛

0

(2 ) ⎞

0

u (4 )

2  − u (2 )

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(A.6)

u 11

⎜ (2 ) ⎟ ⎜ u 12 ⎟ ⎜ (2 ) ⎟ ⎜u ⎟ ⎟ ⎜ Z = −2 ⎜ 13 ⎟ ⎜ u (2 ) ⎟ ⎜ 22 ⎟ ⎜ (2 ) ⎟ ⎝ u 23 ⎠

(A.7)

(2 )

u 33 and

 ⎞ (4 ) (2 ) ( 0) ( 2 ) ( 0) ( 2 ) ( 0) ( 2 ) −3u 1 + u (2) u 1 + 2 u 11 u 1 + u 12 u 2 + u 13 u 3 u1 ⎜  ⎟ 5 ⎜ (2 ) ⎟ ⎜ (4 ) (2 ) ( 0) ( 2 ) ( 0) ( 2 ) ( 0) ( 2 ) ⎟ R= ⎝ u 2 ⎠ + ⎜ −3u 2 + u (2) u 2 + 2 u 12 u 1 + u 22 u 2 + u 23 u 3 ⎟ . ⎝ 3τ  ⎠ (2 ) (4 ) (2 ) ( 0) ( 2 ) ( 0) ( 2 ) ( 0) ( 2 ) u3 −3u 3 + u (2) u 3 + 2 u 13 u 1 + u 23 u 2 + u 33 u 3 ⎛

(2 ) ⎞

⎛

(A.8)

References [1] P. Jenny, M. Torrilhon, S. Heinz, A solution algorithm for the fluid dynamics equations based on a stochastic model for molecular motion, J. Comput. Phys. 229 (2010) 1077–1098. [2] M.H. Gorji, M. Torrilhon, P. Jenny, Fokker–Planck model for computational studies of monatomic rarefied gas flows, J. Fluid Mech. 680 (2011) 574–601. [3] M.H. Gorji, P. Jenny, An efficient particle Fokker–Planck algorithm for rarefied gas flows, J. Comput. Phys. 263 (2014) 325–343. [4] G.A. Bird, Approach to translational equilibrium in a rigid sphere gas, Phys. Fluids 6 (1963) 1518–1519. [5] G.A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Clarendon, Oxford, 1994. [6] E.M. Shakhov, Generalization of the Krook kinetic relaxation equation, Fluid Dyn. 3 (1968) 142–145. [7] L.H. Holway, New statistical models for kinetic theory: methods of construction, Phys. Fluids 9 (1965) 1658. [8] J. Fan, C. Shen, Statistical simulation of low-speed rarefied gas flows, J. Comput. Phys. 167 (2001) 393–412. [9] T.M.M. Homollea, N.G. Hadjiconstantinou, A low-variance deviational simulation Monte Carlo for the Boltzmann equation, J. Comput. Phys. 226 (2007) 2341–2358. [10] S.K. Stefanov, On DSMC calculation of rarefied gas flows with small number of particles in cells, SIAM J. Sci. Comput. 33 (2011) 677–702. [11] L. Pareschi, G. Russo, Time relaxed Monte Carlo methods for the Boltzmann equation, SIAM J. Sci. Comput. 23 (2001) 1253–1273. [12] P. Degond, L. Dimarco, L. Pareschi, The moment-guided Monte Carlo method, Int. J. Numer. Methods Fluids 67 (2) (2011) 189–213. [13] D.B. Hash, H.A. Hassan, Assessment of schemes for coupling Monte Carlo and Navier–Stokes solution, J. Thermophys. Heat Transf. 10 (2) (1996) 242–249. [14] W.L. Wang, I.D. Boyd, Hybrid DSMC–CFD simulations of hypersonic flow over sharp and blunted bodies, AIAA paper 3644, 2003. [15] O. Aktas, N.R. Aluru, A combined continuum/DSMC technique for multiscale analysis of microfluidic filters, J. Comput. Phys. 178 (2) (2002) 342–372. [16] M.A. Gallis, J.R. Torczynski, Investigation of the ellipsoidal-statistical Bhatnagar–Gross–Krook kinetic model applied to gas-phase transport of heat and tangential momentum between parallel walls, Phys. Fluids 23 (3) (2011) 030601. [17] D.I. Pullin, Direct simulation methods for compressible inviscid ideal-gas flow, J. Comput. Phys. 34 (2) (1980) 231–244. [18] J.M. Burt, I.D. Boyd, A low diffusion particle method for simulating compressible inviscid flows, J. Comput. Phys. 227 (2008) 4653–4670. [19] J.M. Burt, I.D. Boyd, Extension of a multiscale particle scheme to near-equilibrium viscous flows, AIAA J. 47 (2009) 1507–1517. [20] B.J. Albright, D.S. Lemons, M.E. Jones, D. Winske, Quiet direct simulation of Eulerian fluids, Phys. Rev. E 65 (2002) 055302. [21] S.V. Bogomolov, On Fokker–Planck model for the Boltzmann collision integral at the moderate Knudsen numbers, Math. Model. 21 (2009) 111–117. [22] S.V. Bogomolov, I.G. Gudich, Diffusion model of gas in a phase space for moderate Knudsen numbers, Math. Models Comput. Simul. 5 (2013) 130–144. [23] A.V. Skorokhod, Stochastic Equations for Complex Systems, Springer, 1988. [24] M.H. Gorji, P. Jenny, A Fokker–Planck based kinetic model for diatomic rarefied gas flows, Phys. Fluids 25 (2013) 062002. [25] M.H. Gorji, P. Jenny, A kinetic model for gas mixtures based on a Fokker–Planck equation, J. Phys. Conf. Ser. 362 (2012) 012042. [26] I. Shevtsova, On the absolute constants in the Berry–Esseen type inequalities for identically distributed summands, arXiv:1111.6554, 2011. [27] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, Cambridge University Press, 1970. [28] F. Fei, J. Fan, A diffusive information preservation method for small Knudsen number flows, J. Comput. Phys. 243 (2013) 179–193. [29] R.F. Pawula, Approximation of the linear Boltzmann equation by the Fokker–Planck equation, Phys. Rev. 162 (1967) 186–188. [30] C.W. Gardiner, Handbook of Stochastic Methods, Springer-Verlag, 1985. [31] G.A. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, J. Comput. Phys. 27 (1978) 1–31. [32] S. Tiwari, A. Klar, S. Hardt, A particle–particle hybrid method for kinetic and continuum equations, J. Comput. Phys. 228 (2009) 7109–7124. [33] J.S. Wu, K.C. Tseng, F.Y. Wu, Parallel three-dimensional DSMC method using mesh refinement and variable time-step scheme, Comput. Phys. Commun. 162 (2004) 166–187.

M.H. Gorji, P. Jenny / Journal of Computational Physics 287 (2015) 110–129

129

[34] K.C. Kannenberg, I.D. Boyd, Three-dimensional Monte Carlo simulations of plume impingement, J. Thermophys. Heat Transf. 13 (1999) 226–235. [35] F. Sharipov, D.V. Kozak, Rarefied gas flow through a thin slit into vacuum simulated by the Monte Carlo method over the whole range of the Knudsen number, J. Vac. Sci. Technol., A, Vac. Surf. Films 27 (3) (2009) 479–484. [36] F. Sharipov, Numerical simulation of rarefied gas flow through a thin orifice, J. Fluid Mech. 518 (2004) 35–60.