A comparison of hyperbolic solvers II: ausm-type and Hybrid Lax-Wendroff-Lax-Friedrichs methods for two-phase flows

PROCESS SYSTEMS ENGINEERING

A comparison of hyperbolic solvers II: ausm-type and Hybrid Lax-Wendroff-Lax-Friedrichs methods for two-phase flows

R. M. L. Coelho^I,*; P. L. C. Lage^II; A. Silva Telles^I

^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

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

ABSTRACT

Riemann-solver based schemes are difficult and sometimes impossible to be applied for complex flows due to the required average state. Other methods that do not use Riemann-solvers are best suited for such cases. Among them, AUSM+, AUSMDV and the recently proposed Hybrid Lax-Friedrichs-Lax-Wendroff (HLFW) have been extended to two-phase flows. The eigenstructure of the two-fluid model is complex due to the phase interactions, leading to numerous numerical difficulties. One of them is the well-posedness of the equation system because it may lose hyperbolicity. Therefore, the methods that are not based on the wave structure and that are not TVNI could lead to strong oscillations. The common strategy to handle this problem is the adoption of a pressure correction due to interfacial effects. In this work, this procedure was applied to HLFW and AUSM-type methods and their results analyzed. The AUSM+ and AUSMDV were extended to achieve second-order using the MUSCL strategy for which a conservative and a non-conservative formulation were tested. Additionally, several AUSMDV weighting functions were compared. The first and second-order AUSM-type and HLFW methods were compared for the solution of the water faucet and the shock tube benchmark problems. The pressure correction strategy was efficient to ensure hyperbolicity, but numerical diffusion increased. The MUSCL AUSMDV and HLFW methods with pressure correction strategy were, on average, the best of the analyzed methods for these test problems. The HLFW was more stable than the other methods when the pressure correction was considered.

Keywords: Hyperbolic conservation laws; AUSM+; AUMSDV; Hybrid schemes; Two-phase flow; Two-fluid model.

INTRODUCTION

The simulation of two-phase flows is important in many industrial applications and it is governed by the conservation of mass, momentum and energy and the associated constitutive equations. The later are macroscopic relations, such as an equation of state (EOS) and interfacial flux equations, which represent the behavior of the associated microscopic phenomena. These relations strongly influence the structure and dynamics of the waves in the flow. Molecular dynamic models have been locally applied to model the interfaces in two-phase flows. However, it is not feasible to apply it to the whole flow domain. Therefore, average Navier-Stokes models have been used in which the singular interface with its appropriate jump conditions are averaged. This approach has been proposed in Ishii (1975) in the classical two-fluid model. The two-fluid model is based on field averaging and accounts for both mechanical and thermal non-equilibrium, being able to deal with a wide range of flow patterns. The two-fluid model includes the equations for mass, momentum and energy conservation for each phase. It is possible to sum the phase momentum equations to obtain a single equation for the whole mixture which requires additional relations for the phase slip velocity. This is the idea behind the drift-flux model. It is possible to further simplify the model by eliminating the acceleration terms (no pressure wave model Masella et al., 1998), in which the flow is driven by the mass transport.

Frequently, the validation limits of the EOS's used are violated at the interfaces and the variables evaluated from them are inaccurate or physically unrealistic (negative pressures, densities and volume fractions). Besides, in the case of multicomponent mixtures, the notion of partial pressure is valid only with the local thermodynamic equilibrium assumption, which is acceptable when all components are perfectly mixed and the number of molecular collisions is so large that temperature is rapidly homogenized. When two fluids are not perfectly mixed, there is an interface which is the locus of the collisions between their molecules. Therefore, at an interface the local thermodynamic equilibrium does not hold anymore. It is expected that the interface propagates with an intermediate velocity and has an intermediate pressure which shall be modeled. This implies the presence of non-conservative terms in the two-fluid model which must be carefully modeled to achieve a stable hyperbolic model. Otherwise, the two-fluid model becomes ill-posed, resulting in numerical instabilities during its solution.

The conservation equations for advection-dominated two-phase flow problems can be written as a set of conservation equations whose solution for the flow variables may present large gradients or even discontinuities. Most of the existing commercial codes, such as OLGA (Bendiksen et al., 1991), are based on very diffusive simple upwind pressure-based solvers which may lead to unrealistic solutions due to their failure in solving the flow wave field. These methods have limited accuracy compared to density-based solvers when strong compressibility is present. 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). These methods are based on characteristic upwind differencing and have been applied from nuclear and oil industries to solid combustion. 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 pioneer method was that of Godunov (Godunov, 1959) which is based on the exact solution of the Riemann problem at the interfaces. Its advantage is the assurance of the density, pressure and volume fraction positivity, but its drawbacks are its large computer cost and the lack of generality due to its requirement of an analytical integration of the Riemann invariants.

Flux difference splitting (FDS) schemes are modifications of classical upwind schemes. Such schemes are accurate for all speed flows and two-phase flows but time consuming due to the repeated evaluation of the Jacobian matrix of the flux with respect to the conservative variables and its eigenstructure and they do not guarantee positivity of density and pressure. These schemes are approximate Riemann solvers and the most popular one is due to Roe (Roe, 1981). The generalization of the FDS scheme is difficult because it is not always possible to analytically obtain the necessary expressions for the average state in which the Jacobian matrix must be evaluated. This class of methods was adopted by several authors (Sainsaulieu, 1995; Toumi, 1996, Toumi and Kumbaro 1996; Tiselj and Petelin, 1997; Fjelde and Karlsen, 2002; Cortes et al., 1998; Romate, 1998; Faille, 1999; Saurel and Abgrall, 1999).

The flux vector splitting schemes (FVS) are simplifications of the FDS schemes and are based on scalar rather than matrix calculations. Therefore, these methods are less time consuming but are more diffusive. Recently, great effort has been spent on developing hybrid schemes in order to achieve FVS efficiency with FDS accuracy. These methods are well-known as Advection Upstream Splitting Methods (AUSM). The AUSM scheme (Liou and Steffen, 1993) uses the idea of splitting the flux into a convective part upwinded in the flow direction and a pressure part upwinded according to the acoustic properties of the flow. The flow direction is determined by the Mach number sign considering information for the left and right state for a given cell face. This approach allows the capture of stationary contact discontinuity with reduced numerical diffusion even with first-order discretization. However, for oblique shocks (not aligned with the mesh), the AUSM method admits highly non-monotone solutions for pressure and density, leading to difficulties in convergence. Thereafter, Wada and Liou (1997) proposed a method called AUSMDV, which mitigates the occurrence of non-monotonic solutions allowing robust resolution of shock contact (steady and moving) waves but introduced the carbuncle effect. Recently, Liou (1996) proposed a modified scheme called AUSM+, which is capable of solving contact discontinuities exactly without non-monotonic problems. Besides, the method preserves the positivity of pressure, density and volume fraction.

The application of the density-based solvers to the two-fluid model is not simple due to the non-hyperbolic nature of the two-fluid model, the presence of non-conservative terms, stiff source terms, complex interfacial relations, phase appearance and disappearance and phase change. These methods also experience loss of accuracy for low Mach number flows. In such cases, AUSM+ behaves more like a central difference discretization scheme. This may lead to odd-even decoupling and it is then necessary to couple velocity and pressure (Liou, 1999). The larger the difference between the speed of sound and the velocity, the stronger the odd-even decoupling effect is. Some strategies have been proposed for dealing with one and two-phase low Mach number flows in the AUSM+ context (Liou, 1999 and Paillère et al., 2003). AUSM+ has been used in many works for two-phase drift flux (Liou, 1999, Paillère et al., 1999, Fjelde, 2002, Evje and Fjelde, 2003) and two-fluid (Niu, 2001, Paillère et al., 2003, Evje and flatten, 2003, De Vuyst, 2004) models.

It is well-known that first-order schemes tend to present numerical diffusion, while the second-order schemes tend to present oscillations near high gradient regions. De Vuyst (2004) presented a novel hybrid method and applied it to some one-phase flow problems and to the water faucet two-phase flow problem. The idea of this method (HLFW) is to evaluate the flux with a weighted combination of second-order Lax-Wendroff and first-order Lax-Friedrichs schemes. The weighting parameter is the key to the success of this type of methods. The author proposed a concept of numerical dissipation by convexity in which the entropy and its flux are not directly used for the dissipation evaluation. This allows the dissipation to be expressed as a simple convex function, but also allows some entropy violation. The proposed dissipation term is evaluated in the whole cell instead of the common approach where the dissipation is evaluated at the interface.

Schemes like hybrid Lax-Friedrichs Lax-Wendroff (De Vuyst , 2004), AUSM (Liou, and Steffen, 1993), AUSM+ (Liou, 1996) and AUSMDV (Wada and Liou, 1997, Liou, 200) are not based on the characteristic analysis of the flow and do not require the model to be hyperbolic. However, as shall be shown in this work, their accuracy is improved if the model is hyperbolic.

Another class of methods not studied here are the relaxation schemes (Evje and Fjelde, 2002), which are capable of dealing with strong non-linearities. However, they require a case-by-case analysis of the sub-characteristic stability based on entropy compatibility or Chapman-Enskog-like expansion (De Vuyst, 2004).

Many other authors adopted different methods in order to solve two-phase flows. Recently, Lorentzen and Fjelde (2005) applied a slope limiter to the finite element predictor-corrector shooting technique for the drift-flux model. Wackers and Koren (2005) applied a HLL Riemann Solver to a new model with five equations. Cascales and Salvador (2006) applied a high resolution method and a slope limiter to the two-fluid model. Guinot (2001a, b) applied the Godunov method to solve a simplified two-phase flow model. Castro and Toro (2006) compared a Riemann Solver and a second-order extension using MUSCL and ADER strategies to the Saurel and Abgrall (1999) five equation model. Bedjaoui and Sainsaulieu (1997) examined the effects of diffusion and relaxation parameters on the stability of the solution of a non-hyperbolic two-fluid model, proving its stability depending upon the initial conditions.

As the ideal EOS is not capable of representing non-classical waves, many numerical methods were recently extended to real materials. In our previous work (Coelho at al., 2006), a comparison was made regarding the accuracy and CPU time consumption of second-order extensions of Roe's Riemann solver, VFRoe, AUSM+ and Hybrid Lax-Friedrichs Lax-Wendroff for ideal and real gases modeled with cubic equations of state, including non-classical wave phenomena.

In this work, the AUSM-type and Hybrid Lax-Wendroff-Lax-Friedrichs methods were compared for two-phase flow simulations of benchmark problems. A stiffened gas EOS was adopted in order to make possible comparison with some results presented in the literature. This paper is organized as follows: section 2 is a summary of the two-fluid model equations. In section 3, a summary of the AUSM+, AUSMDV and Hybrid Lax-Friedrichs-Lax-Wendroff (HLFW) methods and the MUSCL (Monotone Upstream-Centered Scheme for Conservation Laws) strategy due to Van Leer (1977 and 1979) are presented. In section 4, the discretization schemes of additional terms in the two-fluid model are described. In section 5, the results of the application of these methods to two-phase flows are discussed. The methods were compared in terms of accuracy for different mesh sizes. The usage of the pressure correction term that enforces hyperbolicity and of different expressions for the weighting function of the AUSMDV scheme was also analyzed. Additionally, MUSCL strategies based on primitive and conservative variables were compared.

THE TWO-FLUID MODEL

The set of hyperbolic one-dimensional conservation laws is generically written as:

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

where is the Jacobian matrix. The Riemann solvers require the calculation of the eigenvalues and eigenvectors of the Jacobian matrix, which is evaluated with an average state. This average state is difficult and sometimes impossible to be obtained.

The system of equations of the two-fluid model (Ishii, 1975) applied to two-phase flow and with no source/sink terms could be written in the form of Eq. (1) as follows:

where ρ is the density, v is the velocity, is the interface pressure, which represents the effects of surface tension and low wavelength dynamics, α is the volume fraction, represents the interfacial and wall drag forces, g_x is the gravity acceleration component along the x-coordinate and subscripts g and l refer to gas and liquid phases, respectively.

Defining Δ= p_k - the following transformation could be applied:

leading to the following form:

In the present work, the pressures and the interfacial pressures were assumed to be equal for both phases, that is, p = p_g = p_l and pⁱ = = . The authors in the references (Toumi, 1996, Cortes et al., 1998, Evje and Fjelde, 2003) used a perturbation method in order to obtain the approximate eigenvalues of the system. It is important to highlight that, if Δp = 0, the model becomes non-hyperbolic due to the existence of complex eigenvalues.

Let w = [α_g, v_g, v_l, p]^t be the vector of non-conservative variables. It is necessary to obtain the transformation of u into w and vice-versa. The following equation is a relation between u and the phase velocities:

Besides,

These equations lead to the following non-linear relation, which can be solved to give the pressure and, thus, the densities:

The transformation of w into u is trivial. The complete description of the system is accomplished by specifying an EOS. For an isentropic evolution, the following relation was assumed:

where γ is the ratio of the specific heats and A_k and Π_k are constants depending upon the fluid. This EOS reproduces the stiffened gas EOS, p_k = (ρ_k/c)- (/c), if A_k = (1/c) and .

AUSM+, AUSMDV AND THE HYBRID LAX-FRIEDRICHS-LAX-WENDROFF METHODS

We are interested in solving the Finite Volume discretization of Eq. (1) with the following form:

where ε = Δt/Δx, F_j+1/2 and F_j–1/2 are the numerical flux calculated at the interface j-1/2 and j+1/2 of the volume j. For the Euler explicit the numerical flux is evaluated at the t time level and for the Euler implicit it is evaluated at the t+Δt time level.

AUSM Type Methods

Liou and Steffen (1993) proposed a scheme based on the separation of a convective flux and a pressure flux, presenting some drawbacks, especially in the case of strong shocks. Wada and Liou (1997) presented an improved scheme called AUSMDV, which is a blend of the FDS and FVS methods. However, AUSMDV does not have the property of capturing a stationary shock. Thereafter, Liou (1996) presented a novel scheme called AUSM+, which is capable of capturing stationary shock and contact discontinuity waves, preserves the positivity of pressure and density and improves the accuracy of the AUSM method.

The AUSM+ is as accurate as the flux splitting methods due to Roe (1981) without the wave decomposition, leading to a large reduction of CPU time consumption and making the formulation simple and quite general. The AUSM+ is able to solve the contact discontinuities exactly, making it suitable for viscous flow calculations.

The AUSM+ flux for the system given by eqs. (1) and (3) is given by the following expression, where the pressure and mass flux parts are separated:

where the subscript k indicates the k^th phase (gas, g, or liquid, l) and the following upwind formula based on the sign of _k^j+1/2 is used to compute (Paillère et al., 2003):

using the value obtained by the following steps:

(i) Evaluation of left and right Mach numbers calculated by 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:

where the superscripts "+" and "–" are associated with the right and left running waves, respectively.

(iv) finally, the mass flux is given as follows:

The common speed of sound defined in Eq. (12) is the key to the AUSM+ accuracy instead of the simple upwind (

^+1/2 = , for > 0, ^+1/2 = ⁺¹, otherwise). Liou (1996) presented a correct expression for the interface sound speed based on the isoenergetic condition. This involves an EOS dependent expression. In this paper this was avoid by using Eq. (12), as suggested by the author, which does not guarantee the exact capture of stationary shocks. Besides, the first RHS term of Eq. (11) is a weighted Mach number averaged flux. The coefficient is merely a scalar, reducing CPU demand for the flux evaluation and keeping the algorithm free of jacobian evaluations. The original AUSM+ polynomials are:

These polynomials were developed in order to satisfy some properties (see Liou 1996, for details). They are responsible for the more general characteristic of the AUSM+ scheme. Liou (1996) adopted the criterion in which µ^±(M) has an inflection point at M=0 and ^±(M) has two inflection points for M = ±1, leading to A = 3/16 and B = 1/8. These parameter values could also be changed in order to improve AUSM+ results, but they were left fixed as suggested in Liou (1996).

The AUSM+ suffers of odd-even decoupling for low Mach numbers, losing accuracy. This is particularly important for liquid flows where the sound speed is usually three orders of magnitude higher than the velocity. As cited in the introduction, density-based solvers exhibit a stiffness problem for low Mach number where it is necessary to couple pressure and density. For this purpose, a cut-off procedure was adopted by Paillère et al. (2003) in order to limit the sound speed and consequently the Mach number. This procedure was applied to the liquid phase in the present work and is explained in the following.

1) Modify Eq. (16) to add a pressure diffusion term for the liquid phase given by:

2) Redefine the interface numerical speed of sound, Eq. (13), to:

where M₀ is a cut-off Mach number, say 10^-4.

3) compute

4) evaluate pressure diffusion term as:

This last term adds the necessary dissipation scaled for low Mach numbers.

Before the AUSM+ formulation, the authors in Wada and Liou (1997) presented an improvement to AUSM creating the AUSMDV scheme. This method presented some problems with the carbuncle effect which were solved in AUSM+. However, for one-dimensional problems, we believe it remains desirable and the results confirmed this expectation. The scheme is obtained by the substitution of Eq. (16) with the following expression:

where:

Additionally:

Wada and Liou (1997) adopted Y = p/ρ . Thereafter, Liou (2000) suggested Y = 1/ρ in order to eliminate shock instability (carbuncle phenomenon) and Evje and Flatten (2003) adopted Y= ρ/α . One should bear in mind that the shock instability is not important in the one-dimensional problems analyzed in this paper. In the analysis of the results, the difference between these options will be further discussed.

Hybrid Lax-Wendroff-Lax Friedrichs Method

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). The author used the Lax-Wendroff and the Lax-Friedrichs fluxes. The Lax-Wendroff flux with a relaxation parameter χ for each phase k can be written as:

where = f(u_k,j), = f(u_k,j+1) and is a consistent mean state function, that is, it has the property that (u,u) = u. In the present work, a constant value of χ = -10^-2 and (u,v) = 0.5(u + v) were used as suggested by the author. 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 local dissipation function, η , which is defined as:

which, in its turn, depends on the numerical flux given by Eq. (28) and on the final updated field, . The method was implemented by the algorithm suggested in De Vuyst (2004) to achieve η < 0, which calculates as the smallest values that make negative. Let Δ > 0 be a small increment and build the sequences as follows:

let q = 0 and = 0 and

while

> ξ and < 1 do

q = q + 1

set = and calculate

= max (,)

evaluate the numerical flux, Eq.(29)

calculate the update field, Eq. (10)

evaluate

end while.

where ξ is a small positive number. A Newton-Raphson like algorithm could reduce CPU demand. In Eq. (29), S is the convexity function and G its derivative.

Extension of AUSM-Type Methods to Second-Order

The extension of the present methods to second-order follows Van Leer (1977, 1979) and the MUSCL method modified later by Van Leer (1985) (MUSCL-Hancock, see Toro, 1999, for details). 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:

For an approximate Riemann solver, the last step is the solution of the Riemann problem at the "j+1/2" face with the values e as left and right values, respectively. For the AUSM type methods the values of e substitute the values of U_j and U_j+1 in the flux evaluation. The Δ_j values are obtained as follows:

where Γ = 1 and 2 leads to the MINMOD (MINBEE) and SUPERBEE schemes, respectively.

The MUSCL data reconstruction can be applied to primitive variables instead of conservative variables, but not when using the previously shown algorithm, because Eq. (32) requires that the dimension of the reconstructed data vector be equal to the flux vector dimension which, in turn, is equal to the number of conservative variables. However, in some tests, the MUSCL-Hancock applied to AUSM+ presented too much pressure oscillation. Therefore, a primitive variable based data reconstruction algorithm was also tested. A primitive-variable version of the MUSCL-Hancock algorithm can be found in Toro (1999), but it requires the Jacobian matrix evaluation. Therefore, a slope limiter scheme was adopted. It consists of Eq. (31) with W (the vector of primitive variables) in substitution to U. Finally, the following expression for instead of Eq. (33):

where Δ_j–1/2 and Δ_j+1/2 have the same definition as in Eq. (33), Θ(r) is a slope limiter and ϖ ∈[-1,1]. The MINMOD (MINBEE) reconstruction is obtained if the slope limiter is taken as:

DISCRETIZATION OF THE REMAINING TERMS

Interfacial Transfer Terms

The two-fluid model is non-hyperbolic for = p_k. In order to regularize this term, some interfacial pressure models were created. A summary of these models can be found in (Cortes et al., 1998). The following form for the interfacial pressure was adopted in (Bestion, 1990, Cortes et al., 1998, Evje and Flatten, 2003, Paillère et al., 2003) (with = = pⁱ and p_l = p_g = p):

where σ is a positive constant. For σ> 1 the two-fluid model becomes hyperbolic, while for σ< 1, the eigenvalues of the flux jacobian matrix flux become complexes. One should note that the value of σ = 1 (used in CATHARE code) is a limit based on an approximate eigenvalues analysis and values of σ> 1.2 are recommended to achieve hyperbolicity. The effect of this parameter on the contact wave is important and will be presented in the results. The authors in Paillère et al. (2003) presented results for σ = 2 and σ = 5, but the adoption of these values is questionable. As equation (36) is more a mathematical artifact to keep the model hyperbolicity than a physical model, simulation results with \sigma=0 are also included to show how much its usage, which implies in additional numerical diffusion, affects the results. Besides, it is not clear if the methods treated here could keep useful results for reasonable refined meshes even if the hyperbolicity is lost since they do not make directly use of the hyperbolicity.

The interfacial drag force depends upon the flow pattern, which is, in turn, dependent upon the superficial velocities of the phases. One possible form for this term is:

where C_D is a positive drag coefficient that depends upon the phase and interfacial properties and the flow pattern, is the specific area and is the wall drag force. As for all interphase change fluxes, we have:

Phase Appearance and Disappearance

As the formation of vacuum is a challenge to single phase flow simulation, the phase appearance and disappearance are important numerical difficulties in multiphase flow simulation. First, the numerical dissipation of the scheme must be evaluated for conditions where the volume fraction tends to zero or one. The AUSM+ flux, Eq. (11), is still non-singular since the Mach number is bounded. The corresponding velocity bounds imply that the velocity of the disappearing phase must tend to the remaining phase velocity as it disappears. For example, if the gas volume fraction tends to 0, the following expression is used:

where (α) is a positive function equal to 1 if α > α_min, with α_min very close to zero and tending smoothly to zero if α <α_min.

Source Terms Treatment

Following Paillère et al. (2003), the source terms were discretized with a simple central formula as follows:

The remaining gravity and drag terms were evaluated at cell j. One could observe that all discretizations were based on the two-fluid model form given in Eq. (3). Some authors adopted the form given in Eq. (5) where the flux presents Δp instead of p and the non-conservative term includes the pressure gradient instead of the volume fraction gradient. Evje and Flatten (2003) used this form but treated the non-conservative term as part of the AUSM flux instead of considering it a source term. Cascales and Salvador (2006) adopted both central and backward discretization for this term.

For the Hybrid Lax-Wendroff-Lax Friedrichs method, the treatment of the non-conservative terms is included in the flux. The modified Lax-Wendroff flux given by Eq. (27) is substituted by:

and the variables update is also changed to:

where B is the non-conservative term and B_j+1/2=B[(U_j+U_j+1)/2].

NUMERICAL RESULTS

In our comparison study, the accuracy of the methods was compared for different mesh sizes. The shock tube and the water faucet test problems were employed. Particularly, the effects of the pressure correction term and of the choice of the AUSMDV weighting term were examined. Some complex cases such as phase separation were not included in this analysis because a special treatment for the model singularity at each limit of the phase fraction value would have to be included and analyzed.

Shock Tube Problem

We consider a shock tube problem with large relative velocity proposed in Cortes et al. (1998), which is a Riemann problem with the left and right initial states given by w_L=[0.71, 65, 1, 265000] and w_R=[0.7, 50, 1, 265000], respectively. The tube length is 100. Since there is no analytical solution available, the results from the HLFW method using a 10000-cell mesh were adopted as reference values. The MUSCL AUSM+ was adopted in Fjelde (2002) for the drift flux model, but their MUSCL reconstruction used non-conservative instead of conservative variables, as frequently adopted. Therefore, the MUSCL-Hancock strategy was substituted by a primitive (non-conservative) variable algorithm described previously. However, as shown in Figure 1 (a) and (b), the results are slightly worse than the ones obtained using the conservative variable approach. Similar results were observed for the water faucet problem. Therefore, only results for the conservative variable MUSCL scheme are presented in the following. Figure 2 (a) and (b) show the results for different choices of the weighting function and σ for AUSMDV scheme. It can be seen that the use of σ =1.2 led to better results since it produces less oscillation. Figure 2 (c) and (d) show the results for different choices of the weighting function and σ =1.2 for the AUSMDV and MUSCL-AUSMDV for 1000 and 10000 cell meshes. Liou (2000) suggested Y = 1/ρ and one can see that this option gives similar results when compared to Y = p/ρ . For refined meshes, only the Evje and Flatten (2003) suggestion, Y= ρ/α , could represent the two separate waves pattern. However, the Y = ρ/α option applied to AUSMDV produced some oscillations. This effect is even worse for MUSCL-AUSMDV. The options Y = p/ρ and Y = 1/ρ produce too much diffusion for AUSMDV. Since this is a shock tube problem, the less diffusive option is usually considered to be more suitable. Therefore, for this example, the option Y=ρ /α seems to be more adequate for AUSMDV, despite some oscillatory behavior. For the MUSCL-AUSMDV scheme, Y=ρ /α produces unacceptably large oscillations. Therefore, although Y=1/ρ produces too much numerical diffusion, it was adopted as the option for MUSCL-AUSMDV, which was capable of capturing the two separate waves only for the 10000 cell mesh.

Figures 3 to 6 show results for localized regions that have been magnified to highlight the differences among the schemes. Comparing these results, the following remarks can be made:

(i) The AUSMDV scheme was more diffusive than the corresponding AUSM+ scheme and the MUSCL AUSMDV is more diffusive than MUSCL AUSM+ with or without pressure correction. For σ = 0, all methods presented unstable solutions and no solution for a 10000 cell mesh was possible.

(ii) In the results presented by Evje and Flatten (2003) with a Roe Riemann solver, two separate void fraction waves can be seen. In our results, this scenario is observed for meshes up to 1000 cells for all σ. However, the two separate waves are accompanied by overshoots for all AUSM type methods with some advantage for AUSMDV. For the HLFW methods the overshoot was observed for σ = 0, but not for σ up to 1.2. Therefore, the convex function proposed in De Vuyst (2004) was efficient for this case.

(iii) Although the usage of σ = 1.2 reduced the numerical oscillations for both pressure and void fraction, the amplitude of the two separate waves remains slightly higher than the Riemann solver solution used as reference by Evje and Flatten (2003).

(iv) The MUSCL-Hancock strategy presented the expected behavior, favoring numerical oscillations for meshes up to 200 cells. Therefore, the extension to second-order of such methods by MUSCL was suitable only for small meshes for the adopted Y expression.

(v) From the results, the most suitable methods were the MUSCL AUSMDV and HLFW schemes with the pressure correction ( σ = 1.2). The accuracy of HLFW was superior to that of MUSCL AUSMDV for all mesh sizes. The MUSCL AUSMDV presented small oscillations only for the 10000-cell mesh. This result is particularly interesting because this behavior is different from the results for one-phase flows obtained in our previous work (Coelho et al., 2006), in which MUSCL AUSMDV was slightly better than HLFW. The MUSCL AUSM+ presented results with strong oscillatory behavior for all mesh sizes and it is not recommended.

(vi) The authors in Paillère et al. (2003) observed that the value of σ strongly affects a contact discontinuity. In Figure 2 (a), it is possible to observe the effect of σ =2. For all methods and weighting parameter functions, the contact discontinuity wave showed the expected positive pressure gradient at x=0.5 for σ =0 while the results for σ =2 showed a negative pressure slope at this point. Therefore, values of σ higher than 1.2 should be considered with care.

It should be pointed out that the AUSMDV less diffusive weighting function (Y = 1/ρ ) was adopted to generate these results, aiming to represent the two separate void fraction waves. However, it actually added considerable numerical oscillation. The other weighting functions were not tried because they are expected to produce more numerical diffusion and the weighting function used already gave results with larger numerical diffusion than those obtained with HLFW.

Water Faucet

The water faucet problem (Ransom, 1987) became a benchmark and has been extensively used by many authors in order to test numerical schemes in tracking volume fraction fronts (Toumi, 1996, Coquel et al., 1997, Paillère et al., 2003, Evje and Flatten, 2003, De Vuyst, 2004). This problem is a study of the gravity effects on a vertical downward water jet flow. Consider a vertical pipe of length 12 m initially at a uniform state with p = 10⁵ Pa, α_l = 0.80, v_g = 0 and v_l = 10 m/s. The pressure at the outlet is kept constant at the initial value and, at the inlet, the values of the remaining variables are kept constant and equal to their initial values. A scheme of the problem is shown in Figure 7. This problem has the following analytical solution:

The problem was also solved for different choices of the weighing function, Y, in the AUSMDV scheme. The results are presented in Figure 8. The adoption of Y=α/ρ instead of Y=p/ρ or Y=1/ρ produced no improvement in AUSMDV or in MUSCL-AUSMDV when replacing Y= 1/ρ . Liou (2000) has observed that Y= 1/ρ eliminates the carbuncle phenomenon, but this is not important for one-dimensional problems. In our tests, no appreciable differences were observed between Y= 1/ρ and Y=p/ρ for AUSMDV for the water faucet problem. However, Y=p/ρ was unstable for MUSCL-AUSMDV and the usage of Y=1/ρ produced much better results then the other options. Particularly, for refined meshes, the usage of Y=1/ρ showed results without the strong oscillations and with an accuracy equivalent of the HLFW method. This is due to the more diffusive characteristic of the scheme imparted by Y=1/ρ , which was prohibitive for the shock tube problem. Therefore, Y= p/ρ and Y= 1/ρ were considered to be the better choices for the water faucet problem for the AUSMDV and MUSCL-AUSMDV methods, respectively.

Figures 9 (a) to (c) shows the comparison between the solutions obtained for the water faucet problem for 320, 640 and 1280-cell meshes, using AUSM+, AUSMDV (Y=p/ρ ), these two methods with MUSCL-Hancock data reconstruction (Y=1/ρ ) and the HLFW method without interfacial pressure correction, that is, Δp = 0. Figure 10 (a) to (c) shows the same solution but with σ = 1.2 in Eq. (36) for Δp calculation. The following behavior can be observed from Figure 9:

(i) The AUSMDV scheme was more diffusive than the corresponding AUSM+ scheme and the MUSCL AUSMDV is more diffusive than MUSCL AUSM+ with or without pressure correction.

(ii) Since the methods are not TVD, for σ =0 the oscillations tended to grow when the mesh was refined or when the order is extended and this was confirmed for both MUSCL and HLFW schemes. AUSM+ presented some oscillations for the 320 cell-mesh, which increased when the mesh was refined and AUSMDV presented strong numerical diffusion. The MUSCL second-order extended schemes presented strong oscillations for the 320-cell mesh, while the HLFW scheme presented some oscillation for this mesh and both could not converge for the 1280-cell mesh. However, HLFW presented an oscillatory solution for the 640-cell mesh, while the MUSCL-AUSM+ schemes could not be solved for both 640 and 1280-cell meshes. MUSCL-AUSMDV also presented strong oscillations for the 640-cell mesh.

From Figure 10, it can be seen that:

(iii) The pressure correction term increased the accuracy of all schemes. It made the solution with MUSCL schemes possible for all mesh sizes.

(iv) As observed in the shock-tube problem, the accuracy of HLFW was the best among the methods analyzed.

CONCLUSIONS

The AUSM+, AUSMDV, Hybrid Lax-Friedrichs-Lax-Wendroff, and MUSCL-Hancock second-order extensions of AUSM+ and AUSMDV methods were compared in terms of accuracy using two benchmark problems. Three AUSMDV weighting functions were compared and the most suitable one was found to be not only problem dependent, but also order dependent, requiring different expressions for the AUSMDV and the extended second-order MUSCL-AUSMDV schemes. Primitive and conservative MUSCL-Hancock extensions were compared and the conservative variable option was slightly better. All schemes were compared with and without the interfacial pressure correction term.

For the water faucet test problem, the pressure correction term improved considerably the solution. The results without the pressure correction term showed that the HLFW method was stable only for coarse meshes. The same observation is valid for the MUSCL second-order extensions of the AUSM-type methods. AUSMDV was stable for large meshes, but showed more numerical diffusion than the results using the HLFW method for a mesh four times coarser. The AUSM+ method presented high overshoots for all mesh sizes. The results with the pressure correction term showed much less oscillation. This term stabilized the solution for refined meshes for all methods including the MUSCL second-order extensions. However, only the MUSCL-AUSMDV and HLFW methods gave oscillation-free results. The MUSCL second-order extensions were unstable without the pressure correction term. In this case, an oscillation-free solution was obtained only for the AUSMDV method.

In the large relative velocity shock tube test problem, the two separated void fraction waves could be simulated with AUSM type methods only for refined meshes. However, the Hybrid Lax-Friedrichs-Lax-Wendroff method was capable of handling this two separate void fraction scenario with coarser meshes. As in the case of the water faucet, it was not possible to obtain stable solutions for refined meshes for all methods without the pressure correction term. The results with the pressure correction term also showed more numerical diffusion, but with acceptable results for coarse meshes. However, in such meshes the two separate waves were not captured. The results for refined meshes showed that only the MUSCL-AUSMDV and HLFW methods were capable of eliminating the numerical oscillations with high accuracy results. Different from the water faucet problem, the better choices for the AUSMDV weighting function were different for the first-order and second-order schemes. In the water faucet problem, the requirements on numerical diffusion could be relaxed, allowing MUSCL-AUSMDV to give better results with the more diffusive choice for the weighting function.

The results show that it is not possible to generalize which AUSMDV weighting function is the best. This choice is case and order dependent. For all cases, AUSMDV produced additional diffusion and more stable results than the AUSM+ method. For both test problems, the best results were obtained from the MUSCL-Hancock second-order extension of AUSMDV and the Hybrid Lax-Friedrichs-Lax-Wendroff methods, with some advantage for the latter. One should bear in mind that this method has an iterative optimization of a weighting parameter, which increases accuracy but also increases the CPU time consumption. However, unlike the one-phase flow analysis, where the accuracy of the Hybrid Lax-Friedrichs-Lax-Wendroff method was comparable to AUSM type methods but with larger CPU time consumption (Coelho et al., 2006), this method is an alternative for Riemann solvers and AUSM type solvers with reasonable accuracy and CPU time for two-phase flows. As it does not require the jacobian matrix evaluation, it can be easily generalized for complex flow closure relations.

ACKNOWLEDGEMENTS

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

NOMENCLATURE

specific area A Jacobian matrix A_k EOS parameter B vector of non-conservative terms c speed of sound C_d drag coefficient C_p specific heat at constant pressure C_v specific heat at constant volume CFL Courant-Friedrichs-Lewy number, CFL = Δ t (c+max(v))/Δx f flux vector F numeric flux vector f_k^D drag force due to interface and wall on phase k f_k^Di interfacial drag force on phase k f_k^Dw wall drag force on phase k g_x component of the gravity acceleration vector in the x direction g non-homogeneous vector in conservation equations G derivaive of the convex function max flux M Mach number p pressure P_k pressure part of the flux p_i interfacial pressure S convex function t time u vector of conservative variables U vector of reconstructed variables w vector of non-conservatives variables(w = [α_g, v_g, v_l, p]^t) v velocity for one dimensional flows x spatial coordinate Y AUSMDV weighting function

Greek Symbols

α volume fraction η dissipation function σ pressure correction parameter Δ variation ε Δt/Δx ratio γ ratio of specific heats, γ = C_p/C_v Π EOS parameter ψ_k convenction part of the flux ρ density θ weighting parameter of hybrid method Θ flux limiter

Subscripts

j value at cell j j+1/2 value at x_j+1/2 (cell face) 0 initial values k denotes phase k

Superscripts

AUSM related to AUSM method L left state LW related to Lax-Wendroff method LF related to Lax-Friedrichs method n value at time instant n R right state

(Submitted: October 23, 2008 ; Revised: October 30, 2009 ; Accepted: January 4, 2010)

• Bedjaoui N. and Sainsaulieu L., Analysis of a Non-hyperbolic System Modeling Two-phase Flows Part 1: The Effects of Diffusion and Relaxation. Mathematical Methods in the Applied Sciences, 20, 783-803 (1997).
• Bendiksen, K., Malnes, D., Moe, R. and Nuland, S., The dynamic two-fluid model OLGA: theory and application. SPE production Engineering, 6, 171-180 (1991).
• Bestion, D., The physical closure laws in the CATHARE code. Nucl. Eng. Design, 124, 229-245 (1990).
• Cascales, J. R. G. and Salvador, J. M. C., Extension of a high-resolution scheme to 1D liquidgas flow. Int. J. Numer. Meth. Fluids, 50, no. 9 (2006).
• Castro C. E. and Toro E. F., A Riemann solver and upwind methods for a two-phase flow model in non-conservative form, Int. J. Numer. Meth. Fluids. 50, 275-307 (2006).
• Coelho, R. M. L., Lage, P. L.C. and Telles, A. S., A comparison of hyperbolic solvers for ideal and real gas flows. Brazilian Journal of Chemical Engineering, 23, no.3, 301-318 (2006).
• Coquel F., El Amine K., Godlewski E., Perthame B. and Rascale P., A numerical method using upwind schemes for the resolution of two-phase flows. J. Comput. Phys.; 136, no. 2, 272-288 (1997).
• Cortes, J., Debussche, A. and Toumi, I., A density perturbation method to study the eigenstructure of two-phase flow equation system. J. Comput. Phys., 147, 463-484 (1998).
• De Vuyst, F., Stable and accurate hybrid finite volume methods based on pure convexity arguments for hyperbolic systems of conservation law. J. Comput. Phys., 193, 426-468 (2004).
• Evje, S. and Fjelde, K. K., Relaxation Schemes for the Calculation of the two-phase flow in pipes. Math. Comp. Modeling, 36, 535-567 (2002).
• Evje, S. and Fjelde, K. K., On a rough AUSM scheme for a one-dimensional two-phase model. Comp. Fluids, 32, 1497-1530 (2003).
• Evje, S. and Flatten, T., Hybrid flux-splitting schemes for a common two-fluid model. J. Comput. Phys., 192, 175-210 (2003).
• Faille, I. and Heintzé, E., A rough Finite Volume scheme for modeling two-phase flow in pipeline. Comput. Fluids, 28, 213-241 (1999).
• Fjelde, K. K. and Karlsen, K. H., High-resolution hybrid primitive-conservative upwind schemes for drift flux model. Comput. Fluids, 31, 335-367 (2002).
• Godunov, S. K., A difference method for numerical calculation of discontinuous equation of hydrodynamics. Math. Sbornik (in Russian), 47, 217-306 (1959).
• Guinot V., Numerical simulation of two-phase flow in pipes using Godunov method. Int. J. Numer. Meth. Engng., 50,1169-1189 (2001a).
• Guinot V., The discontinuous profile method for simulating two-phase flow in pipes using the single component approximation. Int. J. Numer. Meth. Fluids., 37, 341-359 (2001b).
• Harten, A., High resolution schemes for hyperbolic conservation laws. J. Comput. Phys., 49(3), 357-393 (1983).
• Ishii, M., Thermo-fluid Dynamics: Theory of two-phase flow. Eyrolles, Paris, 1975.
• Liou, M. S., A sequel to AUSM: AUSM+. J. Comput. Phys., 129, 364-382 (1996).
• Liou, M. S., Mass flux schemes and connection to shock instability. J. Comput. Phys., 160, 623-648 (2000).
• 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. and Steffen, C. J., A new flux splitting scheme. J. Comput. Phys., 107, 23-39 (1993).
• Lorentzen R. J. and Fjelde K. K., Use of slopelimiter techniques in traditional numerical methods for multi-phase flow in pipelines and wells. Int. J. Numer. Meth. Fluids, 48, 723-745 (2005).
• Masella, J. M., Tran, Q. H., Ferre, D. & Pauchon, C., Transient simulation of two-phase flows in pipes. Int. J. of Multiphase Flow, 24, 739-755 (1998).
• Menikoff, R. and Plhor, B. J., The Riemann problem for fluid flow of real materials. Reviews of modern physics, 61(1), 75-131 (1989).
• Niu, Y. Y., Advection upwinding splitting method to solve a compressible two-fluid model. Int. J. Numer. Meth. Fluids., 36, 351-371 (2001).
• Paillère, H., Corre, C. and Cascales, J. R. G., On the extension of the AUSM+ scheme to compressible two-fluid models. Comp. Fluids, 32, 891-916 (2003).
• Paillère, H., Kumbaro, A., Viozat, C. and Clerc, S., Development of an AUSM+-type scheme for homogeneous equilibrium two-phase flow, Sèminare "Schemes de flux pour la simulation numèrique des ecoulements diphasiques". Reunion du groupe de Travail CEA-EDF-CMLA-Cargèse, 22-24, Septembre (1999).
• Ransom, V. H., Numerical benchmark tests, Multiphase science and technology, 3 (1987).
• Roe, P. L., Approximate Riemann solvers, parameter vector, and difference schemes. J. Comput. Phys., 43, 357-372 (1981).
• Romate, J. E., An approximate Riemann solver for a two-phase flow model with numerically given slip relation. J. Comput. Phys., 27(4), 445-477 (1998).
• Sainsaulieu, L., Finite Volume approximation of two-phase fluid flows based on an approximate Roe-type Riemann solver. J. Comput. Phys., 121, 1-28 (1995).
• Saurel, R. and Abgrall, R., A multiphase Godunov method for compressible multifluid and multiphase flows. J. Comput. Phys., 150, 425-467 (1999).
• Tiselj, I. and Petelin, S., Modelling of two-phase flow with second-order accurate scheme. J. Comput. Phys., 136, 503-521 (1997).
• Toro, E. F., Riemann Solvers and Numerical Methods for Fluid Dynamics-A practical introduction. 2^nd edition, Springer, Berlin (1999).
• Toumi, I., An upwind numerical method for two-fluid two-phase flow models. Nucl. Sci. Eng., 123, 147-168 (1996).
• Toumi, I. and Kumbaro, A., An approximate linearized Riemann solver for a two-fluid model. J. Comput. Phys., 124, 286-300 (1996).
• Van Leer, B., Towards the ultimate conservative difference scheme IV. A new approach to numerical convection, J. Comput. Phys., 23, 276-299 (1977).
• Van Leer, B., Towards the ultimate conservative difference scheme V. A second-order sequel Godunov method. J. Comput. Phys., 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).
• Wackers J. and Koren B., A fully conservative model for compressible two-fluid flow. Int. J. Numer. Meth. Fluids., 47, 1337-1343 (2005).
• Wada Y. and Liou M. S., An Accurate and robust flux splitting scheme for shock and contact discontinuities. SIAM J. Sci. Comput., 3, 633-657 (1997).

Publication Dates

• Publication in this collection
  15 Apr 2010
• Date of issue
  Mar 2010

History

• Accepted
  04 Jan 2010
• Reviewed
  30 Oct 2009
• Received
  23 Oct 2008