A conservative discrete velocity method for the ellipsoidal Fokker-Planck equation in gas-kinetic theory

TL;DR

This paper introduces a conservative discrete velocity method for the ellipsoidal Fokker-Planck equation, improving accuracy and efficiency in simulating non-equilibrium gas flows, especially at high Mach numbers.

Contribution

A novel conservative discrete velocity method with a finite volume framework and quasi reduction for efficient, accurate simulation of non-equilibrium gas flows in the ellipsoidal Fokker-Planck model.

Findings

Method accurately predicts shock structures at Mach 8.

Conservation properties are maintained in the numerical scheme.

Predictions agree well with benchmark solutions and DSMC results.

Abstract

A conservative discrete velocity method (DVM) is developed for the ellipsoidal Fokker-Planck (ES-FP) equation in prediction of non-equilibrium neutral gas flows in this paper. The ES-FP collision operator is solved in discrete velocity space in a concise and quick finite difference framework. The conservation problem of discrete ES-FP collision operator is solved by multiplying each term in it by extra conservative coefficients whose values are very closed to unity. Their differences to unity are in the same order of the numerical error in approximating the ES-FP operator in discrete velocity space. All the macroscopic conservative variables (mass, momentum and energy) are conserved in the present modified discrete ES-FP collision operator. Since the conservation property in discrete element of physical space is very important for numerical scheme when discontinuity and large gradient exist in flow field, a finite volume framework is adopted for the transport term of ES-FP equation. For $nD$-$3V$ ($n<3$) cases, a $nD$-quasi $nV$ reduction is specially proposed for ES-FP equation and the corresponding FP-DVM method, which can greatly reduce the computational cost. The validity and accuracy of both ES-FP equation and FP-DVM method are examined using a series of $0D$-$3V$ homogenous relaxation cases and $1D$-$3V$ shock structure cases with different $Mach$ numbers, in which $1D$-$3V$ cases are reduced to $1D$-quasi $1V$ cases. Both the predictions of $0D$-$3V$ and $1D$-$3V$ cases match well with the benchmark results such as analytical Boltzmann solution, direct full-Boltzmann numerical solution and DSMC result. Especially, the FP-DVM predictions match well with the DSMC results in the $Mach$ 8.0 shock structure case, which is in high non-equilibrium, and is a challenge case of the model Boltzmann equation and the corresponding numerical methods.

Tables

Table 1. Table 1: macroscopic variables and ε 𝜀 \varepsilon at different settings.

	${50}^3$ mesh 2nd order	${50}^3$ mesh 4th order	${100}^3$ mesh 2nd order
	(base setting)	(high order)	(dense mesh)
$\left|\rho -1\right|$	8.5e-5	8.5e-5	1.3e-5
$\left|u_1\right|$	5.0e-5	5.0e-5	4.7e-6
$\left|T-1\right|$	5.2e-4	5.2e-4	9.5e-5
$\left|{\epsilon }_f-1\right|$	7.7e-5	4.6e-5	4.5e-6
$\left|{\epsilon }_{A,1}-1\right|$	3.0e-5	7.4e-5	1.1e-6
$\left|{\epsilon }_D-1\right|$	2.0e-2	6.4e-4	5.0e-3
$\log \left|{\epsilon }_f-1\right|/\log \Delta \xi$	5.8	6.2	5.3
$\log \left|{\epsilon }_{A,1}-1\right|/\log \Delta \xi$	6.5	5.9	6.0
$\log \left|{\epsilon }_D-1\right|/\log \Delta \xi$	2.4	4.6	2.3

Equations

$$\frac{\partial f}{\partial t}+{\xi }_i\frac{\partial f}{\partial x_i}+a_i\frac{\partial f}{\partial {\xi }_i}=C\left(f,f\right),$$

$$\rho =⟨mf⟩=mn,$$

$$\rho u_i=⟨m{\xi }_{i}f⟩,$$

$$\rho e=⟨\frac{1}{2}m{\xi }_k{\xi }_{k}f⟩=\frac{1}{2}\rho u_k{u}_k+\frac{3}{2}nkT,$$

$${\tau }_{ij}=-⟨m\left(c_i{c}_j-\frac{1}{3}{c}_k{c}_k{\delta }_{ij}\right)f⟩=-nk\left(T_{ij}-T{\delta }_{ij}\right),$$

$$q_i=⟨\frac{1}{2}m{c}_i{c}_k{c}_k⟩,$$

$$⟨\cdot ⟩={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\left(\cdot \right)d{\xi }_{1}d{\xi }_{2}d{\xi }_3.$$

$$\frac{\partial f}{\partial t}+{\xi }_i\frac{\partial f}{\partial x_i}=\frac{1}{{\tau }_{FP}}\left\{\frac{\partial \left(\left({\xi }_i-u_i\right)f\right)}{\partial {\xi }_i}+RT{\delta }_{ij}\frac{{\partial }^{2}f}{\partial {\xi }_i\partial {\xi }_j}\right\},$$

$$T_{ES,ij}=\left(1-\nu \right)T{\delta }_{ij}+\nu T_{ij},$$

$$\nu =\max \left(-\frac{5}{4},-\frac{T}{{\lambda }_{\max }-T}\right),$$

$$\mathrm{Pr}=\frac{3}{2\left(1-\nu \right)},$$

$${\tau }_{ES}=2\left(1-\nu \right)\frac{\mu }{p}.$$

$$\frac{\partial f}{\partial t}+{\xi }_i\frac{\partial f}{\partial x_i}=\frac{1}{{\tau }_{ES}}\left\{\frac{\partial \left(\left({\xi }_i-u_i\right)f\right)}{\partial {\xi }_i}+R{T}_{ES,ij}\frac{{\partial }^{2}f}{\partial {\xi }_i\partial {\xi }_j}\right\}.$$

$$\frac{f^{\ast }-f^n}{\Delta t}+\frac{1}{\Omega }\sum _{a=1}^A\left({\xi }_i{f}_a\right)S_{a,i}=0,$$

$$f_a=f^n+\frac{\partial f^n}{\partial x_j}\left(x_{a,j}-x_{c,j}\right).$$

$$\frac{f^{n+1}-f^{\ast }}{\Delta t}=\frac{1}{{\tau }_{ES}^{\ast }}\left\{\frac{\partial \left(\left({\xi }_i-u_i^{\ast }\right)f^{\ast }\right)}{\partial {\xi }_i}+R{T}_{ES,ij}^{\ast }\frac{{\partial }^2{f}^{\ast }}{\partial {\xi }_i\partial {\xi }_j}\right\},$$

$$\frac{\partial f}{\partial t}=\frac{1}{{\tau }_{ES}}\left\{3f+\left({\xi }_i-u_i\right)\frac{\partial f}{\partial {\xi }_i}+R{T}_{ES,ij}\frac{{\partial }^{2}f}{\partial {\xi }_i\partial {\xi }_j}\right\}.$$

$$\frac{\partial f}{\partial t}=\frac{1}{{\tau }_{ES}}\left\{3{\epsilon }_{F}f+{\epsilon }_{A,i}\left({\xi }_i-u_i\right)\frac{\partial f}{\partial {\xi }_i}+{\epsilon }_{D}R{T}_{ES,ij}\frac{{\partial }^{2}f}{\partial {\xi }_i\partial {\xi }_j}\right\}.$$

$$\left\{{\begin{array}[]{*{20}{c}}{{\varepsilon_{A,i}}\neq{\rm{1}},{\rm{}}{\xi_{i}}<0},\\ {{\varepsilon_{A,i}}{\rm{=1}},{\rm{}}{\xi_{i}}\geq 0}.\end{array}}\right.$$

$${\epsilon }_F\sum 3f+{\epsilon }_{A,1}\sum _{{\xi }_1<0}{A}_1+{\epsilon }_{A,2}\sum _{{\xi }_2<0}{A}_2+{\epsilon }_{A,3}\sum _{{\xi }_3<0}{A}_3+{\epsilon }_D\sum D=-\left(\sum _{{\xi }_1\ge 0}{A}_1+\sum _{{\xi }_2\ge 0}{A}_2+\sum _{{\xi }_3\ge 0}{A}_3\right),$$

$${\epsilon }_F\sum 3{\xi }_{1}f+{\epsilon }_{A,1}\sum _{{\xi }_1<0}{\xi }_1{A}_1+{\epsilon }_{A,2}\sum _{{\xi }_2<0}{\xi }_1{A}_2+{\epsilon }_{A,3}\sum _{{\xi }_3<0}{\xi }_1{A}_3+{\epsilon }_D\sum {\xi }_{1}D=-\left(\sum _{{\xi }_1\ge 0}{\xi }_1{A}_1+\sum _{{\xi }_2\ge 0}{\xi }_1{A}_2+\sum _{{\xi }_3\ge 0}{\xi }_1{A}_3\right),$$

$${\epsilon }_F\sum 3{\xi }_{2}f+{\epsilon }_{A,1}\sum _{{\xi }_1<0}{\xi }_2{A}_1+{\epsilon }_{A,2}\sum _{{\xi }_2<0}{\xi }_2{A}_2+{\epsilon }_{A,3}\sum _{{\xi }_3<0}{\xi }_2{A}_3+{\epsilon }_D\sum {\xi }_{2}D=-\left(\sum _{{\xi }_1\ge 0}{\xi }_2{A}_1+\sum _{{\xi }_2\ge 0}{\xi }_2{A}_2+\sum _{{\xi }_3\ge 0}{\xi }_2{A}_3\right),$$

$${\epsilon }_F\sum 3{\xi }_{3}f+{\epsilon }_{A,1}\sum _{{\xi }_1<0}{\xi }_3{A}_1+{\epsilon }_{A,2}\sum _{{\xi }_2<0}{\xi }_3{A}_2+{\epsilon }_{A,3}\sum _{{\xi }_3<0}{\xi }_3{A}_3+{\epsilon }_D\sum {\xi }_{3}D=-\left(\sum _{{\xi }_1\ge 0}{\xi }_3{A}_1+\sum _{{\xi }_2\ge 0}{\xi }_3{A}_2+\sum _{{\xi }_3\ge 0}{\xi }_3{A}_3\right),$$

$${\epsilon }_F\sum 3{\xi }^{2}f+{\epsilon }_{A,1}\sum _{{\xi }_1<0}{\xi }^2{A}_1+{\epsilon }_{A,2}\sum _{{\xi }_2<0}{\xi }^2{A}_2+{\epsilon }_{A,3}\sum _{{\xi }_3<0}{\xi }^2{A}_3+{\epsilon }_D\sum {\xi }^{2}D=-\left(\sum _{{\xi }_1\ge 0}{\xi }^2{A}_1+\sum _{{\xi }_2\ge 0}{\xi }^2{A}_2+\sum _{{\xi }_3\ge 0}{\xi }^2{A}_3\right),$$

$$A_i=\left({\xi }_i-u_i\right)\frac{\partial f}{\partial {\xi }_i},$$

$$D=R{T}_{ES,ij}\frac{{\partial }^{2}f}{\partial {\xi }_i\partial {\xi }_j}.$$

$$\sum \left(\cdot \right)=\sum _{{\xi }_1\in \left({\xi }_{1,\min },{\xi }_{1,\max }\right)}\sum _{{\xi }_2\in \left({\xi }_{2,\min },{\xi }_{2,\max }\right)}\sum _{{\xi }_3\in \left({\xi }_{3,\min },{\xi }_{3,\max }\right)}\left(\cdot \right)\Delta {\xi }_1\Delta {\xi }_2\Delta {\xi }_3,$$

$$\sum _{{\xi }_1<0}\left(\cdot \right)=\sum _{{\xi }_1\in \left({\xi }_{1,\min },0\right)}\sum _{{\xi }_2\in \left({\xi }_{2,\min },{\xi }_{2,\max }\right)}\sum _{{\xi }_3\in \left({\xi }_{3,\min },{\xi }_{3,\max }\right)}\left(\cdot \right)\Delta {\xi }_1\Delta {\xi }_2\Delta {\xi }_3,$$

$$\sum _{{\xi }_1\ge 0}\left(\cdot \right)=\sum _{{\xi }_1\in [0,{\xi }_{1,\max })}\sum _{{\xi }_2\in \left({\xi }_{2,\min },{\xi }_{2,\max }\right)}\sum _{{\xi }_3\in \left({\xi }_{3,\min },{\xi }_{3,\max }\right)}\left(\cdot \right)\Delta {\xi }_1\Delta {\xi }_2\Delta {\xi }_3.$$

$$\frac{f^{n+1}-f^{\ast }}{\Delta t}=\frac{1}{{\tau }_{ES}^{\ast }}\left\{3{\epsilon }_F^{\ast }{f}^{\ast }+{\epsilon }_{A,i}^{\ast }\left({\xi }_i-u_i^{\ast }\right)\frac{\partial f^{\ast }}{\partial {\xi }_i}+{\epsilon }_F^{\ast }R{T}_{ES,ij}^{\ast }\frac{{\partial }^2{f}^{\ast }}{\partial {\xi }_i\partial {\xi }_j}\right\}.$$

$$T_{ij}^{\ast }=\frac{\sum c_i{c}_j{f}^{\ast }}{R\sum f^{\ast }},$$

$$u_i^{\ast }=\frac{\sum {\xi }_i{f}^{\ast }}{\sum f^{\ast }}.$$

$$\frac{\partial f}{\partial x_2}=0,\frac{\partial f}{\partial x_3}=0,$$

$$u_2=0,u_3=0,$$

$$T_{12}=T_{21}=0,T_{23}=T_{32}=0,T_{13}=T_{31}=0,$$

$$T_{22}=T_{33}=\left(3T-T_{11}\right)/2.$$

$$\frac{\partial f}{\partial t}+{\xi }_1\frac{\partial f}{\partial x_1}=\frac{1}{{\tau }_{ES}}\left\{3f+\left({\xi }_1-u_1\right)\frac{\partial f}{\partial {\xi }_1}+{\xi }_2\frac{\partial f}{\partial {\xi }_2}+{\xi }_3\frac{\partial f}{\partial {\xi }_3}+R{T}_{ES,11}\frac{{\partial }^{2}f}{\partial {\xi }_1^2}+R{T}_{ES,22}\frac{{\partial }^{2}f}{\partial {\xi }_2^2}+R{T}_{ES,33}\frac{{\partial }^{2}f}{\partial {\xi }_3^2}\right\}.$$

$$F={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }mfd{\xi }_{2}d{\xi }_3,$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

The Fokker-Planck (FP) equation with advection-diffusion collision operator is widely used in modeling dynamic systems such as neutral molecule [1, 2, 3], plasma [4, 5, 6], photonics [7, 8], and even biological [9], economic [10], and social [11] systems. The first FP equation for molecule system is derived from Boltzmann equation in gas kinetic theory when counting the gazing effect of molecule collisions [12]. The Prandtl ($Pr$) number yielded from this FP equation is fixed at $3/2$. Since the $Pr$ number of real gas is below unity ($2/3$ for monatomic gas), two types of modified FP equations are proposed. They are cubic-FP equation [2] and ellipsoidal-FP (ES-FP) equation [3]. Recently, by mapping these FP equations to Stochastic Differential Equations (SDEs), the FP equations are solved in a stochastic and particle way [1, 2]. Comparing with other particle methods such as Direct Simulation Monte-Carlo (DSMC) [13], its computational cost is greatly reduced in continuum limit (dissipation limit in FP research). Since the mechanism of drag and diffusion forces holds for both micro and macro scales, then large time step and cell length can be used. This advantage is very important for the prediction of flow fields in either transitional or continuum regimes where the molecular mean free path (m.f.p.) and mean collision time (m.c.t.) are comparable or greatly less than the characteristic length and time, respectively. On the other hand, since the deterministic methods are not affected by statistical fluctuation, they are very useful in the precise computation of multi-scale non-equilibrium flows, and are helpful in investigating the mechanism of such flows. Moreover, the accuracy of modified FP equations, especially their collision operators, should be examined using a deterministic method.

The first deterministic numerical scheme for FP equation, which can preserve equilibrium, is proposed for a $1V$ (one dimensional velocity space) isotropic linear Fokker-Planck-Landau (FPL) system [14]. Then it is extended to mass/energy/equilibrium preserving scheme [15], $2V$ cases in discrete cylindrical velocity space [16], and nonlinear equation [17, 18]. Besides the FPL-type, another form of FP equation is Rosenbluth-Fokker-Planck equation (RFP) [4]. RFP equation has a similar mathematical form of the FP equation that is derived from Boltzmann equation [12]. Its collision operator is written in a differential form with nonlinear advection and diffusion coefficients. By using a finite volume framework in velocity space and extra coefficients on advection terms for conservation purpose, RFP equation is well solved by the deterministic numerical method in Ref. [19, 20].

Plenty of the FP researches focus on the homogenous FP equation where the particle transport in physical space is assumed to be zero [21]. Several works address the non-homogenous FP equation in the topic of FPL-type equation [22, 23]. In flow predictions, the particle transport in physical space is an essential aspect that can not be ignored. The transport term in the FP equation of gas-kinetic theory (similar to RFP-type) should be calculated. By taking the particle transport into consideration, the deterministic numerical framework will be the discrete velocity method (DVM). There are several multi-scale methods in gas kinetic theory, such as Unified Gas-Kinetic Scheme (UGKS) [24], Discrete Unified Gas-Kinetic Scheme (DUGKS) [25], Gas-Kinetic Unified Algorithm (GKUA) [26], using the DVM framework for predictions of flows from continuum regime to rarefied regime. Recently, these methods have been widely used in the prediction of non-equilibrium flows [24], plasmas [27], and photonics [28]. Since they use Bhatnagar-Gross-Krook (BGK)-type model equations, certain degree of deviation (from the Boltzmann equation) exists in the prediction of high non-equilibrium flows [29], such as the shock structure cases with high $Mach$ numbers, which will be calculated in the later section.

In this paper, a novel deterministic method (FP-DVM) is proposed for non-equilibrium flows, which solves the ES-FP equation numerically in the DVM framework. In FP-DVM, the ES-FP collision operator is solved in a deterministic way in discrete velocity space using a framework of Finite Difference Method (FDM). Instead of considering the conservation in discrete element of velocity space, the discrete ES-FP collision operator is treated in a relaxation way, and the conservations of mass, momentum and energy are ensured by coefficients being added to both advection and dissipation terms. The deviations of their values from unity have the same order as the truncation error of the numerical approximation to ES-FP collision operator. The computational complexity of solving the discrete ES-FP collision operator is $O(N)$, here “$N$” is the number of discrete points in velocity space. Four $0D$-$3V$ ($nD$-$mV$ is a denotation of cases whose dimension of physical space is “$n$” and the dimension of velocity space is “$m$” in the scope of FP research) homogenous relaxation cases and three $1D$-$3V$ shock structure cases with different $Mach$ numbers are conducted to examine the validity and accuracy of both ES-FP equation and the present numerical method. Using the dimensional reduction method for ES-FP equation proposed in this paper, $1D$-$3V$ cases are reduced to $1D$-quasi $1V$ cases, and the computational cost is greatly reduced. The remaining of this paper is organized as follows: Sec. 2 is a quick review of gas-kinetic theory and FP equation. Sec. 3 is the construction of FP-DVM method; Sec. 4 is the numerical experiment; Sec. 5 is the discussion and conclusion.

2 Gas kinetic theory and Fokker-Planck equation

2.1 Distribution function and Boltzmann equation

In gas kinetic theory, molecular system is described using distribution function $f(x_{i},\xi_{j},t)$ depending on location $x_{i}$, molecular velocity $\xi_{j}$ and time $t$. It is the number density of molecules that arrived at $x_{i}$ at time $t$ with velocity $\xi_{j}$. For dilute gas, the evolution of $f$ is governed by Boltzmann equation [30]:

[TABLE]

where $a_{i}$ is acceleration of molecule due to body force such as gravity. Einstein summation convention is used throughout this paper if without special statement. The Left Hand Side (LHS) of Eq. 1 is the free transport operator, while the Right Hand Side (RHS) is the collision operator which is mathematically a five-fold nonlinear integral.

Given the distribution function $f$, macroscopic physical variables, such as mass density $\rho$, momentum density $\rho u_{i}$ ($u_{i}$ is macroscopic velocity), energy density $\rho e$ ($e$ is energy per mass), stress $\tau_{ij}$ and heat flux $q_{i}$, can be obtained using their definition in gas-kinetic theory as follows,

[TABLE]

where $n$ is number density, $c_{i}$ is the peculiar velocity defined as $\xi_{i}-u_{i}$, $T$ is thermodynamic temperature, $T_{ij}$ is the temperature tensor in gas kinetic theory whose trace is $3T$, $k$ is Boltzmann constant, $m$ is the mass of molecule, $\delta_{ij}$ is the Kronecker delta, the operator $\left\langle\cdot\right\rangle$ is an integral over the whole velocity space, which can be written as

[TABLE]

2.2 Fokker-Planck equation for gas kinetic theory

In the scope of gas kinetic theory, the original Fokker-Planck equation without body force is derived in Ref. [12], where the grazing effect of binary collisions is considered. By changing its relaxation rate to $\tau_{FP}=2\mu/p$ ($\mu$ is viscosity, $p$ is pressure), the standard Fokker-Planck equation as a model of the Boltzmann equation can be written as follows,

[TABLE]

where $R=k/m$ is the specific gas constant. Since the standard Fokker-Planck equation corresponds to a fix Prandtl number of $3/2$, two types of modified Fokker-Planck equations, the cubic-FP equation equation [2] and ES-FP equation [3], are proposed. In cubic-FP, the advection term is multiplied by a polynomial of molecular velocity $\xi_{i}$, whose coefficients are used to get the right relaxation rate of both stress and heat flux, thus a right Prandtl is realized. In ES-FP equation, the diagonal dissipation coefficient $RT\delta_{ij}$ in the standard FP equation is replaced by $T_{ES,ij}$ which is defined as follows

[TABLE]

as a combination of isotropic temperature $T$ and anisotropic temperature $T_{ij}$ ($T_{ij}=\left\langle{m{c_{i}}{c_{j}}f}\right\rangle/\rho R$), and $\nu$ is defined as

[TABLE]

where $\lambda_{max}$ is the maximum eigenvalue of the positive definite matrix $T_{ij}$. Since $\nu$ and $Pr$ number have the following relation

[TABLE]

the $Pr$ number is $2/3$ except in the extreme condition $\lambda_{\max}>1.8T$. In this extreme condition, $Pr$ varies from $2/3$ to unity. The relaxation time $\tau_{ES}$ in ES-FP collision operator is defined as

[TABLE]

Since the procedure of deriving ES-FP equation from standard FP equation is similar to that of extending BGK equation to ES-BGK equation, it is called ES-FP in Ref. [3]. Similar to standard FP, ES-FP is written in the following form

[TABLE]

3 Deterministic discrete velocity method for ellipsoidal Fokker-Planck equation

In the DVM framework, the physical space $x_{i}$, the velocity space $\xi_{i}$ and the time $t$ are discrete. The ES-FP equation is solved in an operator splitting way. The free transport operator (LHS of Eq. 9) is solved in the physical space first, in order to get the distribution $f^{*}$ at intermediate step in each discrete element in physical space. Given $f^{*}$, the collision operator (RHS of Eq. 9) can be solved in the discrete velocity space, then the distribution can be evolved to the next time step.

3.1 Free transport operator

For a numerical scheme in flow predictions, the conservation property in discrete cell (in physical space) is very important when discontinuity, such as the shock wave, exists in the flow field. So, the transport operator of ES-FP equation is solved in a finite volume framework where the extra numerical viscosity needed by capturing the discontinuity is provided by the slope limiters. In this paper, a Euler method is used for temporal discretization. Second order upwind reconstruction in physical space is used for the flux term. The FVM-type numerical scheme for transport operator can be written as

[TABLE]

where

[TABLE]

In Eq. 10, $S_{a,i}$ is cell interface whose direction is from inside to outside. Its subscript “$a$” is an index of discrete cell interface, and the total number of discrete interfaces in a cell is denoted by “$A$”. $\Omega$ is the volume of cell. $\Delta t$ is the time interval. The superscript “$n$” denotes the n-th iteration step, and “$*$” denotes the intermediate time step between the $n$th and $\left(n+1\right)$th steps in the operator splitting treatment. In Eq. 11, the subscript “$c$” denotes the “cell center”. In the present method, the calculation of slope $\partial f/\partial x_{j}$ is to the second order, and van Leer slope limiter is used.

3.2 Collision operator

The collision operator is solved in a finite difference framework, since it is computational efficient. Theoretically, the evolution equation of collision operator can be directly written using the information at intermediate time step as follows

[TABLE]

where the first and second order slopes in velocity space can be numerically approximated using central difference.

The numerical approximation to the slopes has truncation error related to $\Delta\xi$. If a second order central difference is used, the truncation error is $O\left(\Delta{\xi}^{2}\right)$. The numerical quadratures in velocity space also generate numerical errors. If the above numerical scheme for ES-FP collision operator is directly used without treatment of these numerical errors, aggregate effect will produce undesired variations of mass, momentum and energy which should be zero since the collision operator fulfills the conservation property. As a result, it often leads to a non-convergent and non-conservative numerical scheme. This problem is addressed in several works in the topic of FPL equation [6, 31]. For RFP equation which has a similar mathematical form as the FP and ES-FP equations used in gas-kinetic theory, Ref [19] constructed a conservative finite volume scheme in velocity space.

In this paper, the non-convergence problem for efficient finite difference framework is handled in a similar way as Ref [19] as follows. Because in a finite difference framework in velocity space, the distribution function only lives at discrete nodes, the corresponding discrete collision operator should be slightly different from the continuous one due to the inevitable numerical errors. First, the collision operator is rewritten as follows by decomposing the advection term into a distribution function and a first order slope,

[TABLE]

Then each term in the brace is multiplied by a coefficient $\varepsilon$ which is designed to eliminate the influence of numerical errors on conservation property, and the collision operator turns into

[TABLE]

The values of $\varepsilon_{F}$, $\varepsilon_{A,i}$ and $\varepsilon_{D}$ are very close to unity, and their departures from unity (denoted by $\left|{\varepsilon-1}\right|$) are directly related to the numerical error. Being the similar with Ref. [19], the treatment of $\varepsilon_{A,i}$ is as follows,

[TABLE]

That means that $\varepsilon_{A,i}$ only exerts on half of the velocity space. For continuous velocity space, these $\varepsilon$ become unity, since it is the basic property of ES-FP equation that mass, momentum and energy conservations are fulfilled. For discrete velocity space, these coefficients can be obtained using the conservation of mass, momentum, and energy, and solving the following algebraic equations analytically,

[TABLE]

where $A_{i}$ ($i=1,2,3$) and $D$ are defined as

[TABLE]

In Eq.16, the numerical integrals in whole velocity space and half velocity space are denoted by $\sum{\left(\cdot\right)}$, $\sum\limits_{{\xi_{i}}<0}{\left(\cdot\right)}$, and $\sum\limits_{{\xi_{i}}\geq 0}{\left(\cdot\right)}$ respectively. Take index $i=1$ for example, these integrals are in the following form

[TABLE]

where the subscript “$min$” and “$max$” denote the boundaries in each direction of truncated velocity space.

Finally, the evolution of distribution function according to discrete ES-FP collision operator can be written as

[TABLE]

The calculation process of collision operator can be summed up as follows. First, using the information $f^{*}$ at the intermediate step, the first and second order slopes (in velocity space) in Eq. 19 can be calculated and stored. Then the contribution of each term in the brace of Eq. 19 to mass, momentum, and energy can be calculated using numerical integration in Eq. 18. In the process of numerical integration, $u_{i}$ and $T_{ij}$ can also be obtained as

[TABLE]

Then, using the obtained numerical integrals, the coefficient $\varepsilon^{*}$ can be calculated using Eq. 16. Up to this point, every term in Eq. 19 is obtained, and the distribution function can be updated to the $\left(n+1\right)$th time step.

3.3 Reduced ES-FP equation

Real monatomic gas flows have a three dimensional physical space and a three dimensional velocity space ($3D$-$3V$ case). For $nD$-$3V$ case where $n<3$, the ES-FP equation can be reduced to $nD$-quasi $nV$ case whose computational cost is greatly reduced. The following context takes the reducing process from $1D$-$3V$ case to $1D$-quasi $1V$ case for example (which is also used in the case of shock wave structure calculation in this paper). For $1D$ case in $x_{1}$ direction, there is

[TABLE]

The slopes in $x_{2}$ and $x_{3}$ directions are zero. The off-diagonal elements in temperature tensor $T_{ij}$ are zero because tangential stress is zero (according to Eq. 2).

ES-FP equation can be first reduced to

[TABLE]

Define a mass distribution and an energy distribution in $\xi_{1}$ axis as follows

[TABLE]

Then multiply ES-FP equation by unity and $\xi_{2}^{2}+\xi_{3}^{2}$, and integrate it in both $\xi_{2}$ and $\xi_{3}$ directions. After calculating the integrals, the reduced ES-FP equation becomes a system of two equations about “$F$” and “$G$” as follows

[TABLE]

Then the solving of $f$ in three dimensional velocity space is turned into the solving of $F$ and $G$ in one dimensional velocity space. Here the evolutions of $F$ and $G$ are coupled through a relaxation process from $2G$ to $2R\left({{T_{ES,22}}+{T_{ES,{\rm{33}}}}}\right)F$ in the second equation ($G$ equation) in Eq. 24. If the integrated equation is Eq. 14 (the discrete form), $\varepsilon$ will appear in the corresponding terms in both $F$ and $G$ equations in Eq. 24. Practically, these $\varepsilon$ can only appear in $F$ equation, since it can be seen from later equation (Eq. 27) that both mass and momentum conservations are only involved by $F$, and $F$ also appears in the expression of energy conservation. So the reduced ES-FP equations in discrete velocity space can be written as

[TABLE]

here, in $1D$ case, the “$A$” and “$D$” in Eq. 17 are reduced to

[TABLE]

Using conservations of mass, momentum and energy, the $\varepsilon$ for conservation purpose can be obtained by solving the following linear equations,

[TABLE]

here $RHS_{G}$ is the RHS of the $G$ equation (Eq. 25). The numerical process for $1V$ case is the same with the $3V$ case in Sec. 3.1 and Sec. 3.2, expect that the operation of $f$ is now on $F$ and $G$.

3.4 the numerical error in discrete Fokker-Planck collision operator

The numerical error in discrete velocity space comes from the three items below:

the truncation in velocity space, 2. 2.

the error in numerical integration, 3. 3.

the truncation error in calculating the slopes using discrete points.

For the 1st item, the domain of truncated velocity space should be as large as possible. But in order to achieve high computational efficiency, it can not be too large. The Maxwellian distribution suggests that the domain should be at least larger than $3\sqrt{RT}$, since beyond $3\sqrt{RT}$, the distribution only contributes $0.3\%$ of mass.

For the 2nd item, high order numerical integration can be used, such as Newton-Cotes integration, to suppress the numerical error in this item. While for the sake of clarity, the rectangular integration is used in this paper.

For the 3rd item, high order central difference can be used, which will be analyzed in Sec. 4.1. Much of its influence is on the order of the departure of $\varepsilon$ from unity. A too high order central difference will harm the computational efficiency.

4 Numerical experiment

4.1 Maintain the thermal equilibrium Maxwellian distribution ($0D$-$3V$ case)

In this case, the initial distribution function is the Maxwellian distribution $g$ in the following form,

[TABLE]

If a discrete numerical method is conservative, the distribution function will maintain the Maxwellian distribution.

Using this case, the numerical stability of the present method is examined. The validity of conservative coefficients $\varepsilon$ and the accuracy of the present scheme are also investigated. For Maxwellian distribution, $n=1$, $u_{i}=0$, and $T=1$ are chosen. The domain of truncated velocity space in each direction is $[-5,5]$. Both $50*50*50$ and $100*100*100$ meshes in velocity space are tested. Both 2nd and 4rd order central difference are used in approximating the 1st and 2nd order slopes in ES-FP collision operator. The maintained distribution functions that are predicted using $50*50*50$ and $100*100*100$ meshes are shown in radial direction in Fig. 1, respectively. Although the distribution on $50*50*50$ meshes slightly deviates from the analytical Maxwellian distribution near the zero point, it is stable since the discrete ES-FP is conservative. In Table 1, there is a comparison of integral error and $\varepsilon$ under dense/coarse meshes and using low/high order central difference for calculating slopes. It can be seen that integral error ($\left|\rho-1\right|$, $\left|u_{1}\right|$, $\left|T-1\right|$) is related to the mesh number, and almost has no relation to the order of central difference. Comparing with the coarse mesh, by using a dense mesh, the precision of macroscopic variables will increase, while the computational cost will also increase. To increase the integration precision, high order integration method such as Newton-Cotes can be used without using a dense mesh. While for clarity, the direct rectangular integration is used in this paper. Using either a dense mesh or a higher order difference, will make the deviation of $\varepsilon$ from unity a smaller value. The largest deviation comes from $\varepsilon_{D}$. From the $log$ data, it can be seen that when 4rd order difference is used, the order of $\left|\varepsilon_{D}-1\right|$ is about 4, while for a 2nd order difference, it is about 2. This data shows that the order of $\left|\varepsilon_{D}-1\right|$ is related to the order of numerical difference, and is almost not affected by the mesh number. In the following test cases, the second order central difference is used in the velocity space for efficiency.

4.2 Energy relaxation among directions ($0D$-$3V$ case)

In this case, initially the temperatures in different directions are not the same (anisotropic temperature). Through molecular collisions, these temperatures will gradually achieve equilibrium during several m.c.t.. Here the initial temperatures are set to be $T_{1}=2.0$, $T_{2}=1.0$, $T_{3}=1.0$. The truncated discrete velocity space is $[-7,7]$ in each direction with 70 discrete points (707070 mesh). The iteration time $\Delta t$ is chosen as $0.005\tau$. Here $\tau=\mu/p$ has the same order of magnitude as $\tau_{FP}$ and $\tau_{ES}$, but their values are not the same.

For Maxwell molecule, the relaxation of stress and heat flux from the Boltzmann equation is [12]

[TABLE]

For homogenous case, since the density is a constant, the relaxation of anisotropic temperature can be derived from Eq. 29 as follows

[TABLE]

So, the analytical solution of temperature and heat flux can be obtained as

[TABLE]

The relaxation process of distribution in $3V$ space predicted by the present method is shown in Fig. 2, where the iso-surface of distribution gradually transforms from an ellipsoid to a sphere during several $\tau$. In Fig. 3, the relaxation process of anisotropic temperatures predicted by the present method matches precisely with the analytical solution (Eq. 31).

4.3 Relaxation of bi-model distribution function ($0D$-$3V$ case)

In this case, the distribution function is composed of two Maxwellian distributions determined by the physical variables before and after the shock wave respectively. According to the Rankine-Hugoniot relation for a $Mach$ 8.0 shock wave, the physical variables before the shock are $u_{a,1}=8.0$, $u_{a,2}=u_{a,3}=0$, $T_{a}=1.0$, and the physical variables after the shock wave are $u_{b,1}=2.09$, $u_{b,2}=u_{b,3}=0$, $T_{b}=20.87$. The weights of two Maxwellian distributions are chosen as $\rho_{a}=0.9$ and $\rho_{b}=0.1$, in order to mimic the distribution function in the front of shock wave, where high non-equilibrium exists. This case investigates the relaxation of this highly non-equilibrium distribution function. The truncated discrete velocity space is $[-26,26]$ in each direction with 260 discrete points (260260260 mesh). The iteration time $\Delta t$ is chosen as $0.001\tau$. The time evolution of this initial bi-model distribution is shown in Fig. 4, where two Maxwellian distributions merge into a single one during about $10\tau$. The evolutions of anisotropic temperatures and heat flux predicted by the present method are shown in Fig.5, and they match with the analytical solution (Eq. 31) precisely.

4.4 Relaxation of discontinuous distribution function ($0D$-$3V$ case)

The discontinuous distribution function in this case mimics the non-equilibrium distribution at the gas-solid boundary or in Knudsen layer. It is composed of two half Maxwellian distributions. The interface of two half Maxwellian distributions in velocity space is the face $\xi_{1}=0$. Across the interface, the distribution is discontinuous. The Maxwellian distribution on left is determined from $\rho_{a}=1.0$, $u_{a,i}=0$, $T_{a}=2.0$, while the Maxwellian distribution on the right is determined from $\rho_{b}=1.0$, $u_{b,i}=0$, $T_{b}=1.0$. This setting mimics the situation that the temperature of fluid is different from the temperature of the solid wall. The truncated discrete velocity space is $[-8,8]$ in each direction with 80 discrete points (808080 mesh). The iteration time $\Delta t$ is chosen as $0.005\tau$. The time evolution of the initial discontinuous distribution is shown in Fig. 6, where the discontinuity disappears during only one $\tau$, and gradually achieves equilibrium during several $\tau$. The evolutions of anisotropic temperatures and heat flux predicted by the present method are shown in Fig.7. Due to the exponent term in Eq. 31, at first the distribution function approaches the equilibrium in a fast rate. This phenomenon can also be seen from the quick disappearance of discontinuity. Then the rate slows down when the distribution is near equilibrium.

4.5 Normal shock wave structure ($1D$-quasi $1V$ case reduced from $1D$-$3V$ case)

Shock structure prediction is a benchmark test case for non-equilibrium flow models and corresponding numerical methods. In macroscopic point of view, normal shock wave is a discontinuity in space, across which physical variables change suddenly. While, in microscopic point of view (zoom into the thin shock wave), the physical variables in the shock wave changes smoothly from the front to the back of shock wave. Physically, the molecules in the shock wave are composed of the molecules before the shock (super/hyper-sonic, low temperature) and the molecules after the shock (subsonic, high temperature). When shock $Mach$ number is high, the separation of distribution functions before and after the shock in velocity space is large. Since the molecular collisions in the thin shock wave (about twenty m.f.p.) are insufficient, the distribution function will be far from equilibrium (high non-equilibrium).

Variable Soft Sphere (VSS) model is used in this case since it can be reduced to Hard Sphere (HS) model and inverse power potential model directly by using their scattering factor $\alpha$ and heat index $\omega$. The m.f.p. of VSS model is defined as

[TABLE]

where $\beta$ is defined as

[TABLE]

In shock structure case, the upstream and downstream conditions are determined by Rankine-Hugoniot relation. The iteration time step of FP-DVM is chosen as

[TABLE]

here the subscript “TP” stands for “transport” and corresponds to the free transport operator. $\Delta t_{FP}$ and $\Delta t_{TP}$ can be calculated using the following equation,

[TABLE]

4.5.1 Mach 1.2

In this case, HS molecule model ($\omega=0.5$, $\alpha=1.0$) is used which is the same with the deterministic solution of full-Boltzmann equation in Ref. [32]. For full-Boltzmann solution, x coordinate is non-dimensionalized using the m.f.p. of HS molecule [32]. The density, temperature, stress, and heat flux in the shock wave are non-dimensionalized using

[TABLE]

here subscripts “up” and “down” are used to indicate the variables in the upstream and downstream of shock wave, respectively. The truncated discrete velocity space is $[-7\sqrt{RT_{up}},7\sqrt{RT_{up}}]$ with 70 points. The cell Knudsen number $Kn_{cell}={m.f.p.}/{\Delta x}$ is chosen as 4.0 (the cell length in physical space is a quarter of m.f.p.). $CFL_{FP}$ and $CFL_{TP}$ are set 1.0 and 0.9, respectively. The density/temperature profile and stress/heat flux profile are illustrated in Fig. 8(a) and Fig. 8(b). The FP-DVM predictions match well with the full-Boltzmann result in Ref. [32]. In this case, the $Mach$ number is low, and the non-equilibrium is not too strong.

4.5.2 Mach 3.0

The same HS model ($\omega=0.5$, $\alpha=1.0$) as $Mach$ 1.2 case is used here along with the same non-dimensionalized x coordinate and physical variable (Eq. 36). Comparing with the $Mach$ $1.2$ case, the degree of non-equilibrium increases in this case. The truncated discrete velocity is $[-10\sqrt{RT_{up}},10\sqrt{RT_{up}}]$ with $100$ points. $Kn_{cell}=4.0$ is used. $CFL_{FP}$ and $CFL_{TP}$ are set 1.0 and 0.8, respectively. The density/temperature profile and stress/heat flux profile are illustrated in Fig. 9(a) and Fig. 9(b). The FP-DVM predictions match well with the full-Boltzmann result in Ref. [32], expect that the temperature profile rises a little earlier, so are the stress and heat flux profiles.

4.5.3 Mach 8.0

The working gas is Argon in this case. When $Mach$ number is 8.0, the flow inside the shock wave is in high non-equilibrium. Being the same with Ref. [33], a 11th power inverse power potential model is used, whose model coefficients can be calculated from Ref. [34] as $\omega=0.68$, $\alpha=1.4225$. The density, temperature, stress and heat flux profiles are calculated and compared with the DSMC results in Ref. [33]. Instead of using the m.f.p. of inverse power potential model and being the same with the setting in Ref. [33], the x coordinates in the profiles are non-dimensionalized using the m.f.p. of HS molecule. The density and temperature are normalized using

[TABLE]

The stress and heat flux are non-dimensionalized using

[TABLE]

The truncated discrete velocity space is $[-30\sqrt{RT_{up}},30\sqrt{RT_{up}}]$ with 300 points, and the cell Knudsen number $Kn_{cell}=4.0$. $CFL_{FP}$ and $CFL_{TP}$ are set 1.0 and 0.6, respectively. The density/temperature profile and stress/heat flux profile are illustrated in Fig. 10(a) and Fig. 10(b) along with the numerical prediction using BGK-type Shakhov model in Ref. [35]. The temperature, stress and heat flux profiles predicted by Shakhov model deviate from the DSMC results in the front of the shock wave, while the FP-DVM predictions match well with the DSMC results in such a high $Mach$ number and high non-equilibrium case.

Since the aim of the above cases is examining the validity of the present FP-DVM method, then the scope of the truncated velocity space and the amount of the discrete velocity points are set to be large. For $Mach$ 8.0 case, the velocity scope $[-25\sqrt{RT_{up}},25\sqrt{RT_{up}}]$ is sufficient. With different amount of discrete velocity points in $\xi_{1}$ direction ($300$, $200$, $100$ and $50$ points, respectively), the numerical results predicted by FP-DVM are examined in Fig. 11. It can be seen that the results with different amount of discrete velocity points almost coincide with each other, expect in the front of the shock wave where the results obtained using $50$ discrete velocity points deviate sightly from the others. The distribution functions $F$ (mass distribution along $\xi_{1}$) and $G$ (energy distribution along $\xi_{1}$) at different locations inside the shock wave are shown in Fig. 12. They are predicted using $300$ points and $50$ points in $\xi_{1}$ direction, respectively. Since this case is a high non-equilibrium one, the distribution functions inside the shock wave deviate much from the Maxwellian distribution. It can be seen from Fig. 12 that the positive property of distribution function is fulfilled in this high non-equilibrium case. It can also be seen in Fig. 12 that distributions predicted using $50$ discrete velocity points match well with that predicted using $300$ points. Since the discrete velocity space with $50$ points is very coarse, its resolution for precipitous distribution is low. For example, the setting of $50$ discrete velocity points only has $7$ points for approximating the peak of $F$ at $x=-5$ (Fig. 12(a)), and may be the reason of slight deviations in the front of the shock wave (Fig. 11).

5 Conclusion

In this paper, a deterministic FP-DVM method is proposed for the non-equilibrium flow simulations. The conservation problem of the discrete ES-FP equation is resolved by multiplying conservative coefficients whose differences with unity are small and have the same orders with the truncation error of the difference scheme for ES-FP collision operator. Using four $0D$-$3V$ cases which mimic different types of distributions that exist in real flow fields, the validity of FP-DVM method and ES-FP model for homogenous cases are proved. In these cases, the evolution of non-equilibrium anisotropic temperatures and heat flux match with the analytical Boltzmann solution precisely. To further extend the scope to ordinary inhomogeneous cases, a $nD$-quasi $nV$ reduction for $nD$-$3V$ ($n<3$) ES-FP equation is developed, which can greatly reduce the computational cost. Using the reduced $1D$-quasi $1V$ FP-DVM method, the shock structure cases from low to high $Mach$ numbers are calculated. The negative distribution function and early rise of temperature profile for high $Mach$ number cases do not appear in the present FP-DVM predictions. All the density, temperature, stress and heat flux profiles match well with the direct full-Botlzmann results and DSMC results. The validity and accuracy of both FP-DVM and ES-FP model for non-equilibrium flow simulation are proved, and the shock structure profiles predicted by the FP-DVM method are probably the best numerical prediction using model Boltzmann equations up to now to the best of our knowledge. Since the FP-type model equations have stiffness problem, the penalty and implicit treatments used in the previous researches should be consider in the further works of the FP-DVM method in order to further increase its iteration time step. The time integral solution of BGK-type equations can also be used in the flux calculation of the present method to extend its scope to whole flow regime (make it a UGKS-type method).

Acknowledgements

The authors thank Prof. Kun Xu in Hong Kong University of Science and Technology for discussions of gas kinetic theory and multi-scale flow mechanism. Sha Liu thanks Prof. Jun Zhang in Beihang University and Dr. Fei Fei in Huazhong University of Science and Technology for discussion of Fokker-Planck equation and its numerical method. This work is supported by Fundamental Research Funds for the Central Universities (No. G2018KY0302) and 111 Project of China (No. B17037).

Reference