An accurate and robust genuinely multidimensional Riemann solver for Euler equations based on TV flux splitting
AAn accurate and robust genuinely multidimensional Riemannsolver for Euler equations based on TV flux splitting
S. Sangeeth and J. C. Mandal
Department of Aerospace Engineering,Indian Institute of Technology Bombay,Mumbai, 400076, India
Abstract
A simple robust genuinenly multidimensional convective pressure split (CPS) , contact preserving, shockstable Riemann solver (GM-K-CUSP-X) for Euler equations of gasdynamics is developed. The convective andpressure components of the Euler system are seperated following the Toro-Vazquez type PDE flux splitting [Toroet al, 2012]. Upwind discretization of these components are achieved using the framework of Mandal et al [Man-dal et al, 2015]. The robustness of the scheme is studied on a few two dimensional test problems. The resultsdemonstrate the e ffi cacy of the scheme over the corresponding conventional two state version of the solver. Re-sults from two classic strong shock test cases associated with the infamous Carbuncle phenomenon, indicatethat the present solver is completely free of any such numerical instabilities albeit possessing contact resolutionabilities.Such a finding emphasizes the preexisting notion about the positive e ff ects that multidimensional flowmodeling may have towards curing of shock instabilities. Keywords : Genuinely multidimensional, Riemann solver, Contact preserving, Convective Pressure Split, numer-ical shock instability, Carbuncle phenomenon.
Past few decades have witnessed commendable advancement in the computation of high speed flows. Birth ofthe upwind schemes for finite volume methods is a milestone in this regard. Among these characteristics basedschemes, those based on the solution of two state ‘Riemann problems’at cell interface are the most popular class ofschemes lately. Godunov [4] proposed the first Riemann problem based algorithm (also called Riemann solvers)that used the exact solution of a Riemann problem to construct a first order accurate scheme. However, this exact1 a r X i v : . [ m a t h . NA ] M a r olver was critized for being prohibitively expensive because of the iterative technique involved [5]. To mitigate thisshortcoming, a class of approximate Riemann solvers have been proposed and continue to be developed. Detailsabout some famous approximate Riemann solvers can be found in [5, 6, 7, 8]. Such algorithms have brought forthan era of much cheaper yet accurate computation of gasdynamic flows.Since most of these approximate schemes are inherently one dimensional in their framework, they enjoy accu-rate and robust resolution of both linear and nonlinear wave fields arising in one dimensional problems. However,accurate resolution of flow features in practically relevant multidimensional problems where important flow fea-tures are oblique to the grid still pose a major challenge. This is mainly because the wave system in such problemspossess infinitely many propagation direction as compared to the limited ones in a one-dimensional problem. Fur-ther vorticity wave enters the formulation and has to be dealt with appropriately [9]. Currently, the standard wayof extending the one-dimensional schemes to multidimensions is through a dimensionally split operator approachmotivated in [10] where local one dimensional Riemann problems are solved in a direction normal to each cellinterface. Such a strategy has been criticised because of its unnatural selection of the interface normal as the onlywave propagation direction even for problems with infinitely many propagation directions. This lacuna is attributedto cause for example, pressure disturbances across a grid oblique shear wave [11].Another disturbing problem encountered by dimensionally split approximate Riemann solvers is the occurrenceof various forms of numerical shock instabilities when simulating strong normal shocks. A catalouge of many suchfailures can be found in [12]. It has been observed that only those schemes that have exact shear wave resolutionabilities are known to produce these instabilities. There is increasing evidence [13, 14] to consider that suchfailures occur due to lack of adequate dissipation across cell faces which are normal to the shock wave front.These problems solicit a strong need for a genuinely multidimensional formulation that resolves all characteristicfields while introducing appropriate dissipation along necessary directions.Probably the earliest attempt at creating a multidimensional Riemann solver was by Raithby [15] who sug-gested discretizing the convective terms along flow dominant directions. Davis [16] and others [17, 18, 19] intro-duced the idea of rotated Riemann solvers that involves identifying grid oblique shock orientations and solvingRiemann problems across them.Roe and others [9, 20] suggested using extrapolated data from the left and right states across a Riemann solverto construct simple wave solutions for a multidimensional Riemann problem. This technique is not generally validfor arbitrary governing equations and demands substantial reformulation for three dimensional problems.Collela [21] is credited with proposing the CTU method that uses the characteristics of the system to incorporatediagonal cell contirbutions to interface fluxes in a predictor corrector framework.Leveque [22] proposed a multidimensional scheme in which cross derivative terms that constitute the transversecontributions are included by interpreting the usual one dimensional Gudonov method in a fluctuation splitting2ramework.Ren et al [23] constructed an operator split predictor corrector scheme based on CTU and Leveque’s wavepropagation method for both Euler and Navier Stokes systems. While the predictor step solves linearized Eulerequations in characteristic variables, the corrector step adds viscous contributions.Wendro ff [24] introduced a multistate Riemann solver that attempted the extension of one dimentional HLLEscheme [6] into several dimensions. The scheme used an expensive nine point stencil and su ff ered unacceptabledissipation due to unnatural wavespeed selection [26].Another multistate Riemann solver is from the work of Brio et al [25]. Multidimensional e ff ects are incorpo-rated by adding correction terms to the standard face normal fluxes at every computational cell interface. Thesecorrection terms are obtained by solving a three state Riemann problem using Roe’s FDS at the corners of therespective cells.In a similar spirit, Balsara [26] presented a generic multidimensional HLLE solver (GM-HLLE) with simpleclosed form expressions that can be extended to any hyperbolic system. This method too relies on construction ofmultidimensional correction terms like [25] wherein these terms are obtained using a four state HLLE Riemannsolver at interface corners.By building upon the wave model introduced in GM-HLLE, Balsara [27] further proposed a multidimensionalHLLC solver for Euler and MHD systems. To deal with contact discontinuities in two dimensions, a set of twelvepossible contact orientations on a given cell are included in the wave model. Although posessing closed formexpressions in two dimensions, such a wave model would not be easily tractable in three dimensions.Improvements to the model was proposed in [28] by reformulating the scheme in terms of characteristic vari-ables but this too remains complicated to implement on a computer code.Mandal et al [2] proposed a multidimensional convective-pressure split scheme (GM-HLLCPS-Z) based onZha-Bilgen [29] way of splitting Euler fluxes. The scheme consists of a wave speed averaged upwinding for theconvective part and GM-HLLE type discretization for the pressure part. This scheme is basically a multidimen-sional extension of HLL-CPS strategy [30] that uses a HLL type discretization for distinctly treating convectiveand pressure flux parts. The splitting of full Euler flux into convective and pressure parts can be achieved by adopt-ing AUSM type or Zha-Bilgen type PDE splitting [8, 29]. Di ff usion control is achieved by careful tuning of thedissipation vector of the pressure flux. This renders the scheme with stationary and moving contact preserving abil-ities in addition to capability of surviving the most stringent test problems. Toro has corroborated such a methodby showing that exact contact preserving ability can be incorporated into convective-pressure split framework bydiscretizing the pressure fluxes using a Riemann solver [1]. Surprisingly, although contact preserving, HLL-CPSis found to evade most common forms of numerical shock instabilities, particularly the carbuncle phenomenon.Recently Xie et al [31] has proposed a contact capturing convective-pressure split scheme named K-CUSP-X.3he scheme uses exactly similar discretization as HLL-CPS for both convective and pressure fluxes but di ff ers onlyin the method of splitting the Euler fluxes into these components; instead of using Zha-Bilgen or AUSM way ofPDE splitting, it adopts Toro-Vazquez method wherein the pressure terms embedded in the energy is also seperated[1]. However unlike HLL-CPS, the present investigations reveal that this scheme is found to su ff er from numericalshock instabilites.In this paper, a new genuinely multidimensional scheme based on the conventional K-CUSP-X is proposed. Amultidimensional correction term similar to [25] and [26] is constructed by solving appropriate four state Riemannproblem at the corners of each interface. The multidimensional terms are incorporated such that the final fluxescan be calculated at interfaces with the familiar ease of pre-existing dimensional split methods.Adhereing to the existing convention, this new scheme may be called GM-K-CUSP-X. It will be the purposeof this paper to demonstrate that such a construction is not only as accurate as GM-HLL-CPS-Z, but also curesthe numerical shock instability that plagued the corresponding two state conventional K-CUSP-X Riemann solver.The positive e ff ect of genuinely multidimensional flow modeling on K-CUSP-X was first demonstrated in [33]wherein authors demonstrated that shock instability associated with standing shock problem [32] was completelyremoved by extending the original two state solver into its genuinely multidimensional variant. Such a findingcorroborates the pre existing opinion that multidimensional dissipation acts as a reliable cure for numerical shockinstability problems and incentivizes the need for extending the accurate Riemann solvers in literature into theirrobust genuinely multidimensional counterparts.This paper is organized as follows. In the next section, the governing equations and the type of PDE fluxsplitting used is detailed. In Section 3, the second order version of the new formulation is described. In Section4, results for some complex flow problems like double Mach reflection and multidimensional Riemann problemare discussed. Further, two classic shock instability test problems, the odd-even decoupling problem and standingshock problem are used to study instability behaviour of the newly developed scheme. Section 5 contains someconcluding remarks. Consider the two-dimensional Euler equations in di ff erential form ∂ U ∂ t + ∂ F ∂ x + ∂ G ∂ y = U = [ ρ ρ u ρ v ρ e ] T is the vector of conserved quantities. F and G are the flux vectors in the x and y directionsrespectively given by F = ρ u ρ u + p ρ uvu ( ρ e + p ) G = ρ v ρ uv ρ v + pv ( ρ e + p ) (2)In the present work, the above flux vectors are split into corresponding convective and pressure parts following theapproach of Toro-Vazquez [1] . Using ideal gas law, the split flux vectors can be written as F = u ρρ u ρ v ρ ( u + v ) F = p γγ − pu G = v ρρ u ρ v ρ ( u + v ) G = p γγ − pv (3)where F and G are the convective fluxes while F and G are the pressure fluxes. γ represents the ratio ofspecific heat capacities. As suggested by Toro et al [1] the convective fluxes F and G can be interpreted assimple advection of mass, momentum and kinetic energy along the x and y directions respectively, while thepressure fluxes F and G are interpreted to be sonic impulses that spread in all directions with reference to theseconvecting particles. This type of PDE splitting primarily di ff ers from those that already exist in the literature like[29] and [8], in terms of the quantity concerning energy that is being advected. Such di ff erences may have strongbearing on the robustness of these schemes. Figure 1: Cell ( i , j ) of area ∆ x × ∆ y and corners c , c , c and c in a Cartesian grid.5onsider the integral form of equation (1) ∂∂ t (cid:90) Ω U d Ω + (cid:73) d Ω ( F + G ) · ˆ n ds = n = ( n x , n y ) is the unit vector along the face normal. Consider a Cartesian cell of area ∆ x × ∆ y as shown inFigure 1. The finite volume discretization for the cell ( i , j ) can be written as¯ U n + i , j = ¯ U ni , j + ∆ t ∆ x ( ¯ F (cid:48) i − , j − ¯ F (cid:48) i + , j ) + ∆ t ∆ y ( ¯ G (cid:48) i , j − − ¯ G (cid:48) i , j + ) (5)where n is the time level, i and j are cell indices. The ¯( . ) quantities depict the respective averaged quantities. Thetotal fluxes at the cell interfaces ( i + , j ) and ( i − , j ) are denoted as ¯ F (cid:48) i + , j and ¯ F (cid:48) i − , j respectively. Similarly ¯ G (cid:48) i , j + and ¯ G (cid:48) i , j − are the y directional total fluxes in corresponding y interfaces. In the proposed scheme the interface flux will comprise of both a two state conventional Riemann flux and amultidimensional flux. These multidimensional fluxes are sought from the solution of four state Riemann problemthat occurs at each corner of a Cartesian cell. To see this more clearly, consider a typical cell ( i , j ) as shown inFigure 2. The four constant states that come together at corner c at t = t = T , for sake of simplicity, the simple wave model consistingof only four waves as introduced in [26] is used in this work. Accordingly the wave propagation at the corners isassumed to span, at any time T, a rectangular region. Figure 3 represents a three dimensional view of these wavesas they evolve in time where the shaded region depicts the domain of influence of this multidimensional Riemannproblem.In principle, the Equation 1 can be integrated along the space-time volume of Figure 3 appropriately to obtainthe time averaged multidimensional fluxes F ∗ in x-direction and G ∗ in y-direction. For a typical cell interface( i + , j ) the total flux will be an ensemble of the two multidimensional fluxes from corners c and c denoted as F ∗ i + , j + , F ∗ i + , j − and the single mid-point flux F midi + , j . This is shown in Figure 4. A conservative ensemble of thecorner and mid point fluxes are achieved by using Simpson’s rule of intergration [26] along the interface as givenby, 6 F (cid:48) i + , j = F ∗ i + , j + + F midi + , j + F ∗ i + , j − (6)Figure 2: Figure showing a typical cell ( i , j ) in a Cartesian computational grid.Figure 3: Figure showing the evolution of four waves making up the simple wave model used in this work.Figure 4: Figure showing the various components of an interface flux in the present genuinely multidimensionalscheme. This section deals with determination of the two state Riemann flux at the mid point of a typical interface ( i + , j )denoted as F midi + , j . Following the Convective Pressure Split (CPS) philosophy the total flux at this interface is7rst split into convective and pressure parts. This work uses the Toro-Vazquez type flux splitting as mentioned inEquation (3), F midi + = F midi + , j + F midi + , j (7)These convective and pressure parts are discretized independently following the original HLL-CPS strategy[30]. It may be noted that the interface ( i + , j ) admits only the x-directional Riemann flux F midi + , j while they-directional flux, G midi + , j is zero on it. midi + , j ) The upwind discretization of the convective flux denoted by F midi + , j at the mid point of the cell interface ( i + , j ),following the strategy of HLL-CPS method [30], is given by, F midi + , j = M k ρρ u ρ v ρ ( u + v ) k a k (8) k = L if ¯ u ≥ R if ¯ u < M k = ¯ u ¯ u − S cL if ¯ u ≥ u ¯ u − S cR if ¯ u < a k = u L − S cL if ¯ u ≥ u R − S cR if ¯ u < u = u L + u R S cL and S cR are8arefully selected wave speeds which are discussed in section 3.1.3. midi + , j ) The pressure flux at the midpoint of the cell interface is obtained by applying a HLL type discretization of thepressure flux vector [30]. F midi + , j = S cR S cR − S cL F L − S cL S cR − S cL F R + S cR S cL S cR − S cL ( U R − U L ) (12)The above equation can be rewritten as F midi + , j =
12 ( F L + F R ) + δ U (13)where F L and F R are the left and right side x-directional pressure fluxes normal to the interface ( i + , j ) and δ U is the numerical di ff usion given by δ U = S cR + S cL S cR − S cL ) ( F L − F R ) − S cL S cR S cR − S cL ρ L − ρ R ( ρ u ) L − ( ρ u ) R ( ρ v ) L − ( ρ v ) R ( ρ e ) L − ( ρ e ) R (14)It is clear from Equation 14 that the second term will give rise to numerical di ff usion across a contact. Thus inorder to capture the contact discontinuity accurately, density terms in δ U are replaced by the pressure terms byusing isentropic assumption ¯ a c = δ p δρ as described in the reference [30]. δ U = S cR + S cL S cR − S cL ) ( F L − F R ) − S cR S cL ¯ a c ( S cR − S cL ) p L − p R ( pu ) L − ( pu ) R ( pv ) L − ( pv ) R ¯ a c ( p L − p R ) + [( pq ) L − ( pq ) R ] (15)where ¯ a c is the average speed of sound at the interface given by ¯ a c = a L + a R q = u + v is twice the localkinetic energy per unit mass. 9 .1.3 Selection of wave speeds for the 1D Riemann problem at the interface (i + , j) The wave speeds are selected according to conventional HLL-CPS method [30] S L = min (0 , u L − a L , u ∗ − c ∗ ) S R = max (0 , u r + a r , u ∗ + c ∗ ) (16)where u ∗ and c ∗ are given by [30], u ∗ = u L + u R + a L − a R γ − c ∗ = a L + a R + γ −
14 ( u L − u R ) (17)Supersonic conditions are taken care by including ’0’ in the above expressions. For a stationary flow the wavespeeds are modified as S cL = − ¯ a c S cR = ¯ a c (18)It is to be noted that midpoint fluxes at other interfaces can be obtained in the similar manner. This section will detail how to evaluate the fluxes F ∗ i + , j + and G ∗ i + , j + that results from the interaction of fourRiemann states at a representative corner c . These fluxes forms the multidimensional component of the interfacefluxes and are shown in Figure 5. Once again, resorting to the (CPS) philosophy, the total flux at this corner is splitinto convective ( F ∗ i + , j + and G ∗ i + , j + ) and pressure fluxes ( F ∗ i + , j + and G ∗ i + , j + ) originating at the corner asper Equation 3. Analogous to the procedure for evaluation of the two state mid point flux at the interface, the splitconvective and the pressure parts at the corners will also be evaluated using di ff erent upwind strategies. c of the cell interface (F ∗ i + , j + ) Following [2], the x-directional convective flux is evaluated as,10igure 5: Figure showing showing the multidimensional fluxes F ∗ i + , j + and G ∗ i + , j + at corner c of ( i + , j )interface of a typical cell ( i , j ) [2]. F ∗ i + , j + = ¯ u S U ρρ u ρ v ρ ( u + v ) k − S D ρρ u ρ v ρ ( u + v ) k S U − S D (19)where ¯ u is the wave speed averaged x-directional local fluid speed at the corner defined as,¯ u = u LU S U − u LD S D + u RU S U − u RD S D S U − S D ) (20)Based on the direction of the average x- directional flow, the upwind states k , k are chosen as,If ¯ u > k = LU and k = LD (i.e upwinding is done from the left states).If ¯ u < k = RU and k = RD (i.e upwinding is done from the right states).In similar spirit, the y directional convective flux is evaluated as, G ∗ i + , j + = ¯ v S R ρρ u ρ v ρ ( u + v ) k − S L ρρ u ρ v ρ ( u + v ) k S R − S L (21)where ¯ v is the the wave speed averaged y-directional local fluid speed at the corner defined as,¯ v = v RU S R − v LU S L + v RD S R − v LD S L S R − S L ) (22)11he states k , k are chosen accordingly as,If ¯ v > k = RD and k = LD (i.e upwinding is done from the down states).If ¯ v < k = RU and k = LU (i.e upwinding is done from the up states).Since the above strategy is developed for a subsonic case, slight modification ¯ u and ¯ v is done to extend theabove formulation to supersonic cases:1. If the flow is supersonic in positive x-direction, then at the corner upwinding is done from left states. Thereforefor the evaluation of x-directional convective flux ¯ u is taken as¯ u = u LU S U − u LD S D S U − S D
2. If the flow is supersonic in negative x-direction the upwinding is done from right states. Therefore for theevaluation of x-directional convective flux ¯ u is taken as¯ u = u RU S U − u RD S D S U − S D
3. If the flow is supersonic in positive y-direction the upwinding is done from down states. Therefore for theevaluation of y-directional convective flux ¯ v is taken as¯ v = v RD S R − v LD S L S R − S L
4. If the flow is supersonic in negative y-direction the upwinding is done from upper states. Therefore for theevaluation of y-directional convective flux ¯ v is taken as¯ v = v RU S R − v LU S L S R − S L c of the cell interface (F ∗ i + , j + ) Following [2], the x-directional convective flux is evaluated as, F ∗ i + , j + = F LU S R S U + F RD S L S D − F LD S R S D − F RU S L S U ( S R − S L )( S U − S D ) − S R S L ( S R − S L )( S U − S D ) ( G RU − G LU + G LD − G RD ) + S R S L ( S R − S L )( S U − S D ) ( S U ( U RU − U LU ) − S D ( U RD − U LD )) (23)12imilarly the y-directional pressure flux is evaluated as, G ∗ i + , j + = G RD S R S U + G LU S L S D − G RU S R S D − G LD S L S U ( S R − S L )( S U − S D ) − S U S D ( S R − S L )( S U − S D ) ( F RU − F LU + F LD − F RD ) + S U S D ( S R − S L )( S U − S D ) ( S R ( U RU − U RD ) − S L ( U LU − U LD )) (24)The above equations can be rewritten as, F ∗ i + , j + =
12 ( F L + F R ) + δ U − S R S L ( S R − S L )( S U − S D ) ( G RU − G LU + G LD − G RD ) (25) G ∗ i + , j + =
12 ( G D + G U ) + δ U − S U S D ( S R − S L )( S U − S D ) ( F RU − F LU + F LD − F RD ) (26)where F L = F LU S U − F LD S D S U − S D (27) F R = F RU S U − F RD S D S U − S D (28) G D = G RD S R − G LD S L S R − S L (29) G U = G RU S R − G LU S L S R − S L (30)13he δ U and δ U terms in Equations (25,26), are the numerical di ff usion terms in x and y-directions respec-tively. To remove the numerical dissipation across a contact wave, the dissipation terms in x and y-directions areremodeled using isentropic expression as, δ U = S R + S L S R − S L ) ( F L − F R ) − S R S L ( S R − S L )( S U − S D )(¯ a ) S U ( p LU − p RU ) − S D ( p LD − p RD ) S U (( pu ) LU − ( pu ) RU ) − S D (( pu ) LD − ( pu ) RD ) S U (( pv ) LU − ( pv ) RU ) − S D (( pv ) LD − ( pv ) RD ) S U ( e ∗ LU − e ∗ RU ) − S D ( e ∗ LD − e ∗ RD ) (31) δ U = S U + S D S U − S D ) ( G D − G U ) − S U S D ( S R − S L )( S U − S D )(¯ a ) S R ( p RD − p RU ) − S L ( p LD − p LU ) S R (( pu ) RD − ( pu ) RU ) − S L (( pu ) LD − ( pu ) LU ) S R (( pv ) RD − ( pv ) RU ) − S L (( pv ) LD − ( pv ) LU ) S R ( e ∗ RD − e ∗ RU ) − S L ( e ∗ LD − e ∗ LU ) (32)where e ∗ k is given as e ∗ k = ¯ a γ − p k + p k ( u + v ) k and ¯ a as ¯ a = a LU + a RU + a LD + a RD c of theinterface can be obtained in a similar manner. Most importantly, it must be noted that while F ∗ i + , j ± contributes tothe total interface flux at the interface ( i + , j ), G ∗ i + , j + and G ∗ i + , j − contributes to the total interface flux at theinterfaces ( i , j + ) and ( i , j − ) respectively. As previously mentioned, the present work adopts a simple wave model as proposed by [26] to represent the wavesemerging from the four state Riemann problem at the corners of every interface. A top view of the area swept14y these four waves S L , S R , S U , S D is shown in Figure 6. The rectangle ABCD in Figure 6 depicts the domain ofinfluence of the four state Riemann problem at corner c at time T on x-y plane. An estimate for these wavespeedsFigure 6: Figure showing the rectangular area swept by the four waves that constitute the simple wave model usedin this work at any time t = T. [2]can be obtained as, S R = max (0 , λ Nx ( U RU ) , λ Nx ( U RD ) , ¯ λ Nx ( U LU , U RU ) , ˜ λ Nx ( U LD , U RD )) S L = min (0 , λ x ( U LU ) , λ x ( U LD ) , ¯ λ x ( U LU , U RU ) , ˜ λ x ( U LD , U RD )) S U = max (0 , λ Ny ( U RU ) , λ Ny ( U LU ) , ¯ λ Ny ( U RD , U RU ) , ˜ λ Ny ( U LD , U LU )) S D = min (0 , λ y ( U RD ) , λ y ( U LD ) , ¯ λ y ( U RD , U RU ) , ˜ λ y ( U LD , U LU )) (33)where typically, λ x ( U k ) denote smallest x-directional wave speed in the state U k , λ Nx ( U k ) denote largest x-directional wave speed in the state U k ,˜ λ x ( U k , U l ) denotes smallest x-directional wave speed from Roe averaged state between U k and U l ,˜ λ Nx ( U k , U l ) denotes largest x-directional wave speed from Roe averaged state between U k and U l such that k , l ∈{ LU , LD , RU , RD } and k (cid:44) l .If ¯ u = v (cid:44)
0, then x-directional wave speeds are modified as S L = − ¯ a and S R = ¯ a .If ¯ u (cid:44) v =
0, then y-directional wave speeds are modified as S D = − ¯ a and S U = ¯ a .If ¯ u = v =
0, then wave speeds are modified as S L = − ¯ a , S R = ¯ a , S D = − ¯ a and S U = ¯ a .It should be noted that zero has been added in the above expressions in order to ensure fully one sided flux insupersonic flow. 15 Results
A second order accurate version of the present solver is developed using the SDWLS strategy [34] whose details areomitted here for brevity. Formal order of accuracy of the second order version of GM-K-CUSP-X is investigatedusing the isentropic vortex problem [2]. The problem involves an isentropic vortex centered initially at (0,0) andmade to traverse the domain diagonally under a periodic boundary condition. A Cartesian domain of size [-5,5] × [-5,5] is used. The initial conditions consists of a unperturbed state given by ( ρ, p , u , v ) = (1 . , . , . , . T = p ρ . A perturbation is added to the flow as,( δ u , δ v ) = (cid:15) π e . − r ) ( − y , x ) δ T = − γ − γπ (cid:15) e (1 − r ) δρ = δ T γ − δ p = δρ γ Here, (cid:15) which defines the strength of the vortex is taken as 5.0. r denotes the Cartesian distance from vortex center.The accuracy of the scheme is measured in L and L ∞ norms of the density variable. The formal order of accuracy O is obtained by the formula, O = log ( η ) − log ( η ) log ( ∆ x ) − log ( ∆ x ) (34)where, η and η depict consecutive norms ( L or L ∞ ) for progressively refined grids with dimensions ∆ x and ∆ x respectively. The results are tabulated in table 1.Mesh size L error L order L ∞ error L ∞ order64 ×
64 0.007724 - 0.145543 -128 ×
128 0.001949 1.9866 0.044705 1.7029256 ×
256 0.000510 1.9341 0.013459 1.7318512 ×
512 0.000110 2.2129 0.002416 2.4778Table 1: Second order accuracy analysis of GM-K-CUSP-X using SDWLS reconstruction strategy.It is observed from the analysis that GM-K-CUSP-X scheme is able to achieve second order accuracy onsu ffi ciently refined grids in both L and L ∞ norms. 16 .2 Two dimensional Riemann problems Two dimensional Riemann problems provides an excellent test case to assess the qualitative improvement of gen-uinely multidimensional formulation of GM-K-CUSP-X over the corresponding conventional two state K-CUSP-Xscheme. Second order versions of both solvers are used for these test cases. The first two dimensional Riemannproblem investigated in the present work consists of a double Mach reflection and an oblique shock wave propa-gating at an angle to the grid. The initial conditions of the 2D Riemann problem are given by [2],Zone ρ p u vx > y > x > y < x < y > x < y < × − , × [ − ,
1] is used. The solution is evolved up to a timeof 1.05 units with a CFL number of 0.95 as suggested in reference [27]. The density contours at the final time areshown in Figure 7. Result for K-CUSP-X solver under identical conditions is given for comparison. It is observedthat the mushroom cap structure is well resolved by GM-K-CUSP-X as compared to the original K-CUSP-X.Further GM-K-CUSP-X is able to resolve the Kelvin-Helmholtz roll up much better than K-CUSP-X.17 a) GM-K-CUSP-X(b) Original K-CUSP-X
Figure 7: Results for the Riemann problem 1 showing 30 iso density contours from 0.15 to 1.7The second Riemann problem consist evolution of two weak shock waves and two contact waves and can beset using conditions [2], Zone ρ p u vx > y > x > y < x < y > x < y < × [-1,1] domain is used. CFL number of 0.95 is used and thesolution is evolved for a time of 0.5 units. The results obtained are represented using iso density contours. For18 a) GM-K-CUSP-X(b) Original K-CUSP-X Figure 8: Results for the Riemann problem 2 showing 30 iso density contours from 0.5 to 1.7comparison, results obtained under identical conditions by corresponding two state conventional K-CUSP-X is alsoprovided in Figure 8. Although both solvers are able to resolve the resultant contact waves and Mach reflections[26], only GM-K-CUSP-X is able to resolve the Kelvin Helmholtz instability along the Mach stems.
This problems deals with a Mach 10 shock inclined at 60 o with the x-axis and propagating downstream of arectangular duct of size [0,4] × [0,1] and interacting with a reflective bottom boundary wall. Detailed initial andboundary conditions can be found in [2]. The problem is solved on a Cartesian mesh of size 1920 ×
480 with aCFL of 0.7 and evolved up to time of 0.2 units. The second order accurate results obtained for GM-K-CUSP-X are represented using iso density contours and are compared with that obtained using K-CUSP-X under same19onditions. It is visible from the result that the present solver is able to resolve the complex shock structures crisplyand also the slipping contact line emerging from the triple point. (a) GM-K-CUSP-X(b) Original K-CUSP-X
Figure 9: Results for double Mach reflection problem showing 25 iso density contours from 1.77 to 22.44
One of the objectives of the present work is to demonstrate the e ff ect of genuinely multidimensional formulation onthe numerical shock instability characteristics of a Riemann solver. Two standard test cases namely, the odd-evendecoupling problem [12] and the standing shock problem [32] will be used to carry out the investigations. Sinceshock instability is most explicitly observed in the first order simulations, first order version of solvers will be usedfor these tests. Odd-even decoupling test consists of a M = =
140 units. It is clearly observed that the genuinelymultidimensional extension is able to preserve the shock profile without any instabilities. (a) K-CUSP-X(b) GM-K-CUSP-X
Figure 10: Results for odd-even decoupling problem showing 40 iso density contours from 1.4 to 7.37
To clearly di ff erentiate the behavior of K-CUSP-X and GM-K-CUSP-X schemes, a standing shock instability testcase is used. The details of the test are available in [32]. As seen in the figure 11a, K-CUSP-X fails miserably inthis test case. However, from figure 11b it is very evident that GM-K-CUSP-X scheme is free of this instability.21 a) K-CUSP-X(b) GM-K-CUSP-X Figure 11: Results for odd-even decoupling problem showing 50 iso density contours from 1.0 to 5.4The reason for the stabilizing nature of GM-K-CUSP-X can be discerned from observing Equations (19), (21),(25), (26). Discretization of the convective terms using Equations (19), (21) employs a wave weighted averagingthat strives to accurately model the underlying multidimensional wave evolution phenomenon. Specifically, it iseasy to note that Equations (20) and (20) defines wave averaged convection velocities in x and y directions thatdepends on all the adjoining states at the corners. Such a formulation admits information from transverse cellswhich would contribute to dissipation that would help supressing the instabilities. A similar observation can befound from Equations (25), (26) which describe the discretization of the pressure terms. For example, in Equation25 for multidimensional flux F ∗ i + , j + , the third term on right hand side consist of terms G RU , G LU , G LD and22 RD clearly indicating the influence of y-directional flux terms in the evaluation of x-directional flux. Equation 26depicts a vice versa situation for the flux G ∗ i + , j + . These multidimensional coupling terms along with the waveaveraged convective terms may be providing the additional cross dissipation that damps any unprecedented growthin instabilities thereby making the present scheme immune to shock instabilities. These equations reveal that thediscretization of convective and pressure fluxes in the corner of interfaces using GM-HLLE strategy may have apositive impact in curing such instabilities. The present work introduces a new genuinely multidimensional contact preserving Riemann solver GM-K-CUSP-X. This scheme is based upon Toro-Vazquez type PDE level flux splitting and the ensuing convective-pressurefluxes are discretized following [30]. While the convective fluxes are upwinded based on appropriate wave speedsthat emerge from the interacting states, the pressure fluxes are treated in a HLLE framework. Restoration ofstationary contact preservation ability is improved by explicitly reducing the numerical dissipation in the pressureflux discretization. The resulting solver is found to produce improved results as compared to the original K-CUSP-X solver on standard test problems. Particularly, interesting flow features like Kelvin-Helmholtz roll up andmushroom cap structure in two dimensional Riemann problem and complex shock interactions in double Machreflection problem are well resolved. Further, the present genuinely multidimensional solver is able to mitigatevarious forms of shock instabilities that plagued the corresponding conventional two state Riemann solver, K-CUSP-X. Such a finding reassures the pre existing notion that multidimensional dissipation is one of the mostpromising methods to cure shock instabilities. Due to the simplicity of the formulation, the present solver can beeasily extended to unstructured framework and three dimensional problems in principle.
References [1] E.F. Toro, V. Cendon Flux splitting schemes for the Euler equations, Comput Fluids, 70(2012) 1-12.[2] J.C. Mandal , V. Sharma A genuinely multidimensional convective pressure flux split Riemann solver forEuler equations, J. Comput. Phys., 297(2015) 669-688.[3] J. C Mandal, J. Subramanian, On the link between weighted least-squares and limiters used in higher-orderreconstructions for finite volume computations of hyperbolic equations, Appl Numer Math, 58(2008) 705-725.[4] S. K. Godunov, A finite di ff erence method for the computation of discontinuous solutions of the equations offluid dynamics, Mat. Sb., 47(1959) 357-393. 235] P. L. Roe, Approximate Riemann solvers, parameter vectors, and di ff erence schemes, J. Comput. Phys.,43(1981) 357-372.[6] A. Harten, P. D. Lax, B. van Leer, On upstream di ff erencing and Godunov-type schemes for hyperbolicconservation laws, SIAM Rev., 25(1983) 35-61.[7] E. F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann Solver, ShockWaves, 4(1994) 25-34[8] M. S. Liou, C. J. Ste ff en, A new flux splitting scheme, J. Comput. Phys., 107(1993) 23-39.[9] P. L. Roe, Discrete models for the numerical analysis of time dependent multidimensional gas dynamics, J.Comput. Phys., 63(1986) 458-476.[10] G. Strang, On the construction and comparison of di ff erence schemes, SIAM J. Numer. Anal., 5(1968) 506-517.[11] B. van Leer, Progress in multidimensional upwind di ff erencing, Proceedings of Thirteenth International con-ference on Numerical Methods in Fluid Dynamics, Lecture Notes in Physics, vol 414 (1993), Springer.[12] J. J. Quirk, A contribution to the great Riemann solver debate, Int. J. Numer. Methods Fluids, 18(1994)555-574.[13] Z. Shen, W. Yan, G. Yuan A robust HLLC-type Riemann solver for strong shock, JCP, 2016.[14] M. Pandolfi, D. D’Ambrosio, Numerical instabilities in upwind methods: Analysis and cures for the “Car-buncle” phenomenon, J. Comput. Phys., 166(2001) 271-301.[15] G. D. Raithby, Skew upstream di ff erencing schemes for problems involving fluid flow, Comput. Methods.Appl. Mech, 9(1976), 153-164.[16] S. F. Davis A rotationally biased upwind di ff erence scheme for the Euler equations, J. Comput. Phys.,56,65-92,1984[17] D. W. Levy, K. G. Powell, B. van Leer, Use of a rotated Riemann solver for the two dimensional Eulerequations, J. Comput. Phys.,106, 201-214,1993.[18] H. Nishikawa,K. Kitamura Very simple carbuncle free boundary layer resolving rotated hybrid Riemannsolvers, J. Comput. Phys.,227, 2560-2581,2008[19] Y. X. Ren, A robust shock-capturing scheme based on rotated Riemann solvers, Comput Fluids, 32(2003)1379-1403. 2420] C. L. Rumsey, B. v. Leer, P. L. Roe A grid independent approximate Riemann solver with applications to theEuler and Navier Stokes equations, 29th Aerospace sciences meeting, AIAA, Nevada, 1991.[21] P. Collela, Multidimensional upwind methods for hyperbolic conservation laws, JCP, 87(1990), 171.[22] R. J. Leveque, Wave propagation algorithms for multi-dimensional hyperbolic systems, JCP, 131(1997), 327-353.[23] Y. X. Ren, Y. Sun, A multi dimensional upwind scheme for solving Euler and Navier-Stokes equations, J.Comput. Phys., 219(2006) 391-403.[24] B. Wendro ff , A two-dimensional HLLE Riemann solver and associated Godunov-type di ffff