An accurate and robust genuinely multidimensional Riemann solver for   Euler equations based on TV flux splitting

S. Sangeeth; J. C. Mandal

arXiv:1703.09621·math.NA·March 29, 2017

An accurate and robust genuinely multidimensional Riemann solver for Euler equations based on TV flux splitting

S. Sangeeth, J. C. Mandal

PDF

Open Access

TL;DR

This paper introduces a new multidimensional Riemann solver for Euler equations that is robust, shock-stable, and preserves contact discontinuities, demonstrating superior performance on complex shock problems.

Contribution

A novel genuinely multidimensional Riemann solver based on TV flux splitting that effectively handles shock stability and contact preservation in Euler equations.

Findings

01

Solver is free of Carbuncle instability.

02

Demonstrates superior shock stability in 2D tests.

03

Effectively preserves contact discontinuities.

Abstract

A simple robust genuinely multidimensional convective pressure split (CPS) , contact preserving, shock stable Riemann solver (GM-K-CUSP-X) for Euler equations of gas dynamics is developed. The convective and pressure components of the Euler system are separated following the Toro-Vazquez type PDE flux splitting [Toro et al, 2012]. Upwind discretization of these components are achieved using the framework of Mandal et al [Mandal et al, 2015]. The robustness of the scheme is studied on a few two dimensional test problems. The results demonstrate the efficacy of the scheme over the corresponding conventional two state version of the solver. Results from two classic strong shock test cases associated with the infamous Carbuncle phenomenon, indicate that the present solver is completely free of any such numerical instabilities albeit possessing contact resolution abilities.Such a finding…

Figures29

Click any figure to enlarge with its caption.

Tables3

Table 1. Table 1: Second order accuracy analysis of GM-K-CUSP-X using SDWLS reconstruction strategy.

Mesh size	$L_{1}$ error	$L_{1}$ order	$L_{\infty}$ error	$L_{\infty}$ order
64 $\times$ 64	0.007724	-	0.145543	-
128 $\times$ 128	0.001949	1.9866	0.044705	1.7029
256 $\times$ 256	0.000510	1.9341	0.013459	1.7318
512 $\times$ 512	0.000110	2.2129	0.002416	2.4778

Table 2. Table 2: Initial condition for Riemann problem 1.

Zone	$ρ$	$p$	$u$	$v$
$x > 0$ , $y > 0$	1.5	1.5	0.0	0.0
$x > 0$ , $y < 0$	0.5323	0.3	0.0	1.206
$x < 0$ , $y > 0$	0.5323	0.3	1.206	0.0
$x < 0$ , $y < 0$	0.1379	0.029	1.206	1.206

Table 3. Table 3: Initial condition for Riemann problem 2.

Zone	$ρ$	$p$	$u$	$v$
$x > 0$ , $y > 0$	0.5313	0.4	0.0	0.0
$x > 0$ , $y < 0$	1.0	1.0	0.0	0.7276
$x < 0$ , $y > 0$	1.0	1.0	0.7276	0.0
$x < 0$ , $y < 0$	0.8	1.0	0.0	0.0

Equations101

\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = 0

\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = 0

F = ρ u ρ u^{2} + p ρ uv u (ρ e + p) G = ρ v ρ uv ρ v^{2} + p v (ρ e + p)

F = ρ u ρ u^{2} + p ρ uv u (ρ e + p) G = ρ v ρ uv ρ v^{2} + p v (ρ e + p)

F_{1} = u ρ ρ u ρ v \frac{1}{2} ρ (u^{2} + v^{2}) F_{2} = 0 p 0 \frac{γ}{γ - 1} p u G_{1} = v ρ ρ u ρ v \frac{1}{2} ρ (u^{2} + v^{2}) G_{2} = 00 p \frac{γ}{γ - 1} p v

F_{1} = u ρ ρ u ρ v \frac{1}{2} ρ (u^{2} + v^{2}) F_{2} = 0 p 0 \frac{γ}{γ - 1} p u G_{1} = v ρ ρ u ρ v \frac{1}{2} ρ (u^{2} + v^{2}) G_{2} = 00 p \frac{γ}{γ - 1} p v

\frac{\partial}{\partial t} \int_{Ω} U d Ω + \oint_{d Ω} (F + G) \cdot \hat{n} d s = 0

\frac{\partial}{\partial t} \int_{Ω} U d Ω + \oint_{d Ω} (F + G) \cdot \hat{n} d s = 0

\overset{ˉ}{U}_{i, j}^{n + 1} = \overset{ˉ}{U}_{i, j}^{n} + \frac{Δ t}{Δ x} (\overset{ˉ}{F}_{i - \frac{1}{2}, j}^{'} - \overset{ˉ}{F}_{i + \frac{1}{2}, j}^{'}) + \frac{Δ t}{Δ y} (\overset{ˉ}{G}_{i, j - \frac{1}{2}}^{'} - \overset{ˉ}{G}_{i, j + \frac{1}{2}}^{'})

\overset{ˉ}{U}_{i, j}^{n + 1} = \overset{ˉ}{U}_{i, j}^{n} + \frac{Δ t}{Δ x} (\overset{ˉ}{F}_{i - \frac{1}{2}, j}^{'} - \overset{ˉ}{F}_{i + \frac{1}{2}, j}^{'}) + \frac{Δ t}{Δ y} (\overset{ˉ}{G}_{i, j - \frac{1}{2}}^{'} - \overset{ˉ}{G}_{i, j + \frac{1}{2}}^{'})

\overset{ˉ}{F}_{i + \frac{1}{2}, j}^{'} = \frac{1}{6} F_{i + \frac{1}{2}, j + \frac{1}{2}}^{*} + \frac{4}{6} F_{i + \frac{1}{2}, j}^{mi d} + \frac{1}{6} F_{i + \frac{1}{2}, j - \frac{1}{2}}^{*}

\overset{ˉ}{F}_{i + \frac{1}{2}, j}^{'} = \frac{1}{6} F_{i + \frac{1}{2}, j + \frac{1}{2}}^{*} + \frac{4}{6} F_{i + \frac{1}{2}, j}^{mi d} + \frac{1}{6} F_{i + \frac{1}{2}, j - \frac{1}{2}}^{*}

F_{i + \frac{1}{2}}^{mi d} = F_{1}_{i + \frac{1}{2}, j}^{mi d} + F_{2}_{i + \frac{1}{2}, j}^{mi d}

F_{i + \frac{1}{2}}^{mi d} = F_{1}_{i + \frac{1}{2}, j}^{mi d} + F_{2}_{i + \frac{1}{2}, j}^{mi d}

F_{1}_{i + \frac{1}{2}, j}^{mi d} = M_{k} ρ ρ u ρ v \frac{1}{2} ρ (u^{2} + v^{2})_{k} a_{k}

F_{1}_{i + \frac{1}{2}, j}^{mi d} = M_{k} ρ ρ u ρ v \frac{1}{2} ρ (u^{2} + v^{2})_{k} a_{k}

k = {L \mbox i f \overset{u}{ˉ} \geq 0 R \mbox i f \overset{u}{ˉ} < 0

k = {L \mbox i f \overset{u}{ˉ} \geq 0 R \mbox i f \overset{u}{ˉ} < 0

M_{k} = ⎩ ⎨ ⎧ \frac{u ˉ}{u ˉ - S _{L}^{c}} \mbox i f \overset{u}{ˉ} \geq 0 \frac{u ˉ}{u ˉ - S _{R}^{c}} \mbox i f \overset{u}{ˉ} < 0

M_{k} = ⎩ ⎨ ⎧ \frac{u ˉ}{u ˉ - S _{L}^{c}} \mbox i f \overset{u}{ˉ} \geq 0 \frac{u ˉ}{u ˉ - S _{R}^{c}} \mbox i f \overset{u}{ˉ} < 0

a_{k} = {u_{L} - S_{L}^{c} \mbox i f \overset{u}{ˉ} \geq 0 u_{R} - S_{R}^{c} \mbox i f \overset{u}{ˉ} < 0

a_{k} = {u_{L} - S_{L}^{c} \mbox i f \overset{u}{ˉ} \geq 0 u_{R} - S_{R}^{c} \mbox i f \overset{u}{ˉ} < 0

F_{2}_{i + \frac{1}{2}, j}^{mi d} = \frac{S _{R}^{c}}{S _{R}^{c} - S _{L}^{c}} F_{2}_{L} - \frac{S _{L}^{c}}{S _{R}^{c} - S _{L}^{c}} F_{2}_{R} + \frac{S _{R}^{c} S _{L}^{c}}{S _{R}^{c} - S _{L}^{c}} (U_{R} - U_{L})

F_{2}_{i + \frac{1}{2}, j}^{mi d} = \frac{S _{R}^{c}}{S _{R}^{c} - S _{L}^{c}} F_{2}_{L} - \frac{S _{L}^{c}}{S _{R}^{c} - S _{L}^{c}} F_{2}_{R} + \frac{S _{R}^{c} S _{L}^{c}}{S _{R}^{c} - S _{L}^{c}} (U_{R} - U_{L})

F_{2}_{i + \frac{1}{2}, j}^{mi d} = \frac{1}{2} (F_{2 L} + F_{2 R}) + δ U_{2}

F_{2}_{i + \frac{1}{2}, j}^{mi d} = \frac{1}{2} (F_{2 L} + F_{2 R}) + δ U_{2}

δ U_{2} = \frac{S _{R}^{c} + S _{L}^{c}}{2 ( S _{R}^{c} - S _{L}^{c} )} (F_{2}_{L} - F_{2}_{R}) - \frac{S _{L}^{c} S _{R}^{c}}{S _{R}^{c} - S _{L}^{c}} ρ_{L} - ρ_{R} (ρ u)_{L} - (ρ u)_{R} (ρ v)_{L} - (ρ v)_{R} (ρ e)_{L} - (ρ e)_{R}

δ U_{2} = \frac{S _{R}^{c} + S _{L}^{c}}{2 ( S _{R}^{c} - S _{L}^{c} )} (F_{2}_{L} - F_{2}_{R}) - \frac{S _{L}^{c} S _{R}^{c}}{S _{R}^{c} - S _{L}^{c}} ρ_{L} - ρ_{R} (ρ u)_{L} - (ρ u)_{R} (ρ v)_{L} - (ρ v)_{R} (ρ e)_{L} - (ρ e)_{R}

δ U_{2} = \frac{S _{R}^{c} + S _{L}^{c}}{2 ( S _{R}^{c} - S _{L}^{c} )} (F_{2 L} - F_{2 R}) - \frac{S _{R}^{c} S _{L}^{c}}{a ˉ _{c}^{2} ( S _{R}^{c} - S _{L}^{c} )} p_{L} - p_{R} (p u)_{L} - (p u)_{R} (p v)_{L} - (p v)_{R} \frac{a ˉ _{c}^{2}}{2} (p_{L} - p_{R}) + \frac{1}{2} [(p q^{2})_{L} - (p q^{2})_{R}]

δ U_{2} = \frac{S _{R}^{c} + S _{L}^{c}}{2 ( S _{R}^{c} - S _{L}^{c} )} (F_{2 L} - F_{2 R}) - \frac{S _{R}^{c} S _{L}^{c}}{a ˉ _{c}^{2} ( S _{R}^{c} - S _{L}^{c} )} p_{L} - p_{R} (p u)_{L} - (p u)_{R} (p v)_{L} - (p v)_{R} \frac{a ˉ _{c}^{2}}{2} (p_{L} - p_{R}) + \frac{1}{2} [(p q^{2})_{L} - (p q^{2})_{R}]

S_{L}

S_{L}

S_{R}

u^{*}

u^{*}

c^{*}

S_{L}^{c} = - \overset{a}{ˉ}_{c} S_{R}^{c} = \overset{a}{ˉ}_{c}

S_{L}^{c} = - \overset{a}{ˉ}_{c} S_{R}^{c} = \overset{a}{ˉ}_{c}

F_{1}_{i + \frac{1}{2}, j + \frac{1}{2}}^{*} = \overset{u}{ˉ} \frac{S _{U} ρ ρ u ρ v \frac{1}{2} ρ ( u ^{2} + v ^{2} ) _{k_{1}} - S _{D} ρ ρ u ρ v \frac{1}{2} ρ ( u ^{2} + v ^{2} ) _{k_{2}}}{S _{U} - S _{D}}

F_{1}_{i + \frac{1}{2}, j + \frac{1}{2}}^{*} = \overset{u}{ˉ} \frac{S _{U} ρ ρ u ρ v \frac{1}{2} ρ ( u ^{2} + v ^{2} ) _{k_{1}} - S _{D} ρ ρ u ρ v \frac{1}{2} ρ ( u ^{2} + v ^{2} ) _{k_{2}}}{S _{U} - S _{D}}

\overset{u}{ˉ} = \frac{u _{LU} S _{U} - u _{L D} S _{D} + u _{R U} S _{U} - u _{R D} S _{D}}{2 ( S _{U} - S _{D} )}

\overset{u}{ˉ} = \frac{u _{LU} S _{U} - u _{L D} S _{D} + u _{R U} S _{U} - u _{R D} S _{D}}{2 ( S _{U} - S _{D} )}

G_{1}_{i + \frac{1}{2}, j + \frac{1}{2}}^{*} = \overset{v}{ˉ} \frac{S _{R} ρ ρ u ρ v \frac{1}{2} ρ ( u ^{2} + v ^{2} ) _{k_{1}} - S _{L} ρ ρ u ρ v \frac{1}{2} ρ ( u ^{2} + v ^{2} ) _{k_{2}}}{S _{R} - S _{L}}

G_{1}_{i + \frac{1}{2}, j + \frac{1}{2}}^{*} = \overset{v}{ˉ} \frac{S _{R} ρ ρ u ρ v \frac{1}{2} ρ ( u ^{2} + v ^{2} ) _{k_{1}} - S _{L} ρ ρ u ρ v \frac{1}{2} ρ ( u ^{2} + v ^{2} ) _{k_{2}}}{S _{R} - S _{L}}

\overset{v}{ˉ} = \frac{v _{R U} S _{R} - v _{LU} S _{L} + v _{R D} S _{R} - v _{L D} S _{L}}{2 ( S _{R} - S _{L} )}

\overset{v}{ˉ} = \frac{v _{R U} S _{R} - v _{LU} S _{L} + v _{R D} S _{R} - v _{L D} S _{L}}{2 ( S _{R} - S _{L} )}

\overset{u}{ˉ} = \frac{u _{LU} S _{U} - u _{L D} S _{D}}{S _{U} - S _{D}}

\overset{u}{ˉ} = \frac{u _{LU} S _{U} - u _{L D} S _{D}}{S _{U} - S _{D}}

\overset{u}{ˉ} = \frac{u _{R U} S _{U} - u _{R D} S _{D}}{S _{U} - S _{D}}

\overset{u}{ˉ} = \frac{u _{R U} S _{U} - u _{R D} S _{D}}{S _{U} - S _{D}}

\overset{v}{ˉ} = \frac{v _{R D} S _{R} - v _{L D} S _{L}}{S _{R} - S _{L}}

\overset{v}{ˉ} = \frac{v _{R D} S _{R} - v _{L D} S _{L}}{S _{R} - S _{L}}

\overset{v}{ˉ} = \frac{v _{R U} S _{R} - v _{LU} S _{L}}{S _{R} - S _{L}}

\overset{v}{ˉ} = \frac{v _{R U} S _{R} - v _{LU} S _{L}}{S _{R} - S _{L}}

F_{2}_{i + \frac{1}{2}, j + \frac{1}{2}}^{*}

F_{2}_{i + \frac{1}{2}, j + \frac{1}{2}}^{*}

- 2 \frac{S _{R} S _{L}}{( S _{R} - S _{L} ) ( S _{U} - S _{D} )} (G_{2}_{R U} - G_{2}_{LU} + G_{2}_{L D} - G_{2}_{R D})

+ \frac{S _{R} S _{L}}{( S _{R} - S _{L} ) ( S _{U} - S _{D} )} (S_{U} (U_{R U} - U_{LU}) - S_{D} (U_{R D} - U_{L D}))

G_{2}_{i + \frac{1}{2}, j + \frac{1}{2}}^{*}

G_{2}_{i + \frac{1}{2}, j + \frac{1}{2}}^{*}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsComputational Fluid Dynamics and Aerodynamics · Gas Dynamics and Kinetic Theory · Fluid Dynamics and Turbulent Flows

Full text

An accurate and robust genuinely multidimensional Riemann solver 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, shock stable Riemann solver (GM-K-CUSP-X) for Euler equations of gasdynamics is developed. The convective and pressure components of the Euler system are seperated following the Toro-Vazquez type PDE flux splitting [Toro et al, 2012]. Upwind discretization of these components are achieved using the framework of Mandal et al [Mandal et al, 2015]. The robustness of the scheme is studied on a few two dimensional test problems. The results demonstrate the efficacy of the scheme over the corresponding conventional two state version of the solver. Results from two classic strong shock test cases associated with the infamous Carbuncle phenomenon, indicate that the present solver is completely free of any such numerical instabilities albeit possessing contact resolution abilities.Such a finding emphasizes the preexisting notion about the positive effects that multidimensional flow modeling may have towards curing of shock instabilities.

Keywords : Genuinely multidimensional, Riemann solver, Contact preserving, Convective Pressure Split, numerical shock instability, Carbuncle phenomenon.

1 Introduction

Past few decades have witnessed commendable advancement in the computation of high speed flows. Birth of the upwind schemes for finite volume methods is a milestone in this regard. Among these characteristics based schemes, those based on the solution of two state ‘Riemann problems’at cell interface are the most popular class of schemes 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 exact solver was critized for being prohibitively expensive because of the iterative technique involved [5]. To mitigate this shortcoming, a class of approximate Riemann solvers have been proposed and continue to be developed. Details about some famous approximate Riemann solvers can be found in [5, 6, 7, 8]. Such algorithms have brought forth an 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 accurate 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 features are oblique to the grid still pose a major challenge. This is mainly because the wave system in such problems possess infinitely many propagation direction as compared to the limited ones in a one-dimensional problem. Further vorticity wave enters the formulation and has to be dealt with appropriately [9]. Currently, the standard way of extending the one-dimensional schemes to multidimensions is through a dimensionally split operator approach motivated in [10] where local one dimensional Riemann problems are solved in a direction normal to each cell interface. Such a strategy has been criticised because of its unnatural selection of the interface normal as the only wave propagation direction even for problems with infinitely many propagation directions. This lacuna is attributed to cause for example, pressure disturbances across a grid oblique shear wave [11].

Another disturbing problem encountered by dimensionally split approximate Riemann solvers is the occurrence of various forms of numerical shock instabilities when simulating strong normal shocks. A catalouge of many such failures can be found in [12]. It has been observed that only those schemes that have exact shear wave resolution abilities are known to produce these instabilities. There is increasing evidence [13, 14] to consider that such failures 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 characteristic fields while introducing appropriate dissipation along necessary directions.

Probably the earliest attempt at creating a multidimensional Riemann solver was by Raithby [15] who suggested discretizing the convective terms along flow dominant directions. Davis [16] and others [17, 18, 19] introduced the idea of rotated Riemann solvers that involves identifying grid oblique shock orientations and solving Riemann problems across them.

Roe and others [9, 20] suggested using extrapolated data from the left and right states across a Riemann solver to construct simple wave solutions for a multidimensional Riemann problem. This technique is not generally valid for 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 incorporate diagonal cell contirbutions to interface fluxes in a predictor corrector framework.

Leveque [22] proposed a multidimensional scheme in which cross derivative terms that constitute the transverse contributions are included by interpreting the usual one dimensional Gudonov method in a fluctuation splitting framework.

Ren et al [23] constructed an operator split predictor corrector scheme based on CTU and Leveque’s wave propagation method for both Euler and Navier Stokes systems. While the predictor step solves linearized Euler equations in characteristic variables, the corrector step adds viscous contributions.

Wendroff [24] introduced a multistate Riemann solver that attempted the extension of one dimentional HLLE scheme [6] into several dimensions. The scheme used an expensive nine point stencil and suffered unacceptable dissipation due to unnatural wavespeed selection [26].

Another multistate Riemann solver is from the work of Brio et al [25]. Multidimensional effects are incorporated by adding correction terms to the standard face normal fluxes at every computational cell interface. These correction terms are obtained by solving a three state Riemann problem using Roe’s FDS at the corners of the respective cells.

In a similar spirit, Balsara [26] presented a generic multidimensional HLLE solver (GM-HLLE) with simple closed form expressions that can be extended to any hyperbolic system. This method too relies on construction of multidimensional correction terms like [25] wherein these terms are obtained using a four state HLLE Riemann solver at interface corners.

By building upon the wave model introduced in GM-HLLE, Balsara [27] further proposed a multidimensional HLLC solver for Euler and MHD systems. To deal with contact discontinuities in two dimensions, a set of twelve possible contact orientations on a given cell are included in the wave model. Although posessing closed form expressions 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 variables 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 on Zha-Bilgen [29] way of splitting Euler fluxes. The scheme consists of a wave speed averaged upwinding for the convective part and GM-HLLE type discretization for the pressure part. This scheme is basically a multidimensional extension of HLL-CPS strategy [30] that uses a HLL type discretization for distinctly treating convective and pressure flux parts. The splitting of full Euler flux into convective and pressure parts can be achieved by adopting AUSM type or Zha-Bilgen type PDE splitting [8, 29]. Diffusion control is achieved by careful tuning of the dissipation vector of the pressure flux. This renders the scheme with stationary and moving contact preserving abilities in addition to capability of surviving the most stringent test problems. Toro has corroborated such a method by showing that exact contact preserving ability can be incorporated into convective-pressure split framework by discretizing the pressure fluxes using a Riemann solver [1]. Surprisingly, although contact preserving, HLL-CPS is 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. The scheme uses exactly similar discretization as HLL-CPS for both convective and pressure fluxes but differs only in the method of splitting the Euler fluxes into these components; instead of using Zha-Bilgen or AUSM way of PDE 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 suffer from numerical shock instabilites.

In this paper, a new genuinely multidimensional scheme based on the conventional K-CUSP-X is proposed. A multidimensional correction term similar to [25] and [26] is constructed by solving appropriate four state Riemann problem at the corners of each interface. The multidimensional terms are incorporated such that the final fluxes can 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 purpose of this paper to demonstrate that such a construction is not only as accurate as GM-HLL-CPS-Z, but also cures the numerical shock instability that plagued the corresponding two state conventional K-CUSP-X Riemann solver. The positive effect 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 completely removed by extending the original two state solver into its genuinely multidimensional variant. Such a finding corroborates the pre existing opinion that multidimensional dissipation acts as a reliable cure for numerical shock instability problems and incentivizes the need for extending the accurate Riemann solvers in literature into their robust genuinely multidimensional counterparts.

This paper is organized as follows. In the next section, the governing equations and the type of PDE flux splitting used is detailed. In Section 3, the second order version of the new formulation is described. In Section 4, results for some complex flow problems like double Mach reflection and multidimensional Riemann problem are discussed. Further, two classic shock instability test problems, the odd-even decoupling problem and standing shock problem are used to study instability behaviour of the newly developed scheme. Section 5 contains some concluding remarks.

2 Preliminaries

2.1 Governing equation

Consider the two-dimensional Euler equations in differential form

[TABLE]

where $\mathbf{U}=[\rho$ $\rho u$ $\rho v$ $\rho e]^{T}$ is the vector of conserved quantities. $\mathbf{F}$ and $\mathbf{G}$ are the flux vectors in the x and y directions respectively given by

[TABLE]

In the present work, the above flux vectors are split into corresponding convective and pressure parts following the approach of Toro-Vazquez [1] . Using ideal gas law, the split flux vectors can be written as

[TABLE]

where $\mathbf{F_{1}}$ and $\mathbf{G_{1}}$ are the convective fluxes while $\mathbf{F_{2}}$ and $\mathbf{G_{2}}$ are the pressure fluxes. $\gamma$ represents the ratio of specific heat capacities. As suggested by Toro et al [1] the convective fluxes $\mathbf{F_{1}}$ and $\mathbf{G_{1}}$ can be interpreted as simple advection of mass, momentum and kinetic energy along the x and y directions respectively, while the pressure fluxes $\mathbf{F_{2}}$ and $\mathbf{G_{2}}$ are interpreted to be sonic impulses that spread in all directions with reference to these convecting particles. This type of PDE splitting primarily differs from those that already exist in the literature like [29] and [8], in terms of the quantity concerning energy that is being advected. Such differences may have strong bearing on the robustness of these schemes.

2.2 Finite volume discretization

Consider the integral form of equation (1)

[TABLE]

where $\hat{\mathbf{n}}=(n_{x},n_{y})$ is the unit vector along the face normal. Consider a Cartesian cell of area $\Delta x\times\Delta y$ as shown in Figure 1. The finite volume discretization for the cell ( $i,j$ ) can be written as

[TABLE]

where n is the time level, i and j are cell indices. The $\bar{(.)}$ quantities depict the respective averaged quantities. The total fluxes at the cell interfaces ( $i+\frac{1}{2},j$ ) and ( $i-\frac{1}{2},j$ ) are denoted as $\bar{\mathbf{F}}^{\prime}_{i+\frac{1}{2},j}$ and $\bar{\mathbf{F}}^{\prime}_{i-\frac{1}{2},j}$ respectively. Similarly $\bar{\mathbf{G}}^{\prime}_{i,j+\frac{1}{2}}$ and $\bar{\mathbf{G}}^{\prime}_{i,j-\frac{1}{2}}$ are the y directional total fluxes in corresponding y interfaces.

3 Formulation

In the proposed scheme the interface flux will comprise of both a two state conventional Riemann flux and a multidimensional flux. These multidimensional fluxes are sought from the solution of four state Riemann problem that occurs at each corner of a Cartesian cell. To see this more clearly, consider a typical cell ( $i,j$ ) as shown in Figure 2. The four constant states that come together at corner $c_{1}$ at $t=0$ and forms a two-dimensional Riemann problem are marked as LU (Left-Upper), LD (Left-Down), RU (Right-Upper) and RD (Right-Down). Although realistically, the waves pertaining to this multistate Riemann problem will travel in infinitely many directions and swap a curvilinear area in space at any later time $t=T$ , for sake of simplicity, the simple wave model consisting of only four waves as introduced in [26] is used in this work. Accordingly the wave propagation at the corners is assumed to span, at any time T, a rectangular region. Figure 3 represents a three dimensional view of these waves as they evolve in time where the shaded region depicts the domain of influence of this multidimensional Riemann problem.

In principle, the Equation 1 can be integrated along the space-time volume of Figure 3 appropriately to obtain the time averaged multidimensional fluxes $\mathbf{F}^{*}$ in x-direction and $\mathbf{G}^{*}$ in y-direction. For a typical cell interface ( $i+\frac{1}{2},j$ ) the total flux will be an ensemble of the two multidimensional fluxes from corners $c_{1}$ and $c_{4}$ denoted as $\mathbf{F}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ , $\mathbf{F}^{*}_{i+\frac{1}{2},j-\frac{1}{2}}$ and the single mid-point flux $\mathbf{F}_{i+\frac{1}{2},j}^{mid}$ . This is shown in Figure 4. A conservative ensemble of the corner and mid point fluxes are achieved by using Simpson’s rule of intergration [26] along the interface as given by,

[TABLE]

3.1 Evaluation of flux at the midpoint of the cell interface

This section deals with determination of the two state Riemann flux at the mid point of a typical interface $(i+\frac{1}{2},j)$ denoted as ${\mathbf{F}_{i+\frac{1}{2},j}^{mid}}$ . Following the Convective Pressure Split (CPS) philosophy the total flux at this interface is first split into convective and pressure parts. This work uses the Toro-Vazquez type flux splitting as mentioned in Equation (3),

[TABLE]

These convective and pressure parts are discretized independently following the original HLL-CPS strategy [30]. It may be noted that the interface $(i+\frac{1}{2},j)$ admits only the x-directional Riemann flux ${\mathbf{F}_{i+\frac{1}{2},j}^{mid}}$ while the y-directional flux, ${\mathbf{G}_{i+\frac{1}{2},j}^{mid}}$ is zero on it.

3.1.1 Evaluation of convective flux at the midpoint of the cell interface ( ${\mathbf{F}_{1}}_{i+\frac{1}{2},j}^{mid}$ )

The upwind discretization of the convective flux denoted by ${\mathbf{F}_{1}}_{i+\frac{1}{2},j}^{mid}$ at the mid point of the cell interface ( $i+\frac{1}{2},j$ ), following the strategy of HLL-CPS method [30], is given by,

[TABLE]

where, the average local x-directional velocity at the interface is taken as $\bar{u}=\displaystyle\frac{u_{L}+u_{R}}{2}$ . Depending on the sign of the average local x-directional velocity, left (L) or right (R) states are selected for upwinding. $S_{L}^{c}$ and $S_{R}^{c}$ are carefully selected wave speeds which are discussed in section 3.1.3.

3.1.2 Evaluation of pressure flux at the midpoint of the cell interface ( ${\mathbf{F}_{2}}_{i+\frac{1}{2},j}^{mid}$ )

The pressure flux at the midpoint of the cell interface is obtained by applying a HLL type discretization of the pressure flux vector [30].

[TABLE]

The above equation can be rewritten as

[TABLE]

where $\mathbf{F}_{2L}$ and $\mathbf{F}_{2R}$ are the left and right side x-directional pressure fluxes normal to the interface $(i+\frac{1}{2},j)$ and $\delta\mathbf{U}_{2}$ is the numerical diffusion given by

[TABLE]

It is clear from Equation 14 that the second term will give rise to numerical diffusion across a contact. Thus in order to capture the contact discontinuity accurately, density terms in $\delta\mathbf{U}_{2}$ are replaced by the pressure terms by using isentropic assumption $\bar{a}_{c}^{2}=\displaystyle\frac{\delta p}{\delta\rho}$ as described in the reference [30].

[TABLE]

where $\bar{a}_{c}$ is the average speed of sound at the interface given by $\bar{a}_{c}=\displaystyle\frac{a_{L}+a_{R}}{2}$ and $q^{2}=u^{2}+v^{2}$ is twice the local kinetic energy per unit mass.

3.1.3 Selection of wave speeds for the 1D Riemann problem at the interface ( $\mathbf{i+\frac{1}{2},j}$ )

The wave speeds are selected according to conventional HLL-CPS method [30]

[TABLE]

where $u^{*}$ and $c^{*}$ are given by [30],

[TABLE]

Supersonic conditions are taken care by including ’0’ in the above expressions. For a stationary flow the wave speeds are modified as

[TABLE]

It is to be noted that midpoint fluxes at other interfaces can be obtained in the similar manner.

3.2 Evaluation of flux at the corner of the cell interface

This section will detail how to evaluate the fluxes $\mathbf{F}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ and $\mathbf{G}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ that results from the interaction of four Riemann states at a representative corner $c_{1}$ . These fluxes forms the multidimensional component of the interface fluxes and are shown in Figure 5. Once again, resorting to the (CPS) philosophy, the total flux at this corner is split into convective ( $\mathbf{{F_{1}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ and $\mathbf{{G_{1}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ ) and pressure fluxes ( $\mathbf{{F_{2}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ and $\mathbf{{G_{2}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ ) originating at the corner as per Equation 3. Analogous to the procedure for evaluation of the two state mid point flux at the interface, the split convective and the pressure parts at the corners will also be evaluated using different upwind strategies.

3.2.1 Evaluation of convective flux at the corner $c_{1}$ of the cell interface ( $\mathbf{{F_{1}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ )

Following [2], the x-directional convective flux is evaluated as,

[TABLE]

where $\bar{u}$ is the wave speed averaged x-directional local fluid speed at the corner defined as,

[TABLE]

Based on the direction of the average x- directional flow, the upwind states $k_{1},k_{2}$ are chosen as,

If $\bar{u}>0$ , $k_{1}=LU$ and $k_{2}=LD$ (i.e upwinding is done from the left states).

If $\bar{u}<0$ , $k_{1}=RU$ and $k_{2}=RD$ (i.e upwinding is done from the right states).

In similar spirit, the y directional convective flux is evaluated as,

[TABLE]

where $\bar{v}$ is the the wave speed averaged y-directional local fluid speed at the corner defined as,

[TABLE]

The states $k_{1},k_{2}$ are chosen accordingly as,

If $\bar{v}>0$ , $k_{1}=RD$ and $k_{2}=LD$ (i.e upwinding is done from the down states).

If $\bar{v}<0$ , $k_{1}=RU$ and $k_{2}=LU$ (i.e upwinding is done from the up states).

Since the above strategy is developed for a subsonic case, slight modification $\bar{u}$ and $\bar{v}$ is done to extend the above formulation to supersonic cases:

If the flow is supersonic in positive x-direction, then at the corner upwinding is done from left states. Therefore for the evaluation of x-directional convective flux $\bar{u}$ is taken as

[TABLE]

If the flow is supersonic in negative x-direction the upwinding is done from right states. Therefore for the evaluation of x-directional convective flux $\bar{u}$ is taken as

[TABLE]

If the flow is supersonic in positive y-direction the upwinding is done from down states. Therefore for the evaluation of y-directional convective flux $\bar{v}$ is taken as

[TABLE]

If the flow is supersonic in negative y-direction the upwinding is done from upper states. Therefore for the evaluation of y-directional convective flux $\bar{v}$ is taken as

[TABLE]

3.2.2 Evaluation of pressure flux at the corner $c_{1}$ of the cell interface ( $\mathbf{{F_{2}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ )

Following [2], the x-directional convective flux is evaluated as,

[TABLE]

Similarly the y-directional pressure flux is evaluated as,

[TABLE]

The above equations can be rewritten as,

[TABLE]

where

[TABLE]

The $\delta\mathbf{U_{2x}}$ and $\delta\mathbf{U_{2y}}$ terms in Equations (25,26), are the numerical diffusion terms in x and y-directions respectively. To remove the numerical dissipation across a contact wave, the dissipation terms in x and y-directions are remodeled using isentropic expression as,

[TABLE]

where $e^{*}_{k}$ is given as

[TABLE]

and $\bar{a}$ as

[TABLE]

The flux contributions due to the genuinely multidimensional Riemann problem at the other corner $c_{4}$ of the interface can be obtained in a similar manner. Most importantly, it must be noted that while $\mathbf{F}^{*}_{i+\frac{1}{2},j\pm\frac{1}{2}}$ contributes to the total interface flux at the interface $(i+\frac{1}{2},j)$ , $\mathbf{G}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ and $\mathbf{G}^{*}_{i+\frac{1}{2},j-\frac{1}{2}}$ contributes to the total interface flux at the interfaces $(i,j+\frac{1}{2})$ and $(i,j-\frac{1}{2})$ respectively.

3.3 Selection of wave speeds for the multidimensional Riemann problem at the corners

As previously mentioned, the present work adopts a simple wave model as proposed by [26] to represent the waves emerging from the four state Riemann problem at the corners of every interface. A top view of the area swept by 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 of influence of the four state Riemann problem at corner $c_{1}$ at time T on x-y plane.

An estimate for these wavespeeds can be obtained as,

[TABLE]

where typically,

$\lambda_{x}^{1}(\mathbf{U}_{k})$ denote smallest x-directional wave speed in the state $\mathbf{U}_{k}$ ,

$\lambda_{x}^{N}(\mathbf{U}_{k})$ denote largest x-directional wave speed in the state $\mathbf{U}_{k}$ ,

$\tilde{\lambda}_{x}^{1}(\mathbf{U}_{k},\mathbf{U}_{l})$ denotes smallest x-directional wave speed from Roe averaged state between $\mathbf{U}_{k}$ and $\mathbf{U}_{l}$ ,

$\tilde{\lambda}_{x}^{N}(\mathbf{U}_{k},\mathbf{U}_{l})$ denotes largest x-directional wave speed from Roe averaged state between $\mathbf{U}_{k}$ and $\mathbf{U}_{l}$ such that $k,l\in\{LU,LD,RU,RD\}$ and $k\neq l$ .

If $\bar{u}=0$ and $\bar{v}\neq 0$ , then x-directional wave speeds are modified as $S_{L}=-\bar{a}$ and $S_{R}=\bar{a}$ .

If $\bar{u}\neq 0$ and $\bar{v}=0$ , then y-directional wave speeds are modified as $S_{D}=-\bar{a}$ and $S_{U}=\bar{a}$ .

If $\bar{u}=0$ and $\bar{v}=0$ , then wave speeds are modified as $S_{L}=-\bar{a}$ , $S_{R}=\bar{a}$ , $S_{D}=-\bar{a}$ and $S_{U}=\bar{a}$ .

It should be noted that zero has been added in the above expressions in order to ensure fully one sided flux in supersonic flow.

4 Results

4.1 Isentropic vortex problem

A second order accurate version of the present solver is developed using the SDWLS strategy [34] whose details are omitted here for brevity. Formal order of accuracy of the second order version of GM-K-CUSP-X is investigated using the isentropic vortex problem [2]. The problem involves an isentropic vortex centered initially at (0,0) and made to traverse the domain diagonally under a periodic boundary condition. A Cartesian domain of size [-5,5] $\times$ [-5,5] is used. The initial conditions consists of a unperturbed state given by $(\rho,p,u,v)$ = $(1.0,1.0,1.0,1.0)$ . The temperature is defined as $T=\frac{p}{\rho}$ . A perturbation is added to the flow as,

[TABLE]

Here, $\epsilon$ 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_{1}$ and $L_{\infty}$ norms of the density variable. The formal order of accuracy $O$ is obtained by the formula,

[TABLE]

where, $\eta_{1}$ and $\eta_{2}$ depict consecutive norms ( $L_{1}$ or $L_{\infty}$ ) for progressively refined grids with dimensions $\Delta x_{1}$ and $\Delta x_{2}$ respectively. The results are tabulated in table 1.

It is observed from the analysis that GM-K-CUSP-X scheme is able to achieve second order accuracy on sufficiently refined grids in both $L_{1}$ and $L_{\infty}$ norms.

4.2 Two dimensional Riemann problems

Two dimensional Riemann problems provides an excellent test case to assess the qualitative improvement of genuinely multidimensional formulation of GM-K-CUSP-X over the corresponding conventional two state K-CUSP-X scheme. Second order versions of both solvers are used for these test cases. The first two dimensional Riemann problem investigated in the present work consists of a double Mach reflection and an oblique shock wave propagating at an angle to the grid. The initial conditions of the 2D Riemann problem are given by [2],

A Cartesian grid of size $2000\times 2000$ spanning $[-1,1]\times[-1,1]$ is used. The solution is evolved up to a time of 1.05 units with a CFL number of 0.95 as suggested in reference [27]. The density contours at the final time are shown in Figure 7. Result for K-CUSP-X solver under identical conditions is given for comparison. It is observed that 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.

The second Riemann problem consist evolution of two weak shock waves and two contact waves and can be set using conditions [2],

A 2000X2000 Cartesian grid spanning [-1,1] $\times$ [-1,1] domain is used. CFL number of 0.95 is used and the solution is evolved for a time of 0.5 units. The results obtained are represented using iso density contours. For comparison, results obtained under identical conditions by corresponding two state conventional K-CUSP-X is also provided 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.

4.3 Double Mach reflection problem

This problems deals with a Mach 10 shock inclined at $60^{o}$ with the x-axis and propagating downstream of a rectangular duct of size [0,4] $\times$ [0,1] and interacting with a reflective bottom boundary wall. Detailed initial and boundary conditions can be found in [2]. The problem is solved on a Cartesian mesh of size 1920 $\times$ 480 with a CFL 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 same conditions. It is visible from the result that the present solver is able to resolve the complex shock structures crisply and also the slipping contact line emerging from the triple point.

5 Numerical shock instability tests

One of the objectives of the present work is to demonstrate the effect of genuinely multidimensional formulation on the numerical shock instability characteristics of a Riemann solver. Two standard test cases namely, the odd-even decoupling problem [12] and the standing shock problem [32] will be used to carry out the investigations. Since shock instability is most explicitly observed in the first order simulations, first order version of solvers will be used for these tests.

5.1 Odd-Even decoupling problem

Odd-even decoupling test consists of a M = 6 shock propagating down a computational duct whose centerline grid is slightly perturbed to induce random oscillations into the initial conditions [12]. It has been long argued that odd-even decoupling and the Carbuncle have the same origin [35]. Much alike the Carbuncle phenomenon, the moving shock profile in the odd-even decoupling problem also deteriorates over time (producing the recognizable bulge at the centerline) and pollutes the after shock flow field. The original K-CUSP-X scheme was found to have slight after shock perturbations in this test. In comparison, the behavior of GM-K-CUSP-X scheme on this problem is shown in 10b after the solution has evolved for t = 140 units. It is clearly observed that the genuinely multidimensional extension is able to preserve the shock profile without any instabilities.

5.2 Standing shock instability

To clearly differentiate the behavior of K-CUSP-X and GM-K-CUSP-X schemes, a standing shock instability test case is used. The details of the test are available in [32]. As seen in the figure 11a, K-CUSP-X fails miserably in this test case. However, from figure 11b it is very evident that GM-K-CUSP-X scheme is free of this instability.

The 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 averaging that strives to accurately model the underlying multidimensional wave evolution phenomenon. Specifically, it is easy to note that Equations (20) and (20) defines wave averaged convection velocities in x and y directions that depends on all the adjoining states at the corners. Such a formulation admits information from transverse cells which would contribute to dissipation that would help supressing the instabilities. A similar observation can be found from Equations (25), (26) which describe the discretization of the pressure terms. For example, in Equation 25 for multidimensional flux $\mathbf{F_{2}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ , the third term on right hand side consist of terms $\mathbf{{G_{2}}}_{RU}$ , $\mathbf{{G_{2}}}_{LU}$ , $\mathbf{{G_{2}}}_{LD}$ and $\mathbf{{G_{2}}}_{RD}$ clearly indicating the influence of y-directional flux terms in the evaluation of x-directional flux. Equation 26 depicts a vice versa situation for the flux $\mathbf{G_{2}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ . These multidimensional coupling terms along with the wave averaged convective terms may be providing the additional cross dissipation that damps any unprecedented growth in instabilities thereby making the present scheme immune to shock instabilities. These equations reveal that the discretization of convective and pressure fluxes in the corner of interfaces using GM-HLLE strategy may have a positive impact in curing such instabilities.

6 Conclusions

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-pressure fluxes are discretized following [30]. While the convective fluxes are upwinded based on appropriate wave speeds that emerge from the interacting states, the pressure fluxes are treated in a HLLE framework. Restoration of stationary contact preservation ability is improved by explicitly reducing the numerical dissipation in the pressure flux 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 and mushroom cap structure in two dimensional Riemann problem and complex shock interactions in double Mach reflection problem are well resolved. Further, the present genuinely multidimensional solver is able to mitigate various 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 most promising methods to cure shock instabilities. Due to the simplicity of the formulation, the present solver can be easily extended to unstructured framework and three dimensional problems in principle.

Bibliography35

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] E.F. Toro, V. Cendon Flux splitting schemes for the Euler equations, Comput Fluids, 70(2012) 1-12.
2[2] J.C. Mandal , V. Sharma A genuinely multidimensional convective pressure flux split Riemann solver for Euler equations, J. Comput. Phys., 297(2015) 669-688.
3[3] J. C Mandal, J. Subramanian, On the link between weighted least-squares and limiters used in higher-order reconstructions for finite volume computations of hyperbolic equations, Appl Numer Math, 58(2008) 705-725.
4[4] S. K. Godunov, A finite difference method for the computation of discontinuous solutions of the equations of fluid dynamics, Mat. Sb., 47(1959) 357-393.
5[5] P. L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys., 43(1981) 357-372.
6[6] A. Harten, P. D. Lax, B. van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev., 25(1983) 35-61.
7[7] E. F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann Solver, Shock Waves, 4(1994) 25-34
8[8] M. S. Liou, C. J. Steffen, A new flux splitting scheme, J. Comput. Phys., 107(1993) 23-39.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Taxonomy

An accurate and robust genuinely multidimensional Riemann solver for Euler equations based on TV flux splitting

Abstract

1 Introduction

2 Preliminaries

2.1 Governing equation

2.2 Finite volume discretization

3 Formulation

3.1 Evaluation of flux at the midpoint of the cell interface

3.1.1 Evaluation of convective flux at the midpoint of the cell interface (F1i+12,jmid{\mathbf{F}_{1}}_{i+\frac{1}{2},j}^{mid}F1​i+21​,jmid​)

3.1.2 Evaluation of pressure flux at the midpoint of the cell interface (F2i+12,jmid{\mathbf{F}_{2}}_{i+\frac{1}{2},j}^{mid}F2​i+21​,jmid​)

3.1.3 Selection of wave speeds for the 1D Riemann problem at the interface (i+12,j\mathbf{i+\frac{1}{2},j}i+21​,j)

3.2 Evaluation of flux at the corner of the cell interface

3.2.1 Evaluation of convective flux at the corner c1c_{1}c1​ of the cell interface (F1i+12,j+12∗\mathbf{{F_{1}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}F1​i+21​,j+21​∗​)

3.2.2 Evaluation of pressure flux at the corner c1c_{1}c1​ of the cell interface (F2i+12,j+12∗\mathbf{{F_{2}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}F2​i+21​,j+21​∗​)

3.3 Selection of wave speeds for the multidimensional Riemann problem at the corners

4 Results

4.1 Isentropic vortex problem

4.2 Two dimensional Riemann problems

4.3 Double Mach reflection problem

5 Numerical shock instability tests

5.1 Odd-Even decoupling problem

5.2 Standing shock instability

6 Conclusions

3.1.1 Evaluation of convective flux at the midpoint of the cell interface ( ${\mathbf{F}_{1}}_{i+\frac{1}{2},j}^{mid}$ )

3.1.2 Evaluation of pressure flux at the midpoint of the cell interface ( ${\mathbf{F}_{2}}_{i+\frac{1}{2},j}^{mid}$ )

3.1.3 Selection of wave speeds for the 1D Riemann problem at the interface ( $\mathbf{i+\frac{1}{2},j}$ )

3.2.1 Evaluation of convective flux at the corner $c_{1}$ of the cell interface ( $\mathbf{{F_{1}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ )

3.2.2 Evaluation of pressure flux at the corner $c_{1}$ of the cell interface ( $\mathbf{{F_{2}}}^{*}_{i+\frac{1}{2},j+\frac{1}{2}}$ )