How to avoid very large discrete velocity grids in deterministic - - PowerPoint PPT Presentation

How to avoid very large discrete velocity grids in deterministic simulations of rarefied gas flows?

Luc Mieussens

Institut de Mathématiques de Bordeaux, Université de Bordeaux, France

Issues in Solving the Boltzmann Equation for Aerospace Applications (June 3-7, 2013)

collaborators: C. Baranger - J. Claudel - N. Hérouard (CEA)

• S. Brull - L. Forestier-Costes (IMB)

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 1 / 49

Introduction

flows of rarefied gases Deterministic simulations based on a discrete velocity approximation

• f the Boltzmann (BGK) equation (discrete ordinate method)

A global discrete velocity grid, generally Cartesian, is commonly used for the whole computational domain (Kyoto group, Aristov et al., etc.) For practical applications in aerodynamics (atmospheric re-entry problems), the grid is so large that the computational ressources (memory storage and CPU time) required by the simulation are huge two projects: automatic generation of a locally refined velocity grid (steady flows), and construction of local velocity grids (unsteady flows)

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 2 / 49

Outline of the talk

1

Basics for discrete velocity approximations of the models of RGD

2

Steady flows: local refinement of the velocity grid

3

Unsteady flows: a method to use local velocity grids

4

Conclusions

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 3 / 49

Outline

1

Basics for discrete velocity approximations of the models of RGD

2

Steady flows: local refinement of the velocity grid

3

Unsteady flows: a method to use local velocity grids

4

Conclusions

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 4 / 49

Discrete velocity approximation of a kinetic equation of RGD

Physical model : steady equation v · ∇xf = Q(f ) + BCs Velocity space R3 V =

• vk, k ∈ N3

discrete velocity set Discrete kinetic equation: vk · ∇xfk = Qk(f ) + BCs

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 5 / 49

How to design a discrete velocity set?

Generally: V is a Cartesian grid

V =   vl,m,n = c +   l∆vx m∆vy n∆vz   ,   l m n   =   0 : lmax 0 : mmax 0 : nmax     

Constraints: bounds and grid step : the grid must be

large enough to contain the main part of the distribution functions for every position in the computational domain fine enough to capture the core of the distribution functions for every position in the computational domain

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 6 / 49

How to design a Cartesian velocity grid?

At each space position x, we have information on the support of the distribution function f (x, .): assume f is close to its local Maxwellian, then it is centered on u(x) with standard deviation

• RT(x)

u(x) − c

• RT(x)

u(x) u(x) + c

• RT(x)
• RT(x)

The grid contains all the distributions if the bounds are at least vmin = min

x∈Ω{u(x) − c

• RT(x)}

vmax = max

x∈Ω{u(x) + c

• RT(x)}

(c around 4)

There should be at least 3 points into the “core” of each distribution. The grid step should be: ∆v ≤ min

x∈Ω{

• RT(x)}

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 7 / 49

Estimation of the bounds and the step of the grid

Temperature and velocity fields are a priori unknown! Several solutions :

Use the boundary values: upstream flow (u∞, T∞) and velocity and temperature of the body (uwall, Twall) More information with flow parameters inside the shock: ushock and Tshock can be estimated by Rankine-Hugoniot relations. Even better: a compressible Navier-Stokes pre-simulation:

gives uNS and T NS in Ω set vmin = min

x∈Ω{uNS(x) − c

• RT NS(x)}

vmax = max

x∈Ω{uNS(x) + c

• RT NS(x)}

∆v = a min

x∈Ω{

• RT NS(x)}

simple, fast and efficient

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 8 / 49

Drawbacks of a Cartesian grid

For re-entry problems: the velocity is very large (≈ 6.000 m/s), the temperature is very large in the shock (∼ 105K), and very small in the upstream and at the boundary (∼ 102K) Consequently: very large grid bounds, and very small step, so very large number of discrete velocities Example: Mach 20, altitude 90 km: V contains 52 × 41 × 41 points. Around 350 GB memory requirements with a coarse 3D mesh in space!...

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 9 / 49

What can be done?

existing works

for unsteady flows: variable velocity grids (all the grids have the same bounds, and the grids are refined or coarsened by AMR techniques), see Kolobov et al. (2011) and K. Xu et al. (2012) for steady flows: use the unsteady solution, but expensive

this talk:

for steady flows: estimation of an “optimal” locally refined global grid for unsteady flows: variable velocity grids, with variable bounds and steps

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 10 / 49

Outline

1

Basics for discrete velocity approximations of the models of RGD

2

Steady flows: local refinement of the velocity grid

3

Unsteady flows: a method to use local velocity grids

4

Conclusions

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 11 / 49

Idea in 1D

1D view: several distributions to be represented on the same grid Natural idea: refine the grid in the main parts of the distributions: Main contribution of this work: automatic generation of a locally refined velocity grid based on a CNS pre-simulation

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 12 / 49

Automatic generation of a locally refined velocity grid

Main “ingredient”: construction of a “support” function

Assume we have an estimation of the macroscopic velocity and temperature fields (compressible Navier-Stokes simulation) to build a Cartesian velocity grid V At each point x ∈ Ω, the “main” support of the distribution function f (x, .) in the velocity space is centered at u(x) with a radius 4

• RT(x)

We define the support function φ(v) by

φ(v) = min

• RT(x), where x is such that ||v − u(x)|| ≤ 4
• RT(x)
• .

Take a velocity v: among all the distribution functions whose support contains v, one of them has a smallest support, and φ(v) is the standard deviation of this distribution. φ computed on the thin Cartesian grid defined by u and T.

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 13 / 49

Automatic generation of a l. r. velocity grid: algorithm

Example : flow around a cylinder Tells us how the Cartesian grid around vl,m,n should be refined, or coarsened Algorithm to generate an AMR velocity grid:

recursive algorithm that cuts a cell if it is larger than the minimum of φ in the velocity cell property of the final grid: every cell has a size smaller than the minimum of the support function in the cell.

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 14 / 49

Automatic generation of a l. r. v. grid: application

Computational code

2D steady rarefied flows, monoatomic gases, BGK model finite volume scheme on curvilinear structured grids entropic and conservative velocity discretization

Why BGK?

simpler than Boltzmann (parametric model)

• ften sufficiently accurate in transition regime and slow microflows

(see "A Unified Gas-Kinetic Scheme for Continuum and Rarefied Flows II: Multi-Dimensional Cases" Hung et al. Commun. Comput. Phys. 12 (3) 2012)

more physics might be included (polyatomic, multispecies, reactive)

Code of the CEA

3D flows polyatomic gases

2D test case: steady flow over a cylinder of radius 0.1m, Mach 20, altitude 90

km, space mesh of 50 × 50 cells

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 15 / 49

Automatic generation of a l. r. v. grid: application

Compressible Navier-Stokes pre-simulation gives a uniform Cartesian grid of 44 × 45 = 1980 points, total memory required: 1 Gb Support function and induced AMR velocity grid:

AMR grid of 316 points, memory required: 140 Mb

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 16 / 49

Automatic generation of a l. r. v. grid: application

Gain in memory: factor 7 Gain in CPU time: factor 7 Accuracy:

average error less than 2% for macroscopic quantities maximum error on the normal heat flux along the boundary smaller than 3%

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 17 / 49

Automatic generation of a l. r. v. grid: application - 2D axisymmetric

air flow around a sphere (radius=1cm), Ma=20, altitude=90km gain: factor 10 (CPU and memory) accuracy: less than 1% error

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 18 / 49

Automatic generation of a l. r. v. grid: application - 2D axisymmetric

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 19 / 49

Automatic generation of a l. r. v. grid: application - full 3D

air flow, cone with spherical nose, (radius= 1 cm), Ma= 20, altitude= 90km gain: factor 27 (CPU and memory)

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 20 / 49

Automatic generation of a l. r. v. grid: application - full 3D

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 21 / 49

Modifications to use our l. r. v. grid on an existing discrete velocity code

for a code that uses a Cartesian grid, the 3D array-data structure must be replaced by a 1D array (or octree) the quadrature formulae (moments of the distribution function) must be modified (new weights): here we use quadrature obtained through constant or bilinear interpolation collision operator: the discretization must be modified (easy for BGK)

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 22 / 49

Outline

1

Basics for discrete velocity approximations of the models of RGD

2

Steady flows: local refinement of the velocity grid

3

Unsteady flows: a method to use local velocity grids

4

Conclusions

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 23 / 49

Needs for local grids

previous approach: still a global velocity grid if the bounds are large, we still need a lot of points even difficult to know a priori the correct size of the grid : how to estimate the bounds for an unsteady flow?

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 24 / 49

Needs for local grids: example

“two interacting blast waves”: at t = 0 ρ = 1, u = 0, T = [1000, 0.01, 100] max (u + 4 √ RT) = 126 at t = 0.02 (after interaction): max (u + 4 √ RT) = 236: larger bounds! the optimal global velocity grid has 2551 points!

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 25 / 49

Local velocity grids

idea: define a different velocity grid for every t and x no global bounds (different from Kolobov et al.) different from the strategy of rescaled velocity v → (v − u)/ √ RT

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 26 / 49

Local velocity grids: problems

1 how to define the local velocity grids (LVG): bounds and step? 2 how to exchange information between two grids? Example:

∂xf (t, x, v) ≈ f (t, x + ∆x, v) − f (t, x, v) ∆x but f (t, x, .) and f (t, x + ∆x, .) are not defined on the same velocity grid... this work: preliminar results (in 1D and 2D)

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 27 / 49

Definition of a LVG: conservation laws

at tn: space cell [xi− 1

2 , xi+ 1 2 ]

f n

i,k approximates f (tn, xi, vn i,k), where

vn

i,k ∈ Vn i =

• vn

i,1, vn i,2, . . . , vn i,K n

i

• local velocity grid (with K n

i points)

the discrete moment vector Un

i approximates U(tn, xi) with a

quadrature on the LVG:

U(tn, xi) = (ρ, ρu, E)(tn, xi) =

• R

m(v)f (tn, xi, v) dv = mf (tn, xi, .) ↓ Un

i = (ρn i , ρn i un i , E n i ) = K n

i

• k=1

m(v n

i,k)f n i,kωn i,k = mf n i Vn

i

where m(v) = (1, v, 1

2|v|2). Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 28 / 49

Definition of a LVG: conservation laws

approximation of U(tn+1, xi) by the discrete conservation laws: conservation laws ∂tU + ∂xvmf = 0 finite volume upwind scheme Un+1

i

− Un

i

∆t + Φn

i+ 1

2 − Φn

i− 1

2

∆x = 0 where the numerical fluxes are defined by Φn

i+ 1

2 =

• v+mf n

i

• Vn

i +

• v−mf n

i+1

• Vn

i+1 Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 29 / 49

Definition of a LVG

then we have approximations of u(tn+1, xi) and T(tn+1, xi): un+1

i

and T n+1

i

we set the bounds of Vn+1

i

to vn+1

min,i = un+1 i

− 4

• RT n+1

i

and vn+1

max,i = un+1 i

+ 4

• RT n+1

i

remark: these bounds are correct if f is not too far from its local Maxwellian new dicrete velocity grid Vn+1

i

: so far, we define uniform grids with a constant number of points (K n+1

i

= 10 to 30) remark: variable step and number of points could be necessary (discontinuous f ).

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 30 / 49

Computation of f n+1

i

kinetic equation: ∂tf + v∂xf = Q(f ) finite volume - upwind scheme (continuous v):

f n+1

i

(v) − f n

i (v)

∆t + v + f n

i (v) − f n i−1(v)

∆x + v − f n

i+1(v) − f n i (v)

∆x = Q(f n

i )(v)

problem: in our discrete velocity method, f n+1

i

, f n

i , f n i±1 are not

defined on the same grids.

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 31 / 49

Computation of f n+1

i

reconstruction procedure: for each i, f n

i

is reconstructed on its LVG grid Vn

i

¯ f n

i (v) =

• R(f n

i )(v)

if v n

min,i ≤ v ≤ v n max,i

else,

¯ f n

i (v)

f n

i

vn

i,k

vn

min,i

vn

max,i

numerical scheme: f n+1

i

(v) − f n

i (v)

∆t + v + f n

i (v) − f n i−1(v)

∆x + v − f n

i+1(v) − f n i (v)

∆x = Q(f n

i )(v)

written for every v n+1

i,k

• f Vn+1

i

to get: f n+1

i,k

− ¯ f n

i (v n+1 i,k )

∆t + v n+1

i,k + ¯

f n

i (v n+1 i,k ) − ¯

f n

i−1(v n+1 i,k )

∆x + v n+1

i,k − ¯

f n

i+1(v n+1 i,k ) − ¯

f n

i (v n+1 i,k )

∆x = ¯ Q(f n

i )(v n+1 i,k ) Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 32 / 49

Advantages of this scheme

1 advantage 1: the local velocity grids adapt in time and space to the

local temperature T and velocity u

2 advantage 2: only initial values for u and T are required Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 33 / 49

Properties of the numerical scheme

1 if R is positive, then the scheme is positive, provided that

∆t satisfies a standard CFL condition T n+1

i

is positive

2 under this slight modification on the quadrature:

f n

i Vn

i := ¯

f n

i

T n+1

i

can be proved to be positive

3 in practice: even with non positive R, the scheme generally preserves

the positivity

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 34 / 49

Illustration 1: Sod test case

velocity profile:

0,1 0,2 0,3 0,4 0,5 0,1 0,2 0,3 0,4 0,5 0,6

... : global velocity grid (wrong bounds, with initial data)

• : global velocity grid (correct bounds, 100 points: bounds and steps estimated

by an Euler computation, then a grid convergence study) .- : local velocity grid (10 points)

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 35 / 49

Illustration 2: two interacting blast waves

before the interaction (top), after (bottom)

T u p

0,2 0,4 0,6 0,8 1 0,5 1 1,5 2 2,5 3 0,2 0,4 0,6 0,8 1

• 5

5 10 15 20 0,2 0,4 0,6 0,8 1 100 200 300 400 0,2 0,4 0,6 0,8 1 1 2 3 4 0,2 0,4 0,6 0,8 1 8 10 12 14 16 18 20 0,2 0,4 0,6 0,8 1 300 350 400 450 500 550 600

global velocity grid: 2531 points (bounds and steps estimated by an Euler computation, then grid convergence study) local velocity grid: 30 points

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 36 / 49

Reconstruction procedure

we use piecewise polynomial interpolation high order interpolation is necessary (while very low order interpolation is sufficient for quadratures) in case of discontinuous distributions (large Knudsen numbers), a non

• scillatory method is necessary (ENO interpolation used)

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 37 / 49

Influence of the reconstruction step: Sod test

ρ u T

• 1
• 0,5

0,5 1 0,2 0,4 0,6 0,8 1

• 1
• 0,5

0,5 1 0,2 0,4 0,6 0,8

• 1
• 0,5

0,5 1 0,6 0,8 1 1,2 1,4 1,6 1,8 2

• 1
• 0,5

0,5 1 0,2 0,4 0,6 0,8 1

• 1
• 0,5

0,5 1 0,2 0,4 0,6 0,8 1

• 1
• 0,5

0,5 1 0,6 0,8 1 1,2 1,4 1,6 1,8 2 2,2 2,4

Free transport regime (no collisions): top: global grid (30 points) bottom: local velocity grid (30 points): linear interpolation, ENO3, ENO4

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 38 / 49

A 2d view of the LVG: Sod test case, free transport

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 39 / 49

Boundary conditions

diffuse reflection: f (t, x = 0, v > 0) = − v−f

v+MwMw

ghost cell approach Mw is approximated on its own local velocity grid (based on the wall temperature and velocity) however, such reflection can generate very discontinuous and non-symmetric distributions!

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 40 / 49

Boundary conditions

example: heat transfer problem for small t, the support of f (t, x ≈ 1, .) is non symmetric:

TR TL t = 0   ρ0 u0 = 0 T0 = TL   x = 0 x = 1

• RTL
• RTR

the local velocity grid (symmetric by construction) is not large enough for v < 0 (−1900 instead of −2200).

0,1 0,2 0,3 0,4 0,5 0,6

• 50
• 40
• 30
• 20
• 10

(Kn = 0.01, global velocity grid=100 points, LVG=30 points)

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 41 / 49

Extension of the LVG: algorithm

compute f n+1

i,k

for every vn+1

i,k

• f Vn+1

i

set w = vn+1

i,1

(leftmost point) loop

w = w − ∆v n+1

i

compute f n+1

i

(w) by the numerical scheme:

f n+1

i

(w) − ¯ f n

i (w)

∆t +w + ¯ f n

i (w) − ¯

f n

i−1(w)

∆x +w − ¯ f n

i+1(w) − ¯

f n

i (w)

∆x = ¯ Q(f n

i )(w)

if f n+1

i

(w) is too large then add w to the grid and continue the loop else stop

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 42 / 49

Extension of the LVG: illustration

heat transfer problem: T u p

0,1 0,2 0,3 0,4 0,5 0,6 300 400 500 600 700 800 900 1000 0,1 0,2 0,3 0,4 0,5 0,6

• 50
• 40
• 30
• 20
• 10

0,1 0,2 0,3 0,4 0,5 0,6 1 1,5 2 2,5 3 0,1 0,2 0,3 0,4 0,5 0,6 300 400 500 600 700 800 900 0,1 0,2 0,3 0,4 0,5 0,6

• 50
• 40
• 30
• 20
• 10

0,1 0,2 0,3 0,4 0,5 1 1,5 2 2,5 3

no extension, symmetric LVG (top) with extension (bottom) Kn = 0.01, global velocity grid=100 points, LVG=30 points

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 43 / 49

2D LVG

Heat transfer problem between two cylinders

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 44 / 49

2D LVG

Heat transfer problem between two cylinders

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 45 / 49

2D LVG: problems

systematic interpolation: for each time step, in every space cell, for every velocity very large computational cost slower than the Global Grid DVM for some cases!

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 46 / 49

2D LVG: developpements

embedded local grids: two neighbouring space cells i and j have local grids Vi and Vj such that:

∆vi = k∆vj, k integer their bounds might be shifted of l∆vi

• 4
• 2

2 4 50 100 150 200 250 300

very few interpolation= much faster still performance issues: slower memory acces? tests in progress

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 47 / 49

2D LVG: developpements

embedded local grids: two neighbouring space cells i and j have local grids Vi and Vj such that:

∆vi = k∆vj, k integer their bounds might be shifted of l∆vi

• 4
• 2

2 4 50 100 150 200 250 300

very few interpolation= much faster still performance issues: slower memory acces? tests in progress

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 47 / 49

2D LVG: developpements

embedded local grids: two neighbouring space cells i and j have local grids Vi and Vj such that:

∆vi = k∆vj, k integer their bounds might be shifted of l∆vi

• 4
• 2

2 4 50 100 150 200 250 300

very few interpolation= much faster still performance issues: slower memory acces? tests in progress

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 47 / 49

2D LVG: developpements

embedded local grids: two neighbouring space cells i and j have local grids Vi and Vj such that:

∆vi = k∆vj, k integer their bounds might be shifted of l∆vi

• 4
• 2

2 4 50 100 150 200 250 300

very few interpolation= much faster still performance issues: slower memory acces? tests in progress

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 47 / 49

2D LVG: developpements

embedded local grids: two neighbouring space cells i and j have local grids Vi and Vj such that:

∆vi = k∆vj, k integer their bounds might be shifted of l∆vi

• 4
• 2

2 4 50 100 150 200 250 300

very few interpolation= much faster still performance issues: slower memory acces? tests in progress

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 47 / 49

Outline

1

Basics for discrete velocity approximations of the models of RGD

2

Steady flows: local refinement of the velocity grid

3

Unsteady flows: a method to use local velocity grids

4

Conclusions

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 48 / 49

Conclusions

two different methods (steady/unsteady) encouraging results efficient methods for problems with large variations of temperature (re-entry, shock interactions) local refinement of a global grid:

simple and efficient algorithm to generate a locally refined velocity grid for discrete velocity approximations of the BGK equation important gains both on memory storage and CPU time validated on 2D cases (both for plane and axisymmetric configurations) for re-entry problems perspectives: 3D validation (in progress, CEA code)

local velocity grids:

still some problems perspectives: more tests, better performance, analysis

Luc Mieussens (Bordeaux, France) Variable velocity grids for DVM ICERM 2013 49 / 49