IMplicit-EXplicit schemes for BGK kinetic equations Sandra - - PowerPoint PPT Presentation

IMplicit-EXplicit schemes for BGK kinetic equations

Sandra Pieraccini, Gabriella Puppo Politecnico di Torino Dipartimento di Matematica

HYP06, Lyon, July 17-21 – p.1/22

Boltzmann equation

Rarefied gas flow obeys Boltzmann equation which describes the evolution of the probability density f = f(x, v, t)

• f finding a particle at position x with microscopic velocity v at

time t: ∂tf + v · ∇xf = Q(f, f) where Q(f, f) is the collision integral which accounts for particle interactions.

HYP06, Lyon, July 17-21 – p.2/22

Collision integral

Particles interact through collisions in which a particle located at the point (x, v, t) in phase space collides with a particle at (x, v∗, t). The collision integral accounts for all possible interactions: Q(f, f) = σ

• ℜ3
• B−

(f(v′)f(v′

∗) − f(v)f(v∗)) |(v − v∗) · n|dS dv∗

During collisions, particles exchange momentum, so their location immediately after impact will be given by (x, v′, t) and (x, v′

∗, t). The relation between v, v∗ and v′, v′ ∗ contains the

physics of the interaction between particles.

HYP06, Lyon, July 17-21 – p.3/22

Adding chemistry...

The model described so far governs the motion of a monoatomic gas. To include more phenomena, the model becomes more complicated: For a polyatomic gas, the internal energy of the molecule must also be added as an independent variable. Modelling chemical reactions requires to consider distribution functions for all species involved and their reciprocal interactions in suitable collision integrals. Clearly the equations become more and more complicated, and their numerical solution extremely expensive. For this reason, simplified models for Boltzmann equation are introduced.

HYP06, Lyon, July 17-21 – p.4/22

Motivation for BGK model

The BGK model (Bhatnagar-Gross-Krook ’54) approximates Boltzmann equation for the evolution of a rarefied gas for small and moderate Knudsen numbers: Kn = mean free path characteristic length of the problem Lately, interest in this model has increased because: Several desirable properties have been shown to hold for the BGK model and its variants, such as BGK-ES, (Perthame et al. from 1989 on)

HYP06, Lyon, July 17-21 – p.5/22

Motivation for BGK model

The BGK model (Bhatnagar-Gross-Krook ’54) approximates Boltzmann equation for the evolution of a rarefied gas for small and moderate Knudsen numbers: Kn = mean free path characteristic length of the problem Lately, interest in this model has increased because: The BGK model has been extended to include more general fluids and can now be applied to the flow of a polytropic gas (Mieussens) and to mixtures of reacting gases (Aoki et al.)

HYP06, Lyon, July 17-21 – p.5/22

Motivation for BGK model

The BGK model (Bhatnagar-Gross-Krook ’54) approximates Boltzmann equation for the evolution of a rarefied gas for small and moderate Knudsen numbers: Kn = mean free path characteristic length of the problem Lately, interest in this model has increased because: New applications of kinetic models have appeared. For instance, fluid flow in nanostructures is described by Boltzmann equation, since it involves also moderate Knudsen numbers

HYP06, Lyon, July 17-21 – p.5/22

Numerical schemes for the BGK model

The development of numerical methods for the BGK model has started only recently. Yang, Huang ’95 This scheme is high order accurate in space, but only first

• rder accurate in time

HYP06, Lyon, July 17-21 – p.6/22

Numerical schemes for the BGK model

The development of numerical methods for the BGK model has started only recently. Aoki, Kanba, Takata ’97 This is a second order scheme, designed for smooth solutions

HYP06, Lyon, July 17-21 – p.6/22

Numerical schemes for the BGK model

The development of numerical methods for the BGK model has started only recently. Mieussens, ’00 A second order scheme, where conservation is exactly enforced

HYP06, Lyon, July 17-21 – p.6/22

Numerical schemes for the BGK model

The development of numerical methods for the BGK model has started only recently. Andries, Bourgat, le Tallec, Perthame ’02 A stochastic Monte Carlo scheme

HYP06, Lyon, July 17-21 – p.6/22

Numerical schemes for the BGK model

The development of numerical methods for the BGK model has started only recently. Pieraccini, Puppo ’06 A non oscillatory high order scheme in space and time

HYP06, Lyon, July 17-21 – p.6/22

BGK model

The main variable is the probability density f that a particle be in the point x ∈ Rd with velocity v ∈ RN at time t, thus f = f(x, v, t). The evolution of f is given by: ∂f ∂t + v · ▽xf = 1 τ (fM − f) , with initial condition f(x, v, 0) = f0(x, v) ≥ 0. Here τ is the collision time τ ≃ Kn, with τ > 0 and close to the hydrodynamic regime τ can be very small.

HYP06, Lyon, July 17-21 – p.7/22

The Maxwellian

fM is the local Maxwellian function, and it is built starting from the macroscopic moments of f: fM(x, v, t) = ρ(x, t) (2πRT(x, t))N/2 exp

• −||v − u(x, t)||2

2RT(x, t)

• ,

where ρ and u are the gas macroscopic density and velocity and T is the temperature. They are computed from f as:     ρ ρu E     =

• f

    1 v v2    

• where

g =

• RN g dv

E is total energy, and the temperature is: NT/2 = E − 1/2ρu2.

HYP06, Lyon, July 17-21 – p.8/22

The Maxwellian

Thus the BGK equation ∂f ∂t + v · ▽xf = 1 τ (fM − f) , describes the relaxation of f towards the local equilibrium Maxwellian fM. The local equilibrium is reached with a speed that is inversely proportional to τ. Thus the system is stiff for τ << 1.

HYP06, Lyon, July 17-21 – p.8/22

The Maxwellian

Note that the Maxwellian fM has the same moments of f, namely:     ρ ρu E     =

• f

    1 v v2    

• =
• fM

    1 v v2    

• HYP06, Lyon, July 17-21 – p.8/22

Conservation

As in Boltzmann equation, the macroscopic moments of f are conserved: ∂t f + ∇x · fv = 0, ∂t fv + ∇x · v ⊗ vf = 0, ∂t 1

2v2f

• + ∇x ·

1

2v2vf

• = 0.

Moreover, for τ → 0 the macroscopic solution converges to the gas dynamic solution of Euler equations.

HYP06, Lyon, July 17-21 – p.9/22

Conservation

As in Boltzmann equation, the macroscopic moments of f are conserved: ∂t f + ∇x · fv = 0, ∂t fv + ∇x · v ⊗ vf = 0, ∂t 1

2v2f

• + ∇x ·

1

2v2vf

• = 0.

Thus a numerical scheme for the BGK model must be con- servative, and its numerical solution must converge to the Eu- ler solution as τ → 0.

HYP06, Lyon, July 17-21 – p.9/22

Entropy principle

The BGK model satisfies an entropy principle: ∂t f log f + ∇x vf log f ≤ 0, ∀f ≥ 0 where equality holds if and only if f = fM. Thus the Maxwellian fM is the equilibrium solution of the system. The macroscopic entropy is: S = f log f Note that as τ → 0, entropy is conserved on smooth solutions, as for Euler solutions.

HYP06, Lyon, July 17-21 – p.10/22

Splitting schemes

To circumvent the stiffness of the source term, the BGK equation can be integrated through operator splitting.

HYP06, Lyon, July 17-21 – p.11/22

Splitting schemes

First, the convective part of the equations is integrated explicitly in time: ∂f ∂t + v · ▽xf = 0 In this fashion, space cells remain decoupled notwithstanding the term v · ▽x.

HYP06, Lyon, July 17-21 – p.11/22

Splitting schemes

First, the convective part of the equations is integrated explicitly in time: ∂f ∂t + v · ▽xf = 0 In this fashion, space cells remain decoupled notwithstanding the term v · ▽x. Next, the relaxation part of the equations is integrated implicitly in time, to circumvent stability restrictions on the time step induced by the stiff term: ∂f ∂t = 1 τ (fM − f) Note however that this procedure is only first order accurate in time

HYP06, Lyon, July 17-21 – p.11/22

IMEX RK schemes

To achieve the advantages of operator splitting, coupled with high order accuracy in time, we use IMplicit EXplicit (IMEX) Runge -Kutta schemes. To set notation, we consider the following autonomous ODE: y′ = f(y) + 1 εg(y) Let yn denote the numerical solution at time tn = n∆t. IMEX schemes are composed of an explicit and an implicit Runge- Kutta schemes, coupled together.

HYP06, Lyon, July 17-21 – p.12/22

IMEX RK schemes

The updated solution is computed as: yn+1 = yn + ∆t

ν

• i=1

˜ bif(y(i)) + ∆t ε

ν

• i=1

big(y(i)) where the stage values y(i) are given by: y(i) = yn + ∆t

i−1

• l=1

˜ ailf(y(l)) + ∆t ε

i

• l=1

ailg(y(l))

HYP06, Lyon, July 17-21 – p.12/22

IMEX RK schemes

The updated solution is computed as: yn+1 = yn + ∆t

ν

• i=1

˜ bif(y(i)) + ∆t ε

ν

• i=1

big(y(i)) Thus at each stage we group together all known terms: B(i) = yn + ∆t

i−1

• l=1

˜ ailf(y(l)) + ∆t ε

i−1

• l=1

ailg(y(l)) and then we solve implicitly the equation for y(i): y(i) = B(i) + ∆t ε aiig(y(i))

HYP06, Lyon, July 17-21 – p.12/22

An IMEX scheme for the BGK model

We propose a family of numerical schemes for the BGK model (Pieraccini, P ., ’06). The schemes we propose are based on a IMEX RK time advancement: The convective terms which are non stiff but highly non linear are integrated explicitly (thus avoiding coupling between neighboring cells) The source term, which can be very stiff for low Knudsen numbers, is integrated implicitly. We will show that the solution of the implicit step can be computed without solving non-linear equations.

HYP06, Lyon, July 17-21 – p.13/22

Time discretization

We start discretizing the equations in time. Let f n(x, v) = f(x, v, n∆t). Dropping the x and v dependency:

HYP06, Lyon, July 17-21 – p.14/22

Time discretization

f n+1 = f n − ∆t

ν

• i=1

˜ biv · ∇xf (i) + ∆t

ν

• i=1

bi τ (i)(f (i)

M − f (i)).

where f (i)

M (x, v), i = 1, ..., ν, are the Maxwellians computed

from the corresponding stage values f (i), which are: f (i) = f n − ∆t

i−1

• l=1

˜ ailv · ∇xf (l) + ∆t

i

• l=1

ail τ (l)(f (l)

M − f (l)).

˜ bi e ˜ ail define the explicit RK scheme for the convective terms, while bi and ail define the diagonally implicit RK scheme for the source term. Together, they compose the IMEX Runge-Kutta scheme for the time advancement of the solution.

HYP06, Lyon, July 17-21 – p.14/22

Examples

A second order IMEX scheme:

IMEX2 2nd order scheme, (Pareschi, Russo, JSC ’06)

1

1 2 1 2 1 2

− 1

2 1 2 1 2 1 2 1 2 1 2

HYP06, Lyon, July 17-21 – p.15/22

Examples

A third order IMEX scheme:

IMEX3 3rd order scheme, (Pareschi, Russo, JSC ’06)

1

1 4 1 4 1 6 1 6 2 3

α −α α 1 − α α β η

1 2 − β − η − α

α

1 6 1 6 2 3

with α = 0.24169426078821, β = 0.06042356519705, η = 0.1291528696059.

HYP06, Lyon, July 17-21 – p.15/22

Solution of the implicit step

At each stage i, the equation for f (i) can be written as: f (i) = B(i) + ∆t τ (i)aii

• f (i)

M − f (i)

where B(i) contains only previously computed terms. Note that f (i)

M depends non linearly on f (i) through the

moments of f. We want to show how it is possible to compute f without solv- ing non linear equations

HYP06, Lyon, July 17-21 – p.16/22

Solution of the implicit step

Let φ(v) be the vector φ(v) = (1, v, 1

2v2)T. We start

computing the moments of f (i):

• f (i)φ
• =
• B(i)φ
• + ∆t aii

τ (i)

• f (i)

M − f (i)

φ

• .

Since

• f (i)

M − f (i)

φ

• = 0

HYP06, Lyon, July 17-21 – p.16/22

Solution of the implicit step

• f (i)φ
• =
• B(i)φ
• .

and one immediately finds the macroscopic variables ρ(i), u(i) and T (i) corresponding to f (i). With these the updated Maxwellian is obtained: f (i)

M =

ρ(i) (2πRT (i))N/2 exp

• −||v − u(i)||2

2RT (i)

• With this knowledge, the equation for f (i) is linear and can be

solved exactly, for each x, v.

HYP06, Lyon, July 17-21 – p.16/22

Solution of the implicit step

In other words, the Maxwellian can be computed explicitely, and the collision term of the BGK scheme is linear in f, so that the solution is easily computed. (see also Coron Perthame, ’90). Note also that conservation is exact after time discretization.

HYP06, Lyon, July 17-21 – p.16/22

Space discretization

We choose a uniform grid in space, and discretize the equation with a conservative finite difference scheme, based

• n a high order non-oscillatory reconstruction. In particular

the convective term is discretized with upwinding for the first order scheme non oscillatory piecewise linear reconstruction and upwinding for the second order scheme WENO 5th order finite differences for the 3rd order in time scheme

HYP06, Lyon, July 17-21 – p.17/22

Convergence

Errors obtained with schemes of order 1, 2, and 3 on a smooth solution density

Scheme BGK1 Scheme BGK2 Scheme BGK3 L∞-errors Nx error

• rder

error

• rder

error

• rder

20 2.31e-2 1.17e-2 1.15e-2 40 1.24e-2 0.90 2.33e-3 2.33 1.69e-3 2.76 80 6.52e-3 0.93 5.02e-4 2.21 1.39e-4 3.61 160 3.33e-3 0.97 1.20e-4 2.06 5.54e-6 4.65 320 1.68e-3 0.99 2.89e-5 2.06 1.21e-7 5.52 640 8.42e-4 1.00 7.08e-6 2.03 8.72e-9 3.80

HYP06, Lyon, July 17-21 – p.18/22

Convergence

Errors obtained with schemes of order 1, 2, and 3 on a smooth solution Error on conservative variables

Scheme BGK1 Scheme BGK2 Scheme BGK3 Nx = NV = 21 ρ 6.060e-08 1.377e-06 5.680e-08 m 4.078e-10 2.748e-06 1.543e-07 E 9.815e-07 3.462e-06 1.154e-07 Nx = NV = 41 ρ 9.561e-13 5.941e-10 2.193e-10 m 1.365e-12 1.091e-09 3.902e-10 E 1.459e-12 1.106e-09 3.241e-10

HYP06, Lyon, July 17-21 – p.18/22

Smooth solution

The initial data has uniform density and temperature, with a localized perturbation in velocity Density

−0.2 0.2 0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 reference solution BGK1 BGK2 BGK3 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 1.01 1.012 1.014 1.016 1.018 1.02 1.022 1.024 1.026 1.028 1.03 reference solution BGK1 BGK2 BGK3

BGK1, BGK2, BGK3, reference solution.

HYP06, Lyon, July 17-21 – p.19/22

Smooth solution

The initial data has uniform density and temperature, with a localized perturbation in velocity Entropy decay in time

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 −2.8381 −2.838 −2.838 −2.8379 −2.8378 −2.8378 reference solution BGK1 BGK2 BGK3 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 −2.8381 −2.838 −2.838 −2.8379 −2.8378 −2.8378 reference solution BGK1 BGK2 BGK3

Left: τ = 10−2. Right: τ = 10−5. BGK1, BGK2, BGK3, reference solution.

HYP06, Lyon, July 17-21 – p.19/22

Riemann problem

The initial data has an initial step in density and pressure Density

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 reference solution BGK1 BGK2 BGK3 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 reference solution BGK1 BGK2 BGK3

Left: τ = 10−2, Right:τ = 10−5 BGK1, BGK2, BGK3 Reference solution.

HYP06, Lyon, July 17-21 – p.20/22

Riemann problem

The initial data has an initial step in density and pressure Entropy decay

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 −1.1 −1.09 −1.08 −1.07 −1.06 −1.05 −1.04 reference solution BGK1 BGK2 BGK3 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 −1.1 −1.09 −1.08 −1.07 −1.06 −1.05 −1.04 reference solution BGK1 BGK2 BGK3

Left: τ = 10−2, Right:τ = 10−5 BGK1, BGK2, BGK3 Reference solution. Note that entropy decays even for small τ due to the presence

• f a shock wave.

HYP06, Lyon, July 17-21 – p.20/22

Behavior of f

This figure shows the behavior of the probability density f for various values of x in the Riemann problem.

−10 −5 5 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x=0.1 x=0.3 x=0.4 x=0.5 x=0.75 x=0.85

HYP06, Lyon, July 17-21 – p.21/22

Behavior of f

Comparison between f and the corresponding local equilibrium Maxwellian: τ = 1e − 2.

−4 −3 −2 −1 1 2 3 4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 rarefaction contact discontinuity shock HYP06, Lyon, July 17-21 – p.21/22

Behavior of f

Comparison between f and the corresponding local equilibrium Maxwellian: τ = 1e − 5.

−4 −3 −2 −1 1 2 3 4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 rarefaction contact discontinuity shock HYP06, Lyon, July 17-21 – p.21/22

.

Thank you!

HYP06, Lyon, July 17-21 – p.22/22