A new Roe-type scheme for all speeds

Computers and Fluids 121 (2015) 11–25

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

A new Roe-type scheme for all speeds Feng Qua,b, Chao Yanb,∗, Di Sunb, Zhenhua Jiangb a b

Institute of Manned Space System Engineering, China Academy of Space Technology, Beijing 100094, China National Laboratory for CFD, Beihang University, Beijing 100191, China

a r t i c l e

i n f o

Article history: Received 15 July 2014 Revised 26 March 2015 Accepted 6 July 2015 Available online 13 July 2015 Keywords: Shock anomaly Roe Low speeds All speeds RoeMAS Computational fluid dynamics

a b s t r a c t The all-speed Roe schemes still encounter the problem that the original Roe scheme does at high speeds because they are just improved in terms of low speeds. To overcome this problem, we propose a new scheme called RoeMAS (Roe Modified for All Speeds) in this paper. It changes the original Roe scheme’s coefficients DP and DρU to the order of O(c−1 ) and O(c0 ) to improve the level of accuracy at low speeds. To be robust against the shock anomaly and the expansion shock at high speeds, it controls the coefficients of the density perturbation in the numerical mass flux instead of controlling those of the pressure perturbation like the RoeM-type schemes. Moreover, it changes the non-linear wave speeds in rarefactions’ simulations to avoid using an entropy fix which is not satisfied in engineering area. Various numerical tests show that RoeMAS satisfies the following attractive properties independent of any tuning coefficient: (1) robustness against the shock anomaly and high discontinuity’s resolution (2) free from the expansion shock’s appearance (3) low dissipation and high resolution at low speeds. These properties suggest that RoeMAS is promising to be widely used to simulate flows of all speeds accurately and efficiently. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Since being proposed by Roe [1], the original Roe scheme has been widely used due to its high discontinuity and boundary resolutions [2]. However, it suffers many problems in the simulation for both cases of low speed and high speed. At low speed cases, the original Roe scheme encounters the problem of large disparity between the fluid speed and the acoustic speed, which leads to difficult or stalled convergence [3–5]. To overcome this defect, preconditioning matrix methods are widely used [6–11]. However, they need at least one tuning coefficient and how to define this coefficient depends on the user’s experience. In recent years, some improved Roe schemes for all speeds [12–18] are proposed to avoid using such tuning coefficients. However, both the preconditioning methods and the all-speed Roe schemes are inclined to encounter the original Roe scheme’s shortcomings at high speeds. For example, they suffer from the shock-induced anomaly [19], such as the famous carbuncle phenomena [20], in high speeds’ simulations. Also, they need to employ an entropy fix [21,22] in rarefactions’ simulations. But the entropy fix needs a problem-dependent coefficient. Improper coefficient’s definition may deteriorate the resolution [23].

∗

Corresponding author. E-mail address: [email protected] (C. Yan).

http://dx.doi.org/10.1016/j.compfluid.2015.07.007 0045-7930/© 2015 Elsevier B.V. All rights reserved.

Moreover, our test cases in this paper will show that the all-speed schemes’ accuracy will be deteriorated if they are combined with an entropy fix. To avoid the appearances of the shock anomaly, the expansive shock, and the artificial entropy fix, Kim introduced the RoeM-type schemes through a linear perturbation analysis on the odd–even decoupling [24]. However, they still encounter the shock anomaly in some cases [19]. Also, they are incapable of obtaining satisfactory results at low speeds. In this paper, we will propose a new scheme called RoeMAS. It overcomes the Roe schemes’ defects at high speeds by modifying the maximum and minimum wave speeds [25] and improving the density perturbation’s damping rates based on the KIM’s analysis [23]. Also, it can avoid the non-physical behavior, the global cut-off and the checkerboard problem in low speeds’ simulations. All in all, Numerical tests below will show that the RoeMAS proposed can be with a high level of accuracy, robustness, and applicability. This paper is organized as follows. In the second section, we will briefly overview the governing equations. Moreover, we will introduce the RoeMAS’s construction in the 3rd section. To make a thorough analysis, we conduct several numerical tests in Section 4. The last section contains concluding remarks.

12

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

2

10

Roe RoeM2 RoeMAS

Roe RoeM2 RoeMAS

8

6

Density

Pressure

1.5

1

4

0.5 2

0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

X

0.8

1

X

Fig. 1. Pressure and density distributions of a slowly moving contact discontinuity.

2. Governing equations Euler equations are usually used to describe the motion of invscid flows. It can be written in the following form:

∂q + f c (q) = 0 ∂t

(1)

where q is the conservative variable vector, s is the auxiliary variable vector, and f c (q) is the inviscid flux vector. For two-dimensional cases, we have

q = (ρ , ρ u, ρv, ρ e)T

⎛

ρu

(2)

⎞

⎛

(3)

where ρ is density, u and v are velocity vector components, p is static pressure, e is internal energy, and T is temperature. 3. RoeMAS’s construction

1 1 ( f c,i + f c,i+1 ) − R||R−1 (qi+1 − qi ) 2 2

(4)

where R is the right eigenvector matrix of ∂∂fqc with the following form:

⎢ ⎢ ⎣ 1 2



0

1

u˜

−ny

u˜ − nx c˜

v˜

nx

v˜ − ny c˜

⎥ v˜ + ny c˜⎦

(nx v˜ − ny u˜)

H˜ − c˜U˜

H˜ + c˜U˜



1

⎤

1

u˜2 + v˜ 2

(6)

λ1 = U˜ λ2 = U˜ + c˜ λ3

(7)

= U˜ − c˜

1 1 f + δ pRoe d) ( f c,i + f c,i+1 ) − (|λ1 |q + δURoe  2 2

(8)

T

where d is (0, nx , ny , U˜ ) . Specific expressions of δU and δ p for the original Roe scheme can be given as follows:

δURoe = DPRoe

 p (|λ2 | − |λ3 |) U + 2 c˜ ρ˜ c˜2

(|λ2 | − |λ3 |) 2

(9)

ρU  p + DRoe ρ˜ U

(10)

where

The “classical” Roe scheme can be expressed in the following form [1]:

R=⎢

λ2

with

δ pRoe =

3.1. The original Roe scheme

⎡

⎥ ⎥ ⎥ ⎦

λ1 λ3

f Roe c,i+1/2 =

⎟→ ⎜ ρ u2 ⎟ → ⎜ ⎟ ⎜ ⎟ ⎜ f c (q) = ⎜ ⎟ j ⎟ i +⎜ ⎝ ρ uv ⎠ ⎝ ρv2 ⎠ (ρ e + p)u (ρ e + p)v

f Roe c,i+1/2 =

⎤

Through mathematical transformation, Eq. (4) can be written in the following form [26]:

⎞

ρv ρ uv

⎡ λ1 ⎢ ⎢ =⎢ ⎣

u˜ + nx c˜⎥ ⎥

DPRoe =

(|λ2 | + |λ3 |) 2

ρU

− |λ1 | = DRoe

(11)

In addition, the original Roe scheme’s mass flux can be written in the following form if we combine Eq. (4) with Eq. (8):



f ρRoe ,i+1/2 =

1 1 p ρ ( f ρ ,i + f ρ ,i+1 ) − DRoe ρ + DURoe U + DPRoe 2 2 2 ρ˜ c˜

 (12)

ρ

where f ρ is the mass flux and DRoe = |λ1 |. (5)

Here, c is the sound speed, “∼” means Roe average, nx , ny are the unit ˜ x + v˜ ny . normal vectors, and U˜ = un  is a diagonal matrix consisting of the relevant eigenvalues:

3.2. The RoeMAS scheme’s construction 3.2.1. Cure for the defects at low speeds By asymptotic analysis, Guillard found that the original Roe scheme suffers from unphysical behavior problems at low speeds [27,28]. In the Ref. [27], he got the following conjecture for the one-dimensional subgrid model: in the low Mach limit the Roe

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

1.2

1.2

1

1

Roe RoeM2 RoeMAS

Roe RoeM2 RoeMAS

0.8

Density

0.8

Pressure

13

0.6

0.6

0.4

0.4

0.2

0.2

0

0 -5

-4

-3

-2

-1

0

1

2

3

4

5

-5

-4

-3

-2

-1

X

0

1

2

3

4

5

X

Fig. 2. Pressure and density distributions of the Sod shock tube test (t = 1.8 s).

scheme supports pressure fluctuations of the order of the Mach number, while the physical pressures scale with the square of the Mach number. Preconditioning methods slow down the acoustic speed towards the fluid speed and show the correct asymptotic. However, they need to define at least one problem-dependent coefficient. Different coefficient’s definitions may make the solutions different. To avoid using the problem-dependent problem and show the correct asymtotic, Li reached the following conclusions for low speeds through theoretical analysis: (a) DP in the form of Eq. (9) should be of the order between O(c0 ) and O(c−1 ) to suppress the checkerboard problem. (b) DρU in the form of Eq. (10) should be of the order O(c0 ) or lower to be free from the unphysical problem. (c) Cancelling global cut-off in the numerator can overcome the global cut-off problem. To satisfy these requirements, we can keep the term δU the same as the original Roe scheme’s in Eq. (9) while the term δ p is modified as follows:

δp =

(|λ2 | − |λ3 |) 2

ρU  p + DRoeMAS ρ˜ U

(13)

where ρU

DRoeMAS =

(|U˜ + c | + |U˜ − c |) 2

(14)

c = f (M) · c˜

(15)

f (M) = min(M, 1)

(16)

Through mathematical transformation, modifications above can also be written in the form of Eq. (8) as follows:

f c,i+1/2 =

1



1 f + δ pRoe · d ( f c,i + f c,i+1 ) − DρRoeMAS · q + δURoe ·  2 2 1 + [1 − f (M)] f ∗∗ (17) 2

where ρ ˜ |t DRoeMAS = c˜|M

⎛

0

(18)

⎞

⎜ρ cn U ⎟ ⎜ x ⎟ ⎟ ⎝ρ cny U ⎠ ρ cU U

f ∗∗ = ⎜

(19)

3.2.2. Cure for the shock instability Quirk conjectured that any scheme that does not survive the odd– even decoupling is incapable to be free from the shock anomaly’ appearance. Through analyses on the odd–even decoupling, Kim argued that the damping rate of the density perturbation and the feeding rate of the pressure perturbation make the shock unstable simultaneously. By balancing the contribution of the pressure perturbation to the mass flux with that of the density perturbation, the odd–even decoupling could be cured. Based on such idea, Kim introduced the RoeM-type schemes, such as RoeM and RoeM2, by multiplying the coefficients DP and Dρ with functions f and g which are not larger than 1. However, the RoeM-type schemes’ terms DP in the mass flux are of the order O(c−2 ) or smaller and they do not satisfy the requirements for low speeds in Section 3.2.1. Therefore, we will repeat Kim’s analysis in brief and check whether it is necessary to follow the RoeM-type schemes’ ideas to make DP of the order O(c−2 ) or smaller. More details can be found in the Ref. [23]. Assuming that the computational mesh is uniform and the solution at time t n is as follows: If i is even, then

ρin = ρ + ρˆin ,

pni = ρ + pˆni ,

uni = u0 ,

vni = v0

(20)

uni = u0 ,

vni = v0

(21)

else if i is odd, then

ρin = ρ − ρˆin ,

pni = ρ − pˆni ,

where ρˆ , pˆ are the amplitudes of density and pressure perturbations. Thus, we can get the following conclusions:

14

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

3

3

2.9

Pressure

Pressure

2.5

2

2.8

2.7 Roe Roe (with Muller's E-fix) RoeM2 RoeMAS

1.5

Roe Roe (with Muller's E-fix) RoeM2 RoeMAS

2.6

1

2.5 -5

-4

-3

-2

-1

0

1

2

3

4

5

-0.7

-0.6

-0.5

X

-0.4

-0.3

-0.2

-0.3

-0.2

X 3

3

2.9 2.5

Density

Density

2.8 2

2.7 Roe Roe (with Muller's E-fix) RoeM2 RoeMAS

1.5

Roe Roe (with Muller's E-fix) RoeM2 RoeMAS

2.6

1

2.5 -5

-4

-3

-2

-1

0

1

2

3

4

5

-0.7

-0.6

-0.5

X

-0.4

X

Fig. 3. Pressure and density distributions of the Sod shock tube II test (t = 1.5 s).

  2v DP Dρ ρˆin+1 = 1 − 2vy ρˆin − 2y pˆni c˜

c˜

pˆn+1 = (1 − 2vy ) pˆni + 2vy i

(22)

c˜

  γ − 1 u˜2 + v˜ 2 DP + Dρ c˜2

2

c˜

 −1

 − v˜ 2 pˆni (23)

where

vy =

ct y

ρˆin+1 = (1 − 2vy M˜ )ρˆin −

Combing Eqs. (22), (23) with Eq. (11), we can express the original Roe scheme’s amplitude as follows:

ρ

ˆin+1

transverse speed is so small and the pressure perturbation is so large behind a strong shock that the original Roe scheme is incapable of damping the perturbation and inclined to blow up the solution [24]. As mentioned above, the RoeM-type schemes were proposed to control the pressure perturbation near shocks and achieved a considerable success. For example, RoeM can be expressed in the form of Eqs. (24) and (25) as follows:

˜ |)ρ = (1 − 2vy |M

ˆin

2vy ˜ |) pˆn − 2 (1 − |M i c˜

˜ 2 pˆn = (1 − 2vy ) pˆni − 2vy (γ − 1)M pˆn+1 i i

(24)

(25)

Eqs. (24) and (25) show that the density field is perturbed by the pressure perturbation, and the transverse velocity is capable of damping the density perturbation more or less. However, the





2vy ˜ |) pˆn f (1 − |M i cˆ2

˜2 = 1 − 2vy f 1 + (γ − 1)M pˆn+1 i



pˆni

(26)

(27)

However, analyses aforementioned show that the RoeM-type schemes would show the checkerboard problem at low speeds due to their reducement of pressure perturbation. From the right-hand of Eq. (24), we can see the factors that corrupt the density field are comprised of two aspects: the density perturbation and the pressure perturbation at time n. Thus, we can choose to control the density contribution, in place of the

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

15

100

10

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS

-2

Residual

10-4

10-6

10

-8

10

-10

10

-12

0

20000

pressure contribution, to avoid the unphysical phenomenon at low speeds. Because the density damping’s mechanism is different from the pressure damping’s, we modify the Dρ to the following form: ρ

˜ |t Di+1/2, j = c˜|M

t=

(28)

f1

|M˜ |

100000

Fig. 5. Error histories of inviscid flow over NACA0012 airfoil (M∞ = 0.01, 1st order reconstruction and primary grid).

⎡ |λ1 | ⎢ ⎢ || = ⎢ ⎣

⎤ ⎥ ⎥ ⎥ ⎦

|λ1 | |λ2 |

(34)

where

|M˜ | < 1

(30)

u˜2 + v˜ 2 = 0

h1 = min(P1, P2, P3, P4, P5)

P1 = Pi+1/2, j ,

P2 = Pi, j−1/2 ,

P4 = Pi+1, j−1/2 ,



(31)

P3 = Pi, j+1/2 ,

P4 = Pi+1, j+1/2

pi, j pi+1, j , pi+1, j pi, j

|λ1 | = |U˜ | |λ2 | = max(|U˜ + c˜|, |Ui+1 + c˜|)

(35)

|λ3 | = min(|U˜ − c˜|, |Ui − c˜|)

with

Pi+1/2, j = min

80000

(29)

where the function f1 is different from the RoeM-type schemes’ f :

f1 =

60000

|λ3 |

|M˜ |

 ˜ h1 |M|

40000

Iteration Number

Fig. 4. The primary computational grid of NACA0012 (120 × 50).

It should be noticed that the term U is of the order O(c−1 ) or smaller while the term cis of the order O(c) for the low Mach number limit. Thus, the term U’s modification can be disregarded at low speeds, ie: the eigenvalues’ modifications have little influence on the scheme’s performances at low speeds.

(32)

 (33)

Also, in order to capture the contact discontinuity accurately, we make the function t equal to 1 if the function P1 is equal to 1. 3.2.3. Expansion shock, instability in the expansion region It is known that the original Roe scheme is incapable of distinguishing expansion shocks from compression ones. To overcome this defect, an entropy fix is widely used. However, it needs a problemdependent coefficient. Improper coefficient’s definition may deteriorate the resolution. Moreover, if combined with the all-speed schemes, it will deteriorate the schemes’ accuracy at low speeds. Einfeldt showed that underestimating the numerical signal velocities cause the Roe-type scheme’s failure [25]. As a remedy, he computed the numerical signal velocity by incorporating the values at neighboring cells. To avoid the entropy fix’s defects above, we borrow the Einfeldt’s idea and modify the non-linear eigenvalue of Eq. (6) to the following form:

3.2.4. Total enthalpy conservation Another defect of the Roe-type schemes is that they are incapable of preserving total enthalpy in steady inviscid flow. Jameson suggested that this defect can be remedied by introducing a modified state vector which contains the product of density and enthalpy [29]. Following his idea, the q in Eq. (15) can be modified T to (ρ , ρ u, ρv, ρ H ) to satisfy Hanel’s conclusion [30].

3.2.5. The complete form of the RoeMAS scheme Through analysis and constructions above, we obtain the new scheme called RoeMAS. In theory, it is robust against the shock anomaly and satisfies the following property at low speeds: DP is of the order O(c−1 ) to satisfy the property (a) in Section 3.2.1. (2) DρU is of the order O(c0 ) to satisfy the property (b) in Section 3.2.1. (3) Independent of any tunable parameter to satisfy the property (c) in Section 3.2.1. (1)

The RoeMAS scheme for all speeds can be written as follows:

16

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

Roe (with Muller's E-fix)

Roe

F-Roe (with Muller's E-fix)

F-Roe

RoeM2

RoeMAS

Fig. 6. Schemes’ pressure contours of inviscid flows around NACA0012 airfoil with the 1st order reconstruction (Ma = 0.01, primary grid).

RoeMAS f c,i+1/2 =

1 1 f + δp · d ( f c,i + f c,i+1 ) − DρRoeMAS · q∗ + δU ·  2 2

(36)

and

δp =

where T

d = (0, nx , ny , U˜ )

(37)

q = (ρ , ρ u, ρv, ρ H )

δU =

T

DPRoeMAS

 p (|λ2 | − |λ3 | ) U + 2 c˜ ρ˜ c˜2

(38)

(39)

2

ρU  p + DRoeMAS ρ˜ U

(41)

with ρU

∗

(|λ2 | − |λ3 | )

DRoeMAS =

(|λ2 | + |λ3 | ) 2

 |λ2 | = max |U˜ + c |, |Ui+1 + c |  |λ3 | = min |U˜ − c |, |Ui − c |

(42)

(43)

with

DPRoeMAS =

(|λ2 | + |λ3 | ) 2

− |U˜ |

(40)

c = f (M) · c˜

(44)

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

0

17

0

-1

-1

X X

-2

X

log10lnd(P)

log10lnd(P)

-2

-3

X -4

X

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS Theretical Prediction

-5

-6

-7 -3.5

-3

-2.5

-2

-1.5

-1

-0.5

-3

X -4

X

X -5

-6

-7 -3.5

0

-3

-2.5

-2

-1.5

log10M

log10M

(a)

(b) 0

X

-1

-0.5

0

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS X Theretical Prediction

-1

-2

log10lnd(P)

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS Theretical Prediction

-3

X -4

X -5

-6

-7 -3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

log10M

(c) Fig. 7. Pressure fluctuations VS Inflow Mach number (primary grid, (a) 1st order reconstruction, (b) 2nd order reconstruction, (c) 3rd order reconstruction).

f (M) = min(M, 1)

(45)

In addition,



ρ

DRoeMAS =

˜| c˜|M ˜| c˜|M

if P1 = 1 |M˜ | f1

(46)

else if P1 = 1

where the functions f1 and P1 are expressed in Eqs. (30)–(33). 4. Numerical experiments 4.1. 1d Moving contact discontinuity The contact wave is a discontinuous wave across which both pressure and particle velocity are constant but density jumps discontinuously. So as do variables which depend on density,

such as specific internal energy, temperature, sound speed, and entropy. That’s to say, where there is a contact discontinuity, there is a large temperature’s gradient [31]. In supersonic and hypersonic flows, contact discontinuities are so common for the appearance of strong viscosity, shock/shock interaction, and shock/boundary layer interaction. Thus, simulating it with high resolution is a must for a scheme to depict complex flow fields and to predict aerodynamic performance parameters such as the aerodynamic heating. The initial conditions are (ρ , u, p)L = (0.125, 0.1125, 1), (ρ , u, p)R = (10, 0.1125, 1) with a CFL number of 0.7. The iteration count is 500 and the number of grid points is 100. To conduct a fair comparison, all the schemes are combined with the 1st order reconstruction which is constant by cell. For time integration, the explicit Euler scheme is used [32].

18

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

X

-1

log10lnd(P)

-2

-3

-4

0

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS X Theretical Prediction

-2

X X

X -4

X -5

-6

-6

-3

-2.5

-2

-1.5

-1

-0.5

0

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS X Theretical Prediction

-3

-5

-7 -3.5

X

-1

log10lnd(P)

0

-7 -3.5

-3

-2.5

-2

-1.5

log10M

log10M

(a)

(b)

-1

-0.5

0

Fig. 8. Pressure fluctuations VS Inflow Mach number (1st order reconstruction, (a) denser grid, (b) densest grid).

As can be seen in Fig. 1, RoeMAS is capable of capturing contact discontinuity as accurately as Roe and RoeM2. 4.2. Sod shock tube This test is conducted with 400 cells to show RoeMAS’s capability of capturing shocks. The initial conditions are (ρ , u, p)L = (1, 0, 1), (ρ , u, p)R = (0.125, 0, 0.1) with a CFL number of 0.5. The explicit 4th Runge–Kutta time-marching scheme is used. To conduct a fair comparison, all the schemes are combined with the 1st order reconstruction. Fig. 2 shows schemes’ pressure and density distributions at 1.8 s. In it we can see that RoeMAS is capable of producing similar steep gradients near the shock as Roe and RoeM2. 4.3. Sod shock tube II This test is given special initial conditions as (ρ , u, p)L = (3, 0.9, 3), (ρ , u, p)R = (1, 0.9, 1), which is different from the Section 4.2 and a sonic point exits along rarefaction waves. The domain is of length 10 and is discretized using 400 cells. The CFL number is chosen as 0.5 and the explicit 4th Runge–Kutta timemarching scheme is used. It is reminded that in this case only the 1st reconstruction is combined with to compare the schemes fairly. The schemes’ density and pressure distributions at 1.5 s are shown in Fig. 3. The original Roe scheme yields an unphysical expansion shock. Both RoeMAS and RoeM2 prevent the formation of an expansion shock. Also, the entropy fix can make the Roe scheme be free from the entropy violating solution. However, it adds much dissipation. Moreover, RoeMAS is capable of giving much sharper profiles at the head and tail of the rarefaction than the others. 4.4. 2d Inviscid NCA0012 airfoil In this case, we choose an inviscid flow around NACA0012 airfoil to check the RoeMAS’s accuracy at low speeds. The flow conditions are as follows: Ma = 0.001–0.1, no angle of attack.

The grid is C-type. To check the schemes’ grid dependencies, we conduct the computations within three different grids. These grids are 120 (circumferential) × 50 (wall normal) called primary grid (Fig. 4), 240 (circumferential) × 50 (wall normal) called denser grid, and 240 (circumferential) × 100 (wall normal) called densest grid respectively. In addition, we also use different reconstruction technologies for primitive variables to analyze these schemes’ disparities. These reconstruction technologies are as follows: 1st order reconstruction, 2nd MUSCL reconstruction with van Albada limiter [33], 3rd reconstruction with van Leer limiter [34]. For time integration, LU-SGS is used with CFL number equal to 5, and the computations are conducted for 100,000 steps to achieved five orders reduction of the density residual at least. Fig. 5 shows their error histories when they are combined with the 1st reconstruction. In it we can see that all the computations are convergent. RoeMAS and F-Roe converge much faster than the other schemes. Another noticeable phenomenon is that the entropy fix makes F-Roe converge much slower. Fig. 6 shows the schemes’ pressure contours with the same levels when the Mach number is equal to 0.01. It indicates that only RoeMAS and F-Roe are capable of obtaining accurate flowfields as that in the Ref. [35] at low speeds. Also, the entropy fix deteriorates the F-Roe scheme’s accuracy at low speeds. Fig. 7a shows the schemes’ pressure fluctuations Ind(p) = (Pmax − Pmin )/Pmax VS the Mach number when they are combined with the 1st order reconstruction. According to the theoretical asymptotic predictions, the pressure fluctuations scale with the square of the Mach number in the continuous case. Thus, the followings can be drawn: (a) Only RoeMAS and F-Roe are capable of agreeing with the theoretical predictions. (b) The entropy fix deteriorates the F-Roe’s high resolution at low speeds. To make a thorough analysis, Fig. 7b and c display the schemes’ results when they are combined with the 2nd order MUSCL reconstruction with thevan Albada limiter and the 3rd order

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

2

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1

-1.5

-1.5

-2

-2

-2.5

-2.5

-3 -3.5

-3

-2.5

-2

-1.5

-1

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS

1.5

log10Cd

log10Cd

2

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS

1.5

19

-3 -3.5

-0.5

-3

-2.5

-2

log10M

log10M

(a)

(b) 2

-1.5

-1

-0.5

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS

1.5 1

log10Cd

0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5

-3

-2.5

-2

-1.5

-1

-0.5

log10M

(c) Fig. 9. Drag coefficients VS Inflow Mach number (primary grid, (a) 1st order reconstruction, (b) 2nd order reconstruction, (c) 3rd order reconstruction).

reconstruction with the van Leer limiter. Comparing them with Fig. 7a, we can draw conclusions as follows: (a) RoeMAS and F-Roe are nearly independent of the reconstruction method’s choice. (b) Only RoeMAS and F-Roe are capable of agreeing with the theoretical predictions no matter which reconstruction method they are combined with. (c) RoeMAS and F-Roe’s results with the 1st reconstruction are even much better than the other schemes’ with higher order reconstruction methods. Fig. 8 indicates that RoeMAS and F-Roe perform much better than the other schemes no matter which computational mesh they are computed in.

Theoretically, the drag in the flow of this case should be zero. Thus, the computed drag is also a good measure of numerical error [11,36]. Figs. 9 and 10 display the schemes’ drag coefficients when they are combined with different reconstruction methods or different computational meshes. In them we can see that only RoeMAS and F-Roe are with the property of Mach uniformity. Correspondingly, the other schemes’ dissipations increase obviously with the Mach number’s decreasing, even in cases where they are combined with higher order reconstruction methods or computed in denser meshes. Compared with the classic compilation of airfoil data by Abbott and von Doenhoff, we can find that the inviscid drag coefficients predicted by RoeMAS and F-Roe are roughly equal to the measured drag coefficients of real airfoils in viscous flow. This means that if

20

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

2 1.5 1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1

-1.5

-1.5

-2

-2

-2.5

-2.5

-3 -3.5

-3

-2.5

-2

-1.5

-1

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS

1.5

log10Cd

log10Cd

2

Roe Roe (with Muller's E-fix) F-Roe F-Roe (with Muller's E-fix) RoeM2 RoeMAS

-3 -3.5

-0.5

-3

-2.5

-2

log10M

log10M

(a)

(b)

-1.5

-1

-0.5

Fig. 10. Drag coefficients VS Inflow Mach number (1st order reconstruction, (a) denser grid, (b) densest grid).

4.5. Quirk’s test (Odd–Even decoupling)

3.5 3

This test case depicts a planar moving shock in a duct where the centerline grids are perturbed. It is considered that any scheme that does not survive this test case is incapable to be free from the shock anomaly’ appearance. The initial conditions are given as (ρ , u, v, p)L = (1, 6, 0, 1), (ρ , u, v, p)R = (5.25, 0.353, 0, 40.64) to produce a normal shock propagate to the left with a Mach number of 6. The computational domain is [0:1, 0:10] and the shock’s initial location is at x = 8 (Fig. 11). Moreover, the centerline grids are perturbed in the following manner:

2.5

Y

2 1.5 1 0.5 0 -0.5 -2

0

2

4

6

8

10



12

X

Yi, j,mid =

Fig. 11. The initial density contour and axis lines in the primary grid.

Y j,mid + 10−4 −4

Y j,mid − 10

for i even for i odd

For time integration, the explicit 4th Runge–Kutta time-marching scheme is used. At the very start, we just combine these schemes with the 1st order reconstruction to analyze them fairly. In Quirk’s report, a computational mesh with 20 (Y) × 800 (X) was used to examine the schemes’ robustness against the shock anomaly. Also, Kim used this mesh to validate the RoeM-type

these methods were to be employed, then when the viscous drag was added, the drag would be about twice its actual value. Thus, the RoeMAS scheme and the F-Roe scheme are still incapable of obtaining acceptable drag coefficients and a further research on this issue is necessary.

Fig. 12. Aspect ratio of grids.

Level R :

1 2

2 2 .5

3 3

4 3 .5

5 4

6 4 .5

7 5

(47)

8 5 .5

9 6

10 6 .5

11 7

12 7 .5

13 8

Fig. 13. The level of density contours shown in the Quirk’s test case.

14 8 .5

15 9

16 9 .5

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

Roe, delta_x=1

Roe with Muller's E-fix, delta_x=5

F-Roe, delta_x=1

F-Roe with Muller's E-fix, delta_x=5

RoeM2, delta_x=5

RoeMAS, delta_x=5

Fig. 14. Density Contours of the odd–even grid perturbation problem in the primary grid (Ms = 6.0).

Roe with Muller's E-fix, delta_x=5

F-Roe with Muller's E-fix, delta_x=5

RoeM2, delta_x=5

RoeMAS, delta_x=5

Fig. 15. Density Contours of the odd–even grid perturbation problem in the primary2 grid (Ms = 6.0).

21

22

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

Roe with Muller's E-fix, delta_x=1

F-Roe with Muller's E-fix, delta_x=1

RoeM2, delta_x=5

RoeMAS, delta_x=5

Fig. 16. Density Contours of the odd–even grid perturbation problem in the primary3 grid (Ms = 6.0).

2

2

1.8

1.8 RoeM2 RoeMAS

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

RoeM2 RoeMAS

1.6

Pressure

Density

1.6

0 0

0.2

0.4

0.6

0.8

1

Y

0

0.2

0.4

0.6

0.8

1

Y

Fig. 17. Transverse slices of the pressure and density fields for the dotted line in front of the shock shown in Fig. 16.

schemes’ robustness against the shock anomaly. However, a series of researches show that the aspect ratio of cells near the shock is a major factor that influenced the schemes’ performances against the shock anomaly [19,37,38] (Fig. 12). Thus, we use the following three different computational meshes to examine the schemes’ robustness against the shock anomaly: 20 (Y) × 800 (X) with an aspect ratio of 4 called primary grid, 20 (Y) × 200 (X) with an aspect ratio of 1 called primary2 grid, and 80 (Y) × 200 (X) with an aspect ratio of 0.25 called primary3 grid. It is reminded that all the schemes’ resultant density contours are plotted with the same level defined in Fig. 13. Fig. 14 shows that both F-Roe and Roe are more inclined to show the shock anomaly than the others. Moreover, the entropy fix can help them be more robust against the unphysical solution, even in cases where the shock moves five times farther. Also, both RoeMAS and RoeM2 are free from the shock anomaly.

Modifying the aspect ratio from 4 to 1, the entropy fix can still make Roe and F-Roe more robust against the shock anomaly (Fig. 15). However, they all show the unphysical oscillation when the shock moves to x = 5. Correspondingly, RoeMAS and the RoeM2 are still capable of obtaining right solutions. Fig. 16 displays that with a small aspect ratio, even RoeM2 which is robust enough in the other two meshes, shows asymmetry behind the shock (Fig. 17) while RoeMAS does not (Fig. 18). This shows that RoeMAS is even more robust against the shock anomaly than RoeM2. Fig. 19 shows RoeMAS and RoeM2’s density contours when they are computed in the primary3 mesh and combined with the 2nd order MUSCL reconstruction. Moreover, the minmod limiter, which is the most dissipative one among all the limiters, is used to monitor the local gradient of a solution and control spatial order. In this figure we can see that the higher order spatial accuracy can be

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

6

50

RoeM2 RoeMAS

5

RoeM2 RoeMAS

45

Pressure

5.5

Density

23

4.5

40

35

4

30 0

0.2

0.4

0.6

0.8

1

0

0.2

Y

0.4

0.6

0.8

1

Y

Fig. 18. Transverse slices of the pressure and density fields for the dotted line behind the shock shown in Fig. 16.

Roe, delta_x=1

Roe with Muller's E-fix, delta_x=1

F-Roe, delta_x=1

F-Roe with Muller's E-fix, delta_x=1

RoeM2, delta_x=5

RoeMAS, delta_x=5

Fig. 19. Density contours of the odd–even grid perturbation problem in the primary3 grid (Ms = 6.0, 2nd order MUSCL reconstruction with the minmod limiter).

24

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

2

2

1.8

1.8 RoeM2 RoeMAS

RoeM2 RoeMAS

1.6

1.4

1.4

1.2

1.2

Pressure

Density

1.6

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

Y

0.6

0.8

1

Y

Fig. 20. Transverse slices of the pressure and density fields for the dotted line in front of the shock shown in Fig. 19.

6

50

RoeM2 RoeMAS

45

Pressure

Density

5.5

RoeM2 RoeMAS

5

4.5

40

35

4

30 0

0.2

0.4

0.6

0.8

1

0

0.2

Y

0.4

0.6

0.8

1

Y

Fig. 21. Transverse slices of the pressure and density fields for the dotted line behind the shock shown in Fig. 19.

helpful to control the unphysical solution’s appearance more or less. However, RoeMAS is still the only one that can be free from the odd– even decoupling (Figs. 19–21). 5. Conclusion remarks In this paper, we present a new Roe-type scheme called RoeMAS. Theoretical analyses in Section 3 verify that the RoeMAS scheme satisfies the following two merits independent of any tuning coefficient: (a) in low speed cases, it approaches F-Roe’s high resolution with the Mach number approaching 0, (b) in high speed cases, it is robust against the odd–even decoupling and free from the expansion shock’s appearance. Three one dimensional test cases (Sections 4.1–4.4) show that RoeMAS can capture contact discontinuities and shocks as sharply as Roe and RoeM2. Also, it is capable to be free from the expansion shock’s appearance. Moreover, it is capable of giving much sharper profiles at the head and tail of the rarefaction in supersonic/hypersonic flows’ simulations. In Section 4.4,

RoeMAS displays its merits of high resolution and good convergence property at low speeds as F-Roe. In addition, both of their 1st order results are much better than the other schemes’ when they are combined with higher order reconstruction methods at low speeds. Quirk’s test shows that RoeMAS is even with a much higher level of robustness against the shock anomaly than RoeM2. All in all, RoeMAS is attractive because it is robust, accurate, and efficient at all speeds without depending on any problem-dependent parameter. In addition, it can easily be coded and extended to practical simulations of complex flows. Acknowledgements This study was co-supported by National Basic Research Program of China (No. 2009CB724104). We are grateful to Shuhai Zhang, China Aerodynamics Research and Development Center, for his guidance. The first author thanks Ju Wang, University of Michigan, for much helpful advice.

F. Qu et al. / Computers and Fluids 121 (2015) 11–25

References [1] Roe PL. Approximate Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357–72. [2] Qu F, Yan C, Yu J, et al. A new flux splitting scheme for the Euler equations. Comput Fluids. doi: http://dx.doi.org/10.1016/j.compfluid.2014.07.004. [3] Liou MS. A further development of the AUSM_ scheme towards robust and accurate solutions for all speeds. AIAA Paper 2003-4116; 2003. [4] Jack RE, Liou MS. Low-diffusion flux-splitting methods for flows at all speeds. AIAA Paper 1997-1862; 1997. [5] Li XS, Gu CW. An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour. J Comput Phys 2008;227:5144–59. [6] Weiss JM, Smith WA. Preconditioning applied to variable and constant density flows. AIAA J 1995;33:2050–7. [7] Turkel E. Preconditioning technique in computational fluid dynamics. Annu Rev Fluid Mech 1999;31:385–416. [8] Unrau D, Zingg DW. Viscous airfoil computations using local preconditioning. AIAA Paper 1997-2027; 1997. [9] Li XS, Gu CW. Preconditioning method and engineering application of large eddy simulation. Sci China Ser G: Phys Mech Astron 2008. doi:10.1007/s11433-0080054-1. [10] Li XS, Gu CW. Development of Roe-type scheme for all-speed flows based on preconditioning method. Comput Fluids 2009;38:810–17. [11] Kitamura K, Shima E. Parameter-free simple low-dissipation AUSM-family scheme for all speeds. AIAA J 2011;49:1693–709. [12] Li XS, Gu CW. The momentum interpolation method based on the time-marching algorithm for all-speed flows. J Comput Phys 2010;229:7806–18. [13] Thornber B, Drikakis D. Numerical dissipation of upwind schemes in low Mach flow. Int J Numer Meth Fluids 2008;56:1535–41. [14] Rieper F. A low-Mach number fix for Roe’s approximate Riemann solver. J Comput Phys 2011;230:5263–87. [15] Fillion P, Chanoine A, Dellacherie S, et al. A new platform for core thermal– hydraulic studies. Nucl Eng Des 2011;241:4348–58. [16] Thornber B, Mosedale A, Drikakis D, et al. An improved reconstruction method for compressible flows with low Mach number features. J Comput Phys 2008;227:4873–94. [17] Dellacherie S. Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number. J Comput Phys 2010;229:978–1016. [18] Li XS, Gu CW. Mechanism of Roe-type schemes for all-speed flows and its application. Comput Fluids 2013;86:56–70. [19] Kitamura K, Shima E, et al. Evaluation of Euler fluxes for hypersonic heating computations. AIAA J 2010;48:763–76.

25

[20] Peery KM, Imlay ST. Blunt-body flow simulations. AIAA Paper 1988-2904; 1988. [21] Kermani MJ, Plett EG. Modified entropy correction formula for the Roe scheme. AIAA Paper 2001-0083; 2001. [22] Muller B. Simple improvements if an upwind TVD scheme for hypersonic flow. AIAA Paper 1989-1977; 1989. [23] Kim SS, Kim C, Rho OH, et al. Cures for the shock instability: development of a shock-stable Roe scheme. J Comput Phys 2003;185:342–74. [24] Quirk JJ. A contribution to the great Riemann solver debate. Int J Numer Meth Fluids 1994;18:555. [25] Einfeldt B, Munz CD, et al. On Godunov-type methods near low densities. J Comput Phys 1991;92:273. [26] Weiss JM, Smith WA. Preconditioning applied to variable and const density flows. AIAA J 1995;33:2050–7. [27] Guillard H, Viozat C. On the behaviour of upwind schemes in the low mach number limit. Comput Fluids 1999;28:63–86. [28] Guillard H, Murrone A. On the behaviour of upwind schemes in the low mach number limit: II. Godunov type schemes. Comput Fluids 2004;33:655–75. [29] Jameson A. Analysis and design of numerical schemes for gas dynamics 2: artificial diffusion and discrete shock structure. Int J Comput Fluid Dyn 1995;5. [30] Hanel D. On the accuracy of upwind schemes in the solutions of the Navier–Stokes equations. AIAA Paper 87-1105; 1987. [31] Toro EF. Riemann solvers and numerical methods for fluid dynamics. 3rd ed. Springer-Verlag; 2009. [32] Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. J Comput Phys 1997;126:202–8. [33] van Leer B. Towards the ultimate conservation difference scheme V: a secondorder sequal to Godunov’s method. J Comput Phys 1979;32:101–36. [34] Kim KH, Kim C. Accurate, efficient and monotonic numerical methods for multidimensional compressible flows Part II: multi-dimensional limiting process. J Comput Phys 2005;208:570–615. [35] Kitamura K, Fujimoto K, et al. Performance of low-dissipation euler fluxes and preconditioned implicit schemes in low speeds. AIAA Paper 2010-1272; 2010. [36] Qu F, Yan C, Yu J, et al. A study of parameter-free shock capturing upwind schemes on low speeds’ issues. Sci China Tech Sci 2014;57:1183–90. [37] Henderson J, Menart JA. Grid study on blunt bodies with the Carbuncle phenomenon. AIAA Paper 97-3904; 1997. [38] Pandolfi Maurizio, D’Ambrosio Domenic. Numerical instabilities in upwind methods: analysis and cures for the “carbuncle” phenomenon. J Comput Phys 2011;166:271–301.