Healing of the Carbuncle Phenomenon for AUSMDV Scheme on Triangular Grids

1 Introduction

Inviscid numerical flux functions for the compressible Euler equations have been subjected to extensive investigation by many computational fluid dynamics (CFD) researchers during the past several years. This is due to the fact that the flux functions are critical to the robustness and accuracy of numerical schemes. The upwind schemes, in the framework of the Riemann solver by Godunov [1], have proved to be efficient schemes to capture shock wave and contact discontinuities in many compressible inviscid flows. They can be categorized as either flux vector splitting (FVS) or flux difference splitting (FDS). The FVS schemes [2, 3] are known to be fast, simple, and robust to capture strong shock waves and rarefaction waves, especially in hypersonic flows. However, many numerical simulations indicate that these types of schemes are too dissipative and tend to deteriorate the boundary layer profiles. The FDS schemes [4, 5], which employ approximate Riemann solvers, are generally very robust and require no explicit dissipation. But the scheme requires longer computational time compared with the FVS scheme for its matrix calculations. The Roe’s FDS scheme produces unphysical solutions to the system of Euler equations in certain problems such as a strong receding flow (or double rarefaction) that violates the entropy condition and shock anomalies associated with the carbuncle phenomenon on both quadrilateral and triangular grids [6–10]. These two deficiencies are believed to be caused by different mechanisms. The carbuncle phenomenon usually occurs in high-speed aerodynamic simulations, where shock-capturing finite-volume methods are designed to preserve contact discontinuities [6, 11, 12]. In multidimensional simulations, it is difficult or almost impossible to predict what kind of problems are affected by the carbuncle phenomenon. Many researchers have proposed the entropy fix formulations to replace the near-zero eigenvalues by some tolerances in order to alleviate the shock anomalies associated with the carbuncle phenomenon [8, 9, 13, 14]. There are many explanations for the carbuncle phenomenon. Robinet et al. [11] explained that the origin of the carbuncle could lie in the physical instability of the surface of discontinuity (such as shock wave) itself. Pandolfi and D’Ambrosio [14] concluded that the schemes provide insufficient dissipation in the direction parallel to shock wave region, thereby suffering from instabilities that may result in the formation of the carbuncle phenomenon. Recently, Ramalho et al. [10] show that the origin of the carbuncle phenomenon may come from the vorticity generated by the misalignment of pressure gradients across the shock wave with density gradients artificially created within the unphysical numerical shock structure. By introducing unphysical states in the solution, shock-capturing schemes could produce the Richtmyer–Meshkov instability, which is physically and inherently related to the interaction of shock waves with density inhomogeneities.

In 1993, a less dissipative upwind scheme, the so-called advection upstream splitting method (AUSM), was presented by Liou and Steffen [15] as a fast, simple, robust, and accurate method in comparison with the existing numerical schemes. The idea behind the AUSM scheme is to combine the efficiency of the FVS and the accuracy of the FDS schemes together and to give a vanishing numerical diffusivity at stagnation points. The AUSM-family schemes split the flux into a convection and pressure parts that are treated separately. The interface flux functions can be written as a central difference contribution and a diffusive contribution similar to most upwind schemes. Differences among the various approaches are due primarily to behavior of the diffusive contribution to the critical sonic and stagnation points. The accuracy and efficiency of the AUSM scheme are comparable to the Roe’s FDS scheme with less computational resource requirements. In recent years, several attempts have been made to improve the original AUSM scheme [16–22]. Nowadays AUSM-family schemes are often used for hypersonic flow computations due to their robustness and simplicity. The AUSMDV scheme is proposed by Wada and Liou [17] to remove numerical dissipation on contact discontinuities and to conserve enthalpy for steady flows. Unfortunately, it cannot completely remove the carbuncle phenomenon and post-shock oscillations in many numerical tests. Some numerical treatments are needed to stabilize the schemes, especially for higher-order computation. In spite of that, the AUSMDV scheme continues to be used by researchers [23–26] for simulating a wide range of compressible flow problems due to its high accuracy and simplicity.

This paper is aimed to study the numerical instability (carbuncle phenomenon) of the AUSMD, AUSMV, and AUSMDV schemes on two-dimensional structured triangular grids by examining these scheme performances against certain problems suggested by Quirk [6]. Then a multidimensional dissipation technique is proposed to alleviate the numerical instability of these schemes. Solutions are then compared with those obtained from the original schemes. The rest of this paper is organized as follows. In Section 2, a brief description of the governing equations and the numerical flux formulations of the AUSMD, AUSMV and AUSMDV schemes are presented. Some well-known problems that exhibit the numerical shock instability from many shock-capturing schemes are examined. In Section 3, the linearized analysis is applied to the odd–even decoupling problem with the presence of perturbations that may cause the carbuncle phenomenon. The recursive equations are derived in order to study the damping mechanism of these schemes. A more stable numerical flux computation is introduced by implementing the proposed multidimensional dissipation technique to the numerical dissipation term. Finally, the proposed schemes are further extended to achieve a second-order accuracy and assessed against several test cases shown in Section 4.

2 Review of the AUSMDV scheme

Two-dimensional compressible Euler equations can be expressed in a conservation form as

where flow and flux vectors are

A compact vector form of eq. (1) can be written as

where F=Exn⃗x+Eyn⃗y, with n⃗x and n⃗y are unit normal vectors in x- and y-directions, respectively. For a calorically perfect gas, the equation of state is given by

with γ=1.4 for air.

For all AUSM-family schemes, the inviscid flux is explicitly split into convective and pressure fluxes [16] as

The numerical flux at the cell interface of the original AUSMD scheme [17] is calculated from

where

A common mass flux is defined by

with

where V=unx+vny,

The interface pressure is calculated from

with

To improve a shock-capturing capability, the AUSMV scheme is then proposed by modifying the convective flux in the momentum equations of the AUSMD scheme as

Finally, the AUSMDV scheme combines both the AUSMD and AUSMV schemes for the convective flux in the momentum equations

where s is a switching function depending on the pressure gradient given by

and the value of K is fixed at 10 is in this paper.

To investigate the robustness of these schemes, various numerical tests including a Mach 15 flow over a cylinder, a Mach 6 slowly moving normal shock, and a Mach 5 shock reflection over a 46° wedge are used. The imposed boundary conditions are as follows: free-stream value is specified as the inflow condition, an extrapolation from the inner computational domain is used for the outflow condition, and slip condition at a wall is specified for velocity. All the computations are carried out with the first-order scheme in both space and time with an explicit time integration method. The first-order scheme is used because it is meaningful for studying and analyzing the flux function numerically.

2.1 Mach 15 flow over a cylinder

The carbuncle phenomenon refers to a spurious bump on detached shock near the centerline of the flow just ahead a cylinder. The phenomenon is highly grid dependent [14] and can occur even if the number of grid points is not so large. In this case, a Mach 15 flow over a cylinder with a radius of 1.5 is studied. The computational domain consists of a structured triangular grid with 14×318 cells in the radial and circumferential directions, respectively. Figure 1 shows the steady-state density contours obtained from the AUSMD, AUSMV, and AUSMDV schemes. A serious carbuncle phenomenon is produced along the centerline of the cylinder when AUSMD scheme is used, but completely removed by employing the AUSMV and AUSMDV schemes. It should be noted that there is a circulation behind the kinked detached shock wave in the vicinity of the stagnation line that perturbs the solution and may induce some oscillations as shown by streamline contours and velocity vectors in Figure 2. Robinet et al. [11] concluded that the carbuncle phenomenon is a pure numerical mechanism with no connection to the physics.

Figure 1:

A Mach 15 flow over a blunt body: AUSMD (left), AUSMV (middle), and AUSMDV (right).

Figure 2:

AUSMD scheme: Streamline (left) and velocity vectors (right).

2.2 A Mach 6 slowly moving normal shock

The problem of a Mach 6 planar shock wave moving through a two-dimensional duct [6] is chosen to examine the carbuncle phenomenon. The grid at the mid-channel is perturbed alternately at odd and even points with a small magnitude of ±10−6 in order to initiate the instability. The computational domain consists of a uniform triangular grid with 800×20 uniform intervals along x− and y−directions, respectively. For some numerical schemes such as the Roe’s FDS and Harten-Lax-van Leer contact (HLLC) schemes, this grid perturbation is sufficient to trigger the carbuncle phenomenon [6, 8, 9]. Figure 3 shows the density contours of the normal shock wave at the three locations along the duct obtained from the AUSMD, AUSMV and AUSMDV schemes. Carbuncle phenomenon is generated by all schemes but it is less severe than that obtained from the AUSMDV scheme reported in Wada and Liou [17]. This is because the triangular grid has more numerical dissipation than the quadrilateral grid.

Figure 3:

A Mach 6 slowly moving normal shock: AUSMD (top), AUSMV (middle), and AUSMDV (bottom).

2.3 Mach 5 shock reflection over a 46° wedge

The kinked Mach stem generated from a Mach 5 shock wave that moves over a 46° wedge [6] is the third test case for testing the robustness of these schemes. The computational domain consists of a structured triangular grid with 200×200 cells in x- and y-directions, respectively. It is well known that the grid aspect ratio is very significant in the establishment of the carbuncle phenomenon [14] such that very elongated grid sizes along the normal direction to the shock wave promote instability. For this problem, the grid aspect ratio (Δx:Δy) of cells around the right end of the wedge is approximately 5:1, which can potentially triggering carbuncle phenomenon. Figure 4 shows the density contours of double Mach reflection simulation at time t=0.11. The AUSMD scheme clearly exhibits the kinked Mach stem and post-shock oscillation. The carbuncle phenomenon is also present in the solutions by the AUSMV and AUSMDV schemes, but it is much weaker than that produced by the AUSMD scheme.

Figure 4:

A Mach 5 shock reflection over a 46° wedge: AUSMD (top), AUSMV (middle), and AUSMDV (bottom).

3 Numerical analysis

Numerical experiments in the previous section indicated that the AUSMD, AUSMV, and AUSMDV schemes potentially produce numerical instabilities when applied to certain problems. To study the odd–even decoupling numerical instability mechanism, the linearized analysis presented by Pandolfi and D’Ambrosio [14] is applied to the schemes. By assuming uniform base flow conditions with normalized values of ρ0=1, u0≠0, v0=0, and p0=1, a perturbation of (ρˆ,uˆ,pˆ) is then imposed onto the base flow properties, i. e., ρ=1±ρˆ, u=u0±uˆ, and p=1±pˆ at cell (N,M) and (N,M+1). The numerical flux at the cell interface (N,M+1/2) for these three schemes is given in eqs (17)–(19):

It is clear from eqs (17)–(19) that for all schemes, pressure perturbation contributes to every equation except the y-momentum equation and is expected to play an important role to the numerical instability of all these schemes. The odd–even decoupling problems resulting from these numerical fluxes in eqs (17)–(19) are further examined. By starting from the set of initial perturbation (ρˆ0,uˆ0,pˆ0), the result of the first-order scheme at step K+1 is calculated from

By defining the numerical flux at the cell interface as

eq. (20) becomes

where ν is the Courant number.

Perturbations at time step K+1 of the AUSMD, AUSMV, and AUSMDV schemes are given by eqs (23)–(25), respectively:

To obtain a better understanding of damping mechanism of the perturbation from each scheme, the responses to initial values ρˆ0≠0, pˆ0≠0, and both of ρˆ0,pˆ0≠0 are shown in Figures 5–7. All schemes fail to damp density perturbations generated either from an initial density perturbation ρˆ0≠0 or an initial pressure perturbation (pˆ0≠0), i. e., AUSMDV-family schemes are sensitive to density perturbation. However, when both initial density and pressure perturbations are not zero, i. e., ρˆ0≠0 and pˆ0≠0, both density and pressure perturbations are damped.

Figure 5:

Odd–even decoupling problem of the AUSMD scheme: ρˆ0=0.01 (top), pˆ0=0.01 (middle), and ρˆ0=pˆ0=0.01 (bottom).

Figure 6:

Odd–even decoupling problem of the AUSMV scheme: ρˆ0=0.01 (top), pˆ0=0.01 (middle), and ρˆ0=pˆ0=0.01 (bottom).

Figure 7:

Odd–even decoupling problem of the AUSMDV scheme: ρˆ0=0.01 (top), pˆ0=0.01 (middle), and ρˆ0=pˆ0=0.01 (bottom).

From the viewpoint of linearized analysis of the odd–even decoupling problem, eqs (23)–(25) suggest that the AUSMV and AUSMDV schemes are more stable than the AUSMD scheme. According to the previous work by many researches [6, 8, 9, 14, 17, 19, 22], the carbuncle phenomenon seems to be a common problem of multidimensional shock-capturing schemes that require an information exchange between an intermediate point and its neighbor points – also intermediate shock points. The carbuncle appears inside the zone of planar shock wave that aligns with an edge of the grid. And the scheme will suffer from the carbuncle phenomenon if one of these perturbation quantities ρˆMK+1, uˆMK+1, or pˆMK+1 does not suppress to machine zero. To remedy the numerical instability of these schemes, the multidimensional dissipation technique in Phongthanapanich and Dechaumphai [8] is written in a scalar form for the AUSMDV-family schemes. This technique is originally proposed for quadrilateral grid by Sanders et al. [13]. A more efficient version of this technique is presented in Pandolfi and D’Ambrosio [14]. Phongthanapanich and Dechaumphai [8, 9] apply this technique to the triangular grid. This paper proposes a modified version of the multidimensional dissipation technique suitable for applying to |m1/2| term in the AUSMD scheme as

and δ is calculated from

where ηi=|Vi,L−Vi,R|+|ai,L−ai,R| and i=1,…,4 (see Figure 8). A constant κ is usually less than or equal to 1 of the first-order schemes and more than or equal to 1 of the second-order schemes. After tuning κ with some numerical experiments, an appropriate range of κ for the higher-order scheme is found to be in between 1 and 6 (usually between 2 and 5).

Figure 8:

Cell interface of the numerical instability region.

In order to apply eq. (26) to the AUSMV scheme, it is recommended to rewrite eqs (14a) and (14b) in the following form:

where

Then eq. (26) can be applied to the |m1/2| term directly. The convective flux in the momentum equation for the AUSMDV scheme can be determined from eq. (15).

The perturbations at time step K+1 of the modified AUSMD, AUSMV and AUSMDV schemes are derived again as shown in eqs (30)–(32), and the responses to initial non-zero perturbations are also shown in Figures 9–11, respectively. It was found that the dissipation factor (δ) affects all perturbations except density perturbation. However, as pressure disturbance diminishes, the density perturbation damps out. By implementing the multidimensional dissipation technique to these schemes, all perturbations are damped after a few iterations, and recursive equations of perturbations for these schemes are given by

Figure 9:

Odd–even decoupling problem of the AUSMD with δ=0.5: ρˆ0=0.01 (top), pˆ0=0.01 (middle), and ρˆ0=pˆ0=0.01 (bottom).

Figure 10:

Odd–even decoupling problem of the AUSMV with δ=0.5: ρˆ0=0.01 (top), pˆ0=0.01 (middle), and ρˆ0=pˆ0=0.01 (bottom).

Figure 11:

Odd–even decoupling problem of the AUSMDV with δ=0.5: ρˆ0=0.01 (top), pˆ0=0.01 (middle), and ρˆ0=pˆ0=0.01 (bottom).

The modified version of the multidimensional dissipation technique in eq. (26) is tested for problems (2.1)–(2.3). Only solutions of the AUSMD scheme are shown in Figures 12–14, respectively. For Section (2.1), Figure 12 shows that the carbuncle phenomenon ahead of the bow shock wave is completely removed and there is no oscillations downstream of the shock wave. For Section (2.2), Figure 13 shows that even though the dissipation technique cannot completely remove the carbuncle phenomenon, it is much weaker than that in the original AUSMD scheme. The dissipation technique can also recover severely kinked Mach stems of Section (2.3) as shown in Figure 14. However, there are still some oscillations downstream of the incident shock wave.

Figure 12:

A Mach 15 flow over a cylinder: AUSMD with κ=1.0.

Figure 13:

A Mach 6 slowly moving normal shock: AUSMD with κ=0.5.

Figure 14:

A Mach 5 shock reflection over a 46° wedge: AUSMD with κ=1.0.

4 Numerical results

In this section, the robustness against the carbuncle phenomenon of the multidimensional dissipation technique is further examined with a second-order scheme. The second-order spatial discretization is achieved by applying the least-squares method with Venkatakrishnan’s limiter [27], while the second-order Runge–Kutta time-stepping method [28] is used for a temporal discretization. The test cases include: (1) a Mach 3 flow over a forward-facing step, (2) a Mach 6 slowly moving normal shock, (3) a Mach 5.09 shock diffraction around a corner, and (4) a Mach 2 shock reflection over a 46° wedge.

4.1 Mach 3 flow over a forward-facing step

The first problem is a Mach 3 supersonic flow over a forward-facing step [29]. The height of the channel and step are 1 and 0.2 units, respectively, and the step located at 0.6 unit downstream of the inlet. The forward-facing step generates a detached shock wave that then reflects between the upper and lower walls as shown in Figure 15. A structured triangular grid of Δx=Δy=0.0125 is used for this problem. Figure 15 shows the density contours from the AUSMD, AUSMV, and AUSMD schemes. The AUSMD scheme generates a small circulation, which disturbed the solution and may introduce some oscillations to downstream of the detached shock wave. Figure 16 shows the solution from the AUSMD scheme with κ=1.0. It can be seen that the oscillations were now significantly removed.

Figure 15:

A Mach 3 flow over a forward-facing step: AUSMD (top), AUSMV (middle), and AUSMDV (bottom).

Figure 16:

A Mach 3 flow over a forward-facing step: AUSMD with κ=1.0.

4.2 Mach 6 slowly moving normal shock

The problem of a planar shock wave propagating through a two-dimensional duct is tested again. Figure 17 shows the density contours of the original schemes. The AUSMD scheme produces the most serious carbuncle phenomenon and the normal shock wave moves faster than those produced by the AUSMV and AUSMDV schemes. The post-shock oscillations are clearly seen in the second-order computation of all schemes due to low numerical dissipation. The carbuncle phenomenon is totally removed but slight oscillations behind normal shock waves are still present when the multidimensional dissipation technique with κ=3.0 is applied to all schemes as shown in Figure 18. Figure 19 shows the comparison of the computed density distribution along the centerline with the exact solution at time t=0.36. We can see that the post-shock oscillations are significantly improved by introducing the multidimensional dissipation technique. The original AUSMD and AUSMD with κ=3.0 schemes can capture the moving shock with two intermediate points. The original AUSMDV and AUSMDV with κ=3.0 schemes can capture the moving shock with three intermediate points due to their more dissipative flux blending functions. However, post-shock oscillations are common to shock-capturing finite-volume methods at least if one tries to resolve shock waves as sharply as possible. Arora and Roe [30] reported that slowly moving shock waves would generate spurious oscillations in a code employing either the Godunov or Roe fluxes. There is no completely satisfactory solution to this problem. Fedkiw et al. [31] proposed a method to alleviate spurious oscillations in the conservative schemes for multiphase and multispecies Euler equations with a non-conservative modification of the total energy computed by solving a coupled evolution equation for the pressure.

Figure 17:

A Mach 6 slowly moving normal shock: AUSMD (top), AUSMV (middle), and AUSMDV (bottom).

Figure 18:

A Mach 6 slowly moving normal shock: AUSMD (top), AUSMV (middle), and AUSMDV (bottom) with κ=3.0.

Figure 19:

Comparison of density distribution along a centerline at time t=0.36: AUSMD (top) and AUSMDV (bottom).

4.3 A Mach 5.09 shock diffraction around a corner

For a strong shock diffraction around a corner [32], the structure of the diffracted shock wave consists of expansion waves, a slipstream, and a secondary shock wave. Figure 20 shows the density contours of the original schemes. Clearly, the AUSMD scheme still exhibits spurious oscillations in the vicinity of an incident shock wave, while the AUSMV and AUSMDV schemes generate slight oscillations near a tip of an incident shock wave. If the multidimensional dissipation technique is applied with κ=2.0 to all schemes, the carbuncle phenomenon is now removed, but there are still slight post-shock oscillations from the second-order AUSMD scheme as shown in Figure 21.

Figure 20:

A Mach 5.09 shock diffraction around a corner: AUSMD (top), AUSMV (middle), and AUSMDV (bottom).

Figure 21:

A Mach 5.09 shock diffraction around a corner: AUSMD (top), AUSMV (middle), and AUSMDV (bottom) with κ=2.0.

4.4 A Mach 2 shock reflection over a 46° wedge

A classical unsteady shock phenomenon, the so-called shock reflection over a wedge [33], is tested as the last example. A Mach 2 shock wave moves over a 46° wedge from left to right of the domain. The flow phenomenon consists of an incident shock wave, a Mach stem, a reflected shock wave, and a slipstream. Many existing numerical schemes produce severely kinked Mach stem [9] for this test case. The AUSMD scheme provides a kinked Mach stem solution, and a large amount of oscillations behind a Mach stem is exhibited by all schemes as shown in Figure 22. To solve this problem, the multidimensional dissipation technique is applied with κ=4.0 to all schemes. The carbuncle phenomenon and post-shock oscillations are totally removed for the AUSMD scheme, whereas there is still slight oscillations behind a Mach stem for the AUSMV and AUSMDV schemes as shown in Figure 23.

Figure 22:

A Mach 2 shock reflection over a 46° wedge: AUSMD (top), AUSMV (middle), and AUSMDV (bottom).

Figure 23:

A Mach 2 shock reflection over a 46° wedge: AUSMD (top), AUSMV (middle), and AUSMDV (bottom) with κ=4.0.

5 Conclusion

The numerical instability, the so-called carbuncle phenomenon, of the AUSMD, AUSMV, and AUSMDV schemes is investigated. Many numerical experiments show that the schemes have potentials to trigger the carbuncle phenomenon under some certain conditions, especially when the incident shock wave aligns with the grid. Post-shock oscillations are common for all schemes. From the linearized analysis of the odd–even decoupling problem, the AUSMV and AUSMDV schemes are more robust than the AUSMD scheme. The multidimensional dissipation technique adopted in the author’s previous work is modified and applied to the original schemes to heal carbuncle phenomenon and suppress post-shock oscillations. Test cases suggest that the proposed multidimensional dissipation technique is capable of removing the carbuncle phenomenon, especially where strong physical discontinuities exist, in a wide range of high-speed compressible flow problems. However, there exist post-shock oscillations downstream of the moving shock wave, at least in the situation where one tries to resolve shock waves as sharply as possible by means of higher-order schemes.