Use of Jordan forms for convection-pressure split Euler solvers

Use of Jordan forms for convection-pressure split Euler solvers arXiv:1607.00947v5 [math.NA] 21 Sep 2017 Naveen Kumar Garga,1 , S.V. Raghurama Raob , M. Sekharc a Department of Mathematical Sciences, Indian Institute of Science (IISc), Bangalore, India b Department of Aerospace Engineering, IISc, Bangalore, India c Department of Civil Engineering, IISc, Bangalore, India Abstract In this study, we analyze convection-pressure split Euler flux functions which contain genuine weakly hyperbolic convection subsystems. A system is said to be a genuine weakly hyperbolic if all eigenvalues are real with no complete set of linearly independent (LI) eigenvectors. To construct an upwind solver based on flux difference splitting (FDS) framework, we require to generate complete set of LI eigenvectors. This can be done through addition of generalized eigenvectors which can be computed from theory of Jordan canonical forms. Once we have complete set of LI generalized eigenvectors, we construct upwind solvers in convection-pressure splitting framework. Since generalized eigenvectors are not unique, we take extra care to ensure no direct contribution of generalized eigenvectors in the final formulation of both the newly developed numerical schemes. First scheme is based on Zha and Bilgen type splitting approach, while second is based on Toro & Vázquez splitting. Both the schemes are tested on several bench-mark test problems on 1-D and one of them is tested on some typical 2-D test problems which involve shock instabilities. The concept of generalized eigenvector based on Jordan forms is found to be useful in dealing with the genuine weakly hyperbolic parts of the considered Euler systems. Keywords: Convection-pressure splittings, Jordan forms, Upwind schemes 1. Introduction Numerical algorithms based on compressible Euler systems are of great importance and are frequently used in simulations. These algorithms can be broadly divided into two major categories, namely, central discretization methods and upwind discretization methods. In this study we mainly focus on upwind methods and, in particular, on the Flux Difference Splitting (FDS) based upwind schemes. Another class of upwind methods based on Flux Vector Splitting (FVS) schemes, like Steger & Warming [22] and van Leer [27] schemes, are constructed using eigenvector structure of Euler system. But these schemes are quite diffusive and can’t capture isolated contact discontinuities crisply. Similarly, schemes based on FDS framework, such as an approximate Riemann solvers of Roe [20] and Osher [16], also depend on the eigen-structure (both eigenvalues Email addresses: garg.naveen70@gmail.com, naveen@tifrbng.res.in (Naveen Kumar Garg), raghu@aero.iisc.ernet.in (S.V. Raghurama Rao), muddu@civil.iisc.ernet.in (M. Sekhar) 1 Current address: Post Doctoral Fellow, TIFR Center for Applicable Mathematics, Bangalore, India Preprint submitted to Elsevier September 22, 2017 and eigenvectors) but are quite accurate. Osher scheme captures expansion waves well but is computationally expensive. Roe’s approximate Riemann solver is accurate and can capture steady discontinuities either exactly or with a single interior point, without being as expensive. Because of less numerical diffusion, Roe scheme tends to produce unphysical expansion shocks [29], post-shock oscillations [23] and it is non-trivial to avoid instability problems [18] in its application. Similarly, HLL scheme [5] and HLLC scheme [25] are special upwind schemes mainly dependent on structure of eigenvalues. HLLC scheme is a modified version of HLL scheme and can capture an isolated contact discontinuity exactly but suffers from infamous carbuncle phenomena and some other shock instabilities in 2-D, as shown in [6], [21]. There is a different class of upwind schemes based on splitting of the Euler flux function, with several possible splittings. We consider here the popular splittings of the Euler flux function as is usually done in three distinct ways. In the first category, we have the AUSM family of schemes [13], [14] and CUSP scheme of Jameson [8], in which flux vector is split into a convection and a pressure part. Liou and Steffen first introduced the category of AUSM family of schemes, such that in their convectionpressure splitting, the pressure term is alone present in the pressure part of the split flux vector. The Second category of convection-pressure split schemes were initially proposed by Steger and Warming in [22], but were throughly explored by Jameson [9], Zha & Bilgen [32], Balakrishnan & Deshpande [2] and by Raghurama Rao & Deshpande [19]. The structure of this splitting is such that pressure term of momentum equation and the term containing the product of pressure and velocity in the energy equation constitute the pressure part of the convection-pressure split Euler flux function. The main feature of this type of convection-pressure splitting is that the eigenvalues corresponding to Jacobian of pressure subsystem become completely free from fluid velocity u, unlike in the splitting utilized by Meng-Sing Liou and others. Another interesting and more recent convection-pressure splitting is proposed by Toro & Vázquez [26]. In this category, the authors split convection and pressure parts in such a way that convection flux becomes completely free from pressure terms. For all the three splittings, convection part always corresponds to a genuine weakly hyperbolic system. As the each convection part is weakly hyperbolic, we can utilize the theory of Jordan Canonical forms to recover complete set of linearly independent (LI) generalized eigenvectors. In contrast, each pressure part corresponds to strict or non-strict hyperbolic system. In the first category of convection-pressure splitting, although the convection part contains contribution of two different eigenvalues namely, u and γu, with u as repeated eigenvalue, the eigenvalues of pressure part do not contain any contribution of acoustic speed and this may result in an unstable scheme if used in FDS framework. As pressure part of other two splittings contain the contribution of acoustic speeds as well, we propose two numerical schemes namely, Zha and Bilgen split - flux difference splitting scheme (ZBS-FDS) and Toro and Vázquez split - flux difference splitting TVS-FDS scheme. The main idea of each scheme is first to construct traditional FDS scheme for pressure sub-system and then utilize the resulting averaged values of all variables together with theory of Jordan forms for convection part. Our motivation is to develop efficient and workable flux difference split schemes based on convection-pressure splitting, together with the use of Jordan forms for convection subsystems, rather than focusing on reproducing an ideal approximate Riemann solver in this framework. Both schemes are tested on various shock tube problems in 1-D and are found to require no entropy fix for sonic points and in strong expansion regions. Both the schemes capture isolated and steady contact discontinuities exactly. Out of the two, ZBS-FDS scheme is extended to 2 two dimensions and is further tested on a variety of shock instability problems, including shock diffraction around a corner, flow over over a half cylinder and reflection of a plane shock from a wedge. 2. Convection-Pressure splittings for Euler flux function Consider the one-dimensional inviscid Euler system ∂U ∂F (U ) + = 0, ∂t ∂x (x, t) ∈ IR × [0, ∞), (1) where U : (IR × IR+ ) 7−→ Ω ⊆ IR3 , Ω is an open subset, is the conserved variable vector and F : Ω 7−→ IR3 is the flux vector defined by     ρ ρu     U =  ρu  and F (U ) =  p + ρu2  (2) pu + ρuE ρE Here the total energy E is defined as the sum of internal energy (e) and kinetic energy p ( 12 u2 ), given as: E = e + 12 u2 = ρ(γ−1) + 12 u2 . Till now, three distinct convection-pressure splittings have been proposed, which are described in the following. 2.1. Liou and Steffen splitting procedure Liou and Steffen, in formulating their upwind scheme [14], introduced a unique convection-pressure splitting by taking out pressure part from momentum equation of full Euler system. F = F LS + F LS (3) c p where   F LS = c ρu   ρu2  and F LS p ρuE + pu   0   =  p 0 (4) Let us split the system (1) into Liou and Steffen type convection and pressure subsystems, for gaining better insight by analyzing each part separately. ∂U ∂F LS c (U ) + = 0 ∂t ∂x (5) and ∂F LS ∂U p (U ) + = 0 ∂t ∂x Both subsystems can also be written in quasilinear form as follows. (6) ∂U ∂U + ALS = 0 c ∂t ∂x (7) ∂U ∂U + ALS = 0 p ∂t ∂x (8) 3 Here ALS and ALS are Jacobian matrices for convection and pressure parts respectively c p and are given by   0 1 0   −u2 2u 0 ALS = c γE − 23 (γ − 1)u2 −γuE + (γ − 1)u3 and  0 0  ALS =  12 (γ − 1)u2 p 0 −(γ − 1)u 0 0 γu   (γ − 1) 0 LS Eigenvalues corresponding to convective Jacobian matrix ALS are λLS c c,1 = γu, λc,2 = LS λc,3 = u and algebraic multiplicity (AM) of the eigenvalue u is 2. Similarly, eigenvalues LS LS corresponding to pressure Jacobian matrix ALS are λLS p p,1 = −(γ − 1)u, λp,2 = λp,3 = 0. Since AM of u is 2, so we have to find its eigenvector space to see whether ALS has c complete set of linearly independent eigenvectors or not. The analysis of matrix ALS c shows that convective subsystem is weakly hyperbolic as there is no complete set of linearly independent eigenvectors. Indeed, its eigenvectors are     1 0     LS (9) RLS c,1 =  0 and Rc,2 =  u  1 2 1 2u LS LS Similarly, eigenvectors corresponding to λLS p,1 = −(γ − 1)u, λp,2 = λp,3 = 0 of pressure subsystems are:       1 0 0       LS LS LS 0 1 (10) Rp,1 =   , Rp,2 =   , Rp,3 =  1  1 2 0 u −2u Convection subsystem turns out to be weakly hyperbolic, and Jordan theory can be applied to explore it further, whereas pressure subsystem is non-strict hyperbolic. Apart from eigenvectors, traditional FDS solvers depend heavily on eigenvalues also, but for present case all eigenvalues are either u or constant times u. In other words, there is no direct or indirect contribution of acoustic speed a as an eigenvalue for both subsystems. This is a serious issue as u frequently goes to zero or near to zero in a flow field which results in zero or near zero diffusion at some parts of the flow. Thus, the scheme may generate near zero diffusion which effectively reduces the scheme to forward in time and central in space (FTCS) framework, and as FTCS doesn’t preserve the hyperbolicity, the solution blows-up. In fact, we constructed FDS scheme for present splitting but unfortunately, it led to blow-up of the solution for almost all problems. Note that we are only considering the application of flux difference splitting to Liou and Steffen splitting here and not their alternative upwinding procedure. 2.2. Zha and Bilgen splitting procedure Another type of flux splitting is given by Zha and Bilgen [32], in which they split the full Euler flux function into convection and pressure fluxes in such a way that eigenvalues 4 corresponding to Jacobian of pressure flux AZB contains no contribution of fluid velocity p u, unlike in Liou and Steffen splitting. Their convection-pressure splitting is as follows. F = F ZB + F ZB c p where  ρu  (11)  0      F ZB =  ρu2  and F ZB =p c p pu ρuE (12) As done earlier, we split system (1) into convection and pressure subsystems, using Zha and Bilgen type flux splitting, separately as ∂U ∂F ZB (U ) c + = 0 ∂t ∂x (13) and ∂F ZB (U ) ∂U p + = 0 ∂t ∂x Again, both subsystems can also be written in quasilinear form as follows. (14) ∂U ∂U + AZB = 0 c ∂t ∂x (15) ∂U ∂U + AZB = 0 p ∂t ∂x (16) where, AZB and AZB are Jacobian matrices for convection and pressure parts respecc p tively and are given by   0 1 0   2u 0 AZB =  −u2 c −uE E u and   AZB = p 0 0 1 2 2 (γ − 1)u 2 3 − aγu + (γ−1) 2 u −(γ − 1)u a2 γ − (γ − 1)u2 0  (γ − 1)   (γ − 1)u ZB Now, eigenvalues corresponding to convective Jacobian matrix AZB are λZB c c,1 = λc,2 = ZB λc,3 = u, thus algebraic multiplicity (AM) of eigenvalue u is 3. Similarly, eigenvalues q (γ−1) ZB λZB corresponding to pressure Jacobian matrix AZB are λ = − p p,2 = p,1 γ a, q (γ−1) 0 and λZB p,3 = γ a. Since AM of u is 3, so we have to find its eigenvector space to see whether AZB has complete set of linearly independent eigenvectors or not. The c analysis of matrix AZB shows that convective subsystem is weakly hyperbolic as there c is no complete set of linearly independent eigenvectors. Indeed, its eigenvectors are     0 1     ZB ZB Rc,1 =  u and Rc,2 =  0 (17) 0 1 5 Since all eigenvalues for pressure part are real and distinct, this makes pressure subsystem strictly hyperbolic. Analysis of the flux Jacobian matrix for the pressure part shows complete set of eigenvectors, as given below.       0 0 1       ZB 1 1   , RZB   u RZB = (18) , R =   p,1 =  p,2 p,3    a a 1 2 √ √ u− u+ u γ(γ−1) γ(γ−1) 2 2.3. Toro and Vázquez splitting Procedure More recently, Toro & Vázquez-Cendón [26] presented a flux splitting in which convection part contains no pressure term at all, leading to following type of splitting. V V F = FT + FT c p where V FT c (19)     0 ρu     V =  ρu2  and F T = p  p γ 1 3 γ−1 pu 2 ρu (20) Let us split the system (1) into convection and pressure subsystems, for gaining better insight by analysing each part separately. V ∂U ∂F T (U ) c + = 0 ∂t ∂x (21) and V ∂F T (U ) ∂U p + = 0 ∂t ∂x Again, both subsystems can also be written in quasilinear form as follows. (22) ∂U V ∂U + AT = 0 c ∂t ∂x (23) ∂U V ∂U + AT = 0 p ∂t ∂x (24) V V Here AT and AT are Jacobian matrices for convection and pressure parts respectively c p and are given by   0 1 0   V 2u 0 AT =  −u2 c 3 2 0 −u3 2u and   V AT = p 0 0 1 2 2 (γ − 1)u 2 ua − (γ−1) + 21 γu3 −(γ − 1)u 2 a (γ−1) − γu2 0  (γ − 1)  γu V V TV Eigenvalues corresponding to convective Jacobian matrix AT are λT c c,1 = 0, λc,2 = TV λc,3 = u and algebraic multiplicity (AM) of eigenvalue u is 2, so we have to find V its eigenvector space to see whether AT has complete set of linearly independent c 6 V eigenvectors or not. The analysis of matrix AT shows that convective subsystem c is weakly hyperbolic as there is no complete set of linearly independent eigenvectors. Indeed, its eigenvectors are     1 0     TV V RT (25) c,1 =  0 and Rc,2 =  u  1 2 1 2u V Similarly, the eigenvalues corresponding to pressure Jacobian matrix AT p , when 1 1 V TV TV evaluated, are found to be λT p,1 = 2 (u − β), λp,2 = 0, λp,3 = 2 (u + β), where β √ = u2 + 4a2 . All eigenvalues for pressure part are real and distinct and this makes pressure subsystem strictly hyperbolic. Analysis of the flux Jacobian matrix for the pressure part shows complete set of eigenvectors, as given below.       0 0 1       V TV TV 1 1 RT (26)  , Rp,2 =  u  , Rp,3 =   p,1 =  1 u+β 1 2 u + 21 ( u−β ) ( ) u + u γ−1 2 γ−1 2 Since the convective subsystems for both Zha-Bilgen splitting and Toro-Vázquez splitting have incomplete set of linearly independent (LI) eigenvectors, it will be nontrivial to construct an upwind scheme based on eigenvector structure. 3. Addition of generalized eigenvectors V First we consider Jacobian matrix AT corresponding to Toro-Vázquez convective c subsystem, for which, our aim is to get complete set of linearly independent generalized eigenvectors. For this system, we have two different sets of eigenvalues. Here, we briefly discuss a procedure to find generalized eigenvectors for cases where resultant Jordan matrix possess exactly one Jordan block for each set of eigenvalues. Let    J =   J (λ1 ) 0 .. . 0 J (λ2 ) .. . ··· ··· .. . 0 0 .. . 0 0 ··· J (λp )     , where λ1 , λ2 , · · · , λp ∈ σ(A)   are set of distinct eigenvalues, some or all of them with arithmetic multiplicity greater than one. Moreover, assume there exists a single Jordan block for each λi . Let us focus on one such λi , with AM equal to m > 1. Then     J (λi ) =    λi 1 .. . .. . .. .       1 λi m×m In order to find set of generalized eigenvectors corresponding to λi , we need to focus 7   on portion P ∗ = X 1 , X 2 , X 3 , ........., X m of P = [...P ∗ ...] that corresponds to the position J (λi ) in J . Now AP = P J implies AP ∗ = P ∗ J (λi ), i.e.,   λi 1   .. ..  . .       A X 1 , X 2 , X 3 , ........., X m = X 1 , X 2 , X 3 , ........., X m     ..  . 1 λi m×m On equating columns on both sides, we get AX 1 = λi X 1 AX 2 = λi X 2 + X 1 AX 3 = λi X 3 + X 2 .. . (27) AX m = λi X m + X m−1 V Now u is a repeated eigenvalue of matrix AT with AM is equal to two and other c eigenvalue is zero with multiplicity one. First we have to compute ranks of matrices V V V (AT − uI), (AT − uI)2 , · · · . It turns out that rank(AT − uI)2 = 1 = c c c TV rank(Ac − uI)3 , which means there should be a Jordan block of order 2. Therefore, there is single Jordan block of order two corresponding to an eigenvalue u. Thus, a V , i.e., Jordan chain of order two will be formed by matrix AT c V AT X 1 = uX 1 and c (28) V AT X 2 = uX 2 + X 1 c should hold. From first relation we get V X 1 = RT c,2   1   = u  1 2 2u (29) and on using X 1 in the second relation of (28), we can find required generalized eigenvector X 2 which is given below.   x1   V (30) X 2 = RT c,3 =  1 + ux1  1 2 u + 2 u x1 Here, x1 ∈ IR is a real constant and det(P ) is equal to one. If we take P equal to     0 1 x1 0 0 0     V u 1 + ux1  then P −1 AT u 1 = J 1 (31) P = 0 0 c 1 2 1 2 0 0 u 1 u + 2 u x1 2u and if we take P equal to  1 x1  u 1 + ux1  1 2 2u u+ 1 2 2 u x1   0 u   V 0 then P −1 AT P = 0 c 0 1 8  1 0 u  0 = J 2 0 0 (32) Next, we have to find generalized eigenvectors corresponding to Jacobian matrix AZB for c a convective subsystem of Zha and Bilgen type splitting. As explained earlier, eigenvalues for AZB are u, u, u and set of LI eigenvectors are, c     0 1     ZB u (33) and R = RZB =  0   c,2 c,1 1 0 On computing ranks of matrices (AZB − uI), (AZB − uI)2 , · · · , we find rank(AZB − c c c ZB uI)2 = 0 = rank(Ac − uI)3 . Thus there will be one Jordan block of order 2 and since all eigenvalues are equal then there must be another Jordan block of order one. In short, Jordan matrix J corresponding to matrix AZB is made up of two Jordan blocks, c which clearly shows that for the present case there is no single Jordan block for given set of eigenvalues. Thus, earlier theory may not be directly applicable for this case. But a Jordan chain of order two should form by matrix AZB . If possible, without loss of c generality, let first assume ZB AZB RZB c c,1 = uRc,1 AZB X = uX + RZB c c,1 (34) holds. On expanding second relation         0 1 0 x1 x1 1         2 2u 0   x2  = u  x2  +  u  −u −uE E u x3 x3 0 and from first two equations, we get x2 = ux1 + 1 (35) −uEx1 + Ex2 + ux3 = ux3 (36) Similarly, from third equation ⇒ x2 = ux1 . We get two different expressions for real constant x2 , which is a contradiction. Thus ZB eigenvector Rc,1 can’t form a Jordan chain of order two corresponding to matrix AZB . c If possible, let us assume now RZB forms a Jordan chain of order two, i.e., c,2 ZB AZB RZB c c,2 = uRc,2 AZB X = uX + RZB c c,2 (37) holds. Again after expanding second relation, we have         0 1 0 x1 x1 0         2 −u 2u 0 x x = u +   2   2  0 1 −uE E u x3 x3 and further on solving first two equations of expanded system we get x2 = ux1 . 9 (38) From third equation we have −uEx1 + Ex2 + ux3 = ux3 + 1 1 + uEx1 ⇒ x2 = E (39) If we compare both values of x2 , we get 0 = 1 which is impossible, hence a contradiction. ZB Therefore, neither RZB c,1 nor Rc,2 helps in forming a Jordan chain of order two corresponding to matrix AZB . Thus, we need to go more deep into the theory of Jordan c canonical forms to obtain proper generalized eigenvectors. Since there will be a Jordan block of order 2, this means that we need to construct a generalized eigenvector which should help in generating a Jordan chain. Let R(A) denote the space spanned by the columns of matrix AZB − uI. Then c R(A) = x1 A1 + x2 A2 + x3 A3 where, A1 , A2 , A3 are column vectors of AZB − uI. Now c      −u 1      R(A) = x1  −u2  + x2  u  + x3  E −uE or (40)  0  0 0      1 1 1       R(A) = −ux1  u  + x2  u  = (−ux1 + x2 )  u  E E E  (41) (42) Therefore, column vector X = (1, u, E)t is a range space of R(A). Next, let N (AX) be a null space of column vectors (AZB − uI)X. By definition c  N (AX) = v ∈ IRn ; (AX)v = 0 (43) and dimension of v is equal to number of entries in each row of matrix AX. Now       0 1 −u 1 0       ZB 2 u 0  u  =  0 . (44) AX = (Ac − uI)X =  −u 0 E −uE E 0 For present case, AX is just a null vector, therefore v reduces to a scalar coefficient. By definition of null space of column vectors we have,     0 0     (45)  0 v =  0 0 0 which holds for any v ∈ IR. And by definition, Xv which is equal to Xv should − uI) ∩ N (AZB − uI). Thus, X 1 = (1, u, E)t should be a from a basis for R(AZB c c generalized eigenvector and to check that we need to see whether for X = X 1 , relation AZB X = uX holds or not. Now c       0 1 0 1 1       2 −u 2u 0 u u AZB X = (46) = u       = uX 1 1 c −uE E u E E 10 This generalized eigenvector is expected to form a Jordan chain of order two corresponding to matrix AZB , i.e., c AZB X 1 = uX 1 c (47) AZB X 2 = uX 2 + X 1 c Eigenvector X 2 can be find from second relation and in expanded from it is written as,        x1 0 1 0 x1 1        2 x x −u 2u 0 u (48) = u +  2     2  x3 −uE E u x3 E After little algebra, X 2 comes out as  x1    X 2 =  1 + ux1  = RZB c,2 x3 where x1 , x3 ∈ IR.  (49)  1   X 1 =  u  = RZB c,1 E (50) and we can take RZB c,3 either equal to     1 0      0 or  u . 0 1 If we take    0 1    =  0 then P =  u E 1 RZB c,3 (51) x1 0 1 + ux1 x3 and det(P ) = 1. Further,  u  −1 ZB P Ac P =  0 0 Similarly, if we take  0 1  P = 0 u 1 E 1 u 0   0 1 (52)  0  0 = J 1 u   (53)  u 0 0   1 + ux1  then P −1 AZB P = 0 c 0 x3 u 0  1 = J 2 u x1 (54) Similarly, let  then, det(P ) = −E 6= 0. 1  P =  u E x1 1 + ux1 x3 11 1   u 0 (55) 4. Formulation of ZBS-FDS and TVS-FDS schemes 4.1. ZBS-FDS scheme We first consider pressure subsystem of Zha and Bilgen type splitting and on comparing (14) and (16), we get dF ZB = AZB (56) p p dU The finite difference analogue of the above differential relation is, ZB ∆F ZB = Āp p ∆U (57) ZB ZB ZB where Āp is now a function of left and right states, i.e., Āp = Āp (U L , U R ). Since ∆U is a column vector, therefore it can be written as linear combination of LI eigenvectors. 3 X ZB ZB ᾱp,i R̄p,i (58) ∆U = i=1 On using above expression in (57), ZB ∆F ZB = Āp p 3 X ZB ZB ᾱp,i R̄p,i (59) i=1 or ZB ZB ∆F ZB = ᾱp,1 Āp p ZB ZB ZB R̄p,1 + ᾱp,2 Āp ZB ZB ZB R̄p,2 + ᾱp,3 Āp ZB R̄p,3 (60) ZB (61) which is further equal to ZB ZB ZB ZB ZB ZB ZB ZB ∆F ZB = ᾱp,1 λ̄p,1 R̄p,1 + ᾱp,2 λ̄p,2 R̄p,2 + ᾱp,3 λ̄p,3 R̄p,3 p Now, ∆F +ZB must have contribution of positive part of eigenvalues only, i.e., p ZB ZB ZB (62) ZB ZB ZB (63) ZB +ZB ZB +ZB ZB +ZB ∆F +ZB = ᾱp,1 λ̄p,1 R̄p,1 + ᾱp,2 λ̄p,2 R̄p,2 + ᾱp,3 λ̄p,3 R̄p,3 p Similarly, ZB −ZB ZB −ZB ZB −ZB λ̄p,3 R̄p,3 ∆F −ZB = ᾱp,1 λ̄p,1 R̄p,1 + ᾱp,2 λ̄p,2 R̄p,2 + ᾱp,3 p We now define the standard Courant splitting for the eigenvalues as λ̄± p,i = λ̄p,i ± |λ̄p,i | 2 (64) On using standard upwinding along with above definition, we finally get ∆F +ZB − ∆F −ZB = p p 3 X ZB ZB λ̄ZB ᾱp,i p,i R̄p,i (65) i=1 To determine right side of (65) completely, we need to find average values of eigenvalues along with coefficients which are attached with LI eigenvectors. First we consider ZB linearization equation, ∆F ZB = Āp ∆U , of pressure subsystem of ZBS-FDS scheme. p In expanded form it can be written as      0 0 0 0 ∆(ρ)   1    2 −(γ − 1)ū (γ − 1)  (66)  ∆(p)  =  2 (γ − 1)ū   ∆(ρu)  2 2 (γ−1) 3 ūā ā 2 ∆(pu) ∆(ρE) − γ + 2 ū (γ − 1)ū γ − (γ − 1)ū 12 From the second equation, we get ∆p = 1 (γ − 1)ū2 ∆ρ − (γ − 1)ū∆(ρu) + (γ − 1)∆(ρE) 2 (67) p 1 (γ − 1)ū2 ∆ρ − (γ − 1)ū∆(ρu) + (γ − 1)∆( ) 2 γ−1 1 + (γ − 1)∆(ρu2 ) 2 (68) or ∆p = It further reduces to ū2 ∆(ρ) − 2ū∆(ρu) + ∆(ρu2 ) = 0 which gives average value for conserved variable ū as given below. √ √ ρL u L + ρR u R ū = √ √ ρL + ρR (69) (70) Other root is being neglected as it contains negative sign in the denominator. Let us consider the relation ∆(ρu) = ρ̄∆u + ū∆ρ (71) in expanded form it is written as, ρR uR − ρL uL = ρ̄(uR − uL ) + ū(ρR − ρL ) √ On using ū in the above relation we get ρ̄ = ρL ρR . Similarly, third equation can be written as   ūā2 (γ − 1) 3 ∆(pu) = − + ū ∆ρ γ 2   2 ā 2 − (γ − 1)ū ∆(ρu) + (γ − 1)ū∆(ρE) + γ On expanding further, above equation looks like    2  (γ − 1) 3 ā ūā2 2 + ū ∆ρ + − (γ − 1)ū ∆(ρu) ∆(pu) = − γ 2 γ 1 + ū∆(p) + (γ − 1)ū∆(ρu2 ) 2 (72) (73) (74) We next use (69) and (71) in above equation and after cancellations of some terms we get, ∆(pu) = ū∆p + p̄∆u (75) On using p̄ = ā2 ρ̄ , we have γ ∆(a2 ρu) − ū∆(a2 ρ) = ā2 ρ̄∆u (76) Let η be any flow variable. Consider the relation ∆(ρuη) − ū∆(ρη) = ρ̄η̄∆u 13 (77) which can be written as √ ∆(ρuη) − ū∆(ρη) = ρL ρR ( √ √ ρ L η L + ρR η R )∆u √ √ ρ L + ρR (78) We put η = a2 , then above relation becomes an equation and further on comparing it with (76), we get average value of acoustic speed from ā 2 = √ √ ρL a2L + ρR a2R √ √ ρL + ρ R (79) Similarly, we find coefficients which are attached with LI eigenvector of pressure subsystem for ZBS-FDS scheme as follows. ∆U = 3 X ZB ZB ᾱp,i R̄p,i (80) i=1 In expanded form,  ∆(ρ)   0   1   0         ZB  ZB ZB 1 1  + ᾱp,2   ū  + ᾱp,3   ∆(ρu)  = ᾱp,1      ā ā 1 2 √ √ ū − ū + ∆(ρE) ū γ(γ−1) γ(γ−1) 2 (81) On comparing first equation we have, ZB ᾱp,2 = ∆ρ (82) ZB ZB ZB ∆(ρu) = ᾱp,1 + ūᾱp,2 + ᾱp,3 (83) From second equation, we get On using (82) in the above equation we get expression as ZB ZB ᾱp,1 + ᾱp,3 = ρ̄∆u Similarly, from third equation     1 2 ZB ā ā ZB ZB ∆(ρE) = ᾱp,1 ū − p + ū ᾱp,2 + ᾱp,3 ū + p 2 γ(γ − 1) γ(γ − 1) (84) (85) or ∆   p 1 ā ZB ZB ZB ZB (ᾱp,3 − ᾱp,1 ) + ρu2 = ū(ᾱp,1 + ᾱp,3 )+ p (γ − 1) 2 γ(γ − 1) 1 ZB + ū2 ᾱp,2 2 (86) on using (84) in the above equation, we get 1 1 1 ā ZB ZB (ᾱp,3 − ᾱp,1 ) + ū2 ∆ρ ∆p + ∆(ρu2 ) = ρ̄ū∆u + p (γ − 1) 2 2 γ(γ − 1) 14 (87) Now, 12 ∆(ρu2 ) = ρ̄ū∆u + 12 ū2 ∆ρ is an equation for above defined averages values of ρ̄ and ū. Thus we are left with r γ ∆p ZB ZB (88) ᾱp,3 − ᾱp,1 = γ − 1 ā On solving (84) and (88) simultaneously, we get r ρ̄∆u γ ∆p ZB ᾱp,1 = − 2 γ − 1 2ā r γ ∆p ρ̄∆u ZB + ᾱp,3 = 2 γ − 1 2ā and (89) Let us consider convective subsystem of Zha and Bilgen type splitting. On comparing (13) and (15) and writing in finite difference analogue, we have ZB ∆F ZB = Āc c ∆U (90) ZB ZB ZB where Āc is now a function of left and right states, i.e., Āc = Āc (U L , U R ). It is worth nothing that (90) is just a relation and may not become an equation for already defined average values. Further ∆U is a column vector and from the theory of Jordan Forms we are able to get complete set of LI generalized eigenvectors. Thus we can form generalized basis for ∆U , i.e., ∆U = 3 X ZB ZB ᾱc,i R̄c,i (91) i=1 On using above expression in (90), we get ZB ∆F ZB = Āc c 3 X ZB ZB ᾱc,i R̄c,i (92) i=1 or ZB ZB ∆F ZB = ᾱc,1 Āc c ZB ZB ZB Āc R̄c,1 + ᾱc,2 ZB ZB ZB R̄c,2 + ᾱc,3 Āc ZB R̄c,3 (93) which is further equal to ZB ZB ZB
ZB ZB ZB ∆F ZB = ᾱc,1 λ̄c,1 R̄c,1 + ᾱc,2 λ̄ZB c c,2 R̄c,2 + R̄c,1 ZB ZB ZB + ᾱc,3 λ̄c,3 R̄c,3 (94) Now, ∆F +ZB must have contribution of positive part of eigenvalues only, i.e., p ZB ZB ZB ZB +ZB ZB +ZB ZB +ZB ∆F +ZB = ᾱc,1 λ̄c,1 R̄c,1 + ᾱc,2 λ̄c,2 R̄c,2 + ᾱc,3 λ̄c,3 R̄c,3 c ZB ZB + ᾱc,2 R̄c,1 (95) Similarly, ZB ZB ZB ZB −ZB ZB −ZB ZB −ZB ∆F −ZB = ᾱc,1 λ̄c,1 R̄c,1 + ᾱc,2 λ̄c,2 R̄c,2 + ᾱc,3 λ̄c,3 R̄c,3 c ZB ZB + ᾱc,2 R̄c,1 (96) Again, we define the standard Courant splitting for the eigenvalues as λ̄± c,i = λ̄c,i ± |λ̄c,i | 2 15 (97) On using standard upwinding along with above definition, we finally get ∆F +ZB − ∆F −ZB = c c 3 X ZB ZB λ̄ZB ᾱc,i c,i R̄c,i (98) i=1 As we can see the resultant of extra contribution, which is coming because of weak hyperbolicity of convective subsystem, turns out be equal to zero. Unlike pressure subsystem here we don’t need to find wave strengths because all eigenvalues corresponding to convective subsystem are equal, which leads to ∆F +ZB c or − ∆F −ZB c = λ̄ZB c 3 X ZB ZB ᾱc,i R̄c,i (99) i=1 − ∆F −ZB = λ̄ZB ∆F +ZB ∆U c c c (100) Now,  ∆U1   ∆(ρ)      ∆U = ∆U2  =  ∆(ρu)  ∆(ρE) ∆U3 (101) where, ∆U1 = ρR − ρL ∆U2 = ∆(ρu) = ρ̄∆u + ū∆ρ and   p 1 ∆U3 = ∆(ρE) = ∆ + ρu2 (γ − 1) 2 1 1 2 = ∆p + ū ∆ρ + 2ρ̄ū∆u) γ−1 2 (102) 1 Here, we did an experiment to check ∆(ρE) = γ−1 ∆p + 12 ū2 ∆ρ + 2ρ̄ū∆u) holds or not. As we know from the theory of gas dynamics [34], ratio of densities, i.e., ( ρρrl ) attains a constant value of 6 as Mach number M → ∞. Next, we define error3 = ∆(ρE) − 1 2 1 ∆p − ū ∆ρ + 2ρ̄ū∆u) γ−1 2 (103) and we consider a test case with variable Mach number taken from [33], which is given below.   2 − (γ−1) pl 2γM γ+1        1    pl pr γ+1 pr γM 2   γ−1 pl + 1          ρl  =  1.0  and ρr  =  pr γ+1   + γ−1 pl   ul ur 1.0 q  (2+(γ−1)M 2 )pr γ (2γM 2 +(1−γ))ρr As per expectations density ratio approaches to limit 6 as shown in Figure 1(a) and error3 , which is given in Figure 1(b) is not exactly zero and it fluctuates between limits −10−15 to 10−15 , which is anyhow very small. The possible reason of this could be the generation of round-off error and if we take macroscopic scale, error looks almost equal to zero as shown in Figure (2). 16 −15 x 10 6.5 2 ZBS−FDS scheme 6 1.5 5.5 5 1 error3 4 R rho /rho L 4.5 3.5 0.5 3 0 2.5 2 −0.5 1.5 ZBS−FDS scheme 1 −1 0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 Mach number axis 400 500 600 700 800 900 1000 Mach number axis (a) (b) Figure 1: (a) represents ratio of densities results and (b) represents error3 results at microscopic level 1 ZBS−FDS scheme 0.8 0.6 0.4 error3 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0 100 200 300 400 500 600 700 800 900 1000 Mach number axis Figure 2: represents error3 results at macroscopic level The final expressions in the finite volume framework, with ZBS-FDS scheme for equations are as follows. i ∆t h n n n = U − U n+1 F − F 1 1 j j j+ 2 j− 2 ∆x where the cell-interface fluxes, F I = F j± 21 , are defined by i   1 h 1 +ZB −ZB −ZB + ∆F − ∆F ∆F +ZB − ∆F F I = [F L + F R ] − p p c c 2 2 where   ρR − ρ L    −ZB ZB  ρ̄∆u + ū∆ρ ∆F +ZB − ∆F = λ̄   c c c   1 1 2 ū ∆ρ + 2ρū∆u) ∆p + γ−1 2 17 Euler (104) (105) (106) and ∆F +ZB − ∆F −ZB = p p 3 X ZB ZB ᾱp,i λ̄ZB p,i R̄p,i (107) i=1 4.2. Formulation of TVS-FDS scheme Let us consider pressure subsystem corresponding to Toro and Vázquez type splitting and on comparing (22) and (24), and on using the finite difference analogue we have, TV V ∆F T = Āp ∆U p (108) Similar procedure like what we did for pressure subsystem of Zha and Bilgen type splitting is followed here and finally we get V V ∆F +T − ∆F −T = p p 3 X TV TV V ᾱp,i λ̄T p,i R̄p,i (109) i=1 To determine right side of (109) completely, we need to find average values of eigenvalues along with wave strengths. Consider the linearization equation of pressure subsystem TV V for TVS-FDS scheme, ∆F T = Āp ∆U , i.e., p      0 0 0 0 ∆(ρ)     1  2 −(γ − 1)ū (γ − 1) ∆p (110)    ∆(ρu)   =  2 (γ − 1)ū γ ūā2 ā2 ( γ−1 )∆pu ∆(ρE) − (γ−1) + 21 γ ū3 (γ−1) − γ ū2 γ ū From the second equation, we get ∆p = 1 (γ − 1)ū2 ∆ρ − (γ − 1)ū∆(ρu) + (γ − 1)∆(ρE) 2 (111) or 1 (γ − 1)ū2 ∆ρ − (γ − 1)ū∆(ρu) 2 1 p ) + (γ − 1)∆(ρu2 ) + (γ − 1)∆( γ−1 2 (112) ū2 ∆(ρ) − 2ū∆(ρu) + ∆(ρu2 ) = 0 (113) ∆p = which further reduces to ⇒ ū = √ √ ρL u L + ρR u R √ √ ρL + ρR (114) Average density can be found by substituting ū in the relation ρR uR − ρL uL = ρ̄(uR − uL ) + ū(ρR − ρL ) (115) √ and after some simple calculation, we get ρ̄ = ρL ρR . Similarly, third equation can be written as     γ ūā2 1 ∆(pu) = − + γ ū3 ∆ρ γ−1 (γ − 1) 2 (116)   2 ā − γ ū2 ∆(ρu) + γ ū∆(ρE) + (γ − 1) 18 On expanding we have,       γ ūā2 1 3 ā2 2 ∆(pu) = − + γ ū ∆ρ + − γ ū ∆(ρu) γ−1 (γ − 1) 2 (γ − 1)   p 1 + ρu2 + γ ū∆ γ−1 2 (117) After further simplifications, above relation reduces to On using p = γ∆(pu) = ā2 ρ̄∆u + γ ū∆p (118) ∆(a2 ρu) = ū∆(a2 ρ) + ā2 ρ̄∆u (119) a2 ρ , we have γ as explained earlier, ā2 comes out as ā2 = √ √ ρL a2L + ρR a2R √ √ ρL + ρ R (120) Similarly, wave strengths for pressure subsystem of TVS-FDS scheme can be calculated from the relation 3 X TV TV ᾱp,i R̄p,i (121) ∆U = i=1 i.e.,  ∆(ρ)   0 1   1   0 1       TV  TV  TV   ∆(ρu)  = ᾱp,1   + ᾱp,2  ū  + ᾱp,3   1 2 1 ū−β̄ 1 ū+β̄ ∆(ρE) ū + 2 ( γ−1 ) ū + 2 ( γ−1 ) 2 ū (122) From the first equation, we get TV ᾱp,2 = ∆ρ (123) TV TV TV ∆(ρu) = ᾱp,1 + ūᾱp,1 + ᾱp,3 (124) TV TV ⇒ ᾱp,1 + ᾱp,3 = ρ̄∆u (125) Similarly, the second equation gives Third equation implies 1 1 u + β̄ 1 u − β̄ TV TV TV )}ᾱp,1 + u2 ᾱp,2 + {u + ( )}ᾱp,3 ∆(ρE) = {u + ( 2 γ−1 2 2 γ−1 (126) On rearrangement of terms and after some algebra, the above equation reduces to TV TV ᾱp,3 − ᾱp,1 = 1 {2∆p − ūρ̄∆u} β̄ 19 (127) TV TV as and ᾱp,3 Finally on comparing (125) and (127), we get both ᾱp,1 1 1 ∆p ρ̄∆u + ρ̄ū∆u − 2 β̄ 2β̄ 1 1 ∆p = ρ̄∆u − ρ̄ū∆u + 2 2β̄ β̄ TV ᾱp,1 = (128) TV ᾱp,3 (129) Therefore, the wave strengths are finally given by 1 1 ∆p TV ρ̄∆u + , ᾱp,2 = ∆ρ ρ̄ū∆u − 2 β̄ 2β̄ 1 1 ∆p = ρ̄∆u − ρ̄ū∆u + 2 β̄ 2β̄ TV ᾱp,1 = TV and ᾱp,3 (130) We know that the convective subsystem is weakly hyperbolic but we can still form a basis of generalized eigenvectors, i.e., ∆U = 3 X TV TV ᾱc,i R̄c,i (131) i=1 and on comparing (21) and (23), and after writing in a finite difference form, we have the relation 3 X TV TV V TV ∆F T = Ā ᾱc,i R̄c,i (132) c c i=1 or TV V TV ∆F T = ᾱc,1 Āc c TV TV TV R̄c,1 + ᾱc,2 Āc TV TV TV R̄c,2 + ᾱc,3 Āc TV (133) R̄c,3 Since Ā is non-diagonalizable, this means TV Āc TV TV V 6= λ̄T c,i R̄c,i R̄c,i TV TV for some i’s. We know R̄c,3 is a generalized eigenvector and corresponding to R̄c,2 , a Jordan chain of order two will be formed, i.e., TV Āc TV R̄c,2 TV TV V = λ̄T c,2 R̄c,2 and Āc TV R̄c,3 TV TV V = λ̄T c,3 R̄c,3 + R̄c,2 (134) On using above relations in (133), we get TV TV TV TV V TV TV TV TV TV TV TV ∆F T = ᾱc,1 λ̄c,1 R̄c,1 + ᾱc,2 λ̄c,2 R̄c,2 + ᾱc,3 λ̄c,3 R̄c,3 + ᾱc,3 R̄c,2 c (135) V V After using standard Courant splitting for the eigenvalues, ∆F +T and ∆F −T are c c given by TV TV TV V T V +T V T V +T V T V +T V ∆F +T = ᾱc,1 λ̄c,1 R̄c,1 + ᾱc,2 λ̄c,2 R̄c,2 + ᾱc,3 λ̄c,3 R̄c,3 c TV TV R̄c,2 + ᾱc,3 (136) and TV TV TV V T V −T V T V −T V T V −T V ∆F −T = ᾱc,1 λ̄c,1 R̄c,1 + ᾱc,2 λ̄c,2 R̄c,2 + ᾱc,3 λ̄c,3 R̄c,3 c TV TV + ᾱc,3 R̄c,2 20 (137) Therefore, we have V V ∆F +T − ∆F −T = c c 3 X T V +T V ᾱc,i λ̄c,i R̄c,i − 3 X TV ᾱc,i ⇒ V ∆F +T c − V ∆F −T c TV TV TV + ᾱc,3 R̄c,2 i=1 (138) V λ̄−T c,i TV R̄c,i − i=1 = 3 X i=1 or V V ∆F +T − ∆F −T = c c TV ᾱc,3 TV R̄c,2
TV V V TV ᾱc,i λ̄+T − λ̄−T R̄c,i c,i c,i 3 X TV V TV λ̄T ᾱc,i c,i R̄c,i (139) (140) i=1 Like previous case, resultant of extra contribution comes out equal to zero. Above relation can be further written as, V V V ∆F +T − ∆F −T = λ̄T c c c V where λ̄T c TV TV TV TV ᾱc,2 R̄c,2 + ᾱc,3 R̄c,3 TV V = λ̄T c,2 = λ̄c,3 .
(141) This can be further written as V V V − ∆F −T = λ̄T ∆F +T c c c  TV TV ∆U − ᾱc,1 R̄c,1  (142) In order to determine (142) fully, we need to find wave strengths which can be calculated from the relation 3 X TV TV ∆U = ᾱc,i R̄c,i (143) i=1 or   ∆(ρ)   TV  ∆(ρu)  = ᾱc,1 ∆(ρE)       x1 1 0     TV  TV   0 + ᾱc,2  ū  + ᾱc,3  1 + ūx1  1 2 1 ū + 12 ū2 x1 2 ū (144) From the first equation, we get TV TV ∆(ρ) = ᾱc,2 + x1 ᾱc,3 (145) From second equation, we have TV TV ∆(ρu) = ūᾱc,2 + (1 + ūx1 )ᾱc,3 TV TV ⇒ ρ̄∆u + ū∆ρ = ū ᾱc,2 + x1 ᾱc,3
(146) TV + ᾱc,3 (147) On using (145) in the above equation and after cancellation of some terms, we get TV ᾱc,3 = ρ̄∆u (148) TV ᾱc,2 = ∆ρ − x1 ρ̄∆u (149) On using (148) in (145) we get 21 Similarly, from third equation we have TV ∆(ρE) = ᾱc,1 + 1 2 TV ū ᾱc,2 + 2 1
TV ū + x1 ū2 ᾱc,3 2  p 1 1  TV + ρu2 ) = ᾱc,1 + ū2 ∆ρ − x1 ρ̄∆u + γ−1 2 2 ⇒ ∆( After little algebra we get, TV ᾱc,1 = (150) 1
ū + x1 ū2 ρ̄∆u (151) 2 1 ∆p (γ − 1) (152) Therefore, the wave strengths are finally given by TV ᾱc,1 = 1 TV TV ∆p , ᾱc,2 = ∆ρ − x1 ρ̄∆u and ᾱc,3 = ρ̄∆u (γ − 1) TV TV On using above calculated wave strength ᾱc,1 along with eigenvector R̄c,1 (142), we get   ρR − ρ L    V −T V TV  ρ̄∆u + ū∆ρ λ̄ ∆F +T − ∆F =   c c c   1 2 ū ∆ρ + 2ρ̄ū∆u) 2 (153) = e3 in (154) Thus, the final expressions for TVS-FDS scheme is written as, = U nj − U n+1 j i ∆t h n F j+ 1 − F nj− 1 2 2 ∆x (155) where the cell-interface fluxes, F I = F j± 21 , are defined by FI = i   1 1 h V V V V (156) + ∆F +T − ∆F −T ∆F +T − ∆F −T [F L + F R ] − p p c c 2 2 5. Results and discussion We first consider here smooth solution problem with periodic boundary conditions to check experimental order of convergence for both constructed upwind schemes. After that, both ZBS-FDS and TVS-FDS schemes are tested on various one-dimensional Riemann problems of gas dynamics. For most of numerical examples, computational domain lies between 0 and 1, i.e., 0 ≤ x ≤ 1.0 except for Sod’s shock tube and shock-entropy test problems. For each problem computational domain is divided into 100 equally spaced cells, except for the shock-entropy interaction test case and for the blast problem test case in which computational domain is divided into 800 and 3000 equally spaced cells respectively. 5.1. Experimental order of convergence (EOC) Order of accuracy of both constructed upwind schemes can be determined using the EOC analysis, as given below. E = C(∆x)s (157) 22 where, E is an error between the exact solution and the numerical solution using some appropriate norm. In particular, we are taking three norms namely, L1 , L2 and L∞ . Here, C is a constant, ∆x is grid spacing and s is the order of accuracy, which need to be calculated. On taking logarithms on both sides of (157), we get log E = log C + s log ∆x (158) which is the equation of a straight line with slope s. For a given norm, let us initially take ∆x = h1 and on using this in (158), we get log Enorm,h1 = log C + s log h1 Next we take ∆x = h2 , preferably h2 = h1 2 , (159) with same norm in (159), we get log Enorm,h2 = log C + s log h2 (160) and on subtracting (160) from (159), formula for finding experimental order of convergence “s” comes out as follows. s = (log Enorm,h1 − log Enorm,h2 ) (log h1 − log h2 ) (161) To check the performance of both schemes in term of errors associated with each norm for different grid sizes, we choose a test case from [1] with initial smooth conditions, for one dimensional Euler system, which are given below. ρ(x, t) = 1.0 + 0.2 sin(π(x − ut)), u(x, t) = 0.1, p(x, t) = 0.5. For the present case, final solutions remain smooth and periodic boundary conditions are being employed to get meaningful solutions. Computational domain is chosen as [0, 2], i.e., 0 ≤ x ≤ 2 and all solutions are obtained at final time t = 0.5. L1 error norm for both schemes are given in Table 1 and as per expectations, there is reduction in error on refinement of mess size. Similarly, L2 error norm and L∞ error norm are given in Table 2 and Table 3. It is clear from all three tables that performance of both schemes is similar, if no discontinuity is present in the solution. Table 1: L1 error norm for smooth solution problem corresponding to both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.004476 0.002529 0.001258 0.000624 0.000308 EOC 0.82 1.007 1.011 1.018 TVS-FDS scheme 0.004476 0.002529 0.001258 0.000624 0.000308 EOC 0.82 1.007 1.011 1.018 5.2. Sod’s shock tube problem and Lax problem First we consider Sod’s shock tube problem in which, an initial discontinuity in the middle evolves to, going from right to left as we observe, a shock, a contact discontinuity and an expansion fan. The initial conditions [11] are (ρL , uL , pL ) = (1.0, 0.0, 100000.0), (ρR , uR , pR ) = (0.125, 0.0, 10000.0) with initial discontinuity at xo = 0.0 and computational domain lies between −10 to 10. All numerical results are obtained 23 Table 2: L2 error norm for smooth solution problem corresponding to both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.005238 0.003227 0.001707 0.000885 0.000452 EOC 0.6988 0.9187 0.9477 0.9693 TVS-FDS scheme 0.005238 0.003227 0.001707 0.000885 0.000452 EOC 0.6988 0.9187 0.9477 0.9693 Table 3: L∞ error norm for smooth solution problem corresponding to both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.019067 0.013968 0.007847 0.003994 0.002006 EOC 0.44 0.83 0.974 0.9935 TVS-FDS scheme 0.019067 0.013968 0.007847 0.003994 0.002006 EOC 0.44 0.83 0.974 0.9935 at final time t = 0.01. For this test case, both ZBS-FDS scheme and TVS-FDS scheme exhibit almost similar results except near normal shock region, where ZBS-FDS scheme performs slightly better than TVS-FDS scheme. Results corresponding to density variable are presented in Figure 3(a). We also present error analysis of both schemes corresponding to L1 -norm and L2 -norm, and results are given in Table 4, Table 5 respectively. Error analysis indicates that ZBS-FDS scheme is a little more accurate than TVS-FDS scheme. Next we consider Lax test case for which initial conditions are given as (ρL , uL , pL ) = (0.445, 0.698, 3.528), (ρR , uR , pR ) = (0.5, 0.0, 0.571) with xo = 0.5 and all numerical results are obtained at final time t = 0.15. The contact and shock discontinuities here are stronger than those in Sod’s shock tube problem. Results of both ZBS-FDS and TVS-FDS schemes are given in Figure 3(b). For this problem too, ZBS-FDS scheme performs slightly better than TVS-FDS scheme. 1.4 1 ZBS−FDS scheme TVS−FDS scheme Exact plot ZBS−FDS scheme TVS−FDS scheme 0.9 1.2 Exact plot 0.8 1 Density Density 0.7 0.6 0.5 0.8 0.6 0.4 0.3 0.4 0.2 0.1 −10 −8 −6 −4 −2 0 2 4 6 8 0.2 10 0 0.1 0.2 0.3 0.4 x 0.5 0.6 0.7 0.8 0.9 x (a) (b) Figure 3: (a) represents results of density variable for Sod’s shock tube problem and (b) represents density plots for Lax problem. 24 1 Table 4: L1 error norm for the Sod’s shock tube problem for both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.502947 0.352076 0.235865 0.156230 0.101461 TVS-FDS scheme 0.582406 0.397561 0.268590 0.176909 0.114140 Table 5: L2 error norm for the Sod’s shock tube problem for both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.177260 0.134255 0.097736 0.073471 0.055553 TVS-FDS scheme 0.196831 0.144438 0.105535 0.078268 0.058590 5.3. Sonic point problem and strong shock problem Next, we present numerical results of both schemes for a modified version of Sod’s problem. For this problem, solution has a right shock wave, a right travelling contact discontinuity and a left sonic rarefaction wave. This test case is useful in assessing the entropy condition satisfaction property of numerical methods. Initial conditions for this problem are given as (ρL , uL , pL ) = (1.0, 0.75, 1.0), (ρR , uR , pR ) = (0.125, 0.0, 0.1) with initial discontinuity at xo = 0.3 and all numerical solutions are obtained at final time t = 0.2. For this test case, low diffusive schemes like Roe scheme may violate the entropy condition and give unphysical rarefaction shocks in the expansion region at sonic points. To avoid this drawback, additional numerical diffusion is typically required, which is usually introduced as an entropy fix and one such famous fix is given by Harten [4]. Because of sufficient inbuilt numerical diffusion, both ZBS-FDS scheme and TVS-FDS scheme are seen to satisfy the entropy condition, as can be seen in the results shown in Figure 4(a) with no rarefaction shock or non-smoothness being present in the solution. For this problem, error analysis of L1 -norm and L2 -norm show TVS-FDS scheme is a slightly more accurate as given in Table 6 and Table 7. Second test case is taken from [24] and is designed to assess the robustness and accuracy of numerical methods. Its solution consists of a strong shock wave with Mach number 198, a contact discontinuity and a left rarefaction wave. Initial conditions are given as (ρL , uL , pL ) = (1.0, 0.0, 1000.0), (ρR , uR , pR ) = (1.0, 0.0, 0.01) with xo = 0.5 and all solutions are obtained at time t = 0.012. Both schemes work well and results are given in Figure 4(b). Error analysis of L1 -norm and L2 -norm show ZBS-FDS scheme is a little more accurate and results are given in Table 8 and Table 9. 5.4. Stationary contact discontinuity A contact discontinuity occurs when a family of characteristics are parallel to each other in the x − t domain. Since fluid velocity is the same on both sides, contact discontinuities move with fluid. The initial conditions as given in [24] are (ρL , uL , pL ) = (1.4, 0.0, 1.0) and (ρR , uR , pR ) = (1.0, 0.0, 1.0). The initial discontinuity is present at xo = 0.5. Both ZBS-FDS and TVS-FDS schemes capture the steady contact discontinuity exactly, as shown in Figure 5(a). 25 1 6 ZBS−FDS scheme TVS−FDS scheme Exact plot 0.9 ZBS−FDS scheme TVS−FDS scheme Exact plot 5 0.8 4 Density Density 0.7 0.6 0.5 0.4 3 2 0.3 1 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 1 0 0.1 0.2 0.3 0.4 x 0.5 0.6 0.7 0.8 0.9 x (a) (b) Figure 4: (a) represents density plots for sonic point problem and (b) represents density plots for strong shock problem. Table 6: L1 error norm for the sonic point problem for both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.038718 0.028065 0.019058 0.012698 0.008396 TVS-FDS scheme 0.036894 0.026387 0.017863 0.011879 0.007822 Table 7: L2 error norm for sonic point problem for both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.053061 0.043452 0.033069 0.025245 0.019600 TVS-FDS scheme 0.051063 0.040727 0.031264 0.024187 0.018961 Table 8: L1 error norm for the strong shock problem for both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.317106 0.241142 0.180898 0.131432 0.088449 TVS-FDS scheme 0.334709 0.258266 0.192025 0.138044 0.092496 5.5. Strong shock problem with slowly moving contact discontinuity This test case is also devised to test the robustness of numerical methods but the main reason for devising this test case is to assess the ability of numerical methods to resolve slowly-moving contact discontinuities. The exact solution of this test consists 26 1 Table 9: L2 error norm for strong shock problem for both schemes grid points 40 80 160 320 640 ZBS-FDS scheme 0.824979 0.665983 0.558651 0.473423 0.366668 TVS-FDS scheme 0.856110 0.699312 0.574795 0.479052 0.372136 of a left rarefaction wave, a right-travelling shock wave and a slowly moving contact discontinuity. Initial conditions are given as (ρL , uL , pL ) = (1.0, −19.59745, 1000.0), (ρR , uR , pR ) = (1.0, −19.59745, 0.01) with xo = 0.8 and all numerical solutions are obtained at time t = 0.012. In case of ZBS-FDS scheme, numerical solution goes towards the top portion of slowly moving contact wave whereas, for TVS-FDS scheme it is a little below as given in Figure 5(b). 6 ZBS−FDS scheme 1.4 ZBS−FDS scheme TVS−FDS scheme Exact plot TVS−FDS scheme Exact plot 1.35 5 1.3 4 Density Density 1.25 1.2 3 1.15 2 1.1 1.05 1 1 0.95 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 1 0 0.1 0.2 0.3 0.4 x 0.5 0.6 0.7 0.8 0.9 x (a) (b) Figure 5: (a) represents density plot for stationary contact discontinuity problem and (b) represents density plots for strong shock problem with slowly moving contact discontinuity. 5.6. Slowly moving shock Sometimes numerical methods tend to produce oscillations near the shock regions, which are completely unphysical. The oscillations associated with slowly moving shock problems are usually linked with lack of sufficient numerical diffusion in the scheme. We took a test case from [23] with initial conditions as (ρL , mL , EL ) = (3.86, −3.1266, 27.0913) and (ρR , mR , ER ) = (1.0, −3.44, 8.4168), where m = ρu is momentum and E is total energy. Final solutions are obtained at t = 4 units and results are given in Figure 6(a). 5.7. Mach 3 problem The initial conditions for this problem are, (ρL , uL , pL ) = (3.857, 0.92, 10.333) and (ρR , uR , pR ) = (1.0, 3.55, 1.0) with xo = 0.4 and all solutions are obtained at t = 0.1 units. This problem consists of a supersonic flow with Mach number 3 in expansion region and it produces a strong expansion fan. Low diffusive upwind schemes such as 27 1 4 4 ZBS−FDS scheme TVS−FDS scheme Exact plot 3.5 3.5 3 Density Density 3 2.5 2.5 2 2 1.5 1.5 ZBS−FDS scheme TVS−FDS scheme Exact plot 1 0 0.1 0.2 0.3 0.4 1 0.5 0.6 0.7 0.8 0.9 0.5 1 0 0.1 0.2 0.3 0.4 x 0.5 0.6 0.7 0.8 0.9 x (a) (b) Figure 6: (a) represents density plots for slowly moving shock problem and (b) represents density plots for Mach 3 problem. Roe’s approximate solver fail for this problem and require an entropy fix. According to Wesseling [29], even after use of Harten’s entropy fix, Roe scheme still gives sonic glitch. Both ZBS-FDS and TVS-FDS schemes perform well without needing any entropy fix and results are given in Figure 6(b). 5.8. Interacting blast wave problem This is one of the most severe test cases used to assess the numerical algorithm for its performance and is taken from Woodward and Colella [30]. Computational domain is divided into 3000 equally spaced finite volumes. Initial conditions for density and velocity are constants and given by ρ = 1.0, u = 0. For pressure variable, two discontinuities are present at position x = 0.1 and 0.9. Initially, pL = 1000 if x ∈ [0.0, 0.1] , pM = 0.01 if x ∈ [0.1, 0.9] and pR = 100 if x ∈ [0.9, 1.0]. Solution of this problem consists of multiple shocks, contact discontinuities and expansions waves. Results for ZBS-FDS scheme are given at two different time levels, as shown in Figure 7(a) and 7(b). For this test case, TVS-FDS scheme blew up in our simulations. 5.9. Shock-entropy wave interaction The shock-entropy wave interaction test case considered here for testing the present schemes is taken from [3], with computational domain x ∈ [−1, 1] being divided into 800 equally spaced finite volumes and all solutions are obtained at final time t = 0.47. The initial conditions are given below.   (ρL , uL , pL ) = 3.857143, 2.629369, 10.3333 if x < −0.8   (162) (ρR , uR , pR ) = 1 + 0.2 sin(5πx), 0, 1 if x > −0.8. In this problem, a Mach 3 shock wave interacts with density disturbances created by perturbing the initial density. This initial disturbance gives rise to the continuous interaction of smooth flow with the discontinuities. Similar kind of interaction can be observed in compressible turbulence. Therefore, this is a suitable test case to test the scheme for its ability to resolve complex interactions, which can be used in turbulent computations. First order results for ZBS-FDS scheme are presented in Figure 8(a). To 28 1 7 6 ZBS−FDS scheme Reference ZBS−FDS scheme Reference 6 5 5 Density Density 4 3 4 3 2 2 1 0 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 1 0 0.1 0.2 0.3 0.4 x 0.5 0.6 0.7 0.8 0.9 x (a) (b) Figure 7: (a) Density plot, for blast wave problem at time t = 0.026 units and (b) density plot, for same problem at time t = 0.038 units. achieve second order accuracy, we used Venkatakrishnan’s limiter which is a modified version of van Albada limiter [28] and deals with piecewise linear reconstruction of primitive variables. As an example, let us consider a piecewise linear reconstruction for density variable, i.e., to obtain

1 ∆2+ + ǫ2 ∆− + ∆2− + ǫ2 ∆+ ρL = ρ + (163) i i+1/2 2 ∆2+ + ∆2+ + 2ǫ2 where, ∆+ = ρi+1 − ρi (164) ǫ2 = (K∆x)3 (165) ∆− = ρi − ρi−1 and Similarly, other primitive variables can be reconstructed. Here, K is a constant and ∆x is a grid spacing. Large values of K indicate no limiting and in the present case, we take K = 0.1. For this problem, second order results for ZBS-FDS scheme are computed and are given in Figure 8(b). Comparison of both first order results and second order results for ZBS-FDS scheme are given in Figure 9. Results for TVS-FDS scheme are presented in Figure 10(a), 10(b). Both the schemes produce results of nearly similar accuracy for this test case. In both cases, second order accurate results are substantially better, compared to the first order accurate results. 6. Two-dimensional Euler system The 2-D Euler equations form a system of four coupled non-linear hyperbolic PDEs with independent space variables x, y and independent time variable t. In the differential, as well as conservative, form the 2-D Euler system can be written as ∂F 1 ∂F 2 ∂U + + = 0 ∂t ∂x ∂y 29 (166) 1 5.5 5.5 ZBS−FDS scheme Reference 4 3.5 3.5 Density 4.5 4 3 2.5 3 2.5 2 2 1.5 1.5 1 −1 2nd−order−ZBS−FDS scheme Reference 5 4.5 1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 x 0 0.2 0.4 0.6 0.8 x (a) (b) Figure 8: (a) 1st-order results for ZBS-FDS scheme, for shock-entropy wave interaction problem and (b) 2nd-order results for ZBS-FDS scheme for same problem. 5 4.5 4 Density axis Density 5 3.5 3 2.5 2 1st−order ZBS−FDS scheme 2nd−order ZBS−FDS scheme Reference 1.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Position axis Figure 9: Comparison of 1st-order and 2nd-order numerical results of ZBS-FDS scheme for shock-entropy wave interaction problem. where, U is vector of conserved variables and F 1 , F 2 are flux vectors which are given as follows.       ρv ρu ρ        ρuv   ρu2 + p   ρu        (167) U =   , F1 =   and F 2 =  ρv 2 + p     ρuv   ρv  ρvE + pv ρuE + pu ρE Using the divergence form, equation (166) can be written as ∂U + ∇.F = 0 ∂t 30 (168) 1 5.5 5.5 TVS−FDS scheme Reference 4.5 4.5 4 4 3.5 3.5 3 3 2.5 2.5 2 2 1.5 1.5 1 −1 2nd−order−TVS−FDS scheme Reference 5 Density Density 5 1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 x 0 0.2 0.4 0.6 0.8 x (a) (b) Figure 10: (a) Density plot for 1st-order TVS-FDS scheme, for shock-entropy wave interaction problem and (b) density plot for 2nd-order TVS-FDS scheme for same problem. where  ρu⊥     ρuu⊥ + pnx   F =  ρvu + pn  ⊥ y   ρEu⊥ + pu⊥ (169) is the flux vector and the vector u⊥ is defined as the scalar product of the velocity vector u and the unit normal vector n, i.e., u⊥ = u.n = nx u + ny v (170) where nx and ny represent the direction cosines of the unit normal n̂ to the cell-interface and are given by ∆y ∆x nx = , ny = − (171) ∆s ∆s On integrating (168) over domain Ω with boundary ∂Ω and on further using Green’s theorem, we get Z I ∂ U dΩ + F ds = 0 (172) ∂t Ω ∂Ω On calculating average value of U over Ω, the first integral can be re-written and the above equation becomes I ∂ Ū 1 F ds (173) = − ∂t A ∂Ω where A is area enclosed by Ω. For typical two dimensional finite volume, for general quadrilaterals, integral as given by second term of above semi-integral from can be approximated by line integrals. After a little algebra, the discretized finite volume form for 2-D Euler system is given by ∂ Ū m 1 X k k = − F ∆s ∂t Am k 31 (174) 1 k=2 m k= 3 u‫װ‬ y ∆y ∆S k u⊥ k= 1 x k=4 ∆x Figure 11: Schematic representation of general quadrilateral. where subscript m denotes the cell number, k is the cell-interface index of the mth cell. Similarly, F k and ∆sk are the normal flux and the perimeter of the k th face. Further, Am is the area of mth cell. In 1-D, both ZBS-FDS and TVS-FDS schemes performed reasonably well for most of the important test cases except for the blast wave problem (5.8), where TVS-FDS blew up quickly. Otherwise, it is bit difficult to judge which of two is more accurate as ZBS-FDS scheme produces slightly better results for Sod’s shock tube problem, Lax problem and for both strong shock problems, whereas TVS-FDS scores over ZBS-FDS in case of sonic point problem, slowly moving shock problem and for Mach 3 problem. In this work, we opt for Zha and Bilgen type splitting in formulating FDS concept based numerical scheme for 2-D Euler system, while emphasizing that this choice is purely based on our convenience and in future other possibility can be explored. 6.1. Analysis of Zha and Bilgen type splitting in 2-D The flux vector F in the 2-D case is split into a convection part and a pressure part, based on Zha and Bilgen splitting, as follows. F = F ZB + F ZB c p where F ZB c  ρu⊥  (175)  0       ρuu⊥   pnx  ZB     =  and F p =  pn   ρvu⊥   y ρEu⊥ pu⊥ Let AZB denote convection flux Jacobian matrix which is given below. c   0 nx ny 0    −uu⊥ u⊥ + unx uny 0 ZB   Ac =  vnx u⊥ + vny 0   −vu⊥ −Eu⊥ Enx Eny u⊥ 32 (176) (177) Eigenvalues corresponding to matrix AZB are real and equal with set of eigenvalues as c u⊥ , u⊥ , u⊥ , u⊥ . Analysis of AZB shows that it has a defective set of LI eigenvectors, c i.e.,       nx ny 0            u 0 0  ⊥  , RZB   and RZB   RZB (178) c,1 =  c,2 =  c,3 =      0  u⊥   0 0 0 1 On evaluating rank of matrices (AZB − u⊥ I4 ), (AZB − u⊥ I4 )2 ..., we find that there c c ZB will be one Jordan block of order two as rank(Ac − u⊥ I4 )2 = rank(AZB − u ⊥ I4 ) 3 . c ZB Let R(A) denote the space spanned by the columns of matrix Ac − u⊥ I4 . Then, as explained in 1-D case, we have R(A) = x1 A1 + x2 A2 + x3 A3 + x4 A4 where, A1 , A2 , A3 and A4 are column vectors of AZB − u⊥ I4 . Now c        −u⊥ nx ny          −uu⊥   unx   uny        R(A) = x1  + x2  + x3  + x4       −vu⊥   vnx   vny   −Eu⊥ Enx Eny or  1 (179) 0   0  0  0 (180)     u  R(A) = (−u⊥ x1 + nx x2 + ny x3 )   v   E (181) The column vector (1, u, v, E)t , which is a range space of R(A), becomes generalized eigenvector of AZB as null space N (AX) is just a scalar coefficient. Let us take AZB X1 c c which is equal to    0 nx ny 0 1     −uu⊥ u⊥ + unx uny 0   u    (182)  −vu   vnx u⊥ + vny 0 ⊥   v −Eu⊥ Enx Eny u⊥ E which, on solving, is equal to    u⊥    1 u   v (183) E Thus, AZB X 1 = u⊥ X 1 holds. Now, this generalized eigenvector is expected to form c a Jordan chain of order two corresponding to matrix AZB , i.e., c AZB X 1 = u⊥ X 1 c X 2 = u⊥ X 2 + X 1 AZB c 33 (184) Other generalized eigenvector X 2 can be found from second relation from, it is written as      x1 x1 0 nx ny 0       −uu⊥ u⊥ + unx x  uny 0   x2     = u⊥  2  +  x      −vu vnx u⊥ + vny 0   x3  ⊥  3  x4 −Eu⊥ Enx Eny u⊥ x4 and, in expanded  1     u    v   E (185) here, each xi ∈ IR, where i runs from 1 to 4, and on solving all four simultaneous equations, we get nx x 2 + ny x 3 = 1 + u ⊥ x 1 (186) which together with (x1 , x2 , x3 , x4 )t defines a generalized eigenvector. Now, if we take   1 x1 nx ny    u x2 u⊥ 0  P =  (187)  v x3 0 u⊥    E x4 0 0 then P −1 AZB P c  u⊥   0 =   0  0 1 0 u⊥ 0 0 u⊥ 0 0 0   0  holds. 0  u⊥ (188) Let AZB denote the Jacobian matrix corresponding to pressure flux function F ZB p p . After a little algebra AZB comes out equal to p    (γ − 1)    0 0 0 Θ 2 nx −nx u −nx v Φ 2 nx − u ⊥ u Φ 2 ny − u ⊥ v Θ 2 ny
Θ 2 − Φ 2 u⊥ −ny u −ny v 0   nx   ny   u⊥ (189) where we define u2 + v 2 and 2 2 a = γ(γ − 1) Θ2 = Φ2 (190) The eigenvalues of the flux Jacobian matrix AZB are: p λZB p,1 = − r γ−1 ZB ZB a, λZB p,2 = 0, λp,3 = 0, λp,4 = γ 34 r γ−1 a γ (191) Since all eigenvalues are real and distinct, therefore AZB must have full set of LI p eigenvectors and are given by: 0  RZB p,1 RZB p,3    =     nx ny u⊥ − p    =    1 a γ(γ − 1)  uk       uuk + Θ2 ny     =  , RZB p,2  vu − Θ2 n   x k   0   0   nx  nx u ⊥    , RZB =  p,4 ny  ny u ⊥    a u⊥ + p u2⊥ − Θ2 γ(γ − 1) (192)        (193) Now, both convection and pressure fluxes at a cell-interface are calculated by using the following.  1 1  ZB F ZB = F c,L + F ZB |ū⊥ | ∆U (194) c,I c,R − 2 2 where   ρR − ρL     ρ̄∆u + ū∆ρ     (195) ∆U =   ρ̄∆v + v̄∆ρ    
1 1 2 2 ∆p + ū + v̄ ∆ρ + ρ̄(ū∆v + v̄∆v) γ−1 2 and 4 F ZB p,I =  1 X ZB ZB ZB 1  ZB λ̄p,i R̄p,i ᾱ F p,L + F ZB p,R − 2 2 i=1 p,i respectively. Like in 1-D case, the average quantities are defined as √ √ ρL u L + ρR u R √ ρ̄ = ρL ρR , ū = √ √ ρ L + ρR √ √ √ √ ρ L v L + ρR v R ρL a2L + ρR a2R v̄ = , ā2 = √ √ √ √ ρL + ρR ρ L + ρR (196) (197) ū⊥ = ūnx + v̄ny and ūk = −ūny + v̄nx where uk denotes velocity component parallel to cell-interface and is given by uk = −uny + vnx 35 (198) ZB , where i runs from 1 to 4, are given as Similarly, wave strengths ᾱp,i r γ ∆p ρ̄∆u⊥ − 2 γ − 1 2ā ūk ∆ρ + ρ̄∆uk = Θ̄2 − ū2⊥ ZB ᾱp,1 = ZB ᾱp,2 ZB ᾱp,3 = ∆ρ − ZB ᾱp,4 = ρ̄∆u⊥ 2 ū2k ∆ρ + ρ̄ūk ∆uk (199) Θ̄2 − ū2⊥ r γ ∆p + γ − 1 2ā where ū2 + v̄ 2 2 ∆ρ = ρR − ρL Θ̄2 = ∆u⊥ = u⊥ R − u⊥ ∆uk = uk R − uk ∆p = pR − pL L (200) L 6.2. Numerical examples In this subsection, the ZBS-FDS scheme is tested on various well-established benchmark test problems. Special attention is given to problems with complex interactions of strong shocks which leads to various shock instabilities. Many well known upwind schemes are known to produce unphysical features in such cases [18]. 6.2.1. Oblique shock reflection In this test case [31], an oblique shock wave is introduced at the top left corner by means of initial conditions and post shock boundary conditions, at the left and top side of the domain, respectively. The computational domain considered for this test case is [0, 3] × [0, 1]. The initial conditions for this test problem are as given below. (ρ, u, v, p)0,y,t = (1.0, 2.9, 0, 1/1.4) (ρ, u, v, p)x,1,t = (1.69997, 2.61934, −0.50633, 1.52819) The incident shock angle measured from the top side of the domain is 290 and the free stream Mach number M = 2.9. Wall boundary conditions are prescribed at the bottom boundary and supersonic outflow boundary conditions are used at the right side of the computational domain. Both first order and second order results on various grids are presented in Figure 12. 6.2.2. Supersonic flow across a compression ramp in a wind tunnel The computational domain of [0, 3] × [0, 1] is considered for this test problem. Other geometrical features of the problem include a 150 ramp at the lower part of the computational domain. In this two-dimensional steady test case [12], supersonic flow of a Mach number M = 2 encounters a fifteen degree ramp to form an oblique shock wave. This shock wave reflects from the upper wall and interacts with the expansion wave generated at the tip of the ramp corner. The so weakened expansion wave again reflects from the 36 (a) (a) 1 1 0.5 0.5 0 0 0.5 1 1.5 2 2.5 0 3 0 0.5 1 (b) 1 1 0.5 0.5 0 0 0.5 1 1.5 1.5 2 2.5 3 2 2.5 3 (b) 2 2.5 0 3 0 0.5 1 1.5 Figure 12: First order results of ZBS-FDS scheme are presented on left, where second order results are given in right side for shock reflection problem; pressure contours (0.7: 0.1: 2.9) on the grids: (a) 120 × 40 and (b) 240 × 80 top wall and further interacts with the second reflected shock from the ramp surface. Initial conditions are prescribed at the left boundary, wall boundary conditions are used at the top and on the ramp, and supersonic outflow boundary conditions are imposed at the exit boundary. Both first order and second order results are presented in Figure 13. (a) (a) 1 1 0.5 0.5 0 0 0.5 1 1.5 2 2.5 0 3 0 0.5 1 (b) 1 1 0.5 0.5 0 0 0.5 1 1.5 1.5 2 2.5 3 2 2.5 3 (b) 2 2.5 3 0 0 0.5 1 1.5 Figure 13: First order results of ZBS-FDS scheme are given on left and second order results are presented on right for ramp reflection problem; pressure contours (1.1: 0.05: 3.8) on the grids: (a) 120 × 40 and (b) 240 × 80 6.2.3. Reflection of a plane shock from wedge In this is a two-dimensional problem in which the reflection of a plane shock wave from a wedge lies in the double-Mach reflection regime, some Riemann solvers are known to generate kinked Mach stems [18]. In this test case, the kinked Mach stem occurs when a strong normal shock wave moving with Mach 5.5 encounters the 300 ramp to form Mach reflection and represents a typical shock instability phenomenon. Three shocks meet to form a triple point and the computational domain considered for this problem 37 is [0, 2.0] × [0, 1.5] with initial shock location at x0 = 0.25. All computational results are obtained at time t = 0.25. The computational domain to the right of the shock is initialized with a rest fluid of density 1.4 and pressure 1. To the left of the shock, values obtained from the moving shock relations for Mach 5.5 are used to initialize the domain. Figure 14 shows the density contours computed with the present scheme and no kink is observed. 1.5 1.5 1 1 0.5 0.5 0 0 0.5 1 1.5 2 0 0 0.5 1 1.5 2 Figure 14: First (left) and Second (right) order results of reflection of a plane shock from a wedge problem with ZBS-FDS scheme on 400 × 400 grid points 6.2.4. Hypersonic flow past a half-cylinder The hypersonic flow around a half-cylinder is also a well-known test case to examine the capability of numerical methods in resolving complex flow features accurately without giving shock instabilities. In this case, the shock instability is known as ’carbuncle shock’ which was initially reported by Peery and Imlay [17]. This test is computed for Mach 6 and Mach 20 flows on fine and coarse grids in circumferential directions and results are given in Figure 15. Many Riemann solvers generate carbuncle shocks [18, 10] in the their numerical solutions. As an example we presented results of Roe scheme, where carbuncle phenomena can be seen very clearly. The present method did not exhibit any such phenomena. 7. Summary In this study, we attempted to develop flux difference split upwind schemes for convection-pressure splitting frameworks, based on Jordan canonical forms to avoid defective matrices. FDS solver for Liou and Steffen type splitting is not attractive as pressure subsystem doesn’t have any contribution of acoustic signals. Newly constructed ZBS-FDS and TVS-FDS schemes are tested on various benchmark problems and don’t require any entropy fix for sonic point and strong expansion problems. ZBS-FDS scheme and TVS-FDS scheme perform in similar ways as is evident from several 1-D test cases. We further extend ZBS-FDS scheme to 2-D Euler system, specifically to those test problems for which several Riemann solvers generate shock instabilities [18, 15, 6, 21]. The performance of ZBS-FDS scheme for these test cases is impressive and deserves further research. For example, accuracy of ZBS-FDS scheme in multi-dimensions can be further improved using some good diffusion regulator such as [7] and other possible direction is to pursue genuinely multi-dimensional modeling, as the eigenvalues of the 38 (a) Roe (b) (c) (d) 1.5 1.5 1.5 1.5 1.5 1 1 1 1 1 0.5 0.5 0.5 0.5 0.5 0 0 0 0 0 -0.5 -0.5 -0.5 -0.5 -0.5 -1 -1 -1 -1 -1 -1.5 0.5 1 1.5 -1.5 0.5 -1.5 1 1.5 0.5 1 1.5 -1.5 0.5 1 1.5 -1.5 0.5 1 1.5 Figure 15: First order results of ZBS-FDS scheme for half cylinder problem; density contours (2.0: 0.2: 5.0): (a) Mach 6 on 45 × 45 grid, (b) Mach 6 on 20 × 320 grid, (c) Mach 20 on 45 × 45 and (d) Mach 20 on 20 × 320 grid convection and pressure parts of the fluxes neatly reflect uni-directional and multidirectional information propagation respectively. 8. Acknowledgments The authors thank Prof. Michael Junk, Fachbereich Mathematik und Statistik, Universität Konstanz, Germany for very useful discussions. The authors also thank Indian Institute of Science for supporting this research. References [1] K. R. Arun, M. Lukáčová-Medvidová, P. Prasad and S. V. Raghurama Rao (2013). A second order accurate kinetic relaxation scheme for inviscid compressible flows. In Recent Developments in the Numerics of Nonlinear Hyperbolic Conservation Laws (pp. 1-24). Springer Berlin Heidelberg. [2] N. Balakrishnan and S. M. Deshpande (1995). New upwind method exploiting the wave-particle behavior of fluid flow. CFD Journal, 3(4), 433-446. [3] D. S. Balsara and C. W. Shu (2000). Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. Journal of Computational Physics, 160(2), 405-452. [4] A. Harten (1984). On a class of high resolution total-variation-stable finite-difference schemes. SIAM Journal on Numerical Analysis, 21(1), 1-23. [5] A. Harten, P. D. Lax and B. V. Leer (1983). On upstream differencing and Godunovtype schemes for hyperbolic conservation laws. SIAM review, 25(1), 35-61. [6] K. Huang, H. Wu, H. Yu and D. Yan (2011). Cures for numerical shock instability in HLLC solver. International journal for numerical methods in fluids, 65(9), 1026-1038. 39 [7] S. Jaisankar and S.V. Raghurama Rao (2007). Diffusion regulation for Euler solvers. Journal of computational Physics, 221(2), 577–599. [8] A. Jameson (1995). Analysis and design of numerical schemes for gas dynamics, 1: Artificial diffusion, upwind biasing, limiters and their effect on accuracy and multi-grid convergence, International Journal of Computational Fluid Dynamics, 4 (3-4), 171-218. [9] A. Jameson (1995). Analysis and design of numerical schemes for gas dynamics, 2: Artificial diffusion and discrete shock structure. International Journal of Computational Fluid Dynamics, 5(1-2), 1-38. [10] S. S. Kim, C. Kim, O. H. Rho and S. K. Hong (2003). Cures for the shock instability: development of a shock-stable Roe scheme. Journal of Computational Physics, 185(2), 342-374. [11] C. B. Laney (1998). Computational gasdynamics. Cambridge University Press. [12] D. W. Levy, K. G. Powell and B. van Leer (1993). Use of a rotated Riemann solver for the two-dimensional Euler equations. Journal of Computational Physics, 106(2), 201-214. [13] M. S. Liou, Ten Years in the Making - AUSM Family, AIAA Paper no. AIAA-20012521; also NASA/TM 2001-210977, 2001. [14] M. S. Liou and C. J. Steffen (1993). A new flux splitting scheme. Journal of Computational physics, 107(1), 23-39. [15] J. C. Mandal and V. Panwar (2012). Robust HLL-type Riemann solver capable of resolving contact discontinuity. Computers & Fluids, 63, 148-164. [16] S. Osher and F. Solomon (1982). Upwind difference schemes for hyperbolic systems of conservation laws. Mathematics of computation, 38(158), 339-374. [17] K. M. Peery and S. T. Imlay (1988). Blunt-body flow simulations. AIAA paper, 88, 2904. [18] J. J. Quirk (1997). A contribution to the great Riemann solver debate (pp. 550-569). Springer Berlin Heidelberg. [19] S.V. Raghurama Rao (1995). Peculiar velocity based upwind method for inviscid compressible flows. Computational Fluid Dynamics Journal, Vol. 3, No. 4, pp. 415-432. Japanese Society for CFD, Japan. [20] P. L. Roe (1981). Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of computational physics, 43(2), 357-372. [21] Z. Shen, W. Yan and G. Yuan (2016). A robust HLLC-type Riemann solver for strong shock. Journal of Computational Physics, 309, 185-206. [22] J. L. Steger and R. F. Warming (1981). Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods. Journal of computational physics, 40(2), 263-293. 40 [23] Y. Stiriba and R. Donat (2003). A numerical study of post-shock oscillations in slowly moving shock waves. Computers & Mathematics with Applications, 46(5), 719-739. [24] E. F. Toro (2009). Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media. [25] E. F. Toro, M. Spruce and W. Speares (1994). Restoration of the contact surface in the HLL-Riemann solver. Shock waves, 4(1), 25-34. [26] E.F. Toro and M.E. Vázquez-Cendón (2012). Flux splitting schemes for the Euler equations. Computers & Fluids, 70, pp.1-12. [27] B. Van Leer (1979). Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. Journal of computational Physics, 32(1), 101-136. [28] V. Venkatakrishnan (1995). Convergence to steady state solutions of the Euler equations on unstructured grids with limiters. Journal of computational physics, 118(1), 120-130. [29] P. Wesseling (2009). Principles of computational fluid dynamics (Vol. 29). Springer Science & Business Media. [30] P. Woodward and P. Colella (1984). The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of computational physics, 54(1), 115-173. [31] H. C. Yee, R. F. Warming, and A. Harten (1982). A high-resolution numerical technique for inviscid gas-dynamic problems with weak solutions. In Eighth International Conference on Numerical Methods in Fluid Dynamics (pp. 546-552). Springer Berlin Heidelberg. [32] G. C. Zha and E. Bilgen (1993). Numerical solutions of Euler equations by using a new flux vector splitting scheme. International Journal for Numerical Methods in Fluids, 17(2), 115-144. [33] S. Zhang and C. W. Shu (2007). A new smoothness indicator for the WENO schemes and its effect on the convergence to steady state solutions. Journal of Scientific Computing, 31(1-2), 273-305. [34] M. J. Zucrow and J. D. Hoffman (1976). Gas Dynamics, Vol. I. John Wiley and Sons, New York, 112-115. 41