 # Numerical boundary conditions in Finite Volume and Discontinuous Galerkin schemes for the simulation of rarefied flows along solid boundaries

We present a numerical comparison between two standard finite volume schemes and a discontinuous Galerkin method applied to the BGK equation of rarefied gas dynamics. We pay a particular attention to the numerical boundary conditions in order to preserve the rate of convergence of the method. Most of our analysis relies on a 1D problem (Couette flow), but we also present some results for a 2D aerodynamical flow.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The kinetic theory of Rarefied gas flows describes the behavior of a gas in a system which length is of same order of magnitude as the mean free path of the gas molecules. Numerical simulations of these flows are fundamental tools, for instance they are used in aerodynamics to estimate heat fluxes at the surface of a re-entry space vehicle at high altitudes or to estimate the attenuation of a micro-accelerometer by the surrounding gas in a micro-electro-mechanical system. One last example is when one wants to estimate the pumping speed or compression rate of a turbo-molecular pump.

For such problems, it is important to take into account the interactions of the gas with the solid boundaries of the system. They are modeled by boundary conditions like the diffuse or specular reflections. In this work, our first aim is to improve the accuracy of our in-house kinetic solver, which is based on a finite volume approach (FV), in particular for the computation of parietal fluxes. Finite volume schemes are often used in rarefied gas dynamics (RGD), see for instance the work of [20, 15, 1, 9, 22]. Our scheme uses a standard extension to second order accuracy by using the non linear flux limiter approach of Yee . We observe that the accuracy of this scheme decreases to first order at solid boundaries, because the boundary conditions are not discretized properly. The first goal of this paper is to analyze this problem, both analytically and numerically. We show that our flux limiter can hardly be compatible with second order discretization of reflection boundary conditions.

Then we study another kind of second order extension based on linear reconstruction and nonlinear slope limiters . In this case, we are able to propose a discretization of the boundary conditions (inspired by [11, 8]), based on extrapolations that are consistent with the slope limiter technique, which preserves second order accuracy up to solid boundaries. On numerical tests, the improvement of the accuracy of this method is spectacular, for 1D and 2D problems.

Finally, we want to compare this modified finite volume scheme to the popular Discontinuous Galerkin method (DG). Incidentally, this method has originally been made for neutron transport (hence a kinetic equation, like the Boltzmann equation of RGD) in the 1970’s by Reed and Hill . Such schemes have many interesting properties, and they became more and more popular in the past decades in many different fields. In particular, the DG method has been used for charged particles , and more recently for the BGK equation . According to this last paper, the DG method would be more efficient that the FV method, at least for the test cases used by the authors. In this paper, we show that the discretization of reflection conditions with a DG method is very simple. We compare a simple second order DG scheme to our two FV schemes, mainly to investigate their accuracy at solid boundaries. This study shows the order of accuracy of both DG and FV method (with slope limiters) are not decreased at solid boundaries.

The outline of this article is the following. In section 2, we present the basic BGK equation of RGD. In 3, we analyze the two versions of our FV scheme and a simple DG scheme. These schemes are compared for several test cases in section 4.

## 2 The BGK equation

In kinetic theory, a monoatomic gas is described by a particle distribution function defined such as is the mass of molecules that at time are located in an elementary space volume centered in and have a velocity in an elementary volume centered in . The macroscopic quantities associated to such as mass density , momentum and total energy

are defined as the five first moments of

with respect to the velocity variable, namely:

 (ρ(t,x),ρu(t,x),E(t,x))=∫R3(1,v,12|v|2)f(t,x,v)dv. (1)

The temperature of the gas and its pressure are defined by the relation

 E=12ρ|u|2+32ρRT (2)

where is the gas constant defined as the ratio between the Boltzmann constant and the molecular mass of the gas. When the gas is in a thermodynamical equilibrium state, it is well known that the distribution function is a Gaussian function of , called Maxwellian distribution, that depends only on macroscopic quantities and satisfies relations (1):

 M[ρ,u,T]=ρ(2πRT)32exp(−|v−u|22RT). (3)

The evolution of a gas at thermodynamical nonequilibrium is described by the following Boltzmann equation:

 ∂tf+v⋅∇xf=Q(f), (4)

which means that the total variation of (described by the left-hand side) is due to collisions between molecules ( is the collision kernel). The most realistic collision kernel is the Boltzmann operator but its use is still very computationally expensive. Here we choose the simpler BGK model [2, 21]

 Q(f)=1τ(M[ρ,u,T]−f) (5)

which makes relax towards the local equilibrium corresponding to the macroscopic quantities defined by (1). The relaxation time defined as ( being the viscosity of the gas) is chosen to recover the correct viscosity in the Chapman-Enskog expansion for Navier-Stokes equations.

The interactions of the gas with solid boundaries are described with a diffuse reflection model. If the boundary has a temperature and a velocity so that , where is the normal to the wall at point

directed into the gas, a molecule that collides with the boundary is re-emitted with a random velocity normally distributed around

(

being the variance of the normal distribution). This reads

 f(t,x,v)|Γ=σw(f)M[1,uw,Tw](v)  if v⋅n(x)>0. (6)

The parameter is the mass density of reflected particles, and it is defined so that there is no normal mass flux across the boundary (all the molecules are re-emitted), see . Namely, that is,

 σw(f)=−(∫v⋅n(x)<0v⋅n(x)f(t,x,v)dv)(∫v⋅n(x)>0v⋅n(x)M[1,uw,Tw](v)dv)−1. (7)

There are other reflection models, like the Maxwell model with partial accommodation, but they are not used in this work.

The numerical resolution of the BGK equation is done by a deterministic method using a discrete velocity approach (see [15, 16, 17], and the extension to polyatomic gases ). The velocity variable is replaced by discrete values of a Cartesian grid (see  for the use of locally refined velocity grids). The continuous distribution is then replaced by its approximation at each point , and we get the following discrete velocity BGK equation

 ∂fk∂t+vk⋅∇xfk=1τ(Mk[ρ,u,T]−fk), (8)

where . The discrete function is an approximation of the equilibrium functions defined in (3) at point . This discrete equilibrium is defined in order to satisfy the conservation of mass, momentum and energy (see ). The macroscopic quantities are now defined by the quadrature rule

 (ρ(t,x),ρu(t,x),E(t,x))=∑k(1,vk,12|vk|2)fk(t,x)ωk, (9)

where are the quadrature weights: for a Cartesian velocity grid of steps , , , the weights are given by . Moreover, the boundary condition (6) is replaced by

 fk(t,x)|Γ=σw(f)Mk[1,uw,Tw]  if vk⋅n(x)>0, (10)

where is defined as in (7) in which the integrals are replaced by quadratures:

 σw(f)=−⎛⎝∑vk⋅n(x)<0vk⋅n(x)fk(t,x)ωk⎞⎠⎛⎝∑vk⋅n(x)>0vk⋅n(x)Mk[1,uw,Tw]ωk⎞⎠−1. (11)

Again this ensures a zero mass flux across the solid wall.

## 3 Numerical schemes and boundary conditions

For the sake of simplicity, the numerical scheme will be presented in the 1D configuration, in the domain . Equation (8) is then written

 ∂fk∂t+vk∂fk∂x=1τ(Mk[ρ,u,T]−fk), (12)

where is the first component of . Since (at ), the boundary condition (10)–(11) now reads

 fk(t,0)=σw(f)Mk[1,uw,Tw]  if vk>0, (13)

where

 σw(f)=−⎛⎝∑vk<0vkfk(t,0)ωk⎞⎠⎛⎝∑vk>0vkMk[1,uw,Tw]ωk⎞⎠−1. (14)

For the following, it is important to note that boundary condition (13)–(14) implies a zero mass flux at the boundary. Indeed, (13) yields:

 ϕ:=∑kvkfk(t,0)ωk=∑vk<0vkfk(t,0)ωk+∑vk<0vkσw(f)Mk[1,uw,Tw]ωk=0

from (14). In particular, this property ensures a global mass conservation for any internal flow.

### 3.1 Finite Volume schemes

Finite Volume schemes approximate the solution of a given problem by integrating the equation on each cell of a mesh. The integration of the advection term results in a flux at cell interfaces: the numerical flux. The accuracy of the scheme depends on the accuracy of this flux. Moreover, this flux is closely linked with the boundary conditions. This is why, in this section, we study in detail three different finite volume schemes and their properties close to a solid boundary. Note that the collision term is local in space and its discretization is consequently the same for all Finite Volume schemes: therefore, even if we take the BGK equation as an example, the study performed in this section can be applied to every collision operator.

First, we assume we have a mesh of nodes with steps , for to . The discrete time variable is with a time step . We first integrate equation (12) in a cell between and and divide by to obtain

 fn+1i,k−fni,kΔtn+1Δxi(Fni+12,k−Fni−12,k)=Qk(fni), (15)

where is an approximation of the average of in the cell . The numerical flux is an approximation of the integral , and can be interpreted as a flux across the cell interface between cells and .

In the literature, there are many different constructions of this numerical flux, and each one leads to a different scheme. The simplest one is the upwind scheme (see section 3.1.1 below), which gives a first order scheme. There are several ways to increase the order of a this scheme, all of them related to a modification of the numerical flux. One solution is to reconstruct the solution inside the cells of the mesh, in order to define new values on the interfaces. Another one is to limit a high order numerical flux at the cell interfaces. Here, we study one method of each category, and in particular their properties at the solid boundary. We first present the standard first order upwind scheme. Then we study a second order flux limiter method and show it decreases to first order at the boundary. Then we present the linear reconstruction method, which gives a correct second-order approximation of the BGK problem up to the boundary.

Finally, note that, for simplicity, all our schemes are presented with a standard forward (explicit) Euler time discretization. In practice, since we are interested in steady flows only, our numerical tests are made with a linearized backward (implicit) Euler method, which is standard in aerodynamics (see [15, 1]). However, the analysis made below does not depend on the time discretization.

#### 3.1.1 A first order scheme finite volume scheme

The simplest numerical flux is the upwind flux:

 Fni+12,k=v+kfni,k+v−kfni+1,k,

where denotes the positive and negative parts of . It is well known that the corresponding scheme is first order accurate in space, with a strong numerical diffusion on coarse meshes.

Assume that the first cell of the mesh is adjacent to a solid wall. Then the numerical flux approximates the flux across the solid boundary. It requires an artificial value , that is generally interpreted as the value of the distribution function in a ghost cell inside the solid wall, adjacent to the first cell. This value must be defined so as to (a) account for the boundary condition (13)–(14), and (b) satisfy a zero mass flux across the wall, that is to say .

These constraints can be satisfied as follows. First, the value of is defined for outgoing velocities by using a zeroth order extrapolation of the value in , namely

 fn0,k=fn1,k.

Then, the value of for incoming velocities is defined by using the boundary condition (13-14), that is to say:

 fn0,k=σw(fn0)Mk[1,uw,Tw],% where σw(fn0)=−(∑lv−lfn0,lωl)/(∑lv+lMl[1,uw,Tw]ωl).

It can easily seen that this definition ensures the mass conservation. Indeed, the mass flux across the wall is

 ∑kFn1/2,kωk=∑k(v+kfn0,k+v−kfn1,k)ωk=∑k(−v+k∑lv−lfn0,lωl∑lv+lMl[1,uw,Tw]ωlMk[1,uw,Tw])ωk+∑k(v−kfn1,k)ωk=−∑l(v−lfn0,l)ωl+∑k(v−kfn1,k)ωk=0,

since for .

#### 3.1.2 Second order finite volume scheme with flux limiters

In this kind of scheme, a nonlinear limitation is applied to a centered numerical flux (which is second order accurate) to ensure a stability property. Here, we use the Yee limiter, which ensures that the scheme is “Total Variation Diminishing” (TVD) (see ):

 (16)

where the minmod function is defined by

 minmod(x,y,z)={sgn(x) min(|x|,|y|,|z|)  if sgn(x)=sgn(y)=sgn(z),0 in other cases.

Like in the first order scheme, some ghost cell values have to be defined to compute the numerical flux at the wall interface . Indeed, the first order part of the flux requires the value of for incoming velocities (), while the limiter needs and for every velocities. Again, these ghost cell values have to be defined so as to account for the boundary condition (13)–(14) to satisfy a zero mass flux across the wall.

A simple idea that satisfies these constraints is the following. For outgoing velocities, we define and by a zeroth order extrapolation from the first boundary cell:

 fn−1,k=fn0,k=fn1,k.

For incoming velocities, we define like in the first order scheme, that is to say:

 fn0,k=σw(fn0,k)Mk, where σw(fn0,k)=−(∑lv−lfn0,lωl)(∑lv+lMlωl)−1.

And again, is defined by zeroth order extrapolation with .

This definition implies that the numerical flux at the wall interface reduces to the first order flux , since by construction, and hence the limiter is zero. Consequently, the mass conservation is satisfied (see section 3.1.1). At the same time, the drawback of this approach is that the numerical flux at the wall is only first order accurate, while it is second order inside the computational domain. This lower accuracy at the wall can be clearly observed in our numerical tests (see section 4).

To increase the accuracy of the scheme, a natural idea is to use a first order extrapolation to define and for outgoing velocities, and then to use the boundary conditions to define these values for incoming velocities. This might give a second order flux, but we have not been able to find a definition for which the conservation property still holds. In our opinion, it is unlikely that such a definition exist (see an explanation at the end of section 3.1.3).

Consequently, we believe that such a scheme cannot be at the same time second order up to the boundary and conservative. This is why we propose to consider another kind of second order scheme in the following section.

#### 3.1.3 Second order finite volume second order scheme based on linear reconstruction

In the finite volume schemes using a linear reconstruction, the numerical fluxes appearing in relation are still defined as upwind fluxes, but with new values at the cells interfaces, obtained with a linear reconstruction of the distribution inside each cell :

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Fni+12,k=v+kfn,−i+12,k+v−kfn,+i+12,k,fn,−i+12,k=fni,k+Δxi2δni,k,fn,+i+12,k=fni+1,k−Δxi+12δni+1,k, (17)

where is the slope of the linear reconstruction of in cell . The most classical way to compute this slope is a least square method, which gives, for a regular mesh, the centered slope

 δni,k=fni+1,k−fni−1,k2Δxi. (18)

To obtain a TVD scheme, the slope must be limited to avoid the creation of new extrema. Here, we use the MC slope limiter defined on each cell by

 δn,limi,k=minmod(δni,k,2αifni,k−fni−1,kΔxi,2αifni+1,k−fni,kΔxi), (19)

where is a free parameter between 0 and 1 .

When this scheme is applied to the first cell , the numerical flux is

 Fn12,k=v+kfn,−12,k+v−kfn,+12,k.

It uses two values of the distribution at the solid interface: is the value on the right side of the wall, required for outgoing velocities , given by the linear reconstruction

 fn,+12,k=fn1,k−Δx12δn,lim1,k, (20)

while is the value on the left side of the wall, required for incoming velocities , given by the linear reconstruction

 fn,−12,k=fn0,k+Δx02δn,lim0,k, (21)

where the slopes are defined by (18)–(19) with and , respectively. This requires to define two ghost cell values and , and the corresponding cell size . Again, they have to be defined so as to account for the boundary condition and to preserve the conservation property.

First, note that for the conservation property holds, it is sufficient that the value of at the wall interface satisfies

 fn,−12,k=σw(fn,+12,k)Mk. (22)

Consequently, we have to construct the ghost cell values so that relations (20)–(22) hold. This can be done as follows.

The method, inspired by  [11, 8], consists in first defining the wall interface values by extrapolation (for outgoing velocities) and by the boundary condition (for incoming velocities), and then in defining the ghost cell values and so that  (20)–(21) hold. This method is detailed below in three steps, and summarized in figure 1. Note that in the following, we will assume that the ghost cells have the same size as the first cell: .

Step 1: computation of for outgoing velocities .
This value is defined by a linear extrapolation using the values and :

 fn,+1/2,k=32fn1,k−12fn2,k. (23)

Step 2: computation of for incoming velocities .
We use the boundary condition (13)–(14) to get

 fn,−1/2,k=σw(fn,+1/2)Mnk. (24)

Step 3: definition of the values and in the ghost cells.

For outgoing velocities, only the value of is required. We use the same linear extrapolation as used to compute extrapolation (see Step 1) to get

 fn0,k=2fn1,k−fn2,k.

For incoming velocities, we need the values of and . Then we use again a linear extrapolation, but based on the values of the incoming value of the wall interface distribution and on . This gives

 fn0,k =2fn,−1/2,k−fn1,k, fn−1,k =4fn,−1/2,k−3fn1,k.

Now it is not difficult to prove that these definitions satisfy the previous constraints, provided that . First, (22) is imposed at step 2, and hence it is satisfied, which gives mass conservation. Now, it remains to prove that defined by (23) and (24) also satisfy (20) and (21). This is due to the extrapolation procedures that make the points are on the same straight line, and hence the limiters can be computed. Indeed, for instance, we have for incoming velocities :

Consequently, the right-hand side of (21) is

 fn0,k+Δx02δn,lim0,k=2fn,−1/2,k−fn1,k+fn1,k−fn,−1/2,k=fn,−1/2,k

and (21) is satisfied. For outgoing velocities, we have

 δn,lim1,k=minmod(fn2,k−fn0,k2Δx1,2α1fn1,k−fn0,kΔx1,2α1fn2,k−fn1,kΔx1)=minmod(fn2,k−fn1,kΔx1,2α1fn2,k−fn1,kΔx1,2α1fn2,k−fn1,kΔx1)=fn2,k−fn1,kΔx1.

Consequently, the right-hand side of (20) is

 fn1,k−Δx12δn,lim1,k=32fn1,k−12fn2,k=fn,+1/2,k

from (23), and hence (20) is satisfied.

In summary, we have proved that this finite volume scheme is based on slope limiters whose ghost cell values are consistent with the slope reconstruction up to the solid boundary: this implies the second order accuracy up to the boundary. Moreover, we have proved that this scheme also preserves the mass conservation even with solid boundaries.

Finally, note that even if the extrapolation procedure we use at the solid wall might induce non positive values of the distribution functions, we did not observe this problem in our numerical tests. This is probably due to the fact that we use small cells in the Knudsen layer, which make the gradients, and hence the slopes, small enough too avoid the creation of negative values.

###### Remark 3.1.

The fact that scheme with a flux limiter presented in section 3.1.2 cannot be both second order and conservative at the solid wall can be seen as follows.

The previous analysis for the scheme with a slope limiter relies on the fact that the minmod function (19) reduces to a single slope. Consequently, the scheme is linear at the solid wall, and hence is compatible with boundary condition (22) (which is linear too). This implies the conservation property. This reduction is due to the fact that the ghost cell values and are defined by extrapolations that make all the points used in (19) aligned (and hence with three equal slopes).

For the scheme with a flux limiter, this extrapolation cannot make all the points aligned. Indeed, the limiter (16) uses the same stencil for negative and positive velocities (it is a “symmetric” limiter (see ). This implies that for positive velocities, and are defined through extrapolation of and , but they cannot be aligned with , in general (see figure 2). Our previous analysis cannot be applied, and there is no reason for the conservation property holds true.

However, if instead a non symmetric (or “upwind”) flux limiter is used (see ), it is still possible to prove the conservation property. Such is scheme is not used in this paper.

### 3.2 Discontinuous Galerkin scheme

In finite volume methods, the distribution is approximated by a piecewise constant function (on each cell, the discrete distribution is equal to its cell average). Discontinuous Galerkin schemes can be viewed as an extension of the finite volume method in which the distribution is now approximated by piecewise polynomial functions. Here, we will consider a piecewise linear approximation.

#### 3.2.1 Weak form and piecewise linear approximation

On each cell , we use two basis affine functions and such that , , and , . We project equation (8) in this basis and integrate by parts to get for every and :

 ∫Ωi∂fk∂t(t,x)φi,p(x)dx−∫Ωivkfk(t,x)∂xφi,p(x)dx+(vkfk(t,xi+12)φi,p(xi+12)−vkfk(t,xi−12)φi,p(xi−12))=∫Ωi1τi(Mk[ρ,u,T]−fk)φi,p(x)dx, (25)

Now, we assume that can be approximated on each cell by a piecewise linear function defined by (see figure 3). Note that while the components and are also the pointwise values of at the edges of , the whole function itself is not continuous across these edges: indeed, at each cell edge , has two left and right values and (see figure 3).

This approximation has now to be injected in the weak form (25). The first term of the left-hand side is easily computed

 2∑q=1∫Ωi∂fi,k,q∂tφi,q(x)φi,p(x)dx=2∑q=1mipq∂fi,k,q∂t, (26)

where the form the matrix

 (mipq)=¯¯Mi=|Ωi|6(2112). (27)

The second term of the left-hand side is

 −vk∫Ωifi,k(x)∂φi,p(x)∂xdx=2∑q=1dpqfi,k,q, (28)

where the form the matrix

 (dpq)=¯¯Dk=vk2(  1  1−1−1). (29)

For the third term of the left-hand side of (25), we muse take into account that the piecewise linear approximation of is not continuous across cell edges (see the remark above). Consequently, we write this term as , where a single value for has to be defined. A good choice (for accuracy, stability, and treatment of boundary conditions) is to use the upwind value

 ^fk(t,xi+12)={fi,k,2(t)if vk>0,fi+1,k,1(t)if vk<0.

Taking into account that the basis functions take values or at the cell edges of , we get

 vk^fk(t,xi+12)φi,p(xi+12)−vk^fk(t,xi−12)φi,p(xi−12)={−v+kfi−1,k,2(t)−v−kfi,k,1(t)for p=1,v+kfi,k,2(t)+v−kfi+1,k,1(t)for p=2. (30)

The right-hand side of (25) is obtained by using the same piecewise linear approximation for each term, which gives

 ∫Ci1τ(Mk[^U]−^fk)φi,p(x)dx=2∑q=1mipq1τi,q(Mi,k,q−fi,k,q), (31)

where has been defined in (27), , , and the macroscopic quantities are defined by applying (9) in each cell.

Finally, we collect the different terms (2631), and we find that the approximation of (25) can be written in the following vectorial form

 ¯¯Mi ∂fi,k∂t+¯¯Dk fi,k+(¯¯Ak fi−1,k+¯¯Bk fi,k+¯¯Ck fi+1,k)=¯¯Mi 1τi(Mi,k−fi,k), (32)

where we use the following notations for the two-component vectors

and , the 22 matrices , , , while the matrices and have been defined in (27) and (29). Note that the product of vectors must be understood component-wise. Multiplying (32) by the matrix (the inverse of ), we obtain the semi-discrete scheme

 ∂fi,k∂t= −1|Ωi|(  3vk  3vk−3vk−3vk)fi,k (33) −2|Ωi|[(  0−2v+k  0v+k)fi−1,k+(−2v−k−v+kv−k2v+k)fi,k+(−v−k   02v−k   0)fi+1,k] +1τi(Mi,k−fi,k).

#### 3.2.2 Time discretization

As shown in , Discontinuous Galerkin scheme are unstable when used with a forward Euler method. Since we are interested in steady flows only, we use the same linearized implicit scheme as for the finite volume method (see [15, 1]).

#### 3.2.3 Boundary conditions

The Discontinuous Galerkin scheme is very compact, as the value only depends on ,