# Numerical Simulations For Non Conservative Hyperbolic System. Application to Transient Two-Phase Flow with Cavitation Phenomenon.

1 Introduction

This paper is concerned with the transient two-phase flow in pipes. The main interest focuses on the special case of transient flow initiated by a rapid-closing valve located in a pipe containing dissolved gas in water. The pressure fluctuations generated by such closing are often associated with the presence of domains with lower pressures involving the growth of gas bubbles, and subsequent the cavitation phenomena whenever the pressure reaches the vapor pressure. In many cases, the modeling of the problem is more difficult because the flow may turn into complex configurations if the released gas bubbles or the vapor coalesce and form large pockets.

The present work is restricted to the transient flow where only the diffusion of gas occurs, leading to a topology constituted with gas bubbles dispersed in liquid. Numerous studies have been made for this case [2,6,9,10,14,15,16,17,18]. However, for the sake of simplicity, the velocities of the phases are almost always assumed to be equal. Our work aims to investigate this problem by considering a complete modelization taking into account the velocity differences between the phases. The mathematical model presented is an extension of the model established by Biesheuvel and Wijngaarden [3] for a suspension of gas bubbles in liquid. By comparison with this model, the compressibility of liquid and the section of pipe are not constant, but functions of the pressure. Furthermore, the diffusion of gas is taken into account due to the presence of gaseous cavitation. Unlike the previous works in which mathematical models are a set of conservation laws, the model presented has a non conservative form. Consequently, the numerical methods of conservative hyperbolic systems cannot be applied in this case.

To solve the system of equations, we develop a finite difference method including a specific treatment of the boundary conditions. The validity of the model and the computational procedures is demonstrated by the comparison between our numerical results and the experimental results obtained by Wiggert and Sundquist [18].

2 Governing equations

As mentioned above, the transient flow is assumed to be a dispersed two-phase nature. To describe such flow, we adopt the model of Biesheuvel and Wijngaarden [3] which we complete by adding some relationships about the mass transfer, the two-fluids and the pipe. Here the relationships correspond to the gaseous cavitation conditions investigated by Wiggert and Sundquist [18]. The set of equations is formulated as follows:

Conservation of mass

[mathematical expression not reproducible], (2.1)

where [[alpha].sub.k], [[rho].sub.k], [U.sub.k] are the volumetric fraction, the density and the velocity of the phase k = g,l. [GAMMA] and S are the rate of gas release and the pipe cross sectional area.

Conservation of momentum for the mixture: d d dPi

[mathematical expression not reproducible] (2.2)

where [P.sub.l], [T.sup.w.sub.l], g, [theta] are respectively the pressure of the liquid, the term of friction on the pipe wall, the gravitational acceleration and the angle of inclination of the pipe.

Momentum of gas bubble:

[mathematical expression not reproducible], (2.3)

where [tau] and R are respectively the volume and the radius of a bubble. [U.sub.0] is the average velocity of the fluids.

These equations are supplemented by the following constitutive laws.

Volumetric fraction:

[[alpha].sub.g] + [[alpha].sub.l] = 1. (2.4)

Equilibrium of pressures:

[P.sub.l] = [P.sub.g] + [P.sub.v],

where [P.sub.v] is the vapor pressure, the term of surface tension being neglected.

Average velocity:

[U.sub.0] = [[alpha].sub.g] [U.sub.g] + [[alpha].sub.l] [U.sub.l].

Compressibility of fluids:

1/[[rho].sub.k] d[[rho].sub.k]/d[P.sub.k] = 1/[K.sub.k], (2.5)

where [K.sub.k] is the bulk modulus of the phase k = g, l.

Cross sectional area-strain:

DS/SdPl = D/eE, (2.6)

where e, E and D are respectively the thickness of the pipe wall, the Young's modulus for the pipe material and the circle pipe diameter.

Viscous stress:

[T.sup.w.sub.l] = 1/2 [[lambda].sub.[rho]l] [U.sub.l][absolute value of [U.sub.l]]/D,

where [lambda] is the friction coefficient.

Henry's law of diffusion:

[mathematical expression not reproducible]

where [P.sub.s] is the gas saturation pressure, [beta] and [K.sub.r] are two constants.

Gas isothermal behavior:

[partial derivative] ([P.sub.g][tau])/[partial derivative]t = 0. (2.7)

3 Non-conservative formulation of the equations

To analyse the above system of equations, it is appropriate to introduce the unknown vectors:

W = [([W.sub.1], [W.sub.2], [W.sub.3], [W.sub.4]).sup.T] = [([S[alpha].sub.l[rho]l], [S[alpha].sub.g] [[rho].sub.g], [S[alpha].sub.lpl] [U.sub.l], 1/2([U.sub.g] - 3[U.sub.0])).sup.T].

The variables [W.sub.k], k = [bar.1,4] can determine all the other variables in particular the void fraction ag, the pressure of gas P and the velocities of liquid and gas [U.sub.l], [U.sub.g]. There exists an interesting relationship between the pressure P, the void fraction ag and their associated conservative variables [W.sub.1] and [W.sub.2]. We present below some mathematical properties derived from this relationship that can be used in any numerical or theoretical analysis of the equations (numerical scheme, boundary conditions, analytic solution, etc.)

In the previous constitutive laws, the compressibility of the fluids and the cross sectional are not constant and depend mainly on the pressure. Following (2.5) and (2.6), these parameters can be defined according to

[mathematical expression not reproducible]

Here, the subscript refers to an initial state. Based on these functions, the differential functions of [W.sub.1] and [W.sub.2] can be written as follows:

[mathematical expression not reproducible]

where [E.sub.k,S] = -1/[[rho].sub.p] +d[[rho].sub.k]/dp + 1/S ds/dp, k = g, l. Then it comes

[mathematical expression not reproducible]

where [E.sub.M] = [[alpha].sub.g][E.sub.g,s] + [[alpha].sub.l][E.sub.l,s]. This means that

[mathematical expression not reproducible]

Now, it is possible to rewrite the equation (2.3) by using the change of variable W = ([W.sub.1], [W.sub.2], [W.sub.3], [W.sub.4]). By developing this equation, we find that

[mathematical expression not reproducible]

It follows that:

[mathematical expression not reproducible]

Using the isothermal assumption of the gas (2.7), we can write

[mathematical expression not reproducible]

This yields to

[mathematical expression not reproducible]

Consequently, the equation (2.3) with the change of variable W, can be expressed as follows:

[mathematical expression not reproducible]

where

[mathematical expression not reproducible]

On the other hand, in order to write the equations (2.1)-(2.2) as a system of conservation laws, the term S [partial derivative]P/[partial derivative]x appearing in the equation (2.2) needs to be written in the conservative form. To get this, it suffices to introduce the primitive function of S defined by

[mathematical expression not reproducible]

Then the conservative form of (2.2) can be written as:

[mathematical expression not reproducible]

Finally, the general form of (2.1)-(2.3) can be expressed as follows:

[mathematical expression not reproducible], (3.1)

where

[mathematical expression not reproducible]

3.1 Formulation with a temporal non conservative term

The set of equations given by (3.1) is characterized by the presence of a non conservative term described by a temporal derivative. From the expression of the matrix [??], it can be seen that the non-conservative term depends directly on the drift equation (2.3). It is therefore, appropriate to further simplify the system (3.1) by separating the drift equation from the others equations. This yields the following system:

[mathematical expression not reproducible] (3.3)

where the unknown vectors are respectively

[mathematical expression not reproducible]

3.2 Formulation with a spatial non conservative term

An alternative formulation of the system (3.3) can be derived if we express the non-conservative term as a spatial derivative. To achieve this, it suffices to eliminate the temporal derivative by using the conservation laws of the system. Hence, the system (3.3) can be formulated as follows:

[mathematical expression not reproducible] (3.4)

A general form of this system is given by

[mathematical expression not reproducible], (3.5)

where F(W) and [??](W) are defined by (3.2).

4 Numerical method

In this section, we propose two numerical methods for solving the formulations

(3.3) and (3.4). The construction of such numerical methods is based on the extension of the conservation laws schemes as well as the type of non conservative term: spatial or temporal derivative.

4.1 Numerical method for non conservative system with temporal derivative

To compute the system (3.3), we adopt the approach followed in [8]. It consists of using two steps for solving the system (3.3):

* Step 1 : Applying a classical numerical scheme to evaluate the conservation laws.

* Step 2 : Performing a time discretization for the drift equation. The algorithm of such method can be expressed as:

[mathematical expression not reproducible] (4.1)

where

[mathematical expression not reproducible]

Here, the function H represents the numerical flux of the MacCormack scheme.

[mathematical expression not reproducible]

The superscript n refers to the time level and the subscript j denotes the node of space. The non-conservative scheme (4.1) can also be formulated as:

[mathematical expression not reproducible] (4.2)

4.2 Numerical method for non conservative system with spatial derivative

As it has been mentioned above, the system (3.4) provides another formulation of the non conservative term by using a spatial derivative. This type of equations has been investigated by many authors particularly in the domain of multiphase flow where the equations are described by a non conservative system of equations. Following the approach used in [1,4,5,7,11,12,13], the system (3.4) needs to be written in the general form:

[partial derivative] W/[partial derivative]t + [partial derivative]F(W)/[partial derivative]x + B(W) [partial derivative]G(W)/ [partial derivative]x = S((W).

The main difficulty of this non-conservative system is to define the notion of weak solutions and the Rankine-Hugoniot Jump Condition as in the conservative case. To solve this system of equations, an extension of these concepts has been done for the non conservative case, which yields to the following numerical scheme:

[mathematical expression not reproducible]

where

[mathematical expression not reproducible]

Here, the terms [B.sub.j+1/2] and [G.sub.j+1/2] are determined from the state vectors [W.sup.n.sub.j] and [W.sup.n.sub.j+1].

By applying this scheme to the system (3.5), we obtain

[mathematical expression not reproducible]

Using the definitions of the matrix [??] and the function F, we get a formulation similar to that found in (4.2):

[mathematical expression not reproducible]

Here, the terms [A.sup.n.sub.j-1/2], [A.sup.n.sub.j+1/2] and [[??].sup.n.sub.j+1/2] are determined from the state vectors [W.sup.n.sub.j] and [W.sup.n.sub.j+1].

By taking A" = Arr+1/2 = Arr-1/2, the above scheme can be reduced to:

[mathematical expression not reproducible]

Using the approximation [[??].sup.n.sub.j+1/2] = 1/2 ([f.sup.n.sub.j] + [f.sup.n.sub.j+1]), the scheme can also be formulated as

[mathematical expression not reproducible], (4 3)

4.3 Numerical treatment of boundary conditions

In order to determine the physical parameters in the boundary nodes, we must combine the boundary conditions with the compatibility relations given by the method of characteristics presented in Appendix 1.

These relations are

[l.sub.j] (W)A(W)([partial derivative]W/[partial derivative] t + [[lambda].sub.j](W) [partial derivative]W/[partial derivative] x) = [l.sub.j] (W)S(W),

where [[lambda].sub.j] and [l.sub.j] are respectively the eigenvalues and the associated eigenvectors of the system (2.1)-(2.3) determined in Appendix 1.

The boundary conditions considered in this work correspond to a closing valve at upstream and a constant pressure at downstream. In the first case, the velocities of the two phases are taken to be zero. In the second case, the pressure is kept constant equal the reservoir pressure. Following the analysis made in Appendix 1, the computation of the boundary conditions is formulated as follows:

At the upstream:

[mathematical expression not reproducible]

where the subscript 0 refers to the first node, [[lambda].sub.j](W) [member of] {[[lambda].sub.1](W), [[lambda].sub.3](W)}.

At the downstream:

[mathematical expression not reproducible]

where the subscript J refers to the last node, [[lambda].sub.j](W)[member of]{ [lambda].sub.1](W), [[lambda].sub.2](W), [[lambda].sub.4](W)}.

Concerning the stability of the numerical method, the eigenvalues must satisfy the CFL criterion: [DELTA]t/[DELTA]x[absolute value of [[lambda].sup.n.sub.j]] [less than or equal to] 1, j = [bar.1, 4], [for all]n. To ensure the stability of the non conservative schemes, an iterative procedure has been implemented in the computer program in order to reach the CFL number allowing the convergence of the numerical computation.

5 Results and discussion

5.1 Comparison between the numerical methods

To investigate the influence of the time and space discretization, we present in the Figures 1-3 some comparisons between the numerical results obtained with the schemes (4.2) and (4.3). The simulations have been performed for different fluids and initial radius.

As illustrated by Figures 1-3, we observe that the numerical schemes (4.2) and (4.3) provide similar results. This means that the choice of the type of discretization has little influence on the numerical prediction.

5.2 Comparison between simulation and experimental data

To validate the model and the computational procedures, we present some comparisons with the experimental results of Wiggert and Sundquist [18]. All computations have been made with a distance increment [DELTA]x = L/200, L being the pipe length. The time increment is evaluated at each iteration through the CFL criterion.

In this work, we investigate numerically the experiments conducted by Wiggert and Sundquist [18] concerning the gaseous cavitation due to a valve closure. In their theoretical approaches, the authors assumed that liquid and gas move with the same velocity. The aim of the present study is to examine the effects of such cavitation on physical properties for both liquid and gas as void fraction, pressure and especially relative velocity between the phases.

The problem of the cavitation investigated by Wiggert and Sundquist [18] is an interesting case study because its occurring may initiate the propagation of serious pressure waves along the pipeline. It is therefore necessary to analyze the extreme pressure fluctuations caused by these waves in order to understand how the relative velocity is affected, particularly in the critical situations of downstream or upstream pressure.

The present model includes the velocity differences between the phases via a relationship depending on the initial radius of bubbles [R.sub.0]. Because of the lack of experimental values, the computations were made with different values of [R.sub.0]. As shown in Figure 4, there are a few influences of the initial bubble radius on the pressure fluctuations. Concerning the velocities of the two phases, it seems that the velocity of liquid does not vary with [R.sup.0], then only one graph is reported. However, the velocity of gas is very affected by the changes of [R.sub.0], in particular when a strong pressure reduction occurs in the fluid.

It can be seen from the Figure 4, that the region near the closed valve is subject to a strong cavitation, accompanied by a decrease in pressure and an increase of the relative velocity between the phases. The maximum amplitude of relative velocity is observed during the first wave cycle after which it becomes negligible. This indicates a causal relationship between the cavitation and the increase of the relative velocity amplitude. The effect of the cavitation on the pressure transient is quite evident as shown by the numerical simulations presented in Figure 4. Herein, the computational results reproduce the typical events that can be encountered when the cavitation takes place in the fluid. It can be noted for all the fluids considered that:

* The pressure peaks were created essentially at the closed valve.

* There exists many pressure surges that followed the first pressure peak as the consequence of the reflection of the pressure wave between the closed valve and the reservoir.

* The pressure period is longer for downstream surge and shorter for upstream surge.

* Because of the damping effect, the pressure surge magnitudes were progressively reduced with time.

To investigate the influence of the gas content on the pressure transient, a comparative study was made in Figure 4 for different fluid mixture. The numerical simulations were carried out with a fluid mixture consisting of water and air or carbon dioxide. The gas concentrations employed in the simulation are identical to those used by experiments of Wiggert and Sandquist [18].

In the Figure 4, the effect of gas content is clearly demonstrated by the distinct characteristics of pressure surges in terms of amplitude, damping or period.

It is observed that: the first pressure peak is much higher with lower gas concentration, this is the case for air-water mixture; the damping of the pressure surges depends mainly on gas content and increases grossly with it; the pressure surge periods were increased proportionally with gas content, which can be observed in the case of C[O.sub.2]-water mixture with large concentration.

Hence, we can conclude that the dissolved gas type and concentration may affect considerably the behavior of transient pressure, and then the magnitude of relative velocity. Furthermore, the first peak pressure is amplified with small amount of gas content, while the damping and the period of pressure surges are increased with the increasing of gas content. As it has been mentioned previously, the effect of the initial bubble radius [R.sub.0] on the pressure is so small that it can be ignored. It is therefore possible to compare the experiment and numerical results in order to validate the proposed model.

In the case of air-water mixture, the comparison between experiment and numerical results indicts a good concordance for the first pics, and presents some differences on phase and amplitude for the last ones. For the C[O.sub.2]-water mixtures, a consistent agreement has been obtained specially for the larger value of initial dissolved gas. With large concentration of gas, the computed and the experimental results are compatible on amplitude as well as on phase.

In constraint to carbon dioxide, the investigation of the air-water mixture produces some numerical instabilities when the initial radius bubble [R.sub.0] becomes lower or higher than certain order. Here, the simulation was carried out with a range of values between [10.sup.-6] and [10.sup.-4]. This limit can be explained by the non-smooth character of the pressure fluctuations and the presence of severe gradient pressure due to the fast change of fluid transient between upstream and downstream pressure.

6 Conclusions

This work presents a numerical analysis of a system of equations describing the transient cavitating two-phase flow. We propose two non conservative schemes to solve this system of equations. The first scheme has been proposed in a previous work of thesis (see [8]) as an extension of the conservation laws schemes. The principle of this scheme consisted of applying the MacCormack method to solve the conservation laws and make a temporal discretization to the non-conservative equation.

The second scheme is a numerical method that has been subject of many studies. In the literature, this scheme is mostly applied to compute the Riemann problem. The novelty of our work is to apply this non conservative scheme to a complex mathematical model involving the gaseous cavitation, the liquid compressibility, the pipe elasticity and slipping between the gas bubbles and liquid. As shown earlier, the computation given by these schemes is identical for the conservation laws, and different when we treat the non-conservative equation. The comparison made between these schemes gives a similitude between the results obtained for different configurations of two phase flow.

The good agreement found between the computed results and the experimental results of Wiggert and Sundquist [18] demonstrates the validity of the proposed model and the accuracy of the numerical schemes.

In addition, further information was obtained about the evolution of the relative velocity between the two phases. The effect of the cavitation on the relative velocity has been clearly demonstrated by the computational results. It is found that under downstream conditions, the decreasing of pressure yields to the increasing of the relative velocity. Moreover, the maximum magnitude of the relative velocity occurs during the first wave cycle, after which it becomes negligible. On the other hand, the numerical computation illustrates the influence of initial radius of bubble as a parameter which could have a substantial impact on the relative velocity.

It is observed that the relative velocity magnitude vary significantly with this parameter and increases with its higher value. It should also be noted that the dissolved gas type and concentration need to be considered when evaluating the effects of cavitation. As shown in the comparative study, there is a significant influence of these parameters on the behavior of pressure as well as the relative velocity between the two phases.

https://doi.org/10.3846/mma.2019.015

References

[1] R. Abgrall and S. Karni. A comment on the computation of non-conservative products. Comput. Phys., 229(8):2759-2763, 2010. https://doi.org/10.1016/j.jcp.2009.12.015.

[2] A. Bernard-Champmartin, O. Poujade, J. Mathiaud and J.M. Ghidaglia. Modelling of an homogeneous equilibrium mixture modelHEM. Acta Applicandae Mathematicae, 129(1):1-21, 2014. https://doi.org/10.1007/s10440-013-9827-2.

[3] A. Biesheuvel and L.V. Wijngaarden. Two-phase flow equations for a dilute dispersion of gas bubbles in liquid. Fluid Mech., 148(11):301-318, 1984. https://doi.org/10.1017/S0022112084002366.

[4] A. Canestrelli, A. Siviglia, M. Dumbser and E.F. Toro. Well-balanced high-order centred schemes for non-conservative hyperbolic systems. applications to shallow water equations with fixed and mobile bed. Adv. Water Resour., 32(6):834-844, 2009. https://doi.org/10.1016/j.advwatres.2009.02.006.

[5] C.E. Castro and E.F. Toro. A Riemann solver and upwind methods for a twophase flow model in non-conservative form. Numer. Meth. Fluids, 50(3):275-307, 2006. https://doi.org/10.1002/fld.1055.

[6] P.J. Martinez Ferrer, T. Flatten and S.T. Munkejord. On the effect of temperature and velocity relaxation in two-phase flow models. ESAIM: M2AN, 46(2):411-442, 2012.

[7] U.S. Fjordholm and S. Mishra. Accurate numerical discretizations of nonconservative hyperbolic systems. ESAIM: M2AN, 46(1):187-206, 2012. https://doi.org/10.1051/m2an/2011044.

[8] A. Qadi El Idrissi. Ecoulement Transitoire en Conduite avec Prise en Compte de Phenomenes Diphasiques. Thesis, Institut National des Sciences Appliques de Lyon (France), 1996.

[9] A. Jafarian and A. Pisheva. An exact multiphase Riemann solver for compressible cavitating flows. Multiphase Flow, 88:152-166, 2017. https://doi.org/10.1016Zj.ijmuItiphaseflow.2016.08.001.

[10] D. Liuzzi. Two-Phase Cavitation Modelling. Thesis, University of Rome--La Sapienza, 2012.

[11] S.T. Munkejord, S. Evje and T. Flatten. A Musta scheme for a nonconservative two-fluid model. SIAM J. Sci. Comput., 31(4):2587-2622, 2009. https://doi.org/10.1137/080719273.

[12] M.L. Munoz-Ruiz and C. Pares. Godunov method for nonconservative hyperbolic systems. Math. Model. Numer. Anal, 41(1):169-185, 2007. https://doi.org/10.1051/m2an:2007011.

[13] C. Pares. Numerical methods for nonconservative hyperbolic systems: A theoretical framework. SIAM J. Numer. Anal., 44(1):300-321, 2006. https://doi.org/10.1137/050628052.

[14] A. Prosperetti and L.V. Wijngaarden. On the characteristics of the equations of motion for a bubbly flow and the related problem of critical flow. Eng. Math., 10(2):153-162, 1976. https://doi.org/10.1007/BF01535658.

[15] O.V. Vasilyev R.S. Lagumbay and A. Haselbacher. Homogeneous equilibrium mixture model for simulation of multiphase/multicomponent flows. Numer. Meth. Fluids, 00:1-6, 2007.

[16] N.D. Tam. Fluid Transients in Complex Systems with Air Entrainment. Thesis, National University of Singapore, 2009.

[17] H.S. Tang and D. Huang. A second-order accurate capturing scheme for 1D inviscid flows of gas and water with vacuum zones. Comput. Phys., 128:301318, 1996. https://doi.org/10.1006/jcph.1996.0212.

[18] D.C. Wiggert and M.J. Sundquist. The effect of gaseous cavitation on fluid transients. Fluids Eng., 101(1):79-86, 1979. https://doi.org/10.1115/1.3448739.

Appendix 1

1. Characteristics

To obtain the characteristics of the system (2.1)-(2.3), it is convenient to use as dependent variables [W.sub.1] = [S[alpha].sub.lpl], [W.sub.2] = [S[alpha].sub.g][[rho].sub.g], [U.sub.l] and [U.sub.g]. The relationships (2.4)-(2.7), determine all other variables. Based on these variables, the quasilinear system may be expressed as

[mathematical expression not reproducible]

where the matrix functions A and B are defined by

[mathematical expression not reproducible]

with

[mathematical expression not reproducible]

The characteristic roots A are obtained by solving the equation [absolute value of B - [lambda]A] = 0. By developing this, we show that one root is vanished ([lambda] = 0) and the three others are the roots of the algebraic equation:

[([lambda] - [U.sub.l]/[c.sub.f]).sup.3], [eta][psi][([lambda]-[U.sub.l]/[c.sub.f]).sup.2] - ([lambda]- [U.sub.l]/cf) +[eta] = 0

where

[mathematical expression not reproducible]

Using the variable X = [lambda]-[U.sub.l]/[c.sub.f], the previous cubic equation can be rewritten as

[X.sup.3] - [eta][psi][X.sup.2] - X + [eta] = 0.

To solve this equation, it is appropriate to set [xi] = X - [eta][psi]/3. Then, we get

[[xi].sup.3] + a[xi] + b = 0,

a = - 1/3 (3 + [([eta][psi]).sup.2]), b = 1/27 (-2[([eta][psi]).sup.3] - 9[eta][psi] + 27[eta].

The resolution of this equation depends on the sign of the expression [b.sup.2]/4 + [a.sup.3]/27. It is demonstrated that the roots of this polynomial are real and distinct if the following condition is satisfied : [b.sup.2]/4 + [a.sup.3]/27 < 0. In addition, the exact expression of the roots is given by

[mathematical expression not reproducible]

In the present work, the term n is assumed to be neglected because the sound velocity is often larger than the velocities of the two phases. It follows that

[mathematical expression not reproducible]

Then

[b.sup.2]/4 + [a.sup.3]/27 [??] 1/27 [(3 + [([eta][psi]).sup.2]).sup.3] < 0

This inequality implies the existence of three real and unequal roots. An approximation of these roots may be obtained by neglecting the term 3b/2a[square root of -3/a] which yields to the following expressions:

[phi] = [pi]/2, [[xi].sub.1] = [square root of -a], [xi]2 = -[square root of -a], [[xi].sub.3] = 0.

Finally, the characteristic roots can be deduced as

[mathematical expression not reproducible]

For each eigenvalue [[lambda].sub.j], we define an eigenvector [l.sub.j] as [l.sub.j] (B - [[lambda].sub.j]A) = 0. The resolution of such system gives: if [[lambda].sub.j] = 0, then [l.sub.j] = (0,0, 0,1), and if [[lambda].sub.j] [not equal to] 0, then

[mathematical expression not reproducible]

2. Compatibility relations

Unlike the interior nodes for which a numerical method has been established, the extremity nodes need a specific treatment based on the compatibility relations. The formulation of these compatibility relations must take into account the boundary conditions imposed by experiment.

In the present work, we investigate the case of a rapid closing valve situed in the upstream of a pipeline located between two reservoirs. This case study has been studied experimentally and numerically by Wiggert and Sundquist [18].

Generally for such physical experiment, the boundary conditions are expressed as follows: the velocity vanish for the two phases at the upstream side and the pressure is maintained constant at the downstream side. To determine the physical parameters at the extremity nodes, we must combine the above boundary conditions with the following compatibility relations:

[l.sub.j] (W)A(W)([partial derivative]W/[partial derivative]t + [[lambda].sub.j](W) [partial derivative]W/[partial derivative]x) = [l.sub.j] (W)S (W),

where [[lambda].sub.j] and [l.sub.j] are respectively the eigenvalues and the associated eigenvectors of the system(2.1)-(2.3). In practice, we have to use the negative characteristics for the upstream side, and the positive characteristics for the downstream side.

2.1 Boundary conditions at upstream

By definition, the boundary conditions for a rapid closing valve are given by [U.sub.l] = [U.sub.g] = 0. By substituting these values in the polynomial characteristics, we get the following eigenvalues: [[lambda].sub.1] = [[lambda].sub.2] = 0, [[lambda].sub.3] = [c.sub.f], [[lambda].sub.4] = [c.sub.f]. For these eigenvalues, the associated eigenvectors are: [l.sub.1,2] = (0, 0, 0,1),

[mathematical expression not reproducible]

At the upstream, we determine the pressure and the void fraction by solving the compatibility relations using the negative eigenvalues. This leads to resolve the following system:

[mathematical expression not reproducible]

2.2 Boundary conditions at downstream

At the downstream side, the pressure is assumed to be constant. The other physical parameters: [[alpha].sub.g], [U.sub.l], [U.sub.g] are determined by using the positive compatibility relations. As seen before, the comparison between the eigenvalues demonstrates clearly that [absolute value of [[lambda].sub.2]] [much less than] [absolute value of [[lambda].sub.4]]. It is therefore reasonable to simplify the compatibility relations by considering the following approximations:

dx/dt = [[lambda].sub.1] = 0, dx/dt = [[lambda].sub.2] [congruent to] 0, dx/dt = [[lambda].sub.3] < 0, dx/dt = [[lambda].sub.4] > 0.

At the downstream, the boundary conditions can be formulated as follows:

[mathematical expression not reproducible]

Appendix 2: Experiments of Wiggert and Sundquist

The experiments of Wiggert and Sundquist [18] were conducted with installation constituted with a pipeline located between two reservoirs. The experience consists of a rapid closing valve located at the downstream end or the upstream end of the pipe. The Figure 5 shows the pipe system and gives the data of experiments.

Boujemaa Achchab (a), Abdellatif Agouzal (b) and

(a) University Hassan 1, National School of Applied Sciences BP 218, Berrechid, Morocco

(b) CNRS--Institut Camille Jordan Universite Claude Bernard Lyon, France

E-mail: achchab@estb.ac.ma

E-mail: abdellatif.agouzal@univ-lyon1.fr

Received April 25, 2018; revised January 22, 2019; accepted January 23, 2019

Caption: Figure 1. Comparison between the numerical schemes (4.2) and (4.3) for Air/Water mixture with dissolved gas concentration C0=0.02.

Caption: Figure 2. Comparison between the numerical schemes (4.2) and (4.3) for C[O.sub.2]/Water mixture with dissolved gas concentration C0=0.6.

Caption: Figure 3. Comparison between the numerical schemes (4.2) and (4.3) for C[O.sub.2]/Water mixture with dissolved gas concentration C0=1.15.

Caption: Figure 4. Variations of pressure and velocities for an upstream closing valve.
```Figure 5. Pipe system and data of experiments.

Valve             Initial            Gas         Gas Content
location       Velocity (m/s)                    ([C.sub.o])

Upstream            0.77             air            0.02
Upstream             077          C[O.sub.2]         06
Upstream             077          C[O.sub.2]         115
Downstream           077             air             002
Downstream          0.77          C[O.sub.2]         059
Downstream          0.77          C[O.sub.2]         127

Valve             Reservoir           Fluid
location           pressure        Temperature
(kN/[m.sup.2])     ([degrees]C)

Upstream           1.72 (a)             16
Upstream           1.70 (a)             16
Upstream           1.75 (a)             16
Downstream         2.97 (b)             18
Downstream         3.05 (b)             18
Downstream         3.08 (b)             18
```
COPYRIGHT 2019 Vilnius Gediminas Technical University
No portion of this article can be reproduced without the express written permission from the copyright holder.