The simulation of distributed free gas in liquids using Verwey-Yu scheme

 

Sameh Mansour

 

Researcher, HRI, Egypt

Currently, research fellow, TU Braunschweig, Germany

 

 

Abstract

Describing fluid transients in which free air is present has attracted a lot of research effort in the last two decades. The numerical modeling of these flows poses special problems due to numerical dispersion and attenuation. The available methods are the method of characteristics using several kinds of interpolations, implicit methods and discrete methods. In 1993 a numerical scheme was suggested by A.Verwey and J.H. Yu. This scheme was tested to be applicable in water hammer and cavitation as well as fluid-structure interaction simulations. This scheme has the distinct advantage of being stable over wide range of courant applications. It is intended in this paper to use this scheme in the simulation of transients in liquids with free gas. The inherent appeal of this scheme is that it will be possible to carry out the simulation even when the existence of free gas in the liquid leads to reducing the wave propagation velocity, benefitting from the numerical stability of the scheme. In this way it would be possible to avoid interpolation, which introduces numerical damping. Special attention has to be paid to the numerical procedure in order to reach a stable application of this technique. The results show good accordance with those produced by delft hydraulic which were used as a benchmark in order to assess the applicability of the Yerwey-Yu scheme in this type of simulation.

 

Introduction

In liquid in which free gas is present the propagation wave velocity becomes a strong function of the pressure level, and the simulation of the transient response of a system becomes a highly nonlinear problem. Figure 1 shows the relation between the wave propagation celerity and the percentage of released air in a pipe (m) together with the pressure.

 

 

The available simulation methods are the method of characteristics, implicit methods and discrete methods. When the method of characteristics is used with a specified grid in x-t plane, interpolations are necessary either in the spatial or time directions, or both. The numerical modeling of all these methods poses special problems due to numerical dispersion and attenuation. Using Verwey-Yu scheme it will be possible to allow for propagation wave speed changes due to the released air in the fluid and pressure changes. As the scheme allows for courant changes in the range of 0.1 to 1 and gives good and stable results. First the transient equations will be numerically treated, using Taylor's series expansion in order to subtract the third order terms from the truncation error. Then the simulation will be carried out and using an iterative technique, it would be possible to determine the pressure and the amount of released air in the fluid and then calculate its effect on the wave propagation velocity.

 

Governing equations

In a discrete air model isothermal expansion and compression of the air volume is assumed. A mass conservation is used at the computational sections. Together with the momentum equation it would be possible to simulate the transient flow in a fluid filled pipe system. The governing equations can be written as follow:

 

 

in which

 

 

a propagation wave velocity m/s

A cross sectional area m2

D pipe diameter m

e pipe wall thickness m

E Young's modules of elasticity Pascal

g gravitational acceleration m/s2

H pressure head m

k fluid bulk modules Pascal

m percentage of dissolved gas -

p pressure Pascal

Q discharge m3/s

R gas solubility coefficient Pascal.Kelvin

T temperature Kelvin

ρ fluid density kg/m3

 

Methodology

Equation 3 shows that the propagation wave celerity is a function of the pressure in the pipe and the percentage of air released. Due to this fact it would not be possible to maintain a unique value of courant number during the simulation. Instead of using interpolation, as in the method of characteristic, the Verwey-Yu scheme provides a practical solution. Using this scheme it will be possible to use a unique grid of space and time. During a simulation courant number will vary between 1 and 0.1 depending on the amount of released air and the pressure inside the pipe. The scheme can be implemented after treating equations 1 and 2 according to Taylor's series analysis to describe the third order terms from the truncation error and deduct them from the original equations. Then it will be possible to carry out a stable simulation even under a very low courant conditions.

 

Numerical implementation

Using Taylor's series expansion to analyze the continuity equation, the following equation can be written:

 

 

Compare equation 4 with the original continuity equation, and using θ, ψ value equal to 0.5, the second order differential terms vanish. The truncation error remains with the third order differential terms, which are equal to :

 

 

It is clear that this important term will vanish when courant number is equal to unity or when it is possible to deduce this term from the original continuity equation. This is possible when the Verwey-Yu scheme is used as it will be possible to describe this term numerically. As this scheme consists of two steps in space and 3 steps in time.

Then the continuity equation and momentum equations can be written in the following form:

In which:

and

In which:

 

After subtracting these terms from the simulated equations a stable simulation can be carried out even under a very low courant condition.

 

Test and results

A numerical model based on this approach is built and tested against the benchmarks defined by delft hydraulic. This test involves a single horizontal, frictionless pipeline. The data of this test is shown in figure 2. The transient action is introduced by lowering instantaneously the upstream head from 60 m to 10.3 m, being the level of the pipe.

 

 

At the first instance the sharp negative pressure pulse disperses. Due to reflection at downstream boundary a positive pressure wave appears. This wave will steepen until a shock appears. The shock wave will amplify and be transferred to the upstream end of the pipeline etc. The response of the system at a distance x = 1200 m from the LHS boundary is shown in figure 3 . This figure presents the results given by delft hydraulics. Figure 4 presents the results at the same point given by the produced code, based on Verwey-Yu scheme and the concentrated cavity model. The following numerical data were used in both simulations :Δx = 200 m and Δt = 0.2 s. These data means that the simulations will be carried out under courant number equal to unity at the beginning of the simulation.

Comparing these figures, it can be seen from the results that good accordance is reached. The results given by delft hydraulics show numerical oscillations or a numerical damping of the physical oscillation as they concluded in their research. They have recommended to decrease the time step in order to decrease the numerical damping. The implementation of the Verwey-Yu scheme in the simulation of this type of flow eliminates the necessity of interpolation. This result directly in the elimination of numerical damping and keeps the physical oscillations. Using this new technique, it was possible to maintain a predefined computational grid of space and time. This was done by computing a time step from the initial steady state computations and allows the system to oscillate according to the variation in courant number, due to the change in wave speed, tacking advantage of the fact that this new technique maintains its applicability over a wide range of courant numbers.

 

Another test was performed in order to obtain further validation of the produced model.

The second test was performed also by delft hydraulics. This test has the same features as the first test but the under pressure wave was introduced by lowering the pressure head at the LHS boundary from 60 mwc to 30 mwc instantaneously. The results of this test are introduced in the form of time series of pressure head at a point X=1200m from the LHS boundary. Figure 5 presents the results given by the produced model. Some oscillations can be noticed in the time series. It is also possible to conclude that the air entrained in the system acts as damper to the pressure wave. The wave amplitude decays until the system reaches the steady state again, in this case it is 30 mwc along the pipeline. Figure 6 presents the results given by delft hydraulics. Comparing both figures, it is possible to notice some phase difference between the two solutions and to conclude that the two models perform in more or less the same manner. The phase difference can be elaborated due to the fact that in the case of modeling this test using the produced model, the under pressure wave was introduced to the system over 5 time steps. While in the case of the results presented by delft hydraulics, the under pressure wave was introduced over two time steps. This results in an initial phase difference between the two solutions. It is not possible to simulate the under pressure wave, using the produced model, with less than these 5 time steps. As the instantaneous pressure drop simulation will produce numerical oscillation due to the fact of lacking the resolution of the under pressure wave.

 

 

Conclusion and recommendation

The Verwey-Yu scheme can be used successfully for water hammer computations as well as cavitation computations. It is very important to understand the numerical behavior of the scheme before applying it to practical computations. The scheme can be used throughout a wide range of Courant Numbers (0.1<Cr<1.0). It is very important to introduce the changes in boundary conditions through enough point per wave length. This is essential to resolve the changes in boundary conditions. If the under pressure wave, which will produce cavitation, is introduced through less than 20 points per wave length, an instability problem can be produced. This instability is detected through the oscillations that show up in the time series of the pressure head. These oscillations are due to the phase error which is produced during simulations carried out under Courant conditions less than unity. The phase error produces short waves. This short waves are the oscillations that appears in the time series of pressure head. To eliminate these oscillations, two different methods can be applied.

-More weight can be added to the calculated time step. This can be applied by introducing larger values of (θ) than 0.5. The drawback of this method is that more values of (θ) than 0.505 means numerical diffusion to the solution. This method is not recommended due to the produced numerical diffusion.

-The changes in the boundary conditions should be introduced through more points per wave length than 20 . This will eliminate the production of short waves and oscillates. This method is recommended and essential to be noted and respected in all simulations that are carried out under courant number less than unity.

 

Reference

[1]                M.B.ABBOTT,"Computational Hydraulics", Ashgate Publishing LTD, England, 1992

[2]                M.B.ABBOTT and D.R.BASCO,"Computational Fluid Dynamics, An Introduction for Engineers", Longman Group UK, LTD, 1989

[3]                E.S.CLYDE and J.P.TULLIS,"Aeration Scale Effects", ASCE Spcialty Conference, Boston, MA, ASCE, New York, August 1983

[4]                DELFT HYDRAULICS,"Courses in the Hydrodynamics of Pipeline Systems", Delft, The Netherlands, November 1991

[5]                A.HEINSBROEK,"Cavitation Conceptual Model Aspects", Delft Hydraulics, Delft, The Netherlands, 1993

[6]                S.G.S.Mansour,"The Application of a Space-Compact High-Order Implicit Scheme for Cavitation Computation in Waterhammer Simulation",Hydroinformatics'96, ETH Zürich, Switzerland

[7]                S.G.S.Mansour,"The Application Verwey-Yu Scheme in FSI simulations", Hydroinformatics'98, DHI Copenhagen, Denmark

[8]                C.A.PROVOOST and E.B.WYLIE,"Discrete Gas Model to Represent Distributed Free Gas in Liquids", Delft Hydraulics, delft, The Netherlands, April 1982

[9]                J.PAUL TULLIS,"Hydraulics of Pipelines", Wiley-Interscience Publications, New York 1989

[10]           A.Verwey and J.H.Yu,"A Space-Compact High-Order Implicit Scheme for Water Hammer Simulations",XXV IAHR, Tokyo, Japan, August 1993