Adaptive particle–cell algorithm for Fokker–Planck based rarefied gas flow simulations

Accepted Manuscript Adaptive particle cell algorithm for Fokker–Planck based rarefied gas flow simulations M. Pfeiffer, M.H. Gorji

PII: DOI: Reference:

S0010-4655(16)30357-5 http://dx.doi.org/10.1016/j.cpc.2016.11.003 COMPHY 6091

To appear in:

Computer Physics Communications

Received date : 14 September 2016 Revised date : 13 November 2016 Accepted date : 18 November 2016 Please cite this article as: M. Pfeiffer, M.H. Gorji, Adaptive particle cell algorithm for Fokker–Planck based rarefied gas flow simulations, Computer Physics Communications (2016), http://dx.doi.org/10.1016/j.cpc.2016.11.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Adaptive Particle Cell Algorithm for Fokker-Planck Based Rarefied Gas Flow Simulations M. Pfeiffera , M.H. Gorjib a

Institute of Space Systems, University of Stuttgart, Pfaffenwaldring 29, D-70569 Stuttgart, Germany b Department of Mathematics, RWTH Aachen University, Schinkelstrasse 2, D-52062 Aachen, Germany

Abstract Recently, the Fokker-Planck (FP) kinetic model has been devised on the basis of the Boltzmann equation [1, 2]. Particle Monte-Carlo schemes are then introduced for simulations of rarefied gas flows based on the FP kinetics. Here the particles follow independent stochastic paths and thus a spatiotemporal resolution coarser than the collisional scales becomes possible. In contrast to the direct simulation Monte-Carlo (DSMC), the computational cost is independent of the Knudsen number resulting in efficient simulations at moderate/low Knudsen flows. In order to further exploit the efficiency of the FP method, the required particle-cell resolutions should be found, and a cell refinement strategy has to be developed accordingly. In this study, an adaptive particle-cell scheme applicable to a general unstructured mesh is derived for the FP model. Virtual sub cells are introduced for the adaptive mesh refinement. Moreover a sub cell-merging algorithm is provided to honor the minimum required number of particles per cell. For assessments, the 70 degree blunted cone reentry flow [3] is studied. Excellent agreement between the introduced adaptive FP method and DSMC is achieved. Keywords: Rarefied gas flows, Fokker-Planck, DSMC, reentry simulations

Email addresses: [email protected] (M. Pfeiffer), [email protected] (M.H. Gorji)

Preprint submitted to Journal of Computer Physics Communications November 13, 2016

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1. Introduction Rarefied gases are encountered in various phenomena ranging from astrophysical flows to reentry flights. Due to the strong departure of the flow from the thermodynamic equilibrium, very often the conventional Navier-StokesFourier equations fail to provide a satisfactory description of the rarefied gases. A higher level of closure is then accounted by the Boltzmann equation in which the dynamics of the velocity distribution function is devised. As an integro-differential equation for the velocity distribution, the direct solution of the Boltzmann equation is very costly. Therefore approximate simulation methods have to be typically employed. The dominant simulation method for rarefied gas flows is DSMC introduced by Bird [4], where the computational particles evolve in a stochastic manner consistent with the Boltzmann equation. While DSMC has been employed with significant success, its high computational cost at low Knudsen regions of the flow can become a serious deficiency. The very fact that the binary collisions have to be resolved in DSMC leads to stringent spatio-temporal restrictions at low Knudsen numbers. The Fokker-Planck solution algorithm which was originally proposed in [1] and then generalized in [2], provides an approximation of the Boltzmann equation suitable for rarefied gas flow simulations. Since the FP approximation replaces the binary collisions described by the Boltzmann equation by continuous Markov processes, the resulting computational algorithm becomes very efficient especially at low Knudsen numbers. Therefore providing accurate time integration schemes, time steps larger than the mean collision time and computational cells larger than the mean free path can be employed in the FP solution algorithm. At the same time, the FP method is capable of capturing non-equilibrium phenomena such as counter Fourier heat fluxes or Knudsen minima [2, 5] compared to conventional continuum methods e.g. Navier-Stokes-Fourier equations. It is worth mentioning that other FP models have been proposed [6, 7, 8] but their comparison with respect to the present model is beyond the scope of this work. In order to exploit intrinsic computational advantages of the FP type kinetic models properly, the necessary spatio-temporal resolution for an accurate FP simulation has to be investigated. Furthermore, the minimum number of 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

computational particles that leads to an acceptable level of statistical error has to be found. Note that the moments of the distribution function enter the FP based stochastic processes through the drift and the diffusion. Therefore the gradients of the present moments together with the corresponding CFL condition may serve as guiding criteria for the required resolution [9]. In this study, an adaptive algorithm is devised such that the necessary spatiotemporal resolution for the FP model is achieved. The main ingredients of the algorithm are borrowed from mesh refinement schemes already introduced in DSMC and Particle-In-Cell solvers ( see e.g. [10, 11, 12]). Here, virtual sub cells are defined around interpolation points where velocity moments are sampled accordingly. Furthermore, an accurate yet efficient treatment of unstructured meshes is provided through a mapping into and from the Cartesian mesh [13, 14, 15]. Finally, a sub cell merging algorithm is devised in order to fulfill the minimum required number of particles. For validation of the adaptive FP method, the results are compared with respect to the normal FP as well as DSMC. For this, simulations of a 70◦ blunted cone are performed using the mentioned methods, and the comparisons are made in terms of the accuracy and the computational efficiency. It is shown that the introduced adaptive FP scheme leads to very accurate results compared to DSMC for the mentioned test case. 2. Review of the Fokker-Planck model 2.1. Stochastic Processes The Fokker-Planck kinetic model ∂ ∂ D2 ∂ 2 F ∂F + (Vi F) = − (Ai F) + ∂t ∂xi ∂Vi 2 ∂Vi ∂Vi

(1)

describes the evolution of the velocity distribution function F(V , x, t) defined over the velocity sample space V at each point in space x and instant in time t. Here the evolution of F(V , x, t) is governed by the drift A(V , F) and the diffusion D(F). The FP model has been employed successfully as an approximation of the Boltzmann equation for low Knudsen numbers ( see e.g. [1, 2]).The cubic FP model devised in [2] provides closures for A and D such that correct viscosity µ and the Prandtl number of 2/3 can be obtained for monoatomic gases in the Navier-Stokes-Fourier limit. The macroscopic 3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

coefficients c, γ and Λ introduced in the drift Ai

 
1 0 2qi 0 0 0 0 0 0 = − vi + cij vj + γi vj vj − 3θ + Λ vi vj vj − τ ρ

(2)

together with the time scale τ = 2µ/p and θ = kT /m lead to exact relaxations of the viscous stresses σij and heat fluxes qi consistent with the Boltzmann collision operator; assuming Maxwell type molecules. Besides Λ, the coefficients c and γ constitute a 9 × 9 linear system including 6 equations (2)

(2)

cik ukj + cjk uki + γi uj + γj ui

p (2) = − pij − 2Λuij µ

(3)

for the stresses, and 3 equations  
(2) (2) cik uk + 2ckj uijk + γi u(4) − 3θu(2) + 2γj uij − 3θuij   2p qj (4) (2) qi = − qi − Λ 3ui − 4uij − 2u 3µ ρ ρ

(4)

for the heat fluxes. Here v 0 = V − U denotes the fluctuating velocity with respect to the bulk velocity U , p is the pressure given by the ideal gas law, ρ the density and µ the viscosity. Furthermore, for abbreviation we employed Z (n) |v 0 |n vi01 ...vi0k F d3 V . (5) ui1 ...ik = R3

The diffusion coefficient

r

D =

2θ τ

(6)

ensures the relaxation towards Maxwellian with the temperature T , the Boltzmann constant k and molecular mass m. Similar to DSMC, in the FP solution algorithm we do not solve Eq. (1) directly but through the equivalent system of stochastic processes. Suppose M (t) and X(t) represent random variables for particle velocity and position, respectively. Therefore the distribution function generated by the following Stochastic Differential Equations (SDEs) and

dMi = Ai dt + DdWi dXi = Mi dt 4

(7) (8)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

lead to the solution of the FP equation (1) (see [16, 5] for details). In practice, we employ a particle Monte-Carlo scheme in order to approximately solve the above SDEs. The solution algorithm derived in [5] is briefly mentioned in the following section. 2.2. Solution Algorithm Suppose the distribution function is discretized by a finite number of particles with statistical weights w whose velocities M and positions X follow the SDEs (7)-(8). Therefore the discrete system   q 1 n 0,n −∆t/τ n n+1 −∆t/τ n n n n −2∆t/τ )ηi Mi = Mi e + (1 − e )τ Ni + θ (1 − e α +Uin (9) n+1 n n Xi = Xi + Mi ∆t (10) provides velocity and position updates of each particle for the time step size ∆t = tn+1 − tn . Note that all the values with the superscript n and n + 1 are evaluated at tn and tn+1 , respectively. Here M 0 = M − U , and η is a standard normal variate. The above integration scheme employs the analytic solution related to the linear part of (7) together with a first order integration of the remainder, i.e.  
2qi 0 0 0 0 0 0 . (11) Ni = cij Mj + γi Mj Mj − 3θ + Λ Mi Mj Mj − ρ

The conservation of energy is ensured by the correction factor α α2 = 1 +



τ τ (1 − e−∆t/τ )2 hNi Ni i + 2 e−∆t/τ − e−2∆t/τ hM 0 Ni i , 3θ (12)

where h...i is the expectation. All the macroscopic coefficients present in the right hand side of Eq. (9) have to be evaluated along the particle path at tn . In order to cope with that, one has to employ a computational grid such that different macroscopic moments are sampled on the grid and then interpolated to the particle position.

5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3. Spatial, temporal and particle resolution In this section, we study the resolution requirements for the FP based particle Monte-Carlo scheme. First, the spatial discretization is investigated based on the flow gradients. Upon that, the time step prerequisite is derived from the CFL criterion. Finally, the minimum number of particles per cell for unbiased samplings of the moments is analyzed. 3.1. Spatial resolution The velocity update scheme i.e. Eq. (9) presumes that the moments of the distribution function remain constant along the particle path during ∆t. Certainly that assumption is relevant only if the moments are sampled on the fine enough grid in the physical space. That limitation leads us to define the required spatial resolution upon the gradients of the macroscopic variables similar as described in [17, 18]. Let Φ(x, t) represent a scalar flow variable, where Φ(j) and (∂Φ/∂xi )(j) are its value and gradient estimated at the center of the computational cell j, re(j) spectively. Furthermore ∆xi denotes the average extend of the cell j along the coordinate xi . Therefore we define a normalized operator (j) ∇Φ

=



∂Φ ∂xi

(j)

(j)

∆xi Φ(j)

(13)

which merely gives the normalized change of Φ across the cell j, averaged over all directions. Note that the summation convention is suppressed for indices inside brackets. For convenience, the employed flow field variables are the temperature T , the pressure p, the density ρ and the flow speed |U |. To decide whether or not a cell should be adapted, thus an upper limit of ∇max has to be introduced for the maximal change of the mentioned variables. Therefore, the cell j is too coarse and should be refined if (j)

max(∇Φ ) > ∇max ; Φ

Φ ∈ {T, p, ρ, |U |}.

(14)

Based on various sets of simulations ∇max = 0.05 proved to be a good compromise between efficiency and accuracy.

6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3.2. Temporal resolution The explicit integration applied in the position update scheme i.e. Eq. (10) leads to a CFL type criterion for the maximum allowed time step size. Providing the mean speed of molecules in each cell Z 1 c¯ = (Vi Vi )1/2 F d3 V (15) ρ R3 the CFL condition ∆t(j) = min

(j) ∆x1 , c¯(j)

(j) ∆x2 , c¯(j)

(j) ∆x3 c¯(j)

!

Cmax

(16)

can be devised. Here Cmax ≤ 1 is the Courant number and c¯(j) denotes the estimated c¯ in the cell j. While the condition (16) can be directly employed in the simulations using local time stepping techniques [11, 19], it is more convenient to adopt a simpler expression for ∆t. First, let us use an estimate for the average speed based on the thermal component and the mean velocity. Therefore the equilibrium approximation √ (17) ceq = 2 2θ + |U | (j)

is considered together with ceq being the corresponding value in the cell j. Furthermore, an average cell length is introduced ∆x =

1 X  (j) 1/3 Vcell Ncell j

(18) (j)

with the total number of cells Ncell and the cell j volume Vcell . Therefore the following CFL condition can be obtained for the time step size ∆t =

∆x (j)

max(ceq )

Cmax ,

(19)

where Cmax = 0.95 is adopted. Typically in supersonic flow simulations max(|U |) can be estimated based on the flow Mach number. Whereas in subsonic simulations |U | ≈ 0 can be applied.

7

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3.3. Particle resolution In this section, the minimum number of particles required for evaluating the moments involved in the drift and the diffusion is investigated. In general, two approaches should be considered for moment evaluations: instantaneous averaging versus moving time averaging. Let Q(x, t) be an arbitrary velocity moment Z q(V )F (V , x, t)d3 V . (20) Q(x, t) = R3

Therefore it can be approximated according to np

˜ Q(x, t) =

1 X (j) w q(M (j) ) δΩ j=1

(21)

from a finite ensemble of particles j ∈ {1, ..., np } whose positions at t are inside the small volume δΩ around point x. Note that w(j) are the corresponding statistical weights which account for the real mass that is represented by the particle j. The arithmetic averaging (21) can be seen as an approximation of the identity F (V , x, t) =

∞ X j=1

w(j) δ(M (j) − V )δ(X (j) − x)

(22)

with δ() being the Dirac delta. Equation (21) gives a formula for instantaneous averaging. For statistically stationary processes, the realizations at different times can be taken into account in order to improve the accuracy ˜ n (x) at of the moment evaluation. Therefore the moment approximation Q time step n is updated np

1 X (j) ˜ n+1 (x) = µ ˜ n (x) + (1 − µ w q(M (j) ) Q ˜Q ˜) δΩ j=1

(23)

according to the time averaging constant 0 ≤ µ ˜ < 1. Depending on the averaging strategy i.e. Eq. (21) versus (23), the minimum required number of particles per computational cell can differ. Starting with ˜ one can work value µ very close to 1 and appropriate initialization of Q, 8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

with as few as 2 or 3 particles per cell for the FP algorithm. However since Eq. (23) combined with Eq. (9) only ensures the momentum and energy conservation in a time averaged sense, it may slow down the convergence rate severely. Consider the random number corrections described in [5], which simply shifts P (j) P (j) (j) and scales η in Eq. (9) such that j ηi = 0 and j ηi ηi = 1. Therefore using instantaneous averages (21), the conservation for momentum and energy would be honored in each computational cell at each time step. Thus providing enough computational memory, Eq. (21) should be preferred for evaluating the moments involved in the drift and the diffusion. Nevertheless for sampling the final results, the moving time averaging (23) is the method of choice (assuming stationary flows). Adopting the instantaneous averaging (21), minimum number of particles required per cell becomes more crucial. Let us first consider the equilibrium limit of the cubic pmodel as a starting point. The linear closures Ai = 1/τ (Vi −Ui ) and D = 2kT /(τ m) require the mean velocity and the temperature to be sampled from particles. Therefore at least two particles per cell is required. Yet the cubic closure for the drift involves many more velocity moments, i.e. up to fifth order as can be seen from Eqs. (3) and (4). Therefore one would need at least 5 particles per cell in order to have meaningful macroscopic coefficients. In practice, the particles should move with the old velocity and without the FP velocity update if the particle number per cell np ≤ 2, the linear model should be used for 3 ≤ np < 5 and the cubic model then only if 5 ≤ np . Nevertheless typically around 10 particles per cell are required for a low bias simulation. 4. Adaptation of the particle-space resolution The spatial resolution of the moment estimators can be improved by introducing virtual sub cells. The unstructured cells are virtually split into sub cells and subsequently the sub cells are partially merged again in order to guarantee the minimum number of particles per sub cell.

9

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4.1. Defining sub cells Here a nearest interpolation point scheme is introduced in order to define the sub cells. The volume around each interpolation point then becomes identical to the sub cell which only exists virtually. Thereby, the particles are assigned to the sub cells and all particles operations i.e. moment evaluation and velocity updates are performed at the sub cell level. Interpolation points. There exist different types of interpolation points e.g. equidistant and Gauss points, among others. The simple equidistant points are the best choice for our purpose, due to the fact that the sub cells become equally weighted. In order to deal with an unstructured mesh, the calculations should be carried out in a reference element space E = [−1, 1]3 and then mapped back to the real space [14, 13] . Therefore, given Ni,inter interpolation points per direction i, the positions of the interpolation points in the reference space read (j)

ξi

=

2j + 1 − 1 ; i ∈ {1, 2, 3} Ni,inter

(24)

with j ∈ {0, .., Ni,inter − 1}. Using this definition, the directional weight of each interpolation point wi = 2/Ni,inter leads to the volume ref Vinter = w1 w2 w3

(25)

in the reference space. Let h(ξ) be the mapping 1 h(ξ) = { P1 (1 − ξ1 )(1 − ξ2 )(1 − ξ3 ) + P2 (1 + ξ1 )(1 − ξ2 )(1 − ξ3 ) + 8 P3 (1 + ξ1 )(1 + ξ2 )(1 − ξ3 ) + P4 (1 − ξ1 )(1 + ξ2 )(1 − ξ3 ) + P5 (1 − ξ1 )(1 − ξ2 )(1 + ξ3 ) + P6 (1 + ξ1 )(1 − ξ2 )(1 + ξ3 ) + P7 (1 + ξ1 )(1 + ξ2 )(1 + ξ3 ) + P8 (1 − ξ1 )(1 + ξ2 )(1 + ξ3 )}, (26) from the reference space ξ = (ξ1 , ξ2 , ξ3 )T to the physical one x = (x1 , x2 , x3 )T , with Pv the element vertices and v ∈ {1, .., 8} [13]. The sub cell volume in the physical space can be obtained by the Jacobian of the map Jij =

∂hi (ξ) . ∂ξj 10

(27)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The volume of each sub cell (around the interpolation points) in the real space Vinter then becomes ref Vinter = |J (ξ) |Vinter .

(28)

Assigning particles to the interpolation points. The assignment of the particles to the nearest interpolation point is carried out in the reference space. This is simply due to the fact that calculating the distances to each interpolation point would be very costly, especially for large Ninter per cell. The drawback of this practice is that particles may not be assigned to their nearest interpolation point in strongly deformed cells. The first step is to map the particle position X into the reference space [15]. However the mapping from the physical space to the reference space, is the inverse of the vector-polynomial given in Eq. (26). Unfortunately, an analytically straight-forward calculation is not available for the inverse vector-polynomial. Let Ξ be the position of the particle in the reference space. Therefore, a Newton method is employed in order to find Ξ iteratively Ξn+1 = Ξn − L−1 (Ξn )F (Ξn ), (29) where

F (Ξ) = h(Ξ) − X,

Lij =

∂Fi . ∂Ξj

(30)

For the initial guess, only the linear content of h(Ξ) is used. With the mapped particle position Ξ ∈ [−1, 1]3 , the assignment to the interpolation point ind = (ind1 , ind2 , ind3 ) is straight-forward   Ξ(i) + 1 indi = INT N(i),inter , (31) 2 where INT(...) returns the integer part of its argument and the summation over i is suppressed by the brackets. 4.2. Adaptation of virtual sub cells By refining the cells into sub cells as described in the last section, the spatial resolution requirements derived in §3.1 can be met. However, the demand on the minimum number of particles per cell as mentioned in §3.3 is not necessarily guaranteed. To cope with that, an additional step has to be taken by 11

merging the neighbouring sub cells up to a point where the minimum number of particles is honored. To proceed, a space filling curve is defined in each cell connecting the sub cells. The virtual sub cells are then merged along the space filling curve until the number of particles in the merged cell meets the minimum requirement. After that, a new merged cell starts and so forth. Based on several computational experiments, the so-called Peano curve is adopted for our purpose. As depicted in Fig. 1, the Peano curve first connects the sub cells in the planar slices and only after the slices are merged along the third dimension. Therefore a favored direction emerges naturally after planar slices are formed. For example in Fig. 1, x-coordinate is the favored one i.e. it will be dissolved only in the last step. The curves then can be exploited more meaningfully if

2

z

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1 0 0 0

1 1 y

2

2

x

Figure 1: Example of the space filling curve for Ninter = 3. The large box represents one computational cell in which the gray sub cells are connected via the red curve.

the largest gradients always lie on the favored direction. However note that the normalized operators ∇Φ defined in §3.1 are averaged over all directions. Therefore besides ∇Φ , the directional temperatures
m 2 Ti = Mi − hMi i2 . (32) k 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

are considered. The favored direction in cell j then can be found by the temperature normalized gradient (j) ∇Ti

=



∂Ti ∂xl

(j)

(j)

∆xl ; T (j)

i ∈ {1, 2, 3},

(33)

where T is the temperature. Consequently the space filling curve inside cell (j) j is constructed such that the direction with the largest ∇Ti is merged at last. Certainly the above definition of the favored direction is not unique. Yet based on analysis of strong shocks and the fact that directional temperatures can differ significantly from the total one, Eq. (33) is adopted as the basis. 4.3. Algorithm of flow field adaptation A flow chart of the algorithm for flow field adaptation is depicted in Fig. 2. The simulation starts with an initial mesh. In addition to the FP solution (j) algorithm described in sec. 2, the normalized operators ∇Φ and gradients (j) ∇Ti are sampled in all cells. Next, each cell with fewer than 12 particles or with a fine enough spatial resolution (according to Eq. (14)) is assigned as non-adaptive. Others will be refined using sub cells as described in §4.1. The number of interpolation points per direction Ni,inter = (Ninter )1/3 is chosen depending on the particle number in the cell np . A trade-off between the high spatial resolution and computational performance leads to minimum 6 particles per each sub cell. Therefore, Ninter per cell can be calculated using h n i p . (34) Ninter = INT 6

Those sub cells which do not meet the minimum 6 particles, are merged according to §4.2 until the minimum is achieved in the merged sub cell. 5. Results All simulations are performed with the coupled Particle-In-Cell - DSMC code PICLas [20]. The code uses unstructured hexahedral meshes and allows high order simulations using the PIC solver. Instead of an adaptive mesh refinement for the DSMC method, a virtual sub cell method using an octree with a 13

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

sample gradient data

yes

calculate Ninter

define virtual sub cells perform normal FP

no

max(∇Φ ) > 0.05? Φ ∈ T, p, v, ρ

adapt virtual sub cells perform FP in each adapted virtual sub cell

Figure 2: Flowchart of the Adaptive FP algorithm.

subsequent nearest neighbor search is implemented to resolve the mean free path [12]. The structure of the code allows highly parallel simulations with DSMC, PIC and the described adaptive FP method. For the parallelization, a domain decomposition is carried out according to [21]. The 70◦ blunted cone described in [3] is chosen to assess the performance of the FP cell adaptation scheme. The simulations were carried out for atomic nitrogen and using the following inflow conditions: v∞ [ms−1 ] T∞ [K] n∞ [m−3 ] 1502.4 13.58 1.115 · 1021 The angle of attack is set to α = 30◦ . Therefore, a full 3D simulation is set up using a particle weighting factor of wp = 2×1010 resulting in around 3.5×107 particles in the stationary condition. The simulations were performed using 1200 cores on the CRAY XC40 system of the high performance computing center Stuttgart. The end time of all simulations is tEnd = 1 × 10−2 s, for this the normal FP method takes around 25 minutes using 1200 cores. The result of max(∇T , ∇p , ∇v , ∇ρ ) and the consequent classification of the normal FP and the adaptive cells is shown in Fig. 3, for the last adaptation step. Steady state results are shown for the gas temperature and the heat flux on the cone using DSMC and the Adaptive-FP method in Fig. 4. Furthermore, the heat flux as well as the force on the cone surface are adopted in order to compare the accuracy of different methods (DSMC, FP and AdaptiveFP). Note that the both values are important for analyzing re-entry vehicles. 14

0.3 0.1

0.25

8 · 10−2

0.2 z [m]

max(∇T , ∇p , ∇v , ∇ρ ) [−]

0.15

6 · 10−2

0.1

4 · 10−2

5 · 10−2

2 · 10−2

0

0

0

2 · 10−2 4 · 10−2 6 · 10−2 8 · 10−2

0.1

x [m]

0.1

8 · 10−2 z [m]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

6 · 10−2 4 · 10−2 2 · 10−2 0

0

2 · 10−2 4 · 10−2 6 · 10−2 8 · 10−2

0.1

x [m]

Figure 3: Top: Result of max(∇T , ∇p , ∇v , ∇ρ ) for the last adaptation step. Bottom: The Adapt-FP algorithm adjusts the red cells, whereas the blue ones remain intact.

Figure 4 shows a white line along the surface that corresponds to the path variable S which is used in the next plots. Additionally, markers {A,B,C,D} are shown to highlight special points on the path S. Figure 5 exhibits obvious differences between the heatflux of the normal FP algorithm and the adaptive FP algorithm. While the former overestimates the DSMC results, the adaptive FP matches DSMC pretty well. Figure 6 shows the force on the surface along x and z directions (along y is almost zero on the path S). Here, the differences among all three methods are small.

15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Isolines of the temperature and the velocity of the flow field are shown in Fig. 7. Here again, the adaptive FP method matches the DSMC results much better than the normal FP method; especially in the shock region in front of the heat shield. Nevertheless, there are visible differences between the FP methods and the DSMC method in the rarefied region behind the heat shield. Concerning the computational effort of the methods, the adaptive FP method needs ≈ 4 times less CPU time compared to DSMC, while the normal FP method only required ≈ 7.5 times less CPU time compared to DSMC. Note that these CPU times are only a rough measure of the computational effort. Especially, the adaptive FP method is in an early stage of development, and additional optimizations are possible and necessary. Moreover, the gap between the FP methods and DSMC in required computational times becomes more pronounced at larger densities [9]. It is also interesting to compare the time step sizes employed by DSMC in comparison to the Adaptive-FP and normal FP methods. The time step size of ∆tDSM C = 5 × 10−8 s was utilized for DSMC in order to resolve the collision frequency, whereas for the both FP methods a much larger value of ∆tF P = 4 × 10−7 s resulting from CFL condition was employed. Note that this time step difference further increases for denser gases. 6. Conclusion and Outlook In this paper a thorough technical analysis of the FP solution algorithm and in particular the cubic model [2] is delivered. First, different prerequisites for spatio-temporal discretizations along with particle number requirements are derived. Furthermore virtual sub cells are defined for unstructured grids such that local spatial refinements become possible by assigning the computational particles to the introduced sub cells. However since that may lead to a number of particles fewer than the minimum allowed, a merging sub cell algorithm is devised. By using Peano space fitting curves, the particle number requirement is met in each merged sub cell. Therefore an AdaptiveFP algorithm is outlined by incorporating the refinement and merging ideas. For assessment, the 70◦ blunted cone is investigated by DSMC, the regular FP and the Adaptive-FP method. The results of the heat flux and the total force on the cone surface demonstrate a very good agreement between Adaptive-FP and DSMC. In the future works more complex geometries such 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

as double cone structure will be studied. Moreover, further validations for hard sphere as well as diatomic gases will be carried out. 7. Acknowledgements The authors gratefully thank the Deutsche Forschungsgemeinschaft (DFG) for funding this research within the project “Kinetic Algorithms for the Maxwell-Boltzmann System and the Simulation of Magnetospheric Propulsion Systems”. We thank the High Performance Computing Center (HLRS) for granting the computational time that has allowed the execution of the presented simulations. References [1] P. Jenny, M. Torrilhon, S. Heinz, A solution algorithm for the fluid dynamic equations based on a stochastic model for molecular motion, Journal of computational physics 229 (4) (2010) 1077–1098. [2] M. H. Gorji, M. Torrilhon, P. Jenny, Fokker–Planck model for computational studies of monatomic rarefied gas flows, Journal of fluid mechanics 680 (2011) 574–601. [3] J. Allgre, D. Bisch, J. Lengrand, Experimental Rarefied Heat Transfer at Hypersonic Conditions over 70-Degree Blunted Cone, Journal of Spacecraft and Rockets 34 (6) (1997) 724–728. [4] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Clarendon Press, Oxford, 1994. [5] M. H. Gorji, P. Jenny, An efficient particle Fokker–Planck algorithm for rarefied gas flows, Journal of Computational Physics 262 (2014) 325–343. [6] S. Heinz, Molecular to fluid dynamics: the consequences of stochastic molecular motion, Physical Review E. 70 (3). [7] F. Fei, J. Fan, A diffusive information preservation method for small Knudsen number flows, Journal of Computational Physics 243 (2013) 179–193.

17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[8] J. Mathiaud, L. Mieussens, A FokkerPlanck Model of the Boltzmann Equation with Correct Prandtl Number, Journal of Statistical Physics 162 (2016) 397–414. [9] M. H. Gorji, P. Jenny, Fokker-Planck-DSMC algorithm for simulations of rarefied gas flows, Journal of Computational Physics 287 (2015) 110 – 129. [10] G. A. Bird, M. A. Gallis, J. R. Torczynski, D. J. Rader, Accuracy and efficiency of the sophisticated direct simulation Monte Carlo algorithm for simulating noncontinuum gas flows, Physics of Fluids 21 (2009) 1–12. [11] J. S. Wu, K. C. Tseng, F. Y. Wu, Parallel three-dimensional dsmc method using mesh refinement and variable time-step scheme, Computer Physics Communications 162 (3) (2004) 166 – 187. [12] M. Pfeiffer, A. Mirza, S. Fasoulas, A grid-independent particle pairing strategy for dsmc, Journal of Computational Physics 246 (2013) 28 – 36. [13] D. A. Kopriva, Implementing Spectral Methods for Partial Differential Equations, Springer, 2009. [14] F. Hindenlang, G. J. Gassner, C. Altmann, A. Beck, M. Staudenmaier, C.-D. Munz, Explicit discontinuous Galerkin methods for unsteady problems, Computers & Fluids 61 (2012) 69–93. [15] M. Pfeiffer, A. Mirza, C.-D. Munz, S. Fasoulas, Two statistical particle split and merge methods for Particle-in-Cell codes , Computer Physics Communications 191 (0) (2015) 9–24. [16] C. W. Gardiner, Handbook of Stochastic Methods, Berlin-Springer, 1985. [17] W. G. Habashi, J. Dompierre, Y. Bourgault, D. Ait-Ali-Yahia, M. Fortin, M.-G. Vallet, Anisotropic mesh adaptation: towards userindependent, mesh-independent and solver-independent CFD. Part I: general principles, International Journal for Numerical Methods in Fluids 32 (6) (2000) 725–744.

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[18] K. L. Bibb, P. A. Gnoffo, M. A. Park, W. T. Jones, Parallel, gradientbased anisotropic mesh adaptation for re-entry vehicle configurations, AIAA Paper 3579 (2006) 2006. [19] K. Kannenberg, I. D. Boyd, Strategies for efficient particle resolution in the direct simulation monte carlo method, Journal of Computational Physics 152 (2) (2000) 727–745. [20] C.-D. Munz, M. Auweter-Kurtz, S. Fasoulas, A. Mirza, P. Ortwein, M. Pfeiffer, T. Stindl, Coupled particle-in-cell and direct simulation monte carlo method for simulating reactive plasma flows, Comptes Rendus Mcanique 342 (1011) (2014) 662 – 670, theoretical and numerical approaches for Vlasov-maxwell equations. [21] P. Ortwein, T. Binder, S. Copplestone, A. Mirza, P. Nizenkov, M. Pfeiffer, T. Stindl, S. Fasoulas, C.-D. Munz, Parallel Performance of a Discontinuous Galerkin Spectral Element Method Based PIC-DSMC Solver, Springer International Publishing, Cham, 2015, pp. 671–681.

19

Temperature [K]

A

  Heatflux kWm−2

B

D

Temperature [K]

C

Velocity [m/s]

  Heatflux kWm−2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 4: Top: The temperature and the surface heat flux contour plots by DSMC, the white path S and the positions of the markers {A,B,C,D}. Bottom: The temperature and the surface heat flux contour plots by the adaptive FP method. Additionally, the velocity stream lines and magnitude are shown in gray.

20

C

  Heat flux kWm−2

30 20

DSMC FP-Adapt FP

B A

10

D

0 0

2

4 6 S/Rn [-]

8

10

100

Force Z per Area [Nm−2 ]

Figure 5: Comparison of the heat flux along the path S of 70◦ blunted cone with markers {A,B,C,D}. The path S is normalized with the nose radius Rn = 0.0125 m of the cone.

Force X per Area [Nm−2 ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

DSMC FP-Adapt FP

B

50 C

A

D

0 0

2

4 6 S/Rn [-]

8

10

C

40 20

B

0

A

0

DSMC FP-Adapt FP

D

2

4 6 S/Rn [-]

8

10

Figure 6: Comparison of the force per area along the path S of 70◦ blunted cone with markers {A,B,C,D}. Up: The force per area along x direction, Bottom: The force per area along z direction.

21

DSMC FP-Adapt FP

0.1

0.1 8 · 10−2 z [m]

8 · 10−2 z [m]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

6 · 10−2

6 · 10−2

4 · 10−2

4 · 10−2

2 · 10−2

2 · 10−2

0

0

2 · 10−24 · 10−26 · 10−28 · 10−2 0.1 x [m]

0

DSMC FP-Adapt FP 0

2 · 10−24 · 10−26 · 10−28 · 10−2 0.1 x [m]

Figure 7: Isoline comparisons on a flow field slice at position y = 0. Left: the temperature, Right: the velocity magnitude.

22