A meshfree method for the BGK model for rarefied gas dynamics

TL;DR

This paper introduces a meshfree semi-Lagrangian method using MLS interpolation to solve the BGK model for rarefied gas dynamics, demonstrating improved accuracy over traditional spline methods through validation against DSMC and Euler solutions.

Contribution

It presents a novel meshfree semi-Lagrangian approach with MLS interpolation for BGK equations, enhancing solution accuracy in rarefied gas dynamics simulations.

Findings

Meshfree MLS interpolation outperforms spline interpolation.

The method accurately captures shock tube problems across various mean free paths.

Results align well with DSMC and Euler solutions for validation.

Abstract

In this paper we have applied a Semi-Lagrangian schemes with meshfree interpolation, based on a Moving Least Squares (MLS) method, to solve the BGK model for rarefied gas dynamics. Sod's shock tube problems are presented for a large range of mean free paths in one dimensional physical space and three dimensional velocity space. In order to validate the solutions obtained from the meshfree method, we have used the piecewise linear spline interpolation. Furthermore, we have compared the solutions of the BGK model with the solutions obtained from Direct Simulation Monte Carlo (DSMC) method. In the case of a very small mean free path the numerical solutions are compared with the exact solutions of the compressible Euler equations. Overall we found that the meshfree interpolation gives better approximation than the piecewise linear spline interpolation.

Equations

$$\frac{\partial f}{\partial t}+v_x\frac{\partial f}{\partial x}=\frac{1}{\tau }(M-f)$$

$$M=\frac{\rho }{(2\pi RT{)}^{3/2}}\exp \left(-\frac{|v-U{|}^2}{2RT}\right),$$

$$(\rho ,\rho U,E)={\int }_{{\mathbb{R}}^3}\varphi (v)f(t,x,v)dv.$$

$$e(t,x)=\frac{3}{2}RT,\rho e=E-\frac{1}{2}\rho |U{|}^2.$$

$$\tau =\frac{4\lambda }{\pi \bar{C}},$$

$$\lambda =\frac{k_B}{\sqrt{2\pi \rho R{d}^2}},$$

$$\frac{d{f}_{jkl}}{dt}$$

$$\frac{dx}{dt}$$

$$x(t_n)=\tilde{x},f_{jkl}(t_n)=f_{jkl}^n(\tilde{x})={\tilde{f}}_{jkl}^n,t\in [t_n,t_{n+1}]$$

$$M_{jkl}=\frac{\rho }{(2\pi RT{)}^{3/2}}\exp \left(-\frac{(v_j-U_x{)}^2+(v_k-U_y{)}^2+(v_l-U_z{)}^2}{2RT}\right).$$

$$f_{ijkl}^{n+1}$$

$$x_i^{n+1}$$

$${\rho }_i^{n+1}$$

$$(\rho U_x{)}^{n+1}$$

$$(\rho U_y{)}^{n+1}$$

$$(\rho U_z{)}^{n+1}$$

$$E_i^{n+1}$$

$$f_{ijkl}^{n+1}=\frac{{\tau }_i^n{\tilde{f}}_{ijkl}^n+\Delta t{M}_{ijkl}^{n+1}}{{\tau }_i^n+\Delta t}\text{mbox}fori=1,\dots ,N_x+1,j,k,l=1,\dots ,N_v+1.$$

$$M_{\Gamma }^{n+1}=\frac{{\rho }_w}{(2\pi R{T}_w{)}^{3/2}}\exp \left(-\frac{|v-U_w{|}^2}{2R{T}_w}\right),$$

$${\rho }_w=-\frac{{\int }_{(v-U_w)\cdot \nu <0}(v-U_w)\cdot \nu f_{\Gamma }^{n+1}dv}{{\int }_{(v-U_w)\cdot \nu >0}(v-U_w)\cdot \nu \frac{1}{(2\pi R{T}_w{)}^{3/2}}{e}^{-\frac{|v-U_w{|}^2}{2R{T}_w}}dv}.$$

$$f(\tilde{x})=f_k+\frac{f_{k+1}-f_k}{x_{k+1}-x_k}(\tilde{x}-x_k)=\frac{f_k(x_{k+1}-\tilde{x})+f_{k+1}(\tilde{x}-x_k)}{x_{k+1}-x_k}.$$

$$\displaystyle w_{j}=w({x}_{j}-\tilde{x};h)=\left\{\begin{array}[]{l}\exp\left(-\alpha\frac{(x_{j}-\tilde{x})^{2}}{h^{2}}\right),\quad\mbox{if }\frac{|x_{j}-\tilde{x}|}{h}\leq 1\\ 0,\qquad\qquad\mbox{else},\end{array}\right.$$

$$f(x_j)=f(\tilde{x})+(x_j-\tilde{x})\frac{\partial f}{\partial x}+e_j,$$

$$f_2-f_1$$

$$⋮$$

$$f_m-f_1$$

$$e=b-M\frac{\partial f}{\partial x},$$

$$J=\sum _{j=2}^m{w}_j{e}_j^2=(M\frac{\partial f}{\partial x}-b{)}^{T}W(M\frac{\partial f}{\partial x}-b),$$

$$\frac{\partial f}{\partial x}=(M^{T}WM{)}^{-1}(M^{T}W)b=\frac{{\sum }_{j=2}^m{w}_j(x_j-x_1)(f_j-f_1)}{{\sum }_{j=2}^m{w}_j(x_j-x_1{)}^2}.$$

$$f(\tilde{x})=f(x_1)+(\tilde{x}-x_1)\frac{\partial f}{\partial x}$$

$${\rho }^0={\rho }_l,U^0=0,e^0=2.5\text{mbox}for0\le x<0.5$$

$${\rho }^0={\rho }_r,U^0=0,e^0=2.0\text{mbox}for0.5\le x\le 1,$$

$$U(t)=0,e(t)=2.5\text{mbox}atx=0\text{mbox}andU(t)=0,e(t)=2\text{mbox}atx=1.$$

Full text

Keywords. rarefied gas, kinetic equation, BGK model, meshfree method, semi-implicit method

1 Introduction

The Boltzmann equation is an evolution equation of a probability distribution function consisting of transport and collision terms [6]. Due to the high dimensional integral in the collision term, deterministic numerical approaches are complicated and time consuming. Therefore, stochastic numerical methods like DSMC, see [4, 1, 15], have been used extensively for complex applications. DSMC methods are suitable for high Mach number and stationary flows. However, for low Mach number flows, the statistical noise inherent in these methods dominates the flow quantities. Here, our main interest is to develop a numerical method for such low Mach number, time dependent flows in arbitrary geometries. Since last 20 years micro-nano scale rarefied gas flows have attracted many researchers due to the fabricated techniques in Micro-Electro-Mechanical-Systems (MEMS) [10], devices, for examples, micro pump, micro turbines, micro pipes [13]. We consider a simplified model for rarefied gas flows, where deterministic methods can be applied more easily. We choose a simplified model suggested by Bhatnager, Gross and Krook [3], the so called BGK model for the Boltzmann equation, where the collision term is replaced by a relaxation of the distribution function towards a local thermal equilibrium. For deterministic schemes for this model we refer to [14] and other references therein. In the present paper, we apply the Semi-Lagrangian scheme suggested by Russo and Filbet, see [16] for details. In contrast to [16], where one dimensional physical and one dimensional velocity spaces are considered, we consider here a three dimensional velocity space. Moreover, the reconstruction procedure is different compared to the one applied in [16]. Here, we use a meshfree method for the reconstruction. Meshfree methods are suitable for changing computational domains in time or flows in complicated geometries, see [18, 19]. We note that a meshfree method based on Least-squares was applied to solve the compressible Euler equations, see [9] and other references there.

The paper is organized as follows. In section 2 the BGK model for the Boltzmann equation is presented. In section 3 the semi Lagrangian scheme for the model and the boundary conditions are described. In section 4 we present piecewise linear spline interpolation and the moving least squares (MLS) approximation for the reconstruction of the function. In section 5 Sod’s shock tube problem [17] is solved for several range of mean free paths. For larger Knudsen numbers the numerical solutions for the BGK model obtained from the piecewise linear spline and MLS interpolations are compared with the solutions obtained from the DSMC simulations of the Boltzmann equation. For a very small Knudsen number, we have compared the numerical solutions of the BGK model with the exact solutions of the compressible Euler equations. We found that the solutions obtained from the MLS approximation are closer to the DSMC results or the exact solutions than the solutions obtained from the piecewise linear spline interpolation. Moreover, we have compared numerical approaches based on continuous and discrete Maxwellians as suggested in [14]. We found that the use of discrete Maxwellian allows us to reduce the number of velocity grids, which is very important in higher dimensional cases from the memory as well as computational point of view. Finally, in section 6 some conclusions and future works are presented.

2 The BGK model for rarefied gas dynamics

The BGK model is the simplified model of the Boltzmann equation for a rarefied gas dynamics, where the collision term is modeled by a relaxation of the distribution function $f(t,x,v)$ to the Maxwelian equilibrium distribution. The is the evolution equation of the distribution function $f(t,x,v)$ and is given by the following initial boundary value problem

[TABLE]

with $f(0,x,v)=f_{0}(x,v),\;t\geq 0,x\in[a,b]\subset\mathbb{R},\;{v}\in\mathbb{R}^{3}$ and some boundary conditions assigned at $a$ and $b$, as will be described in the next section. We denote by $v=(v_{x},v_{y},v_{z})$ the $3d$ velocity vector. Here $\tau$ is the relaxation time and $M$ is the local Maxwellian given by

[TABLE]

where the parameters $\rho,U,T$ are macroscopic quantities, namely, density, mean velocity and temperature, respectively. Here, $R$ is the gas constant. The macroscopic quantities $\rho,U,T$ are computed from $f(t,x,v)$ as its moments. In this case we have denoted $U=(U_{x},U_{y},U_{z})$. Let ${\phi}(v)=\left(1,v,\frac{|v|^{2}}{2}\right)$ be the collision invariants. The moments are defined by

[TABLE]

Here, $E$ is the total energy density and it is related to the temperature through the internal energy

[TABLE]

The relaxation time $\tau$ and the mean free path $\lambda$ are related according to [7]

[TABLE]

where $\bar{C}=\sqrt{\frac{8RT}{\pi}}$ and the mean free path is given by

[TABLE]

where $k_{B}$ is the boltzmann constant and $d$ is the diameter of the gas.

3 Semi-Lagrangian scheme for the BGK model

To solve the BGK model, we have used the Semi Lagrangian method suggested by Russo & Filbet, see [16] for details. This method is Semi Lagrangian for the advection and implicit in the treatment of collision. In this paper we give a short description of the method. We consider constant time step $\Delta t$, uniform meshes in velocity space with mesh size $\Delta v$ and in physical space not necessarily uniform meshes with average spacing $\Delta x$. Let $t_{\rm final}$ be the final time step of computation. The time steps are given by $t_{n}=n\Delta t,n=0,1,\ldots$. The space discretization is obtained by generating grid points (regular or irregular) ${x}_{i}\in[a,b],i=1,\ldots,N_{x}+1$, where $N_{x}+1$ is the total number of grid points in physical space. We note that the $N_{x}+1$ grid points include interior as well boundary points $x_{1}=a$ and $x_{N_{x}+1}=b$. Consider the $N_{v}$ velocity grid points in each directions, where the uniform velocity grid size is given by $\Delta v=\frac{2v_{\rm max}}{N_{v}}$. The x-component of velocity grids are defined by $v_{j}=-v_{\rm max}+(j-1)\Delta v,j=1,\ldots,N_{v}+1$. Similarly, the $y$- and $z$- components are defined by $v_{k}$ and $v_{l}$ for $k,l=1,\ldots,N_{v}+1$. Assuming that $f$ is negligibly small for $|v|>v_{\rm max}$.

Let $f_{jkl}=f_{jkl}(t,x)=f(t,x,v_{j},v_{k},v_{l})$. The evolution equation of $f_{jkl}(t,x)$ along the characteristics between time steps $n$ and $n+1$ is calculated from the Lagrangian form of the discrete BGK model

[TABLE]

with initial conditions

[TABLE]

together with boundary conditions for $f_{jkl}$ at boundary points.

Here $M_{jkl}$ is still the local Maxwellian having the same moments of $f_{jkl}$ and is re-expressed by

[TABLE]

We solve Eq. (6) by the implicit Euler scheme

[TABLE]

where $\tau_{i}^{n}$ is the relaxation time in grid point $i$ at time level $n$ and the characteristic equation (7) is solved by

[TABLE]

where the initial position $\tilde{x}_{ijkl}$ is given by $\tilde{x}_{ijkl}=x_{i}-v_{j}\Delta t$. In Figure 1 we have given the geometrical interpretation. At the time level $t^{n}$ all values $f^{n}_{ijkl},i=1,\ldots,N_{x}+1$ are known. At the time level $t^{n+1}$ the corresponding values are $f_{ijkl}^{n+1},i=1,\ldots,N_{x}+1$.

The method consists of three steps:

(i) First, we determine $\tilde{x}_{ijkl}$ form the backward characteristics $\tilde{x}_{ijkl}=x_{i}^{n+1}-v_{j}\Delta t$, see Figure 1. Then, we reconstruct (or interpolate) the function $\tilde{f}^{n}_{ijkl}$ at $\tilde{x}_{ijkl}$ from the values of its neighboring grid points. One can use any reconstruction, for example, spline interpolations, least squares interpolations. In this paper we use the piecewise linear interpolation (linear spline) and the linear moving least squares interpolation. Higher order interpolations are also possible, for examples, piecewise cubic spline polynomial [16], higher order MLS. Higher order reconstructions give oscillations if the solutions develop shocks, therefore, WENO reconstruction [5] are necessary to damp the oscillations.

(ii) In the second step we obtain $M_{ijkl}^{n+1}$. Since $M_{i}^{n+1}$ and $f_{i}^{n+1}$ have same conservative moments, we multiply the above discrete equation (10) by the collisional invariants $\phi(v)$ and sum over the velocity getting

[TABLE]

Now from (12 - 16) together with (4) we obtain all five parameters of the Maxwellian and can define $M^{n+1}_{ijkl}$ from (9).

(iii) Finally, we update the density function solving Eq. (10)

[TABLE]

3.1 Boundary conditions

We use the diffuse reflection boundary conditions. This means, when gas molecules hit boundaries, we forget their history. We reflect them according to the half Maxwellian with the wall density $\rho_{w}$, wall temperature $T_{w}$ and wall velocity $U_{w}$. This is the Maxwell boundary condition, see [6] for details. Let $\nu$ be the unit normal on the wall pointing towards the computational domain. In the case of Spline interpolation we need ghost points next to the boundary points $i$ in order to compute $\tilde{f}^{n}_{ijkl}$. In order to apply the boundary condition, we first compute $f^{n+1}_{ijkl}$ in all interior points $i$, then we extrapolate the new distribution function $f^{n+1}_{\Gamma}$ for $(v-U_{w})\cdot\nu>0$ on the boundary points with the help of MLS interpolation. For $(v-U_{w})\cdot\nu<0$ the diffuse reflection boundary conditions are obtained according to

[TABLE]

where

[TABLE]

4 Interpolation methods

As we have already mentioned, our main aim is to develop a method to simulate the interactions between rigid body motion and rarefied gas. Due to the movement of a rigid body the computational domain for a gas changes. Moreover, the intersection of the surface of a rigid body and cells of rarefied gas makes the numerical scheme more complicated. In the vicinity of a moving rigid body the regular grid structure does not exists any more. In this section we present two interpolation methods, which are suitable for irregular grids. In this paper we consider linear interpolation. Higher order interpolations require some stable reconstructions, like WENO, which will be focused in future works.

4.1 Piecewise linear interpolation (or linear Spline $S_{1}$)

This is simple to implement and faster than the MLS method. However, one has to add the ghost points to apply boundary conditions, which could be complicated for complex boundaries. Let $I_{k}=[x_{k},x_{k+1}]\subset[a,b]$ be an arbitrary interval and $\tilde{x}\in[x_{k},x_{k+1}]$ be an arbitrary point. The corresponding function values are $f_{k}=f(x_{k})$ and $f_{k+1}=f(x_{k+1})$. The linear interpolantion at $\tilde{x}$ is given by

[TABLE]

We note that the size of intervals $I_{k}$ need not to be equal.

4.2 Moving least squares (MLS) interpolation

In contrast to $S_{1}$ interpolation, this is a fully meshfree method. In the $S_{1}$ interpolation, only the next left and right grid points are used to interpolate. However, in the MLS approximation, we use the nearest neighbor points inside a radius, which is about 2.5 times the average grid space. Therefore, the computational costs increases slightly in the case of MLS approximation compared to $S_{1}$ interpolation. But the MLS gives better approximation than the $S_{1}$ interpolation. In this case also the distribution of grids need not to be uniform. Another advantage of this method is that it is not required to add the ghost points to apply the boundary conditions and is easy to handle complex geometries.

Let $\tilde{x}\in[a,b]$ be an arbitrary point. We consider the problem to approximate the function $f=f(\tilde{x})$ at $x$ from the values of its neighboring points. We associate a weight function such that the near particles have higher and the far particles have lower influence. In order to limit the number of points the neighboring points are taken those points inside the circle of radius $h$ with center $\tilde{x}$. We choose the radius $h$, for example, some factor of $\Delta x$, such that we have at least minimum number of neighbors for the least squares approximation. Let $P(\tilde{x})=\{x_{j},j=1,\ldots,m\}$ be the set of $m$ neighbor points of $\tilde{x}$ inside the radius $h$. We note that this neighboring list is similar to the central stencils in the sense of the finite difference method. Therefore, if the relaxation time $\tau$ is very small and the solution of the Boltzmann equation develop shocks, we need to sort out the neighbor list according to the sign of the velocity $v_{x}$. The weight function can be quite arbitrary, but in our computations, we consider a Gaussian weight function

[TABLE]

with $\alpha$ a user defined positive constant. In our computation, we have considered $\alpha=6$. Let us sort out the neighboring points from $1$ to $m$ with respect to distance. This means, the neighbor index $1$ is the nearest neighbor of $\tilde{x}$.

In order to approximate the function we consider the $m$ Taylor’s expansions of $f(x_{j})$ around $\tilde{x}$

[TABLE]

for $j=1,\ldots,m$, where $e_{j}$ is the error in the Taylor’s expansion. We first assume that $f$ approximates the nearest point $f_{1}$. In other words, $e_{1}=0$. The unknowns $f,\frac{\partial f}{\partial x}$ are computed by minimizing the error $e_{j}$ for $j=2,\ldots,m$ and setting the constraint $e_{1}$ = 0. To solve this constraint least squares problem, we subtract the first equation with $e_{1}=0$ to all the other equations and the system of equations can be rewritten in the form

[TABLE]

The system of equations can be written in the vector form as

[TABLE]

where ${e}=[e_{2},\ldots,e_{m}]^{T}$, $b=[f_{2}-f_{1},\ldots,f_{m}-f_{1}]^{T}$ and $M=[x_{2}-x_{1},\ldots,x_{m}-x_{1}]^{T}$ . For $m>2$, this system of equations is over-determined for one unknown $\frac{\partial f}{\partial x}$ . The unknown $\frac{\partial f}{\partial x}$ is obtained from the weighted least squares method by minimizing the quadratic form

[TABLE]

where $W=w_{j}\delta_{jk},k=2,\ldots,m$ is the diagonal matrix. The minimization of $J$ formally yields

[TABLE]

Now from the equation (22) with $e_{1}=0$ for the closest point $x_{1}$ we can compute the value of $f(\tilde{x})$ as

[TABLE]

since $\frac{\partial f}{\partial x}$ is now known. We note that the higher order approximations are straightforward. Moreover, the approximation in two and three dimensional physical space is also straight forward. For example, for higher order, say order $p$, one needs to use the Taylor expansion up to $p$. Then one obtains an overdetermined system of $m-1$ equations in $p<m$ unknowns which are the derivatives of order $k,k=1,\ldots,p$. Such a system can be solved in the least squares sense with the technique shown before. In several dimensions, say dimension $d$, one has to use the Taylor expansion in dimension $d$ up to order $p$, obtaining again an overdetermined set of equations that can be solved in the least squares sense. We refer to our earlier papers [18, 19] for higher orders and several dimensions.

5 Numerical results

We have considered the Sod’s shock tube problem [17] as Benchmark to validate our numerical methods. We consider the computational domain $[a,b]=[0,1]$. The initial conditions are

[TABLE]

[TABLE]

where $\rho_{l}$ and $\rho_{r}$ denote the density on the left and right half intervals, respectively. We consider the following boundary conditions

[TABLE]

In the MLS we set $\alpha=6.0$ in (4.2) and $h$ equal to $2.5$ times the initial spacing of the grids. The initial spacing of the grids is given by $\Delta x=1/N_{x}$. We consider the Argon gas with diameter $d=0.368\times 10^{-9}m$, Boltzmann constant $k_{B}=1.3806\times 10^{-23}J~{}K^{-1}$ and the gas constant $R=208J{\rm Kg}^{-1}~{}K^{-1}$. The corresponding initial temperature are $0.008012~{}K$ on the left half of the domain and $0.00641~{}K$ on the right half of domain. The limit of the velocities in all direction is set by $v_{\rm max}=10m/s$. We note that the thermal velocity corresponding to such a temperature is much smaller than the chosen $v_{\rm max}$. So, initially the gas is distributed according to the Maxwellian with these initial parameters. The final time is $t_{\rm final}=0.17$ seconds.

5.1 Test 1: Comparison of solutions with different CFL numbers

One of the main advantages of this scheme is that the Courant number CFL $=v_{\rm max}\Delta t/\Delta x$ can be bigger than 1. In the first example we consider $N_{x}=200$, $\rho_{l}=10^{-4}\,{\rm Kg}\,m^{-1},\rho_{r}=0.125\times\rho_{l}$, which corresponds to the initial mean free path ${\lambda}_{l}=0.001$ m for $x<0.5$ and ${\lambda}_{r}=0.008$ m for $x>0.5$ and the corresponding initial relaxation time can be computed from the Eq. (5). The initial condition is the Maxwellian distribution (2) with parameters given by the initial density, mean velocity and the temperature. We have plotted the temperature for CFL = 1 and 2, where the solutions are the same for both cases, see Figure 2.

For the optimal choice of the CFL number we refer [11]. In the following test cases, we have considered CFL = 1.

5.2 Test 2: Comparison of solutions with constant and variable $\tau$

In most of the DSMC simulations constant mean free paths are considered. The considerations of variable mean free paths may effect too much in the solutions since there is large fluctuation in the density in the DSMC simulations. Since we validate our numerical scheme for the BGK equation with the full Boltzmann equation with the help of the DSMC simulations, we fix the Knudsen number initially and keep it constant until the final time step. In shock tube problems there is jump in the initial mean free path. Therefore, we consider the average value of the mean free path in all time steps. Based on the average mean free path we define the initial relaxation time $\tau_{l}$ and $\tau_{r}$ with the help of Eq. (5) and use these values, for example, $\tau_{l}$ in the domain $0\leq x\leq 0.5$ and $\tau_{r}$ in the rest of the domain throughout the simulations. In [16] the authors have used the constant relaxation time to solve the BGK equation. However, one can use the variable $\tau$ in the time and space, see Eq. (17). For very small $\lambda$ or $\tau$ the solutions obtained from variable and constant $\tau$ do not differ much. In Figure 3 we have plotted the temperature for constant and variable $\tau$. We have considered the density ratios $\rho_{l}/\rho_{r}=8$. Three densities are considered, which are $\rho_{l}=5\times 10^{-6}{\rm Kg}\,m^{-3},~{}10^{-4}{\rm Kg}\,m^{-3},~{}1{\rm Kg}\,m^{-3}$ which correspond to the mean free paths $\lambda_{l}=0.02,~{}0.001,~{}10^{-7}$ meters, respectively. We see that the difference in the solution obtained from the constant and variable $\tau$ for larger mean free paths. When the mean free path is very small there is no difference between the solutions.

5.3 Test 3: Comparison of solutions in regular vs irregular grids

In this test case we present the comparison of the numerical solutions in regular as well as irregular grids. We have considered $N_{x}=200$. The regular grids are generated according to $x_{i}=(i-1)*\Delta x,i=1,\ldots,N_{x}+1$. To create the irregular grids we have moved the regular grids with velocity $\Delta x/4$ times the random number $i=2,\ldots,N_{x}$. This movement is performed for 2 iterations. The densities are $\rho_{l}=5\times 10^{-6}{\rm Kg}~{}m^{-3}$ and $\rho_{r}=0.125\times 10^{-6}{\rm Kg}~{}m^{-3}$. The corresponding initial mean free paths are $\lambda_{l}=0.02$ m on the left and $\lambda_{r}=0.17$ m on the right half of the domain. The corresponding initial relaxation times are $0.01$ on the left half of the domain and $0.0957$ on the right half of the domain. The flow is in transition regime. In Figure 4 we have plotted the zoom of the regular and irregular grids on the left and the densities obtained from the MLS interpolation in regular as well as irregular grids. We observe that the irregular grids make no difference to the solutions obtained from the regular grids.

5.4 Test 4

In this case, we have compared the numerical solutions obtained from the BGK model with the DSMC simulations [1, 15] for the Boltzmann equation since DSMC results are widely used as Benchmark solutions. All parameters are same as in the Test 2. For DSMC simulations we have considered $200$ cells and $400$ gas molecules are initially generated per cell according to the Maxwellian distribution in the velocity, where the initial density, temperature and velocity are its parameters. Notice that DSMC is a method to approximate the solution to the Boltzmann equation for hard spheres, while BGK is a simplified model of the Boltzmann equation, therefore we do not expect to observe the same behaviour. Both models have in common that in the limit of very small mean free path converge to the compressible Euler equations for a monoatomic gas, therefore we expect that the two models provide similar results when adopted with the same small mean free path. Furthermore, we observe that standard BGK fails to correctly capture first order effects in the (small) mean free path, because it is not possible to match at the same time thermal conductivity and viscosity coefficient with the single parameter $\tau$. There are extensions of the BGK model, such as, for example, the so called ES-BGK. It is possible to mention the existence of other BGK models, such as for example the ellipsoidal BGK [12], that are able to capture the correct Navier-Stokes limit. The use of such models is however beyond the scope of the present paper.

The cell size for DSMC simulations must be smaller than the mean free path. Therefore, the time step of the DSMC solver is restricted by the mean free path, see [1, 15] for details. However, the time step for the Semi Lagrangian scheme is not restricted, which has shown in Test 1. Since the flow has low Mach number, the statistical fluctuations dominates the DSMC solutions. Therefore, we have obtained $20$ independent runs. The initial density is chosen same as in the Test 3. In Figure 5 we have plotted the density, velocity and pressure obtained from the DSMC simulations and the BGK model using linear MLS and linear spline interpolations. We observe that the linear MLS interpolation scheme for the BGK model and DSMC results match perfectly, however, the linear spline interpolation gives some deviations from the solutions of the DSMC simulations.

5.5 Test 5

In this test case we have increased the density by factor $2$ and $20$ times compared to the Test case 3 and 4 such that $\rho_{l}=10^{-5}{\rm Kg}m^{-3}$ and $10^{-4}{\rm Kg}m^{-3}$, respectively. The corresponding mean free paths are $0.01$ and $0.001$ meters on the left half of the domain and $8$ times larger on the right half. For $\rho_{l}=10^{-5}{\rm Kg}m^{-3}$ we have considered the $400$ cells in the DSMC simulations and the same number of grids for the BGK model. Similarly, for $\rho_{l}=10^{-4}{\rm Kg}m^{-3}$ we have considered the $1000$ cells in the DSMC simulations and the $800$ grids for the BGK model. The increase of the number of cells is due to the restriction that the DSMC cells must be smaller than the mean free path. Other parameters are same as in the earlier test cases.

For the DSMC simulations we have performed $20$ independent runs. In Figures 6 and 7 we have plotted the numerical solutions from all three methods. We again observe that the DSMC solutions and the solutions of the BGK model obtained by linear MLS interpolation are closer than the linear spline interpolation.

5.6 Test 6

In this test case we consider $\rho_{l}=1{\rm Kg}m^{-3}$ and $\rho_{r}=0.125{\rm Kg}m^{-3}$. The corresponding left and right mean free paths are $10^{-7}$ m and $8\times 10^{-7}$ m, respectively. We note that the size of DSMC cells must be smaller than $\lambda$, so, we need at least $9\times 10^{6}$ cells for such mean free path. The time step also has to be reduced accordingly. Moreover, the number of gas molecules is also very high and the computational time for the DSMC simulations becomes enormously high. In one dimensional case, this it quite a large number of cells. Therefore, we have not performed DSMC simulations in this case. However, there is no restriction of cell size for the Semi Lagrangian scheme for the BGK model. We have again used $800$ grids for the Semi Lagrangian scheme for this smaller mean free paths.

On the other hand for this small mean free path we can solve the continuum equations, for example, the compressible Euler equations. For the shock tube problem, the compressible Euler equations can be solved exactly. In Figure 8 we have again plotted the density, velocity and pressure at final time $0.17$ seconds obtained from the MLS and Spline interpolations together with the exact solutions of the compressible Euler equations. In the shock region all three solutions match perfectly, however, in the contact discontinuity and the rarefaction region, the linear MLS interpolation scheme gives better approximation than the linear spline interpolation scheme.

5.7 Test 7: Comparing solutions from continuous and discrete Maxwellian

In all above examples we have considered the velocity grids $N_{v}=20$ in all three directions the standard Maxwellian given by (2). Increasing the number of velocity grids do not help much in the accuracy of the solutions. However, smaller values of $N_{v}$ affects the accuracy. This means we loose the conservative properties of the scheme. In this test case we have considered all parameters as in the Test 6. We have observed that $N_{v}=20$ gives solutions close as the exact solutions of the compressible Euler equations. So, the solutions obtained from $N_{v}=20$ are our reference solutions. We have decreased the values of $N_{v}$ and compare the solutions with the reference solutions. The smallest one which gives the stable solutions is $N_{v}=13$ for the case of the standard Maxwellian. But the solutions deviate from the reference solutions. For smaller values of $N_{v}$ we loose the conservative properties. In other words, the moments obtained from the discrete summation in Eq. (12 - 16) are not exactly equal to the moments computed from the standard Maxwellian (2). To obtain conservative properties, one uses the discrete Maxwellian suggested by Mieussiens [14]. The discrete Maxwellian depends on five parameters, which are also related to the moments and the parameters are determined by solving the non-linear system of equations. Since the Jacobian matrix has very bad condition number, the standard Newton’s method for solving non-linear system does not work, one has to use the back tracking line search algorithm [8]. A careful comparison between the use of continuous and discrete Maxwellian in Semi Lagrangian schemes for the computation of shock problems is reported in [2].

In the case of three dimensional velocity space, the main drawback of the method is the memory problem as well as long computational time. Therefore, it is important to reduce the computational time as well as memory allocation. Thus we have considered the standard Maxwellian as well as discrete Maxwellian. If the number of velocity grid points is equal to $20$ or above, we do not see much difference of the solutions obtained from the standard as well as the discrete Maxwellian. However, the discrete Maxwellian requires more computational efforts than the standard Maxwellain because of the iterative method for solving nonlinear system of equations. If we choose the proper initial guess, the Newton iteration converges very fast. We observed that if we consider discrete Maxwellian with $13$ velocity grids, the solutions are the same as the ones obtained from the standard Maxwellian with $20$ velocity grids. But with $13$ velocity grids with standard Maxwellian, the solutions deviate from the reference solutions. In Figure 9 we have plotted the density, velocity and pressure obtained from the BGK model with discrete Maxwellian considering $N_{v}=13$ and with the standard Maxwellian considering $N_{v}=13$ together with the reference solutions. We see that the solutions obtained from the standard Maxwellian with $N_{v}=13$ deviates from the reference solutions, while the solutions with the discrete Maxwellain with $N_{v}=13$ are closer with the reference solutions (standard Maxwellian with $N_{v}=20$). The computation is performed in $dualIntelXeonGold6132("Skylake")@2.6GHz$ with intel fortran compiler. The total CPU to time with the standard Maxwellian having $20$ velocity grids is $443$ seconds while with the discrete Maxwellian with $13$ velocity grids is $166$ seconds. In higher dimensional cases the reduction of velocity grid points $N_{v}$ is very important from the computation as well as memory allocation point of view.

6 Conclusion and Outlook

In this paper we have focussed on a meshfree method in the Semi Lagrangian scheme for the BGK model for rarefied gas flows. The meshfree method is applied for the reconstruction steps as well as for the implementation of boundary conditions. The diffuse boundary conditions on the solid wall is applied. The meshfree approximation is based on the moving least squares (MLS) method. We have used a linear approximation. The advantage of the meshfree approximation is that we do not require regular distribution of grid points in the velocity space, which will be very important if the boundary moves in time or interface between gases and other medium changes in time. Another advantage of the meshfree method is that we do not need to add ghost points. We have presented also the linear piecewise interpolations in order to compare the results from the linear MLS interpolations. We observed that linear MLS gives better results than the linear spline interpolation. We have considered a problem in one dimensional physical space and three dimensional velocity space. Sod’s shock tube problem is solved for several ranges of mean free paths. No difference is found in the results obtained from regular and irregular grids in physical space. For larger mean free paths the solutions of BGK model are compared with the solutions obtained from the DSMC method for the Boltzmann equation. The solutions obtained from linear MLS are closer to DSMC results than the linear spline interpolation. For a very small mean free path we have compared the solutions of BGK model with the exact solutions of the compressible Euler equations and the solutions have very good agreement. In this case also the linear MLS approximation gives better results than the linear spline interpolation. Moreover, we have shown that there is no restriction of the CFL number, where the CFL number can be larger than 1. Furthermore, we have studied the difference between the solutions obtained from the continuous and the discrete Maxwellian distribution. We found that the use of discrete Maxwellian allows us to reduce the number of grid points in velocity space without losing the accuracy in the solutions. The reduction of number of grid points is very important in higher dimensional physical spaces from the point of view of computational time and memory allocations.

Future works will be the extension of the method for higher order reconstruction using WENO. Moreover, we are planning the extension of the method in higher dimensional physical spaces as well as the interaction of moving nano rigid particles immersed in a rarefied gas.

Acknowledgment

This work is supported by the DFG (German research foundation) under Grant No. KL 1105/30-1 and by the ITN-ETN Marie-Curie Horizon 2020 program ModCompShock, Modeling and computation of shocks and interfaces, Project ID: 642768.