A MIXED MOC / FDM NUMERICAL FORMULATION FOR HYDRAULIC TRANSIENTS

Original scientific paper A mixed numerical formulation, based on the method of characteristics (MOC) and finite difference method (FDM), is proposed for the computational simulations of hydraulic transients. To include the effects of unsteady friction, i.e. the contribution of unsteady behavior of wall shear stress, the numerical algorithm includes the Brunone friction model. The governing water hammer equations are discretized by MOC and the discretization of the local and convective acceleration term in the unsteady friction model is conducted by FDM. The procedure results in an explicit-implicit time marching scheme used for predicting the unknown values of piezometric head and discharge from the known initial conditions. The resulting system of equations is resolved in an iterative fashion. Due to its relatively low complexity, the obtained numerical algorithm is attractive for computational implementation. At the end of the paper a few illustrative numerical examples are presented, compared with available experimental data, and discussed.


Introduction
It is well known that a change in the velocity field will cause pressure disturbance propagation through the space occupied by the fluid.If a pressurized pipeline system is considered, the requirement for flow manipulation at specific locations (e.g.pump or valve) will result in a pressure gradient that induces pressure wave propagation known as water-hammer [1,2].The obligation to study such phenomenon is obvious since the induced pressure disturbance can cause explosion or implosion of pipes.
With the intention to predict the consequence of velocity variations, many authors have proposed different numerical models for friction.Although the developed models work well for relatively slow transients, by comparing the theoretical predictions with the measurements for rapid transients a disagreement in attenuation, shape and timing of pressure traces was evidenced [3,4].The main reason of these disagreements was the inappropriate definition of the fiction factor that in such rapid events should account for unsteady behavior of wall shear stress.To include such effects, several constitutive models were suggested for quantifying the contribution of frictional forces [5,6,7].The implementation of such constitutive models in a onedimensional flow model was explored for different numerical methods such as: finite difference method (FDM), finite volume method, finite element method, wave plan method and method of characteristics (MOC).However, among all the mentioned, MOC is still the most popular method and this is due to its relatively simple implementation, efficiency and accuracy.
In this paper the water-hammer phenomenon is approached by considering the Brunone unsteady friction model [5] and its specific numerical implementation in the equations obtained by MOC.The procedure results in a relatively simple numerical algorithm.

Governing equations
The classical water-hammer theory [8], developed for a one-dimensional approximation of the flow, can be summarized in terms of piezometric head H (defined with respect to the arbitrary datum level as the sum of pipe elevation z and pressure head h) and discharge Q by the momentum equation [8,9] and the continuity equation where A denotes the cross-section flow area, α is the pipe angle inclination with respect to the horizontal level, g is the acceleration of gravity, a the celerity of the generated pressure wave and J the frictional head loss per unit length (defined by some the constitutive relations for quasi-steady and unsteady friction model).The wave speed is given by where D denotes the internal pipe diameter, E w denotes the Young's modulus of elasticity of pipe-wall material, E f is the bulk modulus of elasticity of liquid, ρ is the density of the liquid and s denotes the pipe-wall thickness.

Constitutive equations
The frictional term J can be defined by the additive decomposition as J q +J u , where the first term represents the contribution of viscous forces related to the quasisteady flow condition and the term J u represents the contribution of unsteady behavior of wall shear stress.

Quasi-steady friction model
The frictional head loss J q associated to quasi-steady flow conditions is quantified according to the well-known Darcy-Weisbach equation where λ denotes the Darcy-Weisbach friction coefficient.In the present numerical formulation λ is obtained from an explicit approximation of the Colebrook-White equation.Accordingly, if the Reynolds number Re defines a laminar flow regime, λ is explicitly obtained as 64/Re and in a turbulent flow regime as [10] , 7 3 , , e λ (5) where ε denotes the pipe wall roughness.

Unsteady friction model
To include the effect of unsteady behavior of wall shear stress, the Brunone friction model is used [5].In this case the friction head loss J u is dependent on instantaneous mean flow velocity and instantaneous flow acceleration (temporal acceleration).The original form of the model is [5] or according to Vitkovsky's modified formulation [6], which additionally takes into account the correct sign of the convective term, λ is defined by the formula [6,11] in which Φ denotes an integer value and k denotes the Brunone friction coefficient The Brunone coefficient in ( 9) is usually defined by using Vardy and Brown's [6]

Numerical formulation
To obtain a computational algorithm that can be used to perform predictions of practical importance, Eq. ( 1) and Eq. ( 2) are transformed according to MOC [9].

Method of characteristics
Since the governing Eqs. ( 1) and ( 2) are of hyperbolic type, the Lagrange multiplier method can be used to obtain their linear combination that leads to ordinary differential equations along two intersecting families of curves in the x-t coordinate plane.By introducing the Lagrange multiplier L, the procedure leads to By invoking the definition of material-time derivative [8], which for an arbitrary quantity q and a velocity field v is defined as D(q)/Dt=(∂q/∂t) + v(∂q/∂x), it follows that the equality (11) will be true only if L=±1/a.The conclusion can be used to rewrite Eq. ( 11) in terms of one equation for a positive characteristic curve [8] along the line with the slope dx/dt=v+a, and one equation for a negative characteristic curve which defines the solution in terms of H and Q on a line with a slope dx/dt=v−a.

Domain discretization
To obtain a numerical approximation of Eq. ( 12) and Eq. ( 13), the involved material-time derivatives DH/Dt and DQ/Dt are approximated at discrete spatial and temporal positions.For this purpose the spatial domain is discretized by a finite number nx of equidistant spatial increments ∆x and the temporal domain is substituted by nt equidistant time steps ∆t.The discrete spatial coordinate i is introduced as an index to denote the spatial position of relevant quantities.On the other hand, the exponent n is introduced to denote the temporal position of quantities.The time step is defined as ∆x/a.According to the introduced notations for spatial and temporal positions, the positive characteristic C + (12) can be approximated as [8] 0 1 1 and the negative characteristic C -(13) as To close the numerical formulation, the friction term J must be also represented in a discrete manner.For this purpose the quasi-steady and unsteady contribution will be separately defined.
Except the discharge Q and head H, the Darcy friction coefficient λ in the quasi-steady friction term J q (4) should be also computed at each computational section i and at every time step n.This requires the computation of Re i .Although this procedure should be iterative, because λ i depends on the yet unknown velocity v i , it is usually calculated from the known velocity at the previous time step.The procedure for computing λ i in a time step n from the velocity v i in the same time step will require an iterative loop.Apart from the coefficient λ i , the discharge Q in (4) is specified depending on the line C + or C - [8,9].
Differently from the quasi-steady friction term, the computation of J u is slightly more complex.The difficulty arises from the unsteady term ∂Q/∂t in (7).Namely, the local acceleration is defined with respect to the change in velocity between time steps n and n+1.To approximate the derivative ∂Q/∂t, for both equations ( 12) and ( 13), in the present numerical formulation the forward difference approximation is used.On the other hand, the convective acceleration ∂Q/∂x is discretized by the central difference approximation.Accordingly, J i in the positive characteristic C + can be defined as and for J i in the negative characteristic C -as Since in this manner the computation of the unsteady term ∂Q/∂t requires the discharge Q at the unknown time step n+1, an iterative algorithm is required to obtain H and Q at the same time step.Before illustrating the structure of the algorithm, the algebraic approximation for H and Q should be defined.The equation for H i at a time step n+1 can be obtained by multiplying Eq. ( 14) and Eq. ( 15) with gA∆t and then by combining the obtained results.In this case H i at a time level n+1 takes the form In a similar manner, the multiplications of Eq. ( 14) and Eq. ( 15) with a∆t will lead to With the given boundary and initial conditions, the variables H i and Q i inside the flow domain, i.e. between the coordinates i=2 and i=nx-1, can be obtained from the known initial conditions at the time level n.Due to the introduced approximation for the unsteady term in Eq. ( 7), the procedure is here defined in an iterative fashion.The criterion for convergence is simply defined by a tolerance value tol used to check the difference between discharge values Q computed in successive iteration steps.By denoting the initially assumed discharge value by , the structure of the numerical algorithm can be illustrated with the next pseudo-code (algorithm 1).algorithm 1 If the starting value for Q is assumed to be the value from the previous time step ( in algorithm 1), the numerical results show that the discharge value in a time step n+1 is reached within a few iterations.

Boundary condition at the reservoir
If one side of the pipe is attached to the reservoir, the specification of the boundary condition at the same location is usually very simple.Namely, on this case the standing water in the reservoir acts as a «wall» from which the pressure wave reflects.However, it should be noted that this reflection is not perfect and that certain wave energy is transferred into the reservoir.
By assuming that the water level in the reservoir does not change during the time of water-hammer propagation, which is reasonable since the event usually lasts a relatively short period of time, in every time step ∆t the piezometric head H at the same computational section i can be assumed equal to H 0 , i.e. equal to the water level in the reservoir.On the other hand, the discharge Q at the same location follows from the negative characteristic C - (15).However, due to the unsteady friction term, the procedure is conducted through iterations.Also, for the upstream boundary condition the convective acceleration is approximated with the forward difference so that J -for i=1 takes the form By assuming that the upstream piezometric head is defined by H 0 , the resulting numerical algorithm is similar to the previous one and is illustrated in the following pseudo code (algorithm 2).algorithm 2 J according to Eq. ( 16)

Boundary condition at the valve
The boundary condition at the valve, here assumed to be at i=nx, needs the specification of a function β(t) which defines the ratio between the current flow area A n and the open cross-section area A of the valve.However, to implement the boundary condition in the numerical scheme, the predefined function β(t) is discretized in equidistant time steps ∆t.
The discharge at the valve is related to the steady state head loss H at the same location at time t=t 0 .By including the discharge coefficient as c c (β), the discharge Q at i=nx for an open valve at the beginning of the computation is defined by the well known equation ( ) For a smaller flow area at the valve, the same principle as (21) can be applied by including the reduced flow area defined as Aβ n+1 .Accordingly Eq. ( 21) can be rewritten in the form ( ) To retrieve the equation for the discharge Q at the time level n+1, it is wise to define the dimensionless valve opening parameter as a ratio Dividing Eq. ( 22) by Eq. ( 21) the discharge Q at the valve obtains the form By appending Eq. ( 14) to the previous one, the real root of the obtained quadratic equation is given by [8] The piezometric head H at the valve at the time level n+1 can be computed by The friction head loss contribution J + in Eq. ( 27) is calculated by replacing the convective acceleration term in the unsteady friction model by its backward difference approximation (29). 1 The related numerical procedure is summarized in the next pseudo code (algorithm 3).However, note that the procedure is not valid if the pressure head at the valve is zero as it can be in the case in which the discharge is under atmospheric pressure and the minor loss across the valve is negligible (H nx =0).In this case the computation should involve the velocity head in the computation.algorithm 3 Obviously, to implement the boundary condition at the valve, the function β(t) should be defined in advance.Moreover, by specifying different functions β(t) the influence of the time dependent flow area reduction at the valve on the generated pressure variation in the pipe system can be predicted and consequently studied.To illustrate the procedure, the following numerical examples will be conducted for three types of function β(t).

Initial conditions
The initial conditions are defined by the piezometric head H i and discharge Q i in each computational section i at the beginning of the computation (n=0).Regardless of the considered pressurized pipeline system, the distribution of H and Q at the beginning of the computation should define a real flow situation.

Numerical examples
The numerical algorithm is used to compute H and Q at the valve and at the midpoint of a simple pressurized pipeline system.The considered pressurized system is made of two reservoirs which are connected by a pipe with a valve at one end.The pipe length is 37.2 m, the internal diameter is 22 mm and the wall thickness is 1.6 mm.The pipe material is copper with E p =100 GPa.The head H in one tank is 32 m and in the other one it is 31.85 m.By including the minor head loss at the entrance, at the exit of the pipe and at the valve, for a steady flow condition with the valve completely open, the given water levels in the tanks define the flow velocity of 0,2 m/s.With the defined pipe and water characteristics the celerity is obtained to be a=1319 m/s.The described pipe system defines the laboratory apparatus used to study water-hammer phenomena at the Robin hydraulics laboratory at the University of Adelaide (for more details see [6]).The valve closure time is set to t c =0,009 s.Within this time interval three different scenarios are numerically simulated: sudden, linear and bilinear reductions of the flow area in the valve.However, the experimental data are given for a sudden valve closing.
The number of computational sections i for all numerical examples is set to nx=100.With the computed celerity a of 1319 m/s (3), the adopted time step is ∆t=0,00028 s.The considered temporal domain is 1 s and is covered by nt=3546 time steps.

Sudden reduction of flow area
The almost sudden reduction of the flow area at the valve starts at t 0 =0,1 second in the time domain.The reason for that is to show that the numerical model can handle steady flow conditions without alteration of the steady flow variable H and Q in time, i.e. during the period in which there are no changes in the flow area.Namely, between t=0 and t=0,1 s the flow conditions are the same as defined by the initial conditions.
For each numerical example the results are given by three diagrams (Fig. 1).The top diagram defines the dynamics of the relative flow area reduction at the valve, i.e. the function β(t).Below this diagram the temporal variation of piezometric head H at the valve and at the midpoint of the pipe are presented.At the bottom, the discharges Q at the valve and at the midpoint of the pipe are given.For all diagrams the time axis t is the same.The results of the numerical example are in good accordance with the available experimental data [6].Apart from the obvious qualitative agreement, the comparison of quantities, such as the reduction of maximal pressure values in time, is also in accordance with the experiments.

Linear reduction of flow area
The second numerical example considers the linear reduction of the flow area.In this case the reduction of the flow area starts at t=0 and ends at the same time t=0,1 s.With respect to the previous numerical example (Fig. 1), in this case it can be expected that the maximal computed pressure will be reduced because 0,1 s > 2L/a.The results are given in Fig. 2 and accordantly it can be seen that the maximal piezometric head H is about 40 m.

Bilinear reduction of flow area
The last numerical example is performed for the usually assumed function β in practical cases.However, the time in which the reduction takes place is the same as in the previous examples, i.e. 0,1 s.This time period is here assumed to enable the direct comparison between the performed numerical analyses.The results are summarized in Fig. 3.
It can be noted that the maximal obtained piezometric value H is slightly greater that the one obtained in the previous example (Fig. 2).This is because the velocity decay is very rapid in the beginning of the flow area reduction process.Namely, as can be noted from Fig. 3c, the reduction to zero of discharge Q at the valve is almost the same as the function β(t).The maximum peak in H(t) at the valve occurs within this time period.

Conclusions
A mixed numerical formulation, which is based on the MOC and FDM, is developed for the computational simulation of hydraulic transients.The mixed numerical approach arises from the fact that the governing equations for water hammer are discretized by MOC and the unsteady friction model is discretized by FDM (Brunone friction model).The temporal acceleration is discretized by the forward difference approximation and the spatial acceleration by the central difference approximation.The resulting system of equations is solved in an iterative fashion.The performed numerical examples show that the proposed numerical model is able to predict the temporal variations of H(t) and Q(t) retrieved from experimental data [6].Due to its relatively low complexity, the numerical algorithm is attractive for computational implementation.Moreover, the presented numerical approach is sufficiently general that it can be easily exploited for analysis of complex pipe systems.

Figure 1
Figure 1 Results for the case of sudden reduction of the flow area at the valve: a) valve function β(t), b) temporal variations of H at the valve and at the pipe midpoint and c) temporal variations of Q at the valve and at the pipe midpoint.

Figure 2
Figure 2 Results for the case of linear reduction of the flow area at the valve: a) valve function β(t), b) temporal variations of H at the valve and at the pipe midpoint and c) temporal variations of Q at the valve and at the pipe midpoint.

Figure 3
Figure 3 Results for the case of bilinear reduction of the flow area at the valve: a) valve function β(t), b) temporal variations of H at the valve and at the pipe midpoint and c) temporal variations of Q at the valve and at the pipe midpoint.