﻿ A Gas-Kinetic Scheme for Six-Equation Two-Phase Flow Model

Applied Mathematics
Vol.5 No.3(2014), Article ID:42778,13 pages DOI:10.4236/am.2014.53045

A Gas-Kinetic Scheme for Six-Equation Two-Phase Flow Model

Saqib Zia, Munshoor Ahmed, Shamsul Qamar

Department of Mathematics, COMSATS Institute of Information Technology, Islamabad, Pakistan

Email: shamsul.qamar@comsats.edu.pk

Received October 23, 2013; revised November 23, 2013; accepted November 30, 2013

ABSTRACT

A kinetic flux-vector splitting (KFVS) scheme is applied for solving a reduced six-equation two-phase flow model of Saurel et al. [1]. The model incorporates single velocity, two pressures and relaxation terms. An additional seventh equation, describing the total mixture energy, is added to the model to guarantee the correct treatment of shocks in the single phase limit. Some salient features of the model are that it is hyperbolic with only three wave propagation speeds and the volume fraction remains positive. The proposed numerical scheme is based on the direct splitting of macroscopic flux functions of the system of equations. The second order accuracy of the scheme is achieved by using MUSCL-type initial reconstruction and Runge-Kutta time stepping method. Moreover, a pressure relaxation procedure is used to fulfill the interface conditions. For validation, the results of suggested scheme are compared with those from the high resolution central upwind and HLLC schemes. The central upwind scheme is also applied for the first time to this model. The accuracy, efficiency and simplicity of the KFVS scheme demonstrate its potential for modeling two-phase flows.

Keywords:Six-Equation Model; Kinetic Flux-Vector Splitting Scheme; Central Upwind Scheme; Non-Conservative System; Shock Solution

1. Introduction

Multiphase flow problems occur in a variety of scientific and technical disciplines ranging from environmental research to modeling of operations in nuclear, chemical or processing engineering installations and combustion systems. The modeling and simulation of these flows are one of the complex and challenging subjects of computational fluid dynamics (CFD). In such problems the flow regime of interest contains two or more non-mixed fluids separated by sharp interfaces. A challenging part in the numerical simulation of these models is the coupling of interface with the fluid flow model, because a coupling miss-match may introduce large errors in the numerical simulations.

Several two-phase flow models exist in the literature for describing the behavior of physical mixtures. These models use separate pressures, velocities, and densities for each fluid. Moreover, a convection equation for the interface motion is coupled with the conservation laws of flow models. In the literature such models are known as seven-equation models. One of such models for solid-gas two-phase flows was initially introduced by Baer and Nunziato [2] and was further investigated by Abgrall and Saurel [3,4], among others. The seven-equation models are considered as well-established two-phase flow models. However, such models inherit a number of numerical complexities. To resolve such difficulties researchers have proposed reduced models ranging from three to six-equation models [5-7].

The Kapila’s five-equation model [5] derived from Baer and Nunziatoseven-equation model [2] is a wellknown reduced model and has been successfully implemented to study interfacing compressible fluids, barotropic and non-barotropic cavitating flows. The model contains four conservative equations, two for mass conservation, one for total momentum and one for total energy conservation. The fifth equation is non-zero convection equation for volume fraction of one of the two phases.

Although, this five-equation model is simple, it involves a number of serious difficulties. For example, the model is non-conservative and, hence, it is difficult to obtain a numerical solution which converges to the physical solution. In the presence of shocks, the volume fraction may become non-positive. Another issue is related to non-conservative behavior of the mixture sound speed which is governed by the Wood formula [8].

To avoid such difficulties of five-equation model, Saurel et al. [1] introduced pressure non-equilibrium terms which are already present in the seven-equation model of Baer and Nunziato [3]. This reduced model is six-  equation model incorporating single velocity, two pressures and relaxation terms of two-phases. An additional seventh equation, describing the total mixture energy, is added to the model which guarantees the correct treatment of the shocks in the single phase limit. Some salient features of this extended model are that it is hyperbolic with only three wave propagation speeds and the volume fraction remains positive. In [1], the authors used Godunov-type method equipped with acoustic solver and HLLC-type solver. In addition a pressure relaxation procedure is implemented to fulfill the inter face conditions.

Kinetic flux-vector splitting (KFVS) schemes have been successfully used for numerically studying flows in gas dynamics. In [9] KFVS scheme is used to solve bump in a channel problem on structured meshes. Numerically, it was found in [9] that the explicit flux function of KFVS scheme, by employing collisionless Boltzmann transport equation, is similar to the flux function of van Leer [10]. The same result was first obtained by [11]. Moreover, a high order KFVS scheme for the simulation of several two-dimensional problems on structured and unstructured meshes is in [12]. Recently, authors in [13] and [14] have constructed an improved BGK-type KFVS scheme which also incorporates particle collisions at the cell interfaces. Furthermore, [15] and [16] used different KFVS schemes for solving shallow water equations. Very recently, in [17] KFVS schemes are employed to solve a reduced five-equation two-phase flow model.

In this work, the six-equation model of Saurel et al. [1] is numerically investigated. Our primary objective is to develop a robust and efficient numerical method which is also easy to implement. A kinetic flux vector splitting (KFVS) scheme is proposed for solving this model. The suggested scheme is based on the splitting of macroscopic flux functions of the equations of two-fluid flow model. The upwinding bias in the numerical flux function can be naturally obtained by considering a fluid as a collection of particles. The movements of particles in forward or backward directions automatically split the fluxes of mass, momentum and energy into forward and backward fluxes across the cell interfaces, i.e.

where represents a vector of mass, momentum and energy densities inside the cell. In this scheme, we start with the cell averaged initial data of conserved variables and get back their cell averaged values in the same cell at the next time step. In order to get second-order accuracy, MUSCAL-type initial reconstruction and Runge-Kutta time stepping method are employed.

For validation, the numerical results of KFVS scheme are compared with those obtained from the centralupwind scheme [6]. The basic idea of such schemes is that they use information of local propagation speeds and estimate the solution in terms of its cell average. Further, such schemes have an upwind nature, since they take care of directions of wave propagation by measuring the one-sided local speeds. Moreover, to further authenticate our results, HLLC-type Riemann solver [18] is also implemented to the current model.

This article is organized as follows. In Section 2, the one-dimensional single velocity and two pressures two-phase flow model is presented. Afterwards, in Section 3 the one-dimensional KFVS scheme is derived for the numerical approximation of this model. In Section 4, the one-dimensional central-upwind scheme is briefly introduced. Section 5 deals with the pressure relaxation procedure and correction step. In Section 6, some numerical test problems are presented to validate the numerical results. Finally, Section 7 gives conclusions and remarks.

2. Mathematical Model

In this section, a single velocity six-equation two-phase flow model is presented [1]. This model was obtained from [2] by considering the asymptotic limit of stiff velocity relaxation. However, the non-equilibrium pressure terms are kept in the model. The six-equation model has some features which are the basic requirements for the numerical approximation of hyperbolic type systems. For example, the positivity of volume fraction and monotonic property of the mixture sound speed. In one space dimension, the six-equation model without heat and mass transfer is given as [1]:

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

Here, is the mixture density, is the velocity, is the mixture pressure and is the volume fraction of phase. The mixture total energy is given as, with mixture internal energy e is defined by

(2)

where is the mass fraction for phase k. Further, the mixture density is define by , where according to stiffened gas equation of state (EOS)

(3)

Here, and are material specific constants of the fluid. From Equation (2) one can find the mixture pressure as under

(4)

Here, is the interfacial pressure which is defined as

(5)

where is the acoustic impedance of phase. The combination of the two internal energy equations with mass and momentum equations in Equation (1) give an additional mixture energy equation additional equation for the mixture energy

(6)

where .

As the model (1) is non-conservative, the non-conventional jump conditions are required for the convergence of numerical methods which are based on Reimann solver. These jump conditions are derived in [1] for HLLC scheme. It is worth mentioning that the proposed KFVS and central-upwind schemes do not require such jump conditions.

As it is mentioned earlier that the model under consideration has an attractive feature regarding the mixture sound speed, the mixture sound speed behaves monotonically against mass and volume fractions and is called as the frozen sound speed. It is defined as

(7)

where the phasic sound speed is given as

(8)

The Equation (1a) can be written in the following form

(9)

In view of Equations (1) and (9), the model (1) can be rewritten as

(10a)

where

(10b)

The term in the volume fraction equation and the terms in the two energy equations for

are treated in the same manner as the flux vector terms.

To close the system (1) and (6), additional equations are needed. In the present case, it is assumed that each phase is described by the stiffened gas EOS given by Equation (3).

3. The KFVS Scheme

In this section, the KFVS scheme is implemented to solve the system of Equations (1) and (6). In gas-kinetic theory, the flux is related to the particle motion across a cell interface. The numerical discretization of the current system corresponds to the evaluation of local macroscopic flux-vector through each boundary of the mesh cell. The particle motion in the x-direction determines the flux function. The remaining quantities, such as densities, pressure and volume fraction can be considered as passive scalars transporting with the particle velocity. Normally, particles are randomly distributed around the average velocity.

In statistical mechanics, the local Maxwellian distribution function is used to describe the distribution of moving particles along coordinate directions. The Maxwellian distribution function in the normal direction is given as (e.g. [14])

(11)

where ρ and p are mixture density and pressure. In the one-dimensional case. The transport of any flow quantity is due to the movement of particles. Thus, in the one-dimensional case, with the distribution function given by Equation (11) one can split the particles into two groups. One group is moving to the right with positive velocity and the other group is moving to the left with negative velocity. Before splitting the fluxes, let us define

(12)

(13)

The above two moments are sufficient to split all the fluxes. The zeroth order moment (12) is used to split scalars while the first order moment (13) is used for vectors. In order to simplify the notation, we define

(14)

(15)

and

(16)

(17)

In the above equations, the positive sign represents those particles moving in the positive (right) direction and the negative sign represents the particles moving in the negative (left) direction. Moreover, the complementary error function is defined as

(18)

With the help of above flux-splitting technique, we can derive a KFVS scheme for solving the system of Equations (1) and (6).

In order to apply a finite volume scheme, the first step is to subdivide the domain of interest into subdomains or mesh cells. Let us define the cell by interval for Therefore,

represents the uniform cells width, the points refer to the cells center and the points

represent the cells faces. We start with a cell averaged initial data at time step and compute the cell average updated solution over the same cells at the next time step. This is performed easily by assuming the CFL condition

(19)

where is the speed of sound defined in Equation (3). With the help of above defined setup we can split the flux functions of the six-equation model in Equation (10) as

(20)

where

(21)

where. Here, the flux-vector at the right interface of the cell is given as

(22)

Analogously, we can define the left interface flux-vector of the cell. The integration of the Equation

(10) over the cell gives the following semi-discrete kinetic upwind scheme

(25)

Here, we have assumed the cell averaged values of the non-differential source term. The cell averaged values are defined as

(24)

Similarly, can be defined and are given by (22). Moreover,

(25)

The splitting procedure in Equation (25) is analogous to that presented in (20)-(22) for cell interface fluxes.

The above scheme is first order accurate in space. To achieve high order accuracy, the initial reconstruction strategy must be applied for interpolating the cell averaged variables. Here, a second order accurate MUSCL-type initial reconstruction procedure is employed. Starting with a piecewise-constant solution, , one can reconstruct a piecewise linear (MUSCL-type) approximation in x-directions by selecting the slope vector (differences). The boundary extrapolated values are given as

(26)

A possible computation of these slopes, is given by family of discrete derivatives parameterized with, for example

(27)

where is a parameter and denotes central differencing

(28)

Here denotes the min-mod non-linear limiter

(29)

On the basis of above reconstruction procedure, a semi-discrete high resolution kinetic solver is given as

(30)

To obtain second order accuracy in time, we use a second order TVD Runge-Kutta scheme to solve (30). Denoting the right-hand side of (30) as, a second order TVD Runge-Kutta scheme update through the following two stages [13]

(31a)

(31b)

where is a solution at previous time step and is updated solution at next time step. Moreover, represents the time step.

4. Central-Upwind Scheme

In the literature central schemes are regarded as suitable numerical tools to approximate the hyperbolic conservation laws. In this work, to validate the numerical results of KFVS scheme, the central-upwind scheme is for the first time applied to solve the given model [19].

This scheme belongs to the new family of central schemes introduced by Kurganov and Tadmor [19]. The core concepts behind the development of such schemes are the use of information of local propagation speeds and estimate the solution in terms of its cell average. Moreover, these schemes have an upwind nature, since they respect the directions of wave propagation by measuring the one-sided local speeds. Due to this feature, they are called as central-upwind schemes. Here, we present the final formulation of the scheme, see [19] for further details regarding the scheme derivation. In semi-discrete form the scheme is given as [19]

(32)

where is the numerical flux defined by

(33)

The intermediate values and are given by Equation (26). Here, represents the maximal local speed which in the generic case could be

(34)

Moreover, and.

5. Pressure Relaxation Procedure

The pressure relaxation plays a vital role to fulfill the interface conditions. Here, the pressure relaxation procedure of Saurel et al. [10] is used. For this purpose, we have to solve the following system, which is obtained by applying the Strang splitting approach [20] to the system (1) and (6), i.e.

(35a)

(35b)

(35c)

(35d)

(35e)

(35f)

(35g)

where.

Using the volume fraction Equation (35) with mass Equations (35b) and (35c), the internal energy Equations (35e) and (35f) can be written as under

(36)

where is the mass fraction of phase. Integration of these equations results in the following approximations

(37)

Here, the super scripts “0” and “*” indicate the values of the variables before and after relaxation process, respectively. To approximate the average interface pressure, the following estimates can be used

(38)

where is the relaxed pressure. However, as mentioned in [10], the choice of this approximation has no significant effect on the computations.

Clearly, Equation (37) comprises of two equations. But, there are three unknowns, and. To close the system we need an extra equation, which could be the saturation constraint

(39)

But

because are constant during the relaxation process. Therefore, Equation (39) can be written as

(40)

Using the EOS in Equation (3) into Equation (37) the specific volume can be expressed as

(41)

Now, using Equation (41) into (40), we get single equation with only one unknown

(42)

To find the relax pressure, the Newton’s method is used. After having the relax pressure, the specific volume and volume fractions for two phases can be found.

Correction Step

Since the volume fractions have been approximated through relaxed pressure. The mixture pressure can be obtained from Equation (4)

(43)

The mixture pressure in Equation (4) is obtained from the evolution of the mixture total energy in Equation (6). As the energy equation is in conservative form, therefore, it is expected to be valid in the entire flow field. Thus, the mixture pressure is used with EOS for each phase to reset the value of the internal energies.

6. Numerical Test Problems

In this section, three one-dimensional numerical test problems are presented to validate the results obtained by the application of KFVS and central-upwind schemes to six-equation two-phase flow model. For comparison, we have used the HLLC-type Riemann solver [1].

6.1. Sod’s Problem

This problem is similar to the Sod’s problem in the single phase gas dynamics. In this problem, the density and pressure ratios are large and the left and right gases have different ratios of specific heats. Both gases are separated by a very thin membrane placed at and are initially at rest. The left side gas has high density and pressure compared to that on the right side of the membrane. After removing the membrane, the gases evolution in time takes place. The initial data are given as

The ratio of specific heats for the left and right side gases are taken as and, respectively. The numerical results at are shown in Figure 1. The solution contains a left-going rarefaction wave, right-going shock wave and the right-moving two-fluid interface. In Figure 1, the solutions of KFVS and central upwind schemes are compared at mesh cells. The reference solution is obtained from the same central upwind scheme at 2000 grid points. Both schemes give correct location of the discontinuities and have comparable accuracy. Moreover, no pressure oscillations are observed in the solution.

Table 1 gives -errors in the density for KFVS and central schemes for different number of grid points. It is clear from the results that KFVS scheme produces less errors in the solutions compared to the central scheme.

6.2. Shock Tube Problem

In this test problem, a shock-tube of is considered. It is filled with a liquid under high pressure on the left and with the vapors at atmospheric pressure on the right. This test problem was also studied in [21]. The initial data are

For numerical reasons, the presence of a small volume fraction of the other fluid on both sides of the diaphragm is allowed, say. The results are computed at 1000 mesh cells with final simulation time

. Here, ratios of specific heats are given as and. Further, ,

Figure 1. Results of problem 6.1 at t = 0.015.

Table 1. L1-Errors in the density for problem 6.1.

, and.

The same problem was also solved by HLLC-type Riemann solver. The numerical results are shown in Figure 2. It can be observed that proposed schemes give comparable results. However, small overshoot can be observed in the results of KFVS and central schemes at the interface in the velocity profile. Overall, the results of all three schemes agree well.

6.3. Two-Fluid Mixture Problem

The initial values are given as

Here, ratio of specific heats are and.

This problem is considered as a hard test problem for a numerical scheme. A left moving rarefaction wave, a contact discontinuity, and a right moving shock wave can be observed in the solution. The right moving shock hits the interface at. The shock continues to move towards right and a rarefaction wave is created which is moving towards left. We choose 400 mesh cells and the final simulation time is taken as. The numerical results are shown in Figure 3. The figures show that both schemes give comparable results. However, KFVS scheme has better resolved discontinuities in the solution.

6.4. Comparison Test Problem

In this section we study a water-air mixture to analyze and compare the results obtained from five-equation [6,18], six-equation [1] and seven-equation [4] models. The initial values of the problem are

Figure 2. Results of problem 6.2 at t = 0.000473.

Figure 3. Results of problem 6.3 at t = 0.012.

Here, , , , , and CFL = 0.5. This problem is simulated for three models. The KFVS is implemented for five-equation model and HLLC-type Riemann solver is used for six-equation and seven-equation models. The results are computed at 200 mesh cells and the final simulation time is. Figure 4 shows the results for comparison of three models.

Here, the five-equation model solution is considered as the reference solution. It is observed that the solutions obtained from six-equation model do not converge while the seven-equation model results are comparatively better. Such situations may arise in extreme initial conditions in mixtures. Therefore, the schemes for six-equation and seven-equation models may give wrong shock speed and wrong jumps. This observation is also reported in [22]. Thus, further work is needed to tackle such extreme conditions.

7. Conclusion

In this article, a reduced six-equation two-phase flow model of Saurel et al. [10] was numerical investigated. The model incorporates single velocity, two pressures and relaxation terms. An additional seventh equation, describing the total mixture energy, is added to the model to guarantee the correct treatment of shocks in the single phase limit. Some salient features of the model are that it is hyperbolic with only three wave propagation speeds and the volume fraction remains positive. The high resolution kinetic flux-vector splitting (KFVS) scheme and central upwind scheme were applied to solve the given model. The KFVS is based on the direct splitting of macroscopic flux functions of the system of equations. The central upwind scheme belongs to the new family of central schemes introduced by Kurganov and Tadmor [7]. The scheme utilizes information about the local propagation speeds and estimates the solution in terms of its cell average. Moreover, the scheme has upwind nature in semi-discrete form, since it respects the direction of wave propagation by measuring the one-sided local speeds. The second order accuracy of the suggested schemes was achieved by using MUSCL-type initial reconstruction and Runge-Kutta time stepping method. Moreover, a pressure relaxation procedure was used to fulfill the interface conditions. The numerical results of the schemes were compared with each other and HLLC scheme. Both schemes showed comparable performance in the selected test problems. The accuracy, efficiency and simplicity of the KFVS and central schemes demonstrate their potential for modeling two-phase flows. This work was a first step towards the approximation of the full seven-equation model by KFVS and central upwind schemes. The seven-equation model is non-conservative and non-strictly hyperbolic. Therefore, this model gives tough time to a numerical scheme. The current experience with the reduced model will be useful to solve the full seven-equation models more efficiently and accurately. Work is in progress in this direction and will be presented soon for publication.

Figure 4. Results of problem 6.4 at t = 0.00002.

Acknowledgements

A partial support from the Higher Education Commission of Pakistan at the COMSATS Institute of Information Technology Islamabad is gratefully acknowledged.

REFERENCES

1. R. Saurel, F. Petitpas and A. Berry, “Simple and Efficient Relaxation Methods for Interfaces Separating Compressible Fluids, Cavitating Flows and Shocks in Multiphase Mixtures,” Journal of Computational Physics, Vol. 228, No. 5, 2009, pp. 1678- 1712. http://dx.doi.org/10.1016/j.jcp.2008.11.002

2. M. R. Baer and J. W. Nunziato, “A Two-Phase Mixture Theory for the Deflagration-to-Detonation Transition (DDT) in Reactive Granular Materials,” International Journal of Multiphase Flow, Vol. 12, No. 6, 1986, pp. 861-889. http://dx.doi.org/10.1016/0301-9322(86)90033-9

3. R. Abgrall and R. Saurel, “Discrete Equations for Physical and Numerical Compressible Multiphase Mixtures,” Journal of Computational Physics, Vol. 186, No. 2, 2003, pp. 361-396. http://dx.doi.org/10.1016/S0021-9991(03)00011-1

4. R. Saurel and R. Abgrall, “A Multiphase Godunov Method for Compressible Multifluids and Multiphase Flows,” Journal of Computational Physics, Vol. 150, No. 2, 1999, pp. 425-467. http://dx.doi.org/10.1006/jcph.1999.6187

5. A. K. Kapila, R. Menikoff, J. B. Bdzil, S. F. Son and D. S. Stewart, “Two-Phase Modeling of Deflagration-to-Detonation Transition in Granular Materials,” Physics of Fluids, Vol. 13, No. 10, 2001, pp. 3002-3024. http://dx.doi.org/10.1063/1.1398042

6. J. J. Kreeft and B. Koren, “A New Formulation of Kapila’s Five-Equation Model for Compressible Two-Fluid Flow, and Its Numerical Treatment,” Journal of Computational Physics, Vol. 229, No. 18, 2010, pp. 6220-6242. http://dx.doi.org/10.1016/j.jcp.2010.04.025

7. J. Wackers and B. Koren, “Five-Equation Model for Compressible Two-Fluid Flow,” Centrum Wiskunde & Informatica, 2004, CWI-Report: MAS-E0414.

8. A. B. Wood, “A Textbook of Sound,” G. Bell Sons Ltd., London, 1930.

9. J. C. Mandal and S. M. Deshpande, “Kinetic Flux-Vector Splitting for Euler Equations,” Computer and Fluids, Vol. 23, No. 2, 1994, pp. 447-478.

10. B. Van Leer, “Flux Vector Splitting for the Euler Equations,” ICASE, 1982, Report No. 82-30.

11. A. Harten, P. D. Lax and B. Van Leer, “On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws,” SIAM Review, Vol. 25, No. 1, 1983, pp. 35-62. http://dx.doi.org/10.1137/1025002

12. N. P. Weatherill, J. S. Mathur and M. J. Marchant, “An Upwind Kinetic Flux Vector Splitting Method on General Mesh Topologies,” International Journal for Numerical Methods in Engineering, Vol. 37, No. 2, 1994, pp. 623-643.

13. K. Xu, “Gas-Kinetic Theory Based Flux Slitting Method for Ideal Magneto-Hydrodynamics,” Journal of Computational Physics, Vol. 153, No. 2, 1999, pp. 334-352. http://dx.doi.org/10.1006/jcph.1999.6280

14. T. Tang and K. Xu, “A High-Order Gas-Kinetic Method for Multidimensional Ideal Magnetohydrodynamics,” Journal of Computational Physics, Vol. 165, No. 1, 2000, pp. 69-88. http://dx.doi.org/10.1006/jcph.2000.6597

15. H. Tang, T. Tang and K. Xu, “A Gas-Kinetic Scheme for Shallow-Water Equations with Source Terms,” Zeitschrift für angewandte Mathematik und Physik, Vol. 55, 2004, pp. 365-382.

16. K. Xu, “A Well-Balanced Gas-Kinetic Scheme for the Shallow-Water Equations with Source Terms,” Journal of Computational Physics, Vol. 178, No. 2, 2002, pp. 533-562. http://dx.doi.org/10.1006/jcph.2002.7040

17. S. Qamar and M. Ahmed, “A High Order Kinetic-Flux-Vector Splitting Method for Reduced Five-Equation Model of Compressible Two-Phase Flows,” Journal of Computational Physics, Vol. 228, No. 24, 2009, pp. 9059-9078. http://dx.doi.org/10.1016/j.jcp.2009.09.010

18. E. F. Toro, “Riemann Solvers and Numerical Methods for Fluid Dynamics—A Practical Introduction,” 2nd Edition, Springer-Verlag, Berlin Heidelberg, 1999. http://dx.doi.org/10.1007/978-3-662-03915-1

19. A. Kurganov and E. Tadmor, “New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection Diffusion Equations,” Journal of Computational Physics, Vol. 160, No. 1, 2000, pp. 241-282. http://dx.doi.org/10.1006/jcph.2000.6459

20. G. Strang, “On the Construction and Comparison of Difference Schemes,” SIAM Journal on Numerical Analysis, Vol. 5, No. 3, 1968, pp. 506-517. http://dx.doi.org/10.1137/0705041

21. N. Andrianov, “Analytical and Numerical Investigation of Two-Phase Flows,” Ph.D. Thesis, Otto-von-Guericke University, Magdeburg, 2003.

22. A. Zein, “Numerical Methods for Multiphase Mixture Conservation Laws with Phase Transition,” Ph.D. Thesis, Fakultat fu r Mathema tik Otto-von-Guericke Universita t Magdeburg, 2010.