1 Introduction

The Roe scheme^[1] is one of the most famous and important shock-capturing schemes because of its high accuracy. This scheme has undergone considerable development, such as its extension to incompressible flows^[2-4], and has been extensively used for flow computation, such as in Euler flows^[5], large eddy simulation^[6-7], and cavitation^[8]. However, the Roe scheme also suffers from a few shortcomings, such as the shock instability and expansion shock^[9].

The shock instability is a well-known defect of supersonic flows with different performances, such as carbuncle, kinked Mach stem, and odd-even decoupling. Several methods were proposed to cure shock instability, and the cures were achieved by adding an entropy fix^[10], combining a dissipative scheme^[9], increasing the basic upwind dissipation^[11-12], and considering multi-dimensional characteristics^[13-14].

The expansion shock is another defect of the Roe scheme, which is an unphysical solution that violates the entropy condition. Moreover, this defect often yields unacceptable values, such as negative pressure and density, and leads to the divergence of computation for the highly energetic flow. The entropy fix is also often adopted to overcome this drawback. However, this approach has limited effects when introducing the large numerical dissipation and unfavorable empirical parameters. Another considerably common curing method introduces a slight modification by redefining the numerical signal velocities to obtain improved results^{[12, 15-16]}.

In Ref. [17], the momentum interpolation mechanism (MIM) in the Roe scheme^[18-20] is considered as the most important reason for shock instability. Thus, an improved Roe scheme is proposed^[17] by removing the MIM for the non-linear flow. This improvement cures the shock instability when removing the problem-independence empirical parameters and decreasing the numerical dissipation. However, in this paper, several numerical results show an unexpected defect, wherein the expansion shock becomes serious, and the traditional curing method is invalid. To broaden the range of applications of the improved Roe scheme, the current study is aimed at curing the expansion shock by identifying the reason for the deterioration and further elucidating the mechanism that produces the expansion shock, and thus proposes an ideal scheme.

The rest of this paper is organized as follows. Section 2 provides the governing equations and the improved Roe scheme for curing shock instability^[17]. Section 3 analyzes the mechanism of the traditional method of curing the expansion shock. Section 4 provides a new approach to cure the expansion shock when maintaining all the advantages of the improved Roe scheme. Section 5 concludes this paper.

2 Governing equations and Roe scheme 2.1 Governing equations

The governing three-dimensional Navier-Stokes equation can be written as follows:

(1)

where Q = [ρ ρu ρv ρw ρE]^T is the vector of the conservation variables, F = [ρu ρu² + p ρuv ρuw ρuH]^T, G = [ρv ρuv ρv² + p ρvw ρvH]^T, and H = [ρw ρuw ρvw ρw² + p ρwH]^T are the vectors of the Euler fluxes, ρ is the fluid density, p is the pressure, E is the total energy, H is the total enthalpy, and u, v, and w are the velocity components in the Cartesian coordinate (x, y, z), respectively.

2.2 Roe scheme and improvement

The classical Roe scheme can be expressed in the following general form as the sum of a central term and a numerical dissipation term:

(2)

where is the central term, and is the numerical dissipation term. For a cell face of the finite volume method,

(3)

(4)

where n_x, n_y, and n_z are the components of the face-normal vector, and U=n_xu+n_yv+n_zw is the normal velocity on the cell face.

According to Ref. [21], a scale uniform framework for the shock-capturing scheme is proposed^[22]. This framework is simple and easy to analyze and improve with low computation cost,

(5)

where the first term on the right side ξ is the basic upwind dissipation, the term δp_p is the pressure-difference-driven modification for the cell face pressure, δp_u is the velocity-difference-driven modification for the cell face pressure, δU_u is the velocity-difference-driven modification for the cell face velocity, and δU_p is the pressure-difference-driven modification for the cell face velocity.

For the classical Roe scheme,

(6)

(7)

(8)

(9)

(10)

where c is the sound speed, and the eigenvalues of the system are defined as follows:

(11)

(12)

(13)

Based on Eqs. (11)-(13), Eqs. (7)-(10) can also be further simplified as

(14)

(15)

(16)

(17)

Reference [17] suggested that the shock instability is mainly attributed to the term δU_p, which plays a role of the MIM to smooth the shock, and can be cured by removing the MIM for high-Mach number flows and multiplying Eq. (17) to the two functions s₁ and s₂ as follows:

(18)

(19)

where Ma is the Mach number. The function s₂, which is a shock detector and can be obtained in Ref.[17], is not presented in this paper because it is relatively complicated and probably unnecessary for general cases. This improvement is simple and effective.

2.3 Two classical numerical tests

Two classical numerical examples are available to test the expansion shock by the proposed scheme. One is a specific one-dimensional shock tube, and the other is the supersonic corner problem. The initial conditions of the shock tube are given as ρ_L=3, u_L =0.9, p_L =3, ρ_R =1, u_R =0.9, and p_R =1 at the x-axis position of 0.3. The supersonic corner problem considers a moving supersonic shock around a 90° corner. The initial conditions are given as ρ_L =7.041 08, u_L =4.077 94, v_L =0, p_L =30.059 42, ρ_R =1.4, u_R =v_R =0, and p_R =1 at the x-axis position of 0.05. In this study, the grid numbers are 200 for the shock tube and 400×400 for the supersonic corner. For the time discretization, the four-stage Runge-Kutta scheme is adopted. For the space discretization, the first-order accuracy is adopted (unless otherwise specified) to discuss the schemes themselves.

3 Analysis of traditional method curing expansion shock 3.1 Traditional curing method

To avoid the expansion shock, the traditional curing method redefines the physical signal velocities (i.e., non-linear eigenvalues λ₄ and λ₅). For example, Ref.[15] proposed that

(20)

(21)

To precisely obtain a contact discontinuity, Ref.[12] suggested that only U in λ₄ and λ₅ should be improved,

(22)

(23)

Equations (22)-(23) can retain the capability of Eqs. (20)-(21) to suppress the expansion shock when having the advantage of obtaining the contact discontinuity. Therefore, only Eqs. (22)-(23) are discussed in this study.

3.2 Performances of schemes

Figures 1 and 2 show the results by the classical Roe scheme, as described by Eqs. (6) and (14)-(17), and the traditional curing method, as described by Eqs. (6)-(11) and (22)-(23).

For the shock tube, the classical Roe scheme evidently produces an expansion shock at the position x=0.3. The traditional curing method demonstrates substantially the improved performance. However, a sight gap also exists with the traditional method.

For the supersonic corner, a series expansion waves exist around the corner. Thus, the numerical computation could produce an expansion shock. The results of the classical Roe scheme are shown in Fig. 2(a). No evident expansion shock is observed. However, the shock instability is expectedly strong. The traditional curing method produces similar results (see Fig. 2(b)).

By employing the improved Roe scheme in Eqs. (6), (14)-(16), and (18)^[17], the results of the one-dimensional shock tube are similar to those with the classical Roe scheme (see Fig. 3(a)), because the improvement of Eq. (18) only affects the multi-dimensional computation as analyzed in Ref.[17]. For the two-dimensional computation of the supersonic corner, results become significantly different from those of the classical Roe scheme (see Fig. 3(b)). The shock instability becomes substantially weak and is nearly cured. However, a strong expansion shock occurs. In the iterative calculation process, the density occasionally becomes negative, and the following limitation is necessary to prevent the computational divergence:

(24)

where ε is a small positive value.

The above unfavorable phenomena can be explained as follows. According to the discussion in Ref.[17], the role of the MIM is to smooth the shock, whether the shock is physical or unphysical. Therefore, the improvement of removing MIM can prevent physical compression shock from being destroyed. However, the adverse side effect is to keep unphysical expansion shock.

The traditional curing method for expansion shock in Eqs. (22)-(23) can also be integrated into the improved Roe scheme. For the shock tube, this method can produce results as Fig. 1(b). However, this approach is invalid for the supersonic corner and even substantially increases the expansion shock. The computation diverges because of negative density even when using Eq. (24). This unexpected problem seems confusing and may hinder the possible extensive application of the improved Roe scheme because of concerns regarding computational robust. Therefore, in the following sections, the mechanism of preventing expansion shock is further analyzed, and the new method satisfies the stringent requirement of simultaneously curing the expansion shock and the shock instability without additional costs.

3.3 Analysis of schemes

To develop the new method, the mechanism of the traditional curing method is first further analyzed. Equations (22)-(23) can be decomposed into the following five conditions:

(ⅰ) |U| < c

(25)

(26)

(ⅱ) U>c and U_L >c (U_R >c or U_R < c)

(27)

(28)

(ⅲ) U > c, U_L < c, and U_R > c

(29)

(30)

(ⅳ) U < -c and U_R < -c

(31)

(32)

(ⅴ) U < -c, U_R > -c, and U_L < -c

(33)

(34)

By considering

(35)

Eqs. (25)-(34) can be summarized as follows:

(36)

(37)

where

(38)

(39)

Therefore, Eqs. (7)-(10) become

(40)

(41)

(42)

(43)

In the preceding equations, U is the average of U_L and U_R, and the reasonable value of U is between the range of U_L and U_R,

(44)

For the general average method, such as a simple average or a Roe average,

(45)

For compression flows, b_R=b_L=0 because U_R < U < U_L. Therefore, only the expansion flows are considered for analyzing the increment terms b_R -b_L and b_R+b_L in Eqs. (40)-(43).

(ⅰ) |U| < c

(46)

(47)

Therefore, for the subsonic expansion flows, δp_u and δU_p are increased by the modification of Eqs. (22)-(23).

(ⅱ) U > c and U_L > c

(48)

(49)

Consequently, δp_p and δU_u increase. However, these results are not suitable. For supersonic flows, all increment terms should be zero because of the upwind characteristics.

(ⅲ) U > c, U_L < c, and U_R > c

(50)

(51)

Therefore, the increment terms are unsatisfactory because they are not equal to zero and not smooth between the conditions (ⅰ) and (ⅱ).

The other two conditions of U < -c are not discussed for simplicity because they produce the same conclusions as the conditions of U>c.

3.4 Further analysis of curing expansion shock mechanism

The preceding discussion reveals a few unsatisfactory features of the traditional curing method in Eqs. (22)-(23). Moreover, the discussion provides clues regarding the mechanism of expansion shock suppression. Two inspirations are obtained as follows.

(ⅰ) An increment factor is designed as follows:

(52)

where U-U_L in Eq. (38) and U_R -U in Eq. (39) are replaced by to make Δs strictly be equal to zero for the compression and supersonic flows |U|>c. For the subsonic expansion shock, a positive increment which is the same as Eq. (46) is produced.

(ⅱ) Eight cases (see Table 1) are designed to test the effect of changing the values of δp_u, δp_p, δU_u, and δU_p,

(53)

(54)

(55)

(56)

Figure 4 shows the effect of changing the terms. The computations of Cases 2, 4, 5, and 7 diverge, and the results before divergence are provided in Fig. 4. The results reveal that the expansion shock can become considerably serious by decreasing the coefficients of ΔU in δp_u and Δp in δU_p and increasing those in δp_p and δU_u. Otherwise, the expansion shock is suppressed to increase the coefficients in δp_u and δU_p and decrease the coefficients in δp_p and δU_u.

According to the preceding discussion, the mechanism of the traditional curing method for Eqs. (22)-(23) and the improvement of Eq. (18) are substantially understood. Equations (22)-(23) increase the coefficients in δp_u and δU_p (see Eqs. (40) and (43)) and cure the expansion shock. Equation (18) decreases δU_p to zero for high Mach number flows, particularly for multi-dimensional calculations when Ma → 1 but U → 0.. Thereafter, the problem of expansion shock deteriorates (see Fig. 3(b)).

4 Simultaneous improvement of curing expansion shock and shock instability

Although Eq. (18) worsens the expansion shock, this condition is reasonable and necessary to suppress shock instability^[17]. The traditional curing method of Eqs. (22)-(23) only increases the coefficients in δp_u and δU_p and does not completely utilize the potential to decrease the coefficients in δp_p and δU_u. Therefore, an improved Roe scheme is proposed to simultaneously cure expansion shock and shock instability as follows:

(57)

(58)

(59)

(60)

(61)

(62)

Although the analyses in Subsections 3.3 and 3.4 seem complex, the present scheme is significantly concise and easy to implement, and the computational cost only has a negligible increase as well. Compared with the scheme of Eqs. (14)-(16) and (18), only |U| is redefined as |U^'| by Eq. (62), which can also be expressed as follows:

(63)

Thus, the value of |U^'| is decreased for subsonic expansion flows but still within a reasonable range as given in Eq. (44). Therefore, δp_p and δU_u decrease, and δp_u and δU_p increase synchronously, which provide sufficient power to cure the expansion shock even δU_p is decreased by the functions s₁ and s₂ in Eq. (19).

Figures 5 and 6 show the numerical results of the present scheme. Higher-order reconstruction methods^[23-27] are generally used for practical problems. Thus, the monotonic upstream-centered scheme for conservation laws (MUSCL) is also adopted to test the higher-order performance of the improved equations (62) or (63). The computational processes are robust, and all results are satisfactory, particularly for the supersonic corner test, where the expansion shock and the shock instability are simultaneously cured. No adverse side effects are reported for the improvement.

5 Conclusions

The performance of several Roe-type schemes is discussed in terms of expansion shock, and the mechanism of curing expansion shock is analyzed based on the traditional method. Several unfavorable features of the traditional curing method are discovered, and the possible curing mechanism is not completely utilized. Therefore, an improved method is proposed to overcome these problems. The present scheme is substantially with low computational costs, concise, easy to implement, and robust. This scheme is particularly well compatible with the improvement to cure the shock instability by removing the MIM. Therefore, the present scheme is simultaneously free from the problem of shock instability and expansion shock without additional expense.