SciELO - Scientific Electronic Library Online

 
 número20Análisis por medio de la normalización de variables para un modelo de planificación ambiental hídrica estacional índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Revista

Articulo

Indicadores

Links relacionados

  • En proceso de indezaciónCitado por Google
  • No hay articulos similaresSimilares en SciELO
  • En proceso de indezaciónSimilares en Google

Compartir


Obras y proyectos

versión On-line ISSN 0718-2813

Obras y Proyectos  no.20 Concepción dic. 2016

http://dx.doi.org/10.4067/S0718-28132016000200007 

 

Wave speed calculation for water hammer analysis

Cálculo de la velocidad de onda para el análisis del golpe de ariete

 

John Twyman

Twyman Ingenieros Consultores, Pasaje Dos # 362, Rancagua, Chile, john@twyman.cl


In order to accurately solve the water hammer problem using the Method of the Characteristics MOC is necessary to fulfil with the so-called Courant condition which establishes mandatorily that Cn = f(a) = 1 in each pipeline of the system, where a is the wave speed. The value of Cn is dependant of a whose value depends in turn on the fluid properties (density, bulk modulus) and physical characteristics of each pipeline (elasticity modulus, diameter, wall thickness, supporting condition). Because water distribution systems usually has many different pipes, and therefore, many different wave speeds, it can be said that fulfil with Cn = 1 in each pipeline is a very difficult task, more when the solution by MOC needs a common time step At for all pipe sections of the system. A way of solution to this problem is applying the method of the wave-speed adjustment that involves modifying the value of a in each pipe section in a certain percentage up to obtain Cn = 1. With this procedure optimum results are guaranteed in numerical terms, but it is possible to say the same in physical terms? The question which arises is: what parameters within the formula of a must (or can) be changed without exceeding the characteristic values of the component material of the pipes?. This work shows that in some cases the wave speed modification can significantly alter the value of the parameters that define a, leading to values that can be physically inconsistent, fictitious or without practical application.

Keywords: wave speed, water hammer, Courant number


Para resolver en forma precisa el problema de golpe de ariete usando el Método de las Características MC es necesario que el número de Courant Cn = 1 en cada tubería del sistema. El valor de Cn depende de la velocidad de la onda a, cuyo valor depende a su vez de las propiedades del fluido (densidad, módulo de compresión) y de las características físicas de cada tubería (módulo de elasticidad, diámetro, espesor, condición de apoyo). Debido a que las redes de distribución de agua generalmente tienen muchas tuberías diferentes, y por tanto, muchas velocidades de onda distintas, cumplir con Cn = 1 en cada tubería se torna una tarea muy difícil, más aún cuando la solución mediante el MC necesita un paso de tiempo At común para todas las tuberías del sistema. Una forma de solución a este problema es aplicar el método de ajuste de la velocidad de la onda que consiste en modificar el valor de a en cada tubería en un cierto porcentaje hasta obtener Cn= 1. Con esto se garantizan resultados óptimos en términos numéricos, pero ¿es posible decir lo mismo en términos físicos?. La pregunta que se plantea es: ¿qué parámetros dentro de la fórmula de a deben (o pueden) ser cambiados sin exceder los valores característicos del material componente de la tubería?. En este trabajo se muestran algunos casos donde la modificación de a puede alterar significativamente la magnitud de los parámetros que definen su valor, dando lugar a valores que pueden ser físicamente incompatibles, ficticios o sin aplicación práctica.

Palabras clave: velocidad de la onda, golpe de ariete, número de Courant


 

Introduction

For many years the Method of the Characteristics MOC has been used for solving the transient flow in pipe networks due to its numerical efficiency, computational accuracy, and programming simplicity. However, one difficulty that arises is the selection of an appropriate time step At to use for the analysis. The challenge of selecting a time step is made difficult in pipeline systems because to calculate head and discharge in many boundary conditions it is necessary that the time step be common to all pipes. Besides, MOC requires that the ratio of the distance step Δx to the time step Δt be equal to the wave speed a in each pipe, or that the Courant number Cn = a Δt/Δx should ideally be equal to one. For most pipeline systems it is impossible to satisfy exactly the Courant requirement with a reasonable (and common) Δt because they have a variety of different pipes with a range of wave speeds and lengths (Kamey and Ghidaoui, 1997). There are two strategies to deal with this problem. The first strategy is apply the method of the wave-speed adjustment MWSA where one of the pipeline properties is altered (usually wave speed) to satisfy exactly the Courant condition. The second strategy is interpolating between known grid points allowing Courant numbers less than one. At first glance the MWSA appears simpler because is non-dissipative and non-dispersive and in theory only consists in modify the value of the wave speed in a certain percentage to meet Cn = 1. Nevertheless, this procedure distorts the physical characteristics of the problem (Ghidaoui and Karney, 1994). In other words, changing a involves altering, in physical terms, the value of one or more of the parameters that are part of its formulation such as fluid density or the elastic modulus of the constituent element of the pipe. More clearly, the modification of a in numerical terms involves altering the initial physical conditions of the system, leading to a solution that may be correct in numerical terms (to meet Cn), but incorrect in physical terms because the problem is solved using parameters with unreal magnitudes.

Governing equations of transient flow

When analyzing a volume control it is possible to obtain a set of non-linear partial differential equations of hyperbolic type valid for describing the one-dimensional 1D transient flow in pipes with circular cross-section (Chaudhry and Hussaini, 1985):

where the partial differential equations (1) and (2) correspond to the continuity and momentum (dynamics), respectively. Besides, H is the piezometric head, a is the wave speed, c = (gA/a), where g is the acceleration of gravity, A is the pipe cross-section, Q is the fluid flow and R = ƒ/2DA, f is the friction factor (Darcy-Weisbach) and D is the inner pipe diameter. The subscripts x and t denote space and time dimensions, respectively. Partial differential equations (1) and (2), in conjunction with the equations related with the boundary conditions of specific devices, describe the phenomenon of wave propagation for a water hammer event.

Wave speed

For water, without presence of free air or gas, the more general equation to calculate the water hammer wave speed magnitude in one-dimensional flows is (Watters, 1984):

with a the wave speed, K the volumetric compressibility modulus of the liquid, ρ the liquid density, e the pipe wall thickness, E the pipe elasticity modulus (Young); Ψ a factor related with the pipe supporting condition which can be calculated from general expressions (see Table 1) being the case 2 more conservative from an engineering point of view. Equation (3) supposes that:

 

Table 1: Expressions for Ψ according to the pipe supporting
condition (Watters, 1984; Pierre, 2009)

 

• Pipe has a thin internal wall, condition which is met when D/e > 40 (Watters, 1984) or when D/e > 25 (Wylie and Streeter, 1978).

• Pipe remains full of water during the transient event; that is, no separation of the water column is generated, which means that at all times the pressure is greater than the vapour pressure.

• Water has small air content, so that the magnitude of the wave speed may be assumed constant.

• The pressure is uniform across any section of the pipe. It means that inertial forces associated with radial motion of the fluid are negligible (Skalak, 1955).

Equation (3) includes Poisson’s ratio effect but neglects the motion and inertia of the pipe. This is acceptable for rigidly anchored pipe systems such as buried pipes or pipes with high density and stiffness, to name only a few.

Examples include major transmission pipelines like water distribution systems, natural gas lines and pressurized and surcharged sewerage force mains. However, the motion and inertia of pipes can become important when pipes are inadequately restrained (unsupported, free-hanging pipes) or when the density and stiffness of the pipe is small (Ghidaoui et al, 2005).

Method of the characteristics

The Method of the Characteristics MOC is an Eulerian numerical scheme (Wood et al., 2005) very used for solving the equations which goveming the transient flow because it works with a constant and, unlike other methodologies based on finite difference or finite element, it can easily model wave fronts generated by very fast transient flows. MOC works converting the computational space x - time t grid (or rectangular mesh) in accordance with the Courant condition. It is useful for modelling the wave propagation phenomena in water distribution systems due to its facility for introducing the hydraulic behaviour of different devices and boundary conditions (valves, pumps, reservoirs, etc.). Among its main advantages it can be highlighted its ease of use, speed and explicit nature, which allows calculate the variables Q and H directly from previously known values (Chaudhry, 1979; Wylie and Streeter, 1978). The main disadvantage of the MOC is that it must fulfil with the Courant stability criterion that can limit the magnitude of the time step Δt common for the entire network. In order to get Cn = 1, some pipe initial properties can be modified (length and/or wave speed). Another way is to keep the initial conditions and apply numerical interpolations with risk of generating errors (numerical dissipation and dispersion) in the solution (Goldberg and Wylie, 1983). The MOC stability criterion states that (Watters, 1984):

where Cn is the Courant number, Δt is the time step and Δx is the sub-section pipe length (Δx = L/N with L the pipe length and N the number of pipe sub-sections). In general, MOC gives exact numerical results when Cn = 1, otherwise, it generates erroneous results in the way of attenuations (when Cn < 1) or numerical instability (when Cn > 1).

Sectioning for piping systems: method of wave-speed adjustment

In piping systems Δt must equal for all pipes. This involves a certain amount of care in its selection. It is quickly realized that (4) probably cannot be exactly fulfilled in most systems. Inasmuch as the wave speed is probably not known with great accuracy, it may be permissible to adjust it slightly, so that integer N may be found. In equation form this can be expressed as (Wylie and Streeter, 1978):

in which Øj is a permissible variation in the wave speed in pipe j, always less than some maximum limit of say 0.15 or 15% (Wylie and Streeter, 1978). In general, a slight modification in wave speed is more preferable than any alteration in pipe length to satisfy the requirement of a common time step size.

Numerical interpolation

When MOC is applied with Cn < 1 some numerical interpolation must be applied in order to obtain Q and H for every pipe inner section. When the interpolation is applied on the x axis, some analytical expressions can be obtained for the state variables Q and H at interior nodes using numerical schemes with different interpolation orders. The most common numerical interpolation methods include linear interpolation at a fixed time level, including both space line interpolation and reach-out in space interpolation, as well as interpolation at a fixed location, such as time line interpolation or reach-back in time interpolation (Karney and Ghidaoui, 1997). There is a tendency among practitioners to think of interpolation as a numerical device with only numerical side effects. In general, all common interpolation procedures result in numerical dissipation and dispersion, and they considerably distort the original governing equations. The interpolation procedures effectively change the wave speed (Ghidaoui and Karney, 1994). In summary, interpolation fundamentally changes the physical problem and must be viewed as a nontrivial transformation of the governing equations. Because this topic is beyond the scope of this paper, more information will not be included here. In the following paragraphs, the main parameters of the wave speed in (3) will be briefly analyzed, showing their characteristic values.

Compressibility is the property of a fluid to change its volume due to the pressure (Del Valle, 2010). For problems involving the effect of water hammer is necessary to take into account the compressibility of water, which is inversely proportional to its bulk modulus of elasticity and is defined mathematically as:

where v is the specific volume and P is pressure modulus of elasticity K is:

The equation (7) represents the relative change in a fluid volume per unit of applied pressure. The negative sign is because as the pressure increases, the volume decreases and vice versa. The e units are the same for pressure. At a temperature of 20°C and atmospheric pressure (1 bar) the bulk modulus of water is K = 2.07 • 109 Pa. The density of water is the weight of the water per its unit volume:

with ρ the density, m the fluid mass and Vthe fluid volume. The fluid density is function of pressure and temperature (especially in gases), it increases with increasing pressure and it decreases with major temperature. At atmospheric pressure and temperature of 4°C the water density is ρ = 1000 kg/m3. The Young’s elasticity modulus E is the relationship between the force increment and the unitary strain (Martínez and Azuaga, 1997). E has the same value for a tension or compression, being a constant as long as the force does not exceed a maximum value called elastic limit (Hooke’s law). The formula for calculating the elasticity modulus is:

where E is the modulus of elasticity, F is the force, A0 is the surface (area) where the force is applied, ΔL is the length variation and L0 is the initial length. Typical values of E for some materials are shown in Table 2.

 

Table 2: Typical values for E (Larock et al, 2000)

 

When a sample of material is stretched in one direction it tends to get thinner in the other two directions (Figure 1). The Poisson’s ratio is the ratio of the relative contraction strain (or transverse strain) normal to the applied load. It can be expressed as:

 

Figure 1: Contraction strain normal to the applied load

 

where u is the Poisson’s ratio, εt is the transverse strain and εL is the longitudinal or axial strain. Strain can be expressed as:

where dL is the change in length and L is the initial length. For isotropic materials the Poisson’s ratio is in the range of 0 to 0.5 (Greaves et al., 2011). Table 3 shows some typical values of u.

 

Table 3: Typical values for Poisson’s ratio (Larock et al., 2000)

 

Example of application

It is considered a simple system composed by a reservoir (upstream end), a steel pipeline of length L = 4800 m carrying water without presence of free air or gas, and a valve (downstream end). The pipe is anchored against any axial movement (Figure 2). The temperature is 4°C. To calculate a0 we take into account the following values: K = 2.07 GPa, ρ = 1000 kg/m3, E = 207.7 GPa, D = 0.3 m, e = 0.00755 m, Ψ = 0.9531 (see Table 1, case 2). Substituting these values in (3) it is obtained a0 = 1225 m/s.

 

Figure 2: Pipe example sketch

 

Assuming that pipe network was discretized using Δt = 0.3725 s and N = 10, according to (4) we have: Cn = a0Δt/Δx = 1225 • 0.3725 / (4800 / 10) = 0.95065. Because Cn < 1 and the application of the numerical interpolation is not an option, it will be necessary to modify the a0 value in order to get Cn = 1.

Assumptions

In order to obtain Cn = 1, the following assumptions will be taken account: i) it is not possible to modify Δt, L and N; ii) water parameters such as density p or bulk modulus K are known and unalterable; iii) there are not availability of schemes different than MOC as those posed by Twyman et al. (1997), which are more stable and accurate when Cn < 1 and that do not require to modify the wave speed in order to get a more accurate solution.

Wave-speed adjustment

Because of the previous assumptions, in order to obtain Cn = 1, a0 value must be incremented slightly up to 5%, that means up to al = 1289 m/s. Now, the question is what values should adopt parameters u or E in (3) to justify the value of a1, under the scenario that they are the only parameters which can be modified?. For example, to obtain al = 1289 m/s, Poisson’s ratio u must be incremented up to 0.660 (ceteris paribus), see Figure 3. This value for u corresponds to an unknown material and it is out of range because it is greater than 0.5, the isotropic upper limit (Greaves et al., 2011). On the other hand, to obtain al = 1289 m/s, the elasticity modulus E must be incremented up to 3.12 • 103 GPa (ceteris paribus), see Figure 4. This value for E also is out of range because it belongs to a material that cannot be efficiently used in the manufacture of pipes for water distribution systems.

Discussion

Figures 3 and 4 show that to allow the variation of a in 5% (and get Cn = 1), u and E must take unrealistic values that are outside the normal range of physical constituent material of pipes. For example in Figure 3, given the range of extreme values of u between fiberglass reinforced plastic with u = 0.22 and polyethylene with u = 0.46, the nonnumeric physical variation range of a should be between -0.5 and +1.7%, corresponding to a range of variation of u between -25 and 55% of its original value, respectively.

 

Figure 3: Variation of wave speed versus Poisson’s ratio

 

Figure 4: Variation of wave speed versus modulus of elasticity (k/ms2 = Pa)

 

Even though, in the analyzed example u should vary 120% to get a value of 0.66 (out of range) for a = 1289 m/s and therefore Cn = 1. Figure 4 shows that +5% variation of the wave speed to obtain Cn = 1 implies to increase the value of elasticity modulus up to E = 3120 GPa, which corresponds to a variation of 50% of its initial value, being out of range because the existence of a material more rigid than steel and equally efficient and useful as a constituent element of a pipe is unknown. Another interesting point from Figure 4 is that the allowed variation of a falls between 0 and -11%, which only allow incorporate pipes of copper (E = 1100 GPa) or bronze (E = 10000 GPa) into the model because their rigidity. That is, by including less rigid pipes into the model (PVC for example, with E = 2.8 GPa), the variation of a would stay outside the allowable range ±15% recommended by authors such as Wylie and Streeter (1978). The modification of u or E values leads to Cn = 1, assuring the optimum results in numerical terms. Nevertheless, in this case the application of the method of wave-speed adjustment had a cost in physical terms because Cn = 1 was obtained from parameters (E or u) out of range. In general, analysts tend to forget such cost because the MWSA has been recommended in the pipeline literature (Karney and Ghidaoui, 1997), without giving further details about its physical limitations. Finally, another point is that physical limitations of the (4) and (5) show up to where it is possible modify a in order to avoid an out of range value. For example, in the case of u (Figure 3), the range of permitted variation of a is very restricted, between -1 and +1%. For the case of E (Figure 4) such range varies between 0 and -10%. This means that the range of variation of a between -15 and +15% recommended by Wylie and Streeter (1978) is only referential, that is, it shows the maximum range of values to take in case of need to modify the a values, taking care to apply an arbitrary percentage change without checking the numerical aspects and physical constraints that are behind (4) and (5).

Conclusion

The MWSA distorts the physical characteristics of the water hammer problem. Due to this, it is recommendable that in the process of discretization (Δx, Δt) of the pipe network, necessary to solve the water hammer in pipe networks by MOC, before deciding to apply the MWSA to obtain Cn = 1, the analyst must see if the final values adopted to calculate a are consistent and appropriate, both in numerical and physical terms. Otherwise it would solve a very different problem originally raised with implications for all stages of design or verification of the system. Before changing the value of a, it is important to check the implications of changing its magnitude. At this point, it is important to know what parameters of its formulation are known and can be considered as unalterable (pipe length, diameter or wall thickness) and check what of the other parameters can be modified by analyzing its variation range and level of reality.

 

References

Chaudhry, M.H. (1979). Applied hydraulic transients. Van Nostrand Reinhold, New York

Chaudhry M.H. and Hussaini M.Y. (1985). Second-order accurate explicit finite-difference schemes for waterhammer analysis. Journal of Fluid Engineering 107(4), 523-529

Del Valle, V (2010). Fluidos. Apuntes editorial Universidad Tecnológica Nacional, Tucumán, Argentina

Ghidaoui M.S. and Karney B.W. (1994). Equivalent differential equations in fixed-grid characteristics method. Journal of Hydraulic Engineering 120 (10), 1159-1175

Ghidaoui, M.S., Zhao, M., McInnis, D.A. and Axworthy, D.H. (2005). A review of water hammer theory and practice. Applied Mechanics Review 58(1), 49-76

Goldberg, D.E. and Wylie, E.B. (1983). Characteristics method using time-line interpolations. Journal of Hydraulic Engineering 109(5), 670-683

Greaves, G.N., Greer, A.L., Lakes, R.S. and Rouxel, T. (2011). Poisson’s ratio and modern materials. Nature Materials 10(11), 823-837

Karney B.W. and Ghidaoui M.S. (1997). Flexible discretization algorithm for fixed-grid MOC in pipelines. Journal of Hydraulic Engineering 123 (11), 1004-1011

Larock, B.E., Jeppson, R.W. and Watters, G.Z. (2000). Hydraulics of pipeline systems. CRC Press, Boca Raton, Florida, USA

Martínez, P. y Azuaga, M. (1997). Medición del módulo de elasticidad de Young. Apuntes laboratorio IV, Departamento de Física, UBA

Pierre, B. (2009). Pressure waves in pipelines and impulse pumping: physical principles, model development and numerical simulation. Doctoral thesis, Norwegian University of Science and Technology, Trondheim

Skalak, R. (1955). An extension of the theory of water hammer. Tech. Report No. 15, Columbia University

Twyman, J., Twyman, C. y Salgado, R. (1997). Optimización del método de las características para el análisis del golpe de ariete en redes de tuberías. XIII Congreso Chileno de Ingeniería Hidráulica, Universidad de Santiago de Chile, 53-62

Watters, G.Z. (1984). Analysis and control of unsteady flow in pipelines. 2nd edition, Butterworth-Heinemann, USA

Wood, D.J., Lingireddy, S., Boulos, P.F., Karney, B.W. and McPherson, D.L. (2005). Numerical methods for modeling transient flow in distribution systems. Journal of the American Water Works Association 97 (7), 104-115

Wylie, E.B. and Streeter, VL. (1978). Fluid transients. McGraw-Hill, USA

 


Fecha de entrega: 13 de mayo 2016
Fecha de aceptación: 7 de noviembre 2016

 

Creative Commons License Todo el contenido de esta revista, excepto dónde está identificado, está bajo una Licencia Creative Commons