SGS on ANMCFMRP, May 21-27 2014, Beijing 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
Outline 1. 1. Intr ntrod oduct uction ion 1. Phenomena of Numerical Shock Instability 2. Advances in Research of Numerical Shock Instability 2. 2. Gove overn rning ing Equ Equati ation ons 3. 3. Stab tabil ility ity Ana Analys lysis is 1. Linear Stability Analysis 2. Stability Analysis with a Matrix-based Method 4. 4. Hybr ybrid id Sc Scheme heme fo for 2D 2D Eule Euler E r Equ quati ations ons 1. Construction of Hybrid Scheme 2. Switch Function 5. 5. Nume umeri rical cal Exp Experi erime ments nts 6. 6. Conc onclu lusio sion an n and d Pro rospe spect ct
1. Introduction 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 often occur near the shock. carbuncle odd-even decoupling
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.
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.
2. Governing Equations Consider the Euler equations in two dimensions, ( ) ( ) ∂ ∂ ∂ F U G U U + + = 0. （ 1 ） ∂ ∂ ∂ t x y with, ρ ρ ρ u v ρ 2 ρ + ρ uv u ( ) u p ( ) = = = U , F U , G U . 2 ρ ρ + ρ v uv v p ( ) ( ) + + E u E p v E p Discretize Eq.(1) with the finite volume method over a structured quadrilateral grid Ω i j , d U Ω + − + − = （ 2 ） ˆˆ ˆˆ F F G G 0. 1 1 1 1 d t i+ ,j i- ,j i,j+ i,j- i,j 2 2 2 2 where, are numerical fluxes in the x- and y-directions, respectively. ˆ ˆ F , G 1 1 i+ ,j i,j+ 2 2
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 , = u M a 0 0 0 where is Mach number, is sound speed. a M 0 0 The physical variables on is , is constant values, is perturbation values. Ω n n = δ + n δ W W W W W i j , i j , i j , 0 i j , 0 We use forward Euler to discretize the semi-discretized equations (2)into full-discretized equations: = − σ − − σ − ˆˆ ˆˆ n+1 n U U F F G G i,j i,j x 1 1 y 1 1 i+ ,j i- ,j i,j+ i,j- （ 3 ） 2 2 2 2 ∆ ∆ t t σ = σ = with , . ∆ ∆ x y x y Under the odd-even form of perturbation in one direction, assuming , and using + = ˆˆ F F U U ( , ) + 1 i i 1 i linearization, we can obtain: 2 （ 4 ） δ = δ n+1 n W A W i,j i,j The spectral radius of augmented matrix determines the stability: A 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.
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: < u a 0 0 (5) If the flow is supersonic , ( ) > u a 0 0 (6) HLLC and FORCE are all stable if time step satisfies the CFL condition.
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: (7) 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.
Stability Analysis with a Matrix-Based Method We define a uniform Cartesian mesh containing cells. The cells are initialized with a vertical × 11 11 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: M = 7 0 (8) The downstream state can be calculated via the Rankine-Hugoniot relations. U 1 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( ). ∧ δ = + δ ˆ W W W (9) q q q where, cells in the whole domain. Substitute (9) into (2), and do linearization about the q = × = 11 11 121 mean value,
Stability Analysis with a Matrix-Based Method after some lengthy complex computations, we can get the error evolution of all cells in the whole domain: q = × = 11 11 121 (10) Considering only the evaluation of initial errors, the solution of time-linear invariant system(10)is: (11) ( ) ( ) ( ) ( ) ( ) ( ) λ < and remains bounded if , where are the eigenvalues of , denote the λ S Re λ S S max Re 0 S ( ) real part of . λ S
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.
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.
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).
Recommend
More recommend
