AP schemes including boundary layers

# Micro-macro schemes for kinetic equations including boundary layers

Mohammed Lemou CNRS and IRMAR, Université de Rennes 1, France  and  Florian Méhats IRMAR, Université de Rennes 1, France
###### Abstract.

We introduce a new micro-macro decomposition of collisional kinetic equations in the specific case of the diffusion limit, which naturally incorporates the incoming boundary conditions. The idea is to write the distribution function in all its domain as the sum of an equilibrium adapted to the boundary (which is not the usual equilibrium associated with ) and a remaining kinetic part. This equilibrium is defined such that its incoming velocity moments coincide with the incoming velocity moments of the distribution function. A consequence of this strategy is that no artificial boundary condition is needed in the micro-macro models and the exact boundary condition on is naturally transposed to the macro part of the model. This method provides an ’Asymptotic preserving’ numerical scheme which generates a very good approximation of the space boundary values at the diffusive limit, without any mesh refinement in the boundary layers. Our numerical results are in very good agreement with the exact so-called Chandrasekhar value, which is explicitely known in some simple cases.

## 1. Introduction

The development of numerical methods to solve multiscale kinetic equations has been the subject of active research in the past years, with applications in various fields: plasma physics, rarefied gas dynamics, aerospace engineering, semiconductors, radiative transfert, … The general problem is to construct numerical schemes that are able to capture the properties of the various scales in the considered system, while the numerical parameters remain as independent as possible of the stiffness character of these scales.

There is a huge literature dealing with the construction of such schemes, and in particular it is now a usual challenge to design the so-called Asymptotic Preserving (AP) [18] numerical schemes for kinetic equations which are consistent with the kinetic model for all positive value of the Knudsen number , and degenerate into consistent schemes with the asymptotic models (compressible Euler and Navier-Stokes, diffusion, etc) when . Among a long list of works on this topic, we can mention [30, 31, 19, 20, 15, 23, 17] for stationary problems and [22, 21, 24, 27, 25, 26, 16, 33, 35, 11, 5, 4, 29] for time dependent problems.

Here, we are interested in boundary value problems and our aim is to develop a strategy to construct numerical schemes that are AP in the usual sense and are also able to deal with the space boundary conditions and kinetic boundary layers. In this paper, we focus on the diffusive scaling but our strategy could be extended to other asymptotics. We emphasize that in general, the known numerical schemes induce at the limit some numerical boundary conditions which are different from the theoretical limiting boundary values. In most cases, this does not affect only the boundary layer, but also the solution inside the domain.

Boundary conditions in the construction of AP schemes for transient kinetic equation in the diffusive scaling have been taken into account in several works. For instance, in [24, 25, 26], Klar uses an approximation of the solution to the Milne problem which is injected as a boundary condition for the evolution problem. The boundary value used by Jin, Pareschi and Toscani in [22] is also derived from the asymptotic analysis of the problem as the Knudsen number goes to zero. Of course, taking the boundary condition from the limiting model cannot provide a good approximation of the problem for all regime, in particular for large values of . Our goal here is to develop a new approach in which the boundary conditions are not extracted from the limiting model but are directly derived from the original kinetic boundary condition, at any time. Our numerical results will be compared with those obtained with the approaches in [24, 22, 33].

To this purpose, our starting idea is to extend the micro-macro decomposition of [2, 33, 34] in order to incorporate in a natural way the exact inflow boundary conditions. The general micro-macro decomposition method consists in writing the kinetic equation as a system coupling a macroscopic part (say a Maxwellian) and a remaining kinetic part. The main advantage of this method is its robustness and easy adaptability, since it can be applied to various asymptotics (fluid, diffusion, high-field, etc) and to a large class of collision operators. However, the general problem of space boundary condition and boundary layers has not been completely solved in this approach. When incoming boundary conditions on the distribution function are prescribed, it is clear that this cannot be translated into a boundary condition on the macro part (say the total density) in an exact way since the distribution function is only known for incoming velocities. Such boundary conditions for both macro and micro parts of the distribution function are needed for the computation of fluxes. In this case, an artificial boundary condition of Neumann type is used in [33] and leads to some unsatisfactory approximation of the Chandrasekhar value at the diffusive limit.

Our first aim is to develop a new micro-macro decomposition of collisional kinetic equations which naturally incorporates the exact space boundary conditions (BC). The second task in this work is to show how to use this new formulation at the boundary in order to capture the right boundary conditions and boundary layers in the limit of a small Knudsen number.

To achieve this goal we proceed as follows. We consider a kinetic equation with a linear collision operator in a diffusive scaling on a bounded domain, subject to inflow boundary conditions. Our first idea is to decompose the distribution function in its domain as the sum of a Maxwellian part adapted to the boundary (which is not the usual Maxwellian associated with ) and a remaining kinetic part. This Maxwellian is defined on the whole domain such that its ’incoming’ velocity moments coincide with the ’incoming’ velocity moments of the distribution function. In this way, the exact boundary condition can be directly injected into the micro-macro system. The second idea is to use a suitable approximation of the space derivative of the density at the boundary which garantees a balance between the incoming and outgoing fluxes at this boundary. Important consequences of this strategy are the following. No artificial boundary condition is needed in the micro/macro models and the exact boundary conditions on are naturally shared by the macro part and the kinetic part. A further fundamental property of the strategy developed here is that it provides AP schemes inside the physical domain and a very good approximation in the space boundary layers. Indeed, at the end of this paper, we provide various numerical tests which illustrate this property. We emphasize that the numerical boundary value given by our scheme when is small is very close to the theoretical boundary value derived from the Chandrasekhar function, without injecting this value in our scheme. The strategy developed in this paper has been announced in a preliminary note [32], where we showed that this method can also be applied in principle for a hydrodynamic scaling. We will explore in the future several extensions of this work, in particular the possibility of dealing with more complex collision operators. Our approach shares some similarities with the method developed in [10, 12] where a partial moment system with entropy closure was used to approximate the original kinetic equation.

The paper is organized as follows. In Section 2, the general strategy is described at the continuous level. We introduce a new micro-macro decomposition which is adapted to boundary value problems for kinetic equations. In Section 3, we use this decomposition to design a numerical scheme which is AP in the diffusion limit. In a first step (Subsection 3.1), we show how this formulation allows to inject directly the inflow boundary condition in the numerical method, avoiding in this way any artificial boundary condition. In a second step (Subsection 3.2), we introduce a specific (still consistent) discretization of the spatial derivative of the density in order to correctly capture the boundary condition and the boundary layer. The properties of our scheme are summarized in Proposition 3.2. In Section 4, we provide various numerical tests to validate our approach. We test different boundary conditions (with or without boundary layers) and different linear collision operator (with local and non local collision kernels). Finally, we end the paper by a conclusion and some perspectives.

## 2. A boundary matching micro-macro decomposition

### 2.1. The diffusion limit of kinetic equations in bounded domains

Let be a domain in the position space and let be a domain in the velocity space , being endowed with a measure . At any point of the boundary of , we denote by the unitary outgoing normal vector to . We then consider the following linear transport equation written in a diffusive scaling:

 ε∂tf+v⋅∇xf=1εLf,t>0,(x,v)∈Ω×V, (2.1)

with initial condition

 f(0,x,v)=finit(x,v),(x,v)∈Ω×V (2.2)

and with inflow boundary conditions

 f(t,x,v)=fb(t,x,v),t>0,(x,v)∈∂Ω×V % such that v⋅n(x)<0. (2.3)

The unknown is the distribution function of the particles that depends on time , on position , and on velocity .

In (2.1), the linear collision operator acts on the velocity dependence of and describes the interactions of particles with the medium. This operator relaxes the system of particles to an equilibrium , which is supposed to be a positive and even function. For all distribution function , we shall denote

 ⟨h⟩V=∫Vh(v)dv∫VE(v)dv. (2.4)

Note in particular that and . We assume that the collision operator is non-positive and self-adjoint in , with null space and range given by

 N(L)=Span{E}

and

 R(L)=(N(L))⊥={f such that ⟨f⟩V=0}.

Hence, an important property of is the fact that it is invertible on , we shall denote its pseudo-inverse by . Additionally, we assume the invariance of under orthogonal transformations of . The parameter measures a dimensionless mean free path of the particles, or also the inverse of the dimensionless observation time.

When goes to 0 in (2.1), goes formally to an equilibrium state . The diffusion limit is the equation satisfied by the limiting density . This equation is a diffusion equation and has been obtained in various situations thanks to asymptotic expansion in (Hilbert or Chapman-Enskog expansion), see e.g. [28, 3, 1, 9, 36]. We have

 ∂tρ0−∇x⋅(κ∇xρ0)=0% withκ=−⟨vL−1(vE)⟩V>0, (2.5)

with a Dirichlet boundary condition

 ρ0(t,x)=ρb(t,x),t>0,x∈∂Ω. (2.6)

The limit value is usually obtained by a boundary layer analysis and its computation involves the resolution of a kinetic half-space problem, the Milne problem associated with the collision operator . For , let be the bounded solution of the half-space problem

 −v⋅n(x)∂yχt,x = Lχt,x,y>0,v∈V, χt,x(0,v) = fb(t,x,v),v⋅n(x)<0, (2.7)

in which and are parameters. Then, under reasonnable conditions of , this problem is well-posed [3, 1, 14, 36] and we have

 limy→+∞χt,x(y,v)=ρb(t,x)E(v). (2.8)

Notice that, in the special case where the incoming data is an equilibrium state, i.e. when , one has and is independent of . In this case, there is no boundary layer near as . For general incoming distribution functions, the question of determining requires the resolution of the half-space problem (2.1), which may be computationaly demanding. Our aim here is to construct a numerical method which gives a good approximation of the boundary layer without solving the half-space problem (2.1).

In order to check the efficiency of our approach, a validation test will be the following. We will compare our numerical results with the exact Dirichlet value defined by (2.8) in the diffusion limit , in special cases where this value is well-known. The most simple case is the one-group transport equation in slab geometry, which is such that , , the velocity set is , and . Then is given by

 ρb(t,0)=∫10K0(v)fb(t,0,v)dv,ρb(t,1)=∫0−1K0(−v)fb(t,1,v)dv, (2.9)

where

 K0(v)=√32vH(v), (2.10)

the function being the so-called Chandrasekhar function [7, 1, 8], see also [6], which can be computed numerically from the following formula:

 H(v)=1+vH(v)2∫10H(w)v+wdw. (2.11)

A good approximation of the function is usually given by

 K0(v)≃K1(v)=32v2+v, (2.12)

see [30] for instance.

### 2.2. A new micro-macro decomposition

Let us first recall the micro-macro decomposition introduced in [33] and [34]. We decompose

 f=ρE+g1, (2.13)

with

 ρ(t,x)=⟨f⟩Vandg1=f−ρE

( is not necessarily small). We denote by the orthogonal projector in onto the nullspace of , i.e.

 Πϕ=⟨ϕ⟩VE.

Inserting the decomposition (2.13) into the kinetic equation (2.1) and applying and successively ( being the identity operator), one gets

 ∂tρ+1ε∇x⋅⟨vg1⟩V=0, ∂tg1+1ε(I−Π)(v⋅∇xg1)=1ε2[Lg1−εEv⋅∇xρ)].

We now emphasize that in this formulation, the space boundary condition on and are not known because they cannot be inferred from the boundary conditions on in general. Indeed, (2.3) only gives the function at the boundary for incoming velocities, i.e. such that , and it is therefore clear that the values of cannot be determined a priori on the boundary with the only knowledge of . Consequently, the micro-macro decomposition (2.13) appears as a method for capturing the diffusion limit inside the domain only. The numerical scheme constructed in [33] is Asymptotic Preserving in a strict sense (i.e. including the boundary) only for incoming distribution data at the equilibrium , when there is no boundary layer near . Indeed, only a rough approximation of the boundary layer is obtained for non equilibrium boundary conditions.

Our aim in this work is to develop a new micro-macro decomposition which is able to incorporate the boundary conditions in an exact way. This will allow us in the next section to construct a numerical method with natural boundary conditions, that approximate very well the diffusion limit even for incoming data that are not an equilibrium state.

To this purpose, let us introduce a few further notations. We consider a function which extends the function (defined on ) to the whole domain . Let us give explicit examples of for specific geometries. For the unit ball centered at the origin, one can take . For a half plane , one can choose . In dimension one, for the interval , one can take . Then, for all , we split the velocity space into two parts according to the sign of this function :

 V−(x)={v∈V, ω(x,v)<0},  V+(x)=V∖V−(x). (2.14)

For all function , we will denote

 ⟨h⟩V−=∫V−h(v)dμ∫V−E(v)dμ,⟨h⟩V+=∫V+h(v)dμ∫V+E(v)dμ

and

 ΠV−h=⟨h⟩V−E.

The idea is now the following. Instead of looking for an equation on as usual, we seek an equation on the following ’boundary matching’ density:

 ¯¯¯ρ(t,x)=⟨f(t,x,⋅)⟩V− (2.15)

and perform the corresponding micro-macro decomposition:

 f=¯¯¯ρE+g=ΠV−f+g. (2.16)

In particular, we have . When , we know that the solution of (2.1) is, at least formally, close to except in initial or boundary layers. Therefore, since we have the moment will also be close to for small . This shows that the new decomposition (2.16) is still a decomposition of into an asymptotic part (macro part) and a kinetic part (micro part). In order to derive the system of equations satisfied by and , we first integrate (2.1) on and get

 (2.17)

Substracting (2.17) from equation (2.1) on , we obtain the equation on :

 ∂tg+1ε(I−ΠV−)(v⋅∇xg)+1ε(I−ΠV−)(vE)⋅∇x¯¯¯ρ=1ε2(I−ΠV−)(Lg). (2.18)

Finally, we remark that the system (2.17), (2.18) can be replaced by a more convenient (and still equivalent) system in terms of and as follows:

 ∂tρ+1ε⟨v⋅∇xg⟩V=0, (2.19) ∂tg+1ε(I−ΠV−)(v⋅∇xg)+1ε(I−ΠV−)(vE)⋅∇x¯¯¯ρ=1ε2(I−ΠV−)(Lg). (2.20)

Note that is linked to and by the relations

 ¯¯¯ρ=ρ−⟨g⟩V,f=ρE+g−⟨g⟩VE. (2.21)

One main interest of this new micro-macro formulation (2.16) is the following: the fluxes involved in (2.19) and (2.20) only concern the quantities or , and not . This means that numerical schemes of this formulation would only need the values of and at the space boundary which are completely known from the original boundary condition (2.3) on :

 ¯¯¯ρ(t,x,v)=⟨fb(t,x,⋅)⟩V−,t>0,(x,v)∈∂Ω, (2.22) g(t,x,v)=fb(t,x,v)−⟨fb(t,x,⋅)⟩V−E(v),t>0,(x,v)∈∂Ω×V−. (2.23)

Note that a similar decomposition as (2.16) can be performed for a fluid scaling as well, see [32].

## 3. The numerical method

In this section, we construct in dimension a numerical scheme for (2.19), (2.20), see Proposition 3.2, which is uniformly stable in the limit and provides a good approximation of the boundary layers, without introducing any artificial boundary condition. To this purpose, we proceed in two steps. First, we present in subsection 3.1 a preliminary numerical scheme, see Proposition 3.1, which has the desired uniform stability and does not use any artificial boundary condition. However, it turns out that this scheme does not give a good approximation of the limiting boundary condition. To remedy this problem, we then introduce our numerical scheme in subsection 3.2 where a specific discretization is done near the boundary. We emphasize that this scheme is able to reproduce an accurate approximation of the boundary layer without any mesh refinement.

Let us introduce a few notations. The domain in space is the interval . The incoming boundary conditions are

 f(t,0,v)=fℓ(v), for v∈V−(0)={v∈V,v>0}, (3.1) f(t,1,v)=fr(v), for v∈V−(1)={v∈V,v<0}, (3.2)

the data and being independent of time here for simplicity. The function is a smooth function satisfying

 ω(0,v)=−v,ω(1,v)=v. (3.3)

Let be a constant time step and be a uniform space step. We set and consider staggered grids for and for . We will construct approximations of the densities and on the grid and approximations of the remainder on the grid . The velocity variable is kept continuous in this section, and its discretization will be made precise later on. At the time , the initial condition (2.2) induces naturally the initialization

 ρ0i=⟨finit(xi,⋅)⟩V,¯¯¯ρ0i=⟨finit(xi,⋅)⟩V−,g0i+1/2=(I−ΠV−)finit(xi+1/2,⋅). (3.4)

Note that, by construction, the function is known at the boundary , and can be deduced from the boundary conditions (3.1), (3.2):

 ¯¯¯ρ0=∫V−(0)fℓ(v)dμ∫V−(0)E(v)dμ,¯¯¯ρN+1=∫V−(1)fr(v)dμ∫V−(1)E(v)dμ. (3.5)

At the boundary, the function is also known for entering velocities. We introduce two additional unknowns and for . The incoming boundary conditions yield

 gn0(v)=fℓ(v)−¯¯¯ρ0E for v∈V−(0) gnN+1(v)=fr(v)−¯¯¯ρN+1E for v∈V−(1). (3.6)

For outgoing velocities, these functions will be provided by our numerical scheme, when their values are required, and no artificial condition will be needed. This is explained in Subsection 3.2.

### 3.1. A first attempt without specific treatment at the boundary

We assume that we know and for , and for . We proceed by recursion and explain how to compute these same quantities at time . First, for we discretize (2.20) by a semi-implicit upwind scheme:

 gn+1i+1/2−gni+1/2Δt + 1ε(I−ΠV−)(v+gni+1/2−gni−1/2Δx+v−gni+3/2−gni+1/2Δx) (3.7) + 1ε(I−ΠV−)(vE)(¯¯¯ρni+1−¯¯¯ρniΔx)=1ε2(I−ΠV−)(Lgn+1i+1/2),

where and . Note that the source term involving is discretized by a centered finite difference scheme. We observe that and for incoming velocities can be extracted from the relations

 gn0=gn−1/2+gn1/22,gnN+1=gnN+1/2+gnN+3/22, (3.8)

since and are known and since and are known for incoming velocities by (3.6). Note also that and are known by (3.5). Therefore, the scheme (3.7) enables to compute all the desired values for (see below in the proof of Proposition 3.1).

It remains to compute and for . To that purpose, we consider the following discretization of (2.19). For , we set

 ρn+1i−ρniΔt+1ε⟨vgn+1i+1/2−gn+1i−1/2Δx⟩V=0. (3.9)

These relations clearly provide all the desired values for . Then we compute all the values according to (2.21), by setting

 ¯¯¯ρn+1i=ρn+1i−⟨gn+1i+1/2+gn+1i−1/22⟩V (3.10)

for . The recursion is now complete and can be iterated. Observe that, in this first attempt, we do not use the outgoing values of and in order to compute and inside the domain.

We summarize in the following proposition some properties of this scheme.

###### Proposition 3.1 (Formal).

The numerical scheme (3.7), (3.8), (3.9), (3.10) with initial and boundary conditions (3.4), (3.5), (3.6) determines all the values and for , and for , at any time interation . This scheme is stable under the CFL condition

 Δt≤C(Δx2+εΔx),

where is a constant independent of , and . For all fixed , this scheme is consistent with the initial and boundary value problem (2.1), (2.2), (2.3). Moreover, for all fixed and , the scheme degenerates as into a discretization of the diffusion equation (2.5). More precisely, for , we have as , and is defined by the scheme

 ˜ρn+1i−˜ρniΔt−κ˜ρni+1−2˜ρni+˜ρni−1Δx2=0,fori=1,…,N,

where is given in (2.5) and where the boundary values are

 ˜ρn0=2⟨vL−1(I−Π)(v+fℓ)⟩V⟨vL−1(vE)⟩V,˜ρnN+1=2⟨vL−1(I−Π)(v−fr)⟩V⟨vL−1(vE)⟩V,

for all . We recall that .

In the particular case where , , and , we obtain and

 ˜ρ0=∫10K2(v)fℓ(v)dv,˜ρN+1=∫0−1K2(−v)fr(v)dv.

with

 K2(v)=3v2. (3.11)

For instance, if , we get the value as boundary value for the limiting diffusion equation, whereas the right value can be computed from (2.9) and (2.10): . The value obtained by this scheme should clearly be improved. This is the goal of the next section. Note that the same boundary value 0.75 is obtained by the schemes proposed in [33] and [22], which both involve some artificial boundary conditions. We emphasize that, by construction, our scheme only involves the natural boundary conditions (3.5) and (3.6) derived from and and does not need any artificial condition.

Proof of Proposition 3.1.

We first check that all the values are uniquely defined by the scheme. The algorithm has been described above and it only remains to show that (3.7) allows to compute in terms of quantities at time . To this aim, we simply remark that the calculation of is done by solving a linear system of the form

 (ε2ΔtI−˜L)g=h,% with˜L=(I−ΠV−)L(I−ΠV−). (3.12)

Since the operator is self-adjoint and nonnegative, so is the operator and then the operator is invertible, for .

For the stability CFL condition, a similar proof as in [33, 35] can be done. We only sketch the main arguments of this proof and refer to these works for the details. Roughly, when , the convection term in the equation (3.7) prevails and induces the standard stability condition . When , the system behaves as the diffusion limiting equation (see below) and this results in the usual stability condition . In the general case, the CFL condition is expressed as a combination of these two limiting situations.

Now, we analyze formally the scheme when , the parameters , being fixed. First, we remark that the operator is invertible from to itself. More precisely, one has

 (I−ΠV−)Lg=h,withg∈{f:⟨f⟩V−=0},

if, and only if,

 g=(I−ΠV−)L−1(I−Π)h,withh∈{f:⟨f⟩V−=0}.

Now, by induction, it is clear from (3.7) that , for and , and then, from (3.9), we deduce that and for and . Reinserting these estimates in (3.7), we get

 gn+1i+1/2=ε(I−ΠV−)L−1(vE)(¯¯¯ρni+1−¯¯¯ρniΔx)+O(ε2) (3.13)

for , and

 gn+11/2=ε(I−ΠV−)L−1[vE(¯¯¯ρn1−¯¯¯ρ0Δx)−2Δx(I−Π)(v+gn0)]+O(ε2)
 gn+1N+1/2=ε(I−ΠV−)L−1[vE(¯¯¯ρN+1−¯¯¯ρnNΔx)+2Δx(I−Π)(v−gnN+1)]+O(ε2).

Inserting these equations into (3.9) gives, after some computations,

 ρn+1i−ρniΔt−κρni+1−2ρni+ρni−1Δx2=O(ε),% fori=2,…,N−1

and

 ρn+11−ρn1Δt−κρn2−2ρn1+~ρn0Δx2=O(ε),
 ρn+1N−ρnNΔt−κ~ρnN+1−2ρnN+ρnN−1Δx2=O(ε),

the quantities and being defined in Proposition 3.1. ∎

### 3.2. The numerical scheme

In this section, we construct a numerical scheme which is Asymptotic Preserving inside the domain and which provides a good approximation of the exact value at the boundary. To this purpose, we proceed as follows. Inside the space domain, we use the same numerical scheme as the one described in the previous section. At the boundary, we propose a specific treatment, without using neither artificial boundary condition nor any mesh refinement.

Our approach is based on the following idea. We consider the spatial derivative of the density at the boundary as an additional unknown in the scheme. It is known indeed that the density becomes stiff (its derivative is of order in the boundary layer) and the analysis of the boundary layer requires usually a rescaling process. Here, our strategy enables to avoid such rescaling. In the boundary cells, we simply determine the slope of by ensuring the balance between the incoming and outgoing half-fluxes of . We recall indeed that, if is a solution of the half-space problem (2.1), then one has and then, from (2.8), we deduce that and then, at , we have . We shall use this idea in the time evolution context.

Let us now describe our numerical scheme. We assume that we know and for , and for . We also assume that we know the approximation of the function at the boundary, i.e. and for all . We will now show how to compute all these quantities at time .

First step: determination of the quantities at time iteration , near the boundaries: , , and , ,

First, we use a centered scheme for the calculation of and :

 gn+11/2−gn1/2Δt + 1ε(I−ΠV−)(vgn1−gn0Δx) + 1ε(I−ΠV−)(vE)(¯¯¯ρn1−¯¯¯ρ0Δx)=1ε2(I−ΠV−)(Lgn+11/2),
 gn+1N+1/2−gnN+1/2Δt + 1ε(I−ΠV−)(vgnN+1−gnNΔx) + 1ε(I−ΠV−)(vE)(¯¯¯ρN+1−¯¯¯ρnNΔx)=1ε2(I−ΠV−)(Lgn+1N+1/2),

where and are determined by interpolation

 gn1=gn1/2+gn3/22,gnN=gnN−1/2+gnN+1/22.

It remains to determine and . For incoming velocities, these functions are prescribed thanks to (3.6):

 1V−(0)gn+10(v) = 1V−(0)(fℓ(v)−¯¯¯ρ0E), (3.16) 1V−(1)gn+1N+1(v) = 1V−(1)(fr(v)−¯¯¯ρN+1E), (3.17)

where is the characteristic function of a set . For outgoing velocities, we discretize the transport equation (2.20) in a suitable way. The specific point here is that we consider the derivative of at the boundary as an unknown. We approximate by the quantity and by the quantity , which are computed below. For only, we write

 gn+10−gn0Δt + λn+1ℓε2(I−ΠV−)(vE) (3.18) + 1ε(I−ΠV−)(vgn1/2−gn−1/2Δx)=1ε2(I−ΠV−)(Lgn+10)

and, for only, we write

 gn+1N+1−gnN+1Δt + λn+1rε2(I−ΠV−)(vE) (3.19) + 1ε(I−ΠV−)(vgnN+3/2−gnN−1/2Δx)=1ε2(I−ΠV−)(Lgn+1N+1).

In these equations, we have used the quantities and which are again defined by extrapolation:

 gn−1/2=2gn0−gn1/2andgnN+3/2=2gnN+1−gnN+1/2. (3.20)

Let us give a sense to these two equations (3.18) and (3.19) which are both written on . We only deal with the case of (3.18), the other equation can be treated similarly. This equation takes the form

 1V+(ε2ΔtI−˜L)gn+10=1V+h,with˜L=(I−ΠV−)L(I−ΠV−).

We recall that is already known on by (3.16). Hence, we have to solve the following linear system on the unknown :

 1V+(ε2ΔtI−˜L)1V+gn+10=1V+h+1V+˜L1V−gn+10,

where the right-hand side is given. To see that this system is well posed, it suffices to recall that the operator being self-adjoint and nonpositive, the bilinear form associated to is symmetric positive definite on the set of functions defined on . Hence this bilinear form is also symmetric positive definite on the set of functions defined on , and consequently the operator is invertible on this set: Eq. (3.18) admits a unique solution . From (3.18), we get

 1V+gn+10 = (1V+(ε2ΔtI−˜L)1V+)−1[(vgn1/2−gn−1/2Δ