Transcript Nessun titolo diapositiva - DISMA Dipartimento di Scienze

Artificial dissipation methods
for Euler equations
V. Selmin
Multidisciplinary Computation and Numerical Simulation
Euler equations properties
Conservation law form:
W
   F (W )  0
t
Unknowns & Flux vector:
w
  






T
W    w  , F    ww  pI 
e 
  wh

 t

t

State equation:
p  f (  , e)
p  (  1)  e ( perfect gas)
Non conservative form:
W
 A  W  0
t
F
A
W
Euler equations
mass density
:
velocity
: w  (u , v, w)
pressure
:p
total energy
: et
total enthalpy : ht  et 
internal energy : e  et 
spec. heats ratio : 
p

ww
,
2
Euler equations properties
Hyperbolicity:
Let F  F  ν be a linear combination of the fluxes where ν is a nonzero vector.
The system is hyperbolic if the the Jacobian matrix A  A  ν has only real eigenvalues.
Eigenvector and eigenvalues:
The eigenvectors form the matrix T that diagonalises A :
T
1
A T   ,  lm  (l ) lm
  
(1)
( nd 1)

( nd )
sound speed : c
 wν ,
( nd  2)
 wνcν , 
spec. entropy : s
 wνcν
Jacobian matrix decomposition:
The following decompositions is also introduced:
A  A  A
A  A  A


The matrices A and A are diagonalised by the transformation T :
A   T  T
1
,

,   diag. max( (l ) ,0)

Flux homogeneity
Euler fluxes are homogeneous of degree one with respect to W
F  AW
Euler equations

,   diag. min( (l ) ,0)

Finite Volume/Finite Element
Node-Centred finite volume
Discretisation of the integral equation:
vol (Ci )
dWi
  F  n d  0
Ci
dt
Discrete equation:
vol (Ci )
dWi

dt
 (Wi ,W j , ij )  F(W i )  ν i  0
jK (i )
ν ij  
Ci C j
νi  
Ci 
n d
Consistency:
(W ,W , ν)  F(W )  ν
Conservative systems:
ν ij   ν ji
 ν ij  0
 i
 ν ij   ν i
 i
jK ( i )
jK ( i )
n
Second-order spatial approximation:
(Wi ,W j , ν ij ) 
Finite Volume/Finite Element Formulation
F(Wi )  F(W j )
2
n
 ν ij
n d
Finite Volume/Finite Element
Positivity
Scalar equation:
w
  f  0
t
Discrete equation:
vol (Ci )
dwi

dt
 (wi , w j , ν ij )  0
jK (i )
Condition on scheme positivity:
 ij  0 , vol (Ci )
dwi

dt
Monotone numerical scheme:
 ( wi , w j , ν ij ) 
f ( wi )  f ( w j )
2
 ν ij 
a( wi , w j , ν ij )
2
Extension to systems of equations:
(Wi ,W j , ν ij ) 
F(Wi )  F(W j )
2
( w j  wi )
 ij (w j  wi )
jK ( i )
 f ( w j )  f ( wi )
 ν ij
 w w
j
i

a( wi , w j , ν ij  
 f  ν 
 w ij  w

i
1
 ν ij  T ij  ij T ij1 (W j  Wi ) ; T  T
2
A (Wi ,W j , ν ij ) (W j  Wi )  F(W j )  F(Wi )  ν ij
Finite Volume/Finite Element Formulation
1
A
; w j  wi
, w j  wi
Entropy Condition
The weak solutions are not univocally defined for the partial differential problem: different solutions may satisfy the
initial condition.
An additional principle is then needed for the settlement of a solution physically acceptable.
Let consider the hyperbolic equation
w f ( w)

0
t
x
where f(u) is a nonlinear function of the uknown.
The eligible solution is the limit when   0 of the solution of the viscous problem:
w f ( w)

w 

    ( w)  ,  ( w)  0,   0
t
x
x 
x 
Oleinik has shown that the eligible solution must satisfy the following inequalities
f ( w2 )  f ( w) f ( w2 )  f ( w1 ) f ( w)  f ( w1 )


w2  w
w2  w1
w  w1
for w1  w  w2. These relations are called entropy condition.
Entropy condition
Flux Difference Splitting
A central role is played in the algorithm by the modelling of the interactions of fluxes at the interfaces of adjacent
grid cells.
Two distinct sets of state quantities must be combined into one single set of numerical fluxes at the interface.
In the following, some upwind related tecnologies which may be used to build the numerical flux will be discussed.
Flux Difference Splitting:
Let consider the Riemann problem governed by the one-dimensional inviscid flow equation:
W F (W )

0
t
x
Wl x  0
W (t  0)  
Wr x  0
Roe’s scheme may be interpreted as an approximate Riemann solver based on the linearised equation:
W
W
 A(Wl ,Wr )
0
t
x
The matrix A(Wl ,Wr ) has to satisfies the following properties:
1- A(W ,W )  A(W )
(Consistency)
2- A(Wl ,Wr ) has a complete set of real eigenvalues and eigenvectors (Hyperbolicity)
3- A(Wl ,Wr ) (Wr  Wl )  F (Wr )  F (Wl )
(Conservation propertty)
4- When the two states are connected through one discontinuity, the discontinuity can be reproduced exactly.
Conditions 1 and 3 have to be satisfied when the Jacobian is computed at an average state:
A(Wl ,Wr )  A(W )
Flux Difference & Vector splittings
Flux Difference Splitting
Roe introduced the following average operator for a polytropic ideal gas:
a () 
 r () r  l ()l
 r  l
 u
 r ht ,r  l ht ,l
 r u r   l ul
, ht 
 r  l
 r  l
W
Wl  Wr
2
The average operator is in general valid for a gas whose pressure law allows the flux F to be an homogeneous
function of degree one in the conservative variables. As the pressure is expressed in the general form
p   (  ,  u,  et )
that means the following condition has to be satisfied:




 u
  et
 

 u
 et
u2
p   g (et  )  cst
2
g : arbitrary function
The numerical flux on this approximate Riemann solver may be expressed by
1
 R (Wl ,Wr )  F (Wl )  F (Wr )  A(Wl ,Wr ) (Wr  Wl ) 
2
Numerical flux not differentiable
Crisp representation of discontinuities
Does not satisfy the entropy condition
Allows expansion shocks
Another approximate Riemann solver , with similar properties, is due to Osher:
W
1
 O (Wl ,Wr )   F (Wl )  F (Wr )   r A(W ) dW 
Wl

2
Flux Difference & Vector splittings
Numerical flux differentiable
Integration path // to the eigenvectors of A
Satisfy the entropy condition
Flux Vector Splitting
In the case of Flux Vector Splitting (FVS), the numerical flux across an interface separating two states Wl and Wr
is considered as the sum of two parts F


1- F (W )  F (W )  F (W )


and F that should possess the following desirable properties:
(consistency to the governing equations),


2- dF / dW must have all the eigenvalues  0 , dF / dW must have all the eigenvalues  0



3- F must be continuous with F (W )  F (W ) for M  1, F (W )  F (W ) for M  1 where M is the
Mach number normal to the interface


4- The components of F and F together must reproduce the symmetry of F with respect to M (all other
quantities held constant), that is
Fk ( M )   Fk ( M ) if Fk ( M )   Fk ( M )

5- dF / dW must be continuous
Restriction (3) ensures that in supersonic regions leads to standard upwind differencing.
F  are homogeneous of degree one.
dF  (W )

, D (W ) 
dW
In general, the flux decomposition is chosen in such a way that
F  (W )  D  (W ) W
The numerical flux may be expressed as
 (Wl ,Wr )  F  (Wl )  F  (Wr )
which can be rewritten to make it to look like the Osher’s formula
W
1
(Wl ,Wr )   F (Wl )  F (Wr )   r D(W ) dW  ,
Wl

2
D(W )  D  (W )  D  (W )
The integration is independent of the path since D(W ) is a perfect gradient
Flux Difference & Vector splittings
Test case: hypersonic flow around a cylinder
The test case corresponds to compute a M   8 flow around a cylinder section.
In this case, a detached choc of high intensity (from M=8 to M=0.3929) is present in front of the body.
The flow, being subsonic behind the shock, becomes supersonic at approximatively 29% of the cylinder radius.
Enlargement of the grid
Mach number contours
Computational domain (65x41) nodes
Flux Difference & Vector splittings
Symmetry axis+body nodes
Roe’s scheme & the entropy condition
As already anticipated previously, the numerical flux of the Roe scheme may be written for multidimensional
problems according to
(Wi ,W j , ν ij ) 
F(Wi )  F(W j )
2
1
 ν ij  T ij ij T ij1 (W j  Wi ) ; T  T
2
1
A
where the quantities in the dissipation term are computed by using the average proposed by Roe and
T
1
A T   ,  lm  (l ) lm
(1)    ( nd )  w  ν ,
( nd 1)  w  ν  c ν , ( nd  2)  w  ν  c ν
To have a zero value of the eigenvalues independently of the difference (W j  Wi ) may lead to a violation of
entropy condition, whereas the dissipation matrix may become singular when
“Carbuncle” effect in Roe scheme
Flux Difference & Vector splittings
Expansion discontinuity
( k )  0 .
Roe’s scheme & the entropy condition
Harten proposed to regularize the dissipation matrix by using the following correction of the eigenvalues:
( k )
H
 ( k )

  (k ) 2

 2


2
( k )  
( k )  
Entropy condition
not satisfied
Entropy condition
satisfied
Flux Difference & Vector splittings
Multidimensional extension
A method to extend the upwind schemes to multidimensional problems is based on the property of the
Euler’s equations to be invariant with respect to a rotation of the coordinate axis.
Let F  F  ν be a linear combination of the fluxes where ν is a nonzero vector of real quantities.
A local rotation of the velocity vector is defined in the following way:
u
 u 
u
 
 
 
ν
w   v   w   v   R  v  , u  w 
ν
 w
w 
 w
 
 
 
The state W  (  ,  u ,  v,  w,  et ) is transformed into the state Ŵ
W  R W  (  ,  u ,  v ,  w ,  et )  Wˆ
Moreover, applying the the rotation R to the flux allows to describe the flux vector in a one dimensional fashion
  u



2
  u  p 


F (W , ν )  RF    u v   F (Wˆ )
  u w 


  u ht 


Finally, the inverse transformation express the different flux in the original frame
F (W , ν)  ν R 1F (Wˆ )
The eigenvalues of the Jacobian matrix are not modified
Flux Difference & Vector splittings
Matrix Dissipation Operator
The matrix dissipation flux
Fmd (Wi ,W j , ν ij )  T ij  T ij1 (W j  Wi )
may be described in a straightforward way for the Euler equations.
In the rotated frame, the matrix dissipation operator takes the form
dp


 du








2


c

du

u
dp
u

du

dp










  a2 
v
dp

a
v

du

3





 c 



w dp
w  du 






 c 2u  du  h dp 
h

du

u
dp

 t



 


t
w w

, dp  (  1)     w  w v  et  , du  u  u 
 2

 

  u
Fmd (Uˆ )  c a1   v

  w

  et
where
q  q j  qi
It may be demonstrated that by using Roe’s averages
w
and
 j w j  i w i
 j  i
, ht 
 j ht , j  i ht ,i
  i  j , we obtain
dp  p ,  du   u
Matrix Dissipation Operator
 j  i
Matrix Dissipation Operator
The matrix dissipation flux becomes
 

  u
Fmd (Wˆ )  c a1   v

  w

  et
p


u








2


c u  u p
 u u  p 

a


 2
v p
 a3 
v u 




 c



w p
w u 






2


h


u

u

p

 t



 c u u  ht p 
The coefficients ak are functions of the eigenvalues of the system and take the form:
a1 
(1)
c
, a2 
( 4)  (5)  2 (1)
2c
, a3 
( 4)  (5)
2c
which may be expressed with respect to the Mach number normal to the interface according to
a1  M

1  M ; M  1
, a2  
; M  1

0
Matrix Dissipation Operator

 M
, a3  

sign( M )
; M  1
; M  1
a3  (1  a2 ) sign( M )
Matrix Dissipation Operator
a1  M

1  M ; M  1
, a2  
; M  1

0

 M
, a3  

sign( M )
; M  1
; M  1
• The first term of the matrix dissipation operator is able alone to maintain the monotonicity of Roe scheme
• The second term is a diffusive operator acting only in subsonic regime.
• The third term is an antidiffusive operator that acts both in subsonic and supersonic regime.
Matrix Dissipation Operator
Matrix Dissipation Operator
Matrix Dissipation Operator
a1 , a 2 , a 3
a1 , a 2  0, a 3  0
a1 , a 2 , a 3  0
a 1 , a 2  0, a 3
Matrix Dissipation Operator
In order to enforce the entropy condition, the coefficients may be regularised according to
a1reg
 M

  M 2  2

1

 21
Matrix Dissipation Operator
1  M
, M  1

2
2
  M  1   M  1 2   2
reg
; a2  
4 2
, M  1

0

a3  (1  a2 ) sign( M )
, M  1   2
, 1   2  M  1   2
, M  1   2
Isoenergenetic flows
A desirable property of a numerical scheme for the solution of the Euler equations is its capability to reproduce an
isoenergenetic steady state flow.
Total enthalpy constant over the entire field
Whether the discretisation of the convective flux might provide a solution which satisfies the solution, this no more
true when the dissipation operator Fmd is added to the numerical scheme.
To obtain a scheme which satisfies this property, the following dissipative flux is introduced

  



 c2
 u 
a 
h
Fmd
(Wˆ )  c a1   v   2 

 c


w

 





h

t 

c a1  et  

p

u





u  u p 
u


u


p
 




v p
v u 
  a3 



w p
w u 




ht p
h


u


t



a2 2
a
c u u  ht p  a3 ht u  u p   c a1  ht   2 ht p   a3 ht u 
c
c
Properties:
h
1- The flux Fmd preserves the monoticity of the numerical flux.
2- Essential property in order to obtain accurate stagnation temperature
3- Enthalpy damping techniques may be used to speed up convergence of transonic flows solution.
Matrix Dissipation Operator
Isoenergenetic flows
h
A more compact expression of Fmd may be obtained by rewritting the scheme with respect to the coefficients
a1 and a 2 only.
a3  (1  a2 ) sign( M )
Characteristic variables:
 I1  p  c 2
 I 2  v
 I 3  w
 I 4  p   c  u
 I 5  p   c  u
  I 4,5  p  sign(M )  c  u
h
Fmd
may be expressed according to
 

 u
h
Fmd
(Wˆ )  c a1   v

  w

  ht
Matrix Dissipation Operator
1



 u


c

u

sign
(
M
)
a

I


2
4,5 
  sign( M )
v

 

c
 w




 ht

0

 

1 
  sign( M )(p  a  I ) 0 

2
4,5

 

0

0
 

Scalar Dissipation
In the context of centered approximations, a simpler approximation for the dissipation operator may be
obtained by replacing in the matrix dissipation operator all the eigenvalues of the Jacobian by its spectral radius:
 (A )  u  c  c  M  1
The modification leads to the following scalar dissipation operator
  



u

 
Fsdh (Wˆ )  c  M  1   v 



w

 


  ht 
which preserves the monoticity of the scheme,but it is more dissipative. In the original frame, the scalar dissipation
operator may be transformed according to:
 
u 


h
1 h ˆ
Fsd (W , ν )  ν R Fsd (W )   c ( M  1) v 


w 
ht 
Scalar Dissipation Operator
Van Leer Flux Vector Splitting
A very promising study for splitting the flux has been introduced by Van Leer. He deduced formulas that satisfy the
basic requirements , but with the following additional properties:

1- dF / dW must have one eigenvalue vanishing for M  1 .

2- F (M ), like F (M ) , must be polynomials in M and of the lowest degree.
Constraint (1), that van Leer considered crucial, makes it possible to build numerical shock structures using two
interior zones. The requirement (2) makes the splitting unique.
The formulas for F  are given in terms of the Mach number.
For supersonic and sonic flows, M  1, we have
M  1 : F   F , F   0 , M  1 : F   0, F   F
and for subsonic flow, M  1 ,


 c

2

(
M

1
)

 F1   4

  
(


1
)
u

2
c




F   F2    F1

  

 F3  

    (  1) u  2c 2 
 F1
2( 2  1) 

Split fluxes continuously differentiable
Valid only for a perfect gas
Flux Vector Splitting
Montagné Flux Vector Splitting
A more physical way to build a flux splitting consists in approaching the conservative laws by a collisionless
Boltzmann equation. Each state is split into group of particles which are distributed according to their velocities.
The density law of particles, which obeys a Maxwell distribution in the kinetic theory of gases, must be approximated
by a simpler distribution law. The fluxes F  and F  are defined by considering the particles with positive and
negative velocities, respectively. In addition, it is requested that F  ( F  ) reduces to F for supersonic flows with a
positive (negative) velocity.
The formulas for F  are given in terms of the Mach number.
For supersonic and sonic flows, M  1, we have
M  1 : F   F , F   0 , M  1 : F   0, F   F
and for subsonic flow, M  1,




p
p




( M  1)
u  ( M  1)


2c
2
c






p
p
 , F  
F  
( M  1) 2
u 2  (1  2 M  M 2 )

2
2




 

2
p
c
p
 


 F1 e  p  c ( M  1) 2  
F1 (e 
)  ( c 2  p ) M 3  ( M  1)3  



2 2 
2



  2 2
1  M  0
0  M  1
F  F  F
Split fluxes not always continuously differentiable
Valid not only for a perfect gas
Property (1) not satisfied
Flux Vector Splitting
Enthalpy preserving Flux Vector Splitting
The new splitting is based on the following heuristic consideration:
1- The fluxes F1 of van Leer scheme are continuously differentiable and are expressed in a general form (OK).
2- The fluxes F 2 of the Boltzmann type scheme are constinuously differentiable and are very close to the
corresponding fluxes of the van Leer type (OK).
3- For the components of the third fluxes, any of the the previous schemes can be retained. In addition, they do
not satisfy the condition of constant total enthalpy for isoenergenetic steady flows. They can be repaired by




taking F3  F1 ht . Moreover, if fluxes F1 are continuously differentiable, it will also the case for F3 .
For supersonic and sonic flows, M  1 , we have
M  1 : F   F , F   0 , M  1 : F   0, F   F
and for subsonic flow, M  1,
F1  
c
4
( M  1) 2
p
( M  1) 2 ; F2  u  F2 ;  1  M  0
2
p
F2  u  F2 ; F2  ( M  1) 2 ; 0  M  1
2
F3  F1 ht
F2 
Split fluxes continuously differentiable
Valid for any gas
Very simple expressions
Flux Vector Splitting
Multidimensional extension
Applying the rotation R to the flux F allows to split the flux in a one-dimensional fashion:
F (W , ν)  RF  F (Wˆ )
After splitting F (Uˆ ), the inverse transformation
F (W , ν)  ν R 1F (Wˆ )
gives the following relations for the flux splitting:
F  (W , ν)  ν R 1F  (Wˆ ) , F  (W , ν)  ν R 1F  (Wˆ )
Then the numerical flux reads
(Wl ,Wr , ν)  F  (Wl , ν )  F  (Wr , ν ) 

 
F(Wl )  F(Wr )
1
 ν  F  (Wr )  F  (Wr )  F  (Wl )  F  (Wl )
2
2
For multidimensional applications, the proposed flux splitting is described by the following formulas.
Indicating by M the quantity u / c , for M  1, we have
M  1 : F   F , F   0 , M  1 : F   0, F   F
and for subsonic flow,
F1    c / 4 ( M  1) 2
F2  p / 2 ( M  1) 2 ;  1  M  0
F2   u  p / 2 ( M  1) 2 ; 0  M  1
F3  F1 v
F4  F1 w
F5  F1 ht
Flux Vector Splitting
F  F  F

Extension to second order spatial accuracy
All upwind schemes (FDS, FVS) may be interpreted as built from a centred discretisation (second-order accurate)
spatial approximation + an artificial dissipation term (first-order accuracy)
(Wi ,W j , ν ij ) 
(Wi ,W j , ν ij ) 
F(Wi )  F(W j )
2
F(Wi )  F(W j )
2
1
 ν ij  Qij (Wi ,W j , ν ij ) (W j  Wi )
2
 ν ij 

FDS
 
1
F  (W j )  F  (W j )  F  (Wi )  F  (Wi )
2

FVS
The schemes are first-order accurate in space.
In general, artificial dissipation is needed only at discontinuities, thus development of sensors that identify
discontinuities (or rapid change of the solution) and monitor the level of dissipation by using the sensor.
The quantities Wi and W j are replaced by the following expressions
Wj
Wj
where s and


si
(W j  Wi )    2Wi
2
sj
 Wij  W j 
(W j  Wi )    2W j
2
 Wij  Wi 

 are parameters,
 2Wi  W j  2Wi  Wi  
 2W
h
2
x 2
 W j  W j   2W j  Wi 
j+

j
i
i-
2
i+2
i+1
i
i-1
and Wi  ,W j  are “extrapolated” values of the unknowns.


N.B.: The values of Wij and Wij may be used in the evaluation of all the terms of the FDS or FVS schemes
or only for the evaluation of the dissipative term.
Second order accuracy
Extension to second order spatial accuracy
Wij  Wij may be expressed according to
s j  si

W j  Wi  Wij  Wij  (1 
) (W j  Wi )  s j 2W j  si 2Wi
2
2
The variation

Thus,
s  0  Wij  Wij  W j  Wi
s 1 
Wij
 Wij

 W
2

2
First-order scheme
  Wi
2
j



3
 h
2 x 3
3
Higher-order scheme
Then, in order to avoid oscillations in the numerical solution, in particular near shocks, the extrapolation has to be
limited by varying the switching factor s , which is controlled by the solution itself.
Typical formulation of such limiters are given by
van Albada : s 
van Leer :
s
WW  WW
(W ) 2  (W ) 2  

2 WW  WW

Wi  W j  Wi , W j  W j   W j
Wi  Wi  Wi  , W j  W j  Wi
(W  W ) 2  
where  is a small number (   106 ) preventing division by zero in regions of null gradients.
These limiter formulations react to the local gradient and to the curvature of the solution, i.e.
2
2
2
2   W / x 

1  s  h 


W
/

x


Thus, in regions of weak changes of W , the limiter remains nearly to one, but for strong changes the value
decreases and the scheme reduces the accuracy generating a stronger numerical simulation.
Second order accuracy
Extension to second order spatial accuracy
Methods to computed extrapolated variables
1- Centred gradient:
Wi   W j  2(x j  xi )  Wi
W j   Wi  2(x j  xi )  W j
vol (Ci ) Wi 

Wi  W j
jK (i )
2
j*
ν ij  B.T .
j+
j
This formulation works well for subsonic and transonic flows, but not for high
supersonic flows, for which it is not able to forecast the correct slope of the
solution upstream or downstream of the nodes i and j, respectively, since the
gradient approximation is centered.
W j   W j  (x j  xi )  We j
Wei 
i
ei
i-
2- Upstream/downstream element gradient:
Wi   Wi  (x j  xi )  Wei
ej
 N k Wk
i*
kei
This formulation works well for subsonic, transonic flows and supersonic flows.
It forecasts the correct slope of the solution upstream or downstream of the node i and j, respectively.
Nevertheless, it needs additional data structures for the computation of the gradient on the element.
3- Upstream/downstream segment variation:
Wi   Wi  (Wi  Wi* ) x i*  x i / x j  x i
W j   W j  (W j*  W j ) x j*  x j / x j  x i
The nodes i* and j* are selected in such a way that
the segments i-i* and j-j* minimise the scalar product
with the segment i-j
This formulation works well for subsonic, transonic flows and supersonic flows.
It forecasts an approximation of slope of the solution upstream or downstream of the node i and j, respectively.
It gives a good information for the computation of the limiter s, but the approximation for the third-order operator is
less accurate. It needs very few additional data structure information to work.
Second order accuracy
Extension to second order spatial accuracy
Alternative methods to compute the third-order term
s
s

Wij  Wi  i (W j  Wi )    2Wi 
Wij  Wi  i (W j  Wi )    2Wi

2
2
 
sj
sj
Wij  W j 
(W j  Wi )    2W j 
Wij  W j 
(W j  Wi )    2W j

2
2





where
 ij (W j  Wi )
 Wi 
2
jK ( i )
 ij
jK (i )
Remarks:
1- The use of the operator

 




 ij  1
, 
 ij  (x j  xi )  ij
 2 leads to a final discretisation similar to those of a directional fourth-order derivate
4
2
4 W
  h
l 4
2- The use of the operator  leads to a final discretisation similar to those of a biharmonic operator
2
2
 h 4 2 ( 2W )
3- The use of the harmonic operator leads to a more robust discretisation for multidimensional grids with deformed
elements.
Second order accuracy
Extension to second order spatial accuracy
Special treatments for the scalar dissipation
 
u 


h
Fsd (W , ν )   c ( M  1) v 



w


ht 
In the classical scalar artificial method, it is preferred to use the following expression for the variation:
W   ij2 (W j  Wi )  2 ij4 ( 2W j   2Wi )
shock capturing
with
 ij( 2)
 min(1, k
( 2)
background dissipation
max( si , s j )) ,
 ij( 4)
 max( 0, k
( 4)

 ij( 2)
2
)
where the discontinuity sensors are expressed in terms of pressure variations
si 
p j  2 pi  pi 
p j  2 pi  pi 
, sj 
p j   2 p j  pi
p j   2 p j  pi
Remarks:
1- The aim of the background dissipation is to improve the ability of the scheme to damp high frequency modes in
smooth part of the flow.
( 2)
( 4)
2- k , k
are adjustable parameters.
( 2)
3- The value of  ij is limited to 1 because we know that with this value the scheme maintains monoticity.
( 2)
( 4)
4- The term  ij / 2 is substract from k
to cut off the fourth order difference operator, which could otherwise
cause an oscillation in the neighbourhood of a shock wave.
Second order accuracy
Modelling of High Speed flows
At hypersonic speed, a strong shock appears in front of an aerospace vehicle.
The compression of the gas through the shock is associated with an increase in temperature.
Through the shock wave , the total enthalpy is constant and along a streamline, we can write
w
h  
2
2
w
h
2
2
h: specific enthalpy
2
If the velocity w  is very large, h  w  .
A strong shock causes a large variation of the fluid velocity. Thus, in front of the vehicle, w  w  .
Consequently, we have
2
h
w
2
This means that almost all the kinetic energy is transformed into heat, resulting in a very high temperature in the
region between the body and the shock. An important characteristic of hypersonic phenomena is then the
transformation of an hypervelocity flow into an hyperenthalpy one.
Another characteristic is that a part of that energy is transformed into chemical energy, as a consequence of
(i) excitation of the internal degree of freedom of molecules,
(ii) dissociation of molecules,
(iii) ionisation.
These phenomena are endothermic and the final temperatures reached are much lower than those obtained for
a non-reactive gas. Moreover, in these conditions, the properties of air are modified.
The chemical composition of air can be approximately subdivided into the following regimes:
1- T<2500K. The chemical composition is substantially that at room temperature.
2- 2500<T<4000K. The oxygen dissociation regime, no significant nitrogen dissociation,slight nitric oxyde formation
3- 4000<T<8000K. The nitrogen dissociation regime, oxygen fully dissociated.
4- t>8000K. Inonization of the atomic constituents.
Modelling of high speed flows
Modelling of High Speed flows
The model of a real gas can be described as a mixture of perfect gases which are subject to chemical reactions.
The state law must take into account , on one hand, the species decomposition and, in another hand, the
exitation of internal degrees of freedom of atoms and molecules (translation, rotation, vibration).
The internal energy takes the form
T = temperature
hi0 = energy of formation
ns
ns
ns
ns
Yi = mass fraction
e   Yi ei   Yi Cv,iT   Yi hi0   Yi eivib
i 1
i 1
i 1
i 1
Cv,i = specific heat at constant volume
eivib= energy of vibration
Note that
3 / 2 Ri , atoms
Cv,i  
, Ri  R / mi
3
/
2
R
,
diatomic
molecules

i
R = universal gas constant
mi = molar mass
Using Dalton’s law, the pressure is given by
p  T
ns
 Yi Ri
i 1
The specific enthalpy is written according to
ns
ns
ns
i 1
i 1
i 1
h   Yi hi   Yi C p ,iT  
Modelling of high speed flows
Yi hi0
ns
  Yi eivib
i 1
C p,i  Cv,i  Ri
Modelling of High Speed flows
The following alternatives for the definition of the ratio of specific heats may be introduced:
1- True ratio of specific heats

Cp
Cv

h / T  p
e / T v
2- Equivalent  relating pressure to internal energy and density
  1
p
e
3- Equivalent  relating sound speed to pressure and density
 c2
~
 
p
4- Ratio of specific heats of the frozen flow
ns
ˆ 
 Yi (dhi / dT )
i 1
ns
 Yi (dei / dT )
i 1
Modelling of high speed flows
Non-reactive flows
A gas formed by standard air (oxygen+nitrogen) is assumed and the vibrational degrees of freedom of the diatomic
molecules are neglected. We obtain a perfect gas with constant ratio of specific heats, i.e.
    ~  ˆ  1.4
The pressure and the sound speed can be expressed in a simple way as a function of the conservative variables:
1
p
2
p  (  1)(  et   w ) , c 2  
2

Remarks:
1- The upwind schemes described previously have been build in such a way they satify the condition of total
enthalpy at steady state. That is true for first-order scheme. But, special care must be taken when employing
the second order one. The standard procedure is to adopt the primitive variables (  , w, p) as sensor for
the limiters, which does not satisfy the property. The total enthapy per unit volume ht has to be used
instead of pressure.
2- The upwind formulas depends of five quantities (  , w, p, c, ht ). The limiting procedure is performed on
only three of them. This is possible because a constant specific heats ratio is assumed.
3- For hypersonic flows, a drawback is the use of the static speed of sound that control the eigenvalues.
In strong expansions, this quantity becomes very small. In this case, it is suggested to use the critical sound
speed instead:
c* ( , ht )  2
(  1)
ht
(  1)
which does not influence the total flux F and approaches the static speed of sound at the sonic points.
4- The critical sound speed is proportional to the square root of the total enthlapy and is therefore constant
for isoenergenetic Euler flows. That leads to a gain in robustness of the solver
Modelling of high speed flows
Reactive flows
The laws governing chemical reactions have been coupled with the Euler solver. At each iteration of the time
integrator, the Euler code provides and update for (  , w , et ) . The mass fractions are computed by solving a
set of nonlinear algebaric equations in the case of chemical equilibrium and a set of partial differencial equations
in the case of chemical non-equilibrium. The temperature is computed by expressing the identity of the internal
specific energy written in terms of “fluid motion quantities” and “chemical ones”:
ns
ns
ns
1 2
w
2
i 1
i 1
i 1
The outcome is used to update the other thermodynamic quantities ( p, ht , c ) in preparation of the next step.
 YiCv,iT  
Yi hi0
  Yi eivib  et 
The desirable properties that a FVS schemes should possess are:


1- F (W )  F (W )  F (W ) (consistency to the governing equations)

2- F must be continuous, with
F  (W )  F (W ) , M  1
F  (W )  F (W ) , M  1

3- dF / dW must have all eigenvalues  0
dF  / dW must have all eigenvalues  0
Remarks
1- In the definition of the upwind we have a switch at M  1 . The efficiency of the scheme depends on the
definition given to the sound speed.
Property ( 1) is always satisfied independently of the definition of the sound speed.
To ensure that in supersonic regions property (2) is satisfied, we might use the sound speed as it appears in
the true eigenvalues of the system:
cm2  ~
Modelling of high speed flows
p

 ~  ˆ
Reactive flows
2- Some authors prefer to employ a “false sound speed” approach based solely on satisfaction of property (1):
c 2ph  
p

 ~  
whose value is closer to that obtained from available experimental data.
3- Other choices for the definition of the sound speed are a priori possible. However property (3) has to be satisfied.
The consequence is that ~ and  can not vary independently. In particular, it can be shown that ~ must be
greater or equal to  .
4- Contrary to non-reactive flows, the specific heats ratios for reactive flows are no longer constant.
Then, it is propose to perform the limiting procedure also to ˆ and  for second-order accurate computations.
5- The concept of critical sound speed has been extended to reactive flows according to
c* ( , ˆ, ht )  2
Modelling of high speed flows
ˆ (  1)
ht
2  ˆ (  1)
Non reactive flows around a blunt body
A first test of the numerical scheme has been made for a flow around a half cylinder at M   25.
The freestream conditions are taken to approximatively correspond to the standard atmosphere
at the altitude of 75 km ( p  2.52 Pa,   0.43104 kg/m 3 , T  205.3 K ).
1- scheme a: first-order scheme
2- scheme b: second-order scheme with limitation on (  , w, p)
3- scheme c: second-order scheme with limitation on (  , w, ht )
Iso-Cp contours
Computational domain (65x41) nodes
Iso-Pt contours
Modelling of high speed flows
Non reactive flows around a blunt body
Entropy generation along stagnation streamline
Better accuracy of the second
order schemes (b & c).
Sheme c allows better accordance
with respect to entropy generation
and stagnation temperatures
levels
Total enthalpy residual versus iterations
Stagnation temperature versus iterations
Modelling of high speed flows
Reactive flows around a blunt body
The results obtained on the previous case are illustrated now with the hypothesis of reactive flow in thermodynamic
equilibrium. Comparisons with non-reactive flow results show that the shock has considerably moved towards the
body. The temperature level is considerably reduced and the stagnation temperature is now approximatively 22%
of the corresponding non-reactive flow stagnation temperature.
Non reactive
flow
Reactive
flow
Iso-Cp contours
Diatomic
Oxygen
Mass fractions for Oxygen and Nitrogen
Non reactive
flow
Iso-Temperature contours
Modelling of high speed flows
Reactive
flow
Diatomic
Nitrogen
Summary
1- Augmented dissipation schemes have been developed for generalised finite volume schemes.
→
extension to finite elements schemes is automatic
2- Flux Difference Splitting and Flux Vector Splitting schemes may be re-interpreted as
a central scheme with augmented dissipation operator, but are first-order accurate schemes
3- Augmented dissipation has to be activated where discontinuities are present, but has to vanish in regions
where the solution is regular
→
introduction of limiters and shock sensors
4- Numerical schemes have been developed in such a way that total enthalpy is constant for iso-energenetic
steady flows (needed in order to obtain accurate stagnation temperature value)
In particular, the variables on which limitation is performed have to be accurately chosen.
(  , w, p)  (  , w, ht )
5- The treatment of high enthalpy non-reactive and reactive flows has been discussed.
Modelling of high speed flows