SLIDE 1 Staggered schemes for compressible flows
with
et⋆, L. Gastaldo†, W. Kheriji⋆†, J.-C. Latch´ e†, T.T. Nguyen⋆†
⋆ Universit´ e de Provence † Institut de Radioprotection et de Sˆ uret´ e Nucl´ eaire (IRSN)
SLIDE 2 Context
brèche circuit primaire cuve coeur brèche cuve corium béton 6 m 1 à 4 m béton phase liquide mouvement de convection
⊲ General context: nuclear safety Examples: Accident in a nuclear reactor, fol- lowing loss of coolant. Interaction corium–concrete → two phase flow Numerical simulation of compressible flows – numerical code: ISIS, developed at IRSN
SLIDE 3 Context
brèche circuit primaire cuve coeur brèche cuve corium béton 6 m 1 à 4 m béton phase liquide mouvement de convection
⊲ General context: nuclear safety Examples: Accident in a nuclear reactor, fol- lowing loss of coolant. Interaction corium–concrete → two phase flow Numerical simulation of compressible flows – numerical code: ISIS, developed at IRSN ⊲ Aim : obtain efficient schemes
◮ stable and precise for all Mach number ◮ sufficiently decoupled so that the
numerical solution is not too difficult. fractional time step methods. Classical for incompressible flows (pressure correction, Chorin 68, Temam 69), (Guermond 06) for a synthesis. Also developed for compressible flows, with either colocated or staggered unknowns..
SLIDE 4 Context
brèche circuit primaire cuve coeur brèche cuve corium béton 6 m 1 à 4 m béton phase liquide mouvement de convection
⊲ General context: nuclear safety Examples: Accident in a nuclear reactor, fol- lowing loss of coolant. Interaction corium–concrete → two phase flow Numerical simulation of compressible flows – numerical code: ISIS, developed at IRSN ⊲ Aim : obtain efficient schemes
◮ stable and precise for all Mach number ◮ sufficiently decoupled so that the
numerical solution is not too difficult. fractional time step methods. Classical for incompressible flows (pressure correction, Chorin 68, Temam 69), (Guermond 06) for a synthesis. Also developed for compressible flows, with either colocated or staggered unknowns.. ⊲ Theoretical proof of stability and consis- tency, confirmed by numerical tests.
SLIDE 5 The drift-diffusion model
Drift-diffusion model for two phase flows
∂t + ∇ · (ρ u) = 0 (mass) ∂ ρ y ∂t + ∇ · (ρ y u) = −∇ · (ρ y (1 − y) ur ) + ∇ · (D∇y) (ffl) ∂ρ u ∂t + ∇ · (ρ u ⊗ u) + ∇p − ∇ · τ(u) = 0 (mom) p = Ψ(ρ, y) Physical properties
◮ Positivity of the density ◮ Gas mass fraction between 0 and 1 ◮ Total mass conservation and fractional mass conservation ◮ Transport of an interface between phases with constant pressure and velocity ◮ Stability: control on the entropy.
SLIDE 6 Staggered discretization
◮ Positivity of density finite volume scheme on ρ:
ρ piecewise constant on the cells.
SLIDE 7 Staggered discretization
◮ Positivity of density finite volume scheme on ρ:
ρ piecewise constant on the cells. Primal mesh : M = { set of control volumes K, L, M... }.
SLIDE 8 Staggered discretization
◮ Positivity of density finite volume scheme on ρ:
ρ piecewise constant on the cells. Primal mesh : M = { set of control volumes K, L, M... }. Scalar variables defined at cell centers: (pK )K∈M, (̺K )K∈M,. . .
◮ Velocity components defined at the (or some of the) edges : (v (i)
σ )σ∈F(i).
SLIDE 9 Staggered discretization
◮ Positivity of density finite volume scheme on ρ:
ρ piecewise constant on the cells. Primal mesh : M = { set of control volumes K, L, M... }. Scalar variables defined at cell centers: (pK )K∈M, (̺K )K∈M,. . .
◮ Velocity components defined at the (or some of the) edges : (v (i)
σ )σ∈F(i).
Dual mesh(es) : Mn = (D(i)
σ )σ∈F(i).
SLIDE 10 Staggered discretization
◮ Positivity of density finite volume scheme on ρ:
ρ piecewise constant on the cells. Primal mesh : M = { set of control volumes K, L, M... }. Scalar variables defined at cell centers: (pK )K∈M, (̺K )K∈M,. . .
◮ Velocity components defined at the (or some of the) edges : (v (i)
σ )σ∈F(i).
Dual mesh(es) : Mn = (D(i)
σ )σ∈F(i).
Normal velocity to the face σ denoted by vσ · nσ. Dσ Dσ′ σ
′
= K|M K L M | σ | σ = K | L ǫ = Dσ|Dσ′
Figure: Primal and dual meshes for the Rannacher-Turek and Crouzeix-Raviart elements.
SLIDE 11 The MAC mesh
: Dσ (or K x
i− 1
2 ,j)
xi− 3
2
xi− 1
2
xi+ 1
2
yj− 3
2
yj− 1
2
yj+ 1
2
yj+ 3
2
ux
i− 1
2 ,j
ux
i− 3
2 ,j
ux
i+ 1
2 ,j
ux
i− 1
2 ,j−1
ux
i− 1
2 ,j+1
: Dσ (or K y
i,j− 1
2
) xi− 3
2
xi− 1
2
xi+ 1
2
xi+ 3
2
yj− 3
2
yj− 1
2
yj+ 1
2
uy
i,j− 1
2
uy
i−1,j− 1
2
uy
i+1,j− 1
2
uy
i,j− 3
2
uy
i,j+ 1
2
The dual mesh for the MAC scheme, x and y-component of the velocity.
SLIDE 12
Finite volume discretization of the mass equation
∂tρ + div(ρ u) = 0, (mass)
SLIDE 13 Finite volume discretization of the mass equation
∂tρ + div(ρ u) = 0, (mass)
◮
(mass) + implicit time discretization
ρn+1 − ρn δt +
(ρn+1 un+1 · nK) = 0.
SLIDE 14 Finite volume discretization of the mass equation
∂tρ + div(ρ u) = 0, (mass)
◮
(mass) + implicit time discretization
ρn+1 − ρn δt +
(ρn+1 un+1 · nK) = 0.
◮ discretization of the fluxes:
|K| δt (ρn+1
K
− ρn
K ) +
F n+1
K,σ = 0,
◮ F n+1
K,σ = |σ| ˇ
ρn+1
σ
un+1
σ
· nK,σ , numerical flux through σ.
◮ ˇ
ρn+1
σ
upwind approximation of ρn+1 at the face σ with respect to un+1
σ
· nK,σ.
SLIDE 15 Finite volume discretization of the mass equation
∂tρ + div(ρ u) = 0, (mass)
◮
(mass) + implicit time discretization
ρn+1 − ρn δt +
(ρn+1 un+1 · nK) = 0.
◮ discretization of the fluxes:
|K| δt (ρn+1
K
− ρn
K ) +
F n+1
K,σ = 0,
◮ F n+1
K,σ = |σ| ˇ
ρn+1
σ
un+1
σ
· nK,σ , numerical flux through σ.
◮ ˇ
ρn+1
σ
upwind approximation of ρn+1 at the face σ with respect to un+1
σ
· nK,σ.
◮ Positive density: ρn+1 > 0 if (ρn > 0 and ρ > 0 at inflow boundary).
SLIDE 16
Discretization of the phase mass fraction balance
∂tρy + div(ρ uy) = −∇ · (ρ y (1 − y) ur ) + ∇ · (D∇y), (ff)
SLIDE 17 Discretization of the phase mass fraction balance
∂tρy + div(ρ uy) = −∇ · (ρ y (1 − y) ur ) + ∇ · (D∇y), (ff)
◮
(mass) + implicit time discretization
ρn+1y n+1 − ρny n δt +
(ρn+1y n+1 un+1 · nK ) = −
ρ y (1 − y) ur · nK ) +
D∇y · nK .
SLIDE 18 Discretization of the phase mass fraction balance
∂tρy + div(ρ uy) = −∇ · (ρ y (1 − y) ur ) + ∇ · (D∇y), (ff)
◮
(mass) + implicit time discretization
ρn+1y n+1 − ρny n δt +
(ρn+1y n+1 un+1 · nK ) = −
ρ y (1 − y) ur · nK ) +
D∇y · nK .
◮ discretization of the left hand side: |K|
δt (ρn+1
K
y n+1
K
− ρn
Ky n K ) +
Fn+1
K,σy n+1 σ
y n+1
σ
upwind approximation of y n+1 at the face σ with respect to F n+1
K,σ .
◮ monotone numerical flux for
◮ two point flux approximation for
SLIDE 19 Discretization of the phase mass fraction balance
∂tρy + div(ρ uy) = −∇ · (ρ y (1 − y) ur ) + ∇ · (D∇y), (ff)
◮
(mass) + implicit time discretization
ρn+1y n+1 − ρny n δt +
(ρn+1y n+1 un+1 · nK ) = −
ρ y (1 − y) ur · nK ) +
D∇y · nK .
◮ discretization of the left hand side: |K|
δt (ρn+1
K
y n+1
K
− ρn
Ky n K ) +
Fn+1
K,σy n+1 σ
y n+1
σ
upwind approximation of y n+1 at the face σ with respect to F n+1
K,σ .
◮ monotone numerical flux for
◮ two point flux approximation for
: 0 < yn+1 < 1 if (0 < y n < 1). Existence and uniqueness of the solution y n+1 for a given y n (Gastaldo-H.-Latch´ e 2009)
SLIDE 20
FV-FE discretization of the momentum equation
∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u)) = 0
SLIDE 21 FV-FE discretization of the momentum equation
∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u)) = 0
◮
(momentum) + implicit time discretization
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nK ) +
(∇pn+1 − div(τ(un+1))) = 0.
SLIDE 22 FV-FE discretization of the momentum equation
∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u)) = 0
◮
(momentum) + implicit time discretization
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nK ) +
(∇pn+1 − div(τ(un+1))) = 0.
◮ Space discretization
|Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
+|Dσ|(∇pn+1)σ − |Dσ|(divτ(un+1))σ = 0,
◮ |Dσ|(∇pn+1)σ,i = −
pn+1 divϕ(i)
σ dx = |σ| (pn+1 L
− pn+1
K
) nK,σ · e(i),
◮ |Dσ|(divτ(un+1))σ,i = −µ
∇un+1 · ∇ϕ(i)
σ − µ
3
div un+1 div ϕ(i)
σ .
◮ ϕ(i)
σ i-th finite element shape function.
SLIDE 23 Discretization of the convection operator
◮ Choice of ρn+1
σ
and F n+1
σ,ǫ in
|Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
?
SLIDE 24 Discretization of the convection operator
◮ Choice of ρn+1
σ
and F n+1
σ,ǫ in
|Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
?
◮ Discretize
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nK ) so as to obtain a discrete kinetic energy balance.
SLIDE 25 Discretization of the convection operator
◮ Choice of ρn+1
σ
and F n+1
σ,ǫ in
|Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
?
◮ Discretize
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nK ) so as to obtain a discrete kinetic energy balance.
◮ Copy the continuous kinetic energy balance:
- ∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u))
=0
2 ρ |u|2) + div
2 ρ |u|2) u
=0, with some formal algebra... using ∂tρ + div(ρ u) = 0.
SLIDE 26 Discretization of the convection operator
◮ Choice of ρn+1
σ
and F n+1
σ,ǫ in
|Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
?
◮ Discretize
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nK ) so as to obtain a discrete kinetic energy balance.
◮ Copy the continuous kinetic energy balance:
- ∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u))
=0
2 ρ |u|2) + div
2 ρ |u|2) u
=0, with some formal algebra... using ∂tρ + div(ρ u) = 0.
◮ Do the same at the discrete problem ?
♭ Momentum on dual cells, mass on primal cells...
SLIDE 27 Discretization of the convection operator
◮ Choice of ρn+1
σ
and F n+1
σ,ǫ in
|Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
?
◮ Discretize
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nK ) so as to obtain a discrete kinetic energy balance.
◮ Copy the continuous kinetic energy balance:
- ∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u))
=0
2 ρ |u|2) + div
2 ρ |u|2) u
=0, with some formal algebra... using ∂tρ + div(ρ u) = 0.
◮ Do the same at the discrete problem ?
♭ Momentum on dual cells, mass on primal cells... ♯ Idea: reconstruct a mass balance on the the dual cells: ∀σ ∈ Eint, |Dσ| δt (ρn+1
σ
− ρn
σ) +
F n+1
σ,ǫ = 0
SLIDE 28 Discretization of the convection operator
◮ Choice of ρn+1
σ
and F n+1
σ,ǫ in
|Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
?
◮ Discretize
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nK ) so as to obtain a discrete kinetic energy balance.
◮ Copy the continuous kinetic energy balance:
- ∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u))
=0
2 ρ |u|2) + div
2 ρ |u|2) u
=0, with some formal algebra... using ∂tρ + div(ρ u) = 0.
◮ Do the same at the discrete problem ?
♭ Momentum on dual cells, mass on primal cells... ♯ Idea: reconstruct a mass balance on the the dual cells: ∀σ ∈ Eint, |Dσ| δt (ρn+1
σ
− ρn
σ) +
F n+1
σ,ǫ = 0
- btained from the primal mass balance for
◮ ρσ =
1 |Dσ|
|DK,σ| ρK + |DL,σ| ρL
- ◮ Fσ,ǫ : linear combination of the primal fluxes (FK,σ)σ∈E(K).
SLIDE 29 Discrete convection operator
◮ Continuous convection operator
C(ρ, u) =
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nσ).
SLIDE 30 Discrete convection operator
◮ Continuous convection operator
C(ρ, u) =
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nσ).
◮ Discrete convection operator
Cd(ρd, ud) = |Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
◮
ρσ Fσ,ǫ
- same as in discrete dual mass balance
◮ uǫ = 1
2 (un+1 σ
+ un+1
σ′ )
SLIDE 31 Discrete convection operator
◮ Continuous convection operator
C(ρ, u) =
ρn+1un+1 − ρnun+1 δt +
(ρn+1 un+1 ⊗ un+1 · nσ).
◮ Discrete convection operator
Cd(ρd, ud) = |Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σ un σ) +
F n+1
σ,ǫ un+1 ǫ
◮
ρσ Fσ,ǫ
- same as in discrete dual mass balance
◮ uǫ = 1
2 (un+1 σ
+ un+1
σ′ )
SLIDE 32 Continuous and discrete kinetic energy balance
◮ Continuous setting: Multiply continuous momentum by u:
- ∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u)) = 0
- · u
SLIDE 33 Continuous and discrete kinetic energy balance
◮ Continuous setting: Multiply continuous momentum by u:
- ∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u)) = 0
- · u
continuous kinetic energy balance: ∂t(1 2 ρ |u|2) + div(1 2 ρ |u|2) u + ∇p · u − divτ(u) · u = 0 (kin.en) ... with some formal algebra... using ∂tρ + div(ρ u) = 0.
SLIDE 34 Continuous and discrete kinetic energy balance
◮ Continuous setting: Multiply continuous momentum by u:
- ∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u)) = 0
- · u
continuous kinetic energy balance: ∂t(1 2 ρ |u|2) + div(1 2 ρ |u|2) u + ∇p · u − divτ(u) · u = 0 (kin.en) ... with some formal algebra... using ∂tρ + div(ρ u) = 0.
◮ Discrete setting: Similarly, multiply discrete momentum by un+1
σ
: |Dσ| δt (ρn+1
σ
un+1
σ − ρn σun σ) +
F n+1
σ,ǫ un+1 ǫ
+ |Dσ|(∇pn+1)σ,i − |Dσ|(divτ(un+1))σ = 0
SLIDE 35 Continuous and discrete kinetic energy balance
◮ Continuous setting: Multiply continuous momentum by u:
- ∂t(ρ u) + div(ρ u ⊗ u) + ∇p − div(τ(u)) = 0
- · u
continuous kinetic energy balance: ∂t(1 2 ρ |u|2) + div(1 2 ρ |u|2) u + ∇p · u − divτ(u) · u = 0 (kin.en) ... with some formal algebra... using ∂tρ + div(ρ u) = 0.
◮ Discrete setting: Similarly, multiply discrete momentum by un+1
σ
: |Dσ| δt (ρn+1
σ
un+1
σ − ρn σun σ) +
F n+1
σ,ǫ un+1 ǫ
+ |Dσ|(∇pn+1)σ,i − |Dσ|(divτ(un+1))σ = 0
discrete kinetic energy balance: 1 2 |Dσ| δt
σ
(un+1
σ
)2 − ρn
σ(un σ)2
+ 1 2
F n+1
σ,ǫ un+1 σ
un+1
σ′
+|Dσ| (∇pn+1)σ,i un+1
σ
−|Dσ|(divτ(un+1))σ un+1
σ
+Rn+1
σ
= 0 with Rn+1
σ
≥ 0) (kin.en)σ ... with some algebra... using |Dσ| δt (ρn+1
σ
− ρn
σ) +
F n+1
σ,ǫ = 0.
SLIDE 36 Other properties
◮ Discrete potential inequality ◮ Discrete entropy inequality ◮ Passage to the limit on the scheme (1D)
SLIDE 37 Full compressible Navier Stokes or Euler equations
Complete Navier-Stokes equations. Monophasic case Objective – derive a scheme for Euler (or Navier-Stokes) equations which is a natural extension of the scheme for low Mach number flows:
◮ staggered discretization, ◮ upwinding with respect to the material velocity, ◮ solution of the internal energy balance, ◮ fractional time step scheme.
By product : positivity of the internal energy.
SLIDE 38
Problems with shocks for non viscous flows
Internal energy equation yields wrong shock velocities Analogy with Burgers equation: for regular positive solutions (B) : ∂tu + ∂x(u2) = 0 ⇐ ⇒ (BS) : ∂tu2 + 4 3 ∂xu3 = 0. No longer true with irregular solutions: Rankine Hugoniot gives σ = uℓ + ur and σ = 4
3(uℓ + ur ).
Weak solutions of (B) = weak solutions of (BS).
SLIDE 39
Explicit upwind scheme for (B) and (BS)
u0 ≥ 0 Godunov ⇐ ⇒ upwind
Figure: Upwind Scheme for (B) (left) and (BS) (right) with different mesh sizes, CFL = 1.
SLIDE 40
How to get the correct shocks starting from the wrong equation
∂tu + ∂x(u2) = 0, (B) + initial condition Add a small diffusion term ∂tu + ∂x(u2) − ε∂xxu = 0. (B)ε Multiplying by 2u: ∂tu2 + 4 3 ∂xu3 − 2uε∂xxu = 0. (BS)ε For ε = 0, ∂tu2 + 4 3 ∂xu3 = 0. (BS)
SLIDE 41 Explicit upwind scheme
The explicit upwind scheme on (BS) is formally equivalent to: ∂tu2 + (4/3)∂x (u3) − ∂x((2hu2 − 4ku3)∂xu) = 0.
∂tu + ∂x(u2)− 1 u ∂x((hu2 − 2ku3)∂xu)
D D : non conservative numerical diffusion.
:-( Non conservative numerical diffusion on (B) yields wrong shock velocity...
but
:-) Could a non conservative numerical diffusion on (BS) yield a correct shock velocity on
(B) ?
SLIDE 42 Non conservative numerical diffusion on (BS)
Discretize (BS)ε instead of (BS): (BS)ε with ε = ε0hα ∂t(u2) + 4 3∂x(u3) − uε0hα∂x(2u∂xu) = 0, Centered finite volume
i
2 =
i
2 + 4k 3h u(n−1)
i−1
+ u(n−1)
i
2 3 − u(n−1)
i
+ u(n−1)
i+1
2 3 + k h2 ε0hα u(n−1)
i
i−1
− 2u(n−1)
i
+ u(n−1)
i+1
SLIDE 43
Numerical results with centered scheme
Figure: Centered Scheme for (BS)ε0hα, α = 1.
SLIDE 44
Figure: Centered Scheme for (BS)ε0hα, α = 1.5.
SLIDE 45
Figure: Centered Scheme for (BS)ε0hα, α = 2.
SLIDE 46 From Burgers to Euler
For regular solutions, ∂tu + ∂x(u2) = 0 ⇐ ⇒ ∂t(u2) + 4 3 ∂xu3 = 0. ∂t̺ + div(̺u) = 0, ∂t(̺u) + div(̺u ⊗ u) + ∇p = divτ, ∂t(̺E) + div
⇐ ⇒ ∂t̺ + div(̺u) = 0, ∂t(̺u) + div(̺u ⊗ u) + ∇p = divτ, ∂t(̺e) + div(̺eu) + p divu = τ : ∇u, Additional difficulty
◮ we had an equation, we now have a system...
Inspiration comes from the formal derivation of the internal energy equation...
SLIDE 47 Euler equations. . . and ”derived” forms
◮ Euler (Navier-Stokes) equations:
∂t̺ + div(̺u) = 0, (mass) ∂t(̺u) + div(̺u ⊗ u)−divτ + ∇p = 0, (mom) ∂t(̺E) + div
(tot.en) p = (γ − 1) ̺e, E = 1 2 |u|2 + e.
◮ For regular functions, (mom) ·u & (mass) (kin.en):
1 2 ∂t(̺|u|2) + 1 2div(̺|u|2u) + ∇p · u = div(τ) · u. (kin.en) Subtracting from (tot.en) yields the internal energy balance: ∂t(̺e) + div(̺eu) + p divu = τ : ∇u, (int.en) which implies e ≥ 0.
SLIDE 48 Strategy
How to obtain the correct weak solutions of Euler equations while solving the internal energy balance ? General idea: (kin.en)d + (int.en)d + correction term (tot.en)d More precisely: 1- Kinetic energy balance (with residual term) obtained from (mom)d and (mass)d. 2- Suppose bounds and convergence for a sequence of discrete solutions, compatible with the regularity of the sought continuous solutions:
◮ control in BV and L∞, ◮ convergence in Lp, for p ≥ 1.
3- Let ϕ a regular function,
◮ interpolate, ◮ test the kinetic energy balance, ◮ test the internal energy balance, ◮ and pass to the limit in the scheme.
. . . and, on the basis of this computation, build corrective terms in the internal energy balance in such a way to recover, at the limit, the weak form of the total energy equation.
SLIDE 49 An implicit scheme
|K| δt (ρn+1
K
− ρn
K ) +
F n+1
K,σ = 0,
(mass)d |Dσ| δt (ρn+1
σ
un+1
σ
− ρn
σun σ) +
E(Dσ)
F n+1
σ,ǫ un+1 ǫ
+ |Dσ| (∇pn+1)σ = 0, (mom)d |K| δt (ρn+1
K
en+1
K
− ρn
K en K ) +
F n+1
K,σ en+1 σ
+ |K| pn+1
K
(divun+1)
K = Sn+1 K
, (int.en)d pn+1
K
= (γ − 1)ρn+1
K
en+1
K
. Existence, stability
◮ If the solutions exists, and if SK ≥ 0, e ≥ 0. ◮ Kinetic energy conservation + internal energy conservation stability estimate
(provided that SK may be absorbed).
◮ ∃ at least a solution (topological degree argument), and e ≥ 0.
SLIDE 50 Choice of SK
◮ (kin.en)σ:
|Dσ| 2δt (̺n+1
σ
|un+1
σ
|2 − ̺n
σ|un σ|2) + 1
2
F n+1
σ′ un+1 σ
un+1
σ′
+ (∇p)n+1
σ
· un+1
σ
= −Rn+1
σ
, with: Rn+1
σ
= |Dσ| 2δt ̺n
σ |un+1 σ
− un
σ′|2
◮ Multiply kinetic egergy balance and internal energy balance by adequate interpolates of
test function ϕ
(kinetic terms)σϕσ +
(energy terms)K ϕK =
δt|SK |ϕK +
δt|Rσ|ϕσ. The terms at the left hand side tend to −
- Ω×(0,T)
- ρ E ∂tϕ + (ρ E + p) u ∇ϕ
- dx −
- Ω
ρ(x, 0) E(x, 0) ϕ(x, 0) dx, ...
δt|Rσ|ϕσ +
δt|SK |ϕK → 0 as mesh size → 0.
SLIDE 51 Choice of the source term
Example for the implicit scheme, with centered choice for uσ: Rn+1
σ
= − 1 2 |Dσ| δt ρn
σ
σ
− un
σ,i
2 ∀K ∈ M, Sn+1
K
= 1 2
|DK,σ| δt ρn
K
un+1
σ
− un
σ
2 With this choice,
δt|Rσ|ϕσ +
δt|SK |ϕK → 0 as mesh size → 0. thanks to the fact that |Dσ|ρσ = |DK,σ| ρK + |DL,σ| ρL
SLIDE 52 A fractional time step scheme
Prediction step: 1 δt (̺n˜ un+1 − ̺n−1un) + div(̺n˜ un+1 ⊗ un)−divτ(˜ un+1) + ̺nt ̺n−1 1/2 ∇pn = 0, Correction step:
δt (un+1 − ˜ un+1) + ∇p − ̺n ̺n−1 1/2 ∇pn = 0, 1 δt (̺n+1 − ̺n) + div(̺n+1un+1) = 0, 1 δt (̺n+1en+1 − ̺nen) + div(̺en+1un+1) + pn+1 divun+1 = τ(˜ un+1) : ∇˜ un+1 + S, pn+1 = ℘(̺n+1, en+1) = (γ − 1) ̺n+1en+1, where S is again computed so that
δt
Rσϕσ −
SK ϕK
→ 0.
SLIDE 53 1D Riemann problem
◮ Scheme 1, scheme 3, and what is obtained without corrective term:
5 10 15 20 25 30 35 40
0.2 0.4 density x scheme 1 scheme 3 S=0 exact
[E. Toro, Riemann solvers and numerical methods for fluid dynamics, third edition, test 5 of chapter 4].
SLIDE 54 2D Riemann problem
SLIDE 55 Contact discontinuities
1D: does the scheme keep the pressure and the velocity constant across contact discontinuities ? . . . i.e. check whether a constant velocity and pressure may be solution to the scheme: ∀σ ∈ E |Dσ| δt (̺n
σ˜
un+1
σ
− ̺n−1
σ
un
σ) +
F n
σ,ǫ ˜
un+1
ǫ
+ ̺n
σ
̺n−1
σ
1/2 (∇pn)
Dσ = 0.
mass balance over Dσ, (∇pn)
Dσ = 0 ˜
u = cste
|Dσ| δt (̺n
σuσ − ̺n−1 σ
˜ un+1
σ
) + (∇pn+1)
Dσ −
̺n
Dσ
̺n−1
Dσ
1/2 (∇pn)
Dσ = 0,
un+1 = ˜ un+1 ∀K ∈ M |K| δt (̺n+1
K
− ̺n
K ) +
F n+1
K,σ = 0
∀K ∈ M |K| δt (̺K en+1
K
− ̺n
Ken K ) +
F n+1
K,σ en+1 σ
+ pn+1
K
|σ|un+1
σ
· nK,σ = SK divu, SK = 0 ̺e = p = cste ∀K ∈ M pn+1
K
= (γ − 1) ̺n+1
K
en+1
K
,
SLIDE 56 Perspectives
◮ Integration in the home made open source software, ISIS.
ISIS: https://gforge.irsn.fr/gf/project/isis
◮ Extensive tests beginning, for problems ranging from hyperbolic to low Mach number
flows.
◮ Explicit variant of the scheme (Ph.D. Trung Tan Nguyen). ◮ Do the schemes converge to an entropy weak solution ? ◮ To develop: less diffusive schemes.
Similar scheme for barotropic flows application to (low Mach number) shallow water problems ?
◮ Extension to reactive flows (deflagration, detonation) starting now .
