Numerical Hydraulics

1. Numerical Hydraulics Lecture 3: Computation of pressure surges W. Kinzelbach with Marc Wolf and Cornel Beffa

2. The phenomenon Valve Reservoir Till time t=0: Steady state flow Q At time t = 0: Instantaneous closing of valve Observation: Sudden pressure rise at valve →pressure surge (water hammer)

3. And some typical damages Sayano-Shushenskaya plant in southern Siberia Water pipe damage due to pressure surge San Bruno: PG&E Power Outage and Pressure Surge Preceded Blast

4. Where pressure surges occur… • water distribution systems • waste water transfers • storm water rising mains • power station cooling systems • oil pipelines • RAS (activated sludge) onsite pipelines • hydropower stations • and any fluid system in which the inertia (mass and velocity) of the fluid is significant.

6. The phenomenon Pressure wave propagates with wave velocity c If valve closing time is smaller than the run time of the wave to the reflection point and back the surge is called Joukowski surge Pressure vs. time at valve Damping of amplitude through friction

7. Negative pressure wave The negative pressure wave can not become lower than the vapour pressure of the fluid. If the pressure falls below the vapour pressure, a vapour bubble is formed. The water column separates from the valve. When the pressure increases again the bubble collapses. This phenomenon is called cavitation.

8. Pressure surge with vapour bubble formation

9. Computed and observed pressure surge at two places along a pipe Today the reliable computation of pressure surges is possible

10. Measures against pressure surges Surge shaft

11. Measures against pressure surges Surge vessels (Windkessel) Special valves

12. The equations of unsteady pipe flow • Continuity • Momentum equation (per unit volume) (a inclination angle of pipe)

13. Further transformations (1) • Continuity: As density and cross-sectional area depend on x and t only via the pressure p, the chain rule can be applied. Using the moduli of elasticity of water EW und of the pipe Epipe e is the pipe wall thickness, E‘ is the combined modulus of elasticity of the system

14. Further transformations (2) • Compressibility: Definition • Using pressure tank formula

15. Further transformations(3) • Momentum equation

16. The equations of unsteady pipe flow • Continuity • Momentum equation (1) (2) 2 PDE with 2 unknown functions p(x,t) und v(x,t) plus initial and boundary conditions

17. Boundary conditions • Pressure boundary condition: p given • e.g. water level in reservoir, controlled pump • Velocity/Flux boundary condition: v given • e.g. flow controlled (v from Q/A) • Combination: Relation between pressure and flux given • Z. B. function of pressure reduction valve, characteristic curve of pump • Closing of a valve at the end of a pipe • Initially flow Q, then according to closing function reduction to zero withing closing time of the valve.

18. Linearised equations • Delete all terms in (1) and (2) which are non-linear (for convenience: a = 0): • General solution by elimination: • Take partial derivative of first equation with respect to t • Take partial derivative of second equation with respect to x • Subtraction yields:

19. Linearised equations • Wave equation (for p, analogously for v) which has general solution • Wave with wave velocity • Example: Modulus of elasticity of steel = 200‘000 MN/m2, Modulus of elasticity of water = 2‘000 MN/m2, wall thickness e = 0.02 m, D = 1 m, r = 1000 kg/m3 yields c = 1333 m/s

20. Joukowski surge • Estimate of surge pressure after instantaneous closing of valve (neglecting friction, linearized equations): „Worst case“ • General solution: Proof by insertion into linearised equations!!

21. Joukowski surge • After t = 0 only the backward running wave F(x+ct) is found in the upstream • v at the valve is 0 • Maximum Dp is given by: • Solution: Example continued: c=1333 m/s, Q0=1 m3/s, L=100 m yields: Dp=1.7E6 N/m2

22. Numerical solution of the complete equations • Difference method • Discretisation of space and time • Dx and Dt • Difference equations for time step t, t+Dt • Problem: Discretisation „softens“ pressure front numerically • Way out: Method of characteristics • Follows the pressure signal in moving coordinate system

23. Method of characteristics • Normal difference method • Softening of pressure front • Method of characteristics • Grid is adapted to frontal velocity (feasible, as v<<c, c+v ≈ c-v ≈ c) Front of pressure wave cDt < Dx Dx Front of pressure wave cDt = Dx Dx

24. Method of characteristics • Replacing equations (1) and (2) by 2 linear combinations yields:

25. Method of characteristics • With total derivative along x(t) the equations have the form: Forward characteristic Backward characteristic (c is actually relative wave velocity with respect to average water movement.)

26. Difference scheme Divide pipe of length L in N sections, length of one section Dx = L/N Node N+1 Node 1 section N section 1 x Upper index time step, lower index node Chose time step such that Dx = c Dt In every time step there are two unknowns at each of the N+1 nodes: To determine these unknowns 2N+2 equations are required. From quantities at time j quantities at time j+1 are computed. The new times j+1 become the old times j of the next time step.

27. 3 3 2 1 2 1 Difference scheme time Using c + v ≈ c - v ≈ c node i communicates within time interval Dt with node i-1 via the forward characteristic and with node i+1via the backward characteristic j +1 Dt j Dx = c Dt space i-1 i i+1 Dx Total derivative or derivative along characteristic line

28. Difference form of equations • Equations for nodes 2 to N: 2N-2 equations forward characteristic backward characteristic The pressure loss term is linearised by evaluating it at the old time j Equations can be solved for The further equations are determined by the boundary conditions and the one characteristic which can be used at the respective boundary

29. Method of characteristics • Example: Reservoir with pipe which is closed instantaneously at t=0 • 2 further equations from boundary conditions In the example: • 2 further equations from characteristic equations In the example: From forward characteristic for i=N+1 From backward characteristic for i=1

30. forward characteristic = 0 backward characteristic = 0 Simplified case for basic Matlab-Program a=0, friction neglected, equations nodes 2 to N Solution by subtracting resp. adding the two equations

31. Simplified case for basic Matlab-Program 2 equations from boundary conditions 2 equations from characteristics for i = 1 and i = N+1 From forward characteristic for i=N+1 = 0 From backward characteristic for i=1 = 0

32. Additions • Formation of vapour bubble • Branching pipes • Closing functions • Pumps and pressure reduction valves • …. • Consistent initial conditions through steady state computation of flow/pressure distribution

33. Example (1) Use Program „Hydraulic System“ Tank 2 Tank 1 connecting pipe valve L=500 m D = 0.2 m, e = 0.01 m, roughness k = 0.0005 m r = 1000 kg/m3, Ew = 2000 MN/m2, Epiper = 210000 MN/m2 pressure downstream reservoir 80 mWS, pressure upstream reservoir 90 mWS closing time of valve1 s, Q before closing: 0.2 m3/s loss coefficient valve 2, time of calculation 60 s, number of pipe sections n = 10

34. Example (2) Tank 3 L=500 m Tank 2 Tank 1 Valve L=500 m D = 0.2 m, e = 0.01 m, roughness k = 0.0005 m r = 1000 kg/m3, Ew = 2000 MN/m2, Epipe = 210000 MN/m2 pressure of both downstream reservoirs 80 mWS, pressure upstream reservoir 90 mWS closing time 1 s, Q before closing of valve: 0.2 m3/s loss coefficient of valve 2, computation time 60 s, number of pipe sections n = 10