A hybrid numerical method and its application to inviscid compressible flow problems

Computer Physics Communications 185 (2014) 479–488

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

A hybrid numerical method and its application to inviscid compressible flow problems Ghislain Tchuen a,∗ , Ferdinand Fogang b , Yves Burtschell c , Paul Woafo b a

University of Dschang, IUT-FV, LISIE, P.O. Box. 134 Bandjoun, Cameroon

b

University of Yaounde I, Faculty of Science, LaMSEBP, P.O. Box 8210, Yaounde, Cameroon

c

Aix-Marseille University, IUSTI-UMR 7343, 5 rue Enrico Fermi, Technopole de château Gombert, 13453 Marseille Cedex 13, France

article

info

Article history: Received 17 September 2012 Received in revised form 31 August 2013 Accepted 3 October 2013 Available online 10 October 2013 Keywords: AUFS Roe AUFSR Hybrid method Shock instabilities Carbuncle

abstract An improved version of the artificially upstream flux vector scheme, is developed to efficiently compute inviscid compressible flow problems. This numerical scheme, named AUFSR (Tchuen et al. 2011), is obtained by hybridizing the AUFS scheme with Roe’s solver. This approach handles difficulties encountered by the AUFS scheme, in the case where the flux vector does not check the homogeneous property. The present scheme for multi-dimensional flows introduces a certain amount of numerical dissipation to shear waves, as Roe’s splitting. The AUFSR scheme is not only robust for shock-capturing, but also accurate for resolving shear layers. Numerical results for 1D Riemann problems and several 2D problems are investigated to show the capability of the method to accurately compute inviscid compressible flow when compared to AUFS, and Roe solvers. © 2013 Elsevier B.V. All rights reserved.

1. Introduction In a progressive effort to solve complex flow problems described by Euler/Navier–Stokes equations, the research on how to maximize both accuracy and efficiency has been the primary goal of designing an algorithm in numerical analysis. The compressible flow problems involve complex flow phenomena, such as strong shock waves, shock–shock interactions and shear layers. A number of numerical flux functions for inviscid fluxes have been devised as approximate solutions to the Riemann problem. Among these formulations, upwind numerical methods have become popular for solving hyperbolic partial differential equations with discontinuous solutions. They take into account the physical properties of the flows into the numerical formulation, and their essential characteristics are their particular treatment of the convective terms into Navier–Stokes equations from a well adapted flux decomposition. These are usually classified as either flux difference splitting (FDS) or flux vector splitting (FVS), each method having its advantages and disadvantages. The FDS scheme is based on the difference between the decomposition of fluxes, constructed on an approximated solution of the

∗

Corresponding author. Tel.: +237 99 93 04 89; fax: +237 33 01 46 01. E-mail address: [email protected] (G. Tchuen).

0010-4655/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cpc.2013.10.002

local Riemann problem between two adjacent states. Several formulations of FDS schemes have been compared in the literature [1–5]. These methods have been proven to be capable of capturing sharply and correctly the discontinuities of all types such as shock waves corresponding to linear as well as non-linear waves. Unfortunately, certain FDS schemes contain subtle flaws and produce spurious solutions like low frequency post shock fluctuations in the case of slowly moving shock, carbuncle phenomena, odd–even decoupling and kinked Mach stem on certain occasions [2,6–10]. Even though the FVS methods [11–13] when modified to improve their capability to capture contact discontinuities are not spared from shock instabilities and carbuncle phenomena [9,14], they rely on a decomposition of the flux vector into the positive and negative components and according to the sign of the propagation of the associated waves. They are found to be very successful in capturing steady discontinuities represented by nonlinear waves, which include shocks. However, they are not effective in capturing the discontinuities represented by linear waves, which result in incorrect diffusion of the contact surface and shear waves and a high dissipation in strong rotational flows [9,15]. However, when compared with Godunov-type schemes, the FVS results are poorer resolutions of discontinuity, particularly contact discontinuity [13]. Most upwind schemes, either Godunov-type or FVS methods, have difficulties in resolving the sonic point, and produce a spurious expansion shock there. Roe’s [2] approach, which is an FDS scheme, is widely used because of its accuracy, quality and mathematical clarity. However,

480

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

the scheme may sometimes lead to unphysical flow solutions in certain problems. They admit rarefaction shocks that do not satisfy the entropy condition. This flaw can be easily handled by simple entropy fix procedures in a one-dimensional case. To improve the accuracy of solution of these problems, Quirk [8] also pointed out that Roe’s original scheme should be modified or replaced by other schemes in the vicinity of a strong shock. The combination of Roe’s scheme with other solvers, would therefore be promising. The AUFS method, which is a special FVS scheme, was proposed by Sun and Takayama [13] for splitting flux vectors of Euler equations. This scheme has recently been extended [16] to calculate two-dimensional hypersonic flow, in thermochemical nonequilibrium. The AUFS scheme introduces two artificial wave speeds into the flux decomposition. The direction of wave propagation is adjusted by these two wave speeds. One part of the flux is numerically obtained following the Steger–Warming approach, which is an FVS method. Steger and Warming [11] were the first to use the homogeneous property of the governing equations of gas, and expressed the inviscid flux vectors in terms of their Jacobian matrices. This method is simple, accurate and robust for the Euler equations but it is not applicable to non-homogeneous equations (e.g., the magnetohydrodynamics (MHD) equations). Hence, the current approach is motivated by the desire to combine the efficiency of FVS schemes and the accuracy of FDS schemes. In this paper, the flux splitting scheme named the AUFSR scheme [17], is obtained by hybridizing the AUFS solver (FVS) and Roe’s scheme (FDS). The motivation in using Roe’s scheme relies on several advantages, which make it well known. The resulting flux function can be implemented in a very simple manner, in the form of Roe’s solver with modified wave speeds, so that converting an existing AUFS flux into the new fluxes is an extremely simple task. The procedure is to avoid difficulties encountered by the AUFS scheme, when the artificial diffusion is computed with the Steger–Warming approach in the case where the flux vector does not check the homogeneous property. The AUFSR scheme should be able to solve flows with non-homogeneous properties, furthermore, this type of hybrid method can be easily extended to a higher order of accuracy through the formal use of MUSCL interpolations. In this paper, a hybridized method to construct a simple, accurate and robust scheme is presented and used for solving Euler equations. Numerical examples are given to demonstrate that the hybrid scheme has a high computational accuracy when compared to AUFS, Roe’s and exact Riemann solvers. 2. Preliminaries

where λ = ∆t /∆x, ∆t is the time step, and ∆x is the grid size. Fi+1/2 is the numerical flux vector defined by two neighboring cells (or to the left or right cell), F (Ui+1/2 ) = F (Ui , Ui+1 ) = F (U L , U R ).

In the finite-volume formulation, the difference among all the numerical schemes lies essentially in the definition of the numerical flux Fi+1/2 evaluated at the cell interface. The numerical flux will vary, depending on the properties it must meet. The system of Eqs. (1) can be expressed in quasi-linear form Ut + AUx = 0

∂F ∂U + = 0 t > 0, −∞ ≤ x ≤ ∞ ∂t ∂x U (x, 0) = Uo (x)

∂F ∂U

A=

F = uU + P

3. General form of the AUFS solver The fundamental idea of the AUFS scheme, is to split the flux vector (2) as follows [13]: F = (1 − M )[(u − s1 )U + P ] + M [(u − s2 )U + P ]

where P = (0, p, pu) . For the numerical solution of (1), we shall consider piecewise a constant approximation of Uin+1 defined by the explicit three-point scheme in conservative form,

(7)

where s1 and s2 are two scalar constants. M = s1 /(s1 − s2 ), and (7) can be rewritten as F = (1 − M )F1 + MF2

(8)

with F1,2 = (u − s1,2 )U + P. These two flux vectors are different from the original by a term −sU. Their Jacobian matrices are A1,2 = ∂ F1,2 /∂ U = A − sI and their corresponding matrices of eigenvalues become Λ1,2,3 = diagonal(u − s − c , u − s, u − s + c ). An excellent merit of the AUFS scheme is that the eigenvalues of λ1,2,3 can be changed by varying the scalar value s. s is an artificially introduced wave speed. Depending upon the value and sign of s1 that one gets 1 L (P + P R ) + δ U and F2 = U α (uα − s2 ) + P α . (9) 2 The performance of the scheme should rely on the way to express δ U which represents the artificial viscosity. In the AUFS scheme, this viscosity is determined by the Steger–Warming formula. Substituting (8) and (9) to (7), the following final intercell flux is obtained F1 =



1 2

 (P L + P R ) + δ U + M [U α (uα − s2 ) + P α ]

(10)

where



(1)

(2)

(6)

is the Jacobian matrix. If A has a complete set of linearly independent eigenvectors, the eigenvalues associated with one subvector can be all positive and those associated with the other all negative. The system (1) is homogeneous because it satisfies F (U ) = A(U )U.

α=

where U and F are the vectors of conserved variables and fluxes, given respectively by, U = (ρ, ρ u, E )T and F = (ρ u, ρ u2 + p, (E + p)u)T . Here ρ is density, p is pressure, u is particle velocity and E is the total energy per unit volume defined as, E = ρ(e + u2 /2). The ideal-gas equation of state is assumed, p = (γ − 1)ρ e, where γ = 1.4. The flux vector can be rewritten as

(5)

where

F = (1 − M )

Consider as an initial-value problem, the one-dimensional (1D) system of conservation laws for ideal-gas flows,

(4)

L for s1 > 0, R for s1 ≤ 0.

(11)

The artificial numerical velocities s1 and s2 make it possible to write the intercell flux. The expressions proposed by Sun and Takayama [13] are given as follows. s1 =

1 2

( uL + uR )

(12)

and

 s2 =

min(0, uL − cL , u∗ − c ∗ ) max(0, u∗ + c ∗ , uR + cR )

for s1 > 0 for s1 ≤ 0

(13)

with the speed u∗ and the sound speed c ∗ given by [13]

T

Uin+1

=

Uin

− λ(Fi+1/2 − Fi−1/2 ) n ∈ N , i ∈ Z

(3)

u∗ = c∗ =

1 2 1 2

(uL + uR ) +

cL − cR

γ −1 1

;

(cL + cR ) + (γ − 1)(uL − uR ). 4

(14)

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

¯ | is a diagonal matrix defined by with L¯ = R¯ −1 , and |Λ

4. Basis of Roe’s approach Among many upwind schemes successfully applied to hyperbolic equations with discontinuous solutions, Roe’s method stands out by the simplicity and clarity of the underlying physical model. Roe’s scheme is the basis of many numerical schemes for compressible flows. However, the scheme often encounters catastrophic solutions such as non-physical expansion shocks, non-positive density, and shock instabilities related to the carbuncle phenomena [7]. Roe solved the Riemann problem (1) approximately by introducing the Jacobian matrix (6). The numerical interface flux is often written in the following form Fin+1/2 =

1

i+1/2

Fin + Fin+1 − Σk |λk

2

i+1/2 i+1/2 Rk

|ηk



i+1/2

i+1/2

(R−1 )k

.

i+1/2

and averaged wave strengths ηk

=

Roe-averaged Jacobian matrix.

When the nonlinear flux vector F (U ) is a homogeneous function of degree one in U, F (U ) can be split into subvectors, each one of which is associated with a tailored set of eigenvalues. The Steger–Warming splitting approach is no longer applicable for the investigation of first-order conservative systems that are nonhomogeneous. The basic idea used here is to linearize the Steger–Warming matrix by another formulation as follows. Firstly, the AUFSR approach is based on a linearization of the intermediate intercell flux F1 = (u − s1 )U + P, using a hybridization with Roe’s solver. Secondly, to resolve contact discontinuities efficiently, the two artificial velocities s1 and s2 take a new formulation, calculated from the Roe-averaged values. Therefore, the resulting flux keeps the properties of the AUFS scheme, which introduces some dissipation to shear waves. The flux F1 constructed from a Roe-averaged state U¯ satisfies the following

∆F1 = A¯1 ∆U

(16)

where ∆F1 = F1 (U ) − F1 (U ), ∆U = U − U , and A¯1 = ∂ F1 /∂ U evaluated by the Roe-averaged quantities, L

R

L

  ρ L uL + ρ R uR   ; ρL + ρR (17)      12 ρLH L + ρRH R 1 2 ¯ ¯   ; c¯ = (γ − 1) H − u¯ . H = 2 ρL + ρR

 ρ¯ = ρ L ρ R ;

u¯ =

The intermediate intercell flux F1 (U ) is given by, F1AUFSR =

 1 1 ¯ |L¯ ∆U F1 (U L ) + F1 (U R ) − R¯ |Λ 2 2

(18) ∂F

where R¯ is the right eigenvector matrix of A¯1 = ∂ U1 ,



1



R¯ 1 =  u¯ − c¯  ; ¯ − u¯ c¯ H



1



R¯ 3 =  u¯ + c¯  ¯ + u¯ c¯ H



F1AUFSR =



1 u¯  ¯R2 =  1  ; u¯ 2 2

3  1 1 L F1 + F1R − |λ¯ k |ηk R¯ k 2 2 k=1

(21)

After having replaced F1 by its expression, we get the following relation

 1 L  1 L P + PR + (u − s1 )U L + (uR − s1 )(U R ) 2 2 −

3 1

2 k=1

|λ¯ k |ηk R¯ k

(22)

The characteristic variables satisfy the relations

η1 = η3 =

5. The hybrid scheme, AUFSR

(20)

This can be rewritten in the following form

(15)

∆U. All these quantities are expressed directly from a

R

¯ | = diagonal(|λ¯ 1 |, |λ¯ 2 |, |λ¯ 3 |) |Λ = diagonal(|¯u − s1 − c¯ |, |¯u − s1 |, |¯u − s1 + c¯ |).

F1AUFSR =

In order for the approximate method with the Roe-type numerical flux (15) to be completely determined, it is necessary to calculate the average eigenvalues λki+1/2 , the corresponding averaged right eigenvectors Rk

481

1 2c¯ 2 1 2c¯ 2

[∆p − ρ¯ c¯ ∆u],

η2 = ∆ρ − ∆p/¯c 2 , (23)

[∆p + ρ¯ c¯ ∆u].

While identifying (9) to (21), the artificial viscosity of this scheme is the following:

δ U AUFSR =

((uL − s1 )U L + (uR − s1 )U R ) 2

−

3 1

2 k=1

|λ¯k |ηk R¯ k .

(24)

This relation (24) is very simple to handle. It introduces another artificial diffusion around a contact discontinuity. The intermediate intercell flux F2 in the proposed scheme has the same formulation as the one described for the relation (9). After substituting (9) and (22) into (7), the final intercell flux can be expressed as follows F = (1 − M )



1 2

(P L + P R ) + δ U AUFSR



+ M [U α (uα − s2 ) + P α ] .

(25)

In order to determine completely the numerical fluxes, the values of the artificial velocities must be known. They are the lower and upper bounds for the physical wave speeds with which the information of the initial discontinuity is transported. To satisfy the entropy and the positivity conditions, Einfeldt [18] suggested adequate bounds by making use of the Roe-averaged eigenvalues. Then, the approximate wave speeds in this work, are simply set to their Roe’s average s1 = u¯ .

(26)

Numerical values of s2 can be computed from min(0, uL − cL , u¯ − c¯ ) for s1 > 0, max(0, u¯ + c¯ , uR + cR ) for s1 ≤ 0.

 s2 =

(27)

If c L = c R and uL = uR , the numerical parameters are reduced to u∗ = u¯ and c ∗ = c¯ , and the artificial velocities (s1 , s2 ) become the same for both schemes (AUFS and AUFSR). The accuracy of this AUFSR solver will be studied in this paper with the help of the tests of applications to 1D and 2D problems. 6. Extension to multi-dimensions and to higher order accuracy

(19)

The catastrophic carbuncle instability for numerical algorithms is a particular problem for multidimensional computation. Therefore, in the construction of numerical fluxes, it is helpful to take the multidimensional effects into consideration.

482

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

The time-dependent Euler equations in two dimensions are given in a very compact notation by Ut + F (U )x + G(U )y = 0

(28)

with

 ρu 2  ρu + p  F = ; ρ uv  u(E + p)

 ρ ρ u U =  ; ρv 



E

(29)

 ρv  ρ uv  G= 2 . ρv + p  v(E + p) 

The numerical scheme employs the finite volume method for any grid system. It is shown in [19] that the Euler equations satisfy the rotational invariance property. Given an interface with normal vector n = (nx , ny ) and the rotational matrix T, it is important to

˜ = G(TU ). calculate U˜ = TU , F˜ = F (TU ) and G ˜ With U = TU = (ρ, ρ q˜ n , ρ q˜ t , ρ e)T , F˜ = F (TU ) = (ρ q˜ n , ρ q˜ 2n + p, ρ q˜ n q˜ t , ρ eq˜ n + pq˜n )T , where q˜ n is the normal velocity through the interface (˜qn = unx + v ny ), and q˜ t the tangential velocity (˜qt = −uny + v nx ). Thus, a numerical flux across each interface may be determined by solving exactly or approximately a Riemann problem based on the one-dimensional form in the direction of the normal face n. U˜ t + F˜n = 0.

(30)

The Eq. (30) contains four wave speeds q˜ n − c¯ , q˜ n , q˜ n , q˜ n − c¯ . After having applied the rotational invariant matrix, and rotating back the flux, the final numerical flux (10) through the interface can be expressed as



F = (1 − M )



1 2



 (P L + P R ) + δ U + M [U α (˜qαn − s2 ) + P α ]

(31)

¯ ¯ where δ U = 12 (˜qLn − s1 )U L + (˜qRn − s1 )U R − 21 k=1 |λk |ηk Rk and T P = (0, pnx , pny , pq˜ n ) . Symbols α, c¯ , s1 and s2 are defined as those in one dimension, (11), (17), (26) and (27), respectively, but using normal velocities. This technique can be readily extended to three dimensions. Eq. (28) is discretized in a finite volume approach and details can be found in [16]. After integrating (28), over an arbitrary and discrete volume, the system becomes a system of ordinary differential equations and takes the following form: 

4 ∂U 1  − →→ + (F − n )k = 0 ∂t Ω i ,j k = 1



4

(32)

In high order, the scheme can be extended following Van Leer’s MUSCL approach [20], which verifies TVD properties. The fluxes at the interfaces are extrapolated between neighboring cells by F (Ui+1/2,j ) = F (Ui,j , Ui+1,j )

(33)

The limiter minimum-modulus (MINMOD) function is used. The time predictor–corrector algorithm is used to obtain second-order time accuracy. 7. Numerical results and discussion The performances of the AUFSR scheme are illustrated on the 1D Riemann problems and 2D-compressible flows for ideal gas with γ = 1.4.

Table 1 Data for seven test problems. Test

ρL

uL

pL

ρR

uR

pR

1 2 3 4 5 6 7

1.0 1.0 1.0 5.99924 1.0 1.4 1.4

0.75 −2.0 0.0 19.5975 −19.5975 0.0 0.1

1.0 0.4 1000.0 460.894 1000.0 1.0 1.0

0.125 1.0 1.0 5.99242 1.0 1.0 1.0

0.0 2.0 0.0 −6.19633 −19.59745 0.0 0.1

0.1 0.4 0.01 46.0950 0.01 1.0 1.0

7.1. Results for 1D problem tests Seven tests were solved for 1D problems, which are generally tested by Toro [19]. These tests have exact solutions. In all chosen tests given in Table 1, data consists of two constant states WL = (ρ, u, p)L and WR = (ρ, u, p)R , separated by a discontinuity at a position x = x0 . The exact and numerical solutions are found in the spatial domain 0 ≤ x ≤ 1 using 100 cells. The Courant number coefficient is taken as Ccfl = 0.9. For each test problem, an initial location of discontinuity, x0 , and the output time are selected, these are stated in the caption of each figure. Boundary conditions are transmissive. For the sake of clarity, the numerical results of the AUFSR scheme are compared with those of the exact solution, of Roe’s scheme, of the Steger–Warming scheme and of the AUFS scheme. Test 1 is a modified version of Sod’s problem [21]. The solution consists of a right shock wave, a right travelling contact wave, and a left expansion wave with a sonic point inside. This test is useful in assessing the satisfaction of the entropy property of numerical methods. The computed results for test 1 are shown in Fig. 1 for the five schemes, where the solid line denotes the exact solution. As expected, Roe’s solver exhibits an expansion shock near the sonic point because an entropy fix is not used. The hybrid scheme (AUFSR) seems to be the most accurate near the sonic point. It resolves the sonic rarefaction point better than the AUFS, Roe’s and the Steger–Warming schemes. It resolves the left travelling waves and the contact wave more smoothly than the other three schemes, because it implicitly introduces a certain amount of important non-linear viscosity by slightly over-estimating the artificial wave speed s2 , with which the Roe-averaged values are associated. The AUFSR scheme even captures the right travelling shock slightly more sharply than the other schemes. The exact solution of test 2 consists of two symmetric expansion waves and a trivial contact wave of zero speed; the star region between the non-linear waves is close to vacuum, which makes this problem a suitable test for assessing the performance of numerical methods for low-density flows. The results from the AUFSR scheme for this test are shown in Fig. 2. Roe’s scheme fails on this test. The evolutions in pressure from the other three schemes are virtually identical. The numerical schemes give slightly more diffusion of internal energy. The AUFSR algorithm is less dissipative than the AUFS scheme. The AUFSR scheme gives more accurate results than those obtained with the AUFS scheme; in the vicinity of the trivial contact, where both density and pressure are close to zero, the results are somewhat erroneous, see internal energy plots. One observes that all numerical fluxes exhibit this error. Test 3 is designed to assess the robustness and accuracy of numerical methods. Test 3 is a very severe test; its solution consists of a strong running shock wave of shock Mach number 198, a contact surface and left expansion waves. Fig. 3 shows the results for five schemes. Compared to the exact solution, density levels, remain relatively low behind the shock wave, as is clearly seen in the density plots. The results from the AUFSR scheme are also more accurate than the results obtained with the other schemes. The Steger–Warming scheme has a non-physical bump around the contact wave, which is clearly seen in velocity distributions.

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

483

Fig. 1. Solutions of test 1 with x0 = 0.3. Exact (solid line) and numerical solutions are compared at the output time 0.2 U.

Fig. 2. Solutions of test 2 with x0 = 0.5. Exact (solid line) and numerical solutions are compared at the output time 0.15 U.

Fig. 3. Solutions of test 3 with x0 = 0.5. Exact (solid line) and numerical solutions are compared at the output time 0.012 U.

Test 4, like test 3, is also designed to assess the robustness of numerical algorithms, data originates from two very strong shock waves travelling towards each other. The solution consists of three strong discontinuities moving to the right. The left shock wave moves to the right very slowly, which adds another difficulty to numerical methods. The results from the AUFSR scheme, shown in Fig. 4, are overall comparable with those obtained with other

schemes. The four schemes behave similarly well, with a slight improvement in results obtained with the AUFSR scheme. The hybrid AUFSR scheme resolves the right-travelling shock slightly more sharply than the AUFS and Roe’s schemes. Test 5 is effectively test 3, with a negative uniform background speed added so as to obtain a stationary contact discontinuity. Test 5 is also designed to test the robustness of numerical methods but

484

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

Fig. 4. Solutions of test 4 with x0 = 0.4. Exact (solid line) and numerical solutions are compared at the output time 0.035 U.

Fig. 5. Solutions of test 5 with x0 = 0.8. Exact (solid line) and numerical solutions are compared at the output time 0.012 U.

the main reason for devising this test is to assess the ability of the numerical methods to resolve slow-moving contact discontinuities. The exact solution of test 5 consists of a left rarefaction wave, a right-travelling shock wave and a stationary contact discontinuity. Results are shown in Fig. 5. The contact discontinuity is heavily smeared in the results of the Steger–Warming scheme. The AUSFR scheme performs very well with better accuracy than the other schemes. The ability of the AUFSR scheme to resolve contact discontinuity is more straightforwardly demonstrated in Fig. 6,

in which initial conditions give an isolated contact wave. It concerns the exact recognition of isolated contacts. Except for the Steger–Warming scheme, all the other schemes of this paper pass the test successfully. The performance of the AUFSR scheme in resolving moving contact waves is accessed with the numerical result given in Fig. 7. All numerical methods presented lead to identical results. In order to verify the performance of the AUFSR scheme, few 2D problems are solved to illustrate its accuracy and robustness.

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

485

Fig. 8. Perturbation at the centerline (zoomed view). Fig. 6. Resolution of stationary contact wave. Exact (solid line) and numerical methods are compared at time 2.0 U.

Fig. 9. The odd–even grid perturbation (density contours) Roe scheme.

Fig. 10. The odd–even grid perturbation (density contours) exact Riemann scheme. Fig. 7. Resolution of moving contact wave. Numerical methods are compared at time 2.0 U.

Numerical results for two-dimensional Euler equations are computed using the same mesh, the same limiter and the same boundary conditions for different schemes. 7.2. Odd–even grid perturbation problem This test case corresponds to a problem first proposed and investigated by Quirk [8], and now becomes one of the significant examples demonstrating the carbuncle phenomenon. It is well known that many upwind schemes, including the exact Riemann solver and the approximate ones, which explicitly capture the 1D contact wave suffer from this problem of odd–even decoupling. This case is used to investigate the performance of AUFSR schemes. A plane shock wave travels downstream in a straight duct at the speed of Mach 6. The computational domain consists of a uniform mesh with 21 × 800 points respectively along the axial and the transverse direction of the duct. The centerline grid points are modified with involvement of the following grid perturbations δ y = ±10−3 , as sketched in Fig. 8. This perturbation promotes odd–even decoupling along the length of the shock. The initial position of the shock wave is at x = 20. To the right of the shock, the initial conditions are (ρ, u, v, p) = (1.4, 0, 0, 1). The flow is inviscid. To the left of the shock, the flow variables are computed using Rankine–Hugoniot relations. The propagating shock is

Fig. 11. The odd–even grid perturbation (density contours) AUFS scheme.

Fig. 12. The odd–even grid perturbation (density contours) AUFSR scheme.

strong enough for the flow to enter the duct and the inlet boundary condition is considered as supersonic. The non-reflecting simple wave boundary conditions are imposed at the outlet. The upper and lower boundaries are considered as solid walls. All computations are carried out in the first-order scheme. Figs. 9–12 show the computed density contours of the normal shock at the three locations along the duct by the solvers of Roe, exact Riemann, AUFS and AUFSR, respectively.

486

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

Fig. 13. M∞ = 10 flow over a cylinder: pressure contours. (a) Roe, (b) AUFS, (c) AUFSR.

As shown in Figs. 9 and 10, the shock was destroyed by Roe’s flux and exact Riemann schemes, and a typical carbuncle was developed. Although a planar shock wave is recovered, many disturbances still prevail behind the shock wave in comparison with the uniform solution given from the AUFSR scheme. The AUFSR and AUFS schemes show in Figs. 11–12 that decoupling was completely eliminated, and successfully preserved the initial shock. In this test case, the result of the AUFSR is almost identical to that obtained with AUFS.

Fig. 14. Convergence history for the blunt body problem.

7.3. M∞ = 10 flow over a circular cylinder This is another well-known test case to examine the catastrophic carbuncle failings of certain Riemann solvers based on the FDS scheme. Such a phenomenon refers to a spurious bump on the bow shock near the flow center line ahead of the blunt body. The free stream Mach number was taken to be M∞ = 10 flow over a circular cylinder. This test problem is computed on a structured grid with 160 cells in the circumferential direction and 80 cells in the radial direction. The boundary conditions used are similar to those used in [9]. The first-order scheme was used with the three solvers presented. As expected, Roe’s flux created a distorted solution as shown in Fig. 13(a). On the other hand, the AUFS fluxes did not produce in Fig. 13(b) such a kind of shock instability. The AUFSR flux did not allow the carbuncle to appear and produced a correct solution with a clean flow field behind the bow shock (see Fig. 13(c)). This indicates that generally the carbuncle phenomenon is suppressed by adding an extra dissipation, but adding too much also provokes the carbuncle [10]. The solutions are fully converged for all cases (CFL = 0.6). The convergence history is displayed in Fig. 14, exhibiting that the residuals of the Roe solution have stagnated after six orders of magnitude reduction, while the AUFS and AUFSR residuals continue to decrease to machine accuracy, following roughly the same history. 7.4. Mach 1.7 shock reflection over a wedge The shock wave phenomena considered here are the reflection of an incident shock wave of Mach number 1.7 over a wedge of 25 degrees to the incident shock. The computational domain is [0, 25.0] × [0, 16.5] on the x–y plane, with the apex of the wedge placed at x = 4.69, and the shock placed initially at x = 4.0. To the right of the shock, the flow field is initialized with the ambient conditions (ρ = 1.225 kg/m3 , p = 101, 325Pa , u = v = 0). To the left of the shock, the flow variables are computed using the Rankine–Hugoniot relations. The computation is carried out on a uniform mesh with 800 × 600 cells using first-order schemes. The computed density contours are shown in Figs. 15 and 16. The AUFS

Fig. 15. Computed density contours with the AUFS scheme.

Fig. 16. Computed density contours with the AUFSR scheme.

solver is used in Fig. 15, and the AUFSR solver in Fig. 16. It is noted that, the proposed scheme has a better prediction of the flow field. The shock wave reflection pattern that emerges corresponds to a single Mach reflection. There is the formation of a triple point composed of three shocks, namely, part of the incident shock, the reflected shock and the Mach stem. From the triple point there emerges a slipstream that joins the wedge surface at a sharp angle. The flow feature is captured very well and produced a perfectly comparable solution to the experimental result [15], shown in Fig. 17. The qualitative agreement between the experimental result and the computations is excellent as demonstrated with the numerical holographic interferogram in Fig. 18, computed with the first order of the AUFSR scheme. The result with the AUSFR scheme matches well with the experimental interferogram.

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

487

Fig. 17. Takayama experimental result, holographic interferogram Ms = 1.7, wedge at 25°.

Fig. 19. Shock diffraction, density contours with first-order exact Riemann scheme.

Fig. 18. Numerical holographic interferogram of Mach reflection of an incident shock of Mach number 1.7 over a wedge at 25°, with the AUFSR scheme.

7.5. Shock diffraction The last problem considered here is the diffraction of a supersonic planar shock moving over a 90° corner. The shock Mach number is 5.09. Quirk [8] has investigated this test case, and has shown the complexity of the flow, which generates a series of the complex shock diffraction, reflection, and interaction patterns. This problem is another test case for which many Godunov type fluxes suffer from the carbuncle. Quirk [8] pointed out that some Riemann solvers encountered difficulties on a fine mesh and thus placed an upper limit on the resolution of simulations. The geometry of the computational domain is a unit square [0, 5] × [0, 5] that is discretized into a 600 × 600 uniform cell. The step is located at x = 0.1 and y = 0.2. The top boundary is taken as a wall, and the right and bottom boundaries are taken as outflow. The solid walls are treated using the reflecting boundary condition. Despite the rather fine mesh, the calculation of the density contours with the exact Riemann solver in first order suffers from slight shock instability as shown in Fig. 19. This slight instability observed on the incident shock disappears as shown in Figs. 20 and 21 for the results obtained respectively by using the second order of AUFS and AUFSR schemes. Fig. 21 shows the density contours obtained by using the AUFSR scheme with the same number of contours as in Figs. 19 and 20. All of the major flow structures named by Skews [22] can be easily identified in Fig. 20. The recompression shock, secondary shock, slipstream, contact surface, incident shock and expansion wave can all be seen as sharply defined density gradients. The condition of the incident Mach number chosen here cannot allow us to visualize the vortex, that is well defined for Ms < 1.5 [22]. The pressure contours are shown in Fig. 22. It can be seen that the shock can be sharply captured with a good resolution with the AUFSR scheme without any anomaly. 8. Conclusion A hybridized flux method is proposed to improve the AUFS scheme. The method is developed by combining Roe’s solver and

Fig. 20. Shock diffraction, density contours with second-order AUFS scheme.

Fig. 21. Shock diffraction, density contours with second-order AUFSR scheme.

the AUFS. The resulting flux functions are not only robust for linear shock instability, but also accurate for resolving shear layers, sonic points smoothly and contact discontinuities. To further improve solution accuracy, a higher-order spatial coupled to a time predictor–corrector algorithm was also implemented. The entire procedure was found to provide more accurate solutions for both steady-state and transient flow test cases.

488

G. Tchuen et al. / Computer Physics Communications 185 (2014) 479–488

Fig. 22. Shock diffraction, pressure contours with second-order AUFSR scheme.

Numerical results for 1D Riemann problems show that the AUFSR scheme has a high computational accuracy when compared to the AUFS scheme and to Roe’s solver. Several 2D problems are computed to show the capability of the method to compute 2D compressible flows. The robustness and accuracy demonstrated by these series of numerical experiments indicate that this hybrid AUFSR solver offers a quick and effective cure to those multi-dimensional finite-volume Euler codes that suffer from the carbuncle phenomenon. Exploring the method for non-homogeneous problems such as MHD flow is left as a future work. References [1] S.K. Godunov, Finite difference methods for the numerical computation of discontinuous solutions of the equations of fluid dynamics, Mat. Sb. 47 (1959) 271. [2] P.L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys. 43 (1981) 357–372.

[3] P. Woodward, P. Collela, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54 (1984) 115–173. [4] A. Harten, B. Endquist, S. Order, S. Chakravarthy, Uniformly high-order accurate essentially non-ascillatory schemes, J. Comput. Phys. 71 (1987) 231. [5] C.-W. Shu, Essentially non-oscillatory and weighted essentially nonoscillatory schemes for hyperbolic conservation laws. ICASE Report, 97–65, 1997. [6] K. Huang, H. Wu, H. Yu, D. Yan, Cure for numerical shock instability in HLLC solver, Internat. J. Numer. Methods Fluids 65 (2011) 1026–1038. [7] Y. Chauvat, J.M. Moschetta, J. Gressier, Shock wave numerical structure and carbuncle phenomenon, Internat. J. Numer. Methods Fluids 47 (2005) 903–909. [8] J.J. Quirk, A contribution to the great Riemann solver debate, Internat. J. Numer. Methods Fluids 18 (1994) 555–574. [9] M. Pandolfi, D. d’Ambrosio, Numerical instabilities in upwind methods: analysis and cures for the carbuncle phenomenon, J. Comput. Phys. 166 (2001) 271–301. [10] H. Nishikawa, K. Kitamura, Very simple, carbuncle-free, boundary-layerresolving, rotated-hybrid Riemann solvers, J. Comput. phys. 227 (2008) 2560–2581. [11] L.J. Steger, F.R. Warming, Flux Vector splitting of the inviscid gasdynamic equations with applications to finite difference methods, J. Comput. Phys. 40 (1981) 263. [12] B. Van Leer, Flux Vector Splitting for the Euler equations, Lecture Notes in Phys. 170 (1982) 507–512. [13] M. Sun, K. Takayama, An artificially upstream flux vector splitting scheme for the Euler equations, J. Comput. Phys. 189 (2003) 305–329. [14] J. Gressiera, J.M. Moschetta, Robustness versus accuracy in shock-wave computations, Internat. J. Numer. Methods Fluids 33 (2000) 313–332. [15] E.F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL–Riemann solver, Shock Waves 4 (1994) 25–34. [16] G. Tchuen, Y. Burtschell, D.E. Zeitoun, Computation of non-equilibrium hypersonic flow with Artificially Upstream Flux vector Splitting (AUFS) Schemes, Int. J. Comput. Fluid Dyn. 22 (4) (2008) 209–220. [17] G. Tchuen, F. Fogang, Y. Burtschell, P. Woafo, Hybrid upwind splitting scheme by combining the approaches of Roe and AUFS for compressible flow problems, Int. J. Eng. Syst. Model. Simul. 3 (2011) 16–25. [18] B. Einfeldt, On Godunov-type methods for gas dynamics, SIAM J. Numer. Anal 25 (1988) 357. [19] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, second ed., Springer, Berlin, 1999. [20] B. Van Leer, Towards the ultimate conservative difference scheme, III. Upstream-centered finite-difference schemes for ideal compressible flow, J. Comput. Phys. 23 (1977) 263. [21] A.G. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation law, J. Comput. Phys. 27 (1978) 1–31. [22] B.W. Skews, The perturbed region behind a diffracting shock wave, J. Fluid Mech. 29 (4) (1967) 705–719.