A comparison of hyperbolic solvers for ideal and real gas flows

COMPUTER AND FLUIDS ENGINEERING

A comparison of hyperbolic solvers for ideal and real gas flows

R. M. L. CoelhoI, ^*; P. L. C. Lage^III; A. Silva Telles^II

^ICurso de Pós-Graduação em Tecnologia de Processos Químicos e Bioquímicos (TPQB), Departamento de Engenharia Química, Escola de Química, Bloco E, Centro de Tecnologia, Universidade Federal do Rio de Janeiro, CEP 21949-900, Rio de Janeiro - RJ, Brazil

^IIResearch and Development Center (PETROBRAS /CENPES/EB/IP), Phone: +(55) (21) 3865-7406, Fax: +(55) (21) 3865-6793, Cidade Universitária, Q.7, Ilha do Fundão CEP 21949-900, Rio de Janeiro - RJ, Brazil. E-mail: [e-mail]

^IIIPrograma de Engenharia Química, COPPE, Universidade Federal do Rio de Janeiro, P.O. Box 68502, CEP 21945-970, Rio de Janeiro - RJ, Brazil

ABSTRACT

Classical and recent numerical schemes for solving hyperbolic conservation laws were analyzed for computational efficiency and application to nonideal gas flows. The Roe-Pike approximate Riemann solver with entropy correction, the Harten second-order scheme and the extension of the Roe-Pike method to second-order by the MUSCL strategy were compared for one-dimensional flows of an ideal gas. These methods require the so-called Roe's average state, which is frequently difficult and sometimes impossible to obtain. Other methods that do not require the average state are best suited for complex equations of state. Of these, the VFRoe, AUSM+ and Hybrid Lax-Friedrich-Lax-Wendroff methods were compared for one-dimensional compressible flows of a Van der Waals gas. All methods were evaluated regarding their accuracy for given mesh sizes and their computational cost for a given solution accuracy. It was shown that, even though they require more floating points and indirect addressing operations per time step, for a given time interval for integration the second-order methods are less-time consuming than the first-order methods for a required accuracy. It was also shown that AUSM+ and VFRoe are the most accurate methods and that AUSM+ is much faster than the others, and is thus recommended for nonideal one-phase gas flows.

Keywords: Hyperbolic conservation laws; Riemann solvers; Roe solver; VFRoe; AUSM+; Hybrid schemes; Real gases; Compressible flow.

INTRODUCTION

The simulation of compressible flows is important in many industrial applications and it is governed by the conservation of mass, momentum and energy and the associated constitutive equations, which are macroscopic relations such as an equation of state (EOS), which represents the behavior of the associated microscopic phenomena. These relations strongly affect the structure and dynamics of the waves in the flow.

The conservation equations for advection-dominated problems can be written as a set of hyperbolic equations whose solution for the flow variables may present large gradients or even discontinuities. Most of the existing commercial codes are based on very diffusive methods which may result in unrealistic solutions due to their failure at solving the flow wave field. On the other hand, there are several successful methods that account for the flow wave structure through a superposition of Riemann solutions. These are local solutions of the hyperbolic equations, which are composed of elementary waves (Menikoff and Plhor, 1989).

The set of hyperbolic conservation laws that models a transient one-dimensional compressible flow is usually written as

where is the physical flux and represents a source/sink term, like friction loss. The linearized form of equation (1) is given by

where is the Jacobian matrix.

Finite volume methods have proven their efficiency in solving hyperbolic initial-value problems. The most popular methods are based on the solution of the Riemann problem (RP) at the cell interfaces, which enables a very accurate estimate of the interface fluxes. The emergence of TVNI (total variation non-increasing, also popularly known as TVD, total variation diminishing) methods (Harten, 1983) made the accurate capture of discontinuities possible. The pioneering method was that of Godunov (1959), which is based on the exact solution of the RP's at the interfaces. Its advantage is the assurance of the density and pressure positivity, but its drawbacks are its high computer cost and the lack of generality due to its requirement of an analytical integration of the Riemann invariants. Therefore, approximate Riemann solvers (RS's) have been proposed of which the most popular one is that of Roe (1981). Generalization of the Roe scheme is difficult because it is not always possible to analytically obtain the expressions necessary for the average state in which the Jacobian matrix must be evaluated. The average state must have the following properties:

However, the properties given by equations (3) and (4) are not sufficient to guarantee accurate results due to a violation of entropy law, which may result in negative values in the density and pressure fields. The condition given by equation (4) is somewhat complex and it is sometimes impossible to enforce in several systems, such as complex turbulent models, two-phase models, reacting flows, etc. This has motivated the development of several alternatives to the original Roe scheme, such as VFRoe (Buffard et al., 2000), Advection Upstream Splitting Method (AUSM+) (Liou and Steffen, 1993; Liou, 1996) and hybrid methods such as the one proposed by De Vuyst (2004). These schemes are not more accurate than the popular TVD and ENO (Essentially Non-Oscillatory, Harten et al., 1987) schemes, but they are suitable for general purpose solvers.

As the ideal EOS is not able to represent nonclassical waves, many numerical methods were recently extended to real materials. Nonclassical wave phenomena are observed near the liquid-vapor saturation region or in phase-change flows. Under such conditions the fundamental derivative of gas dynamics, given by equation (5), becomes negative and the genuinely nonlinear character of the acoustic waves is lost (Guardone and Vigevano, 2002). The waves in those regions are called mixed or negative waves and the flow is said to be in the dense gas regime.

In equation (5) s is the entropy and V is the specific volume. According to thermodynamics, the term ¶P/¶V has to be negative for stability. Therefore, the negative sign of G(s,V) is associated with nonconvex isentropes in the P-V plane. The polytropic EOS is not able to describe these regions because it is convex everywhere (Menikoff and Plohr, 1989; Guardone and Vigevano, 2002).

Exact RS's for real gases were proposed by Collela and Glaz (1985), Saurel et al. (1994) and Quartapelle et al. (2003). However, the exact RS's need to integrate the Riemann invariants along isentropes in order to satisfy simultaneously all jump conditions behind a shock wave, requiring a large amount of CPU time.

A different approach based on finding a generalized Roe average state was taken by Grossman and Walters (1989), Glaister (1988), Liou et al. (1990), Vinokur and Montagné (1990), Abgrall (1991), Toumi (1992), Cox and Cinnella (1994) and Mottura et al. (1997). This average state includes the thermodynamics derivatives as additional variables and is responsible for the uniqueness of the average state. However, as stated by Toumi (1992) and Guardone and Vigevano (2002), the average thermodynamics derivatives do not retain their exact significance, which could result in inconsistencies because they are employed to derive other thermodynamic quantities, such as the speed of sound. Guardone and Vigevano (2002) obtained the average state for the Roe solver for the Van der Waals EOS. As pointed out previously, this average state is valid just for this EOS and the problem of generalization of RS's remains.

Besides the entropy violation of the original Roe solver, the VFRoe scheme (Masella et al., 1999) has an additional one: it is not consistent with the integral form of the conservation laws (Buffard et al., 2000). Buffard et al. (2000) extended the VFRoe scheme to nonconservative variables and compared the results to the original VFRoe scheme, analyzing the convergence rate for different time formulations. The authors also applied the first-order version to the Van der Waals EOS, but not the second-order extension. In the present work, the method was extended to achieve second-order and was applied to the Buffard et al. (2000) example, whose solution showed some numerical oscillations.

The Advection Upstream Splitting Method (AUSM+) for systems of hyperbolic conservation laws was proposed by Liou and Steffen (1993) and Liou (1996). It does not require the evaluation of the Jacobian of the flux, allowing easy generalization to arbitrary EOS. Besides, the method presents a low computational cost and high accuracy in the capture of contact discontinuities and it preserves the positivity of pressure and density. Therefore, AUSM+ is suitable for real gas flows. However, density-based solvers experience stiffness problems and a loss of accuracy for low Mach number flows. In these cases, AUSM+ behaves more like a central difference discretization, which may result in odd-even decoupling, and it is necessary to couple velocity and pressure (Liou and Edwards, 1999). The larger the difference between the speed of sound and the velocity, the stronger the odd-even decoupling effect is.

It is well known that the first-order schemes tend to present numerical diffusion, while the second-order schemes tend to have oscillations near high-gradient regions. De Vuyst (2004) presented a novel hybrid method. The idea of this method (called HLFW) is to evaluate the flux with a weighted combination of first-order Lax-Friedrichs and second-order Lax-Wendroff schemes.

This paper is organized as follows: section 2 is a short presentation of the basic concepts of Godunov's method and the RS's of Roe (1981), Harten (1983) and the MUSCL (Monotone Upstream-Centered Scheme for Conservation Laws) strategy of Van Leer (1977, 1979). In section 3, the AUSM+ and Hybrid Lax-Friedrich-Lax-Wendroff (HLFW) methods, which do not use any RS, are presented. In section 4, the application of these methods to ideal and real gas flows is described and the results compared. The RS methods were applied to ideal gas flows and those methods which do not require any RS were applied to Van der Walls gases. The methods are compared in terms of accuracy for a given mesh size and in respect to FLOP's (floating-point operations), IA's (indirect addressings) and CPU time, adopting different meshes for each method which result in similar accuracy in the solution.

FINITE VOLUME DISCRETIZATION, THE GODUNOV SCHEME AND RIEMANN SOLVERS

Considering a uniform mesh where Dx=x_j+1/2-x_j-1/2, the finite volume solution of the system given by equation (1) for is written as

where is the numerical flux evaluated at the interface j+1/2 and e=Dt/Dx. The numerical flux must be consistent with the physical flux, i.e., .

The idea proposed by Godunov (1959) is to solve an RP with and at the cell interfaces. These local RP's are solved exactly to give , "j where . The is then calculated by

However, this method is not attractive due to the computational cost of the exact solution of the RP.

Roe (1981) proposed the use of the solution of the linearized local RP:

where satisfies equations (3) and (4). The relation given by equation (4) is called Roe's condition, which is consistent with the integral form of conservation laws. However, Roe's solution may violate the entropy condition, also giving negative values of pressure and density. Harten (1983) introduced a class of second-order schemes which satisfy the entropy condition. Thereafter, Harten and Hyman (1983) proposed an entropy correction for Roe's solver. Other common schemes could be found in Harten (1983) and Toro (1999). One of these variants, the Roe-Pike method (Toro, 1999), will be compared to Harten's ULT1 method (Harten, 1983) in the present work.

The final expression for numerical flux in Roe's and Harten's methods is

where are the eigenvectors of ; ; F_i and F_i+1 are the fluxes evaluated at u_i and u_i+1, respectively; for the first-order Roe scheme and for the second-order Harten ULT1 scheme; , with being the eigenvalues of ; and is the k^th component of the wave strength, D_i+1/2u, in the coordinate system:

Q is the viscosity function which defines the amount of numerical viscosity. It was defined by Harten (1983) as a small perturbation of the absolute value function given by Q(x)=|x| for |x|>2z, and Q(x)=(x²/4z)+z for |x|<2z, where z is a small positive number. The remaining terms in Harten's flux definition are as follows:

For the solution based on Roe's scheme, the average state in which the Jacobian matrix is evaluated must be determined. The derivation is precisely the most important drawback of Roe's scheme due to the difficulty of its generalization to other EOS's.

The Harten-Hyman entropy correction (Harten and Hyman, 1983) is applied to transonic rarefaction waves in the Roe-Pike method. First, the eigenvalues (l_L < 0 < l_R) are examined in order to detect the existence of a transonic rarefaction wave. Then the fluxes and the intermediate states are modified (see Toro, 1999 for details).

The VFRoe scheme uses an average state different from Roe's average state with the drawback of lack of conservation. However, Buffard et al. (2000) proved its accuracy for the flow of ideal and real gases and for two-phase flows. The simplest average is the arithmetic mean, which was used here. The remaining steps of the method are the same as those in the original Roe method.

The methods presented were extended to the second-order in accordance with Van Leer's (1977, 1979) MUSCL method, later modified by Van Leer (1985) (MUSCL-Hancock). The first step is called data reconstruction and consists of changing the u value as follows:

The second step is the evolution step given by

The last step is the application of the first-order method with the solution of the RP at the "i+1/2" face with the values and as left and right values, respectively. The TVD version is obtained by setting the D_i values as follows:

where G = 1 and 2 results in the MINMOD (MINBEE) and SUPERBEE schemes, respectively.

THE AUSM+ AND THE HYBRID LAX-FRIEDRICH-LAX-WENDROFF METHODS

AUSM+ is as accurate as the flux-splitting methods of Roe (1981) without wave decomposition, result in a large reduction in CPU time consumption and making the formulation simple and quite general. AUSM+ is able to solve the contact discontinuities exactly, which is suitable for viscous flow calculations. In the case of one-phase flow, the system given by equation (1) is valid with u=[r,rv,E]^t and f=[rv, rv²+P,v(E+P)]^t. In this case the AUSM+ flux can be obtained from the following expression where the pressure and mass flux parts are separated:

where the following upwind formulae based on sign is used to compute Y^*:

using the value obtained by the following steps:

(i) Evaluation of left and right Mach numbers, calculated using a mean sound speed at the interface:

(ii) Then the interface pressure is calculated by averaging the left and right pressures:

where are consistent, differentiable and symmetric polynomial functions.

(iii) the interface Mach number is calculated as a polynomial function of the Mach numbers given in the first step:

(iv) Finally, giving the mass flux as follows:

The original AUSM+ polynomials are

where A_AUSM = 3/16 and B_AUSM = 1/8. AUSM+ suffers of odd-even decoupling for low Mach numbers, thereby losing accuracy. This is particularly important for liquid flows, where the of sound speed is three orders of magnitude higher than the velocity. In the examples studied in this paper, no strategy to deal with odd-even decoupling was necessary.

De Vuyst (2004) proposed a novel approach for numerical solution of hyperbolic systems using a hybrid method that combines first- and second-order numerical fluxes. The author presented the basis for generalizing the calculation of the optimal weighting of such fluxes. The weight is calculated from a proposed concept of local dissipation by a convexity function which differs from the traditional entropy dissipation functions because it does not need an entropy flux expression and allows some relaxation of the entropy inequality (e.g. entropy violation). De Vuyst (2004) used the Lax-Wendroff and the Lax-Friedrichs fluxes. The Lax-Wendroff flux with a relaxation parameter c can be written as

where f_L=f(u_j), f_R=f(u_j+1), and is a consistent mean state function, (u,u)=u. In the present work, constant value of c = -10^-2 and (u,v)=0.5(u+v)were used. The modified Lax-Friedrichs numerical flux is given by

Then the weighted flux at tⁿ is given by

where the weighting factor, , depends upon the dissipation function, h, which is defined as

which, in turn, depends on the numerical flux given by equation (27) and on the final updated field, . The method was implemented by the algorithm suggested by De Vuyst (2004) to achieve h < 0, which calculates as the smallest values that make negative. Let be a small increment and build the sequences as follows: let q = 0 and and

while and do

set and calculate

evaluate the numerical flux, Eq.(27)

calculate the updated field, Eq.(6)

evaluate

end while

where x is a small positive number. A Newton-Raphson-like algorithm could reduce CPU demand.

In order to make fair comparisons, and since the HLFW method is second-order accurate, AUSM+ method was extended to second-order accuracy using the same MUSCL-Hancock strategy (MUSCL AUSM+) as that applied to the VFRoe (MUSCL VFRoe) scheme.

NUMERICAL RESULTS

Application to Ideal and Real Gases

For one-dimensional compressible flow, the system given by equation (1) is valid with

where r is the density, v is the velocity, P is the pressure and E is the total energy. This system of equations must be completed by an EOS. For an ideal gas, the following relations hold:

where c is the speed of sound and g is the ratio of the specific heats (g = C_p/C_v). The Jacobian matrix is as follows:

where H=(E+P)/r is the total enthalpy. Finally, the eigenvectors and eigenvalues are

Therefore, the solution is composed of three waves with velocities l¹, l² and l³. The second wave is a contact discontinuity associated with a linear degenerated field (Ñl²×R² = 0). The other wave fields are genuinely nonlinear rarefaction and shock waves (Ñlⁱ×Rⁱ¹0, i = 1,3).

The Roe-Pike average state for this case is as follows:

from which the eigenvectors, eigenvalues and wave strengths (Eqs. 10 and 32) can be calculated.

The Van der Waals EOS is expressed as follows:

The corresponding Jacobian matrix is given by (Guardone and Vigevano, 2002):

where m=rv and

The treatment of the boundary conditions is described in detail in the ^Appendix.

Comparison of the Methods

In our comparative study, two approaches were used. First, the accuracy of the methods was compared for a given mesh. Then, different meshes were used for different methods in order to have results with approximately the same accuracy. Finally, the methods were compared with respect to FLOPs, IA's and CPU time costs. This is important because the less accurate methods could be used with a refined mesh to give approximately the same results as the more accurate methods. These refined meshes require more CPU time than the coarse meshes. However, it is not clear if the increase in CPU time due to a refined mesh is compensated for by the smaller computational cost per step for the first-order accurate methods. The answer is bound to be method- and scheme-dependent. Some examples solved in this work are fictitious test problems obtained from the literature where the variable units were not specified. For the remaining cases the units are specified in the text. For ideal gas flow, a shock tube problem and the isothermal pipeline shutting-in test were employed.

a) Shock Tube Ideal Gas Flow

Figure 1 shows the comparison between the solutions obtained for the modified Sod's shock tube problem (Toro, 1999) for a 50-cell mesh, using the Roe-Pike method with the Harten-Hyman entropy correction, the Harten method and both with two flux-limited (MINMOD and SUPERBEE) MUSCL reconstruction schemes applied to the Roe-Pike method. The exact solution, obtained using the HE-E1RPEX routine of the NUMERICA library (Toro, 1999), is also shown in Fig. 1. The initial values for density, pressure and velocity are (1, 1, 0.75) if x < 0.3 and (0.125, 0.1, 0) otherwise. These figures show that the two second-order MUSCL schemes give much better results than the first-order methods. Although not shown, it was verified that accuracy similar to that of the MUSCL results could be obtained with the Roe-Pike and Harten methods when meshes with four times more cells were used.

Figure 1 shows that the MUSCL SUPERBEE scheme gives the best solution for this test case, but followed closely by the MUSCL MINMOD scheme. In order of decreasing accuracy, they are followed by the Harten and Roe-Pike methods. It can be seen that SUPERBEE presented small errors on the rarefaction wave front, but higher peaks and some oscillation on its tail.

b) Ideal Gas Flow During Isothermal Pipeline Shutting-In

Consider now the pipeline shutting-in problem, where the friction loss is represented by g(u)=lrv²/2D. Initially, a pipeline 100 m long was filled with a stagnant ideal gas and, at t = 0, a 70 kg/m²s mass flux was imposed at the pipe inlet. When the steady state solution is obtained from a transient simulation, as in the present work, there is a cumulative error in the solution due to mass losses. Using the || uⁿ⁺¹-uⁿ || < 10^-10 criterion for assuming that the steady-state had been reached, The mean residues for each method based on the exact solution, which is the constant imposed mass flux, are presented in. One can observe that the Roe-Pike and Harten schemes achieve results with a similar accuracy. The Roe-Pike scheme required about 1000 cells to achieve accuracy similar to that of the results of the two MUSCL schemes for a 400-cell mesh. Figure 2 shows the simulation results at t = 0.1 s. Similarly to the shock tube problem, the Roe-Pike and Harten methods achieved accuracy comparable to that of the results obtained using the two MUSCL schemes when a mesh with four times more cells was used.

The only difference in the Roe-Pike and Harten methods regarding FLOP's and IA's lies in the b evaluation in accordance with equation (9), and the difference between these methods and the MUSCL schemes is due to both b and data reconstruction steps in accordance with equations (12), (13) and (14). Table 2 compares the numbers of IA's and FLOP's of these steps for each scheme. One should keep in mind that some conditional "if" operations can increase the CPU demand of the MUSCL schemes result in higher CPU times than those expected from Table 2. It is important to count the FLOP's and IA's as a function of the number of variables to be solved for and the number of time steps required for a given CFL number because they depend upon the problem. As shown previously, the Roe-Pike and Harten schemes require four times more cells to achieve accuracy comparable to that of the MUSCL MINMOD or SUPERBEE schemes at the time instant analyzed. Moreover, for a given CFL number and accuracy, the MUSCL schemes can use larger time steps because Dx is also larger. Furthermore, if the calculation is performed until the steady state is reached, the number of cells required by the Roe-Pike and Harten schemes for a given accuracy increases due to error accumulation during time integration.

Table 2 also shows a comparison of the FLOP's and IA's required by each scheme to achieve the steady-state solution for the pipeline shutting-in problem (nv = 2) for a given accuracy. If the same number of time steps were performed for all schemes, the MUSCL number of FLOP's for a 400-cell mesh would be slightly smaller (40%) than those for the Roe-Pike and Harten schemes for the 3200-cell meshes. However, the possibility of larger time steps decreases the FLOP's in the MUSCL scheme much more than those in to the Roe-Pike and Harten schemes, showing the importance of the second-order extension of these methods.

c) Real Gases

The MUSCL AUSM+, HLFW and MUSCL VFRoe methods were compared in terms of accuracy and computational cost for three flow test cases of Van der Waals gases found in the literature. The MINMOD limiter was adopted for the MUSCL scheme because the SUPERBEE limiter presented undesired oscillations near discontinuities. Test cases 1, 2 and 3 are shock tube problems defined as a RP with vector (r, v, P) given for the left and right states. The left and right states are u^L= (100 kg/m³, 0 m/s, 115931328.765 Pa) and u^R= (10 kg/m³, 0 m/s, 6935545.80532 Pa) for test case 1 (Saurel, 1994), u^L= (250, 0, 35966778) and u^R= (166.6, 0, 27114795) for test case 2 (Buffard et al., 2000) and u^L= (1.818, 0, 3) and u^R= (0.275, 0, 0.575) for test case 3 (Guardone, 2002). Additional data for the EOS are defined by the parameter vector (a, b, d) which is given by (0.138, 3.258´10^-5, 0.4), (1684.54, 1.692´10^-3, 0.3292) and (3, 0.333, 1.25´10^-2) for test cases 1, 2 and 3, respectively.

Figures 3, 4 and 5 compare the HLFW, MUSCL VFRoe and MUSCL AUSM+ results for test cases 1, 2 and 3, respectively. A 200-cell mesh was used with all methods, whose results can also be compared to those obtained with MUSCL AUSM+ using a 2000-cell mesh. Since our aim is to compare the accuracy and CPU cost of these methods, the plots show only the regions where the differences between the results are largest. The complete solution can be found in the original references.

For test case 1, the MUSCL VFRoe solution presented large deviations for v and E. It wrongly predicted the wave amplitudes, keeping these deviations in the two adjacent constant value regions. For all variables except r, the HLFW results are only a little worse than the MUSCL AUSM+ results. For r, the HLFW method provided a poor prediction of the contact discontinuity of about x = 0.29. As can be seen in Figure 3, the MUSCL AUSM+ results were the best for this test case.

For test case 2, the MUSCL AUSM+ and MUSCL VFRoe results showed some numerical oscillations near discontinuities, with the AUSM+ oscillations being smaller but having higher frequencies. Strange oscillations were observed for the MUSCL AUSM+ and MUSCL VFRoe solutions in the region 250 < x < 300. It was found that these oscillations were produced by the MUSCL strategy, as they disappeared when it was turned off, which also increased the numerical diffusion as both methods become first-order accurate. The SUPERBEE MUSCL schemes were also tested, but their oscillations were even larger. Although applied to other cases, the second-order MUSCL strategy was not used by Buffard et al. (2000) for the VFRoe solution of test case 2, maybe due to these numerical oscillations. For this test case, all three methods gave good results, with the MUSCL AUSM+ results being somewhat better. However, the HLFW method presented less oscillatory behavior.

In test case 3, for the HLFW method there was a discontinuity in the central part of the rarefaction wave (x » 0.5) due to entropy violation. It also showed large errors in the mixed wave region (x » 0.55-0.75). Numerical oscillations in the discontinuities similar to those observed in test case 2 also appeared in the MUSCL VFRoe solution, whereas there were no such oscillations for the MUSCL AUSM+ numerical solution, which, as can be seen in Figure 5, presented more numerical diffusion than MUSCL VFRoe.

The accuracies of all schemes were evaluated by comparison to the solution for a 2000-cell mesh, for which the MUSCL AUSM+ and MUSCL VFRoe methods obtain the same results. Table 3 contains the mean and maximum value of velocity, density and energy absolute deviations for the solutions of test case 3 by all three methods. The MUSCL VFRoe and AUSM+ schemes presented comparable accuracies for a given mesh, but the MUSCL AUSM+ results were slightly better but with more oscillations than the MUSCL VFRoe results. However, the MUSCL AUSM+ results presented more numerical diffusion with larger errors than the MUSCL VFRoe results in the shock wave region, as can be seen in Figure 5. These regions make the maximum deviations for the MUSCL AUSM+ larger than those for MUSCL VFRoe results, except for density. However, in a wide region near x = 0.5, the MUSCL VFRoe results presented larger residues than the MUSCL AUSM+ results. This region makes the values of mean deviations for velocity and energy larger for the MUSCL VFRoe solution. The HLFW scheme was less accurate than the MUSCL VFRoe and AUSM+ methods. Its results presented more numerical diffusion than the MUSCL VFRoe results but less oscillation than the MUSCL AUSM+ solution. However, the HLFW method presented large deviation in the regions of mixed waves and contact discontinuities.

In Figure 6 the fundamental derivative for test cases 1, 2 and 3 are presented. For test case 3, where the fundamental derivative presents a sign change result in mixed waves, the MUSCL AUSM+ and MUSCL VFRoe schemes were able to follow the density elevation near x = 0.7. However, the HLFW solution shows strong deviations in this region.

In order to compare the CPU time, test case 3 was considered. For refined meshes the methods tends to achieve the same accuracy. In order to compare the methods for a given accuracy, the MUSCL AUSM+ solution using a 200-cell mesh was considered and the meshes for the MUSCL VFRoe and HLFW methods were increased to 230 and 700 cells, respectively. The results for CPU times relative to that spent by MUSCL AUSM+ in the 200-cell mesh are shown in Table 3 for these meshes and for a 5000-cell mesh. For the given refined mesh, where there is no difference in the results, the HLFW method required 5.6 times more CPU time than the MUSCL AUSM+ method, but it was more than twice as fast as the MUSCL VFRoe method. However, for the meshes with a similar accuracy, MUSCL AUSM+ is tremendously faster than the other methods and MUSCL VFRoe became about nine times faster than HLFW. This kind of comparison is important because only small meshes are usually practical for a long time integration in multidimensional problems.

CONCLUSIONS

Some Riemann solver-based methods were compared for two ideal gas flow problems. Due to the difficulty in application of Riemann solvers to real gases, three Riemann solver-free methods were compared for three real gas flow problems taken from the literature. The methods were compared not only for a given mesh, but the performance for a given accuracy with different meshes was also analyzed. This methodology allows the correct evaluation of the advantages of the MUSCL strategy for achieving second-order accuracy.

Several issues were analyzed for the first time in the literature: the comparison of the recent proposed hybrid Lax-Friedrich-Lax-Wendroff (De Vuyst, 2004) method with the most popular hyperbolic solvers, the MUSCL second-order extension of the VFRoe method and the application of these methods to real gases flows. Besides, an accurate procedure for imposing the boundary condition when there are source terms was clearly presented and compared to common procedures in a pipeline shutting-in problem.

The Riemann solvers of Roe-Pike, Harten and the MUSCL-Hancock MINMOD and SUPERBEE second-order extensions of the Roe-Pike method were compared for accuracy and CPU time cost using ideal gas flows in shock tube and pipeline shutting-in test problems. The CPU cost was also evaluated in terms of the number of discretization cells, time steps and variables. The results show that, although the MUSCL-Hancock second-order extension requires more FLOP's for a given mesh and CFL number, it is faster than first-order methods because, for a given accuracy, it requires fewer cells and consequently allows larger time steps than the first-order methods. It was shown that the Roe-Pike solver with MUSCL-Hancock extensions are more accurate and faster than the original Roe-Pike and Harten methods. The SUPERBEE MUSCL scheme presented numerical oscillations near discontinuities common to second-order methods whereas the MINMOD MUSCL scheme did not.

The hybrid Lax-Friedrich-Lax-Wendroff, VFRoe and AUSM+, the latter two in their MUSCL-Hancock second-order extensions, were applied to real gases using the Van der Waals EOS. The accuracy and CPU cost of these methods were compared for three test cases. In one of the tests the MUSCL VFRoe and the hybrid method presented large errors due to entropy violation. The MUSCL AUSM+ and VFRoe schemes showed small oscillations near discontinuities for one of the test cases, also producing strange oscillations in two adjacent constant value regions. This was shown to be caused by the MUSCL second-order extension. The MUSCL AUSM+ scheme was always much faster than the others. The MUSCL VFRoe method was the slowest scheme for large meshes, but it was faster than the hybrid method for small meshes with comparable accuracies.

Therefore, the AUSM+ scheme is generally recommended for compressible gas-phase flow problems, although caution must be taken in the use of its MUSCL extension due to strange nonphysical oscillations for some real gas flows.

ACKNOWLEDGEMENTS

The authors acknowledge the financial support provided by PRH-ANP/MME/MCT – Agência Nacional do Petróleo – Brazil and CNPq (grant no. 351378/1994-4).

NOMENCLATURE

a parameter in Van der Waals EOS A_AUSM

auxiliary AUSM+ constant

Jacobian matrix b

parameter in Van der Waals EOS

B_AUSM

auxiliary AUSM+ constant

c

speed of sound

C_p

specific heat at constant pressure

C_v

specific heat at constant volume

D

diameter

e

specific internal energy

E

total energy per unit volume, E=r(0.5u²+e)

f

flux vector

F

numeric flux vector

g

nonhomogeneous vector in conservation equations

G

fundamental derivative

l

auxiliary function in the HLFW method given in equation (28)

h

specific enthalpy, h=e+p/r

auxiliary function in the Harten method

H total specific enthalpy, H=0.5²+h l

friction factor

m

mass flux

M

Mach number

nc

number of cells

nv

number of variables in u

ns

number of time steps

P

pressure

AUSM+ polynomials for pressure flux part

Q

viscosity function

R

right eigenvector of the Jacobian matrix

Â

universal gas constant

s

entropy

S

auxiliary variable in the Harten method in equation (11)

t

time

u

vector of unknown variables

v

velocity for one-dimensional flows

u _exact

local RP solution

U

vector of initially modified variables in the MUSCL reconstruction

vector of MUSCL reconstructed variables x

spatial coordinate

local spatial coordinate w

vector of primitive variables

Greek Symbols

a

wave strength

auxiliary variable in the Roe-Pike and Harten methods

g

ratio of specific heats, g = C_p/C_v

G

auxiliary variable in MUSCL data reconstruction

d

parameter in Van der Waals EOS, d =Â/C_v

D variation D auxiliary variable in MUSCL data reconstruction e

Dt/Dx ratio

z small parameter in viscosity function Q h

dissipation function

q

weighting parameter in hybrid method

l

Jacobian matrix eigenvalue

AUSM+ polynomials for mass flux part

x

small positive parameter in the HLFW method

P

thermodynamics derivative defined in equation (36)

r

density

auxiliary variable in the Roe-Pike and Harten methods j

mean state function

c

relaxation parameter of Lax-Wendroff flux

Y

auxilary vector in AUSM+ method

V specific volume W

auxiliary function in the Harten method

Subscripts

j

value at cell j

j+1/2

value at x_j+1/2 (cell face)

L

left state

R

right state

0

initial values

Superscripts

AUSM

related to the AUSM+ method

LW

related to the Lax-Wendroff method

LF

related to the Lax-Friedrich method

MLF

modified Lax-Friedrich method

n

value at time instant n

k

related to the k^th wave

*

variable modified by the AUSM+ method

±

related to different AUSM+ polynomials

Overscripts

~

evaluated at mean state

Abbreviations

AUSM

advection upstream splitting method

CFL

Courant-Friedrich-Lewy number, CFL = Dt (c+max(v))/Dx

CPU

central processing unit

EOS

equation of state

ENO

essentially non-oscillatory

FLOP

floating point operation

HLFW

hybrid Lax-Friedrich Lax-Wendroff method

IA

indirect adressing

MUSL

monotone upstream-centered scheme for conservative laws

RP

Riemann problems

RS

Riemann solver

TVD

total variation diminishing

TVNI

total variation nonincreasing

VFRoe

Volume of Fluid Roe method

Received: December 12, 2004

Accepted April 19, 2006

The boundary conditions are imposed at the inlet and outlet using the well-known ghost-cell approach, in which the first and last faces correspond to the domain limits and additional fictitious cells are outside the domain. The problem consists of how to determine the next time step values in the fictitious cells. In common shock problems the states in those cells are equal to those in the adjacent internal cells and there is no difficulty in imposing the boundary conditions. However, generally, as in pipeline flows, the fictitious state must be determined. For two-phase flows, Masella et al. (1998) pointed out that the commonly used extrapolation from the internal values to the unknown variables could result in inconsistent results, making it necessary to take into account the source terms and the differential equations themselves. Usually, the outlet pressure and inlet mass flow rate and temperature are kept constant. Figure A1 shows an example of those boundary conditions for a pipeline with friction losses with the extrapolation of the pressure from the last internal cell to the fictitious cell. The figure is a result of the propagation of waves in a 100-meter pipeline during its shutting-in and a zoom was applied in order to improve visualization. One can observe the change in curvature in the x = 99.5 m region.

The procedure used to impose the boundary conditions is based on the simultaneous solution of the discrete equations and the EOS using the known variable values at the boundaries. For the ideal gas isothermal one-dimensional flow with friction, the system of equations to be solved for the fictitious cell at the flow inlet with the specified mass flow rate is as follows:

where the index 0 corresponds to the fictitious cell, D is the diameter, l is the friction factor and c is the sound velocity which is constant for isothermal flow of an ideal gas. The solution of this system gives the following equation for the density and, consequently, the pressure, at the inlet fictitious cell.

Since the flow is isothermal and pressure is assumed to be constant at the outlet, the density is also constant at the outlet. Therefore, the discrete form of the momentum equation can be solved to obtain the outlet mass flow rate resulting in

where the index N corresponds to the fictitious cell at the outlet. For a general EOS, the procedure is similar, but the resulting system of equations must be solved numerically.

For the nonisothermal flow of an ideal gas it is assumed that the temperature is known at the inlet. In the general case, this implies the need for solving an additional relation involving enthalpy, pressure and temperature. However, for an ideal gas, the enthalpy is not a function of pressure and a constant temperature implies a constant enthalpy. Hence, the discrete form of the momentum equation and the EOS can be solved simultaneously. The system of equations to be solved is

where g'=g/(g-1). Solving the above system, the inlet density r₁ can be obtained from the following quadratic equation:

The temperature is not known at the outlet. Therefore, the discrete forms of momentum and energy equations are solved together with the EOS in order to evaluate the outlet total energy, density and mass flow rate. The final system is given by

Using the technique described above to impose the boundary conditions, the same flow example of the shutting-in of a 100-meter pipeline was solved and the results are shown in Figure A1. The improvement in the solution is quite clear.

• Abgrall, R., An extension of Roe's upwind scheme to algebraic equilibrium real gases, Computer and Fluids, 19, 171-182 (1991).
• Buffard, T., Gallouet, T. and Hérard, J. M., A sequel to a rough Godunov scheme: application to real gases, Computers and Fluids, 29, 813-847 (2000).
• Collela, P. and Glaz, H., Efficient solution algorithms for the Riemann problem for real gases, J. of Computational Physics, 59, 264-289 (1985).
• Cox, C. F. and Cinnella, P., General solution procedure for flows in local chemical equilibrium, AIAA J., 32, 519 (1994).
• De Vuyst, F., Stable and accurate hybrid finite volume methods based on pure convexity arguments for hyperbolic systems of conservation law, J. of Computational Physics, 193, 426-468 (2004).
• Godunov, S. K., A difference method for numerical calculation of discontinuous equation of hydrodynamics, Math Sbornik (in Russian), 47, 217-306 (1959).
• Glaister, P., An approximate linearized Riemann solver for the Euler equations for real gases, J. of Computational Physics, 74, 382-408 (1988).
• Grossman, B. and Walters, R. W., Analysis of the flux-split algorithms for Euler's equations with real gases, AIAA J., 27, 524 (1989).
• Guardone, A. and Vigevano, L., Roe linearization for the Van der Waals gas, J. of Computational Physics, 175, 50-78 (2002).
• Harten, A., High resolution schemes for hyperbolic conservation laws, J. of Computational Physics, 49(3), 357-393 (1983).
• Harten, A. and Hyman, J. M., Self adjusting grid methods for one-dimensional hyperbolic conservation laws, J. of Computational Physics, 50, 235-269 (1983).
• Harten, A., Enquist, B., Osher, S. and Chakravarthy, S., Uniformly high order essentially non-oscillatory schemes III, J. of Computational Physics, 71, 231-303 (1987).
• Liou, M. S. and Steffen, C. J., A new flux splitting scheme, J. of Computational Physics, 107, 23-39 (1993).
• Liou, M. S., A sequel to AUSM: AUSM+, J. of Computational Physics, 129, 364-382 (1996).
• Liou, M. S. and Edwards, J. R., AUSM schemes and extensions for low mach and multiphase flow, VKI LS 1999-03, Computational Fluid Dynamics (1999).
• Liou, M. S., Van Leer, B. and Shuen, J. S., Splitting of inviscid fluxes for real gases, J. of Computational Physics, 87, 1-24 (1990).
• Masella, J. M., Faille I. and Gallouet T, On an approximate Godunov scheme, International Journal of Computational Fluid Dynamics, 12 (2), 133-149 (1999).
• Masella, J. M., Tran Q. H., Ferre, D. and Pauchon, C., Transient simulation of two-phase flows in pipes, International Journal of Multiphase Flow, 24(5), 739-755 1998.
• Menikoff, R. and Plhor, B. J., The Riemann problem for fluid flow of real materials, Reviews of Modern Physics, 61, No. 1, 75-131 (1989).
• Mottura, L., Vigevano, L. and Zaccanti, M., An evaluation of Roe scheme generalization for equilibrium real gas flow, J. of Computational Physics, 138, 354-399 (1997).
• Quartapelle, L., Castelletti, A., Guardone, A. and Quaranta, G., Solution of the Riemann problem of classical gasdynamics, J. of Computational Physics, 190, 118-140 (2003).
• Roe, P. L., Approximate Riemann solvers, parameter vector, and difference schemes, J. of Computational Physics, 43, 357-372 (1981).
• Saurel, R., Larini, M. and Loraud, J. C., Exact and approximate Riemann solvers for real gases, J. of Computational Physics, 112, 126-137 (1994).
• Toro, E. F., Riemann Solvers and Numerical Methods for Fluid Dynamics - A practical introduction, 2^nd edition, Springer, Berlin (1999).
• Toumi, I., A weak formulation for Roe's approximate Riemann solver, J. of Computational Physics, 102(2), 360-373 (1992).
• Van Leer, B., Towards the ultimate conservative difference scheme IV. A new approach to numerical convection, J. of Computational Physics, 23, 276-299 (1977).
• Van Leer, B., Towards the ultimate conservative difference scheme V. A second-order sequel Godunov method, J. of Computational Physics, 32, 101-136 (1979).
• Van Leer, B., On the relation between the upwind-differencing schemes of Godunov, Enguist-Osher and Roe, SIAM J. Sci. Stat. Comput., 5(1), 1-20 (1985).
• Vinokur, M. and Montagné, J. L., Generalized flux-vector splitting and Roe average for an equilibrium real gas, J. of Computational Physics, 89, 276-300 (1990).

Appendix

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  14 Dec 2006
• Date of issue
  Sept 2006

History

• Accepted
  19 Apr 2006
• Received
  12 Dec 2004