SLIDE 1
Analysis of Numerical Shock Instability and Cure by HLLC-FORCE Hybrid Scheme
Li Yuan
Institute of Computational Mathematics and Scientific/Engineering Computing, AMSS, Chinese Academy of Sciences Joint work with Lijun Hu
SGS on ANMCFMRP, May 21-27 2014, Beijing
SLIDE 2 Outline
1.
ntrod
uction ion
- 1. Phenomena of Numerical Shock Instability
- 2. Advances in Research of Numerical Shock Instability
2.
rning ing Equ Equati ation
3.
tabil ility ity Ana Analys lysis is
- 1. Linear Stability Analysis
- 2. Stability Analysis with a Matrix-based Method
4.
ybrid id Sc Scheme heme fo for 2D 2D Eule Euler E r Equ quati ations
- ns
- 1. Construction of Hybrid Scheme
- 2. Switch Function
5.
umeri rical cal Exp Experi erime ments nts 6.
lusio sion an n and d Pro rospe spect ct
SLIDE 3
Phenomena of Numerical Shock Instability
- HLLC(Harten-Lax-Leer-Contact) Riemann solver has been widely used in CFD due to its
ability for capturing shock, rarefaction, and contact discontinuity accurately.
- However, HLLC for multidimensional problems, unphysical phenomena will appear.
Peery and Imlay found the carbuncle phenomenon in two dimensional flows around a blunt body. Quirk reported that if the grid line has a slight perturbation, the phenomenon of odd-even decoupling will arise when some low-dissipative Riemann solvers are used. These phenomena are called Numerical Shock Instability because they
- ften occur near the shock.
carbuncle
SLIDE 4 Advances in Research of Numerical Shock Instability
- Mechanism Analysis
- 1. Quirk believes that the scheme will suffer from odd-even decoupling if the perturbation of
pressure feedbacks to the density.
- 2. Liou thinks that the root of numerical shock instability is because the flux of mass contains
the term of pressure.
- 3. Xu finds that numerical shock instability is caused by two factors: the absence of dissipation
in the tangential of shock, and the effect that pressure imposes on mass flux.
- 4. Dumbser indicates the structure of numerical shock is very imoportant for stability of shock.
- Method of Cure
- 1. Moschetta modifies the eigenvalue related with the shear wave to cure the numerical shock
instability.
- 2. Quirk amd Pandolfi suggest using the dissipative scheme(such as HLL)near the shock but
using the low-dissipative scheme(such as HLLC)in other region.
- 3. Kim (2009) describes a method which combine the HLLC and HLL to cure the instability
phenomena of HLLC.
- 4. Rotated Roe(Davis, Ren) and genuinely multidimensional schemes
However, hybrid methods in 2 & 3 may introduce unnecessary dissipation.
SLIDE 5 Advances in Research of Numerical Shock Instability
- Our Work
- 1. Analyze shock instability of HLLC scheme.
- 2. Construct a hybrid scheme to cure the numerical shock instability.
Compared with the hybrid scheme proposed by Kim (2009), our method is different in: （1）use the FORCE(First-Order Centred deterministic scheme), whose dissipation seems lower than HLL, as the dissipative scheme; （2）blend the HLLC and FORCE only in mass and tangential momentum, and use HLLC in other components of the fluxes; （3）by defining a switch function, we only use the hybrid scheme near the strong shock, and still use the HLLC scheme in the other region.
- 3. Numerical experiments indicate the hybrid scheme can not only cure the shock instability
phenomenon but also retain the high resolution of HLLC.
SLIDE 6
Consider the Euler equations in two dimensions,
（1）
with, Discretize Eq.(1) with the finite volume method over a structured quadrilateral grid （2） where, are numerical fluxes in the x- and y-directions, respectively.
( ) ( )
0. t x y ∂ ∂ ∂ + + = ∂ ∂ ∂ F U G U U
( ) ( ) ( ) ( )
2 2
, , . u v uv u p u uv v p v E u E p v E p ρ ρ ρ ρ ρ ρ ρ ρ ρ + = = = + + +
U F U G U
, i j Ω
d ˆˆ ˆˆ 0. 1 1 1 1 i+ ,j i- ,j i,j+ i,j- d i,j 2 2 2 2 t Ω + − + − =
U F F G G
ˆ ˆ , 1 1 i+ ,j i,j+ 2 2 F G
SLIDE 7
- 3. Linear Stability Analysis
Assume the fluid flows along the X-direction from left to right, and the shock normal is X-direction. The velocity in Y-direction (shock tangent) is zero, and the flow velocity in X-direction is , where is Mach number, is sound speed. The physical variables on is , is constant values, is perturbation values. We use forward Euler to discretize the semi-discretized equations (2)into full-discretized equations: （3） Under the odd-even form of perturbation in one direction, assuming , and using linearization, we can obtain: （4）
The spectral radius of augmented matrix determines the stability: If it is more than 1, the scheme is unstable; If it is equal to 1, the scheme is marginal stable; If it is less than 1, the scheme is stable.
u M a = M a
, i j Ω , , n n i j i j δ = + W W W W , n i j δ W
n+1 n i,j i,j x 1 1 y 1 1 i+ ,j i- ,j i,j+ i,j- 2 2 2 2
ˆˆ ˆˆ σ σ = − − − − U U F F G G
x y
with , . t t x y σ σ ∆ ∆ = = ∆ ∆
n+1 n i,j i,j
δ δ = W A W
A
1 1 2
ˆˆ ( , )
i i i + + =
F F U U
SLIDE 8 Linear Stability Analysis
HLLC and FORCE schemes are analyzed using the linear stability analysis. In order to explain the instability phenomenon, two perturbation models are analyzed.
1.X-direction perturbation. There are odd-even perturbations in X-direction, and the physical quantities are constant in Y-direction.
If the flow is subsonic , the eigenvalues of the augmented matrix are as follows: If the flow is supersonic ,
HLLC and FORCE are all stable if time step satisfies the CFL condition.
(5) (6)
( )
u a <
( )
u a >
SLIDE 9
Linear Stability Analysis
2.Y-direction perturbation. There are odd-even perturbations in Y-direction, and the physical quantities are constant in X-direction.
Because the velocity in Y-direction is zero,the flow is subsonic. The eigenvalues of augmented matrix are as follows:
In this case, the HLLC is marginally stable and the FORCE is stable. We see that two eigenvalues which are equal to 1 in (7) lead to the instability of HLLC scheme.
(7)
SLIDE 10 Stability Analysis with a Matrix-Based Method
We define a uniform Cartesian mesh containing cells. The cells are initialized with a vertical steady shock on the grid line between the 6-th and 7-th column with upstream Mach number . The upstream values are nomalized as follows: The downstream state can be calculated via the Rankine-Hugoniot relations. For the stability analysis of a steady-state field, submitted to small numerical random errors, we are interested in the temporal evolution of these errors. The field is expanded into its steady mean value( )and the error( ). where, cells in the whole domain. Substitute (9) into (2), and do linearization about the mean value,
(8)
11 11 × 7 M =
1
U
(9)
∧ δ
ˆ
q q q
δ = + W W W
11 11 121 q = × =
SLIDE 11 Stability Analysis with a Matrix-Based Method
after some lengthy complex computations, we can get the error evolution
- f all cells in the whole domain:
Considering only the evaluation of initial errors, the solution of time-linear invariant system(10)is: and remains bounded if , where are the eigenvalues of , denote the real part of .
(10)
11 11 121 q = × = S
(11)
( )
λ S
( )
( )
Re λ S
( )
( )
( )
max Re λ < S
( )
λ S
SLIDE 12
Stability Analysis with a Matrix-Based Method
Fig.1. Distribution of the eigenvalues of S in the complex plane for HLLC scheme.
max(Re)=0.7414, unstable.
SLIDE 13
Stability Analysis with a Matrix-Based Method
Fig.2. Distribution of the eigenvalues of S in the complex plane for FORCE scheme.
max(Re)=-0.1063, stable.
SLIDE 14
Stability Analysis with a Matrix-Based Method
In order to comply with the linear stability analysis, we design two hybrid schemes. First, we use HLLC to calculate X-direction flux and FORCE to calculate Y-direction flux (left, Fig.3). Then, we change the sequence (right, Fig.3).
Fig.3. Distribution of eigenvalues in complex plane. left: max(Real) = -0.0935, stable;
right: max(Re) = 0.5267, unstable.
Once again, we verify that the instability of HLLC scheme will occur if it is applied in the Y-direction (tangential to shock surface).
SLIDE 15
- 4. Construction of Hybrid Scheme
From linear stability analysis, it is found HLLC has insufficient dissipation to suppress perturbation in the equations of density and velocity tangential to the cell interface. To cure it, the HLLC flux and the FORCE flux are combined in grid-aligned framework. Take the fluxes of X-direction for example, the hybrid scheme can be written as: The weight coefficient can be chosen as: where, is the shock normal direction, and is the normal direction of the cell interface.
when shock normal direction is aligned with the interface normal direction, the hybrid scheme is HLLC; when shock normal is perpendicular to interface normal, the weight of FORCE scheme reaches the maximum, 0.5.
(12) (13)
s
n
n
SLIDE 16 Construction of Hybrid Scheme
The shock normal direction is unknown in calculation, but one can approximate it with the velocity- difference vector between two adjacent grids:
(14)
4
where 10 U ε
−
=
SLIDE 17 Switch Function
We define a function to identify the location of shock waves( represents strong shock ): where . The shock instability occurred due to the perturbation in the transverse direction
- f the shock. Therefore, all the neighboring cells must be examined. The required surfaces for the
function in the x-and y-direction are calculated as: Now, we define the switch function :
(15)
p
f
16
10 ε
− ∗ =
(16) (17)
1.0
p
f =
p
f f
SLIDE 18
Switch Function
Finally, take the X-direction for example, the flux can be calculated as:
(18)
SLIDE 19
MUSCL Reconstruction
In the calculation, the high order accuracy can be obtained by MUSCL reconstruction: The limiter function is a minmod function, and the differences of characteristic variables are:
(19) (20)
SLIDE 20
- 5. Numerical Experiments
- 1. Odd–even decoupling
Quirk's odd–even decoupling test problem. The centerline of the grid is perturbed from a perfectly uniform grid by ±10-6. A moving shock of Ms = 6.0 propagates down a duct.
We can see clearly hybrid scheme elimilates the instability phenomenon
SLIDE 21 Numerical Experiments
We turn the odd–even decoupling test problem into another version. The perturbation of centerline is elimilated, but on every grid we add random perturbation of 10-6 in the initial conditions.
We can see clearly hybrid scheme elimilate the instability phenomenon
SLIDE 22 Numerical Experiments
- 3. The diffraction of a supersonic shock moving over a 90o corner
The shock diffraction problem is another test for which many Godunov-type schemes are known to fail. The shock Mach number is 5.09 in this problem. The domain is a unit square [0,1]×[0,1], which is discretized into a 400×400 uniform grids.
It is obvious that HLLC scheme suffers from shock instability. Our hybrid scheme does not have such flaw.
SLIDE 23 Numerical Experiments
- 4. Mach 10 hypersonic flow over a cylinder
This is another well-known test to examine the catastrophic carbuncle failings of upwind
- schemes. The numerical simulation of a Mach number 10 inviscid flow around a circular
- cylinder. In this test problem, 40×160 grids are used.
We can see that the result of HLLC scheme is not correct, but the hybrid scheme have good shock picture.
SLIDE 24 Numerical Experiments
- 5. Forward facing step problem
The computational domain is [0,3]×[0,1]. In this test problem, 240×80 grids are used.
Even in this case, the HLLC solver exhibits instability in regions after the normal shock. However, our hybrid scheme works very well.
SLIDE 25 Numerical Experiments
This test problem is proposed by Lax and Liu. Computational region is [0,1]×[0,1] and uses the grids of 400×400. The initial condition is:
When we use the low-dissipative schemes(such as HLLC and Roe ), two pseudo-shocks will appear:
SLIDE 26 Numerical Experiments
The hybrid scheme eliminates the pseudo-shocks, and retain the high resolution of HLLC.
SLIDE 27
- 6. Conclusion and Prospect
- Conclusion
- 1. The use of HLLC scheme in the tangential direction of the shock causes the
numerical shock instability (A little bit more specific, the marginal stability of mass and momentum tangential to cell interface is the root of instability);
- 2. The present hybrid scheme can not only eliminate the numerical shock
instability of HLLC but also retain the high resolution of HLLC.
- Prospect
- 1. The analysis methods and hybrid scheme can implemented to the other scheme
(such as Roe);
- 2. The methods can be implemented on irregular grids(such as triangular grids).
