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

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 [23]. 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 [13]. 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 [18]. 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 [5], and more recently for the BGK equation [19]. 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 $f (t, x, v)$ defined such as $f (t, x, v) d x d v$ is the mass of molecules that at time $t$ are located in an elementary space volume $d x$ centered in $x = (x, y, z)$ and have a velocity in an elementary volume $d v$ centered in $v = (v_{x}, v_{y}, v_{z})$ . The macroscopic quantities associated to $f$ such as mass density $ρ$ , momentum $ρ u$ and total energy $E$

are defined as the five first moments of

f

with respect to the velocity variable, namely:

(ρ (t, x), ρ u (t, x), E (t, x)) = \int_{R^{3}} (1, v, \frac{1}{2} | v |^{2}) f (t, x, v) d v .

(1)

The temperature $T$ of the gas and its pressure $p = ρ R T$ are defined by the relation

E = \frac{1}{2} ρ | u |^{2} + \frac{3}{2} ρ R T

(2)

where $R$ 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 $f$ is a Gaussian function $M [ρ, u, T]$ of $v$ , called Maxwellian distribution, that depends only on macroscopic quantities and satisfies relations (1):

M [ρ, u, T] = \frac{ρ}{(2 π R T)^{\frac{3}{2}}} exp (- \frac{| v - u |^{2}}{2 R T}) .

(3)

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

\partial_{t} f + v \cdot \nabla_{x} f = Q (f),

(4)

which means that the total variation of $f$ (described by the left-hand side) is due to collisions between molecules ( $Q (f)$ 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) = \frac{1}{τ} (M [ρ, u, T] - f)

(5)

which makes $f$ relax towards the local equilibrium corresponding to the macroscopic quantities defined by (1). The relaxation time defined as $τ = \frac{μ}{ρ R T}$ ( $μ$ 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 $T_{w}$ and a velocity $u_{w}$ so that $u_{w} \cdot n (x) = 0$ , where $n (x)$ is the normal to the wall at point $x$

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

u_{w}

(

R T_{w}

being the variance of the normal distribution). This reads

f (t, x, v)_{| Γ} = σ_{w} (f) M [1, u_{w}, T_{w}] (v) if v \cdot n (x) > 0.

(6)

The parameter $σ_{w} (f)$ 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 [4]. Namely, that is,

σ_{w} (f) = - (\int_{v \cdot n (x) < 0} v \cdot n (x) f (t, x, v) d v) {(\int_{v \cdot n (x) > 0} v \cdot n (x) M [1, u_{w}, T_{w}] (v) d v)}_{v \cdot n (x) > 0}^{- 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 [10]). The velocity variable is replaced by discrete values $v_{k}$ of a Cartesian grid (see [1] for the use of locally refined velocity grids). The continuous distribution $f$ is then replaced by its approximation at each point $v_{k}$ , and we get the following discrete velocity BGK equation

\frac{\partial f_{k}}{\partial t} + v_{k} \cdot \nabla_{x} f_{k} = \frac{1}{} τ (M_{k} [ρ, u, T] - f_{k}),

(8)

where $f_{k} (t, x) \approx f (t, x, v_{k})$ . The discrete function $M_{k}$ is an approximation of the equilibrium functions $M$ defined in (3) at point $v_{k}$ . This discrete equilibrium is defined in order to satisfy the conservation of mass, momentum and energy (see [15]). The macroscopic quantities are now defined by the quadrature rule

(ρ (t, x), ρ u (t, x), E (t, x)) = \sum k (1, v_{k}, \frac{1}{2} | v_{k} |^{2}) f_{k} (t, x) ω_{k},

(9)

where $ω_{k}$ are the quadrature weights: for a Cartesian velocity grid of steps $Δ v_{x}$ , $Δ v_{y}$ , $Δ v_{z}$ , the weights are given by $ω_{k} = Δ v_{x} Δ v_{y} Δ v_{z}$ . Moreover, the boundary condition (6) is replaced by

f_{k} (t, x)_{| Γ} = σ_{w} (f) M_{k} [1, u_{w}, T_{w}] if v_{k} \cdot n (x) > 0,

(10)

where $σ_{w} (f)$ is defined as in (7) in which the integrals are replaced by quadratures:

σ_{w} (f) = - ⎛ ⎝ \sum v_{k} \cdot n (x) < 0 v_{k} \cdot n (x) f_{k} (t, x) ω_{k} ⎞ ⎠ {⎛ ⎝ \sum v_{k} \cdot n (x) > 0 v_{k} \cdot n (x) M_{k} [1, u_{w}, T_{w}] ω_{k} ⎞ ⎠}_{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 $[0, 1]$ . Equation (8) is then written

\frac{\partial f_{k}}{\partial t} + v_{k} \frac{\partial f_{k}}{\partial x} = \frac{}{1} τ (M_{k} [ρ, u, T] - f_{k}),

(12)

where $v_{k}$ is the first component of $v_{k}$ . Since $v_{k} \cdot n (x) = v_{k}$ (at $x = 0$ ), the boundary condition (10)–(11) now reads

f_{k} (t, 0) = σ_{w} (f) M_{k} [1, u_{w}, T_{w}] if v_{k} > 0,

(13)

where

σ_{w} (f) = - ⎛ ⎝ \sum v_{k} < 0 v_{k} f_{k} (t, 0) ω_{k} ⎞ ⎠ {⎛ ⎝ \sum v_{k} > 0 v_{k} M_{k} [1, u_{w}, T_{w}] ω_{k} ⎞ ⎠}_{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:

\begin{matrix} ϕ & := \sum k v_{k} f_{k} (t, 0) ω_{k} = \sum v_{k} < 0 v_{k} f_{k} (t, 0) ω_{k} + \sum v_{k} < 0 v_{k} σ_{w} (f) M_{k} [1, u_{w}, T_{w}] ω_{k} = 0 \end{matrix}

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 $Q (f)$ 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 $i_{m a x} + 1$ nodes $x_{i + \frac{1}{2}}$ with steps $Δ x_{i} = x_{i + \frac{1}{2}} - x_{i - \frac{1}{2}}$ , for $i = 1$ to $i = i_{m a x}$ . The discrete time variable is $t_{n}$ with a time step $Δ t_{n}$ . We first integrate equation (12) in a cell $[x_{i - \frac{1}{2}}, x_{i + \frac{1}{2}}]$ between $t_{n}$ and $t_{n + 1}$ and divide by $Δ t_{n}$ to obtain

\frac{f_{i, k}^{n + 1} - f_{i, k}^{n}}{Δ t_{n}} + \frac{1}{Δ x_{i}} (F_{i + \frac{1}{2}, k}^{n} - F_{i - \frac{1}{2}, k}^{n}) = Q_{k} (f_{i}^{n}),

(15)

where $f_{i, k}^{n}$ is an approximation of the average of $f_{k} (t_{n}, x)$ in the cell $i$ . The numerical flux $F_{i + \frac{1}{2}, k}^{n}$ is an approximation of the integral $\int_{t^{n}}^{t^{n + 1}} v_{k} f_{k} (t, x_{i + \frac{1}{2}}) d t$ , and can be interpreted as a flux across the cell interface $x_{i + \frac{1}{2}}$ between cells $i$ and $i + 1$ .

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:

F_{i + \frac{1}{2}, k}^{n} = v_{k}^{+} f_{i, k}^{n} + v_{k}^{-} f_{i + 1, k}^{n},

where $v_{k}^{\pm} = (v \pm | v |) / 2$ denotes the positive and negative parts of $v_{k}$ . 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 $[x_{\frac{1}{2}}, x_{\frac{3}{2}}]$ is adjacent to a solid wall. Then the numerical flux $F_{\frac{1}{2}, k}^{n}$ approximates the flux across the solid boundary. It requires an artificial value $f_{0, k}^{n}$ , that is generally interpreted as the value of the distribution function in a ghost cell $[x_{- \frac{1}{2}}, x_{\frac{1}{2}}]$ 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 $\sum_{k} F_{\frac{1}{2}, k}^{n} ω_{k} = 0$ .

These constraints can be satisfied as follows. First, the value of $f_{0, k}^{n}$ is defined for outgoing velocities $v_{k} < 0$ by using a zeroth order extrapolation of the value in $Ω_{1}$ , namely

f_{0, k}^{n} = f_{1, k}^{n} .

Then, the value of $f_{0, k}^{n}$ for incoming velocities is defined by using the boundary condition (13-14), that is to say:

f_{0, k}^{n} = σ_{w} (f_{0}^{n}) M_{k} [1, u_{w}, T_{w}], % where σ_{w} (f_{0}^{n}) = - (\sum l v_{l}^{-} f_{0, l}^{n} ω_{l}) / (\sum l v_{l}^{+} M_{l} [1, u_{w}, T_{w}] ω_{l}) .

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

\begin{matrix} \sum k F_{1 / 2, k}^{n} ω_{k} & = \sum k (v_{k}^{+} f_{0, k}^{n} + v_{k}^{-} f_{1, k}^{n}) ω_{k} = \sum k (- v_{k}^{+} \frac{\sum_{l} v_{l}^{-} f_{0, l}^{n} ω_{l}}{\sum_{l} v_{l}^{+} M_{l} [1, u_{w}, T_{w}] ω_{l}} M_{k} [1, u_{w}, T_{w}]) ω_{k} + \sum k (v_{k}^{-} f_{1, k}^{n}) ω_{k} = - \sum l (v_{l}^{-} f_{0, l}^{n}) ω_{l} + \sum k (v_{k}^{-} f_{1, k}^{n}) ω_{k} = 0, \end{matrix}

since $f_{0, k}^{n} = f_{1, k}^{n}$ for $v_{k} < 0$ .

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 [23]):

(16)

where the minmod function is defined by

minmod (x, y, z) = {\begin{matrix} sgn (x) min (| x |, | y |, | z |) if sgn (x) = sgn (y) = sgn (z), 0 in other cases . \end{matrix}

Like in the first order scheme, some ghost cell values have to be defined to compute the numerical flux at the wall interface $F_{\frac{1}{2}, k}^{n}$ . Indeed, the first order part of the flux requires the value of $f_{0, k}^{n}$ for incoming velocities ( $v_{k} > 0$ ), while the limiter needs $f_{0, k}^{n}$ and $f_{- 1, k}^{n}$ 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 $f_{0, k}^{n}$ and $f_{- 1, k}^{n}$ by a zeroth order extrapolation from the first boundary cell:

f_{- 1, k}^{n} = f_{0, k}^{n} = f_{1, k}^{n} .

For incoming velocities, we define $f_{0, k}^{n}$ like in the first order scheme, that is to say:

f_{0, k}^{n} = σ_{w} (f_{0, k}^{n}) M_{k}, where σ_{w} (f_{0, k}^{n}) = - (\sum l v_{l}^{-} f_{0, l}^{n} ω_{l}) {(\sum l v_{l}^{+} M_{l} ω_{l})}_{l}^{- 1} .

And again, $f_{- 1, k}^{n}$ is defined by zeroth order extrapolation with $f_{- 1, k}^{n} = f_{0, k}^{n}$ .

This definition implies that the numerical flux at the wall interface $F_{\frac{1}{2}, k}^{n}$ reduces to the first order flux $v_{k}^{+} f_{0, k}^{n} + v_{k}^{-} f_{1, k}^{n}$ , since $Δ f_{- \frac{1}{2}, k}^{n} = f_{0, k}^{n} - f_{- 1, k}^{n} = 0$ by construction, and hence the limiter $Φ_{1 / 2, k}^{n}$ 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 $f_{0, k}^{n}$ and $f_{- 1, k}^{n}$ 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 $f_{k}$ inside each cell [14]:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} F_{i + \frac{1}{2}, k}^{n} = v_{k}^{+} f_{i + \frac{1}{2}, k}^{n, -} + v_{k}^{-} f_{i + \frac{1}{2}, k}^{n, +}, f_{i + \frac{1}{2}, k}^{n, -} = f_{i, k}^{n} + \frac{Δ x_{i}}{2} δ_{i, k}^{n}, f_{i + \frac{1}{2}, k}^{n, +} = f_{i + 1, k}^{n} - \frac{Δ x_{i + 1}}{2} δ_{i + 1, k}^{n}, \end{matrix}

(17)

where $δ_{i, k}^{n}$ is the slope of the linear reconstruction of $f$ in cell $i$ . The most classical way to compute this slope is a least square method, which gives, for a regular mesh, the centered slope

δ_{i, k}^{n} = \frac{f_{i + 1, k}^{n} - f_{i - 1, k}^{n}}{2 Δ x_{i}} .

(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

δ_{i, k}^{n, l i m} = m i n m o d (δ_{i, k}^{n}, 2 α_{i} \frac{f_{i, k}^{n} - f_{i - 1, k}^{n}}{Δ x_{i}}, 2 α_{i} \frac{f_{i + 1, k}^{n} - f_{i, k}^{n}}{Δ x_{i}}),

(19)

where $α_{i}$ is a free parameter between 0 and 1 [14].

When this scheme is applied to the first cell $i = 1$ , the numerical flux $F_{\frac{1}{2}, k}^{n}$ is

F_{\frac{1}{2}, k}^{n} = v_{k}^{+} f_{\frac{1}{2}, k}^{n, -} + v_{k}^{-} f_{\frac{1}{2}, k}^{n, +} .

It uses two values of the distribution at the solid interface: $f_{\frac{1}{2}, k}^{n, +}$ is the value on the right side of the wall, required for outgoing velocities $v_{k} < 0$ , given by the linear reconstruction

f_{\frac{1}{2}, k}^{n, +} = f_{1, k}^{n} - \frac{Δ x_{1}}{2} δ_{1, k}^{n, l i m},

(20)

while $f_{\frac{1}{2}, k}^{n, -}$ is the value on the left side of the wall, required for incoming velocities $v_{k} > 0$ , given by the linear reconstruction

f_{\frac{1}{2}, k}^{n, -} = f_{0, k}^{n} + \frac{Δ x_{0}}{2} δ_{0, k}^{n, l i m},

(21)

where the slopes are defined by (18)–(19) with $i = 1$ and $i = 0$ , respectively. This requires to define two ghost cell values $f_{0, k}^{n}$ and $f_{- 1, k}^{n}$ , and the corresponding cell size $Δ x_{0}$ . 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 $f$ at the wall interface satisfies

f_{\frac{1}{2}, k}^{n, -} = σ_{w} (f_{\frac{1}{2}, k}^{n, +}) M_{k} .

(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 $f_{\frac{1}{2}, k}^{n, \pm}$ by extrapolation (for outgoing velocities) and by the boundary condition (for incoming velocities), and then in defining the ghost cell values $f_{0, k}^{n}$ and $f_{- 1, k}^{n}$ 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: $Δ x_{- 1} = Δ x_{0} = Δ x_{1}$ .

$∙$ Step 1: computation of $f_{1 / 2, k}^{n, +}$ for outgoing velocities $v_{k} < 0$ .
This value is defined by a linear extrapolation using the values $f_{1, k}^{n}$ and $f_{2, k}^{n}$ :

f_{1 / 2, k}^{n, +} = \frac{3}{2} f_{1, k}^{n} - \frac{1}{2} f_{2, k}^{n} .

(23)

$∙$ Step 2: computation of $f_{1 / 2, k}^{n, -}$ for incoming velocities $v_{k} > 0$ .
We use the boundary condition (13)–(14) to get

f_{1 / 2, k}^{n, -} = σ_{w} (f_{1 / 2}^{n, +}) M_{k}^{n} .

(24)

$∙$ Step 3: definition of the values $f_{0, k}^{n}$ and $f_{- 1, k}^{n}$ in the ghost cells.

For outgoing velocities, only the value of $f_{0, k}^{n}$ is required. We use the same linear extrapolation as used to compute extrapolation $f_{1 / 2, k}^{n, +}$ (see Step 1) to get

f_{0, k}^{n} = 2 f_{1, k}^{n} - f_{2, k}^{n} .

For incoming velocities, we need the values of $f_{0, k}^{n}$ and $f_{- 1, k}^{n}$ . Then we use again a linear extrapolation, but based on the values of the incoming value of the wall interface distribution $f_{1 / 2, k}^{n, -}$ and on $f_{1, k}^{n}$ . This gives

	$f_{0, k}^{n}$	$= 2 f_{1 / 2, k}^{n, -} - f_{1, k}^{n},$
	$f_{- 1, k}^{n}$	$= 4 f_{1 / 2, k}^{n, -} - 3 f_{1, k}^{n} .$

Now it is not difficult to prove that these definitions satisfy the previous constraints, provided that $α_{i} \geq \frac{1}{2}$ . First, (22) is imposed at step 2, and hence it is satisfied, which gives mass conservation. Now, it remains to prove that $f_{\frac{1}{2}, k}^{n, \pm}$ defined by (23) and (24) also satisfy (20) and (21). This is due to the extrapolation procedures that make the points $(x_{i}, f_{i, k}^{n})$ are on the same straight line, and hence the limiters can be computed. Indeed, for instance, we have for incoming velocities $v_{k} > 0$ :

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

\begin{matrix} f_{0, k}^{n} + \frac{Δ x_{0}}{2} δ_{0, k}^{n, l i m} & = 2 f_{1 / 2, k}^{n, -} - f_{1, k}^{n} + f_{1, k}^{n} - f_{1 / 2, k}^{n, -} = f_{1 / 2, k}^{n, -} \end{matrix}

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

\begin{matrix} δ_{1, k}^{n, l i m} & = m i n m o d (\frac{f_{2, k}^{n} - f_{0, k}^{n}}{2 Δ x_{1}}, 2 α_{1} \frac{f_{1, k}^{n} - f_{0, k}^{n}}{Δ x_{1}}, 2 α_{1} \frac{f_{2, k}^{n} - f_{1, k}^{n}}{Δ x_{1}}) = m i n m o d (\frac{f_{2, k}^{n} - f_{1, k}^{n}}{Δ x_{1}}, 2 α_{1} \frac{f_{2, k}^{n} - f_{1, k}^{n}}{Δ x_{1}}, 2 α_{1} \frac{f_{2, k}^{n} - f_{1, k}^{n}}{Δ x_{1}}) = \frac{f_{2, k}^{n} - f_{1, k}^{n}}{Δ x_{1}} . \end{matrix}

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

\begin{matrix} f_{1, k}^{n} - \frac{Δ x_{1}}{2} δ_{1, k}^{n, l i m} & = \frac{3}{2} f_{1, k}^{n} - \frac{1}{2} f_{2, k}^{n} = f_{1 / 2, k}^{n, +} \end{matrix}

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 $f_{0, k}$ and $f_{- 1, k}$ 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 [23]). This implies that for positive velocities, $f_{- 1, k}^{n}$ and $f_{0, k}^{n}$ are defined through extrapolation of $f_{\frac{1}{2}, k}^{-, n}$ and $f_{1, k}^{n}$ , but they cannot be aligned with $f_{2, k}^{n}$ , 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 [14]), 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 $Ω_{i} = [x_{i - \frac{1}{2}}, x_{i + \frac{1}{2}}]$ , we use two basis affine functions $φ_{i, 1}$ and $φ_{i, 2}$ such that $φ_{i, 1} (x_{i - \frac{1}{2}}) = 1$ , $φ_{i, 1} (x_{i + \frac{1}{2}}) = 0$ , and $φ_{i, 2} (x_{i - \frac{1}{2}}) = 0$ , $φ_{i, 2} (x_{i + \frac{1}{2}}) = 1$ . We project equation (8) in this basis and integrate by parts to get for every $i$ and $p = 1, 2$ :

\begin{matrix} \int_{Ω_{i}} \frac{\partial f_{k}}{\partial t} (t, x) φ_{i, p} (x) d x - \int_{Ω_{i}} v_{k} f_{k} (t, x) \partial_{x} φ_{i, p} (x) d x + (v_{k} f_{k} (t, x_{i + \frac{1}{2}}) φ_{i, p} (x_{i + \frac{1}{2}}) - v_{k} f_{k} (t, x_{i - \frac{1}{2}}) φ_{i, p} (x_{i - \frac{1}{2}})) = \int_{Ω_{i}} \frac{1}{τ_{i}} (M_{k} [ρ, u, T] - f_{k}) φ_{i, p} (x) d x, \end{matrix}

(25)

Now, we assume that $f_{k} (t, x)$ can be approximated on each cell by a piecewise linear function ${^f}_{k} (t, x)$ defined by ${^f}_{k} (t, x) |_{Ω_{i}} = f_{i, k} (t, x) = f_{i, k, 1} (t) φ_{i, 1} (x) + f_{i, k, 2} (t) φ_{i, 2} (x)$ (see figure 3). Note that while the components $f_{i, k, 1}$ and $f_{i, k, 2}$ are also the pointwise values of $f_{i, k}$ at the edges of $Ω_{i}$ , the whole function ${^f}_{k}$ itself is not continuous across these edges: indeed, at each cell edge $x_{i + \frac{1}{2}}$ , ${^f}_{k}$ has two left and right values $f_{i, k} (t, x_{i + \frac{1}{2}}) = f_{i, k, 2}$ and $f_{i + 1, k} (t, x_{x_{i + \frac{1}{2}}}) = f_{i + 1, k, 1}$ (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 \sum q = 1 \int_{Ω_{i}} \frac{\partial f_{i, k, q}}{\partial t} φ_{i, q} (x) φ_{i, p} (x) d x = 2 \sum q = 1 m_{p q}^{i} \frac{\partial f_{i, k, q}}{\partial t},

(26)

where the $(m_{p q}^{i})$ form the $2 \times 2$ matrix

(m_{p q}^{i}) = ¯_{i} = \frac{| Ω_{i} |}{6} (\begin{matrix} 2 & 1 1 & 2 \end{matrix}) .

(27)

The second term of the left-hand side is

- v_{k} \int_{Ω_{i}} f_{i, k} (x) \frac{\partial φ_{i, p} (x)}{\partial x} d x = 2 \sum q = 1 d_{p q} f_{i, k, q},

(28)

where the $(d_{p q})$ form the $2 \times 2$ matrix

(d_{p q}) = {¯ ¯ D}_{k} = \frac{v_{k}}{2} (\begin{matrix} 1 & 1 - 1 & - 1 \end{matrix}) .

(29)

For the third term of the left-hand side of (25), we muse take into account that the piecewise linear approximation of $f_{k}$ is not continuous across cell edges (see the remark above). Consequently, we write this term as $v_{k} {^f}_{k} (t, x_{i + \frac{1}{2}}) φ_{i, p} (x_{i + \frac{1}{2}}) - v_{k} {^f}_{k} (t, x_{i - \frac{1}{2}}) φ_{i, p} (x_{i - \frac{1}{2}})$ , where a single value for ${^f}_{k} (t, x_{i + \frac{1}{2}})$ 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 $φ_{i, p}$ take values $0$ or $1$ at the cell edges of $Ω_{i}$ , 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

\int_{C_{i}} \frac{1}{τ} (M_{k} [^U] - {^f}_{k}) φ_{i, p} (x) d x = 2 \sum q = 1 m_{p q}^{i} \frac{1}{τ_{i, q}} (M_{i, k, q} - f_{i, k, q}),

(31)

where $m_{p q}^{i}$ has been defined in (27), $M_{i, k, q} = M_{k} [ρ_{i, q}, u_{i, q}, T_{i, q}]$ , $τ_{i, q} = τ (ρ_{i, q}, T_{i, q})$ , and the macroscopic quantities $ρ_{i, q}, u_{i, q}, T_{i, q}$ are defined by applying (9) in each cell.

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

¯_{i} \frac{\partial f_{i, k}}{\partial t} + {¯ ¯ D}_{k} f_{i, k} + ({¯ ¯ A}_{k} f_{i - 1, k} + {¯ ¯ B}_{k} f_{i, k} + {¯ ¯ C}_{k} f_{i + 1, k}) = ¯_{i} \frac{1}{τ_{i}} (M_{i, k} - f_{i, k}),

(32)

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

f_{i, k} = (f_{i, k, 1}, f_{i, k, 2})^{T}

and

M_{i, k} = (M_{i, k, 1}, M_{i, k, 2})^{T}

, the 2

\times

2 matrices

{¯ ¯ A}_{k} = (\begin{matrix} 0 & - v_{k}^{+} 0 & 0 \end{matrix})

{¯ ¯ B}_{k} = (\begin{matrix} - v_{k}^{-} & 0 0 & v_{k}^{+} \end{matrix})

{¯ ¯ C}_{k} = (\begin{matrix} 0 & 0 v_{k}^{-} & 0 \end{matrix})

, while the matrices

¯_{i}

and

{¯ ¯ D}_{k}

have been defined in (27) and (29). Note that the product of vectors

\frac{1}{τ_{i}} (M_{i, k} - f_{i, k})

must be understood component-wise. Multiplying (32) by the matrix

{¯_{i}}^{- 1} = \frac{2}{| Ω_{i} |} (\begin{matrix} 2 & - 1 - 1 & 2 \end{matrix})

(the inverse of

¯_{i}

), we obtain the semi-discrete scheme

$\frac{\partial f_{i, k}}{\partial t} =$	$- \frac{1}{\| Ω_{i} \|} (\begin{matrix} 3 v_{k} & 3 v_{k} - 3 v_{k} & - 3 v_{k} \end{matrix}) f_{i, k}$	(33)
	$- \frac{2}{\| Ω_{i} \|} [(\begin{matrix} 0 & - 2 v_{k}^{+} 0 & v_{k}^{+} \end{matrix}) f_{i - 1, k} + (\begin{matrix} - 2 v_{k}^{-} & - v_{k}^{+} v_{k}^{-} & 2 v_{k}^{+} \end{matrix}) f_{i, k} + (\begin{matrix} - v_{k}^{-} & 0 2 v_{k}^{-} & 0 \end{matrix}) f_{i + 1, k}]$
	$+ \frac{1}{τ_{i}} (M_{i, k} - f_{i, k}) .$

3.2.2 Time discretization

As shown in [7], 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 $f_{i, k}^{n + 1}$ only depends on $f_{i, k}^{n}$ ,

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

Authors

Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions for the compressible Navier-Stokes equations

Extension and comparison of techniques to enforce boundary conditions in Finite Volume POD-Galerkin reduced order models for fluid dynamic problems

Discontinuous Galerkin Approach to Large Bending Deformation of a Bilayer Plate with Isometry Constraint

One-sided GRP Solver and Numerical Boundary Conditions for compressible fluid flows

One- and multi-dimensional CWENOZ reconstructions for implementing boundary conditions without ghost cells

Mathematical aspects relative to the fluid statics of a self-gravitating perfect-gas isothermal sphere

Discontinuous Galerkin method for macroscopic traffic flow models on networks

1 Introduction

2 The BGK equation

3 Numerical schemes and boundary conditions

3.1 Finite Volume schemes

3.1.1 A first order scheme finite volume scheme

3.1.2 Second order finite volume scheme with flux limiters

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

Remark 3.1.

3.2 Discontinuous Galerkin scheme

3.2.1 Weak form and piecewise linear approximation

3.2.2 Time discretization

3.2.3 Boundary conditions

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

Authors

Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions for the compressible Navier-Stokes equations

Extension and comparison of techniques to enforce boundary conditions in Finite Volume POD-Galerkin reduced order models for fluid dynamic problems

Discontinuous Galerkin Approach to Large Bending Deformation of a Bilayer Plate with Isometry Constraint

One-sided GRP Solver and Numerical Boundary Conditions for compressible fluid flows

One- and multi-dimensional CWENOZ reconstructions for implementing boundary conditions without ghost cells

Mathematical aspects relative to the fluid statics of a self-gravitating perfect-gas isothermal sphere

Discontinuous Galerkin method for macroscopic traffic flow models on networks

This week in AI

1 Introduction

2 The BGK equation

3 Numerical schemes and boundary conditions

3.1 Finite Volume schemes

3.1.1 A first order scheme finite volume scheme

3.1.2 Second order finite volume scheme with flux limiters

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

Remark 3.1.

3.2 Discontinuous Galerkin scheme

3.2.1 Weak form and piecewise linear approximation

3.2.2 Time discretization

3.2.3 Boundary conditions