Water Hammer Simulations by S. Mambretti

Water Hammer Simulations

WITPRESS

WIT Press publishes leading books in Science and Technology. Visit our website for the current list of titles. www.witpress.com

WITeLibrary

Home of the Transactions of the Wessex Institute, the WIT electronic-library provides the international scientific community with immediate and permanent access to individual papers presented at WIT conferences. Visit the WIT eLibrary athttp://library.witpress.com

This page intentionally left blank

Water Hammer Simulations

S. Mambretti

Universidade Estadual de Campinas, Brazil

S. Mambretti Universidade Estadual de Campinas, Brazil

Published by WIT Press Ashurst Lodge, Ashurst, Southampton, SO40 7AA, UK Tel: 44 (0) 238 029 3223; Fax: 44 (0) 238 029 2853 E-Mail: [email protected] http://www.witpress.com For USA, Canada and Mexico WIT Press 25 Bridge Street, Billerica, MA 01821, USA Tel: 978 667 5841; Fax: 978 667 7582 E-Mail: [email protected] http://www.witpress.com British Library Cataloguing-in-Publication Data A Catalogue record for this book is available from the British Library ISBN: 978-1-84564-680-6 eISBN: 978-1-84564-681-3 Library of Congress Catalog Card Number: 2013938170 No responsibility is assumed by the Publisher, the Editors and Authors for any injury and/or damage to persons or property as a matter of products liability, negligence or otherwise, or from any use or operation of any methods, products, instructions or ideas contained in the material herein. The Publisher does not necessarily endorse the ideas held, or views expressed by the Editors or Authors of the material contained in its publications. © WIT Press 2014 Printed by Lightning Source, UK. All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the Publisher.

Preface

Waterhammer is a pressure surge or wave caused when a fluid (usually a liquid but sometimes also a gas) in motion is forced to change its velocity. It commonly occurs when a valve closes or opens at the end of a pipeline system, or a pump starts or stops; as a consequence, a pressure wave propagates in the pipe. This pressure wave can cause problems, and normally the efforts of the designers are aimed at reducing its effects. However, in some cases the pressure pulses are deliberately caused in order to achieve particular results, such as the pumping of a fluid or the mapping of a network. In all cases, the study of the consequences of this phenomenon is required, and this might be performed with rough formulas as those developed between the end of 19th Century and the beginning of the 20th, or through more advanced models. Nevertheless, all the models have limitations and might be criticized. During the years I spent as a professional consultant and researcher in the field, the need for a reference able to guide through the several models that have been developed and the analysis of their results clearly emerged. Hence the reason to write this book. In the book, chapter 2 presents the governing equations and the hypotheses involved; chapter 3 shows the simplified solutions that have been developed and used before the arrival of computers; in chapters 4 and 5 two numerical integration methods are shown. In chapter 6 a number of devices are presented with the methods that can be used for their modeling. In chapter 7, problems related to the implementation of the model are discussed, giving some suggestions for their solution. In chapter 8 an important phenomenon, the presence of air and cavitation, is analyzed, while in chapter 9 the most advanced models are presented and discussed. Chapter 10 presents some

actual hydraulic plants whose behavior has been analyzed and studied with the methods and models presented in the book. Writing a book is always a challenging task, and many people contributed in some way to making this volume a reality. All the case studies have been carried out with the invaluable help of my colleague and friend Dr. Paola Pianta, with whom I discussed many developments and applications of the models presented in this book. Most of my research has been carried out at Politecnico di Milano, where all the experimental tests have been performed, thanks to Prof. Enrico Orsi, head of the Laboratory of Hydraulics, and Prof. Enrico Larcan, head of the Department of Hydraulics. After moving to Brazil, invited by Prof. Paulo Barbosa, head of the Faculty of Civil Engineering, I found an excellent welcome and all the freedom and help I needed to finish the book, thanks to the kindness of Prof. José Geraldo Pena de Andrade, head of the School of Technology, University of Campinas. I cannot forget Prof. Daniele De Wrachien, who has always been a support and a help in my academic life, and Prof. Carlos Brebbia, who not only is the Director of WIT, but he is also an example to be followed both from the scientific and the human point of view. Finally, I have to thank my wife Grazia, both because she allowed me to work overtime and because she drew all the figures for the book. All the mistakes that can be found in the book have to be attributed to my limitations. I would be grateful if the Readers would report any errors and suggest any improvements for the subsequent editions of the book. Stefano Mambretti São Paulo, 2013

Contents

1 An old topic still not completely solved

2

3

1

1.1

An old topic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2

Applications and problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.3

This book . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.4

The computer programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

Compressible flow theory: basic concepts

9

2.1

Instantaneous operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2

Wave celerity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3

Velocity of operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4

Non-negligible headlosses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.5

Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . 17 2.5.1

Continuity equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.5.2

Momentum equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.5.3

The governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Simplified solutions

25

3.1

The Allievi’s method (1913) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.2

The non-elastic hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.3

3.2.1

Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.2.2

Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

Graphical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4

5

6

Numerical solution of the governing equations: The method of characteristics

41

4.1

Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4.2

Initial and boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.2.1

Reservoir . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.2.2

Valve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.2.3

Junction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . 47

4.3

The computer code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.4

First simple application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Numerical solution of the governing equations: finite difference methods

55

5.1

The Courant–Friedrichs–Levy stability condition . . . . . . . . . . . . . . . . . . 57

5.2

The Lax–Wendroff method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.3

Solving the governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.4

Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.4.1

Asymmetrical schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.4.2

Ghost cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.5

The computer code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

5.6

Again the simple application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

Devices – Boundary conditions 6.1

71

Surge tanks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.1.1

Simple surge tanks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

6.1.2

Different types of surge tanks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

6.2

Air chambers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

6.3

Relief valves and rupture disks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

6.4

Centrifugal pumps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

6.5

Other methods for controlling the pressures . . . . . . . . . . . . . . . . . . . . . . . 84

6.6

7

6.6.1

The simple surge tank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

6.6.2

The simple air chamber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

6.6.3

Air chamber with headlosses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

6.6.4

Air chamber and valve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

6.6.5

Valve modelling: an example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

6.6.6

Pumps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

Instabilities 7.1

7.2

7.3

8

The computer codes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

101

Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.1.1

General remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

7.1.2

Computer program for oscillating velocity . . . . . . . . . . . . . . . . 103

Transfer matrix method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.2.1

General remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

7.2.2

Application to hydraulic transients . . . . . . . . . . . . . . . . . . . . . . . 106

7.2.3

Description of simple system: pipes . . . . . . . . . . . . . . . . .. . . . . 108

7.2.4

Description of simple system: valves and effluxes . . . . . . . . . 109

7.2.5

Global matrix of a system . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . 111

7.2.6

A simple application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

Numerical instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7.3.1

Changing CFL number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

7.3.2

First order methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

7.3.3

Flux-limiters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

7.3.4

Artificial dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . 122

Effects of air and cavitation 8.1

127

Cavities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 8.1.1

Formation of the cavities . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . 127

8.1.2

Collapse of the cavities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

8.1.3

Description of the motion in the presence of cavities . .. . . . . 130

8.2

Changing of celerity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

8.3

Water column separation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

8.4

Additional resistance terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

8.5

8.6

8.4.1

Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

8.4.2

Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

Laboratory experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 8.5.1

Experimental set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

8.5.2

Experimental tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

Computer code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

9 Advanced models 9.1

9.2

9.3

10

145

2D models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 9.1.1

Continuity equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146

9.1.2

Momentum equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146

Headlossess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 9.2.1

Pezzinga’s model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . 147

9.2.2

k–ε model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . 149

Cavitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . 150 9.3.1

Release gaseous cavitation model . . . . . . . . . . . . . . . . .. . . . . 150

9.3.2

Second viscosity cavitation model . . . . . . . . . . . . . . . . . . . . . 152

9.4

Numerical schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

9.5

Further problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

Case studies

157

10.1 Simple pressure pipe for petroleum products in Djibouti . . . . . . . . . 157 10.1.1 Plant characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 10.1.2 Expected scenarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 10.1.3 Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 10.1.4 Case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 10.1.5 Case 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 10.1.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161

10.2 A more complex example for seawater treatment plant in Tanzania . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 10.2.1 Plant characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 10.2.2 Unsteady flow simulations: existing plant . . . . . . . . . .. . . . . 162 10.2.3 Plant to be designed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 10.2.4 No air chambers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 10.2.5 Air chamber 3 m3 volume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 10.2.6 Air chamber 5 m3 volume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 10.2.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 10.3 A very complex example for seawater treatment plant in Algeria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 10.3.1 The plant to be modelled . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . 169 10.3.2 A peculiar device: energy recovery PX . . . . . . . . . . . . . . . . . 172 10.3.3 Laboratory plant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 10.3.4 Laboratory tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . 173 10.3.5 Model of the seawater plant in Algeria . . . . . . . . . . . . . . . . . 175 10.3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 10.4 Final remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 References

181


Chapter 1 An old topic still not completely solved

1.1 An old topic The high pressures generated by a sudden variation of the velocity in closed pipes have been studied for centuries. Italian books of the eighteenth century describe the problem of “percossa idraulica” (hydraulic stroke) even if without being able to quantify its value. The pressures generated gave the idea that they might be used to produce work: in 1772 John Whitehurst, in the United Kingdom, invented a manually controlled precursor of the hydraulic ram called the “pulsation engine”, which in its first installation raised water to a height of 4.9 m. The inventor did not patent his idea and the details are obscure. The first self-acting ram pump was invented by the Frenchman Joseph Michel Montgolfier (best known as a co-inventor of the hot air balloon) in 1796. The principle is quite simple, although ingenious; a simplified hydraulic ram is shown in Figure 1.1. The device has to be positioned where fluent water is available, as the ram will use its energy. Initially, the waste valve (3) is open, and the delivery valve (4) is closed. The water in the drive pipe (1) starts to flow under the force of gravity and picks up speed and kinetic energy until the increasing drag force closes the waste valve. The momentum of the water flow in the supply pipe against the now closed waste valve causes a water hammer that raises the pressure in the pump, opens the delivery valve (4), and forces some water to flow into the delivery pipe (2). Because this water is being forced uphill through the delivery pipe farther than it is falling downhill from the source, the flow slows; when the flow reverses, the delivery check valve closes. If the water flow stops, the loaded waste valve reopens against the now static head, which allows the process to begin again. A pressure vessel (5) containing air cushions the hydraulic pressure shock when the waste valve closes, and it also improves the pumping efficiency by allowing a more constant flow through the delivery pipe. Although, in theory, the pump could work without it, the efficiency would drop drastically and the pump would be subject to extraordinary stresses that could shorten its life considerably. One problem is that the pressurized air will gradually dissolve into the water until none remains. One solution to this problem is to have the air separated from the water by an elastic diaphragm (similar to an expansion tank); however, this solution can be problematic in developing

2 WATER HAMMER SIMULATIONS

Air

(5)

(2) Outlet-delivery pipe

Weight (3)

(1)

Waste valve v⬘

(4) Delivery check valve

Waste valve

A

Figure 1.1: Hydraulic ram.

Countries where replacements are difficult to procure. Another solution is to have a mechanism such as a snifting valve that automatically inserts a small bubble of air with each pump cycle. Yet another solution is to insert an inner tube of a car or bicycle tyre into the pressure vessel with some air in it and the valve closed. This tube is in effect the same as the diaphragm, but it is implemented with more widely available materials. The air in the tube cushions the shock of the water the same as the air in other configurations does. The efficiency of this device is obviously quite low, but it has the advantage to be independent from any external source of energy, using only that of the inflowing water, and for this reason it is still used in small plants installed to provide water in villages in the developing countries. The hydraulic ram was developed on the basis of empirical observations, as analytical evaluations of the pressures during transients have been carried out at the end of the nineteenth century. The first equations which provided a quantitative assessment of the pressures were produced by Menabrea in 1858 [1] and Michaud in 1878 [2]; afterwards, Joukowsky in 1898 [3], and, independently, Allievi in 1903 [4] completed Michaud’s work, correcting his results and developing a more comprehensive theory. The paper of Joukowsky was the first that introduced the term “waterhammer”, which later on is used worldwide to refer only to the elastic model, while the phenomenon related to the rigid column model, i.e., mass oscillations, took the name “surge”. Both the Italian “Colpo d’Ariete” and French terms “Coup de Bélier” mean literally “Stroke of a Ram”, and probably they have beed derived from the above-described device “hydraulic ram”. These theories allowed to compute only the maximum value obtained at the beginning of the transients; thenAllievi in 1913 [5] was able to compute the different phases, under simplifying assumptions, with the so-called chained equations.

AN OLD TOPIC STILL NOT COMPLETELY SOLVED

3

Evangelisti, in the sixties of twentieth century [6, 7], developed a numerical method based on the characteristics, which is still used worldwide and implemented in computer programs. In the thirties of the same century graphical methods appeared [8, 9] that allowed, still under simplifying hypotheses, to compute the value of the pressures during the development of the transient. Since then, a really large number of books and papers have appeared on the scientific literature, and among them the works of Streeter and Wylie [10, 11] and Chaudry [12] have to be highlighted. Moreover, many computer programs are already available to study this phenomenon.

1.2 Applications and problems The phenomenon is typically destructive, but some authors [13–16] tried to use the impulsive wave generated during transients to perform network analysis studying the subsequent reflections and in particular to detect the position of the leakages in the water supply systems. As the pressure impulse is reflected by any variation in diameter, material or in presence of a leakage, the analysis of the response to that impulse can show the dissimilarities between the expected and the real networks. As mere evaluation of the position of the leak is not sufficient to evaluate the efficiency of the system, this method has also been developed in order to allow the determination of the amount of the leakage. These authors agree on the following: • An higher leakage implies a greater decrease of the pressure; • The variation of the pressure does not significantly depend on the shape of the orifice or on the outflow conditions; • The friction plays a significant role on the deformation of the pressure wave. In Figure 1.2, a very simple network is shown, where water is provided from reservoir A and discharged to the reservoirs B and C; the leakage is positioned 25 m upstream the reservoir B (position L) and the discharge is abruptly interrupted at the valve V, where the pressures are also recorded during the transients. Actually, the analysis has been performed only with numerical methods and not on physical models or real cases. The comparison between the responses with and without the leakage is reported in Figure 1.3. As the celerity of the wave has been imposed equal to 1000 m/s, the first reflection is recorded after 0.18 s because of the junction J; for the second reflection, the two cases with and without leakage differ: if there are no leakages, the reflection is recorded after 0.28 s while, with the leakage, it is recorded after 0.23 s. The subsequent oscillations could be followed, but they become confuse and lose significance. Moreover, the recorded reflection at 0.23 s means that at a distance equal to 115 m from the valve there is a dissimilarity from the expected network,


70

m

C

A

B

V J

L 25 m

90 m

25 m

Figure 1.2: Simple network: losses are positioned 25 m upstream the reservoir B. 625

Head (m asl)

620

No leakage Leakage

615 610 605 600 595 1.8

2

2.2

2.4 Time (s)

2.6

2.8

3

Figure 1.3: Response of the simple network – comparison between the responses with and without leakages. but, for instance, this piece of information alone does not allow to detect if this dissimilarity is located on the pipe J–C or on the pipe J–B. In depth analysis is required for a better understanding of the actual system. This analysis is surely valid and effective in the localization and assessment of small leaks, but it is sometimes difficult to apply in real cases since it is difficult to estimate the initial conditions and the values of some parameters of the model. Moreover, the reverse procedure that is used to evaluate the extent of the leakages is often burdensome. Finally, this is surely an interesting field of research in the years to come, especially considering the importance of reducing the losses from the water supply networks, both economically and ecologically. However, the methods of modelling the networks and analysing its response need a great improvement. This is not the only field where this phenomenon needs further research: in fact, despite the efforts made to explain the phenomena involved in the transients, many aspects are still to be clarified, and some of them will be discussed in this


5

book. For instance, another important problem still unresolved is represented by headlosses: in the simplified models that have been developed during the twentieth century these have often been neglected; this assumption is surely on the safe side, but does not allow to follow the evolution of the different phases during the transient. Even more complex models available nowadays normally evaluate the headlosses with the expressions carried out for steady flow; as it will be shown in the appropriate chapters of this book, the use of steady flow models does not allow a good description of the unsteady flow phenomena. Therefore, many authors dedicated their time in studying the best methods to improve the models; however, when air is mixed with water, and especially when cavitation occurs, the behaviour of the mixture is very difficult to predict. Finally, as it has been mentioned that the phenomenon is destructive, methods and devices to reduce its effects have been studied and modelled. For instance, surge tanks are very large structures designed since the end of nineteenth century to protect the first hydropower plants, and their dimension together with their costs and their strategic importance justified a very accurate design. However, as will be shown in this book, it is very difficult to accurately model the behaviour of even the simplest devices, and some of them can even produce instabilities that worsen the condition of the plant where they are installed. Pumps, on the other hand, are particular devices that generate unsteady flow and are affected by the additional pressures obtained during transients; their behaviour during transients is very difficult to be reproduced. Again, for the interested researchers there is a lot of work to do.

1.3 This book As the literature on this topic is very extensive, before starting this job the question was whether a new book was really newsworthy and what it would have added to the general knowledge in this area. The main target of this book is providing the methods for modelling transients in closed pipes to the reader. When is a computer model, i.e., a deep analysis necessary? ASCE [17, 18] developed guidelines to identify pipelines at risk; however, Ellis [19] correctly stated that caution is exercised in all cases where some form of study is carried out, even if not a detailed computer analysis. In the first chapters the old formulae are shown, because they are normally used for first evaluations, as a “rule of thumb”; moreover, some ancient methods are also presented to give an idea of the simplifications they required and the differences that can be expected using more advanced models. The governing equations are also derived, because an expert engineer should not use the computer programs as simple tools, always remember in a computer the old adage “trash in – trash out” applies; therefore, the assumptions underlying the equations should be known as well. This book is also devoted to the computer model development. To help the reader to build a model, a methodology is shown. To this end, not all the devices that can be found to reduce the effects of waterhammer have been described; instead, are

6 WATER HAMMER SIMULATIONS described, the most widely used devices and especially the assumptions and the methods to implement them in a computer program in addition to the problems found in their implementation and the reliability of the carried out results. Also, methods available to integrate the governing equations are nowadays in a large number: in the book only the method of characteristics and one based on finite differences (Lax–Wendroff) are shown. However, the methodologies shown in this book can be used as a guide in case the reader is interested in the implementation of a different scheme. In most cases, when a method or a model is described, it is also implemented in a computer program, enclosed with the book: this means seven computer programs are enclosed. The programs developed are not meant to be used for professional purposes: they have been tested and all the attention has been placed in their development, and they can be used for educational purposes or even for small real plants. However, these programs are very simple and do not allow to model complex plants, for which the reader should turn to commercial codes. Anyway, the presence of a running computer program shows that all the methods described have been tested and they work. Both the executable file and the source codes are enclosed, therefore the programs can be used as they are, in order to reproduce the examples presented in the book or to try new cases; but the main purpose is for the reader to study the list code to see how the different methods have been implemented, learn the technique and, if possible, improve the outcome.

1.4 The computer programs The computer programs enclosed with the book have been developed in Delphi [20], which is an object-oriented Pascal. Most of the code is devoted to input/output of data and for managing the hydraulic engine which is constituted by: • setting initial conditions; • starting the loop and call alternatively the procedures to solve: – internal points; – boundary conditions. The procedures to compute internal points and boundary conditions will be discussed in each case through the book, while the initial conditions are normally those of steady flow conditions. The steady flow conditions are computed using the Darcy–Weisbach formula, i.e.: J =λ·

V · |V | 2·g·D

(1.1)

where J is the distributed headlosses per unit length (or slope of the hydraulic grade line), λ is a dimensionless coefficient called the Darcy friction factor, V is the flow velocity, D is the pipe diameter and g is the acceleration of gravity. As


7

normally the roughness of the conduit is expressed with parameters as Manning’s or Hazen-Willimans’, in the programs developed for the book it has been decided to let the User introduce the former, which is easily found in all engineering practice manuals (e.g., Stephenson [21], who also compare the different formulae available in the literature). The Manning’s formula is: 1 2/3 1/2 ·R ·J (1.2) n where n is the Manning’s coefficient dependent on the pipe characteristics, and R is the hydraulic radius. As the pipe and flow characteristics are expressed in diameter D and discharge Q, in the computer programs the headlosses J are computed rearranging eqn (1.2) and obtaining: n2 · Q 2 (1.3) J = 10.29 · 16/3 D Combining eqns (1.1) and (1.3), the values of λ can be computed from the Manning’s coefficient with the formula: V =

n2 (1.4) D1/3 All the units are expressed in SI. The use of practical formulae like Manning’s implicitly assume the flow is turbulent (while the formula of Darcy–Weisbach can be used in all flow regimes), but this hypothesis is strictly verified in most practical cases, and however acceptable. Obviously the application of the formula in the cases presented in the book is trivial, as in the examples there will be one pipe only and because a reservoir is always present (upstream or downstream), the head in each required point in the pipe is computed starting from the reservoir and adding or subtracting the headlosses. For complex networks, a nonlinear system has to be written and many methods have been developed for its solution, which are not part of the topic of this book (for instance see Ref. [22] for a solid handbook or Ref. [23] for a scientific approach and a selected list of references). Results are printed in two text files: one reports the heads and the other the velocities computed during the simulations; these files can be imported for analysis in any spreadsheet. Each of the file has data printed in seven columns: the first column reports the time of simulation, the other six columns the required data (heads or velocities) located at the upstream of the pipe (second column) and at 20%, 40%, 60%, 80% and 100% of the length of the pipe – the last column is obviously related to the downstream condition. As the programs have been developed only for demonstration purposes, they have no checks or alerts in case of overflows or errors occurred during simulations: that mean they can abruptly crash without notice. On the other hand, they are sufficiently light and fast to be used practically in all computers and require few seconds to accomplish the average simulations. λ = 124.528 ·


Chapter 2 Compressible flow theory: basic concepts

Transients in pipes could be analyzed with two different approaches, depending on the hypothesis related to the flow compressibility. In the elastic theory approach (waterhammer) the liquid compressibility is computed, while in the anelastic approach (surge) the liquid compressibility is neglected. Although in most of the hydraulic phenomena liquid compressibility may be neglected, in the analysis of transients this is often a key parameter which cannot be ignored, as its mistreatment would bring completely mistaken results, except in very peculiar cases. On the other hand, a model which take into account the liquid compressibility is more general and obviously allows the analysis of the phenomena that could have it ignored. Therefore, in this book the liquid is considered compressible, even if in few cases comparisons with the anaelastic approach will be performed.

2.1 Instantaneous operations Let us consider the simple sketch shown in Figure 2.1, where for the time being the kinetic energy and the headlosses are neglected. The initial velocity in the pipe of area A is equal to V0 and L is the length of the pipe. At the upstream boundary there is a reservoir, which can be considered having a constant head, while at downstream boundary there is a valve, open in steady flow conditions and instantaneously closed at the beginning of the operations. When the velocity is stopped at the downstream end, under the hypothesis of uncompressible liquid, the whole column stops: the momentum of the water column (a finite quantity) goes to zero in an infinitesimal time, and therefore the force exerted by the liquid on the pipe is infinitive, which is neither true nor acceptable from a designer point of view. Thus, the hypothesis of uncompressible flow has to be discarded, accepting its elasticity: in this case, when the valve closes, the liquid is still flowing in the upstream part of the pipe, and, therefore, the quantity of water which actually stops in the infinitesimal time dt is an infinitesimal volume of length ds, while the remainder of the water column moves with the initial velocity. As the stopped volume is equal to A· ds, the stopped mass is equal to ρ · A · ds, and therefore, the variation of the momentum is equal to ρ · A · ds · V0 . This


Upstream reservoir

Downstream valve

Figure 2.1: Sketch of a simple plant. variation must be equilibrated from the impulse of the forces operating on the same mass; these forces are given by the increase of the pressure p which is found at the downstream boundary, and therefore the force is A · p · dt. As a consequence, the equilibrium can be expressed by: ρ · A · ds · V0 = A · p · dt

(2.1)

and therefore the increase of pressure due to the abrupt closure of a downstream valve is: p = ρ ·

ds · V0 dt

(2.2)

Let c = ds/dt be the wave celerity1 , eqn (2.2) can be written as: p = ρ · c · V0

(2.3)

See Figure 2.2 for a better understanding of the event and of the adopted symbols. The increment of pressure given by eqn (2.3) can be expressed in terms of water column, as: c h = · V0 (2.4) g Equation (2.3) or (2.4) is known as Allievi–Joukowsky’s formula, even if the literature shows Menabrea [1] to have made calculations with these equations. In the following instants, the subsequent volumes get stopped. The pressure wave propagates with celerity c towards the upstream reservoir, which is reached

1

In this book, the word velocity refers to the mass transportation, while the word celerity refers to the velocity of propagation of a wave.

COMPRESSIBLE FLOW THEORY: BASIC CONCEPTS

11

ds

V0

A

V

P0

0

A

c ds dt

P0 p

Figure 2.2: Forces equilibrium – sketch and symbols. at time: t1 =

L c

(2.5)

When t < t1 , the pipe is divided in two halves: the upstream part where the liquid is still moving and the downstream part where the liquid is stopped and pressures have increased of the quantity given by eqn (2.3) (Figure 2.3). As mentioned, when the time is equal to t1 the pressure wave reaches the upstream tank where the head is constant. As a consequence, at time t1 there is a difference in the pressures of the liquid at the upstream end of the pipeline, while the velocity of the liquid is null everywhere. Due to this pressure differences, water starts to flow from the pipeline towards the tank, with velocity −V0 , while pressures are back to those of steady flow. This second phase lasts until time t2 = 2L/c when the pressure wave reaches again the downstream valve. In the time t1 < t < t2 , the pipeline is still divided in two halves: in the upstream part, the liquid is still moving with velocity −V0 and pressure is equal to that of steady flow; in the downstream part, the liquid is stopped and there is high pressure. At time t2 , which is the end of the second phase, the pressure wave reaches again the downstream valve where obviously the velocity must be zero. With the same reasoning carried for the initial phase, but considering that now in the pipe the velocity is equal to −V0 , the wave reflection induces pressures equal to p = −ρ · c · V0 . The third and fourth phases are very similar to the first and second, respectively, but in this case the pressures have opposite sign. At time t4 = 4L/c, the initial conditions are reached again and, therefore, the cycle starts again. These four phases reproduce continuously and, theoretically, never end; in the real world, the waves smooth because of the headlosses and therefore the phenomenon ends in a finite time, which can be very short.

12 WATER HAMMER SIMULATIONS ∆h c V0 g

∆h c V0 g

First phase

Second phase

∆h c V0 g

∆h c V0 g

Third phase

Fourth phase

Figure 2.3: Four phases of transient.

2.2 Wave celerity As it has been observed, the characteristics of the event depend substantially on the wave celerity. In the following, the expression of the celerity is reported only in the case of non-deformable pipes, while in the Section 2.5.1 this value will be carried out for more general cases. The bulk modulus of elasticity of the fluid ε is defined as: ε=ρ·

dp dρ

(2.6)

Let us consider again Figure 2.2 which reproduce the conditions in a generic time t in a pipe after the instantaneous closure of the downstream valve. In the subsequent time dt the section that separates the stopped fluid (downstream) from that in motion (upstream) moves upstream from A−A , together with the pressure wave; instead, the section A−A itself moves downstream of a distance V0 · dt because the volume W between that section and the downstream valve is compressed; the volume decreases of a value equal to: dW = −A · V0 · dt

(2.7)


13

At the same time on the same volume W acts an increased pressure given by eqn (2.3) to which, because of eqn (2.6), corresponds a decrement in volume: dW = −

W · p A · ρ · c · V0 =− · ds ε ε

(2.8)

Equalling eqns (2.7) and (2.8) and reminding that c = ds/dt, it follows that: c=

ε ρ

(2.9)

This is the value of sound celerity in the fluid. For water with a temperature equal to 8◦ C, this value is equal to c ∼ = 1425 m/s, which is an increase of about 3 m/s for each increased degree of temperature. If the liquid is more compressible, as it happens for certain oils, the celerity c decreases. The wave celerity also decreases when the elasticity of the pipeline wall is considered. For a pipe with diameter D, thickness e and modulus of elasticity E, eqn (2.9) becomes (see Section 2.5.1): c=

ε ρ

1+λ·

ε·D E·e

(2.10)

λ being a coefficient which considers the effects of the junctions between pipes and can be assumed equal to 1 for a single pipe. For steel pipes, a widespread assumption is c = 1000 m/s, while for plastic pipes this value is lower because their elasticity is higher. As will be discussed in Chapter 8, the presence of air may have a significant effect on this parameter. Example 1 Let us compute the wave celerity in a steel pipe with diameter D = 500 mm and thickness e = 6.3 mm, knowing that the steel elasticity modulus is equal to E = 2.0 × 1011 N/m2 , the bulk elasticity modulus for water is equal to ε = 2.14 × 109 N/m2 and water density is ρ = 1000 kg/m3 . Wave celerity is equal to: c=

ε ρ

D·ε 1+ e·E

=

2.14 × 109 1000

1462.87 ∼ = √ = 1075 m/s 1.85 500 × 2.14 × 10 1+ 6.3 × 2.0 × 1011 9

14 WATER HAMMER SIMULATIONS About the dimensions, it is to be noted that: kg ε ε m m2 m · s2 = 2 ⇒ kg s ρ s ρ m3 while:

1+

D·ε e·E

is non-dimensional.

2.3 Velocity of operations Let us consider now, instead of an instantaneous closure of the valve, the case when a valve is closed in a finite time Tc > 0. For the sake of simplicity let us evaluate the velocity outflowing from the valve with the following linear expression, which, however, can be easily generalized: ⎧ ⎪ ⎨V (t) = V0 · 1 − t when t ≤ Tc Tc (2.11) ⎪ ⎩ V (t) = 0 when t > Tc In the first instants of the transient, downstream pressures can be computed by the already-known formula p = ρ · c · V . As the velocity linearly decreases because of eqn (2.11), pressure linearly increases. Let Tc the time when the valve is completely closed, for a given time t > Tc , in the section where the valve is positioned the pressure due to the transient can still be computed by eqn (2.3); the front of the wave which moves upstream is inclined, as shown in Figure 2.4. As already mentioned, when the wave reaches the upstream reservoir, it is reflected downstream, where it propagates as descendant perturbation. If the pipe is relatively short, the valve may close in such a long time that the closure is still in progress when the reflected wave is back; this time τ0 is called “critical phase duration” and it is equal to: 2·L τ0 = (2.12) c where L is the length of the pipe. When Tc > τ0 , the pressure in the section of the valve cannot completely develop and the pressure in the section where the valve is positioned is equal to: p = ρ · c · (V0 − Vf ) (2.13) where Vf is the velocity at the valve when t = τ0 .


15

Tc

∆h

c V g 0

Figure 2.4: Pressure wave moving upstream when the downstream valve closes linearly. The maximum value of pressure that is reached in this case can be computed considering Figure 2.5 and applying simple geometric considerations, obtaining: c · V0 p g = 2 ·L Tc c

(2.14)

Rearranging and simplifying the equation, the so-called Michaud’s formula [1] is carried out: 2 · ρ · L · V0 p = (2.15) Tc The increment of pressure due to the transient can again be expressed as water column and, therefore: 2 · L · V0 h = (2.16) g · Tc Example 2 Let us consider a pipe with length L = 1000 m where the wave celerity is equal to c = 1000 m/s and the flow has a velocity equal to V0 = 5 m/s. Pressures have to be computed in both cases of valve closure with time Tc = 0.8 s and 5.0 s. The phase duration of the system is τ = 2L/c = 2 s. In the former case, we have τ0 < Tc and, therefore, pressures have to be computed with the Allievi–Joukowski formula, obtaining: p = ρ · c · V0 = 1000 kg/m3 · 1000 m/s · 5 m/s = 5 × 106 kg/m2 m/s2 = 5 × 106 Pa

16 WATER HAMMER SIMULATIONS p

ρcV0

2L c

Tc

t

Figure 2.5: Computation of waterhammer pressure for long operations. or, in terms of water column: h =

1000 m/s c · 5 m/s ∼ · V0 = = 510 m g 9.806 m/s2

In the latter case, instead, we have τ0 > Tc and, therefore, pressures have to be computed with the Michaud formula: p =

2 × 1000 kg/m3 · 1000 m · 5 m/s 2 · ρ · L · V0 = = 2 × 106 Pa Tc 5s

or, in terms of water column: h =

2 · 1000 m · 5 m/s ∼ 2 · L · V0 = = 204 m g · Tc 9.806 m/s2 · 5 s

2.4 Non-negligible headlosses Now let us consider a pipe where headlosses cannot be neglected, as in the case shown in Figure 2.6. In the time interval t1 , the mass 1 stops and the water column increases to a quantity equal to c · V0 /g. The increased head propagates to the mass 2, which, however, has an initial (steady flow) pressure larger than the mass 1, the difference being equal to J · x. For this reason, the mass 2 is pressed by a downstream head equal to cV0 /g − J · x which is obviously smaller than the value c · V0 /g and is necessary to stop mass 2 completely.


17

J ⋅∆x

c ⋅ V0 g Steady flow Hydraulic grade line

3

2

1

Valve completely and instantaneously closing

Figure 2.6: Pressures due to transients when headlosses cannot be neglected. In other words, the pressure generated at the valve section is not sufficient to stop the mass 2, which, therefore, has a residual velocity equal to V = J · x · g/c. Because of this residual velocity, mass 2 compresses again the mass 1, further increasing its pressure, until the mass 2 stops completely; this happens when for that mass the increased head is equal to c · V0 /g. The increased head on the mass 1 is, therefore, equal to c · V0 /g + J · x. Repeating this procedure for the subsequent temporal steps (and analogously for the following masses), it can be shown that the increased head at the downstream valve is equal to c · V0 /g, to be added to the static head.

2.5 Governing equations The governing equations are the classic continuity and momentum; they are usually written in one dimension, along the pipe, and in this form they will be derived. In the following, the flow will be characterized by its average velocity V (x, t) = Q(x, t)/A(x, t) and pressure p(x, t) computed in the centre of mass of the pipe cross section. 2.5.1

Continuity equation

The continuity equation has to be written considering that during transients the pipeline wall deforms because of the pressure variations. This deformation is a


V

∂V ∂x

∂u ∂x

V u x

Figure 2.7: Control volume for determining the continuity equation. quite complex phenomenon, as it consists not only of radial shape variations, but also elongations or shortening of the pipeline in the longitudinal direction. The continuity equation will be carried out referring to a control volume considered fixed with respect to the pipe, which therefore follows all deformations. The sketch of the phenomenon and the adopted symbols are reported in Figure 2.7. Let u be the velocity of the pipe wall in the direction of the pipe axes; the mass that enters the control volume is equal to ρ · A(V − u), while the mass outflowing downstream is (except for infinitesimal of higher order): ρ · A · (V − u) +

∂[ρ · A · (V − u)] · δx ∂x

(2.17)

The difference between the masses inflowing and outflowing the control volume must be equal to the increase of the mass inside the volume itself; therefore, still neglecting infinitesimal of higher order: −

∂[ρ · A · (V − u)] D · δx = (ρ · A · δx) ∂x Dt

(2.18)

where the symbol D /Dt denotes the total derivative2 with respect to the axial motion of the pipe.

2

The total derivative of a function f with respect to the pipe motion, which has velocity u, is ∂f + u · ∂f . ∂t ∂x


Developing and keeping into account that

19

D (δx) ∂u = · δx, it yields: Dt ∂x

∂(ρ · A · V ) ∂(ρ · A) + =0 ∂x ∂t

(2.19)

Rearranging again and denoting with a dot positioned above the variable the substantial derivative3 with respect to the average flow velocity, it yields: A˙ ρ˙ ∂V + + =0 A ρ ∂x

(2.20)

So far the most general case has been handled without making any hypothesis neither on the shape or the material of the pipe nor on the transported fluid. Now few hypothesis must be introduced, which are quite common in the practical cases of hydraulic engineering. Hypothesis 1 (on the fluid). The fluid is Newtonian, i.e., denoted with ε the bulk elasticity modulus, its rehology can be described by: ρ˙ p˙ = ρ ε

(2.21)

Hypothesis 2 (on the pipe geometry). The pipe is cylindrical with circular base. Let ξD be the specific extension of the diameter, it yields: ˙ A˙ D = 2 · = 2 · ξ˙D A D

(2.22)

Hypothesis 3 (on the material of the pipe). The material that constitutes the pipe is mechanically homogeneous and isotropic, with linear elastic behaviour and the wall thickness e is negligible when compared with the diameter. Therefore, the Hook’s law applies for two dimensions stress states, i.e.: ξC = ξD =

ξx =

σC − µ · σx E

σx − µ · σC E

(2.23)

(2.24)

where ξC is the specific extension of the wall in the radial direction, ξx is the specific extension of the wall in the axial direction, σC is the normal stress in the

3

The substantial derivative of a function f with respect to the average flow velocity V is + V · ∂f . ∂x

∂f ∂t

20 WATER HAMMER SIMULATIONS radial direction, σx is the normal stress in the axial direction, µ is the Poisson’s coefficient (0 ≤ µ ≤ 0.5) and E is the bulk elasticity modulus of the material of the pipe. If the pressure is large enough, its variations on the generic section of the pipe can be neglected, i.e., the pressure p can be assumed as constant along the wall of the pipe and therefore the Mariotte’s4 formula applies: σc =

p·D 2·e

(2.25)

The axial stress depends on the constraints to the pipe. If the pipe is free to expand or reduce in the axial direction, then the axial stress is always zero (this is the case of a pipe with many expansion joints) so that σx = 0. On the contrary, if the pipe is constrained so that ant axial deformation is prevented, it is obviously ξx = 0. Therefore, the expression (2.22) for the former case (axial stress equal to zero) becomes: •

A˙ p·D = A E·e

(2.26)

while in the latter case (axial deformation equal to zero) it becomes: A˙ D = A Dt

p·D · 1 − µ2 E·e

(2.27)

Both eqns (2.26) and (2.27) can be written in the general case as: •

A˙ p·D =λ· A E·e

(2.28)

where λ depends on the constraint conditions and varies within the range 0.75 and 1. This parameter depends only on µ that is constant in space and time because of the hypothesis on the pipe material. As mentioned, the wall thickness is negligible if compared with the diameter, and because the stress can be considered twodimensional, the thickness e can also be considered constant, both in space and time. As a consequence, eqn (2.28) can be simplified as: A˙ λ ˙ = · (˙p · D + p · D) A E·e

4

(2.29)

Mariotte [24] was able to show the strength of the weakest material element in a structure is likely to decrease with increasing structure size.


˙= Considering D

21

D A˙ ˙ finally yields: · , substituting and extracting A/A 2 A A˙ = A

p˙ · D E·e p·D 1−λ· 2·E·e λ·

(2.30)

In the usual applications, the denominator on the right hand equation is very close to one, which allows a drastic simplification of eqn (2.30)5 . This simplified expression, substituted in the continuity equation (2.20) and also considering the eqn (2.221) allows to obtain the following: λ·

p˙ · D p˙ ∂V + + =0 E·e ε ∂x

(2.31)

that can also be written as: ε·D p˙ 1 + λ · E · e ∂V · =0 + ε ρ ∂x ρ

(2.32)

Keeping in account the celerity term defined in eqn (2.10) in Section 2.2, eqn (2.32) becomes: ∂V =0 (2.33) ∂x and finally, developing the substantial derivative, the expression of the continuity equation is: p˙ + ρ · c2 ·

∂p ∂V ∂p +V · + ρ · c2 · =0 ∂t ∂x ∂x 2.5.2

(2.34)

Momentum equation

The momentum equation can be written balancing all the forces applied to a part of the flow between two pipe sections that have distance δx. In Figure 2.8, those sections are shown with all the forces applied. The different forces are as follows. Forces due to the pressures on the upstream and downstream areas are, respectively: m = p · A

(2.35)

In a normal steel pipe, the ratio between the diameter and the wall thickness is D/e = 50 ÷ 100 and the elasticity modulus of steel is ε = 2 × 1011 N/m2 . In this case with a pressure equal to a water column of 1000 m the mentioned denominator is equal to 0.997. 5


Πl Πv Mv I

Πm Mm x

x

Figure 2.8: Control volume to carry out the momentum equation. v = p · A +

∂(p · A) ∂p ∂A · δx + · · δx2 ∂x ∂x ∂x

(2.36)

The difference between eqns (2.36) and (2.35), neglecting higher infinitesimals, is: ∂(p · A) · δx ∂x

(2.37)

The axial component of the force applied by the pipe wall is: ∂p δx ∂A l = p + · · · δx ∂x 2 ∂x

(2.38)

The axial component of the weight of the fluid volume is: G · sin α = γ · A · δx · sin α

(2.39)

where sin α = dz/dx. For the Newton’s second law the addition of all the forces acting on a volume is equal to the inertia force due to the acceleration of the mass in the volume. Accelerating, the mass moves in the axis direction, and, therefore, the acceleration has to be computed as substantial derivative of the velocity V . Therefore: I = ρ · A · δx ·

∂V ∂V +V · ∂t ∂x

(2.40)


23

Note that, because of the assumption of the control volume, in this term the differences between the incoming and outcoming momentum flux is already computed6 . The resulting shear force, τ being the shear stress at the pipe wall, is equal to τ · π · D · δx. Under the hypothesis that the headlosses in unsteady flow conditions are similar to those in steady flow, it yields: τ=

γ ·D·J 4

(2.41)

J being the slope of the hydraulic grade line. Adding all the terms: ∂(p · A) ∂A dz ∂V · δx − p · · δx + γ · A · δx · + ρ · A · δx · ∂x ∂x dx ∂t ∂V γ ·D·J + · π · D · δx = 0 ∂x 4

(2.42)

∂p ∂V ∂V dz +ρ· +ρ·V · +γ ·J +γ · =0 ∂x ∂t ∂x dx

(2.43)

+ρ · A · δx · V · which can be simplified as:

2.5.3 The governing equations Finally, the governing equations can be written together as: ∂p ∂p ∂V +V · + ρ · c2 · =0 ∂t ∂x ∂x

(2.34)

∂p ∂V ∂V dz +ρ· +ρ·V · +γ ·J +γ · =0 ∂x ∂t ∂x dx

(2.43)

where eqn (2.34) is the continuity equation and eqn (2.43) is the momentum equation. Being two partial differential equations, these must be integrated to both initial and boundary conditions.

the local inertias have to be computed, which have to be written as I = Differently, ∂(ρ·V ) ∂(ρ·V ) · dV = · A · dx and the momentum flux are m = ρ · V 2 · A and v = ρ · V 2 · ∂t ∂t

6

V

·A) · dx + ρ · ∂A · V · dx. In the latter equation, the second term is the momentum A + ∂(ρ·V ∂x ∂t flux through the plane section, while the third term is the momentum flux through the lateral surface. Adding up and simplifying these terms, keeping in account that ρ · ∂A · V =V · ∂t ∂(ρ·A) ∂(ρ·A·V ) ∂(ρ·A) ∂ρ − V · A · and that (because of the continuity equation) + = 0, the same ∂t ∂t ∂x ∂t result is achieved. 2

24 WATER HAMMER SIMULATIONS These equations have been carried out using the unknown parameters the pressure p and the velocity V . It is quite obvious that similar results could have been carried out using unknowns the head h and the discharge Q (or any of their combination) being, respectively: p γ

(2.44)

Q =V ·A

(2.45)

h=z+

Independent variables are time t and the longitudinal abscissa x. As mentioned earlier, the other variables c, ρ (or γ ), D and the pipe roughness f do not vary significantly during the transient and therefore they can be let constants. Rigorously, the hydraulic grade line does not depend only on the roughness f (which in the following will be let constant), but also from the Reynolds number Re; however, in practical cases these variations are negligible. In matrix form, the governing equations are written as: V p ∂ +1 ∂t V ρ

ρ · c2 0 p ∂ · = dz V ∂x V −g · −g·J dx

(2.46)

or: ∂U · ∂ U +B =E ∂t ∂x where: p = U V

V = 1 B ρ

ρ · c2 V

0 = E dz −g · −g·J dx

(2.47)

(2.48)

can be found solving the equation (in the The eigenvalues of the matrix B unknown λ): (V − λ)2 = c2

⇒

λ=V ±c

(2.49)

As both eigenvalues are real and different, the system of differential equation is said to be hyperbolic. The governing equations can be solved in different ways, as will be seen in Chapters 4 and 5, but some preliminary consideration has to be drawn beforehand and these will be presented in Chapter 3 together with few simplified methods to carry out approximate solutions.

Chapter 3 Simplified solutions

With the availability of computer the governing equations carried out in Chapter 2 are quite easily numerically integrated. In the past, if only the peak of pressure was needed the Allievi–Joukowsky formula (2.3) or (2.4) or the Michaud formula (2.15) or (2.16) were used. If the values of pressure and velocity were required for the whole transient, simplified approaches were used. In this chapter, three methods are presented. The first is the Allievi’s method, which uses simplified equations and finds the solution iteratively; the second is based on the hypothesis that the fluid is non-elastic, again to find a simplified solution and the third is the graphical method.

3.1 The Allievi’s method (1913) This method has been presented in 1913 [5] and the equations obtained are known as the Allievi chained equations. The governing equations are simplified as follows. Continuity equation is written as: ∂p ∂V + ρ · c2 · =0 ∂t ∂x

(3.1)

∂p ∂V +ρ· =0 ∂x ∂t

(3.2)

while momentum equation is:

These equations are valid under the following hypothesis: 1. velocity V is negligible when compared with the wave celerity c and, therefore, also terms like |V · ∂F/∂x| are negligible when compared with terms like |∂F/∂t|; 2. headlosses are negligible; 3. density variations, due for instance to elevation variations, are negligible. The usual sketch of simple plant is considered (Figure 3.1). Having the abscissa starting from the downstream valve (x = 0) and assuming the x-axis positive in the


Upstream reservoir

Downstream valve

x x

0

Figure 3.1: Sketch of a simplified plant. upstream direction, the governing equations can be written as: ∂V ∂p = ρ · c2 · ∂t ∂x

(3.3)

∂p ∂V =ρ· ∂x ∂t

(3.4)

where eqn (3.3) is the continuity equation and eqn (3.4) is the momentum equation. The initial conditions, which define the pressures and the velocities before the beginning of the transient, are: p(x, t = 0) = p0

(3.5)

V (x, t = 0) = V0

(3.6)

Boundary conditions are given at the upstream reservoir: p(x = L, t) = p0

(3.7)

and at the downstream valve, where the law of velocity variation is known: V = V (x = 0, t)

(3.8)

Deriving eqn (3.3) with respect to t and eqn (3.4) with respect to x, and eliminating the term in V , a single second-order differential equation can be carried out: 2 ∂ 2p 2 ∂ p = c · ∂t 2 ∂x2

(3.9)

Similarly, deriving eqn (3.3) with respect to x and eqn (3.4) with respect to t and eliminating the term in p, it yields: 2 ∂ 2V 2 ∂ V = c · ∂t 2 ∂x2

(3.10)

SIMPLIFIED SOLUTIONS

27

Equations (3.9) and (3.10) are the D’Alembert equations [25] and their solution is well-known in Mathematical Analysis [26]. The solution of eqn (3.9) can be written as: p = p − p0 = F (x − c · t) + f (x + c · t)

(3.11)

where F and f are functions to be defined on the basis of boundary conditions. This equation describes the trend of the pressures during transients. With regards to the velocity, its derivative can be computed deriving eqn (3.11) with respect to x and substituting the result in the continuity equation (3.3), obtaining1 : −c · ρ ·

∂F(x − c · t) ∂f (x + c · t) ∂V = − ∂t ∂t ∂t

(3.12)

This equation, integrated with respect of time, yields: c · ρ · V = c · ρ · (V − V0 ) = F(x − c · t) − f (x + c · t)

(3.13)

With eqns (3.11) and (3.13), it is possible to compute velocities and pressures in the pipe. For each couple of values (xi , ti ) for which the following relation is valid: x1 − c · t1 = x2 − c · t2

(3.14)

the function F must have the same value. The same is for the function f , for which the relation that ties the couple of coordinates is: x1 + c · t1 = x2 + c · t2

(3.15)

That means the values the function F assumes at time t1 and at time t2 given by: t2 = t1 +

x2 − x1 c

(3.16)

must be equal or, in other words, that the function F describes a wave which propagates with celerity c. As mentioned, for the function f the same reasoning can be done. Function F is, therefore, related to a perturbation that propagates from the valve to the reservoir, while function f is a wave which has opposite direction (i.e., from the reservoir towards the valve). Few simple considerations allow to further simplify the problem, letting it to be solved with only one equation.

1

It can be demonstrated that:

∂p ∂x

=

∂F ∂x

+

∂f ∂x

= − 1c ·

∂F ∂t

−

∂f ∂t

.

28 WATER HAMMER SIMULATIONS At the section located near the reservoir, i.e., x = L, the pressure p must be constant, and, therefore: p = F (L − c · t) + f (L + c · t) = 0

(3.17)

At the junction between the reservoir and the pipe, the two functions have the same value but opposite sign. In other words, when from the downstream pipe a wave arrives to the reservoir, a new wave is generated, with opposite sign, which propagates downstream. This result has been already found in Section 2.1. Therefore, if at time t1 a disturbance f · (L − c · t1 ) is generated at the reservoir, it propagates downstream with opposite sign; at a subsequent time t > t1 , the disturbance is located in the section x, according to the expression t = t1 + L −c x and due to the relation L + c · t1 = x + c · t. This relation together with eqn (3.17) of the reservoir, it yields: f (x + c · t) = f (L + c · t1 ) = −F(L − c · t1 ) A more general expression can be written as follows: 2 · (L − x) f (x + c · t) = −F x − c · t − c

(3.18)

(3.19)

which connects the two functions. Expressions (3.11) and (3.13) then become much simpler: 2 · (L − x) p = p − p0 = F(x − c · t) − F x − c · t − (3.11 ) c 2 · (L − x) c · ρ · V = c · ρ · (V − V0 ) = F(x − c · t) + F x − c · t − c (3.13 ) As in the first phase of the transient only the ascending wave is present, the equations simplify to: p = p − p0 = F(x − c · t)

(3.11 )

c · ρ · V = c · ρ · (V − V0 ) = F(x − c · t)

(3.13 )

Removing F from both, the well-known formula p = ρ · c · V is derived again. Let τ = 2·L , for the section x = 0 the eqns (3.11 ) and (3.13 ) can be finally c written as: p = p − p0 = F(t) − F(t − τ )

(3.11 )

c · ρ · V = c · ρ · (V − V0 ) = F(t) + F(t − τ )

(3.13 )


29

In order to study the second phase of the transient, when the ascending waves adds with those descending, eqns (3.11 ) and (3.13 ) are used, choosing a sequence of times ti that: 0 < t1 < τ t2 = t1 + τ t3 = t2 + τ .. . ti = ti−1 + τ so that, for t = t1 : p(t1 ) − p0 = F(t1 )

c · ρ · [V0 − V (t1 )] = F(t1 )

(3.20)

being still in the first phase. For the following instants: t = t2 : p(t2 ) − p0 = F(t2 ) − F(t1 )

c · ρ · [V0 − V (t2 )] = F(t2 ) + F(t1 )

(3.21)

.. . t = ti : p(ti ) − p0 = F(ti ) − F(ti−1 )

c · ρ · [V0 − V (ti )] = F(ti ) + F(ti−1 )

(3.22)

Here again, removing F from the two series of equations, it yields: p(t1 ) − p0 = c · ρ · [V0 − V (t1 )] p(t2 ) + p(t1 ) − 2 · p0 = c · ρ · [V (t1 ) − V (t2 )] .. .

(3.23)

p(ti ) + p(ti−1 ) − 2 · p0 = c · ρ · [V (ti−1 ) − V (ti )] With eqns (3.23), the study of the pressures at the valve during the transient is possible, at any time. In these equations, in fact, the time step has to be equal to τ : but, as the initial time t1 is chosen arbitrarily, under the only condition t1 < τ , the values can be actually computed for any time. For a generic section with abscissa x the procedure is similar; the computations start from the values p and V at the downstream valve. Neglecting all the passage for the sake of brevity, the final equations are: x x − F0 tn−1 + x p = γ · F0 tn − c c

(3.24)

30 WATER HAMMER SIMULATIONS x V =

g x x · F0 tn − + F0 tn−1 + c c c

(3.25)

In eqns (3.24) and (3.25), x p and x V refers to pressure and velocity variations, respectively, at the section x. With F0 (t) the value that the function F has at the section x = 0 at time t is indicated. Values x p and x V can be iteratively computed with the following simple procedure: 1. pressure and velocity values of at the downstream valve are computed with eqn (3.23) at different times; 2. values of the function F0 (t) are computed iteratively from: F0 (ti ) = F0 (ti−1 ) + [p(ti ) − p0 ]; 3. values of x p and x V are computed with eqns (3.24) and (3.25).

3.2 The non-elastic hypothesis 3.2.1 Governing equations In Figure 3.2, the classic scheme of an hydropower plant is shown. This is the simplest scheme of an hydropower plant, with an upstream reservoir; a tunnel, which is usually quite large and where velocities are relatively low, in order to reduce headlosses; a surge tank, the use of which will be discussed in Chapter 6; a penstock, where velocities are very high, as at the end there is a turbine and the target is the production of the maximum possible amount of electricity. This plant can be studied as a whole with the mathematical methods described later in this book. However, in the past this was not possible and the only way to study the behaviour of the plant was splitting it in two parts: first the penstock, which was studied with theAllievi or Michaud’s formula or theAllievi’s method; then the tunnel was studied, pretending the downstream pipeline was nonexistent and therefore positioning the downstream valve immediately downstream the surge tank; moreover, the valve was supposed to close instantaneously.

A1 A R A

B Downstream valve

Tunnel Reservoir

Penstock Turbine T

Figure 3.2: Sketch of a classic hydropower plant.


31

These hypotheses mean the two phenomena (in the tunnel and in the penstock) are not connected or, better, that they can be studied separately; this is not so senseless, as the phenomenon in the penstock depends on wave celerity (around 1000 m/s) and it is so fast that it can be considered finished when the transient in the tunnel starts; the latter depends on the flow velocity (few meters per second), therefore being much slower. However, in this section the accuracy of these hypotheses will be tested. In this simplified system, where only the tunnel is considered, the hypotheses of non-elasticity of the fluid and of the non-deformability of the pipe can be accepted. In essence, the water flowing in the tunnel cannot enter the penstock and therefore flows in the surge tank, rising its water level and consequently the piezometric head at the bottom of the tank itself. The discharge reduces until the maximum level in the tank is reached, which is higher than the level in the upstream reservoir because of the water inertia; when that maximum point is reached, the velocity of the water in all points of the system is zero. At this point, the system has not reached the equilibrium, and therefore the level in the tank decreases while the flow in the pipe has negative velocity, which becomes zero when the level in the tank is minimum. The transient continues with a series of oscillations that reduce their amplitude because of the headlosses. This phenomenon can be someway considered similar to the oscillations that can be recorded in an U-shaped pipe, shown in Figure 3.3. As the fluid is considered uncompressible and without headlosses, Bernoulli’s theorem is: ∂H β ∂V + · =0 ∂s g ∂t

(3.26)

where the Coriolis’ coefficient β ≈ 1. Because of the non-compressibility of the fluid and the non-deformability of the pipe walls, the continuity equation becomes V (s) = constant, and therefore: ∂V (s) dV = ∂t dt

(3.27)

As a consequence, the Bernoulli’s theorem can be integrated, resulting: H2 − H1 = −

1 dV · · g dt

2 1

ds = −

L dV · g dt

(3.28)

Let z the elevation of the free surface, because of the continuity equation, the following equation can be written: A · dz = −Q · dt ⇒

dz d dV 1 d2 = −V ⇒ = −2 · V ⇒ =− · 2 dt dt dt 2 dt

(3.29)

in level ∆

Difference


Area A

Length L

Figure 3.3: Mass oscillation in a U-shaped pipe. being the difference in elevation between the two free surfaces, as shown in Figure 3.3. Finally, Bernoulli’s theorem becomes: d2 + ω2 · = 0 dt 2 where:

ω=

2·g L

(3.30)

(3.31)

Imposing the initial conditions: (t = 0) = 0 V (t = 0) = 0

(3.32)

a periodical motion is obtained, i.e., oscillations, described by: = 0 · cos(ω · t)

(3.33)

If the headlosses cannot be neglected, Bernoulli’s theorem has to be modified as: ∂H β ∂V + · +J =0 ∂s g ∂t

(3.34)

which, upon integration, yields: L dV · −+L·J =0 g dt

(3.35)


33

This equation, together with the continuity equation, form a set which in general does not allow analytical integration. It is not necessary to describe too much in detail these equations, but it has to be said that in turbulent flow and in some cases in laminar flow, the motion is oscillatory and it is smoothed because of the flow resistances. Particular cases may be faced when the flow is laminar and the viscosity of the fluid is very high: in these cases, the initial different elevation asymptotically tends to zero with no changes in sign. However, generally speaking and assuming the presence of concentrated headlosses in the connection between the tunnel and the surge tank (see Chapter 6), the governing equations are: ⎧ L dV ⎪ ⎨ · +Z ±H ±K =0 g dt ⎪ ⎩ Asurge tank · dZ = Q − Qds dt

(3.36)

where Z is the water elevation in the surge tank, computed from the static level; H is the value of the headlosses in the tunnel; K the concentrated headlosses in the connection between the tunnel and the surge tank; Q is the discharge in the tunnel at time t and Qds is the discharge at the same time in the penstock. 3.2.2 Comparisons It is now interesting to compare the results that can be obtained with this simplified model with those that can be carried out with a numerical model that uses the complete equations described in Chapter 2. To do that, let us consider a plant with the following characteristics: Q0 = 15 m3 /s discharge; diameter of the tunnel; Dtunnel = 2.80 m L = 4200 m length of the tunnel; Dsurge tank = 11.50 m diameter of the surge tank. The surge tank has circular cross section and no concentrated headlosses at the junction with the tunnel. With the mentioned hypotheses and integrating with finite differences, the governing equations (3.36) become: ⎧ L V ⎪ ⎨ · + Z ± β · V 2 (t) = 0 g t ⎪ ⎩ Asurge tank · Z = Q(t) = Atunnel · V (t) t

(3.37)

In eqn (3.37), the discharge in the penstock is Qds = 0, that means the transient is finished. Distributed headlosses in the tunnel are computed as: β=

0.018 × 4200 λ·L = = 1.377 2·g·D 2 × 9.806 × 4200

34 WATER HAMMER SIMULATIONS In this case, the resistance term is λ = 0.018, which means full turbulent flow and a Manning resistance term equal to n = 0.143 s · m−1/3 . As the initial velocity in the pipe is equal to 2.44 m/s, the steady flow headlosses are equal to H = β · V 2 = 8.17 m = −Z0 . The method of integration is explicit and really simple: given a timestep t, from the continuity equation, the difference in the water elevation z is computed by: Z =

Atunnel · Vi · t Asurge tank

(3.38)

and therefore the elevation at the new time ti+1 is given by: Zi+1 = Zi + Z

(3.39)

Then the velocity variation is computed: V = −

g · (Zi+1 ± β · Vi2 ) · t L

(3.40)

and finally the velocity at the new time: Vi+1 = Vi + V

(3.41)

Results of the computations, performed with a spreadsheet (but they can be carried out by hand), are reported in Figure 3.4 and in Table 3.1. 10

Velocity (m/s) Elevation (m)

5 0 5 10 15

Water elevation in the surge tank Velocity in the tunnel

20

0

100

200

300 Time (s)

400

500

600

Figure 3.4: Oscillations in the pipe and in the surge tank, computed with the nonelastic model.


35

Figure 3.5 reports the water elevations in the surge tank carried out with this model and those obtained with the more complex simulation models described in the following. As can be seen, the agreement between the two sets of data is very good, at least at the beginning, while in the subsequent instants of the simulation, the anelastic model tends to underestimate the headlosses, thus being on the safe side. Table 3.1: Numerical results of the non-elastic model: water levels and velocities in the surge tank. t (s) 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300

Z (m) 1.44 1.42 1.39 1.34 1.28 1.21 1.13 1.04 0.96 0.87 0.77 0.67 0.58 0.48 0.37 0.27 0.17 0.07 −0.03 −0.13 −0.23 −0.33 −0.43 −0.52 −0.61 −0.70 −0.79 −0.87 −0.95 −1.02

Z (m) −8.17 −6.73 −5.30 −3.92 −2.58 −1.30 −0.10 1.03 2.08 3.03 3.90 4.67 5.34 5.92 6.39 6.77 7.04 7.21 7.28 7.25 7.12 6.89 6.56 6.14 5.62 5.00 4.30 3.51 2.64 1.69 0.67

V (m/s) −0.034 −0.062 −0.085 −0.103 −0.119 −0.131 −0.141 −0.148 −0.155 −0.160 −0.163 −0.166 −0.168 −0.170 −0.171 −0.171 −0.171 −0.171 −0.169 −0.168 −0.166 −0.163 −0.160 −0.156 −0.151 −0.146 −0.139 −0.131 −0.122 −0.111

V (m/s) 2.44 2.40 2.34 2.26 2.15 2.03 1.90 1.76 1.61 1.46 1.30 1.14 0.97 0.80 0.63 0.46 0.29 0.12 −0.05 −0.22 −0.39 −0.55 −0.72 −0.88 −1.03 −1.18 −1.33 −1.47 −1.60 −1.72 −1.83

t (s)

Z (m)

Z (m)

V (m/s)

V (m/s)

310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600

−1.09 −1.15 −1.19 −1.23 −1.26 −1.28 −1.27 −1.25 −1.21 −1.15 −1.06 −0.95 −0.82 −0.66 −0.47 −0.27 −0.06 0.16 0.37 0.57 0.76 0.92 1.05 1.16 1.23 1.28 1.31 1.31 1.29 1.26

−0.42 −1.56 −2.76 −3.99 −5.25 −6.53 −7.80 −9.05 −10.26 −11.42 −12.48 −13.43 −14.25 −14.90 −15.38 −15.65 −15.71 −15.55 −15.18 −14.60 −13.84 −12.92 −11.87 −10.71 −9.48 −8.20 −6.90 −5.59 −4.29 −3.03

−0.098 −0.084 −0.066 −0.046 −0.023 0.004 0.034 0.068 0.105 0.145 0.188 0.231 0.272 0.309 0.339 0.359 0.366 0.361 0.342 0.311 0.271 0.224 0.176 0.128 0.082 0.041 0.005 −0.026 −0.053 −0.075

−1.93 −2.02 −2.08 −2.13 −2.15 −2.15 −2.11 −2.05 −1.94 −1.80 −1.61 −1.38 −1.11 −0.80 −0.46 −0.10 0.27 0.63 0.97 1.28 1.55 1.78 1.95 2.08 2.16 2.20 2.21 2.18 2.13 2.05

36 WATER HAMMER SIMULATIONS 10

5

Elevation (m)

0

5 10 15 20

Anelastic model Elastic model 0

100

200

300 Time (s)

400

500

600

Figure 3.5: Water elevation in the surge tank: comparison between elastic and nonelastic model results. Actually, the hypothesis of non-compressibility of water moving in the tunnel-surge tank is surely acceptable; however, the differences in the effort of implementation of the two models is considerable. Design of simple surge tanks is surely possible even with the simplified model. The disadvantage of these models is not given by their precision, which is proven to be sufficiently good, but instead on their poor versatility. Simplified models cannot keep in due account the devices that can be installed in the surge tanks, or the ramification the pipelines may have, or other problems. With complex models, all these devices can be implemented, even if no model is exempt from errors, being always necessary to critically evaluate the coherence of the carried out results and to contrast the model instabilities.

3.3 Graphical method This is probably the most visual representation of the effects of waterhammer and it should be studied for a deeper understanding of the phenomenon, even if it is not used anymore for the evaluation of real cases [8, 9, 27, 28]. As the theoretical description might be confused, the method will be explained solving a simple problem. Again, the usual very simple plant reproduced in Figure 3.6 is considered, the characteristics are those used in the example in Chapter 2: pipe length L = 1000 m where the wave celerity is equal to c = 1000 m/s and the flow has a velocity equal to V0 = 5 m/s. Pressures have to be computed in both cases of valve closure with


37

Upstream reservoir

Downstream valve

Figure 3.6: A simple plant to be solved with the graphical method.

7 V1

6 5 4 3 2 R2

1 0 0

0.2

0.4

0.6

0.8

1

V3

Figure 3.7: Graphical solution for sudden valve closure.

time Tc = 0.0 s (sudden valve closure) and Tc = 5.0 s. Headlosses are considered negligible. The dimensionless variables are defined as: h = H /H0 , v = V /V0 , t = T /Tc . For the sudden valve closure, the following points have to be computed and placed on a graph which x-axis is v and y-axis is h, as reported in Figure 3.7. The point V0 refers to the steady flow conditions (time t = 0) at the valve V , and therefore has co-ordinates (v = 1; h = 1). Then we have to move on that graph

38 WATER HAMMER SIMULATIONS following two lines whose slope S is equal to those of the characteristics, i.e.: S=±

c · V0 1000 × 5 = ±5.1 =± g · H0 9.81 × 100

(3.42)

Starting from the point at the valve V0 and moving with slope –S (negative characteristic line) we touch the axis v = 0 (which is reached at time t = 1) with the co-ordinate h = 6.1; this is the point V1 (0, 6.1) related to the time when the valve is closed. Because of the hypothesis of negligible headlosses, the point R0 at the reservoir R at time t = 0 is V0 ≡ R0 . The next point is related to the upstream reservoir: starting from V1 and moving on the positive characteristic line, i.e. with slope S, until the intersection h = 1, because the reservoir keeps the head constant; therefore we reach the point R2 (−1, 1). Continuing this procedure the co-ordinates of V3 (0, −4.1) and V4 ≡ V0 are computed. As the headlosses are considered negligible, the oscillations do not reduce and the cycle continues as it has been computed so far. With regards to the slow closure, a timestep of 1 s is chosen, which, being Tc = 5.0 s, in dimensionless form is t = 0.2. In the case of slow closure, the law to be used is: v = (1 − t)n

(3.43)

where n is an exponent that can be assumed equal to 1 for simplicity, but reminding that this choice does not change the application of the procedure. The method is reported in Figure 3.8: again the starting point is V0 (1, 1), and again we have to move on the negative characteristic line to find the point V1 . This point is found at the position computed with the (3.43) and therefore v = 0.8. From the point V1 , moving on the negative characteristic we arrive to the point V2 to be reached at time t = 0.4 and therefore v = 0.6. At the same time, moving on the positive characteristic, the point R2 is found, which has also to be at the intersection of h = 1 because the point R is at the reservoir, and therefore, as known, with constant head. Continuing this procedure, point V3 can be positioned and so on until the requested time of simulation is reached. The points computed with this procedure are shown in Figure 3.8. This solution indicates that the maximum peak of pressure is reached at time t = 0.4 (in the dimensional values this is equal to 2 s) and it has a value that can be estimated equal to 3; dimensionally, this means that the peak has to be expected around 300 m. After this high peak, the oscillations are repeated with smaller values, graphically estimated equal to 2, and therefore around 200 m. Moreover, after the complete closure of the valve, t > 1 and the values of R are on the horizontal line h = 1 as usual, and oscillate between v = −0.2 and 0.2; the values of V are on the vertical line v = 0 and oscillate between h = 0.0 and 2.0. With regard to the effect of the headlosses, they can be considered but using an expedient, which consists in applying them at the upstream end of the pipe. Therefore, the boundary upstream condition is not represented by an horizontal

39

SIMPLIFIED SOLUTIONS 7

t 1.0

t 0.4

t 0.6

t 0.8

t 0.0

t 0.2

6

Ch

ara

5

h

4

cte

rist

ic l

ine

cV2

3 V5 V6

V3

2 R R 6 7 1

V1

V4 R3 R4 R5 h 1

V0 R0 R1

R2

0 V7 1 0.2

0

0.2

0.4 v

0.6

0.8

1

Figure 3.8: Graphical solution for slow valve closure. straight line h = 1.0 but by a parabola given by:

λ·L h=1− 1+ D

·

V02 · v 2 2 · g · H0

Under all the other aspects, the procedure remains unchanged.

(3.44)


Chapter 4 Numerical solution of the governing equations: The method of characteristics

There are many methods to numerically integrate differential equations, which can be roughly divided in three main families: finite differences, finite elements and characteristics. One of the first integration methods used to solve waterhammer problems is the method of characteristics applied by Evangelisti to the case of fluid transients [5, 6] following a more general development [29, 30]. This method has proved to be easy, fast and reliable. Actually, many problems cannot be faced with this integration method, but it surely deserves a description, especially because its worldwide diffusion. The method of characteristics is a technique for solving partial differential equations. Typically, it is applied to solve first-order equations, although more generally it is valid for any hyperbolic partial differential equation. The method is to reduce a partial differential equation to a family of ordinary differential equations along which the solution can be integrated from some initial data given on a suitable hypersurface.

4.1 Numerical solution A linear combination of the partial derivative of a given function u (x, t), for instance: ∂u ∂u +b· (4.1) a· ∂x ∂t is the derivative of the function in the direction: a dx = dt b

(4.2)

The problem is the identification, on the plane (x, t), of particular directions along which in the governing equations (continuity: eqn (2.34) and motion: eqn (2.43)) the functions p(x, t) and V (x, t) are derived, in order to reduce the partial derivative differential equations to ordinary differential equations. These directions are called characteristics and, if they are real and distinct, the system is said to be hyperbolic.

42 WATER HAMMER SIMULATIONS To identify these directions, the two equations are linearly combined, multiplying the momentum equation by a constant K and summing the result to the continuity equation [10, 11]. The result is: ∂p ∂V ∂V dz +K ·ρ· +K ·ρ·V · +K ·γ ·J +K ·γ · ∂x ∂t ∂x dx ∂p ∂V +V · + ρ · c2 · =0 ∂x ∂x

K·

(4.3)

upon collecting the common terms, it yields: (K + V ) ·

∂p ∂p ∂V dz ∂V + + ρ · (K · V + c2 ) · +ρ·K · +γ · K · J + =0 ∂x ∂t ∂x ∂t dx (4.4)

The coefficients of the derivative of p and V with respect to x and t are: ∂p ∂t ∂p ∂x ∂V ∂t ∂V ∂x

→ 1

(=a)

(4.5)

→ K +V

(=b)

(4.6)

→ K ·ρ

(=c)

(4.7)

→ K · ρ · V + ρ · c2

(=d)

(4.8)

Therefore, p and V are derivative along the same direction if a/b = c/d, and consequently if: 1 ρ·K = K +V ρ · (K · V + c2 )

⇒

K +V =V +

c2 K

(4.9)

i.e., if K = ±c. As mentioned, the system is hyperbolic because the two characteristics are real and distinct, and they are: dx =V +c (4.10) dt dx =V −c dt

(4.11)

It can be noted that the velocity V is normally much lower than 10 m/s, while the wave celerity is normally in the range 800 ÷ 1200 m/s. As a consequence, an acceptable approximation is V + c ≈ c and V − c ≈ −c.

NUMERICAL SOLUTION OF THE GOVERNING EQUATIONS

43

The governing equations can be rewritten from eqn (4.4), inserting K = +c and −c. Reminding that: dp =

∂p ∂p · dt + · dx ∂t ∂x

(4.12)

dV =

∂V ∂V · dt + · dx ∂t ∂x

(4.13)

and being in this case:

dp = dV =

∂p ∂p ±c· · dx ∂t ∂x

∂V ∂V ±c· ∂t ∂x

(4.14)

· dx

Writing again eqn (4.4) yields: ∂p ∂V dz ∂V ∂p ±c· ±ρ·c· ±c· ±c·γ ·J ±c·γ · =0 ∂t ∂x ∂t ∂x dx

(4.15)

(4.16)

or: dV dp dz ±ρ·c· ±c·γ ·J ±c·γ · =0 dt dt dx Writing the piezometric head h instead of the pressure p, being: p h=z+ γ it finally results: ∂p =γ · ∂x

∂h ∂z − ∂x ∂x

+

∂γ · (h − z) ∂x

(4.17)

(4.18)

(4.19)

In eqn (4.19) a good approximation is: ∂γ ≈0 ∂x

(4.20)

Moreover: ∂p ∂h =γ · (4.21) ∂t ∂t Inserting these expressions in eqn (4.16), the following system of two ordinary differential equations is obtained, which is equivalent to eqns (2.34) and (2.43): ⎧ dh c dV ⎪ ⎪ ⎨ dt + g · dt + c · J = 0 (4.22) ⎪ ⎪ ⎩ dh − c · dV − c · J = 0 dt g dt

44 WATER HAMMER SIMULATIONS where the first equation is valid along the positive characteristic line (4.10) and the second is valid along the negative characteristic line (4.11). Expressing the headlosses with the uniform flow formula: J =λ·

V · |V | 2·g·D

and writing eqn (4.22) with finite differences, it yields: ⎧ h c V V · |V | ⎪ ⎪ ⎪ ⎨ t + g · t + c · λ · 2 · g · D = 0 ⎪ h c V V · |V | ⎪ ⎪ − · −c·λ· =0 ⎩ t g t 2·g·D

(4.23)

(4.24)

which can be written using an explicit method: ⎧ Vi−1, j−1 · Vi−1, j−1 hi, j − hi−1, j−1 c Vi, j − Vi−1, j−1 ⎪ ⎪ + · +c·λ· =0 ⎪ ⎨ t g t 2·g·D (4.25) ⎪ Vi+1, j−1 · V ⎪ h − h − V V c i+1, j−1 i, j i, j i+1, j−1 i+1, j−1 ⎪ ⎩ =0 − · −c·λ· 2·g·D t g t This is an algebraic system of two equations in the two unknown hi, j and Vi, j that can be easily solved, giving: 1 c hi, j = · hi+1, j−1 + hi−1, j−1 − · Vi+1, j−1 − Vi−1, j−1 + 2 g c·λ · Vi−1, j−1 · Vi−1, j−1 − Vi+1, j−1 · Vi+1, j−1 · t − 2·g·D g c Vi, j = (4.26) · hi+1, j−1 + hi−1, j−1 + · Vi+1, j−1 + Vi−1, j−1 + 2·c g c·λ − · Vi−1, j−1 · Vi−1, j−1 + Vi+1, j−1 · Vi+1, j−1 · t 2·g·D Graphically, the characteristic lines may be represented on a (x, t) plane, where a point has co-ordinated x = i · x and t = j · t and the variables h and V in this point are indicated as hi, j and Vi, j . In Figure 4.1 the (x, t) plane with the symbols mentioned is shown.

4.2 Initial and boundary conditions With eqns (4.26) it is possible to compute the values of the required parameters in the internal points of the field, when the initial and boundary conditions are known. As a differential equations problem, in fact, these conditions must be defined.


45

t

x

Figure 4.1: (x, t) plane with the representation of a point. As mentioned in Chapter 1, initial conditions are normally computed with the usual methods of steady flow in water supply networks, under the hypothesis that before the transient these methods can be applied. With regard to the boundary conditions, these depend on the device positioned at the extreme ends of the pipe, which determine the behaviour of the system. A specific boundary condition settles one of the parameters (head h or velocity V ), or ties these two parameters with a definite function. The two parameters h and V at each end of the pipe are computed solving the equation(s) which describe the boundary condition together with the characteristic line, for downstream or upstream conditions, respectively. 4.2.1

Reservoir

The simplest case that can be handled is probably the reservoir, under the hypothesis that it is large enough to keep its level unchanged during the transient: in this case, the head can be considered constant. To this very simple condition, a characteristic equation has to be associated. Let us assume the pipe is divided into n + 1 sections, from 0 to n and the reservoir is positioned upstream (i.e., in the section 0), it yields: h(0, t) = h0 V (0, t + 1) = V (1, t) +

∀t

the head in the reservoir is constant

(4.27)

g V (1, t) · |V (1, t)| · [h0 − h(1, t)] − λ · · t c 2·D

(4.28)

Very similar is the case of a reservoir positioned downstream (i.e., in the section n), so that to the condition that impose the constant head has to be associated the positive characteristic line; it yields: h(n, t) = h0

∀t

the head in the reservoir is constant

(4.29)

46 WATER HAMMER SIMULATIONS c

Reservoir upstream

c

Reservoir downstream

Figure 4.2: Scheme for the reservoir (upstream and downstream) on the (x, t) plane. g V (n − 1, t) · |V (n − 1, t)| · [h0 −h(n − 1, t)] − λ · · t c 2·D (4.30) The scheme of the method is reported in Figure 4.2. V (n, t +1)=V (n−1, t)−

4.2.2 Valve Valves are very complex devices, but in this section the simplest description is given and implemented: the hypothesis is that to a certain opening of the valve corresponds an univocal value of the flow velocity. In this case, eqn (2.11) seen in Section 2.3 can be used again, even if given more generally as follows, to allow its use even for nonlinear operations: ⎧ t n ⎨ V (0, t) = V0 · 1 − Tc ⎩ V (0, t) = 0

when t ≤ Tc

(4.31)

when t > Tc

In this case, the valve is positioned upstream, i.e., in the section 0. Therefore, the characteristic equation to be associated is negative and yields: c c V (1, t) · |V (1, t)| · [V (0, t + 1) − V (1, t)] + λ · · · t g g 2·D (4.32) Again, very similar is the case of the valve positioned downstream, where the condition on the velocity has to be associated to the positive characteristic line. h(0, t + 1) = h(1, t) +

⎧ t n ⎨ V (n, t) = V0 · 1 − Tc ⎩ V (n, t) = 0

when t ≤ Tc when t > Tc

(4.33)


h(n, t + 1) = h(n − 1, t) − −λ · 4.2.3

47

c · [V (n, t + 1) − V (n − 1, t)] g

c V (n − 1, t) · |V (n − 1, t)| · · t g 2·D

(4.34)

Junction

This case is quite interesting as it allows the understanding of the characteristic lines, visualizing schematically the computing diagram. Let us suppose the simplest case, when two pipelines meet in a junction. The indexes 1 and 2 will be used for the upstream and downstream pipes, respectively. The upstream pipe is divided into n + 1 sections, from 0 (upstream) to n (downstream); the downstream pipe is divided in m + 1 sections numbered from n to n + m. Obviously, the section n, where the junction is positioned, belongs to both the upstream and downstream pipes; however, it will be divided in nup and ndown as the flow parameters h and V are, generally speaking, different. Should the concentrated headlosses in the junction be negligible, the heads have to be equal, and therefore: h(nup , t + 1) = h(ndown , t + 1)

(4.35)

In general cases, headlosses are not negligible, and therefore a further term is to be added to eqn (4.35), normally a function of the velocity V . The continuity equation is: V (nup , t + 1) ·

2 π · Dup

4

= V (ndown , t + 1) ·

2 π · Ddown 4

(4.36)

To these two equations, the two characteristic lines (one positive for the downstream pipe and one negative for the upstream) have to be associated, forming a system of four equations in the four unknown flow parameters hup , hdown , Vup and Vdown . As working with this numerical method the condition t = x/c must be satisfied, the assumption of generic spatial discretization steps xup and xdown is shown in Figure 4.3. The use of arbitrary spatial steps brings to different time steps, which force the user to interpolate the carried out values at the junction to compute the flow parameters h and V at the same time. This is possible, but it introduces approximation errors that decrease the precision of the results; moreover, interpolation algorithms must be implemented in the computer code, slowing the solution: with this solution, the advantages given by the method of characteristics, precision and velocity, have to be partially discarded.

48 WATER HAMMER SIMULATIONS t

∆t2 ∆t1

x 0 nm

n

Figure 4.3: x–t plane for two pipes with different characteristics and “problem” at the junction. As a consequence, working with the method of characteristics the condition that the timestep t be the same for all the pipes is normally imposed. That means: t =

x2 x1 = c1 c2

(4.37)

which is: l1 l2 = n · c1 m · c2

(4.38)

Again a new problem arises: as the parameters n and m must be integer, it may be necessary to slightly change the geometric characteristics of one of the pipes (upstream or downstream). If a new length for the downstream pipe l2 is to be computed, that can be done with the following: l2 = m · c2 ·

l1 n · c1

(4.39)

This new value is obviously different from the real length, but simple computations show that if the pipelines are long few thousands meters, and if their lengths do not differ too much, the difference between the original and the new computed lengths is equal to few tens of centimeters, and therefore completely insignificant on the computed results. A non-negligible problem can arise when a network with very different pipe lengths has to be simulated. In this case, the condition (4.38) may drive to a very small timestep and a huge number of sections in the longest pipes, which are not necessary and that highly increase the computational effort. In these cases, the use of an interpolation algorithm can be suggested, because the slowing down of the program is matched by the possibility to have an higher timestep and a reduced section number in the longest pipes. The generalization for a higher number of pipes is possible, as for each pipe two new unknowns are set up (the head and the velocity in the junction for the pipeline


49

ReadNetworkCharacteristics; ReadInitialConditions; t=0 while t < tmax do begin t: = t + dt; ComputeInternalPoints; ComputeBoundaryConditions; end; Figure 4.4: Main procedure of the solution program.

Procedure ComputeInternalPoints for i:=1 to n-1 do begin vp[i]:=0.5*(v[i-1]+v[i+1]+(g/cel)*(h[i-1]-h[i+1])(Lambda*dt/(2*diam))*(v[i-1]*abs(v[i-1])+v[i+1]*abs(v[i+1]))); hp[i]:=0.5*(h[i-1]+h[i+1]+(v[i-1]-v[i+1])/(g/cel)(Lambda*dt/(2*diam))*(v[i-1]*abs(v[i-1])-v[i+1]*abs(v[i+1]))/(g/cel)); end; Figure 4.5: Core procedure of the solution program: calculation of the unknowns (head and velocity) for the points internal to the pipe. end) and two new equations are added to the system: the equality of the heads, similar to eqn (4.35) and the characteristic line; moreover, the continuity equation (4.36) must be changed accordingly.

4.3 The computer code In this chapter, a simple computer program is shown and applied to simulate a pipe with a reservoir positioned upstream and a simple valve positioned downstream; the program has the simple structure reported in Figure 4.4. The core procedures have the structures reported in Figures 4.5 and 4.6; as it can be seen, they exactly match eqns (4.26) for the internal points (Figure 4.4), eqns (4.27) and (4.28) for the upstream reservoir and eqns (4.33) and (4.34) for the downstream valve (Figure 4.5). The parameters to be used (heads and velocities) are stored in different arrays: h[1 . . . NumberOfSections] and v[1 . . . NumberOfSections]. The unknowns, i.e., the parameters at time t + dt, are stored in the arrays hp[1 . . . NumberOfSections] and vp[1 . . . NumberOfSections].

50 WATER HAMMER SIMULATIONS Procedure ComputeBoundaryConditions; // Reservoir positioned upstream hp[0]:=h0; vp[0]:=v[1]+g*(h[0]-hp[1])/cel-Lambda*v[1]*abs(v[1])*dt/(2*diam); // Valve positioned downstream if t
4.4 First simple application The simplest example is given by a pipe which connects an upstream reservoir with a downstream valve. The pipe has a diameter equal to 75 mm, length of 800 m and the headlosses are negligible. The wave celerity has been computed equal to 1000 m/s and the initial discharge is equal to 2.5 l/s. The downstream valve has a sudden closure. Using the described computer program, the inserted parameters are presented in Figure 4.7, which shows the user’s interface. As for the simulation parameters, the number of computational sections has been let equal to 1000, which is very high, and the effect will be discussed later; the checkbox “Ignore Headlosses” has been thicked and the computation time set equal to 15 s. Results are printed in two text files and imported in a spreadsheet. As the program can produce a very large amount of data, they might be reduced by printing them less frequently, fixing a larger time step, in the example let equal of 0.1 s; moreover, results are printed only in few sections, as mentioned in Chapter 1, at the upstream and downstream ends, and at the 20, 40, 60 and 80% of the length of the pipe, and, therefore, in this case at 160, 320, 480 and 640 m from the edges of the pipe.


51

Figure 4.7: User’s interface of the program. Very simple hand calculation show that the spatial discretization step and the time step are: ds =

PipeLength 800 = = 0.8 m NSections 1000

dt =

ds 0.8 = = 0.001 s celerity 1000

The velocity in the pipe, in steady flow (initial) conditions is equal to 0.566 m/s. As the valve suddenly closes, the surcharge can be computed with the Allievi– Joukowsky’s formula (2.4), which yields: h =

c 1000 · V0 = · 0.566 ∼ = 57.7 m g 9.81

As the wave celerity is equal to 1000 m/s the phase time of the system, i.e., the time, the wave requires to reach the reservoir and to be back to the valve is equal to 2 · L/c = 1.6 s. The results of the simulation carried out with the computer program are shown in Figure 4.8; as can be seen, these results are perfectly consistent with the hand calculations.

52 WATER HAMMER SIMULATIONS 180 160 140

Head (m)

120 100 80 60 40 20 0 0.0

1.5

3.0

4.5

6.0

7.5 9.0 Time (s)

10.5

12.0

13.5

15.0

Figure 4.8: Results of the first simulation: sudden closure and negligible headlosses. When a slower operation is performed at the downstream valve, results are expected to change. In particular, the second simulation is performed considering the valve closes in 3 s. In this case, the Michaud’s formula (2.16) has to be used, which yields: h =

2 × 800 × 0.566 2 · L · V0 = = 30.8 m g · Tc 9.81 × 3

As described in Chapter 2, the maximum head is reached at time 2 · L/c = 1.6 s; results are shown in Figure 4.9. These results can be also compared with those carried out with the graphical method in Chapter 3. In order to check the effect of the headlosses, a third simulation has been performed, letting again the valve closure time equal to zero and the roughness n = 0.0125 m−1/3 s. Results are shown in Figure 4.10 and are interesting because of two aspects. The former is that the maximum value is not reached immediately. The reason of this phenomenon can be found in Section 2.4, where it is explained that at the first timestep only the first element is stopped, but the surcharge reached is not enough to stop the whole water column, and therefore it stops in a larger time. The latter aspect is due to the reduction of the amplitude of the oscillations: again this is due to the headlosses, that induce a loss of energy (head) in the water column, so that the surcharge fades with time. The last simulation is carried out with the same parameter of the first, i.e., valve closure time equal to zero and negligible headlosses, but with a smaller number of section, let equal to 10; moreover, the time of simulation is equal to 30 s. Results are


53

140 120

Head (m)

100 80 60 40 20 0 0.0

1.5

3.0

4.5

6.0

7.5

9.0

10.5

12.0

13.5

15.0

Time (s)

Figure 4.9: Results of the second simulation: slow closure and negligible headlosses.

180 160 140

Head (m)

120 100 80 60 40 20 0 0.0

1.5

3.0

4.5

6.0

7.5

9.0

10.5

12.0

13.5

15.0

Time (s)

Figure 4.10: Results of the third simulation: sudden closure and not negligible headlosses.


Head (m)

120 100 80 60 40 20 0

0

5

10

15

20

25

30

Time (s)

Figure 4.11: Results of the fourth simulation: sudden closure and negligible headlosses; the number of computational sections is quite low and the effects of numerical diffusion are significant. reported in Figure 4.11. As can be seen, the angles of the square wave obtained with the first simulation (in Figure 4.7) become smoother and the amplitude decreases; this phenomenon is called diffusion. In this case, the phenomenon is not physical, and it is due only by numerical errors (for this reason it is called numerical diffusion), which will be discussed in detail in Chapter 5, being introduced here only because the reader could find it by “playing” with the enclosed computer code, or with other codes available.

Chapter 5 Numerical solution of the governing equations: finite difference methods

More flexible methods are those based on finite differences [31]. The idea beyond these methods is very simple, being based on the substitution of the differential ratios with finite differences. Obviously, while analytical solutions are defined continuously on the whole domain, numerical solutions are defined only in a finite number of points. The values computed on this finite number of points can be different from those which would be carried out with an analytical solution, but they should converge when the spatial integration step x decreases. The unknown parameters are computed starting from those already known, i.e., from the values the parameters have at a former time. To this end, these methods can be divided into explicit and implicit. Explicit methods solve the differential equations at time t + 1 using data related only to previous times (Figure 5.1). In Figure 5.1, the derivatives ∂/∂x and ∂/∂t are computed as function of the unknown variable values at time t + 1 and position i and of the known variable values at time t and positions i − 1, i and i + 1. The solution of the equations, therefore, advances one point at a time for each timestep, and consequently each equation contains only one unknown, which can be easily computed. On the other hand, implicit methods solve the differential equations at time t + 1 using data related both to previous times and to the same time t + 1 (Figure 5.2). As the latter data are unknown, the equations have to be solved simultaneously, and this requires the solution of a system which is usually quite large; for this reason implicit methods are computationally more demanding. The advantage of the implicit methods is that they are unconditionally stable, and that means the timestep may be larger. Explicit methods, instead, have to satisfy the (necessary, but not sufficient) stability criterion: to this end, small timesteps have to be chosen and, therefore, the time of simulation increases. However, in having a too large timestep other problems may arise, because the truncation error increases and the precision of the solution diminishes. The approximation of a partial derivative with an algebraic ratio is possible as it can be shown starting from the Taylor series, which allows the approximation of any function with a polynomial.

56 WATER HAMMER SIMULATIONS t

1

2

3

4

5

6

7

1

2

3

4

5

6

7

x

Figure 5.1: Explicit method – example.

t

1

2

3

4

5

∆t

6

7

6

7

∆x 1

2

3

4

5

x

Figure 5.2: Implicit method – example. Let f (x), a real numerical function, defined in [a, b] and let x0 ∈ [a, b]. If f can be derived, the Taylor’s series is: ∂f ∂ 2 f (x − x0 )2 f (x) = f (x0 ) + · (x − x ) + · + ... 0 ∂x x0 ∂x2 x0 2! ∂ n f (x − x0 )n + n · +H ∂x x0 n!

(5.1)

where H is the truncation error. The first derivative in a point is the slope of the tangent to the curve that can be approximated by the straight line passing for x0 and x0+1 : in this case, the approximation is called forward difference; when the points are x0−1 and x0 , the approximation is called backward difference. Different combinations are possible, for instance when the derivative is approximated using points x0−1 and x0+1 , the approximation is called central difference.


57

Expressing these concepts with formulae, the result is: f (xi+1 ) − f (xi ) ∂f = ∂x xi x ∂f f (xi ) − f (xi−1 ) = ∂x xi x ∂f f (xi+1 ) − f (xi−1 ) = ∂x xi 2 · x

forward difference

backward difference

central difference

For the second-order derivative, the following expression can be used: ∂ 2 f f (xi+1 ) − 2 · f (xi ) + f (xi−1 ) = 2 ∂x xi (x)2 The terms which have been neglected are called truncation error, which is a measure of the accuracy of the solution and which diminishes when the distance between points decreases.

5.1 The Courant–Friedrichs–Levy stability condition Very often in the context of hyperbolic systems, the solutions are required for a long duration; in these cases, a strong stability is required, in order to guarantee the solution is limited for any value at any time. In order to have a stable solution, it is necessary the errors due to the discretization are not amplified from each temporal line to the following. The stability criterion, which, as mentioned, is necessary but not sufficient, is the following [32]: t =λ (5.2) x This condition is known as Courant–Friedrichs–Levy (CFL) and it is based on the concept of “dependency domain”, which is defined as the set of points that influence the value of the solution in a generic point (xj , tn ). Referring to Figure 5.3, the area inside the triangle with vertices are (xj , tn ), (xj−n , t0 ) and (xj+n , t0 ) is the dependency domain. The lines drawn in the same figure are the characteristics of the differential equations, as described in Chapter 4. The stability condition requires the physical dependency domain to be inside the numerical dependency domain. If this condition is not satisfied, the numerical solution cannot converge towards the physical solution, because the change of a condition within the physical dependency domain but not within the numerical dependency domain would not have effect on the numerical solution. From the mathematical point of view, thus, it is required that the spatial and temporal steps respect the relation (5.2).

58 WATER HAMMER SIMULATIONS (xj , tn)

tn

t0 (xjn, t0)

(xjn, t )

xj

0

Figure 5.3: Dependency domain. t

t xj, tn

xj, tn

0

0 x

x

Figure 5.4: Physical and numerical dependency domain. On the left, the CFL condition is respected, as the numerical domain contains the physical one; on the right, the CFL condition is not satisfied. As the characteristic lines are the following (cfr. Chapter 4): dx = V ± c ≈ ±c dt the limitations to the timestep is given by: c · t x ≤ 1

(5.3)

(5.4)

This condition is shown in Figure 5.4. As mentioned, this condition may be not sufficient to assure stability, especially when the boundary conditions are improperly-posed.

5.2 The Lax–Wendroff method This method [33] is explicit and it is widely used for the solution of linear hyperbolic differential equations with central differences. It is second order accurate both in space and time [34] and it is obtained taking into account the first three terms of the Taylor’s series. Being an explicit method, the CFL condition has to be


59

respected: in this case, it can be demonstrated that the method converges towards the exact solution. A linear system of differential equations can be written as: →

→

∂u ∂u +A· =0 ∂t ∂x

(5.5)

ut + A · u x = 0

(5.6)

or, in compact form:

As mentioned, at the base of the Lax–Wendroff method is the Taylor’s series: u(x, tn+1 ) = u(x, tn ) + t · ut (x, tn ) +

1 · (t)2 · utt (x, tn ) + · · · 2

(5.7)

Equation (5.6) yields: ut = −A · ux

(5.8)

and deriving again with respect to time: utt = −A · uxt = A2 · uxx

(5.9)

uxt = utx = (−A · ux )x = −A · uxx

(5.10)

Equation (5.6) implies:

Substituting these expressions of ut and utt into eqn (5.7) it yields: u(x, tn+1 ) = u(x, tn ) − t · A · ux (x, tn ) +

1 · (t)2 · A2 · uxx (x, tn ) + O(x3 , t 3 ) 2 (5.11)

Neglecting the truncation error O(x3 , t 3 ) and approximating ux and uxx with central differences, the Lax–Wendroff method is obtained, which can be written as: ujn+1 = ujn − A ·

λ λ2 n n n n · (uj+1 · A2 · (uj+1 − uj−1 )+ − 2ujn + uj−1 ) 2 2

(5.12)

where j and n indicate the spatial and temporal positions, respectively, and λ = t/x. This method is defined as three-point support scheme: examining the points of the grid which have been used in the computation of the new value, the stencil of the method can be drawn (Figure 5.5) and, as can be seen, to compute the values in the point ( j, n + 1) it is necessary the knowledge of the values in the nodes ( j, n), ( j + 1, n) and ( j − 1, n).

60 WATER HAMMER SIMULATIONS Time

n1

n

j1

j1

j

Space

Figure 5.5: Stencil for Lax–Wendroff method. Unfortunately, waterhammer is not a linear problem, and the system to be solved is not precisely defined with eqn (5.5) but instead with: →

→

→ ∂u ∂u + A (u) · =E ∂t ∂x

(5.13)

However, Lax–Wendroff method can be extended to quasi-linear problems [35] with the expression: ujn+1 = ujn − Anj ·

λ λ2 n n n n · (uj+1 − uj−1 )+ · (Anj )2 · (uj+1 − 2ujn + uj−1 ) 2 2

(5.14)

5.3 Solving the governing equations The method has to be applied to the equations governing the waterhammer phenomenon, written in conservative form as follows: ∂V ∂p ∂p +V · + ρ · c2 · =0 ∂t ∂x ∂x

(5.15)

∂V γ ·J γ ∂z 1 ∂p ∂V + + ·V · + + · =0 ρ ∂x ∂t ∂x ρ ρ ∂x

(5.16)

Substituting in eqns (5.15) and (5.16) the following expression: ∂h ∂p =γ · ∂t ∂t ∂p =γ · ∂x

∂h ∂z − ∂x ∂x

(5.17) (5.18)


61

Neglecting the term V · ∂z/∂x as it is very small, the equations become: ∂h ∂h c2 ∂V +V · + · =0 ∂t ∂x g ∂x

(5.19)

∂V ∂h ∂V +g· +V +J ·g =0 ∂t ∂x ∂x

(5.20)

or, in matrix form: ∂ h V + ∂t V g

∂ h 0 · = V ∂x V −g · J c2 g

(5.21)

The system to be solved can be finally written as: →

→

→ ∂u ∂u + A(u) · =E ∂t ∂x

where:

h u= V

V A= g

V c2 g

0 E= −g · J

(5.22)

(5.23)

Again, the headlosses can be computed with the Chézy’s formula: J =λ·

V · |V | 2·g·D

(5.24)

The application of the Lax–Wendroff ’s formula (5.14) yields: λ λ2 n n n n − uj−1 )+ − 2ujn + uj−1 ) + t · Ejn · (uj+1 · (Anj )2 · (uj+1 2 2 (5.25) which can be written as: n hn+1 hn λ V c2 n hj+1 − hnj−1 j j g n+1 · uj = n+1 = n − · n n Vj Vj 2 g V Vj+1 − Vj−1 j ⎛ ⎞2 n c2 n 0 n hj+1 − 2hnj + hnj−1 λ2 ⎝ V g ⎠ + t · · n + · n n 2 −J · g Vj+1 − 2Vj + Vj−1 g V j ujn+1 = ujn − Anj ·

j

(5.26) where:

2 2·V ·c2 2 + c ) (V g 2 A = 2 · V · g (V 2 + c2 )

(5.27)

62 WATER HAMMER SIMULATIONS Separating the two unknown parameters V and h, the following solving equations are obtained: Vjn+1 = Vjn −

λ n n · (g · (hnj+1 − hnj−1 ) + Vjn · (Vj+1 − Vj−1 )) 2

λ2 · ((2 · Vjn · g) · (hnj+1 − 2 · hnj + hnj−1 ) 2 n n + (c2 + (Vjn )2 ) · (Vj+1 − 2 · Vjn + Vj−1 )) − g · t · Jjn +

λ c2 n n Vjn · (hnj+1 − hnj−1 ) + − Vj−1 ) · (Vj+1 2 g 2 λ + · ((Vjn )2 + c2 ) · (hnj+1 − 2 · hnj + hnj−1 ) 2 2 · Vjn · c2 n n n + · (Vj+1 − 2 · Vj + Vj−1 ) g

(5.28)

hn+1 = hnj − j

(5.29)

Equations (5.28) and (5.29) are those implemented in the computer program.

5.4 Boundary conditions As can be easily seen from the stencil of the Lax–Wendroff method shown in Figure 5.5, at boundaries this method cannot be applied. In fact, upstream, to compute the values at time n + 1 in the point j = 0, the solver should use the values at time n in the point j = −1, j = 0 and j = 1; the point j = −1 does not exist; the same applies for the downstream end of the pipe. There are two main ways to overcome this problem: (1) the use of a different scheme and (2) the definition of “virtual” points (called ghost cells) outside the integration space, the characteristics of which will be discussed in section 5.4.2. 5.4.1 Asymmetrical schemes A number of schemes useful to solve the problem at boundaries can be found in the literature; not to be forgotten is the possibility to use the same method of characteristics which has been presented in Chapter 4. However, as the goal of this chapter is the use of finite differences, the very simple method called Upwind is presented and implemented in the computer code. The stencils for the method are shown in Figure 5.6. According to the stencils shown in Figure 5.6, the expressions for the Upwind method upstream and downstream are, respectively: ux (xj , tn ) =

1 · [u(xj+1 , tn ) − u(xj , tn )] + O(x; t) 2 · x

(5.30)

63

NUMERICAL SOLUTION OF THE GOVERNING EQUATIONS n-time

n-time

1

1

0

0

1

1 2

1

0 i-space

1

2

2

1

0 i-space

1

2

Figure 5.6: Stencil for the Upwind method. Upstream (left) and downstream (right). ux xj , tn =

1 · [u(xx , tn ) − u(xj−1 , tn )] + O(x; t) 2 · x

(5.31)

Substituting in eqn (5.25) it yields: ujn+1 = ujn −

A n − ujn ) · λ · (uj+1 2

(5.32)

ujn+1 = ujn −

A n · λ · (ujn − uj−1 ) 2

(5.33)

Substituting in eqns (5.32) and (5.33) the expressions to be inserted in the computer program can be obtained, as follows: = hn0 − hn+1 0

V0n+1 = V0n −

λ · [g · (hn1 − hn0 ) + V0n · (V1n − V0n )] − g · t · J0n 2

hn+1 = hnn − n

Vnn+1 = Vnn −

λ c2 · [V0n · (hn1 − hn0 ) + · (V1n − V0n )] 2 g

λ c2 n · [Vnn · (hnn − hnn−1 ) + · (Vnn − Vn−1 )] 2 g

λ n )] − g · t · Jnn · [g · (hnn − hnn−1 ) + Vnn · (Vnn − Vn−1 2

(5.34)

(5.35)

(5.36)

(5.37)

Equations (5.34) and (5.35) are valid for the upstream conditions, while eqns (5.36) and (5.37) are valid for the downstream conditions. Obviously, more complex schemes can be developed. As example, in Figure 5.7 the Beam–Warming scheme is presented.

64 WATER HAMMER SIMULATIONS n-time

n-time

1

1

0

0

1

1 2

1

0 i-space

1

2

2

1

0 i-space

1

2

Figure 5.7: Stencil for the Beam–Warming method. Upstream (left) and downstream (right). For the Beam–Warming method, the derivatives, according with the stencil, for the upstream condition are: ux (xj , tn ) =

1 · [3u(xx , tn ) − 4u(xj+1 , tn ) + u(xj+2 , tn )] + O(x2 ) 2 · x

(5.38)

1 · [u(xx , tn ) − 2u(xj+1 , tn ) + u(xj+2 , tn )] + O(x) (x)2

(5.39)

uxx (xj , tn ) =

while for the downstream condition are: ux (xj , tn ) =

1 · [3u(xx , tn ) − 4u(xj−1 , tn ) + u(xj−2 , tn )] + O(x2 ) 2 · x

(5.40)

1 · [u(xx , tn ) − 2u(xj−1 , tn ) + u(xj−2 , tn )] + O(x) (x)2

(5.41)

uxx (xj , tn ) =

Again the insertion of these expressions into the eqn (5.25) yields: ujn+1 = ujn + A ·

λ λ2 n n n n + uj+2 )+ + uj+2 ) (5.42) · (3ujn − 4uj+1 · A2 · (ujn − 2uj+1 2 2

ujn+1 = ujn − A ·

λ λ2 n n n n + uj−2 )+ + uj−2 ) (5.43) · (3ujn − 4uj−1 · A2 · (ujn − 2uj−1 2 2

The equations to be inserted in the computer program can be obtained as follows: c2 n n λ n n n n n n hn+1 + · V · 3V = h + · 3h − 4h + h − 4V + V 0 0 0 1 2 0 1 2 0 2 g λ2 2 · V0n · c2 n 2 2 n n n n n n + · ((V0 ) + c ) · (h0 − 2h1 + h2 ) + · (V0 − 2V1 + V2 ) 2 g (5.44)


V0n+1 = V0n + +

65

λ · [g · (3hn0 − 4hn1 + hn2 ) + V0n · (3V0n − 4V1n + V2n )] 2

λ2 · [2 · g · V0n · (hn0 − 2hn1 + hn2 ) + ((V0n )2 + c2 ) · (V0n − 2V1n + V2n )] 2

− g · t · J0n

(5.45)

λ c2 n n + Vn−2 ) · Vnn · (3hnn − 4hnn−1 + hnn−2 ) + · (3Vnn − 4Vn−1 2 g λ2 2 · V0n · c2 + · ((Vnn )2 + c2 ) · (hnn − 2hnn−1 + hnn−2 ) + 2 g n n · (Vnn − 2Vn−1 + Vn−2 ) (5.46)

hn+1 = hnn − n

Vnn+1 = Vnn − +

λ n n · [g · (3hnn − 4hnn−1 + hnn−2 ) + Vnn · (3Vnn − 4Vn−1 + Vn−2 )] 2

λ2 · [2 · g · Vnn · (hnn − 2hnn−1 + hnn−2 ) + ((Vnn )2 + c2 ) 2 n n + Vn−2 )] − g · t · Jnn · (Vnn − 2Vn−1

(5.47)

Unfortunately, this method, probably because of the presence of the second-order derivatives, seems to be very unstable for the case of waterhammer and, therefore, will not be implemented in the computer code.

5.4.2

Ghost cells

The idea of the method [36] is the addition of a further cell beyond the physical domain of the problem, i.e., the numerical domain is defined for the points in the set [−1 . . . n + 1], while in the real domain the points belong to the set [0 . . . n]. The values of the parameters in the ghost cells have to be defined accordingly, and a number of methods exists for this purpose. A widely used method is based on internal solutions, i.e., the definition of a cell n n Qj+1 through the values computed inside the domain, as Qjn or Qj−1 . The simplest approach is defined zero-order extrapolation which defines the external parameters as: n Qj+1 = Qjn

(5.48)

The insertion of this condition inside the computer program is very easy, as the Lax–Wendroff equations can be used and simplified. For the upstream condition

66 WATER HAMMER SIMULATIONS these become: V0n+1 = V0n −

λ · (g · (hn1 − hn0 ) + Vjn · (V1n − V0n )) 2

λ2 · ((2 · V0n · g) · (hn1 − hn0 ) + (c2 + (V0n )2 ) · (V1n − V0n )) − g · t · Jjn 2 (5.49) 2 λ c = hn0 − · V0n · hn1 − hn0 + · V1n − V0n 2 g λ2 2 · V0n · c2 n 2 2 n n n n + · ((V0 ) + c ) · (h1 − h0 ) + (5.50) · (V1 − V0 ) 2 g +

hn+1 0

while for downstream conditions the equations become: λ n )) · (g · (hnn − hnn−1 ) + Vnn · (Vnn − Vn−1 2 2 λ2 n + )) · ((2 · Vnn · g) · (−hnn + hnn−1 ) + (c2 + Vnn ) · (−Vnn + Vn−1 2 − g · t · Jjn (5.51)

Vnn+1 = Vnn −

λ c2 n · (Vnn · (hnn − hnn−1 ) · (Vnn − Vn−1 )) 2 g λ2 2 · Vnn · c2 n 2 2 n n n n + · ((Vn ) + c ) · (−hn + hn−1 ) + · (−Vn + Vn−1 ) 2 g

hn+1 = hnn − n

(5.52) It is obviously possible the use of the so-called first-order extrapolation based on the filtering of the internal solutions through a linear function: n n Qj+1 = 2 · Qjn − Qj−1

(5.53)

Many schemes are available, and the “best” scheme, which can be successfully used for any conditions unfortunately does not exist. However, the methods presented in this chapter are normally valid for most real cases.

5.5 The computer code The code for using the Lax–Wendroff scheme is as simple as that shown in Chapter 4 to implement the characteristic method. The core procedure, to solve internal points (which is the actual Lax–Wendroff scheme), is shown in Figure 5.8. Again, it can be seen that this procedure uses the exact equations (5.28) and (5.29). In order to compare the carried out results with those obtained in Chapter 4, the same simple plant is simulated: with an upstream reservoir and a downstream valve which closes with assigned law. The UpWind and the Ghost Cell Zero Order to solve the boundary conditions are implemented, and the related codes are shown in Figures 5.9 and 5.10.


67

procedure ComputeInternalPoints; var i:integer; begin for i:=1 to n-1 do begin vp[i]:=v[i]-0.5*(dt/ds)*((h[i+1]-h[i-1])*g+v[i]*(v[i+1]-v[i-1]))+ 0.5*(dt/ds)*(dt/ds)*(2*g*v[i]*(h[i+1]-2*h[i]+h[i-1])+ (cel*cel+v[i]*v[i])*(v[i+1]-2*v[i]+v[i-1]))Lambda*v[i]*abs(v[i])*dt/(2*diam); hp[i]:=h[i]-0.5*(dt/ds)*(v[i]*(h[i+1]-h[i-1])+(Cel*Cel/g)*(v[i+1]-v[i-1]))+ 0.5*(dt/ds)*(dt/ds)*((v[i]*v[i]+Cel*Cel)*(h[i+1]-2*h[i]+h[i-1])+ (2*Cel*Cel*v[i]/g)*(v[i+1]-2*v[i]+v[i-1])); end; end; Figure 5.8: Core procedure of the solution program: calculation of the unknown (head and velocity) for the points internal to the pipe – Lax–Wendroff scheme. procedure UpWind; begin (* Reservoir positioned upstream *) hp[0]:=h0; vp[0]:=v[0]-0.5*(dt/ds)*(g*(h[1]-h[0])+v[0]*(v[1]-v[0]))Lambda*v[0]*abs(v[0])*dt/(2*diam); (* Valve positioned downstream *) if t
68 WATER HAMMER SIMULATIONS procedure GhostZero; begin (* Reservoir positioned upstream *) hp[0]:=h0; vp[0]:=v[0]-0.5*(dt/ds)*((h[1]-h[0])*g+v[0]*(v[1]-v[0]))+ 0.5*(dt/ds)*(dt/ds)*(2*g*v[0]*(h[1]-h[0])+ (cel*cel+v[0]*v[0])*(v[1]-v[0]))-Lambda*v[0]*abs(v[0])*dt/(2*diam); (* Valve positioned downstream *) if t
5.6 Again the simple application The simulated plant is exactly that described in section 4.4 and therefore its characteristics will not be repeated here. The user’s interface is very similar to the previous program, with slight differences; it is reported in Figure 5.11. As mentioned, in this case both the spatial interval and the time step can be decided by the user, and the computer program computes the CFL condition. As can be seen (and the reader is invited to try), when the CFL number is higher than 1, the code crashes and it is impossible to see the results. However, some interesting comparisons can be performed for CFL numbers smaller than 1. In Figure 5.12, two simulations are presented: in both cases the time step is very small (t = 0.0008 s) but in one case the number of computational

Figure 5.11: User’s interface of the program. 200 180

CFL 1.0 CFL 0.1

160

Head (m)

140 120 100 80 60 40 20 0 0.0

1.5

3.0

4.5

6.0

7.5

9.0

10.5

12.0

13.5

15.0

Time (s)

Figure 5.12: Results of the simulation of a sudden closure and negligible headlosses, with Lax–Wendroff scheme, UpWind scheme for boundaries and different CFL numbers.

70 WATER HAMMER SIMULATIONS section is let equal to 1000 (and this implies CFL = 1.0) while in the latter case, in order to increase the velocity of computations the section number is equal to 100, and therefore CFL = 0.1. As can be seen, even if the code does not crash, when CFL = 0.1 noticeable instabilities appear and the carried out solution can be said to be less accurate, while with CFL = 1.0, the solution is identical to that of the theoretical. The shown simulations have been performed with the UpWind scheme for the boundary conditions, but results carried out with the Ghost Cells are very similar.

Chapter 6 Devices – Boundary conditions

In a complex plant many devices can be inserted. In this chapter few of them are presented, and in particular those purposely studied for reducing the pressures generated during transients, because, as it has been shown in the past chapters, these pressures can be very high and dangerous. Moreover, a section discussing the possibility to use a simple model for pumps is also presented. As waterhammer is known since many decades (or centuries), a large number of devices have been studied and implemented, and the research is still going on. The devices presented in this chapter are the most widespread, and their discussion and implementation in the computer programs attached to the book should allow the reader to generalize the method and to apply it to any technical case. These devices are modelled as boundary conditions and positioned upstream or downstream the pipes; pipes are, therefore, solved with the methods developed for internal points (Chapters 4 and 5) and boundary conditions are solved in a different computer procedure, as described in the following.

6.1 Surge tanks 6.1.1 Simple surge tanks Hydroelectric plants are normally composed by an upstream reservoir with the intake, a tunnel where water flows with relatively low velocity (to minimize the headlosses) and a penstock, where water flows at high velocity and the differences in elevations are maximized, in order to retrieve the maximum power (Figure 6.1). Between the tunnel and the penstock is normally positioned a surge tank, i.e., a reservoir with the aim to receive the inflowing water when the valve situated downstream the penstock is closed and the turbines are turned off. This has already been mentioned in Chapter 3. When the flow in the penstock is abruptly stopped, the water flowing in the tunnel does not stop immediately but, due to its inertia, flows into the surge tank, increasing its level up to a maximum, that occurs when all the kinetic energy of the water in the tunnel has been transformed, except for the headlosses, into potential energy. After that, the flow reverses and the level in the surge tank decreases. The process continues with a series of oscillations, damped by the headlosses of the system.

72 WATER HAMMER SIMULATIONS 45 m

45 m A Upper expansion chamber

A

2,5

0

Cross section AA

Cross section CC

1,75

0 2,5

3,44

29 m

B

∅ 100

Cross section BB

B

4

3.4

2,10

2,10

Lower expansion chanber

C

D C

Tunnel

Cross section DD

D

70 m

Figure 6.1: Hydroelectric plant – the Plima-Lasa plant, North-East of Italy. The model implemented in a computer code using the method of characteristics is very simple, and it uses the following equations: 1. 2. 3. 4. 5. 6.

positive characteristic line (inside the tunnel); negative characteristic line (inside the penstock); continuity equation for the three pipes: tunnel, surge tank and penstock; equality of the heads for tunnel and surge tank; equality of the heads for surge tank and penstock; equality of the head at the base of the surge tank and the water elevation in the tank; 7. variation of the elevation z within the surge tank, computed as mass conservation, evaluating the inflow discharge Qtank = Atank · z/t where Atank is obviously the cross section of the surge tank. Equations are, respectively: htunnel,j − htunnel−1,j−1 c Vtunnel,j − Vtunnel−1,j−1 + · t g t Vtunnel−1,j−1 · |Vtunnel−1,j−1 | +c · λ · 2·g·D

(6.1)

73

DEVICES – BOUNDARY CONDITIONS

Q surge tank

Downstream tunnel section Upstream penstock section Surge tank base section

Tunnel section 1 Penstock section 1 Q tunnel

Q penstock

Figure 6.2: Sketch of the notation used for developing the computer code.

hpenstock,j − hpenstock+1,j−1 c Vpenstock,j − Vpenstock+1,j−1 − · t g t Vpenstock+1,j−1 · |Vpenstock+1,j−1 | −c · λ · 2·g·D

(6.2)

Qtunnel = Qtank + Qpenstock

(6.3)

htunnel = htank

(6.4)

htunnel = hpenstock

(6.5)

Ztank = htank Qtank Z = · t Atank

(6.6) (6.7)

In the above reported description the subscripts tunnel, penstock and tank have been used to indicate the same point positioned at the junction with the three elements but related to each of the element (Figure 6.2). With the subscripts tunnel − 1 and penstock + 1 the sections in the tunnel immediately upstream and in the penstock immediately downstream the surge tank are designated, respectively. The above equations act as downstream boundary condition for the tunnel and as upstream boundary condition for the penstock. These equations, together with those for the internal points and the initial conditions, allow the solution to the problem.

74 WATER HAMMER SIMULATIONS It is to be observed that no characteristic equation is written for the surge tank, for which eqn (6.6) is written instead: this implies the static pressure distribution and the non-elastic behaviour of the water in the surge tank. As can be seen in the example reported in Figure 6.1 and in other examples that can be found in the literature, these devices are very large and, therefore, are very expensive: any reduction of these dimensions is appreciated. To this end, it is possible to shape the base of the surge tank in order to generate a further (concentrated) headloss, which is active only during transients and it has the purpose to reduce the oscillation amplitudes (restricted orifice). The effect of this orifice can be modelled with a slight change in eqn (6.4), with the more general: htank = htunnel + k · Vtank · |Vtank |

(6.4 )

The coefficient k is chosen to represent the effect of the restricted orifice, noting that during steady flow (Vtank = 0) and when the orifice is not present (k = 0) htank = htunnel [37]. At the beginning of transients, the headlosses due to the orifice produce an high value of pressure at the base of the surge tank, which fosters the deceleration of the velocity in the tunnel and, therefore, damps the oscillations in the tank. If the value of k is very low, the effect of the restricted orifice is negligible; on the other hand, if it is too high, the pressures at the beginning of the transient are very high: this reduces the oscillations inside the surge tank and also its effect in damping the pressure peaks in the tunnel. It is, therefore, possible to devise an optimal dimension of the restricted orifice, so that the pressures reach immediately their maximum value, and they stabilize at that level for the whole duration of the transient. In other words, the tunnel is subjected for the whole transient duration to the maximum value of pressure, which is slightly smaller than the pressure it would reach without this added headloss and, therefore, it reduces the oscillations inside the surge tank. To design a simple surge tank, the Thoma’s stability condition [38–40] has to be satisfied, that, for small oscillations, the minimum cross sectional area of the surge tank is: Atunnel · Ltunnel Amin tank = (6.8) 2 2 · g · (Htunnel /Vtunnel ) · Hnet where Htunnel is the headloss in the tunnel, Hnet is the difference in elevation (so neglecting the headlosses) between the water elevation in the surge tank in static conditions and the elevation of the downstream turbine. Normally, this minimum value is multiplied for a safety factor that can vary between 1.2 and 1.5. 6.1.2

Different types of surge tanks

As mentioned above, dimensions of surge tanks are very large. Moreover, in the practical engineering other targets are required, as the reduction of the oscillation period and the duration of the transient; furthermore, if a surge tank is connected to


75

Riser Upper expansion chamber

Lower expansion chamber

Tunnel

Penstock

Figure 6.3: Sketch of a differential surge tank. a reservoir which water level is not constant but shows wide variations, this surge tank would be oversized for most of the period of use when the water level in the reservoir is low. To overcome these drawbacks various types of surge tanks have been studied [41, 42]. In this section the basic types are described, but the creativity of designers has suggested several others, according to the needs of the specific projects. One of the widespread types is built with expansion chambers, as sketched in Figure 6.3. In general, the higher expansion is kept above the highest level of water elevation in the reservoir, while the lower expansion is built below the level of water obtained with maximum discharge and lowest water elevation inside the upstream reservoir. Under these conditions, during steady flow, the water surface is always contained within the riser of the surge tank. If, during transients, water is contained in the riser, the surge tank behaves as in the simple case (Section 6.1.1). When, instead, the expansion chambers are also affected by the motion, the level in the tank rises (or lower) quickly until it is contained in the riser, to rise (or lower) much more slowly when water reaches one of the expansions. With these surge tanks, then, an high (or low) pressure quickly establish at the base of the riser while the expansions limit the amplitude of the oscillations. This system, surely advantageous, may introduce complications to the analysis, as waves can be produced in the galleries during emptying and filling [43]. In order to slightly increase the pressures at the base of the riser, sometime it is continued beyond the bottom of the chamber, an overflow is created, as sketched in Figure 6.4. A similar device has been studied by Kammuller [44] to take advantage of the effect of rapid lowering of the level in the riser; the device is sketched in Figure 6.5.

76 WATER HAMMER SIMULATIONS Upper expansion chamber

Riser

Overflow

Figure 6.4: Overflow in the upper expansion chamber of the surge tank.

Riser

a Lower expansion chamber

c

o Tunnel

o Penstock

Figure 6.5: Lower expansion chamber with Kammuler’s device. Because of this device, the lower expansion chamber starts to empty only when the water level is below the section 0–0, i.e., when from the conduit c the air can enter the chamber. The pipe a, much smaller than pipe c, serves only to expel the air when the chamber fills. A different idea to reach the same target is the use of a differential surge tank, or Johnson type [45, 46], equipped with an optimal orifice for the inflowing water, but also with an internal riser that overflows in the outer tank when the maximum water level is reached. This type of surge tank is sketched in Figure 6.6. When the valve positioned downstream the penstock opens, and, therefore, water starts flowing in the penstock, the decrease of the pressure under the orifice is equal to the maximum lowering of the water level in the riser, as the orifice is optimal for this operation.


77

Orifice

Q1

Q2 Q Tunnel

Q Penstock

Figure 6.6: Differential surge tank, or Johnson’s type. When the valve closes, water stops in the penstock and flows into the surge tank; the orifice is small and, without the internal riser, pressures at the base of the surge tank would be too high. The velocity inside the riser (related to discharge Q1 ) is high and the water surface reaches the maximum level when the level in the outside tank (given by discharge Q2 ) is still very low; when the water level reaches the top of the riser, it overflows in the outside tank thus stabilizing the value of the pressure. The surge tank is, therefore, optimized for both the operations of opening and closing of the plant.

6.2 Air chambers These devices are the most common in pumping stations, where an abrupt stoppage of the flow induces in the downstream pipe low pressure values that can have significant impact on the system. These very low pressures, in the worst cases, produce cavitation, an undesirable phenomenon that will be described later in the book. Air chambers are closed reservoirs, positioned downstream the pump, filled by both water and pressurized air [47, 48]. When the pump stops, due to this device the flow in the downstream pipe continues, fed by the air chamber, and its velocity slowly decreases, therefore limiting the waterhammer effects. During the transient, the volume of air increases and its pressure decreases until it reaches a minimum; afterwards the flow inside the pipe reverses, and in this second phase, the air volume decreases and the pressure increases up to its maximum value. The process progresses with a series of oscillations, damped by headlosses.

78 WATER HAMMER SIMULATIONS Air chamber

Section “0”

Valve

Figure 6.7: Sketch of an air chamber. Variables that govern the boundary conditions are the head in the air chamber h (or the pressure p), the velocity at the inlet of the chamber Vs and the volume of air in the chamber itself W . Equations that solve the system (in addition to the usual characteristic equation) are the mass conservation for the system air chamber-pipe and the polytropic law that describes the behaviour of the air inside the chamber. Considering the air chamber positioned downstream, the pumps or the valve which suddently stops the incoming discharge, and therefore upstream the subsequent pipe, i.e., in the section 0 of the computational domain (see Figure 6.7) the negative characteristic equation is: h0,t+1 − h1,t −

c V1,t · |V1,t | · (V0,t+1 − V1,t ) − c · λ · · t = 0 g 2·g·D

(6.9)

The continuity equation for the air inside the chamber is written as: Wt+1 = Wt + A · Vo,t · t

(6.10)

The polytropic equation that describe the behaviour of air is: m = c.te ho,t · W0,t

(6.11)

where m depends on the polytropic process and can vary between 1.0 and 1.4 for the isothermal and adiabatic (considering dry air at 20◦ C) processes, respectively [49]; this parameter is normally let equal to 1.4, the differences in the carried out results being very limited and this value being on the safe side. Without the additional headloss, as those described for surge tanks in Section 6.1, the assumption is that the head in the first section of the downstream pipe is equal to that inside the chamber. Actually, even if it is possible to design orifices for these devices as well as for the surge tank, they are not normally used, as the economic saving is negligible and an error in design can worse the normal


79

conditions. However, as will be seen in Section 6.4.3, the implementation of an headloss model is very easy. In the thirties of the twentieth century, Evangelisti [50, 51] devised a graphical method to dimension air chambers, but only under the hypotheses of no headlosses and non-elastic behaviour of the system. With numerical models as those presented in this book, the dimensioning of these devices is carried out by trial and errors, selecting a volume for the air chamber and numerically testing its effect in the system. With regard to the installation of air chambers, it is to be observed that they are normally equipped with a compressor in order to re-establish the mass of air within the chamber, as, because of the high pressures, air tends to dissolve in water and then to be convoyed away. For some installations there is the possibility of having the aeriform contained in a pouch which isolates it from the liquid, so that no dispersions occurs and, therefore, no compressor is required. However, the possibility to recharge that pouch should always be granted; a possible drawback of this system is given by the life expectancy of the pouch.

6.3 Relief valves and rupture disks Let us consider again the simple plant of Figure 2.1. As already mentioned, when the downstream valve closes, the pressure increases up to a value than can be computed with eqn (2.1) or (2.9) if the closure is abrupt or slow, respectively. Relief valves work with the principle sketched in Figure 6.8. The valve is closed by a spring, which is calibrated to keep the valve closed under static and steady flow conditions, with a tolerance of approximately 5%. When, during a transient, the pressure increases and exceeds the calibrated value of the spring, the valve opens and water outflows. The discharge outflowing from the valve can be computed with: Q = µ(x) · A ·

"

2·g·H

(6.12)

where x is the distance between the disk and the outflow of the valve. Obviously, if that distance is zero, the disk prevents the discharge to outflow and, therefore, µ(x) = 0; for small distances, Cozzo [52] experimentally carried out the values of the efflux coefficient µ. This coefficient is affected from the disk because of the turbulence and its value cannot be determined theoretically. The use of these typology of valves requires attention because of the possibility of their resonant behaviour. Tanda and Zampaglione [53] presented a real case where the installation of these valves rose instabilities in the system: they oscillated with a frequency close to the natural frequency of the plant, inducing resonance. In other words, not only the pressures did not decrease, but they actually increased, with the risk of driving the system to the collapse. The plant the authors studied is sketched in Figure 6.9: at the end of each of the four branches, a regulation valve is positioned; immediately upstream of the latter, in order to reduce the pressures generated during


Relief valve

Valve

Figure 6.8: Sketch of a relief valve.

Regulation valves

5

6

7

1

12

8

2

13

9

3

14

10

4

15

11

Relief valves

Figure 6.9: Plant studied by Tanda and Zapaglione [53] which showed instabilities when the anti-waterhammer valves were installed. transients, not only one, but two relief valves were positioned. This decision was taken because of an error in the dimensioning of the relief valves, which were not sufficient to reduce the pressures as it was desired. However, the positioning of two valves for each branch was probably the cause of the generated instabilities. Some authors [19] point out that the discharge capacity of the valve should not be too large, because if that discharge is insufficient to keep the valve fully open, then hammering and oscillation may occur. Manufacturers develop valves that use damping mechanisms which allow a quick opening, as requested in order to control the phenomenon, but a slow closure, to eliminate any risk of resonance. Bianchi and Mambretti [54] showed that the use of different valve typology which adopts the same principle (they open when a given pressure is exceeded) may solve the problem. The valve studied in the mentioned paper (which will be presented in Section 6.5.5), in fact, opens quite quickly, being its velocity similar to the relief valves; the difference with the latter is that once it is completely open, it remains in this position until the transient has ended. The use of rupture disks is now quite widespread, which work using the same principle: they open (break) when the threshold pressure is reached, to allow the


81

water outflow. Also in this case instabilities are not possible, as obviously the disk does not close again at the end of the transient; it is also obvious that the installation of these devices should be designed only for emergencies, as they require their manual substitution and that means the plant has to be stopped, at least partially. Moreover, when the rupture threshold is close to the working pressure, fatigue has been detected and after some time from the installation the disks can leak.

6.4 Centrifugal pumps Pumps have always been part of complex plants and therefore a large number of models have been implemented to describe their behaviour. In this book only centrifugal pumps are covered, as they are the most commonly installed pumps; they produce transients only when they change the pumped discharge. Other pumps, as reciprocating pumps or even the hydraulic ram (see Chapter 1) can produce waterhammer as by-product of their operations. The first and simplest idea was that the pump inertia is so low that the time of stopping of the pumps is negligible, and therefore no models are needed. For instance, all the charts prepared in the first mid of twentieth century to design air chambers adopted this hypothesis, which is surely on the safe side, but that may be too demanding. Therefore, during the years many methods for pump modelling have been presented in the scientific literature. In this paragraph, according with the aim of the book, not all these models are covered, but only a simple model, quite widespread, is shown, discussed and implemented in a computer code. When a pump is stopped, the rotor slows down with a trend that can be computed equalling the kinetic energy variation of the rotor and the power of the hydraulic flow [10, 19]: −

d I · ω2 γ · Q · H = dt 2 η

(6.13)

where I is the rotor inertia, that can be computed by: I=

G · D2 4·g

(6.14)

G being the weight of the rotor and D their inertial diameter. The value of G · D2 is normally provided by the companies producing the pumps, but there have been studies about its computation [55]. In (6.13), ω is the angular velocity of the pump, η is the efficiency and γ is the specific weight of the fluid. As usual, Q is the discharge and H is the head of the pump. Equation (6.13) and the characteristic equation of the system allow the solution of the problem.

82 WATER HAMMER SIMULATIONS Actually, eqn (6.13) can be generalized to model the starting phase of the pump [56]: −

γ · Q · H d I · ω2 = −P dt 2 η

(6.13 )

where P is the power provided to the pump. When the pump is stopping, P = 0 and eqn (6.13 ) is obviously equal to (6.13). In steady flow, γ · Q · H =P η that yields to d I · ω2 =0 dt 2 that means the angular velocity is constant. The problem in modelling the starting phase of the pump is the determination of the variation of the power P when the pump starts, which is normally not known and that makes difficult the construction of a reliable model. However, differently from a pump trip that cannot be controlled, the start-up of the pump can always be planned, and it should be to avoid undesiered surges. It is quite common nowadays to provide the pumps with electrical devices able to variate their speed in steady flow conditions and that are also able to “soft” start them at the beginning of the operations. Cheaper devices are able to control only the start up of the pumps (the so-called “soft-starter”) but not to change their velocity in normal operations. In any case, the manufacturers are able to produce a large variety of methods and devices to “ramp up” the pumps and the value of P in eqn (6.13 ) should be modelled accordingly. A further difficulty is given by the incomplete knowledge of the characteristic curves of the pumps when they are not in steady flow conditions. Therefore, when these curves are not provided by the producers of the pumps, the modeller is forced to accept the similitude hypotheses [57]: n1 Q1 = Q2 n2 2 H1 n1 = H2 n2

(6.15) (6.16)

where n is the rotations per minute of the pump, tied to the angular velocity ω by the: 2·π ·n (6.17) 60 Equations (6.15) and (6.16) are acceptable only when the velocity of the pump is not too far from that of the steady flow, and therefore their use is suggested only when no better information are available. ω=


83

From (6.13 ), as the inertia I is constant, the variation of the rotation with time dn/dt is equal to: −

I · π2 dn γ · Q · H ·n· = −P 900 dt η

(6.18)

and, therefore, the pump rotation at time t + 1 can be computed as: γ · Q · H 900 nt+1 = nt + P − · dt · η nt · I · π 2

(6.19)

The curve of the pump is normally given as second-order polynomial, i.e., in the form: h = a · V2 + b · V + c

(6.20)

where h is the head of the fluid and V its velocity. The parameters a, b and c characterize the behaviour of the pump. The problem is how to scale these parameters to produce new curves to be used in unsteady flow in accordance with the similitude laws (6.15) and (6.16). The following assumptions are needed to produce the required rules [56, 58]. Assumption 1. The curve is scaled so that when V = 0 it results ht+1 = ct+1 = nt+1 2 ht · ; as a consequence: nt ct+1 = ct ·

nt+1 nt

2 (6.21)

Assumption 2. The derivatives of the curves must be equal when V = 0; as a consequence: nt+1 bt+1 = bt · (6.22) nt Assumption 3. When vt+1 = vt ·

nt+1 nt

it must be ht+1 = 0, and therefore:

at+1 = a = c · te

(6.23)

An example of the curve scaling is shown in Figure 6.10. This model can be simply inserted in a computer code that solves the internal points with one of the methods presented in Chapters 4 and 5; the pump is considered as a boundary condition that can be positioned upstream or downstream, even if the former is the normal position for a pump.

84 WATER HAMMER SIMULATIONS n 100% n 80% n 50%

1.0

H/Hmax

0.8

0.6

0.4

0.2

0.0 0.0

0.2

0.4

0.6 V/Vmax

0.8

1.0

Figure 6.10: Curves of the pump for different rotations. Assuming the pump is positioned upstream, the numerical solution is carried out coupling the negative characteristics and the equation of the pump: ht+1,0 − ht,1 −

c Vt,1 · |Vt,1 | · (Vt+1,0 − Vt,1 ) − c · λ · · t = 0 g 2·g·D

(6.24)

at+1 · Vt+1,0 · |Vt+1,0 | + bt+1 · Vt+1,0 + ct+1 = ht+1,0

(6.25)

where at+1 , bt+1 and ct+1 are computed by eqns (6.21)–(6.23). Equation (6.19) shows that a larger pump inertia implies a slower stoppage of the pump itself. This suggests the idea to artificially increase the pump inertia in order to lengthen the time of pump stopping and therefore reducing (according to Michaud’s formula) the pressures related to transients [59]. This target is normally reached designing “flywheels”, i.e., metal disks rotating with the mobile part of the pump, at the same velocity. This solution is surely feasible, especially when the pumps are continuously working. Limitations are given by the capability of the motor and of the equipments and, for large installations, by the increased loads and vibrations. Moreover, when pumps are often turned on and off, the drive shaft might result over-stressed by the increased mass of the rotor. However, the actual effectiveness of this solution will be discussed in Section 6.5.6.

6.5 Other methods for controlling the pressures The reader can find in the literature more methods to try to reduce the pressures generated during the transients, but, as the parameters that govern the phenomenon are the difference in the initial and final velocities and the celerity of the wave, these are also the ways in which actions can be taken.


85

To reduce the wave celerity, it is possible to reduce the density of the fluid or to increase the deformability of the pipe. The reduction of the density of the fluid can be obtained bleeding air in the liquid: but probably this method produces more disadvantages than advantages and it is better to avoid it. Ramenieras [60] developed a device which consisted a small flexible hose in the pipe where the air is trapped, so that the effective bulk elasticity modulus of the system is reducd. Ghilardi and Paoletti [61] studied the possibility to add a flexible viscoelastic pipe to be used as pressure surge suppressor; this pipe deforms during transients so that the celerity and the consequent pressures decrease. The exact behaviour of this added pipe is still under investigation [62, 63]. Moreover, it is still to be investigated how the repeated pressure fluctuations influence the life expectancies of these pipes.

6.6 The computer codes Again the computer codes are based on the method of characteristics, and only the boundary conditions are changed in order to model the devices presented in the chapter. The routine for solving internal points is presented in Chapter 4 together with those for solving the boundary conditions related to the reservoir and the valve. 6.6.1 The simple surge tank Boundary conditions for the simple surge tanks are those described in Section 6.1.1 and in the computer program these are reported as shown in Figure 6.11. As can be seen with very simple calculations, the equations have been solved algebraically SA:=f1*cel1*v[n1-1]*abs(v[n1-1])/(2*g*d1)*dt; SB:=f2*cel2*v[n1+1]*abs(v[n1+1])/(2*g*d2)*dt; SC:=h[n1-1]+cel1*TankArea*ZTank/(g*Area1*dt)+ +cel1*Area2*h[n1+1]/(Area1*cel2)+cel1*Area2*SB/(Area1*cel2)+ -cel1*Area2*V[n1+1]/(g*Area1)+cel1*V[n1-1]/g-SA; SD:=1+cel1*TankArea/(g*Area1*dt)+Area2/Area1; ZpTank:=SC/SD; hpd:=ZpTank; hps:=ZpTank; VpTank:=(ZpTank-ZTank)/dt; vpd:=(ZpTank-h[n1+1]-SB)*g/cel2+V[n1+1]; vps:=((ZpTank-ZTank)*TankArea/dt+Vpd*Area2)/Area1; Figure 6.11: Boundary conditions for the solution of a simple surge tank.

86 WATER HAMMER SIMULATIONS and then implemented. Four service variables (SA, SB, SC and SD) have been used to ease the formulas. A simple example can be studied, still referring to a plant with an upstream reservoir, a valve downstream and a simple surge tank to divide the tunnel and the penstock. The characteristics of the plant are as follows. The discharge in steady flow conditions is Q0 = 15 m3 /s and the head of the upstream reservoir is 100 m. The tunnel has a diameter of Dtunnel = 2.80 m, length Ltunnel = 4200 m, Manning’s roughness n = 0.0143 s m−1/3 , celerity ctunnel = 1000 m/s. The characteristics of the penstock are: diameter equal to Dpenstock = 1.90 m, length Lpenstock = 1400 m, Manning’s roughness n = 0.0111 s m−1/3 , celerity cpenstock = 1000 m/s. The surge tank has no orifice and it is simple in type; its diameter is Dtank = 11.50 m and therefore the area of the surge tank is Atank ∼ = 104 m2 . The downstream valve closes instantaneously. The steady flow water velocity in the penstock is equal to 5.29 m/s, the phase time is 2.8 s and therefore to compute the maximum pressure, the Michaud’s formula (Chapter 2) has to be used, which yields: h =

2 · L · V0 2 × 1400 × 5.29 = 302 m = g · Tc 9.806 × 5

Results from the computer programs are saved in the files heads.txt and velocities.txt in a simple text format that can be imported in any spreadsheet. The plot of the heads computed at the downstream section, where the valve is positioned, are shown in Figure 6.12, and a very good agreement can be observed with the Michaud’s result. The oscillations in the surge tank are reported in Figure 6.12. Results are reported in two different figures because of the large differences between the two phenomena, both in terms of pressure values and in terms of time. As can be seen, the high pressures generated inside the penstock are strongly damped by the surge tank, and in the tunnel only the pressures given by the oscillations within the surge tank propagate. The results presented in Figure 6.13 are those already shown in Chapter 3, when this solution, based on elastic hypothesis, was compared with the solution carried out with the simplified hypothesis of anelastic fluid. 6.6.2 The simple air chamber Boundary conditions reported in Section 6.2 are quite simple, but in this section a number of models are implemented, with increasing complexity. The simplest model is mentioned in section 6.2 and its implementation is shown in Figure 6.14; to call this routine in the computer program attached to the book the option “Simple Air Chamber” has to be selected. With the same program is also possible to simulate

DEVICES – BOUNDARY CONDITIONS 400 350

Head (m)

300 250 200 150 100 50 0 0

10

20

30 Time (s)

40

50

60

500

600

Figure 6.12: Results at the downstream valve. 110 108

Water elevation (m)

106 104 102 100 98 96 94 92 90 0

100

200

300 Time (s)

400

Figure 6.13: Oscillations of the water in the surge tank. Wtp:=Wt+Area1*V[0]*dt; hp[0]:=h[0]*power(Wt/Wtp,m); Vp[0]:=V[1]+g*(hp[0]-h[1])+V[1]*abs(V[1])*dt/(2*d1); Figure 6.14: Boundary conditions for the solution of an air chamber.

87


Head (m)

120 100 80 60 40 20 0 0

15

30

45

60

75

90

105

120

135

150

Time (s)

Figure 6.15: Results of an instantaneous valve closure in a plant without the air chamber. the closure of an upstream valve without air chamber, and therefore comparisons are possible. For instance, let us consider the case of a pipe with diameter of 0.3 m, length 1000 m, Manning’s roughness coefficient 0.0125 s m−1/3 and celerity 1000 m/s. The steady flow discharge is equal to 0.050 m3 /s and downstream is positioned a reservoir and the fixed head is equal to 100 m. Upstream there is a valve which closes completely and instantaneously. The results computed in the upstream section without the air chamber (Figure 6.15) and with an air chamber with 1 m3 volume (Figure 6.16) are shown, with the same scale on the xy axes to allow a better comparison. As can be seen the air chamber has the double effect of reducing the pressures generated during the transient and of increasing the time of the oscillations. 6.6.3 Air chamber with headlosses As mentioned above, this very simple boundary condition can be complicated by inserting an headloss in the passage from the chamber to the conduit. In Section 6.2 it has also been mentioned that, in the real world, it is not usual to insert a valve in the air chambers which has the same aim of an orifice in the surge tank; however, it has been implemented and the results are commented, as it is a quite simple improvement, it is useful for understanding the behaviour of such a device and allows an interesting discussion.


89

180

160 140

Head (m)

120 100 80 60 40 20 0 0

15

30

45

60

75 Time (s)

90

105

120

135

150

Figure 6.16: Results of an instantaneous valve closure in a plant with an air chamber of 1 m3 volume.

Wtp:=Wt+Area1*V[0]*dt; hpChamber:=hChamber*power(Wt/Wtp,m); hp[0]:=hpChamber-KOrifice*V[0]*abs(V[0]); Vp[0]:=V[1]+g*(hp[0]-h[1])/cel1-f1*V[1]*abs(V[1])*dt/(2*d1); Figure 6.17: Boundary conditions for the solution of an air chamber with headlosses (orifice or valve). According to the description of an orifice, the difference between the heads inside the chamber (hChamber ) and in the pipe (hPipe ) are given by the headlosses, so that: hChamber = hPipe + k · V · |V | The constant k has to be appropriately selected and depends on the geometry of the orifice (or valve). The code is shown in Figure 6.17 and the results for different values of the parameter k are shown in Figure 6.18, for the section where the valve is positioned (section 0). As it can be seen (Figure 6.18), the increase of the parameter k produces a sudden diminution of the pressure in the section 0, followed by the usual oscillations which have a reduced amplitude.

90 WATER HAMMER SIMULATIONS 120 K=0 K=5 K = 10 K = 20

115

Head (m)

110 105 100 95 90 85 80 0

5

10

15 Time (s)

20

25

30

Figure 6.18: Results of an instantaneous valve closure in a plant with an air chamber of 1 m3 volume for different values of the parameter K for the headlosses. However, the very simple code here presented shows to be unstable, as can be seen when looking at the data related to k = 20. This result is very interesting and deserves a deeper analysis. At the beginning of the transient, because of the concentrated headloss, a low value of pressure at the base of the air chamber is generated; this consequently reduces the flow velocity and, as the headlosses depend on this latter parameter, the pressure at the base of the air chamber increases again. This rapid variation of the pressure values propagates downstream, producing instabilities. The method of characteristics, in fact, is very stable by itself, but it is strongly affected by instabilities coming from the boundary conditions and this is precisely the case. Actually, the instability is surely numerical, as it is not known to the author any literature case reporting instabilities in the use of air chambers. However, the reduction of the time step do not produce any improvements. The solution, in this case, would be the change of the integration method and the insertion of limiters, as will be discussed in Chapter 7, where the problem will be taken up again and solved. 6.6.4 Air chamber and valve Finally, it is possible to study the case of a valve that closes in finite time and it is positioned together with an air chamber. Equations in this case are the following: 1. Negative characteristic equation (6.9); 2. Head in the pipe at section 0 is equal to head at the valve;


91

if t
3. 4. 5. 6.

Head in the pipe at section 0 is equal to head in the chamber; Continuity equation for the node; Continuity equation for the chamber (6.10); Polytropic equation (6.11).

As can be seen with very simple calculations, the equations can be solved algebraically and they have been implemented in the computer code. Four service variables (SA, SB, SC and SD) have been used to ease the formulas. The computer program is shown in Figure 6.19; the results carried out for different time of closure of the valve are reported in Figure 6.20.

6.6.5 Valve modelling: an example The description of all the valves available on the market is obviously impossible, thus in the following only one analysis is reported in order to show how this problem may be faced. The complete analysis is reported in [54].


100

Head (m)

80

60

40

20

Tc 0 s Tc 10 s

0 0

15

30

45

60

75 Time (s)

90

105

120

135

150

Figure 6.20: Effects of time closures of a valve in a plant with an air chamber with volume equal to 1 m3 .

The valve is sketched in Figure 6.21 and it is set to work as relief valve. It has two compartments: the lower (I) is positioned in direct contact with the pipe where the discharge flows and the upper (II) has the function of controlling the opening. As relief valve, it should be installed on a derivation of the main pipe in order to put in communication, when it opens, the water in the pipe with the atmosphere, thus reducing the pressures. The lower compartment is divided by a breechblock (1) in two sections: downstream (a) and upstream (b). The former is directly connected to the atmosphere so that, during the efflux, here the pressure is atmospheric but for negligible headlosses. The upper compartment is also divided in two sections: the former (c) above the diaphragm and the latter (d) between the diaphragm and the body of the valve. The diaphragm of the upper compartment and the breechblock of the lower compartment are connected by a rigid rod (2); this connection is not sealed and, therefore, the same pressure can be found in the sections (a) and (d). Actually, the valve is for multipurpose and its use depends on the device called “pilot” (4) that controls its functioning. To use the valve to reduce waterhammer effects, the pilot has to put in communication, either the main pipe or the atmosphere. During steady flow, the upper section is in contact with the main pipe. The pressure in the section is therefore equal to that of the main pipe; the same value of pressure is also in the lower part of the diaphragm (3) and on the upper surface of the breechblock (1) because of the mentioned hydraulic continuity between the


93

(4) II d

c (3)

II a b

(2)

(1)

Diversion to the relief valve

Main pipe

Figure 6.21: Sketch of the valve modelled.

sections (a) and (d), while on the lower surface of the breechblock is the atmospheric pressure. The valve is closed due to the force directed downwards which is applied on the breechblock. When in the pipe the pressure increases, for instance because a downstream valve closure that starts a transient, the pilot, which has minimum inertia, impose the almost instantaneous interruption of the communication between the section (c) and the main pipe, connecting the section (c) with the atmosphere. In this condition, the force on the upper surface of the diaphragm is very small as it is the force applied on the lower surface of the breechblock (1); on the contrary, the force acting on the lower surface of the diaphragm and directed upwards is larger than the force acting on the upper part of the breechblock (1) and directed downwards because of the larger area of the diaphragm (3). Therefore the valve opens. The time required to the valve to be completely open is not negligible because the upper section of the valve is in communication with the atmosphere through a

94 WATER HAMMER SIMULATIONS Hd Qd X1

P3A2

Finertia Diaphragm

P2A2 P2A1 Breechblock

Figure 6.22: Installation of the valve and symbols for modelling.

very small pipe, which creates a resistance to the emptying of this volume. As a consequence, the opening of this valve is not as fast as those purposely designed to this aim. On the other end, its closure is remotely controlled and therefore oscillations (and the possible instabilities of the system, mentioned in Section 6.3) are avoided. The regulation of the opening pressure threshold is obtained by acting on a spring of the pilot. The unknown parameters are the head hc and the discharge outflowing Qc in the upper section of the upper compartment, and the displacement x of the “connection”, i.e., the entity made by the diaphragm, the breechblock and the rod which connects the two. To solve the problem, equations are (i) the dynamic equilibrium of the connection, (ii) the continuity equation for the upper section of the upper compartment and (iii) the equation that ties the outflowing discharge Qc and the internal head hc . With the symbols used in Figure 6.22, the equations above-mentioned are: 2 · λ · (p1 + p) · A1 + p2 · A2 = p2 · A1 + p3 · A2 + FInertia " Qc = µoref · Aoref · 2 · g · hc x = x + dt ·

Qc A2


95

In the equations: p2 · A2 is the force acting under the diaphragm of area A2 and tends to open the valve; p2 · A1 is the force acting on the upper surface of the diaphragm and tends to close the valve; p3 · A2 is the force due to the pressure p3 , which corresponds the head hc , acting on the upper surface of the diaphragm and tends to close the valve; 2 · λ · (p1 + p) is the force due to the pressure p1 + p acting on the lower surface of the breechblock during transients (the meaning and values of the parameter λ will be discussed in the following); FInertia is the inertia force, obtained multiplying the mass of the entity made by the diaphragm, the breechblock and the rod which connects the two by its acceleration; µoref , Aoref are the efflux coefficients and the area of the very small pipe (diameter equal to 5 mm) that connects the upper section and the atmosphere. It is to be considered that p2 = 0 as it is the atmosphere. These equations allow to compute the value of the head in the section (c) and consequently the displacement of the diaphragm and the breechblock, i.e., of the opening of the valve. Once this opening is known, it is possible to compute the effluent discharge by means of the efflux law. Coefficients λ and µoref have to be determined through laboratory tests. Coefficient λ is used to reduce the force on the diaphragm in dependence of the distance between the diaphragm itself and the outlet. It has been estimated in steady flow conditions measuring the force exerted by the jet on the diaphragm at different degrees of opening of the valve. Similar tests have been carried out for the determination of the coefficient µoref , and finding it is also a function of the distance x. Numerical tests have been performed on the three plants shown in Figure 6.23, with different configurations and characteristics, for details refer [54]. An example of the carried out results is reported in Figure 6.24. Generally speaking, the efficacy of these valves when they are installed in the plants has been demonstrated. However, their effect is limited when a very fast response is required; for instance, when the downstream valve closes very fast or when the phase of the plant is very short (low lengths, high celerity, etc.) or even when the required reduction is great. The analysis performed to characterize the valve allowed the authors to build a numerical model of the valve itself, which has been inserted in a number of plants that could not be studied on physical models and that can be used for other simulations. Obviously, all results have to be intended as general orientation for understanding the plant and, moreover, many features cannot be theoretically derived, and experimental tests are very often required, especially when the systems are very complex.


Relief valve

Regulation valve

Regulation valves 1

5

6

7

12

8

2

13

9

3

14

10

4

15

11

Relief valves

Figure 6.23: Plants numerically analysed for determining the characteristics of the valve.

6.6.6 Pumps In the computer code just turning off the pump has been implemented and, therefore, eqn (6.19) where P = 0. The first step for the solution of the problem is the evaluation of the new number of rotation of the pump, which is computed as shown in Figure 6.25; then the working point of the system is carried out using the same procedure adopted for the steady flow, but with the new values of the pump parameters. If the velocity becomes negative, in the computer code it is equal to zero, in order to avoid reverse flow; this position is quite close to the reality, as the reverse flow should be avoided also in the real world and appropriate valves (checkvalves) are installed downstream the pumps in order to prevent it. However, the behaviour of the system since the flow velocity is zero as the section of the pump strongly


97

Upstream valve Downstream valve

160

Heat (m)

120 80 0

0

20

40

60 Time (s)

80

100

Figure 6.24: Numerical results carried out for the plant “A” shown in Figure 6.23.

np:=nn-900*9806*v[0]*Area1*h00/(Efficiency*Inertia*Pi*Pi*nn)*dt; PumpATrans:=PumpA; PumpBTrans:=PumpB*(np/NumRot); PumpCTrans:=PumpC*(np/NumRot)*(np/NumRot); Figure 6.25: Boundary conditions for the pump.

depends on the rest of the plant; therefore, the code that has been developed is to be considered valid up to this point. The plant to be simulated is with a pump positioned upstream; then there is a single pipe whose length is equal to 3000 m, diameter equal to 0.45 m and Manning’s roughness n = 0.0125 m−1/3 · s. Downstream a reservoir whose water elevation is equal to 10 m is positioned. The coefficients of the pump are: a = −12.72 m−1 · s2 b = −1.32 s c = 53 m Inertia of the pump is equal to G · D2 = 12 N · m2 , speed is equal to 2980 rpm, efficiency η = 0.88. The value of inertia is computed with the pump filled with water and also adding the inertia of the motor. The curve for the pump is shown in Figure 6.26, together with the curve of the headlosses of the plant and the working point of the system.

98 WATER HAMMER SIMULATIONS 60 50

Head (m)

40 Working point 30 20 10

Curve of the pump Curve of the plant

0 0

50

100

150 Discharge (l/s)

200

250

300

Figure 6.26: Head balance curve for the pump, coupled with the plant headloss curve. The working point is highlighted. 3500

1.6

Flow velocity GD2 12 Flow velocity GD2 24 Pump rotation GD2 12 Pump rotation GD2 24

1.4

2500

1

2000

0.8 1500 0.6 1000

0.4

500

0.2 0

Rotation (rpm)

Velocity (m/s)

1.2

3000

0 0

0.2

0.4

0.6 Temp (s)

0.8

1

Figure 6.27: Flow velocity and pump rotation for the simulated plant and two different values of inertia. As can be easily seen, in the pipe a discharge of 221 l/s is flowing, with a velocity equal to 1.39 m/s. The simulation is related to the sudden shut off of the energy of the pump, while the pump stops in a time dependent on the inertia and the system. Simulations have been compared with those obtained doubling the inertia.


99

As can be seen in Figure 6.27, the increase of the value of inertia has a visible effect on both the flow velocity and the rotation of the pump; however, the flow velocity reaches zero in a very short time, while the pump continue to rotate, but with no effect. In the figure, continuous lines are related to computed velocities, while dashed lines to the estimated values. Actually, larger pumps can continue their rotation for minutes, but this does not solve the problems due to the waterhammer which effects are created by the velocity of the water inside the conduit. Therefore, it is difficult to say whether the increase of the inertia of the pumps has appreciable effects and can be suggested instead of the most common devices as the air chambers or others, especially because of the problems that they can create; probably no general answer exists.


Chapter 7 Instabilities

As mentioned in the previous chapters, changes from the steady flow conditions can be considered as disturbances that can cause serious problems; with time, these “disturbances” dampen and the flow reaches a new steady regime. Theoretically, this new regime is reached asymptotically, in practice after a finite time. However, in some cases, not only the oscillations does not attenuate, but also do they increase and the system collapse. In the previous chapters such instabilities of the system have been encountered: in Chapter 6 they were related to surge tanks, which require a minimum dimension in order not to give rise to undamped oscillations; then, a real case of resonance was presented, due to valves inappropriately installed. Other cases are presented in the scientific literature, as those found in leaking valves [64] or in the behaviour of check valves, if not correctly installed [43]. All these cases depend on both the system and the external acting force characteristics. The instabilities may not be in the system, but in the computer code, which crashes for numerical problems not related to the actual plant. Again, this case has been encountered in Chapter 6. In this chapter the two cases are addressed.

7.1 Vibrations 7.1.1 General remarks The problems related to vibrations are very complex and textbooks are available to deal with them [65, 66], and therefore in the following only few reminders of theory, which will be used later in the book, are presented. For a pipe with length L made with only one material and with constant diameter and thickness, the fundamental period Tt is [67]: Tt =

4·L c

(7.1)

where, as usual, c is the celerity of the wave. This is also the highest harmonic; for instance, Tt /3 is the period of the third harmonic, Tt /7 of the seventh harmonic and so on. The frequency that corresponds to the fundamental period, or highest harmonic, is called resonance frequency, which refers to an oscillation that leads to a pressure


k1

k2

k3

m1

m2

x1

m3

x2

x3

Figure 7.1: System with three degrees of freedom. amplification. If T is the period, the corresponding frequency ω (rad/s) is equal to: ω=

2·π T

(7.2)

while the circular frequency f (cycle/s) is equal to: f =

1 T

(7.3)

The number of coordinates required to define a system is said degrees of freedom. Figure 7.1 shows a system made of three masses connected in series with three springs. As can be seen, the status of this system can be described with the three coordinates x1 , x2 and x3 and therefore it has three degrees of freedom. A system made of a mass and a spring has a natural frequency ωn equal to: ωn =

" k/m

(7.4)

where, with the notation used in Figure 7.1, k is the elastic constant of the spring and m is the mass. A system with n degrees of freedom has n natural frequencies, and each of them defines a normal mode of vibration. When an oscillatory force with frequency ωf is applied at this system, it oscillates with the frequency of the forcing momentum and an amplitude that depends on the ratio ωf /ωn . If the two frequencies have the same value and no friction is applied, the amplitude of the oscillations tends to infinitive, because at each cycle the system increases its energy, which cannot be dissipated. If friction is considered, the amplitude of the oscillations increases until the energy acquired in each cycle

INSTABILITIES

103

vp[n]:=v0*cos(t*Period); hp[n]:=h[n-1]-cel*(vp[n]-v[n-1])/gcel*Lambda*v[n-1]*abs(v[n-1])*dt/(2*g*diam); Figure 7.2: Code implemented in the computer program to simulate an oscillating downstream velocity. is equal to that dissipated by the friction: therefore the system does not increase its energy and the oscillations have finite amplitude. In the description of hydraulic systems, the fluid is described as made of an infinite number of masses and springs in series, with the springs representing the fluid elasticity. As a consequence, an hydraulic system has infinitive degrees of freedom. This assumption can be noted observing the equations that describe the movement: if a system is governed by an ordinary differential equation, then in this system the number of degree of freedom is equal to the equation order; while, if a system is governed by partial differential equations, the number of degree of freedom is infinitive.

7.1.2 Computer program for oscillating velocity To check and to understand the theory shown in Section 7.1.1, a simple computer program is developed, very similar to that presented in Chapter 4, but having as downstream boundary condition the following: 2·π ·c V [n] = V0 · cos t · N · 4·L

(7.5)

where t is time, 2πc/4L is the fundamental frequency of the system and N is a parameter to adjust the period of the oscillating velocity. The negative characteristic line is used to compute the head. Obviously, the condition does not represent an actual system or valve, but its function is purely educational. The code is shown in Figure 7.2. As application, let us consider the case of a plant with an upstream reservoir (the head is constant and equal to 100 m) a pipe 10,000 m length and 0.3 m diameter, celerity equal to 1,000 m/s discharge equal to 70 l/s; downstream a valve is positioned. Headlosses are negligible. Two simulations have been carried out. In the former, the frequency of the oscillations ω is let equal to the natural frequency of the system, i.e.: ω = 2·π·c = π5 . 4·L In the latter simulation, the frequency of the downstream valve is equal to the half of the above computed, i.e.: ω = π/10.

104 WATER HAMMER SIMULATIONS 3000 Force in resonance Force in non resonance

2000

Head (m)

1000

0 0

50

100

150

200

250

300

1000

2000

3000

Heads below 10.33 m are not physically possible because of column rupture. They are shown only to demonstrate the effect of the resonance. Time (s)

Figure 7.3: Results of simulations with oscillating downstream head. Results of simulations are shown in Figure 7.3; obviously, when the pressure drops below −10.33 m the water column breaks and the actual behaviour is different from that seen in the figure (see Chapter 8), but the simulation has not been interrupted because the target was the visualization of the behaviour of a system in resonance. As can be seen, when the forcing momentum has a frequency of ω = π/10 (this is obtained inserting in the user interface the parameter N Periods = 0.5), oscillations are limited and have a period equal to 40 s; when the forcing momentum has a frequency equal to the natural frequency of the system, i.e., ω = π/5 (and this is obtained inserting in the user interface the parameter N Periods = 1.0), the pressures increase more and more. In the simulations, the oscillations are not damped having accepted the hypothesis of negligible headlosses. However, if the forces that can contrast the increase of the oscillations are limited, heads increase and the system eventually breaks.

7.2 Transfer matrix method The case presented in the Section 7.1.2 is very simple and the result of the simulation was known in advance from vibration theory. Obviously, and unfortunately, real cases are much more complex. Different methods have been developed to study this kind of problems. Numerical models as revised in this book may not be able to clarify the problem, because

INSTABILITIES

105

it is not always clear whether unexpected oscillations that are noted during simulations have physical or numerical origin (in other words, if it is a problem of the computer program, or the plant can actually develop unwanted oscillations). When the origin of the oscillations is not clear, a different method is used to evaluate the system stability, which develops overall matrices for a fluid system [68]; this method is widely used also in mechanics, structural and electrical engineering. This method works in the frequency domain, linearizing the equations that govern the phenomenon and considering sinusoidal both the motion and the pressure fluctuations. 7.2.1 General remarks Generally speaking, a linear system which input variables x1 , x2 , …, xn are tied to the output variables y1 , y2 , . . ., yn by the equations: ⎧ ⎪ ⎪ u11 · x1 + u12 · x2 + · · · + u1n · xn = y1 ⎨ u21 · x1 + u22 · x2 + · · · + u2n · xn = y2 ··· ⎪ ⎪ ⎩ un1 · x1 + un2 · x2 + · · · + unn · xn = yn

(7.5)

in matrix notation can be written as: · x y = U

(7.6)

where x and y are the vectors that contain the variables, and they are called state contains the parameters of the system and it is called transfer vectors; the matrix U matrix. In the hydraulic systems the variables to be transferred are, as mentioned, velocities (or discharges) and heads (or pressures). As a consequence, in this case the input data are the upstream velocity and head, while the output data are the downstream velocity and head, i.e.: V V zinput = upstream zoutput = downstream hupstream hdownstream

(7.7)

and, therefore, the system is governed by the equation: · z zoutput = U input

(7.8)

Obviously, systems can be composed by combinations of series or parallel pipes, valves, junctions, pumps and all the devices are only partially described in the book. For each of these subsystems, an equation has to be written based on the transfer matrix; then, the description of the whole system is carried out joining all these subsystems in a larger matrix.

106 WATER HAMMER SIMULATIONS 7.2.2 Application to hydraulic transients The governing equations have been shown in Chapter 2; in the following they are written again for the sake of simplicity and in the form which will be used in this chapter [11]. Continuity equation: g·A · Ht + Q x = 0 c2

(7.9)

Momentum equation: Hx +

1 Qn =0 · Qt + f · g·A 2 · g · D · An

(7.10)

In eqn (7.10) the exponent n depends on the headlosses. As can be seen, in these equations the unknown variables are the heads H (instead of the pressures p) and the discharges Q (instead of the velocities V ); moreover, with a more concise notation, the generic partial derivative has been indicated: Ab =

∂A ∂b

(7.11)

The hypothesis that the transient is developed around average values of head H and discharge Q is introduced. Therefore the following equations can be written: H = H +h

(7.12)

Q = Q+q

(7.13)

where h and q are the oscillation of the head and discharge values around the steady flow values. Deriving with respect to the space and the time, and remembering that in steady flow the discharge is constant in the space and time, while the head is constant only with respect to the time since it changes in the space (because of the headlosses), it yields: H x = Hx + h x = − f ·

Qn + hx 2 · g · D · An

(7.14)

Ht = Ht + ht = ht

(7.15)

Qx = Qx + qx = qx

(7.16)

Qt = Qt + qt = qt

(7.17)

The term related to the headlosses can be written as follows: n

f ·

n−1

Q Q (Q + q)n ∼ Qn =f · +f ·n· ·q =f · n 2·g·D·A 2 · g · D · An 2 · g · D · An 2 · g · D · An (7.18)

107

INSTABILITIES

Substituting these equations in eqns (7.9) and (7.10) yields: g·A · ht + qx = 0 c2

(7.19)

n−1

hx +

Q 1 · qt + f · n · ·q = 0 g·A 2 · g · D · An

(7.20)

Introducing the terms: 1 g·A g·A C= 2 c L=

inertance

(7.21)

capacitance

(7.22)

resistance

(7.23)

n−1

R = f ·n·

Q 2 · g · D · An

The system finally becomes: qx + C · ht = 0

(7.24)

hx + L · qt + R · q = 0

(7.25)

This system can be solved with the method of separation of the variables that consists in deriving the first equation with respect to the time t and the second with respect to the space x, and then to subtract them from each other; afterwards, the process is repeated inverting the order, i.e., deriving the first equation with respect to x and the second with respect to t. At the end of the whole process the following equations are obtained: qxx = C · L · qtt + C · R · qt

(7.26)

hxx = C · L · htt + C · R · ht

(7.27)

A new hypothesis is to be introduced: that the variable h can be computed as the product of two independent functions X = X (x) and T = T (t); if this hypothesis is verified, the integral of the system can be computed. Substituting, the following equations can be carried out: 1 d2 X = γ2 · X dx2 1 d2 T dT · C·L· 2 +C·R· = γ2 T dt dt

(7.28) (7.29)

where γ 2 is a complex value independent from x and t. Integrating it yields: X = A1 · eγ ·x + A2 · e−γ ·x

(7.30)

T = A3 · es·t

(7.31)

108 WATER HAMMER SIMULATIONS where s is again a complex value independent from x and t, while A1 , A2 and A3 are the integration constants. Now it is possible to tie the variables γ and s: γ 2 = C · s · (L · s + R)

(7.32)

Considering that h(x, t) = X (x) · T (t), it yields: h = es·t · (C1 · eγ ·x + C2 · e−γ ·x ) = H (x) · es·t

(7.33)

and finally it is possible to define the expression: q=−

C · s st · e · (C1 · eγ ·x − C2 · e−γ ·x ) = Q(x) · es·t γ

(7.34)

The variable γ is called propagation constant and represents the complex frequency, i.e.: s=σ +i·ω

(7.35)

where σ is referred to the damping of the system and ω is the oscillation circular frequency. Finally, the characteristic impedance can be defined as: Zc =

γ C·s

(7.36)

This impedance is a complex function independent of space and time; this function and the propagation constant are functions of fluid and pipe characteristics and the complex frequency. The functions H (x) and Q(x) are the complex head and the complex discharge, respectively, and depend on the spatial variable only. With these equations it is possible to define the values of the integration constants function of the boundary conditions. Once these constants are known, transfer matrix can be built for each element of the system. 7.2.3

Description of simple system: pipes

In the simple case of a single pipe, represented in Figure 7.4, the hydraulic variables have to be transferred from upstream to downstream. As when x = 0 it must be H = Hupstream and Q = Qupstream , from the definition of H (x) and Q(x), given in eqns (7.33) and (7.34), the values of the integrations constants C1 and C2 are carried out: 1 · (Hupstream − Zc · Qupstream ) 2 1 C2 = · (Hupstream + Zc · Qupstream ) 2

C1 =

(7.37) (7.38)

INSTABILITIES

109

H downstream

H upstream Q upstream

Q downstream

1

Figure 7.4: Sketch of a simple pipe and related symbols. Finally, substitution yields: H (x) = Hupstream · cosh(γ · x) − Zc · Qupstream · sinh(γ · x) Q(x) = −

Hupstream · sinh(γ · x) + Qupstream · cosh(γ · x) Zc

(7.39) (7.40)

To compute the values of the complex head and discharge downstream the pipe, the simple substitution of x = l allows are to obtain: Hdownstream = Hupstream · cosh(γ · l) − Zc · Qupstream · sinh(γ · l) Qdownstream = −

Hupstream · sinh(γ · l) + Qupstream · cosh(γ · l) Zc

(7.41) (7.42)

These relations allow to write the transfer matrix for the element “pipe”, according to the notations in Figure 7.4, indicating as positive the discharges outflowing from each element and introducing the functions: −1 Zc · tanh(γ · l) −1 S= Zc · sinh(γ · l)

T =

(7.43) (7.44)

The final compact transferring relation is: T S Hupstream Qupstream · = S T Hdownstream Qdownstream 7.2.4

(7.45)

Description of simple system: valves and effluxes

Transfer matrixes in these cases are obtained linearizing the efflux equations, which does not imply a severe error if the pressure differences induced by the valve are small compared with the static head.

110 WATER HAMMER SIMULATIONS If valves have a dynamic behaviour, this is assumed oscillatory sinusoidal; this method is obviously possible to study non-periodical behaviours too, decomposing the motion in a series of harmonics through the Fourier analysis, which, however, will not be described in this book [69]. In the case of a free efflux (i.e., in the atmosphere), say H the head upstream the efflux and Q the outflow discharge, the well-known formula can be used: Q = µ · Aefflux ·

"

2·g·H

(7.46)

where µ is the efflux coefficient. Recalling eqns (7.12) and (7.13) and letting τ = µ · A, it yields: Q Q that becomes:

=

Q+q Q

τ +τ = · τ

H +h

1/2

H

τ h 1/2 1+ = 1+ · 1+ τ Q H q

(7.47)

(7.48)

Explaining h, the formula becomes: h=

τ2 2·q q2 · 1+ + 2 −1 ·H τ2 + 2 · τ · τ + τ2 Q Q

(7.49)

2

Neglecting the terms τ 2 and q2 /Q and rearranging, eqn (7.49) becomes: h=

2·τ 2·τ q ·H − · ·H τ +2·τ Q τ +2·τ

(7.50)

When τ << τ , eqn (7.50) can be simplified as: h=2·

H Q

·q−

2·τ ·H τ

(7.51)

In Section 7.2.2, eqns (7.33) and (7.34) have been developed; considering only their real parts the following are derived: h = H (x) · sin(ω · t)

(7.52)

q = Q(x) · sin(ω · t)

(7.53)

The hypothesis that the valve has an oscillatory sinusoidal behaviour means that: τ = k · sin(ω · t)

(7.54)

INSTABILITIES

111

where k is the amplitude of the valve motion. With due simplifications the following is obtained: Hupstream = 2 ·

H Q

· Qupstream −

2·k ·H τ

(7.55)

Because of the hypothesis of free efflux (i.e., directly in the atmosphere), Hdownstream = 0, and therefore: Hdownstream = Hupstream − 2 ·

H Q

· Qupstream +

2·k ·H τ

(7.56)

Finally, from the continuity equation: Qupstream = Qdownstream

(7.57)

In conclusion, denoting with: U = −2 ·

V = −2 · the transfer matrix is:

7.2.5

H Q

H ·k τ

Qdownstream 1 0 0 Qupstream Hdownstream = U 1 V · Hupstream 0 0 1 1 1

(7.58)

(7.59)

(7.60)

Global matrix of a system

Let us consider the usual very simple system, as in Figure 7.5, composed by three elements in series: • an upstream reservoir; • a pipe; • a valve with free efflux. For the pipe (Section 7.2.3): T S H A QA S T · H B = QB pipe pipe and for the valve (Section 7.2.4): valve Q downstream 1 valve H downstream = U 0 1

valve 0 Qupstream valve 1 V · Hupstream 0 1 1

(7.61)

0

(7.62)


Upstream reservoir

Downstream valve

Figure 7.5: Sketch of a simple plant. It has to be noted that, in the case of free efflux, obviously is: valve Hdownstream = Hatm = 0

(7.63)

Moreover, the upstream head can be denoted as: valve B H B = Hupstream = Hpipe

(7.64)

Finally, from the continuity equation: B valve Qpipe = Qupstream

(7.65)

Developing the matrixes contained in eqn (7.63) and (7.64), keeping in account the (7.65) and collecting the common terms, the following equation is carried out: V U HB = − S · HA + · (7.66) U T ·U +1 With the matrix notation, eqn (7.66) can be expressed as: B −S·U −V A H T ·U +1 T ·U +1 H = · 1 0 1 1

(7.67)

In order to study the free vibrations, then, being negligible the external causes, it is necessary the analysis of eqn (7.67). The represented system allows non-trivial solutions if the determinant of the matrix is zero; therefore, it is to be computed the complex variable s = σ + i · ω such that: S ·U =0 T ·U +1

(7.68)

The purpose of the analysis is to identify if there are values of the variable σ function of ω with positive values: in this case, the system is unstable.

INSTABILITIES

113

7.2.6 A simple application The plant to be studied is that described and tested in Section 7.1.2, i.e., a pipe 10,000 m length and 0.3 m diameter, celerity equal to 1000 m/s where flows a discharge equal to 70 l/s. Headlosses are negligible and the upstream head is equal to 100 m. The solutions of eqn (7.68) have to be studied. First, the variables inertance L (eqn (7.21)), capacitance C (eqn (7.22)) and resistance R (eqn (7.23)) are computed (the latter term is zero because of the hypothesis of negligible headlosses): L=

1 = 1.442 g·A

g·A = 6.934 · 10−7 c2 To compute the characteristic impedance of the system, eqns (7.32) and (7.36) have to be used: " γ = C · s · (L · s + R) = 9.978 · 10−4 · s C=

γ = 1438.975 C·s It is to be reminded that in the above equations, s is the free variable. To study the behaviour of the pipe, eqns (7.43) and (7.44) are needed: Zc =

S=

−1 −1 = Zc · sinh(γ · l) 1438.975 · sinh(9.978 · s)

T =

−1 −1 = Zc · tanh(γ · l) 1438.975 · tanh(9.978 · s)

Because of the valve, eqn (7.58) is to be introduced: U = −2 ·

H Q

= −2857.143

Eqn (7.68) in the complex variable s = σ + i · ω can finally be written as: 1.986 + tanh(9.978 · s) −

1.986 =0 cosh(9.978 · s)

The study of this function is possible with any spreadsheet and little effort. It can be shown that no positive values of σ exists to set the equation equal to zero: this result means that the system is stable. As mentioned, this is true when no external

114 WATER HAMMER SIMULATIONS force are considered; on the contrary, the presence of an external force oscillating with the natural frequency of the system would bring it in resonance. In conclusion, to decide whether a system is stable, different methods should be applied. However, a good practice is always the installation of slow-response devices, which reduce the risk of undesired instabilities.

7.3 Numerical instabilities In the case of waterhammer, it may happen that the computer program crashes even when the real plant under analysis is not problematic; this is often to be ascribed to numerical instabilities. Numerical schemes based on the series of Taylor, as the Lax–Wendroff ’s, can be analysed developing the series. As it has been shown in Chapter 5, the solution for a linear system like the following: →

→

∂u ∂u +A· =0 ∂t ∂x

(7.69)

ut + A · u x = 0

(7.70)

or, in compact form:

can be expressed as: u(x, tn+1 ) = u(x, tn ) + t · ut (x, tn ) +

1 · (t)2 · utt (x, tn ) 2

(7.71)

Equation (7.71) shows that the series of Taylor has been stopped at the second order and this implies that the truncation error is given by the following neglected term, which is: uttt = −A3 · uxxx

(7.72)

To study the behaviour of the solution, it is possible to look closely to the original equation (7.70) with a slight modification in order to correlate it with the truncation error. Hedstrom [70] shows that the modified equation for the Lax– Wendroff method is: ut + A · ux = µ · uxxx − ε · uxxxx

(7.73)

where: (x)2 µ= ·A· 6 ε=

(t)2 · A2 − 1 (x)2

(x)3 (t)2 2 ·A· 1− · A 6 (x)2

(7.74) (7.75)

INSTABILITIES

115

3 2,5 2 1,5 1 0,5 0 0,5 1 1,5 2 2.5 2 1.5 1 0.5

0

0.5

1

1.5

2

2.5

Figure 7.6: Example of the instabilities. In eqn (7.73), the two addends have different physical meaning. Generally speaking, even spatial derivatives as ε · uxxxx are related to dissipation (or amplification) terms, while odd terms are related to dispersion. The term ε · uxxxx modifies the trend of the solution when |c · λ| < 1 and that confirms the need to let CFL < 1 and at the same time close to one, in order to have stable and non-dissipative solutions. If the spatial term can be neglected, eqn (7.73) is reduced to: ut + A · ux = µ · uxxx

(7.76)

which is called dispersive equation and allows to explain the formation of oscillations close to the points where the presence of large gradients of the dependable variables is noted, as shown in the Figure 7.6 where a square wave is integrated using Lax–Wendroff method. The effects of diffusion have been encountered and highlighted in Chapter 4 (see Figure 4.11), but without an appropriate analysis, which will be performed in the following. In this book the analysis of the reasons why the instabilities occur is not performed, as more specialized books are available in the scientific literature [31]; however, as this problem is found quite frequently, some methods to solve it, or at least reduce it, are shown. In particular, the following will be analysed: • • • •

the change of the CFL parameter; the use of a first order code; the use of flux limiters; the use of artificial viscosity.

116 WATER HAMMER SIMULATIONS Diffusion errors Lax–Wendroff 1.2 1 0.8

CFL 0.25

0.6

CFL 0.50

0.4

CFL 0.75

0.2

CFL 1.00

0 0

20

40

60

80

100

120

140

160

180

140

160

180

Dispersion errors Lax–Wendroff

1 0.5 0 0

20

40

60

80

100

120

k ∆x

Figure 7.7: Diffusion and dispersion errors, function of CFN and angle of phase parameters (after Ref. [71]). 7.3.1 Changing CFL number As discussed in Chapter 5 and mentioned above, the easiest way to try to reduce the numerical oscillations is to keep the CFL number close but lower than 1; the consequences of the variations of this number have already been shown in Figure 5.12. In the paragraph above (Section 7.3) a brief analytical description of the oscillation problem has been given through the modified equation (7.73), introducing a dissipative (or amplifying) and a dispersive term. Changing the CFL number implies a change in the numerical value of these parameters with obvious implications on the stability and convergence of the solution. Figure 7.7 shows [71] the relative values of the parameters ε and µ when the CFL number changes, for the Lax–Wendroff method; on the x-axis the phase angle of the solution1 is reported. As can be seen, in the case of the Lax–Wendroff method, the only value that guarantees no dispersion along the whole axes of the pipe is CFL = 1, as already mentioned. When the phase angle increases, for CFL values lower than 1, the dissipative and dispersive terms become more important.

1

The phase angle of the kth harmonic is a real number expressed in radiants and equal to k · x.

INSTABILITIES

117

150 125

Head[s] (m) Head

100

4000

75 50

3000

25 0

2000

0.0

0.5

1.0

1.5

2.0 2.5 Time (s)

3.0

3.5

Head (m)

1000 0 0

5

10

15

20

25

30

35

40

45

50

1000 2000 3000 4000 Time (s)

Figure 7.8: Results with boundary conditions solved with the Beam–Warming method and internal points with Lax–Wendroff (CFL = 1). It has also been already mentioned that the condition on the CFL number is not sufficient to guarantee stability: Figure 7.8 allows to focus on a particular case. Results have been carried out with the usual simple plant where the valve is instantaneously closed; the pipe length is equal to 10 000 m and the celerity is 1000 m/s. Therefore, the perturbation starts at time 0 s and it is back every 20 s; as can be seen from Figure 7.8, when the perturbation is back to the valve, the instabilities increase until the computer program crashes. The problem in this case is obviously given at the boundary conditions, which have been solved with the Beam–Warming scheme and cannot be solved operating on the CFL number. 7.3.2 First order methods As already described, the numerical instabilities are due to the truncation errors and are part of the intrinsic diffusive properties of the selected numerical method of solution. Usually, if the main terms of the error are odds, as it happens in the second order schemes, these dispersive errors generate the oscillations. A viable

118 WATER HAMMER SIMULATIONS 1.5

1

0.5

0

0.5 0.2

0

0.4

0.6

0.8

1

Figure 7.9: Numerical diffusion in the integration through first-order schemes. alternative is therefore the reduction of the order of the numerical scheme, using a first order method, for instance the mentioned Lax–Friedrichs. LeVeque [36] illustrates the results carried out with first order schemes although they do not show non-physical oscillations in the presence of discontinuities, and they are less accurate with respect to second order methods. As can be seen in Figure 7.9, in fact, where a generic case of first-order integration is shown, the solution is approximated, but the numerical solution tends to be smoother than the analytical because of the numerical dispersion. In this case, the equivalent equation at the third order is [71]:

where:

ut + A · ux = D · uxx − ν · uxxx

(7.77)

(t)2 (h)2 · A2 · 1− D= 2 · t (h)2

(7.78)

(h)3 (t)2 · A2 · 1− 3 (h)2

(7.79)

ν=

The term D · uxxx is dissipative; in other words, increasing the positive parameter D, the solution tends to flatten, becoming increasingly smoothed where the first derivative is discontinuous. The term µ · uxxx is normally not negligible, although the dispersive effect is usually not really noticeable, because hidden by the dissipative, much larger. As the numerical dispersion is quite small, no oscillations are recorded; however, an increase of the celerity when µ > 0 and a diminution when µ < 0 can be noted. As a consequence, observing the position of the wave crests of Figure 7.9 compared with those of the exact solution, it can be seen that the Lax– Friedrichs method has, generally speaking, a positive dispersive error: as µ > 0 the celerity changes and the carried out solution is faster than the exact solution.

INSTABILITIES

119

Diffusion errors Lax–Friedrichs 1.2 1 0.8 0.6 0.4 0.2 0 0

20

40

60

80

100

120

140

160

180

140

160

180

Dispersion errors Lax–Friedrichs 5

CFL 0.25

4

CFL 0.50

3

CFL 0.75

2

CFL 1.00

1 0

20

40

60

80

100

120

k ∆x

Figure 7.10: Amplification and dispersion errors, function of CFN and angle of phase parameters (after Ref. [71]). In Figure 7.10 for the method of Lax–Friedrichs the values of the coefficients D and ν are reported, when changing the phase angle and the CFL number. Letting the CFL number equal to 1, the centered Lax–Friedrichs method can be considered a good choice; also in this case, the CFL condition is not sufficient to guarantee a solution free of errors of amplification or dispersion, and therefore, the use of a first-order scheme is to be deeply evaluated. The absence of spurious oscillations in fact, in general could be compensated by the presence of numerical errors of different nature such as dispersion. 7.3.3 Flux-limiters A very interesting option to be used to deal with these problems is called “FluxLimiter”; this expression joins a number of methods that generally guarantee good results. The common ground of all these methods is the insertion in the numerical code of a new term which allows to obtain a second-order accuracy when it is possible, but that smooths it when the numerical solution becomes unstable. The general idea is the construction of an hybrid numerical flux, carried out combining a second-order scheme with a monotonic first-order scheme so that the former or the latter prevails when the solution is regular or irregular, respectively.

120 WATER HAMMER SIMULATIONS One of these methods is the so-called total variation diminishing (TVD) [72]; the discretized total variation of a solution u for a scalar conservation law is: # TV (u) = |ui+1 − ui | (7.80) i

and a numerical scheme is called TVD if for any time and any u the following inequality is satisfied: TV (un+1 ) ≤ TV (un )

(7.81)

When the inequality (7.81) is met, the scheme respects the monotonic trend of a function, or, in other words: • in the computational domain no new extremes can be created; • the value of a local minimum is not decreasing, while the value of a local maximum is not increasing. Obviously, within a scheme with such properties, starting from a correct solution, non-physical oscillations cannot be created. Godunov [73] demonstrated that linear first-order TVD schemes do not exist, and therefore in a second order scheme it is necessary to introduct non-linear correction terms. These terms are known as “limiters”, as they “limit” the flux between the cells. Generally speaking, a flux-limiter scheme can be expressed as: H = HL + (HH − HL )

(7.82)

where HL is a first order monotonic flux (“Lower”), HH is a second order flux (“Higher”), is a flux limiter. It is desirable that the function tends to 1 when the solution is regular and that tends to 0 when the solution is discontinuous; moreover, it must guarantee that the scheme converge to the right solution. After having arbitrarily defined HL and HH , the shape of the flux limiter is defined accordingly. As example, HL can be defined as the Upwind first-order numeric flux: n ujn+1 = ujn − Anj · λ · (ujn − uj−1 )

(7.83)

while as HH , the second-order Lax–Wendroff scheme, can be used: ujn+1 = ujn −

1 1 n n n n n + · A2 · λ2 · uj+1 − uj−1 − 2 · ujn + uj−1 · Aj · λ · uj+1 2 2 (7.84)

The parameters have the usual meaning and in particular λ = t/x. The hybrid flux is, therefore: Q(u, j) = Anj · ujn +

1 n n n · A · (1 − λ) · (uj+1 − uj−1 ) · (ϑjn ) 2 j

(7.85)

INSTABILITIES Φ

121

Φ 2θ

Φ2

2

θ

Figure 7.11: Admissibility region for a flux-limiter method. where ϑjn is an index related to the regularity of the solution and it is defined as: ϑjn =

n ujn − uj−1

(7.86)

n uj+1 − ujn

Therefore, when ϑjn → 1 data are presumably regular, while when ϑjn → 0 is likely the presence of discontinuities. LeVeque [36] and Hirsch [31] show that, for the function (ϑjn ) to be consistent and to allow the convergence towards the real solution it is necessary the function is confined in an admissibility region, defined by eqn (7.87) and shown in Figure 7.11. (ϑ n ) j n − (ϑ ) j−1 ≤ 2 ϑjn

n ∀ϑjn , ϑj−1

(7.87)

The inequality (7.87) can be also expressed as: 0 ≤ (ϑjn ) ≤ 2 · ϑjn ∧ 0 ≤ (ϑjn ) ≤ 2

(7.88)

According with Ref. [31, 36], for the method to be second-order accurate the following conditions have to be satisfied as well: • the gradient of the function (ϑjn ) is limited; • (1) = 1; • is continuous in ϑjn = 1. Combining Lax–Wendroff method with the boundary conditions given by Beam–Warming, for which (ϑ) = ϑ, it is possible to definite a new, second-order, admissibility region, shown in Figure 7.12.

122 WATER HAMMER SIMULATIONS Φ

Φ 2θ

Φθ

2

Φ2

1

Φ1

θ

Figure 7.12: Second-order region. Any curve chosen for the function (ϑ) which falls in this region guarantees that the carried out second-order method is converging towards the exact solution. Solutions carried out, even when discontinuities of the first derivative are present, are oscillations free and give a good approximation of the real value. The most common limiters are the following, represented also in Figure 7.13 to show their position inside the above defined second-order region. • Superbee limiter: (ϑ) = max(0, min(1; 2ϑ), min(ϑ; 2))

(7.89)

• Van Leer limiter: (ϑ) =

|ϑ| + ϑ |ϑ| + 1

(7.90)

• MinMod limiter: (ϑ) = max(0, min(1; ϑ))

(7.91)

The described limiters, only a sample from those available in literature, allows to modify the solution depending on the problem and on the implemented integration scheme. Therefore, it is not possible an a priori suggestion of the “best” limiter that guarantee proper solutions for each analyzed case, but it is necessary to make that choice as a result of an analysis of the problem. 7.3.4 Artificial dissipation To improve their solution, Lax–Wendroff [33] developed this method, which consists in the insertion in the numerical code of terms that model an artificial viscosity,

INSTABILITIES Φ

Φ 2θ

123

Φθ

2

Φ2

1

Φ1

θ

Figure 7.13: Flux-limiters: Superbee (dash dot), Van Leer (continuous) and MinMod (dashed). or dissipation, directly derived from the discretization of the equations. These terms have the effect to increase the viscosity around the critical points, although they are negligible in the regions where no oscillations are produced and the solution is regular. Reminding the base scheme of the Lax–Wendroff method: ujn+1 = ujn − A ·

λ n λ2 n n n )+ − 2ujn + uj−1 ) (uj+1 − uj−1 · A2 (uj+1 2 2

(7.92)

which can be written in the general form as: n n − Hj−1/2 ) ujn+1 = ujn − λ(Hj+1/2

(7.93)

n n = h(ujn ; uj+1 ) Hj+1/2

(7.94)

n n Hj−1/2 = h(ujn ; uj−1 )

(7.95)

where

are the numerical fluxes, and in this case they can be expressed as: n = Hj+1/2

n n · A · (uj+1 − uj ) − λ · (A2 · (uj+1 − uj ))

(7.96)

n Hj−1/2

n n − uj ) − λ · (A2 · (uj−1 − uj )) · A · (uj−1

(7.97)

1 2 1 = 2

124 WATER HAMMER SIMULATIONS vp[i]:=v[i]-0.5*(dt/ds)*((h[i+1]-h[i-1])*g+v[i]*(v[i+1]-v[i-1]))+ 0.5*(dt/ds)*(dt/ds)*(2*g*v[i]*(h[i+1]-2*h[i]+h[i-1])+ (cel*cel+v[i]*v[i])*(v[i+1]-2*v[i]+v[i-1]))Lambda*v[i]*abs(v[i])*dt/(2*diam); if abs(h[i+1]+h[i]+h[i-1])>0.01 then begin CorrV:=eps*(dt/(2*ds))*(abs(v[i])+cel)*abs(v[i+1]-2*v[i]+v[i-1])* (v[i+1]-v[i-1])/(h[i+1]+h[i]+h[i-1]); vp[i]:=vp[i]-CorrV; end; hp[i]:=h[i]-0.5*(dt/ds)*(v[i]*(h[i+1]-h[i-1])+(Cel*Cel/g)*(v[i+1]-v[i-1]))+ 0.5*(dt/ds)*(dt/ds)*((v[i]*v[i]+Cel*Cel)*(h[i+1]-2*h[i]+h[i-1])+ (2*Cel*Cel*v[i]/g)*(v[i+1]-2*v[i]+v[i-1])); if abs(h[i+1]+h[i]+h[i-1])>0.01 then begin CorrH:=eps*(dt/(2*ds))*(abs(v[i])+cel)*abs(h[i+1]-2*h[i]+h[i-1])* (h[i+1]-h[i-1])/(h[i+1]+h[i]+h[i-1]); hp[i]:=hp[i]-CorrH; end; Figure 7.14: Computer code for Lax–Wendroff method with artificial viscosity. The addition of the viscosity coefficient changes the eqns (7.96) and (7.97) as follows: n Hj+1/2 =

n Hj−1/2 =

1 n n n n · [A · (uj+1 − uj ) − λ · (A2 · (uj+1 − uj ))] − D(ujn ; uj+1 ) · (uj+1 − uj ) 2 (7.98) 1 n n n − uj ) − λ · (A2 · (uj−1 − uj ))] − D(ujn ; uj−1 ) · (ujn − uj−1 ) · [A · (uj−1 2 (7.99)

The function D represents the artificial viscosity. Inserting that term in the general Lax–Wendroff scheme, the expression is: ujn+1 = ujn − A · +

λ λ2 n n n n − uj−1 )+ − 2 · ujn + uj−1 ) · (uj+1 · A2 · (uj+1 2 2

λ n n n ) · (uj+1 − uj ) − D(ujn ; uj−1 ) · (ujn − uj−1 )] · [D(ujn ; uj+1 2

(7.100)

In order to have a stabilizing effect, function D has to be positive. In practice, n n such term is expressed through a polynomial function of (uj+1 − ujn ) and (ujn − uj−1 ).

INSTABILITIES

125

120

100

Head (m)

80

60

40 Air chamber without orifice Air chamber with orifice

20

0 0

10

20

30

40

50

60

Time (s)

Figure 7.15: Effects of artificial viscosity for the case of air chamber with orifice. Many alternatives are possible [74–76]. In particular, MacCormack and Baldwin [76] proposed the following definition: |V | + c ∂ 2 u (7.101) · 2 D = ε · x2· · p ∂x Finally, the Lax–Wendroff scheme becomes: λ λ2 n n n n − uj−1 )+ − 2 · ujn + uj−1 ) · (uj+1 · A2 · (uj+1 2 2 n n |uj+1 − 2 · ujn + uj−1 | λ n n − · ε · (|V | + c) · · (uj+1 − uj−1 ) (7.102) n n n 2 hj+1 + hj + hj−1

ujn+1 = ujn − A ·

Again there is the problem of a correct determination of the coefficient ε, which if too low the artificial viscosity has no effect and if too high the computer program crashes; therefore, a fine tuning is required and depends on the case to be solved. Modelling transients is never a trivial problem, but when oscillations are present it becomes almost an art. However, the correction has been inserted in the computer program, as shown in Figure 7.14, and the results for the case of air chamber with and without orifice are shown in Figure 7.15. As can be reminded, in Section 6.6.3 (see Figure 6.18), instabilities were found using the method of characteristics to try to simulate the presence of an air chamber with orifice. With the use of the Lax–Wendroff ’s method and limiters, the problem is completely solved.


Chapter 8 Effects of air and cavitation

During transients it may happen that pressures fall below atmospheric, especially when steady flow pressures are low; sometimes pressures can fall even down to the vapour pressure related to the local temperature of the fluid with consequent formation of cavities1 . The pressure inside these cavities is not larger than the vapour pressure mentioned above, if no air is artificially inserted. The formation of those cavities can bring two types of danger to the plant. In the former, the danger of crushing for excess of external pressure and the rupture of parts particularly delicate (for instance, in some treatment water plants, some kinds of membranes). In the latter, in the subsequent phases of the phenomenon, the cavities can abruptly collapse, which can produce very high pressures. Moreover, the effect of the presence of air is noted in the transients for two different cases: the change of the value of the celerity of the propagation wave and an increased flow resistance. Finally, the cavitation can become so large to produce the so-called column rupture, which effects are significant and very difficult to be reproduced with a numerical model. All these effects are briefly discussed in this chapter. However, this topic is still subject to research and many aspects need to be disclosed.

8.1 Cavities 8.1.1 Formation of the cavities As mentioned, when the pressures goes below the value of vapour pressure, first we have the formation of bubbles and then the phenomenon of water column rupture.

1

Theoretically, the absolute pressure of the fluid can fall to zero and only then the rupture of the column should happen; in fact, absolute pressure equal to zero means that normal stresses are zero; positive absolute pressure means that normal stress are of compression and negative pressures means that the stresses are of traction. As water (and all fluids, generally speaking) breaks under traction, negative absolute pressures are not possible. However, in the real world, before reaching this theoretical limit, the dissolved air in the water frees and, therefore, cavities can form much earlier.

128 WATER HAMMER SIMULATIONS In case of equilibrium, the Henry’s law [77] applies: C/p = constant

(8.1)

where C is the gas concentration and p, is the pressure. Obviously, during transients the equilibrium is never reached: gas tends to develop preferentially from free surfaces (in this case, for instance, from existing bubbles) and a diffusive phenomenon starts, where concentration is function of position and time. The governing equation is, therefore, of the type2 : →

←

D · ∇ 2 C + v × grad C =

∂C ∂t

(8.2)

where D is the diffusivity constant of the liquid. Boundary conditions are quite difficult to be set and the equation itself is normally not analytically integrable. Usually two types of conditions are placed where the water column breaks. If no air is inserted, the simplest condition is to consider the cavity as a point where the pressure is constant and equal to vapour pressure. A more sophisticated way to deal with the bubbles is the hypothesis that the gas inside them follows the polytropic law: p · W m = constant

(8.3)

where W is the volume of gas and m is the exponent of the law, which falls in the range from 1.0 (for isotermic transformations) to 1.4 (for adiabatic). This way it is possible to study the bubbles stability in function of their dimensions. Because of the surface tension and imposing the same pressure inside and outside a spherical bubble, the equation that provides the pressure p of the liquid outside the bubble with radius r is: p = p0 ·

r 3 m 0

r

−

2·A + pv r

(8.4)

where pv is the partial pressure of the vapour inside the bubble, A is the constant of surface tension and the subscript 0 refers to the steady conditions. Equation (8.4) is represented in Figure 8.1. As can be seen in Figure 8.1, pressure decreases to a minimum given by the condition: dp p0 r0 3m 2 · A = 0 ⇒ −3m · · + 2 =0 (8.5) dr r r r

2

It is to be reminded that ∇ 2 C =

∂2C ∂x2

+

∂2C ∂y2

+

∂2C ∂z 2

→

and grad C =

∂C ∂C ∂C , , ∂x ∂y ∂z

.

EFFECTS OF AIR AND CAVITATION

129

P

r 2A P Po ( ro )3m – r Pv rmax r

Pmin

Figure 8.1: Equilibrium between external pressure and bubble radius.

which allows to compute the maximum dimension of the bubble as: rmax =

1 ( 1−3m ) 2A 3m · p0 · ro3m

(8.6)

To this maximum value of the bubble radius rmax corresponds the minimum pressure: pmin = 2 · A ·

1 ( 1−3m ) 1 3m + pv −1 · · p0 · r03m 3m 2A

(8.7)

As m > 1, pmin < pv . Looking again at Figure 8.1, the stable part of the function (r < rmax ) is related to the gaseous cavitation. If the pressure of the environment p becomes lower than pmin , the bubble becomes instable and tends to expand; as this expansion is quite fast, at least with respect to the diffusion time, the liquid becomes vapour, which remains inside the bubble: this phase is called vaporous cavitation. This exchange between vaporous and liquid phases is much faster than the emission and reabsorption of air. The former phenomenon is quite frequent and cavities can be considered almost always filled with vapour; the latter phenomenon might be neglected. As a consequence, the current practice to verify the presence of cavitation, which consists in comparing the absolute pressure of the liquid with the vapor pressure is fully justified; moreover, if no external air is introduced, the cavities are filled only with vapour which is produced and absorbed during the different phases of the transient.

130 WATER HAMMER SIMULATIONS Finally, if air is introduced into the cavities, obviously their absolute pressure will not be constant, but function of both the volume of the cavity and the mass of air that is introduced (and therefore of the device which allows the entry). 8.1.2

Collapse of the cavities

When the pressure in the pipe increases, the two “stubs” of the water column, upstream and downstream, are quite close, being separated by a cavity which volume tends to zero. Indicating with vupstream and vdownstream the velocities of the front of the upstream and downstream stubs, respectively, in the moments immediately before the collapse is vm > vv and the piezometric head is the same for the two parts and the cavity as well. During the moments immediately following the collapse, there is an unique value not only of the piezometric head but also of the velocity vc . It can be demonstrated that: vm + vv vc = (8.8) 2 and, therefore, the pressure due to the transient can be computed as: H = 8.1.3

c c c v m − vv · (vc − vv ) = · (vm − vc ) = · g g g 2

(8.9)

Description of the motion in the presence of cavities

A problem that might be faced is the evaluation of the position of the cavity: if its position can be assumed known with a good approximation (because it can be supposed positioned close to the regulatory devices, which originate the waterhammer transients), the equations to be written have to allow the determination of the upstream and downstream fronts of the two parts in which the water column is divided. To this end, it is necessary to write an equation that governs the volume of the cavity, function of its pressure. Therefore, it is possible to: 1. suppose that the pressure inside the cavity is equal to the vapour pressure and therefore constant; 2. write the characteristic equation for the “stub” of interest; 3. compute the volume of the cavity W evaluating the differences between the velocities of the upstream and downstream stubs and considering that W = A · (vm − vv ) · t (being A the section of the pipe). It is to be highlighted that in the past the description of the water column separation was performed under the non-elastic hypothesis. However, for this phenomenon more studies, both numerical and experimental, are necessary as the descriptions developed so far are very approximate and there is a need to a better understanding of the physics of the system.


131

8.2 Changing of celerity As mentioned in the introduction to this chapter, the dimensions of gas bubbles have two important consequences related to the celerity of the pressure waves and to the headlosses. With regards to the celerity, qualitatively the significant decrease of the density and the correspondent increase of the fluid elastic bulk modulus can be considered: in other words, a fluid composed of water and vapour bubbles is less dense and more compressible than a fluid composed of water alone. As a consequence, it can be verified that the value of celerity presents impressive decreases even for very low values of vapour concentration. A fast computation can be performed as follows. Let us have, within a unit volume of liquid, a number Nr of bubbles all with the same radius r. The so-called vacuum fraction α is therefore equal to: α=

4 · π · r 3 · Nr 3

(8.10)

The density of the mixture ρm is equal to: ρm = α · ρv + (1 − α) · ρl

(8.11)

where ρv and ρl are the densities of the pure components vaporous and liquid, respectively. The volumetric compressibility of the mixture εm is therefore equal to: dp εm = − dV = V

−dp dα dr

·

dr dp

· dp +

1−α ε

· dp

(8.12)

If, for the sake of simplicity, the hypothesis of isothermal transformation is introduced (i.e. m = 1), eqn (8.12) becomes: εm =

1 1−α ε

+

3 · α · r3 3 · p0 · r03 − 2 · A · r 2

(8.13)

It can be easily demonstrated that εm < ε when r < rmax and has no meaning when r > rmax . If, as usual, the mixture is inside a deformable pipe with diameter d, thickness e and elastic modulus E, the value of celerity to be used for transient analysis is given by: $ % 1 % ρ cm = & 1 m d (8.14) + E·e εm Therefore, cm can be significantly lower than the usual value computed with the formula (2.10) carried out in Chapter 2, even for low values of α. During transients, therefore, significant reductions of the celerity have to be expected when values of α = 0 are expected.


Valve

Flow meter

Pressure transducer

Figure 8.2: Experimental set-up.

8.3 Water column separation As described in Section 8.1, when the pressure in the pipe goes below that of saturated vapour, cavities start to form. When the stress is high, the cavities enlarge to a point which is called water column separation, when the cavity actually separates the two sections of the liquid pipe. For each of the two sections, which are separated by the cavity, the cavity itself is a condition of constant head, equal to the vapour pressure if no air is inserted, artificially or through devices. Supposing a plant as that sketched in Figure 8.2 (which has been installed in laboratory and where tests have been performed, as will be mentioned in the following), where the valve is positioned upstream and the reservoir downstream, a sudden closure of the upstream valve produces an immediate pressure drop; in certain conditions, the pressure drops down to −10.33 m when the water column separates3 . The upstream valve is normally modelled imposing the velocity and retrieving the head from the equation of the system. When cavitation occurs, the boundary conditions change: in the position of the valve a constant pressure is imposed, equal to the mentioned vapour pressure, and the velocity has to be computed from the equation of the system.

As known, the value of −10.33 m corresponds to the absolute pressure equal to zero and therefore the threshold under which the stresses would be of traction; as the liquid do not stand the traction stresses, the water column separates.

3

EFFECTS OF AIR AND CAVITATION Direction of the flow in steady conditions

133

Valve

Air

Water

Figure 8.3: Cavitation recorded during laboratory tests. Concluding, until the water column is separated, the condition to impose upstream is that of a constant reservoir, with head h0 = −10.33 m; when the two sections reattach, the upstream condition again becomes on the velocity. If in the cavity (upstream section) air is artificially inserted, the absolute pressure is not zero but results to be function of (i) the volume of the cavity and (ii) the amount of mass of air inserted. A picture of the phenomenon is shown in Figure 8.3, taken through a transparent pipe positioned downstream the valve. The whole cavitation phenomenon is shown in a short movie enclosed to the book.

8.4 Additional resistance terms 8.4.1 Models It has already been stated that the presence of the bubbles induces additional energy dissipations compared to ordinary hydraulic resistances because of (i) the conversion of mechanical work into heat in the cycles of compression and decompression of the bubbles and (ii) the energy dissipated by viscosity in the liquid surrounding the bubbles. Moreover, some authors (e.g. [78]) believe that the headlosses during transients are higher than those computed in steady flow not only because of the presence of the bubbles but also because of the differences in the velocity profile. These effects are very difficult to quantify. On the other hand, referring to a constant celerity (and equal to that calculated in the absence of bubbles) and to “conventional” headlosses can lead to results too much on the safe side.

134 WATER HAMMER SIMULATIONS However, it is clear that new resistance terms have to be inserted; to evaluate these new terms, many models are available in the literature, which have been classified into six categories [79]. This new term may depend on: • the instantaneous variations of the average flow velocity V ; • the variations of both the average flow velocity V and the local acceleration ∂V /∂t; • the variations of the average flow velocity V , the local acceleration ∂V /∂t and the convective acceleration ∂V /∂x; • the variations of both the average flow velocity V and the diffusive term ∂ 2 V /∂x2 ; • the variations of the average flow velocity V and of the “flow history”, i.e., the record of the past accelerations; • the variations of the instantaneous velocities along the longitudinal and transversal directions. As it is not in the aim of the book to give a complete exposition of all the methods, but only to show how the problem can be dealt, the most common models will be shown and implemented, presenting their critical points also. The models presented in the following are grouped under the name of instantaneous acceleration based (IAB) in which the headlosses are evaluated adding a new term to that related to steady flow, i.e.: J = JS + JU

(8.15)

where JS is the usual steady flow term and JU is the new unsteady formulation. Brunone and Greco [80] evaluated this new term as: ∂Q ∂t

JU = k1 ·

(8.16)

where k1 is an appropriate damping parameter, with the condition

∂Q if sign(Q) = sign ∂t

then k1 = 0

(8.17)

This model has been changed [81, 82] as the need for keeping in account also the convective acceleration has emerged; therefore, the new formulation is: JU = where:

∂V k2 c · · 1− g V ∂t

(8.18)

∂V ∂t V = ∂V

(8.19)

∂x

and k2 is again an appropriate damping parameter, normally different from k1 .


135

Pezzinga [83] commented that the sign minus in eqn (8.18) is not justified and in certain cases it can bring to an unrealistic increase of oscillation, and therefore proposed the following: ∂V k3 · c ∂V k3 ∂V JU = + sgn V · · (8.20) g ∂t ∂x g ∂x Again with k3 a damping parameter is indicated. Starting with this model, Ramos [84] modified again the expression in order to improve the model performances and to match experimental results. All these models can be decomposed in order to separate the two accelerations (local and convective), as follows: kp ∂V ∂V ka · c ∂V JU = + sgn V · · (8.21) g ∂t ∂x g ∂x According to Ramos, the term g1 ∂V does not produce any additional damp ∂t to the signal, but it is responsible only of a change in the value of the celerity, which implies a variation of the frequency of the signal. On the other hand, the term gc · ∂V is responsible only for the signal damping during the transient. There∂x fore, as the physical effects associated to the two terms are different, it is possible to separate them and to use two different parameters, which in eqn (8.21) are indicated as kp and ka and that are generally different. However, in the numerical tests we carried out, the two effects cannot be completely divided, as, even when changing the value of one parameter and keeping constant the other, the carried out solution changes in both the dissipative effects and the variation of velocity; therefore, the performances of this last model were found not to differ too much from those of Pezzinga. 8.4.2

Parameters

The main problem in the use of these models is the choice of the appropriate value for the damping parameter(s). To this end, in the scientific literature many sets of values that reflects the different conditions are reported. For instance, Brunone et al. [81, 82] suggest a single value k1 = 0.085; Bergant and Simpson [79] found values within the range [0.085, 0.330]; Bughazem and Anderson [85] within the range [0.065, 0.150] and Wylie [86] within the range [0.02, 0.03]. The usual suggestion, however, is to calibrate this parameter using experimental data, as no reliable analytical or graphical procedures to estimate a parameter are available yet. Comparing the results of the above-mentioned models it has been demonstrated that: k1 ≈ 100 · k2 ≈ 100 · k3

(8.22)

Expression (8.22) tries to tie the parameters of the mentioned formulations; in the following some empirical estimations of the required parameters is presented, with

136 WATER HAMMER SIMULATIONS the note that the values that can be carried out with these procedure are strongly approximated and should be verified. Pezzinga [87] tried to build an abacus tying the parameter k to the plant characteristics. As Vardy and Brown [88] have shown, the value of the parameter k decreases when the Reynolds number Re increases; however, a number of different parameters have to be introduced to properly describe the behaviour of k, which are ε/D (the relative roughness of the pipe, where ε is the roughness and D is the pipe diameter) and especially: y0 =

g · J0 · L c · V0

(8.23)

where J0 and V0 are the head grade line and the velocity in steady flow, respectively. Varying these non-dimensional parameters, a number of abacus have been drawn, similar to that presented in Figure 8.4, shown as an example. With only one experimental test performed, Carravetta et al. [89] claims there is the possibility to give an appropriate estimation of the required parameter. The base of this method is the consideration that [90] the term related to steady flow headlosses JS is negligible when compared with JU , at least in the cases of complete closure that ends in few phases. Therefore, Carravetta et al. [89] tied the damping parameter k1 to the peak values recorded in two consecutive periods ht and ht+1 as follows: 2 ht+1 1 = (8.24) ht 1 + k1 k2

ε/D

0.05 0.03 0.015 0.008 0.003 0.001 0.0003 0.00005 0.00001

Re0 102

103

104

105

106

107

108

Figure 8.4: Abacus to estimate the damping parameter k (after Ref. [87]).


137

8.5 Laboratory experiments 8.5.1 Experimental set-up In order to verify the effects of the sudden closure of a valve that creates depression in a closed pipe, a very simple experimental plant has been built in the Hydraulic Laboratory of the Politecnico di Milano, as sketched in Figure 8.2: it is a single pipe that links two reservoirs, the elevation of the upstream reservoir being equal to 17.7 m, while the elevation of the downstream reservoir is equal to 5 m. A pipe which diameter is 300 mm is connected to the upstream reservoir; the diameter of this pipe becomes 150 mm and the final part of the plant is constituted of an iron pipe with diameter D = 52 mm and thickness e = 5.0 mm; this latter pipe is wound in a roll of 26 coils for a total length of 90 m, and then connected to the downstream reservoir. The transient is recorded only in this latter pipe: the flow is interrupted by means of a very fast electro-valve positioned on this pipe and pressures are recorded through a transducer positioned immediately downstream the valve. The discharge is measured with an electromagnetic device. Tests carried out in steady flow showed a Strickler ks resistance parameter equal to 110 m1/3 /s, i.e., Manning equal to 1/110 = 0.009091 m−1/3 × s. The value of the celerity can be computed with the usual formula (2.10) given the pipe characteristics and known that the iron elasticity modulus is E = 2.0 × 1011 N/m2 , water compressibility ε = 2.14 × 109 N/m2 and water density ρ = 1000 kg/m3 . It resulted c = 1388 m/s. 8.5.2

Experimental tests

During the years many tests have been performed, and it is obviously not possible to present all of them here. Among the many tests, five have been selected and discussed in this chapter, which are useful for a further discussion. The characteristics of the five tests are shown in Table 8.1. Moreover, in Figures 8.5–8.9 the recorded pressures are shown for each test. Table 8.1: Characteristics of the experimental tests discussed in the book. Test no. (–) 1 2 3 4 5

Discharge (l/s) 0.112 0.228 2.5 3.0 3.3

Velocity (m/s) 0.0527 0.1074 1.18 1.41 1.55

Head U/S (m) 5 5 8.1 9.7 10.5


15

Head (m)

10

5

0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

2.5

3.0

3.5

4.0

5 10 Time (s)

Figure 8.5: Test #1. 30 25 20

Head (m)

15 10 5 0 0.0

0.5

1.0

1.5

2.0

5 10 Time (s)

Figure 8.6: Test #2. As can be seen, tests 1 and 2 do not reach the value of zero absolute pressure (−10.33 m), but the pressure drops under the atmospheric, and therefore bubbles are expected. Tests 3–5, instead, are characterized by water column separation (cavitation): in this case the increase of the velocity cannot force an increase of the depressions as obviously absolute pressure cannot fall below zero; however, there is an increment of the time during which the water column separation occurs.


139

80 70 60

Head (m)

50 40 30 20 10 0 0.0 10

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

2.5

3.0

3.5

4.0

20 Time (s)

Figure 8.7: Test #3. 100

80

Head (m)

60

40

20

0 0.0

0.5

1.0

1.5

2.0

20 Time (s)

Figure 8.8: Test #4. When the water column reattaches, a large turbulence is recorded and then the oscillations around the mean value start. In test #1, the temporal distance between the first two pressure peaks is equal to 0.15 s. It is therefore possible to compute the celerity: τ=

2·L c

⇒

c=

2·L 2 × 90 ∼ = = 1200 m/s τ 0.15

(8.25)


80

Head (m)

60

40

20

0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

20 Time (s)

Figure 8.9: Test #5. As can be seen, the experimental value of the celerity recorded in test #1 is much lower than that theoretically computed above, as expected because of the presence of air bubbles.

8.6 Computer code A computer code has been developed with different options in order to check whether it is possible to simulate these cases with the methods described so far. The computer program has been developed starting from the code that uses the Lax–Wendroff method for the internal points and the Ghost-Cell Zero-Order for the boundary conditions. It implements the Pezzinga’s method [83] as presented in Section 8.4.2 and in the following better specified. Moreover, the possibility of evaluating the cavitation is given. For the sake of brevity, in this book not all the passage are reported to implement the method of Pezzinga with the Lax–Wendroff ’s integration scheme, but only the final formulae: g λ k3 · c n n Vjn+1 = Vjn − · · (hnj+1 − hnj−1 ) ± · (Vj+1 − Vj−1 ) 2 1 + k3 1 + k3 2 V n · g λ k3 · c · g j · (hnj+1 − 2hnj + hnj−1 ) + ± · 1 + k3 2 (1 + k3 )2 k3 · c 2 g c2 n n n · (Vj+1 − 2Vj + Vj−1 ) − + · t · Jjn + 1 + k3 1 + k3 1 + k3 (8.26)


141

Figure 8.10: Users’ interface for the program to simulate unsteady flow headlosses and cavitation.

hn+1 j

λ c2 n n n n n = − · Vj · (hj+1 − hj−1 ) + · (Vj+1 − Vj−1 ) 2 g c2 λ2 n 2 · (hnj+1 − 2hnj + hnj−1 ) + · (Vj ) + 2 1 + k3 Vjn · c2 k3 · c 3 n n n + ± · (Vj+1 − 2Vj + Vj−1 ) g g · (1 + k3 ) hnj

(8.27)

As can be seen in Figure 8.10, which reports the users’interface, in the computer program there is the possibility to insert two coefficients k: the former, named Pezzinga’s coefficient, is used when cavitation does not occur; if there is cavitation, then a different coefficient is used. In the latter case, however, it is not possible to estimate this parameter through relations as those of Carravetta and Vardy and Brown (see Section 8.4) as the flow characteristics are quite different. To reproduce test #1, the first method which has been applied is the standard one, with no rupture of the water column and no added resistance terms. Then, the method of Pezzinga is used, letting k = 1.5. The results from the application of these two methods are compared with the results of test #1 and shown in Figure 8.11. As can be seen, in the classical model, as presented in Chapters 4 and 5, the oscillations are not appropriately damped by the steady flow resistance term; moreover, there is a slight difference in the celerity of the waves. With the added resistance

142 WATER HAMMER SIMULATIONS 20 Experimental data “Classical” model Resistance terms, k 0.1

15

Head (m)

10

5

0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

5 10 Time (s)

Figure 8.11: Test #1 compared with “classical” model (no added resistance terms) and with the model where the Pezzinga’s resistance term is added.

term, the improvement is highly visible, and the simulated pressures are very close to those recorded. However, it can be observed as the presence of the resistance term is also influencing the celerity of the wave, as it has been already discussed. It should also be noted that the maxima values recorded from the experiments are higher than the theoretical values. This phenomenon is probably due to the turbulence, as will be discussed in Chapter 9. Obviously, the simulation of the experiments where cavitation occurs is more difficult; to do that with the computer program, it is necessary to tick the option Water Column rupture in the users’interface and to choose an appropriate parameter for the model to be used in case of cavitation. Different parameters have been used for the simulations as reported in Figure 8.12. The first simulation has been performed using the same values for the two parameters and equal to that used for the above reported simulation of test #1 (kPezzinga = kCavitation = 0.1). As can be seen, the time when the column reattaches is perfectly estimated, but unfortunately the peak is highly overestimated. Moreover, while in the real data the headlosses due to the turbulence that follows the water column reattachment are enough to avoid subsequent separations, the computer model does not dissipate sufficient energy and therefore behaves as if it was a series of ruptures. As mentioned above, the increase of the value of k brings to a larger damping of the oscillations but also influences the celerity. In the second simulation the values kPezzinga = 0.1 and kCavitation = 1.5 have been used. Here, the time when the water column reattaches is delayed; however, the value of the peak is decreased and,

EFFECTS OF AIR AND CAVITATION 140

143

Experimental data kPezzinga kcavitation = 0.1 kPezzinga 0.1; kcavitation 1.5 kPezzinga kcavitation 1.5

120 100

Head (m)

80 60 40 20 0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

20 Time (s)

Figure 8.12: Test #3 compared with the model where the Pezzinga’s resistance term is added, for different parameters. therefore, more similar to that recorded. In any case, however, the energy dissipations that happens when the water column reattaches are much higher than those simulated by the model. The need to differentiate the damping parameter k in case of cavitation is clear observing the third simulations in the same figure: the results carried out leaving kPezzinga = kCavitation = 1.5 are too damped towards the end of the simulation, i.e., when the cavitation phenomena are concluded. However, even if “playing” with the parameters the results may be improved (and the reader is invited in trying), it is quite clear that this topic is still very uncertain and researchers should investigate it more deeply.


Chapter 9 Advanced models

As it has been seen in Chapter 8, classic models fail when applied to complex cases. To this end, more complex models have been developed, even if the research in this field is still open and the efforts of the researchers working all around the world are yet to be considered concluded. In this chapter two-dimensional (2D) models are shown, and the need to keep in due account the effects of turbulence are also shown, which demonstrates in some cases there might be the need to move towards three-dimensional (3D) models. Obviously, these models have drawbacks which will be described briefly.

9.1 2D models The evaluation of the headlosses during transients can be more adequately performed through 2D or quasi-2D models, even if the computational times increase. In these models it is possible to keep in due account the non-uniformity of the velocity distribution in the transversal section. Bratland [91] proposed the so-called cylinder model for non-laminar flow, in which the continuity equation is solved in one dimension and the momentum equation in two dimension. This formulation has been extended to the turbulence flow by Modica and Pezzinga [92] and Pezzinga [93], who allowed to clarify the nature of the simplifications. The reference system in two dimensions is shown in Figure 9.1; the equations that describe the phenomenon have to be rewritten using this system.

r x

Figure 9.1: Two-dimensional reference system.

146 WATER HAMMER SIMULATIONS 9.1.1 Continuity equation The axially symmetric two-dimensional field in an elastic circular pipe can be described by means of a quasi-two-dimensional model based on the continuity equation written as [94]: ∂ρ ∂(ρu) 1 ∂(ρrv) + + · =0 ∂t ∂x r ∂x

(9.1)

where x is the distance along the pipe, r is the distance from the central axis and u and v are the velocities along x and r, respectively. Analysing the numerical results carried out by Vardy and Hwang [94], it resulted that the velocity v along the radius and its derivative take very low values both in laminar and turbulent flow, and that justifies the choice of Bratland [91] to write the continuity equation neglecting the term 1r · ∂(ρrv) ; eqn (9.1), therefore, becomes: ∂x ∂ρ ∂u +ρ =0 ∂t ∂x

(9.2)

This equation can be integrated along the transversal section, neglecting the derivative of ρ with respect to x and, therefore, obtaining: ∂H c2 ∂Q + =0 ∂t gA ∂x 9.1.2

(9.3)

Momentum equation

The 2D momentum equations can be written as [93]: ∂u ∂u ∂u ∂H 1 ∂σx 1 ∂(rτ ) +u +v = −g − − ∂t ∂x ∂r ∂x ρ ∂x ρr ∂r ∂v ∂v ∂u ∂H 1 ∂τ 1 ∂(rσr ) +u +v = −g − − ∂t ∂x ∂r ∂r ρ ∂x ρr ∂r

(9.4a) (9.4b)

where σx and σr are the longitudinal and radial components, respectively, of the normal stresses (pressures), and τ is the shear stress. Also in this case it is possible to introduce some simplifications in the system. In particular, Pezzinga [93] posed the following hypotheses: 1. As described in Section 9.1.1, the transversal velocity v can be neglected together with all its derivatives; 2. In longitudinal direction the derivative of the normal stress σx is also neglected, together with the residual convective term; 3. For each time step a single value of the piezometric head is assumed in each section, which means the variations of the stresses in the radial direction is negligible, i.e.: ∂H = 0. ∂r

ADVANCED MODELS

147

The momentum equation can, therefore, be written as: ∂u ∂H 1 ∂(rτ ) +g + =0 ∂t ∂x ρr ∂r

(9.5)

Finally, considering that dA = 2π · r dr, eqn (9.6) becomes: ∂u ∂H 2π ∂(rτ ) +g + =0 ∂t ∂x ρ ∂A

(9.6)

The model is called quasi 2D because the dependant variable u is function of x, r and t, while H is a function of x and t.

9.2 Headlossess There is the need to develop a model for evaluating the shear stress τ inside the pipe. In this section two families of model are presented. The former is called “two-layer turbulence model”: in the generic cross section two layers are identified, with different hydrodynamic characteristics, and for each of them a flux law is developed. Closest to the pipe wall there is a delimited field where viscous stresses are higher than those turbulents and, therefore, the flow is laminar. Outside that zone, where turbulent stresses are prevailing, the flow is obviously turbulent. In the scientific literature, a number of methods have been developed, which divide this area into three or even five subareas, where different flux laws are defined, in order to get as close as possible to the actual conditions. The latter family of model presented is called “low Reynolds number k–ε model” that allows the further generalization of the shear stress τ evaluation but requires higher computational times. 9.2.1

Pezzinga’s model

Referring to the works of Bratland [91], Vardy and Hwang [94] and then of Pezzinga [93] who considered turbulent flow, the shear stress τ can be evaluated through a model based on the Newton’s law inside the laminar layer; outside this layer, the hypothesis of the Prandtl’s mixing length is used: ∂u ∂u ∂u τ = −ρ · υ · − ρ · l 2 · · (9.7) ∂r ∂r ∂r where l is the mixing length and υ is the kinematic viscosity. To evaluate the length l the following expression has been used [49, 95]: l y = κ · · e(−y/R) R R

(9.8)

where y is the distance from the pipe wall, R is the radius of the circular pipe and κ is a parameter dependent on the steady flow Reynolds number Re0 .

148 WATER HAMMER SIMULATIONS 0.42

k R/ε 507 R/ε 252

0.41

R/ε 126 R/ε 60

0.40

R/ε 30.6

0.39

R/ε 15 smooth

0.38

Eq. (9.9)

0.37 R

0.36 103

104

105

106

107

108

Figure 9.2: Parameter κ function of Re0 (after [93]). Pezzinga [93] computes the parameter κ interpolating the experimental results carried out by Nikuradse [57, 96, 97], as shown in Figure 9.2 and, therefore, obtaining the equation: 83,100 κ = 0.374 + 0.0132 · ln 1 + Re0

(9.9)

In order to compute the thickness of the laminar sublayer δ, it is convenient to make the hypothesis that inside this layer the velocity is linear and outside it the velocity profile follows a logarithmic law. Let the speed friction u∗ tied to the shear stress at the pipe walls τw from the: τp ∗ u = (9.10) ν Inside the laminar sublayer the linear equation to describe the velocity is: u∗ u∗ y = υ ν

(9.11)

The logarithmic profile of velocity outside the sublayer is defined by: u∗ δ = 2.5 · ln +B υ ε

(9.12)

ε being the roughness and B is a parameter dependent on the parameter: Re∗ =

u·ε ν

(9.13)

ADVANCED MODELS

149

The thickness of the sublayer is computed requiring that the two profiles match, i.e.: u∗ δ δ = 2.5 · ln +B υ ε

(9.14)

The parmeter B is dependent on Re∗ by the expression: 3.32 B = 8.5 − 2.5 · ln 1 + Re∗

(9.15)

9.2.2 k–ε model The model presented in the Section 9.2.1 can be generalized with the so-called low Reynolds number k–ε model. These families of models have been developed [98] to describe the turbulence and integrate the Navier–Stokes equations simplified to the Reynolds equations. Eichinger and Lein [99] applied this method to describe the pressure fluctuations recorded during transients; this procedure has been lately used, with modifications, by Pezzinga and Brunone [100] and Zhao and Ghidaoui [101]. This method can be applied in a quite wide range of Reynolds numbers, from laminar to purely turbulent flow. The shear stress τ is evaluated as: τ = −ρ · (υ + υt )

∂u ∂r

(9.16)

where υt is the value of turbulent viscosity, computed with: υt = fu · cu ·

k2 ε

(9.17)

k being the kinetic energy for mass unit and ε is a parameter for the dissipation effects. These variables are governed by the following transport relations [102]: ∂k 1 ∂ − · r· υ+ ∂t r ∂r ∂ε 1 ∂ · r· υ+ − r ∂r ∂t

2 ∂k ∂u +ε =0 − υr · ∂r ∂r 2 ∂ε υt ε2 ∂u · − ce1 · f1 · υr · + cε1 f2 = 0 σk ∂r ∂r k υt σk

·

(9.18) (9.19)

Pezzinga and Brunone [100] showed that in case of transients, the formulation that produce best results is that of Lam and Bremhorst [103]. The parameters to be applied with this method have been evaluated through semi-empirical formulations and they are: 20.5 fu = [1 − exp(−0.0165 · y∗ )] · 1 + (9.20) Rt

150 WATER HAMMER SIMULATIONS f1 = 1 +

0.05 fu

3

f2 = 1 − exp(−R2t )

(9.21) (9.22)

√ where: Rt = υ·ε and y∗ = υk·y . Boundary conditions for k and ε are: k = 0 and 2 ε = ν · ∂∂rk2 at the pipe walls; ∂k = 0 and ∂ε = 0 on the central axis. Numerical ∂r ∂r parameters have the following values [103]: cu = 0.09; σk = 1; σε = 1; σε1 = 1.44; K2

σε2 = 1.92. Finally, the initial values of k and ε are obtained imposing the steady flow conditions.

9.3 Cavitation In Chapter 8 the possibility to have vaporous or gaseous cavitation has been mentioned, reminding that, as water is normally saturated with air, mainly the latter case happens in the real plants; moreover, Wylie [104] showed that gaseous cavitation models are able to reproduce quite well the vaporous cavitation phenomenon. As a consequence, in the following only that kind of model is described. The presence of air, as shown in Chapter 8, increase the headlosses, already larger than those expected using the steady flow formulas. In the scientific literature two families of models can be found: those that analyse the effects of the free gas in the fluid (“release gaseous cavitation model”) and those that analyse the thermic exchanges between the free gaseous and the liquid phases (“second viscosity cavitation model”). 9.3.1 Release gaseous cavitation model Under this type fall those models that consider constant the mass of gas [11] and keep in account the whole gas release process [105]. The density of the mixture can be evaluated as: mRT mRT ρ = 1− · ρW + · ρg p p

(9.23)

where m is the mass for specific volume; R is the constant of the gases; T is the temperature and p is the pressure. Considering the density of the mixture as function of the pressure and of the mass of the gas free, the continuity equation in an elastic pipe can be written as [105]: ∂p ∂m ρc2 mRT ρc2 ∂Q ρc2 RT + 1+ 0 2 · c02 − 0 + 0 · =0 ∂t p ∂t p A ∂x

(9.24)

where c0 is the celerity for an homogeneous liquid. To apply a consistent numerical scheme, the equation is to be written in conservative form. Introducing the variable

ADVANCED MODELS

151

, defined as: =

c2 mRT c2 m p − 0 + 0 ρg gp ρg

(9.25)

The continuity equation becomes: ∂ ρc02 ∂Q + · =0 ∂t A ∂x

(9.26)

Neglecting the convective term, the variation in time of the air mass can be computed as constant, as mentioned, but the carried out results show that this model is not able to correctly reproduce the experimental data. On the other hand, Zielke et al. [106] introduced a model where the air mass is variable, proving the hypothesis that there is an initial quantity of mass m0 and the water is air saturated; then the gas is released following the law: ∂m 1 β = · · (ps − p) ∂t θ RT

(9.27)

where β is the solubility constant of the gas, θ is the relaxation time, ps the saturated pressure and p the pressure of the gas. Generally speaking, the models that consider the air mass as constant are able to represent the main characteristics of the phenomenon and in particular the effects of propagation of the increased deformability of the mixture; the models that try to use a release law are potentially able to describe the energy dissipations not due to the friction. Obtaining the value of the free mass from eqn (9.27) and inserting it in eqn (9.25) it is possible to get the value of the pressure p:

p=

+

' + 4ρc02 mRT 2

(9.28)

Finally, the piezometric head is given by: H =z+

p + p v − pa ρ·g

(9.29)

where pv is the vapour pressure and pa is the atmospheric pressure. It is to be highlighted that to the gas pressure is added also the vapour pressure, as both are present inside the bubbles, even if the vapour mass might be considered negligible [107]. The problem is the determination of the initial mass m0 and of the relaxation time θ. The former parameter influences the oscillation phases, while the latter their dampening. The calibration of these parameters is very difficult and analytical solutions that allow an empirical evaluation do not exist yet. So far, only estimations based on experimental data are provided [105, 108].

152 WATER HAMMER SIMULATIONS Finally, the momentum equation can be written under the hypothesis of a twodimensional field with axial symmetry in a circular pipe [93]: ∂H 1 ∂(rτ ) ∂u +g + =0 ∂t ∂x ρr ∂r

(9.30)

Again, it is necessary to introduce a model that allows the definition of the shear stress τ ; the models described in Section 9.2 can be used. 9.3.2

Second viscosity cavitation model

Some authors ascribed non-negligible effects to dissipations due to the thermic exchanges between the free gaseous and liquid phases [109]. In order to evaluate these effects, in the momentum equation a new term is added, called pseudovelocity or volumetric velocity ξ . This term allows to merge the mechanical and thermodynamic effects in only one expression, which can be easily inserted in the momentum equation. In the following model [110], the free gas for unit volume has been considered as constant [11] for the sake of simplicity. The continuity and momentum equations, written under the usual hypotheses of a velocity profile in axial symmetry are: ∂p ρc2 ∂Q + =0 ∂t A0 ∂x

(9.31)

∂u ∂H 1 ∂(rτ ) 1 ∂σx +g + − =0 ∂t ∂x ρr ∂r ρ ∂x

(9.32)

In this case only the term dependent on the volumetric viscosity ξ is added, thus obtaining [110]: ∂H 1 ∂(rτ ) 1 ∂ ∂u ∂u +g + − ξ =0 ∂t ∂x ρr ∂r ρ ∂x ∂x

(9.33)

Referring to [110] for details, the involved parameters can be computed as: 2 ρ · m · R · T · c4 ρ 2 · m · R · T · c4 k −1 ·θ · =· ξ= θ+ 2 k p p2 ⎛ ⎞ ⎜ c = c0 · ⎝ '

p p2 + ρ · c02 · m · R · T

⎟ ⎠

(9.34)

(9.35)

where c0 is the celerity evaluated when m = 0, ρ is the density of the mixture water-vapour, which can be computed with eqn (9.23), θ is the relaxation time, k is a coefficient dependent on the gas transformation, m is the mass of the gas

ADVANCED MODELS

153

for unitary volume, R is the constant of the gases, T is the temperature, p is the pressure and is the relaxation time of the equivalent gas. The governing equations are therefore the continuity equation, as expressed in Section 9.1.1, and the momentum equation which finally is: ∂u ∂H 1 ∂(r · τ ) ∂ ρ · (c0∗ )2 · m · R · T ∂p +g· + +· =0 · ∂t ∂x ρ · r ∂r ∂x p2 + ρ · (c02 )2 · m · R · T ∂t (9.36) With a procedure similar to that used for the continuity equation, the variable can be defined as: ' p = ρ · (c0∗ )2 · m · R · T · arctg " (9.37) ρ · (c0∗ )2 · m · R · T Moreover, imposing the pipe is circular and therefore dA = 2π · r · dr and the final equation is: ∂u ∂H 2π ∂(r · τ ) ∂ 2 +g· + · + · =0 ∂t ∂x ρ ∂A ρ ∂x · ∂t

(9.38)

Again, the parameters needed to be used in the model should be calibrated on experimental data, as no empirical relations are available yet. It is to be reminded that in this expression the hypothesis also of constant gas/vapour mass is made, while Cannizzaro and Pezzinga [111] extended the model to be used also when this mass is variable, and therefore improving the results.

9.4 Numerical schemes The integration of the equations is carried out using a cylindrical grid where width is constant in the longitudinal direction, while the area is constant in the radial direction, as shown in Figure 9.1. Velocities are defined on the median circumference, while the shear stresses are defined on the internal and external circumferences. The celerity of the mixture in a deformable pipe is computed by: c= '

c0 1+

ρ·c02 ·m·R·T p2

(9.39)

The integration of the equations can be carried out by means of the method of characteristics [94] although the dependence of the celerity on the pressure would suggest more appropriate and different methods. As in the momentum equation for a given grid the velocities in the neighbouring grids are present, for each transversal section a non-linear equation system has to be solved. This implies the methods like the Lax–Wendroff cannot be used, as they are valid for quasi-linear systems. Hirsch [31] shows that it is still possible to extend the Lax–Wendroff scheme to non-linear systems too, if the final equation is considerably different from that presented in

154 WATER HAMMER SIMULATIONS the preceding chapters. Alternatively, to solve the problem the McCormack method [76] has been successfully used; this method is second order accurate both in space and time; the matrix to be solved is tri-diagonal and dependent on the unknowns vector given by the velocities. The solution is obtained through an iterative method, coefficients are computed at each iteration and the first trial values are let equal to those computed at the preceding temporal step [112]. Generally speaking, the results carried out are significantly better than those obtained with one dimensional models; unfortunately, if a one-dimensional model takes few seconds to produce a solution, two-dimensional ones may take, for the same problem, several hours. The time required to give a solution can even increase if complex turbulence models, like k–ε, are applied. An idea for further research could be the development of one-dimensional code able to shift to a two-dimensional when and where it is needed. This involves the determination of a threshold of validity of one-dimensional models, or acceptability of its results, which so far is not yet quantified in the scientific literature.

9.5 Further problems A simple experimental plant has been set up in the Hydraulic Laboratory at Politecnico di Milano being very similar to that presented in Chapter 8, the valve to interrupt the flow is positioned downstream the plant and the pressures are recorded immediately upstream the valve with a very sensitive pressure transducer, able to sample at 1000 Hz. The recorded pressures are shown in Figure 9.3: as can be seen, 4.0 Recorded data 3.5

Theoretical value

Pressure (bar)

3.0 2.5 2.0 1.5 1.0 τ0

0.5 0.0 4500

4600

4700

4800

4900

5000

5100

5200

5300

5400

5500

Time (ms)

Figure 9.3: Experimental data carried out with the closing valve positioned at 0.10 m from the pressure transducer.

ADVANCED MODELS

155

3.5 Recorded data 3.0

Theoretical value

Pressure (bar)

2.5 2.0 1.5 1.0 0.5 0.0 17.1

17.2

17.3

17.4

17.5

17.6

17.7

17.8

17.9

18.0

18.1

Time (s)

Figure 9.4: Experimental data carried out with the closing valve positioned at 1.50 m from the pressure transducer. the peaks reach higher values than the expected, theoretical, computed with the formulas presented in Chapter 2. This is unexpected and unjustifiable considering even different profile of velocities that could be reproduced using the models described in this chapter. The high turbulence recorded is attributed to its proximity to the valve, which in this configuration was equal to 0.10 m; therefore the distance between the valve and the transducer has been increased to 1.5 m and the test repeated. Values are slightly different because the flowing discharge (which has been measured with an electromagnetic device) was not perfectly the same in the two tests. However, when the valve is moved, the recorded data show an almost perfect agreement with the expected theoretical value and the turbulence practically disappeared (Figure 9.4). That means, on one hand, that to reproduce the effects of turbulence, the models that have been shown in this book, although complex, are still not appropriate, and probably the interested user should try the most complex fully three-dimensional codes, accepting to face all the related problem – for instance, the choice of the appropriate turbulence model, the complexity in building the model, the time required for simulation that would practically prevent its use for complex networks. On the other hand, the user should be aware of these effects and these peaks, but also consider their limited extension in time and space, and therefore whether they are significant for the plant to be designed or verified.


Chapter 10 Case studies

All the theories and formulas presented in this book can be applied to real cases; in the following, studies of existing plants are presented in order to give an idea of the approach that can be used. Three cases are presented with increasing complexity.

10.1 Simple pressure pipe for petroleum products in Djibouti In this case the object of the study was the evaluation of the unsteady flow pressures that can be found in a terminal pipe conveying petroleum products in a plant located in Oriental Africa. The study has been performed with an elastic model integrated using the method of characteristics. 10.1.1 Plant characteristics The plant to be modelled consisted in a simple pipe, with the following boundary conditions: upstream one or two pumps operating in parallel and downstream a motorised valve. System characteristics are as follow. Pipe diameter is D = 16 = 400 mm (but in the last 20 m where D = 12 = 300 mm); total length is equal to L = 2200 m; pipe thickness is equal to e = 7.92 mm. Pipe is built in carbon steel, with elasticity modulus equal to E = 2.0 × 1011 N/m2 and its junctures are welded. Transported fluid is hydro carbonic combustible derived by petroleum, with density equal to ρ = 850 kg/m3 and compressibility coefficient equal to ε = 1.5 × 109 N/m2 . In steady flow, discharge Q and velocity V in the pipe are, respectively, equal to Q = 750 m3 /h ⇒ V ∼ = 1.66 m/s (where D = 16 ) when only one pump is working. When two pumps are working these values can be simply doubled1 , obtaining Q = 1500 m3 /h ⇒ V ∼ = 3.32 m/s. 1

Note that this simplification is not strictly correct if the pumps are working with fixed velocity; however, it was discussed with the client, considered on the safe side and finally accepted, especially because the headlosses were not known.

158 WATER HAMMER SIMULATIONS As in the last end of the pipe the section is smaller, outfall velocity is equal to V = 2.95 m/s (when only one pump is working) and V = 5.89 m/s (when the two pumps are working together). Pressure downstream the pump is p = 1.1 × 106 Pa. 10.1.2

Expected scenarios

Expected movements of the devices in the plant that have been simulated are: 1. the downstream motorised valve gets closed and upstream pumps are still working (they will be switched off at the end of the transient); 2. pumps are switched off and then the motorised valve is closed (case of pumps failure); 3. pumps are switched on and then the motorised valve is open. Due to the lack of more precise information on the valve closure, the linear closure hypothesis is adopted, using the law mentioned in Chapter 4. This hypothesis, however, is quite cautious and thus can be adopted. 10.1.3

Case 1

The first scenario is with the upstream pumps working and the downstream valve getting close. This is the normal procedure adopted to switch off the pumps. The static head (around which the unsteady flow pressures oscillate) is equal to the head of the pump when the discharge is zero. From pump curves it is deduced that this value is equal to 12 atm. Headlosses in the pipe have been computed so that at the outfall the pressure is equal to 1 atm, i.e., it is assumed that the upstream pressure (of the pumps) is necessary and sufficient to convey the discharge Q = 1500 m3 /h. In Figure 10.1 the results of this simulation have been reported; it can be seen that the unsteady flow pressures oscillate, as expected, around the static head equal to 12 atm. On the other hand, the maximum pressure reached is lower than that computed with the Michaud’s formula in Chapter 2, which brought 15.45 atm; instead, the maximum pressure value reached with the computer program is equal to 14.34 atm. This result is ascribed to the headlosses in the downstream end of the plant, where the pipe has smaller diametre. 10.1.4

Case 2

The second scenario described in the Section 10.1.2 is when the upstream pumps are turned off while the downstream valve is still open. The computer program is built modelling the upstream pumps as if they were a valve that closes in fixed times, selected in Tc = 15 and 30 s, which is the range of time in which the discharge is expected to stop. The model was simplified because no precise information

CASE STUDIES

159

16 14

Pressure (atm)

12 10 8 6 4 2 0 0

10

20

30

40

50

60

Time (s)

Figure 10.1: Results of simulations for case 1. were provided. Downstream the liquid outfalls in the atmosphere with concentrated headloss equal to α · V 2 /2 · g; the α coefficient has been calibrated in order to have atmospheric pressure in the downstream section of the plant, in steady flow conditions. The value of the coefficient α carried out in steady flow conditions has also been used in unsteady flow simulations. Simulations show that in both cases (Tc = 15 and 30 s) pressures fall below absolute zero. As largely discussed these are non-physical values: after reaching the points shown in Figure 10.2 indicated as “Vacuum”, the numerical model does not produce realistic results. It has to be noted that this case should be avoided and that the correct sequence for turning off the pumps is that described in case 1. This second scenario happens only in case of pumps failure, and therefore not very frequently. The curves in Figure 10.2 are related to the case of both pumps failure (power failure); however, simulations also showed that when only one pump is working, the pump failure brings the system to cavitation, even if of less important magnitude. 10.1.5

Case 3

The third scenario described in Section 10.1.2 is the opening of the downstream valve while the upstream pumps are already working. This is the correct starting procedure. The downstream valve is modelled with linear opening law, considering it is completely open in 36 s. In Figure 10.3 the results are shown; it can simply be observed that the pressures decrease, never reaching vacuum, and then there are oscillations around the steady flow value.

160 WATER HAMMER SIMULATIONS 14 12

Closure: 15 s

10

Closure: 30 s

8 Pressure (atm)

Vacuum 6 4 2 0 2

0

10

20

30

40

50

60

50

60

4 Pressures below 1 atm are not possible. Left as model output.

6 8

Time (s)

Figure 10.2: Results of simulations for case 2.

14 12

Pressure (atm)

10 8 6 4 2 0 0

10

20

30

40

2 Time (s)

Figure 10.3: Results of simulations for case 3.

CASE STUDIES

161

10.1.6 Conclusions The problem of transients in a very simple plant has been studied by means of the elastic model equations integrated with the method of characteristics. Different scenarios have been drawn and then evaluations have been carried out. Results are quite consistent with those obtained with empirical methods (slight differences are noted); carried out results show pressures are within the range of acceptability for the construction plant standard, if the computational hypothesis assumed can be considered valid and if the correct sequences of starting and stopping of the plant are applied. A problem of cavitation might arise in case of pumps failure, which is obviously not the correct sequence of stopping the system and therefore that happens quite rarely. However, according to the client, this outcome can be also accepted. In normal operations, carried out results show pressures are within the acceptable limits and therefore the plant could be considered safe.

10.2 A more complex example for seawater treatment plant in Tanzania For the design of a seawater treatment plant in Tanzania, we have been asked to verify the steady and unsteady flow conditions of the pumping station both in existing and design conditions. Moreover, we have been asked to determine whether any remedial is required for the pumping station and, in this case, to compute its type and dimensions. In this paragraph, after describing the sketch of the two systems as they have been implemented in the computing code, some results of unsteady flow simulations are presented. The mathematical model used for unsteady flow simulations is based on the method of characteristics.

10.2.1 Plant characteristics The two plants (existing and to be designed) are shown in Figures 10.4 and 10.5 with their computational model. As can be seen, some pipes and junctions had to be inserted into the model in order to keep in account particular devices as the air chambers. Pumps characteristics are reported in Table 10.1, while pipes are described in Table 10.2 and nodes in Table 10.3. In Table 10.2 celerity has been computed with the usual formula (2.10) where the parameter λ is equal to 1. Finally, note that the existing air chamber volume is equal to 6.37 m3 : but this is the total volume of the chamber, while the volume of the air inside the chamber


P1

J1

P2

AC1 P3

01 02

LEGEND: Tank

Junction

Pump

Air chamber 8

1

4

10

18

12 15

19 3

14

6

17

7

9

20

21 11

Figure 10.4: General sktetch and numerical implementation – existing plant.

is smaller. As we do not know that volume, we suppose it is about he half, so equal to 3.2 m3 .

10.2.2

Unsteady flow simulations: existing plant

In the case of the existing plant, it can be easily seen that the times in which the pumps stop (approximately equal to 20 s) is much higher than the phase time of the system, so that unsteady flow pressures are quite low. In Figure 10.6 is shown the static head (equal to the elevation in the downstream tank) and the pressures in the upstream end of pipe 15 (see Figure 10.4, lowerpart). At the beginning (time t = 0) the head is higher than the static value; then as the pumps stop, there is a sudden drop in the pressure; afterwards, the head oscillates around the static value, as expected. As can be seen, oscillations are very small. Other pipes have similar values. During the simulations the clapets

CASE STUDIES

163

P4

J2 P5

AC2

P6

J3

Line J1(a)–J3(a)

P4 P5

03

J1

AC1

Line J1(b)–J3(b)

P6

LEGEND: Tank

Junction

Pump

Air chamber

29

22

30

25

36 33 38

24

32

27

35

37

28

8

48

39

40

18 1

12

4

15

19 7

9

44 45

43 46

41

42

49

47

Figure 10.5: General sktetch and numerical implementation – plant to be designed.

positioned downstream the pumps close in about 6 s, even if the pumps continue their rotating motion. Finally, we checked these results against the Michaud’s formula. As well known, this formula is rigorously valid only for a single pipe; however, in this case we let

164 WATER HAMMER SIMULATIONS Table 10.1: Pumps characteristics.

Pump P1 P2 P3 P4 P5 P6

nreg (rpm) 1450 1450 1450 1450 1450 1450

Qreg (m3 /h) 1500 1500 400 1500 1500 1500

Hreg (m) 78 78 85 85 85 85

Table 10.2: Pipes characteristics. From node P1 P2 P3 AC1 AC1 J1 J1 P4 P5 P6 AC2 AC2 J2 J2 J1 (a) J1 (b) J3

To node AC1 AC1 AC1 Air cham. J1 O1 O2 AC2 AC2 AC2 Air cham. J2 J1 J3 J3 (a) J3 (b) O3

Diameter (m) 0.40 0.40 0.35 0.55 0.60 0.60 0.35 0.60 0.60 0.60 0.60 0.80 0.60 0.60 0.60 0.35 0.90

Length (m) 3 3 3 12 35 525 525 3.5 3.5 3.5 10 15 14 540 525 525 30

n Thick. (m−1/3 /s) (mm) 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6 80 6

Cel. (m/s) 939 939 976 850 826 826 976 826 826 826 826 745 826 826 826 976 713

L = 550 m, V = 1.6 m/s and obviously Tc = 6 s, obtaining: h =

2 · 550 m × 1.6 m/s ∼ 2 · L · V0 = = 30 m g · Tc 9.806 m/s2 × 6

As expected, the carried out value does not perfectly agree with those obtained through the simulation code.

CASE STUDIES

165

Table 10.3: Nodes characteristics. Node P1 P2 P3 AC1 J1 O1 O2 P4 P5 P6 AC2 J2 J3 O3

Elevation (m) Notes 1137 1137 1137 1133.75 Air chamber: 6.37 m3 1136.62 1205.147 1205.147 1136 1136 1136 1135.69 1135.69 1135.69 1204.78

1230 1225 1220

Head (m asl)

1215 1210 1205 1200 1195 1190

Pipe 15–upstream

1185

Static head (downstream tank)

1180 0

10

20

30

40

50

60

70

80

90

100

Time (s)

Figure 10.6: Existing plant: pressures. 10.2.3 Plant to be designed In the following, results of simulations for the plant to be designed are shown having considered three cases: no air chamber, an air chamber of 3 m3 and an air chamber of 5 m3 .


Head (m asl)

1215 1210 1205 1200 1195 1190

Pipe 44–upstream

1185


1180 0

10

20

30

40

50

60

70

80

90

100

Time (s)

Figure 10.7: Pressures in the existing pipes, having eliminated the air chambers.

At first, results presented are related to pipes 44 and 45 because they are already existing and the new upgrade must not worsen the present situation (in the figures only the values for pipe 44 appear as pipe 45 values overlap). Then results for the pipes 38, 39 and 48 are also presented, as they are the main pipes in the new plant. Note that the referred air chamber volume is related to steady flow conditions, while it obviously changes during transients. Also in this case simulations have been carried out considering the presence of a clapet positioned downstream the pumps, which closes when water velocity gets negative. It has to be noted that in some points the piezometric line intersects the pipe: this means that sometime, during transients, pressures can fall below atmospheric. In these cases, vent valves open and let air entry the pipe; this air can change the system’s behaviour, as discussed in Chapter 8.

10.2.4 No air chambers As expected, increasing the discharge, and so the velocity, brings to higher pressures to be recorded during transients. However, changes are not dramatic, and probably the key factors are the short length of the pipe (less than 600 m) and the time of the valve closure (more than 15 s). For the new pipes, similar considerations may be drawn (Figure 10.8). When the velocity becomes negative, the clapet positioned downstream the pump closes and waterhammer pressure oscillations become regular.

CASE STUDIES

167

1230 1225 1220

Head (m asl)

1215 1210 1205 1200 1195 Pipe 38–upstream Pipe 39–downstream Pipe 48–downstream Static head (downstream tank)

1190 1185 1180 0

10

20

30

40

50

60

70

80

90

100

Time (s)

Figure 10.8: Pressures in the new pipes, having eliminated the air chambers. 1230 1225 1220

Head (m asl)

1215 1210 1205 1200 1195 1190 Pipe 44–upstream 1185


1180 0

10

20

30

40

50

60

70

80

90

100

Time (s)

Figure 10.9: Pressures in the existing pipe with 3 m3 air chamber. 10.2.5 Air chamber 3 m3 volume As can be seen in Figure 10.9, due to the distance between the new air chamber and the old pipe, the effect of the former is small. Anyhow, the presence of an air chamber helps in slowing down the system response and in cutting the higher frequency waves.


Head (m asl)

1215 1210 1205 1200 1195 1190 1185 1180 0

10

20

30

40

50

60

70

80

90

100

Time (s)

Figure 10.10: Pressure values in the existing pipe with 5 m3 air chamber. 1230 1225 1220

Head (m asl)

1215 1210 1205 1200 1195 Pipe 38–upstream Pipe 39–downstream Pipe 48–downstream Static head (downstream tank)

1190 1185 1180 0

10

20

30

40

50

60

70

80

90

100

Time (s)

Figure 10.11: Pressure values in the new pipes with 5 m3 air chamber. 10.2.6 Air chamber 5 m3 volume For the reason mentioned in Section 10.2.4 (the distance between the new air chamber and the old pipe) an increase in the air chamber volume has little effect in reducing the overpressure peaks; that can be seen in Figure 10.10 for the existing pipe and in Figure 10.11 for the new pipes.

CASE STUDIES

169

The reduction of the pressures is not dramatic, probably because waterhammer effects are already quite small. The main effect in installing an air chamber in this position is probably the increase of the response time of the system, which become less sensitive to higher frequency perturbations.

10.2.7 Conclusions In this Section unsteady flow simulations are performed with different hypothesis, related to the possible scenarios for a seawater plant to be potentiated in Tanzania. Results show pressure values are within the acceptable range and no cavitation occurs; anyway, in some points of the pipe pressures drop under atmospheric during transients. In these points the positioning of an air valve was designed. Note that the presence of air inside the pipelines changes remarkably the system behaviour in a way that cannot be properly computed with current knowledge. The implementation of the new plant obviously worsen the existing pipe conditions, but it seems this variation is not dramatic. Moreover, it seems that the size of the new air chamber does not affect much the pressure values in the existing pipe, probably because of the distance between these two elements. Anyway, a new air chamber installation is suggested (volume equal to 3 m3 ), in order to smooth the higher frequency perturbations. It has also been shown that the analysis of this plant, performed quite easily with the computer programs developed in this book, could not have been performed with empirical methods.

10.3 A very complex example for seawater treatment plant in Algeria 10.3.1 The plant to be modelled The seawater plant to be modelled is located in Algeria and it is intended to supply potable water to 25% of the Algerian capital city’s population. The project is the first private reverse osmosis potable water project in Algeria and the largest membrane desalination plant in Africa. The plant is fed with open intake Mediterranean seawater. As can be seen in Figure 10.13, the plant comprises nine first-pass trains, each capable of continuously producing 25,100 m3 /day of permeate. In this figure, the intake is depicted by the Storage tanks (S), positioned downstream of which are the main pumps, which guarantee the upstream minimum pressure (NPSH) required for the functioning of the following pumps.

WHS

P

J

J

J

J ERI

J

J

J

J

J

J

J

J

J

J

H

H O

O

J

RO

J

O

J

J

12: 11

J

B RO

H O

ERI

B RO

H O

J ERI

B RO

H O

J ERI

B RO

H

J ERI

B RO

O

J ERI

B

H O

J ERI

RO

H O

FILTER

J

B RO

H

I

J

ERI

B RO

P

J ERI

B

I

J

24/8/2013

I

J

FILTER

CH010.tex

P

170 WATER HAMMER SIMULATIONS REMINERALIZER VESSEL

P

J

J

Junction

P

Main Pump

B

Booster

H

High Pressure Pump

O

Outfall

I

Intake Generic Valve

J

J

J

J

J

J

J

J

O

Page 170

Figure 10.12: Sketch of the seawater treatment plant.

S I

WHS CH010.tex

11

I

15

J 2

12

I

23

19

6

P

21

25

J 36

16

27

J

30

J 52 46 B

J

3

13

I 4

P 14

17

P

10

24

20

26

J

32

31 41

40

57

398

399

402

395

J 401

400

76

88

72 82

B

111

73

114

150

98

99

130

390 387

388

H

393

95

J

392

391

389

J

65

131

B

116

157

193

173

142

407

405 411

410

403

J 408

406

J

108

174

B

416

414

H

181

159

200

J 184

420 415

J

151

210 201 211

J

219

B

419

J 417

228

225 218

202

H

224

227

429 424

J 194

J

243

J

J

276

281

279

428

J 426

262

260

B J

268 261

H

267

245 270

263 271

438 433

J 237

296 287 297

J 308 302

305

B

437

430

J 435

303

314 441

H

431 310

288

329

J

313

447 442

J

280

344 334 343

J

339 330 340

J 351 345

348

B J

443 439

371

O

444

347

J 333

331

J

356

355 349 357

354 446

450

H

440 353

O

436

282

335 326

352

306

J

311 304

373

362

322

J 290

312

309 434

O

427

301 291 300

J

286

J

269

432

422

239

253 244 254

283

J

J

292

J 247

319

324

359

361

J

J

258 248 257

265 259

421

O

418

240

266 425

423

413

196

249

J 204

226 220

217 J

236

316

318

J

J 215 205 214

J

222 216

412

O

409

153

197

223

185

182 175

161

183 177

J

404

O

386

176

180

134

H

167 158 168

179

206

J

233

238

273

275

J

J

172 162 171

J

J

J

140

J

138

110

154

141

139 132

163

J 118

133

383

O

63

124 115 125

136 137

91

190

195

J 129 119 128

J

J

J

97

J

51

67

J 75

90

96 89

120

85

93 94

394

397

J

81

J

147

152

J 86

J

71

87

56

O

22

77 68

J 49

54

J

8

107

230

232

J

456

O

445

451

452 448 455

453

449 454

J

325

323

18 61

105

J

102

J

148

145

J

191

188

J

234

231

J

277

274

J

320

317

J

363

360

J

369

366

374

Page 171

I

7

34

44

396

H

104

109

J 35

45

55 48

64

187

189

J

12: 11

53 47

62

66

144

146

J

24/8/2013

9

101

103

J

5

P

372

O

J Junction P Main Pump B Booster H High Pressure Pump O Outfall I

Intake Generic Valve (Concentrated Headloss)

CASE STUDIES

171

Figure 10.13: Numerical implementation of the seawater treatment plant.

60 1

172 WATER HAMMER SIMULATIONS The membrane elements (RO) of each train are fed by a dedicated high pressure (HP) pump, a booster (B) pump, and an array of 32 ERI®2 model PX-220 Pressure Exchanger energy recovery devices. All nine pumps are fed from a common supply manifold, which is maintained at a pressure greater than 3.5 bar in order to meet the minimum NPSH requirement of the HP pump. The PX device arrays are fed from a second supply manifold, which is maintained at a pressure of greater than 1.3 bar in order to meet the minimum discharge-pressure-requirement of the PX devices. Because the PX device arrays and HP pumps are supplied separately, the PX device arrays are isolated from flow and pressure variations that may occur as individual HP pumps go on- and off-line. The HP pumps selected for the plant operate at 1084 m3 /h and approximately 88% efficiency; the booster pumps operate at 1351 m3 /h and approximately 89% efficiency. The RO system divides the flow into two parts, one of which (permeate) is collected to the remineralizer vessel and then to an outfall storage, while the other (brine) is conveyed through the ERI device to recover energy and then disposed. 10.3.2 A peculiar device: energy recovery PX In a reverse osmosis (RO) water desalination system, a high percentage of water (brine) has to be discarded, consuming a great deal of energy used to obtain the high pressure necessary in these plants. As a consequence, many systems have been developed in order to recover this energy. One of the most advanced devices is implemented in the ERI PX developed for specific use in sea water reverse osmosis plants and transfers pressure from a high pressure line to a low pressure one [113]. The Algerian desalination plant, designed by General Electric – Water section, was at that time the second biggest seawater reverse osmosis desalination plant in the world, and it used one array of 32 PX-220 devices on each of the nine first-pass trains for a total of 288 PX units, the purpose being to minimize the foreseen high costs of processing the expected 200,000 m3 /day of water. The first pass train size is to be 25,100 m3 /day, so that the plant is able to operate at full capacity with just eight trains running [114]. In order to build a complete model and to simulate waterhammer processes, a laboratory investigation of the PX behaviour during flow transients was performed, as its behaviour in unsteady flow conditions was never been tested before, and these results were used for the subsequent simulations. 10.3.3 Laboratory plant In order to verify the behaviour of PX during unsteady flow, several experimental tests were performed in the Hydraulic Laboratory of the Politecnico di Milano [115]. 2

PX, PX Pressure Exchanger, ERI and the ERI logo are registered trademarks of Energy Recovery, Inc.

CASE STUDIES

173

Low pressure line Opening / closing valve

Tank Flowmeter

Pressure transducer ERI-PX Flowmeter

Pressure transducer

Pumps High pressure line

Outlet

Figure 10.14: Sketch of a configuration of the experimental plant (valve on the high pressure line).

A system of pumps brought water to the upstream tank positioned about 9 m above the underground floor. A pipe with diameter D1 = 300 mm was derived from this tank; a reduction brought the diameter to D2 =150 mm and a third reduction brought the diameter to the final one D3 = 2 . This was the low pressure circuit. A similar high pressure circuit was regulated by two high pressure pumps (see Figures 10.14 and 10.15). The final parts of the plant consisted of two iron pipes with diameter D3 = 2 of spiral shape with 26 coils for a total length L of about 90 m. These pipes were installed in order to allow the development of the waterhammer waves in all reflection phases. Celerity of the wave was evaluated as equal to 1388 m/s. There followed two electromagnetic flowmeters, an electro-valve able to close and open almost instantaneously, and the pressure transducers connected to the data collector and the computer. Tests were performed with the valve situated in different positions on the two lines, in order to fully investigate the device’s behaviour.

10.3.4

Laboratory tests

Numerous tests were performed by opening and closing the valve in different positions: on the high pressure or low pressure lines, and upstream or downstream of the coils (in these cases always downstream of the ERI device). Moreover, the valve

174 WATER HAMMER SIMULATIONS High pressure pumps Opening/closing valve

Pressure transducer

Flowmeters

ERI–PX 220

Figure 10.15: Laboratory plant. The valve is located on the low pressure line upstream of the coil and downstream of the EXI – PX. was also positioned upstream of the PX. In what follows, only the tests judged most significant are described. During the tests, discharges were set the same for the two lines and equal to 3 l/s, which gave a velocity equal to 1.41 m/s. Figure 10.16 shows the results of a test with the closure valve positioned downstream of the spiral pipe of the high pressure line. As expected, a large development of pressures on the high pressure line can be observed. On the other hand, overpressures on the low pressure line are also noticeable, but they are almost negligible. The transient ends in a few seconds before the pumps reach the static head. The reliability can be verified by comparing the peak recorded (Figure 10.16) with the expected value computed with the Allievi–Joukowski formula:

h =

c 1388 ·V = · 1.41 = 200 m ≈ 20 bar g 9.806

Similar tests were performed with the closure valve positioned on the low pressure stream. Moreover, the valve was opened on the two lines so as to have initial depression values. The results show that in this case too the disconnection was almost complete. When the valve on the high pressure line was opened, a depression occurred in the same line, but the effect on the low pressure line was practically negligible. When the valve on the low pressure line was opened, the effect on both lines was insignificant. Repeatability tests were performed in order to verify the coherence between results. In fact, in Figure 10.16 two tests are reported, performed with the same conditions and which prove that the results are equivalent.

CASE STUDIES

175

25 HP line - Test 1 LP line - Test 1 HP line - Test 2 LP line - Test 2

Pressure (bar)

20

15

10

5

0 0

1

2

3

4

5

6

7

8

9

10

Time (s)

Figure 10.16: Pressure values recorded in the high and low pressure lines. Instantaneous closure on the high pressure line. 10.3.5 Model of the seawater plant in Algeria Analysis of the laboratory test results showed that, during transients, high pressure and low pressure lines could be considered completely disconnected. This simplified the implementation of the mathematical model to simulate the actual plant. In this case, the plant schematically represented in Figure 10.12 can be divided into two independent parts where the ERI – PX is located [115]. The resulting computing scheme is shown in Figure 10.13. Twenty tests, as listed in Table 10.4, were performed with this mathematical model. The results obtained for all tests showed that the plant was at no risk of dangerous pressures or depressions even in the worst situations which could be foreseen. For instance, simulation 5 modelled the stoppage of both the main pumps positioned in the “lower” part of the plant (see Figure 10.13). As can be seen in Figure 10.17, the pressure dropped and then started to oscillate in typical fashion, but no problematic values were recorded. A similar test (simulation 3) was performed by modelling the stoppage of the two main pumps in the “upper” part of the plant (see again Figure 10.13). Owing to the complete separation of the two parts of the plant produced by the ERI device, there was no need to perform crossed stoppage (stoppage of one pump in the upper part and one in the lower part; stoppage of both main pumps in the upper part and both in the lower, and so on). Simulation 10 modelled the whole train stoppage which means that the high pressure and booster pumps stopped and the valve positioned on the RO system is closed. A sudden fall in the pressure downstream of the RO could be observed,

176 WATER HAMMER SIMULATIONS Table 10.4: List of all tests performed on the Hamma plant with the mathematical model. Device whose trip Test # was modelled 1 Remineralizer vessel 2 Upstream 3 pumps 4 5 6 7

Train 1

8 9 10

11

12 13 14

Train 4

Details about the model’s rules Closure valve Vrem : Tc = 5 s Stopping of pump P1,sup Stopping of pumps P1,sup and P1,inf Stopping of pump P2,sup Stopping of pumps P2,sup and P2,inf Stopping of pumps P1,sup – P1,inf – P2,sup – P2,inf Stopping of pump PHP The following operations were simultaneously performed: a. Stopping of pump PHP b. Closure of valve VRO,PERM : Tc = 5 s Stopping of pump PBOOST RO’s breakage, VRO,PERM : Tc = 5 s Whole train stopping The following operations were simultaneously performed: a. Stopping of pump PHP b. Stopping of pump PBOOST c. Closure of valve VRO,PERM : Tc = 5 s Stopping of pump PHP The following operations were simultaneously performed: a. Stopping of pump PHP b. Closure of valve VRO,PERM : Tc = 5 s Stopping of pump PBOOST RO’s breakage VRO,PERM : Tc = 5 s Whole train stopping The following operations were simultaneously performed: a. Stopping of pump PHP b. Stopping of pump PBOOST c. Closure of valve VRO,PERM : Tc = 5 s

Nodes of the mathematical model involved Node 366 Node 5 Nodes 5 and 6 Node 7 Nodes 7 and 8 Nodes 5 – 6 – 7–8 Node 48 Node 394

Node 46 Node 394

Node 48 Node 46 Node 394

Node 175 Node 177 Node 173 Node 177

Node 175 Node 173 Node 177 (Continued)

CASE STUDIES

177

Table 10.4: Continued. Device whose trip Test # was modelled Details about the model’s rules 15 Train 8 Stopping of pump PHP The following operations were simultaneously performed: a. Stopping of pump PHP b. Closure of valve VRO,PERM : Tc = 5 s 16 Stopping of pump PBOOST 17 RO’s breakage RO VRO,PERM : Tc = 5 s 18 Whole train stopping The following operations were simultaneously performed: a. Stopping of pump PHP b. Stopping of pump PBOOST c. Closure of valve VRO,PERM : Tc = 5 s 19 All trains Stop of all trains at the same time The following operations were simultaneously performed: a. Stopping of pump PHP b. Stopping of pump PBOOST c. Closure of valve VRO,PERM : Tc = 5 s d. Stopping of pump PHP e. Stopping of pump PBOOST f. Closure of valve VRO,PERM : Tc = 5 s g. Stopping of pump PHP h. Stopping of pump PBOOST i. Closure of valve VRO,PERM : Tc = 5 s j. Stopping of pump PHP k. Stopping of pump PBOOST l. Closure of valve VRO,PERM : Tc = 5 s m. Stopping of pump PHP n. Stopping of pump PBOOST o. Closure of valve VRO,PERM : Tc = 5 s p. Stopping of pump PHP q. Stopping of pump PBOOST r. Closure of valve VRO,PERM : Tc = 5 s s. Stopping of pump PHP t. Stopping of pump PBOOST u. Closure of valve VRO,PERM : Tc = 5 s

Nodes of the mathematical model involved

Node 347 Node 349 Node 345 Node 349

Node 347 Node 345 Node 349

Node 48 Node 46 Node 394 Node 89 Node 87 Node 91 Node 132 Node 130 Node 134 Node 175 Node 173 Node 177 Node 218 Node 216 Node 220 Node 261 Node 259 Node 263 Node 304 Node 302 Node 306 (Continued)

178 WATER HAMMER SIMULATIONS Table 10.4: Continued. Device whose trip Test # was modelled Details about the model’s rules v. Stopping of pump PHP w. Stopping of pump PBOOST x. Closure of valve VRO,PERM : Tc = 5 s 20 Electricity failure– stopping of all pumps and VRO,PERM closure. The following operations were simultaneously performed: a. Stopping of pumps P1,sup – P1,inf – P2,sup – P2,inf b. Stopping of pumps PHP

c. Stopping of pumps PBOOST

d. Closure of valves VRO,PERM : Tc = 5 s

Nodes of the mathematical model involved Node 347 Node 345 Node 349

Nodes 5 – 6 – 7–8 Nodes 48 – 89 – 132 – 175 – 218 – 261 – 304 – 347 Nodes 46 – 87 – 130 – 173 – 216 – 259 – 302 – 345 Nodes 394 – 91 – 134 – 177 – 220 – 263 – 306 – 349

70 60

Pressure (m)

50 40 30 20 10 0 0

5

10

15

20

Time (s)

Figure 10.17: Pressure values inside the pipe 17 (pipe 18 values overlap) simulation 5.

CASE STUDIES

179

700 600

Pressure (m)

500 400 Pipe 452 300 200

Pipe 353

100 0 0

2

4

6

8

10

12

14

16

18

20

Time (s)

Figure 10.18: Pressure values inside the plant, simulation 19.

but even upstream of that system, where overpressures were expected, values were within the acceptable range. Figure 10.18 reports the results of simulation 19 when all trains were closed (and hence was probably one of the most dangerous situations that may occur).

10.3.6 Conclusions The Algerian desalination plant designed by General Electric – Water section, the largest seawater reverse osmosis desalination plant in Africa and the second largest in the world, has been simulated using the methods described in this book. Such a large plant requires careful design, and a hydraulic model for steady and unsteady flow simulations has been developed accordingly. A laboratory investigation was performed in order to investigate the ERI – PX’s behaviour during transients, and its results have been reported in this book. The laboratory results clearly show that the PX produced an almost perfect hydraulic disconnection between high and low pressure lines (although from a quality point of view the two lines were obviously not disconnected), and this simplified the mathematical model implementation. A sketch of the Algerian plant and its mathematical representations have also been provided, with some of the results obtained from thorough investigation of the problems most likely to occur. The results show that the plant can be considered safe, from the hydraulic point of view, in case of both steady and unsteady flows.


10.4 Final remarks This paragraph presented the applications of the methods to analyze complex plants during transients, as described in this book. Simplest plants, as the first presented in this paragraph for petroleum product in Djibouti, show that the simplest models can be applied and the differences with the results carried out with most complex models are negligible. However, when the complexity of the plant increases, the oldest and simplest models are not yet sufficient to describe the behavior of the plant, and the carried out results are not appropriate to establish the measure to contrast the effects of waterhammer. As it has been shown for the plant in Tanzania, the design of some devices, in that case it was an air chamber, has to be performed through more complex models: with them, different solutions must be tried in order to see the effect of the introduced modifications and to select the most appropriate. Some very complex plants are very difficult to be modeled. The Algerian plant, shown as third example, included devices which behavior was not known in unsteady flow conditions; this required the performance of laboratory tests. After the tests were carried out, the numerical model has been built and the carried out results allowed the evaluation of the plant performances. However, field tests should always be scheduled and performed before starting the normal activities of the plant because, as we tried to highlight in the whole book, the topic is quite complex and not still resolved, and therefore a deep analysis should be performed in most practical cases.

References

[1] Menabrea L.F., 1858, “Note sur les effects de choc de l’eau dans les conduits.” C. R. Acad. Sci. Paris, 47, 221–224. [2] Michaud J., 1878, “Coup de bélier dans les conduits – Etudedes moyens employés pour en atténuer les effets.” Bulletin de la Société Vaudoise des Ingénieurs et des Architects, 4(3), 56–64; 4(4), 65–77. [3] Joukowsky N., 1898, “Waterhammer.” Mem. Imp. Acad. Soc. St. Petersburg – translated by O. Simin, 1904, Proc. Am. Waterworks Assoc., 24, 341–424. [4] Allievi L., 1903, “Teoria generale del moto perturbato dell’acqua nei tubi in pressione.” Annali della Società degli Ingegneri ed Architetti italiani, Milano. [5] Allievi L., 1913, “Teoria del colpo d’ariete” Atti del Collegio degli Ingegneri e degli Architetti italiani. [6] Evangelisti G., 1965, “Teoria generale del colpo d’ariete col metodo delle caratteristiche” L’Energia Elettrica, Vol. XLII, Milano. [7] Evangelisti G., 1969, “Waterhammer Analysis by the Method of Characteristics”, L’Energia Elettrica Vol. XLVI, Milano. [8] Angus R.W., 1935, “Simple Graphical Solution for Pressure Rise in Pipes and Pump Discharge Lines,” J. Eng. Inst. Canada, 72–81, February. [9] Bergeron L., 1935, “Etude des variations de régime dans les conduits d’eau: Solution graphique générale,” Rev. Gen. Hydraul., Paris, 1, 12–25. [10] Streeter V.L., Wylie E.B., 1967, Hydraulic Transients, McGraw Hill, New York. [11] Streeter V.L., Wylie E.B., 1993, Fluid Transients in Systems, Prentice Hall, New Jersey. [12] Chaudry M.H., 1979, Applied Hydraulic Transients, Van Nostrand Reinhold, New York. [13] Brunone B., Ferrante M., 2001, “Detecting leaks in pressurised pipes by means of transients.” Journal of Hydraulics Research, IAHR, 39(5). [14] Mpesha W., Chaudhry M.H., Gassman S.L., 2002, “Leak detection in pipes by frequency response method using a step excitation.” Journal of Hydraulic Research, 40, 11.1 (55–62). [15] Wang X.W., Lambert M.F., Simpson A.R., Liggett J.A., Vitkovsky J.P., 2002, “Leak Detection in Pipelines using the Damping of Fluid Transients.” Journal of Hydraulic Engineering. ASCE, 128(7), 697–711. [16] Kapelan Z.S., Dragati S.A., Godfrey W.A., 2003, “A hybrid inverse transient model for leakage detection and roughness calibration in pipe network.” Journal of Hydraulic Research, 41(5), 481–492. [17] American Society of Civil Engineers ASCE, 1975, Pipeline design for water and wastewater, Report of the Task Committee on the Engineering Practice in the Design of Pipelines, Philadelphia, PA.

182 REFERENCES [18] American Society of Civil Engineers ASCE, 1992, Pressure Pipe Design for Water and Wastewater, Committee on Pipeline Planning of the Pipeline Division of ASCE, New York, NY. [19] Ellis J., 2008, Pressure transients in water engineering, Thomas Telford Publishing, London. [20] Cantù M., 2001, Mastering Delphi 6, Sybex Inc., Alameda, CA. [21] Stephenson D., 1989, Pipeline design for water engineers, Elsevier Science Amsterdam. [22] Mays L.W., 2000, Water distribution system handbook, The McGraw-Hill Companies, New York. [23] Todini E., 2006, “On the convergence properties of the different pipe network algorithms”, 8th Annual Water Distribution Systems Analysis Symposium, Cincinnati, OH, August 27–30. [24] Mariotte E., 1686, Traité du mouvement des eaux et autres corps fluids, Paris, chéz Estienne Michallet. Edme Mariotte died in 1684, and the book has been edited by M. de La hire. [25] D’Alembert J.R., 1747, Recherches sur les Cordes Vibrantes, Paris. [26] Evans L., 1998, Partial differential equations. American Mathematical Society, Providence, 740 pp. [27] Bergeron L., 1961, Water hammer in hydraulics and wave surges in electricity (translated under the sponsorship of the ASME), Wiley, New York. [28] Parmakian J., 1963, Water-hammer analysis, Dover, New York. [29] Lister M., 1960, “The numerical solution of hyperbolic partial differential equations by the method of characteristics”, in A. Ralston and H.S. Wilf (Eds.) Numerical methods for digital computers, Wiley, New York. [30] Abbott M.B., 1966, An introduction to the method of characteristics, Thames and Hudson, London. [31] Hirsch C., 2007, Numerical computation of internal and external flows, 2nd Edition, Elsevier, Oxford, UK. [32] Godunov S.K., Ryabenkii V.S., 1987, Difference schemes, Elsevier, Amsterdam. [33] Lax P.D., Wendroff B., 1960, “System of conservation laws”, Comm. Pure Appl. Math., 13, 217–237. [34] Thomas J., 1995, Numerical partial different equations: finite difference methods, Springer-Verlag, New York. [35] Roe P.L., 1981, “Approximate Riemann solvers, parameter vectors and difference schemes”, J. Comput. Phys, 43, 357–372. [36] LeVeque R.J., 2002, Finite volume methods for hyperbolic problems, Cambridge University Press. [37] Maione U., 1961, “Perdite di carico nella strozzatura di un pozzo piezometrico”, L’Energia Elettrica, 38, 330–338. [38] Thoma D., 1910, Zur Theorie des Wasserschlosses bei Selbsttaetig Geregelten Turbinenanlagen, Oldenburg, Munchen, Germany. [39] Jaeger C., 1958 “Contribution to the stability theory of systems of surge tanks”, Trans. ASME, 80, 1574–1584. [40] Jaeger C., 1960, “A review of surge-tank stability criteria”, J. Basic Eng., Trans. ASME, December, 765–783. [41] Pickford J., 1969, Analysis of water surge, Gordon and Breach Publisher, New York.

REFERENCES

183

[42] Pearsall I.S., 1963, “Comparative experiments on surge tank performance”, Proc. National Engineering Laboratory, 177, 951–970. [43] Ellis J., Khairulla L.M., 1974, “Oscillations in a surge tank with upper and lower expansion galleries,” Water Power, November, 359–364. [44] Noseda G., sine data, Problemi di moto vario, Ed. Istituto di Idraulica e Costruzioni Idrauliche, Politecnico di Milano. [45] Johnson R.D., 1915, “The differential surge tank”, Trans. Am. Soc. Civ. Eng., 78, 760–784. [46] Gibson W.L., Shelson W., 1956, “An experimental and analytical investigation of a differential surge-tank installation”, Trans. A. Soc. Mech. Eng., 78, 925–938. [47] Angus R.W., 1937, “Air-chambers and valves in relation to water hammer”, Trans. Am. Soc. Mech. Eng., 59, 661–668. [48] Paoletti A., 1972, “Il transitorio negli impianti elevatori muniti di casse d’aria”, L’Energia Elettrica, 6. [49] White F.M., 2003, Fluid mechanics, 4th edition, McGraw Hill Series in Mechanical Engineering, New York. [50] Evangelisti G., 1935, “Sul calcolo delle oscillazioni di carico nelle condotte degli impianti di sollevamento” L’Elettrotecnica, June. [51] Evangelisti G., 1938, “Il colpo d’ariete nelle condotte elevatorie munite di camere d’aria” L’Energia Elettrica, 9. [52] Cozzo G., 1973, “Azione di una piastra piana sull’efflusso da un boccaglio”, Memorie e studi dell’Istituto di Idraulica e Costruzioni Idrauliche del Politecnico di Milano, 1973. [53] Tanda M.G., Zampaglione D., 1991, Analisi di fenomeni di risonanza in un sistema distributore al termine di un’adduttrice a gravità Dipartimento di Ingegneria Idraulica, Ambientale e del Rilevamento, Politecnico di Milano. [54] Bianchi A., Mambretti S., 2000, “Le idrovalvole di scarico rapido per l’attenuazione delle sovrapressioni di colpo d’ariete”, L’Acqua, 6. [55] De Martino G., 1973, “Sul calcolo del GD2 neli impianti di sollevamento”, L’Energia Elettrica, 8. [56] Mambretti S., Orsi E., 2007, “Un modello di simulazione dei transitori di un sistema di pompaggio.” Proc. of the Third Seminar “La ricerca delle perdite e la gestione delle reti di acquedotto”, 185–192, Morlacchi editore, Perugia. [57] Streeter V.L., Wylie E.B., 1985, Fluid mechanics, 8th Edition, McGraw-Hill, New York. [58] Mambretti S., 2004, Fenomeni di moto vario nelle correnti in pressione, pp. 128, Aracne editrice, Roma. [59] Miller D.S., 1978, Internal flow systems, Volume 5 in the BHRA Fluid Engineering Series, Cranfield. [60] Remenieras G., 1952, “Dispositif simple pour réduire la célérité des ondes élastiques dans le conduits en charge”, La Houille Blanche, Special A, 172–196. [61] Ghilardi P., Paoletti A., 1986, “Additional viscoelastic pipes as pressure surge suppressors”, Paper E2, Proc. 5th Int. Conf. Pressure Surges, BHRA, Hannover, Germany. [62] Pezzinga G., 2002, “Unsteady Flow in Hydraulic Networks with Polymeric Additional Pipe”, J. Hydraulic Eng., 128(2), 238–244.

184 REFERENCES [63] Meniconi S., Brunone B., Ferrante M., Massari C., 2012, “Transient hydrodynamics of in-line valves in viscoelastic pressurized pipes: long-period analysis”, Exp. Fluids, 52, 265–275. [64] Kolkman P.A., 1984, “Gate vibrations” in P. Novak (Ed.) Developments in hydraulic engineering – 2, Elsevier Applied Science Publishers, Barking, UK. [65] Tongue B.H., 2002, Principles of vibration, 2nd edition, Oxford University Press. [66] Inman D.J., 2008, Engineering vibration, 3rd Edition, Pearson Prentice Hall, New Jersey. [67] Jaeger C., 1963, “The theory of resonance in Hydropower systems. Discussion of incidents and accidents occurring in pressure systems.” Trans. ASME, Ser. D, 85, 631. [68] Chaudry M.H., 1979, “Resonance in pressurized piping systems”, J. Hydraulic Div. ASCE, 96, HY 9, 1819–1839. [69] Howell K.B., 2001, Principles of Fourier analysis, CRC Press, USA. [70] Hedstrom G., 1975, “Models of difference schemes for ut + ux = 0 by partial differential equations”, Math. Comp., 29, 969–977. [71] Quarteroni A., 2008, Modellistica numerica per problemi differenziali, SpringerVerlag, Milano. [72] Harten A., Lax P.D., 1984, “On a class of high resolution total variation stable finite-difference schemes”, SIAM J. Numer. Anal., 21(1), 1–23. [73] Godunov. S.K., 1959, “A finite difference method for the computation of discontinuous solutions of the equations of fluid dynamics”, Matematicheskii Sbornik, 47, 357–393. [74] Jameson A., 1982, “Transient aerofoil calculations using the Euler equation”, in P.L. Roe (ed.), Numerical methods in aeronautical fluid dynamics, New York. [75] Pulliam T.H., 1984, “Artificial dissipation models for the Euler equation”, AIAA Paper 85-0438, AIAA 23rd Aerospace Sciences Meeting, Reno, Nevada. [76] McCormack R. W., Baldwin B.S., 1975, “A numerical method for solving the Navier-Stokes equation with application to shock-boundary layer interaction”, AIAA Paper, 1–75. [77] Narayanan K.V., 2006, A Textbook of Chemical Engineering Thermodynamics, 8th Printing, Prentice-Hall, New Delhi. [78] Van De Sande E., Belde A.P., Hamer B.J.G., Hiemstra W., 1980, “Velocity profiles in Accelerating pipe Flows started from the Rest”, Proceeding of the 3rd International Conference on Pressure Surges, BHRA Fluid Engineering, Canterbury, England, A1, 1–14. [79] Bergant A., Ross Simpson A., Vitkovsky J., 2001, “Developments in unsteady pipe friction modelling”, J. of Hydraulic Res., 39(3), 249–257. [80] Brunone B., Greco M., 1990, “Un modello per la ricostruzione di fenomeni di colpo d’ariete anche in presenza di cavitazione. Riscontro sperimentale”, Proceedings of 22nd Convegno di Nazionale di Idraulica e Costruzioni Idrauliche, Cosenza, Italy, 4, 114–160. [81] Brunone B., Golia U.M., Greco M., 1991, “Some remarks on the momentum equation for the fast transient”, Proceedings of International Meeting on Hydraulic Transient and Water Column Separation, Valencia, IAHR, E. Cabrera, M.A. Fanelli (eds), 201–205. [82] Brunone B., Golia U.M., Greco M., 1991, “Modelling of fast transient by numerical methods”, Proceedings of International Meeting on Hydraulic Transient and Water Column Separation, Valencia, IAHR, E. Cabrera, M.A. Fanelli (eds), 273–282.

REFERENCES

185

[83] Pezzinga G., 2000, “Affidabilità di modelli semi-empirici di turbolenza per la valutazione delle resistenze di attrito in condizioni di moto vario”, L’Acqua, 78(1), 29–38. [84] Ramos H., Covas D., Borga A., Loureiro D., 2004, “Surge damping analysis in pipe systems: modelling and experiments”, J. Hydraul. Res., 42(4), 413–425. [85] Bughazem M. B., Anderson A., 1996, “Problems with simple models for damping in unsteady flow”, Proceedings of the 7th International Conference on Pressure Surges and Fluid Transients in Pipelines and Open Channels, BHR Group, Harrogate, UK, 537–548. [86] Wylie E. B., 1996, “Frictional effects in unsteady turbulent pipe flow”, Applied mechanics in the Americas, M. Rysz, L.A. Godoy, L.E. Suarez (eds), Vol. 5, The University of Iowa Press, Iowa City, Iowa, 5, 29–34. [87] Pezzinga G., 2000, “Evaluation of Unsteady Flow Resistances by Quasi-2D or 1D Models”, J. Hydraul. Eng., ASCE, 126(10), 778–785. [88] Vardy A.E., Brown J.M.B., 1996, “On turbulent, un steady, smooth-pipe flow”, International Conference on Pressure Surges and Fluid Transients, BHR Group, Harrogate, England, 289–311. [89] Carravetta A., Golia U.M., Greco M., 1992, “Sull’attenuazione spontanea delle fluttuazioni di pressione durante i transitori di colpo d’ariete”, Atti del XXIII Convegno di Idraulica e Costruzioni Idrauliche, Firenze, E67–E79. [90] Golia U.M., 1992, “Sulla valutazione delle Forze Resistenti nel Colpo d’Ariete”, Pubblicazione n .639 del Dipartimento di Idraulica, Università degli Studi di Napoli “Federico II”. [91] Bratland O., 1986, “Frequency-dependent friction and radial kinematic energy variation in transient pipe flow”, Proceedings of the 5th International Conference On Pressure Surge, BHRA, 95–101. [92] Modica C., Pezzinga G., 1992, “Un modello quasi bidimensionale per il moto vario elastico in regime turbolento”, Atti del XXIII Convegno di Idraulica e Costruzioni Idrauliche, vol. 4, E191–E205, Firenze. [93] Pezzinga G., 1999, “Quasi-2D model for Unsteady Flow in Pipe Networks”, Journal of Hydraulic Engineering, 125(7), 676–685. [94] Vardy A.E., Hwang K.L., 1991, “A characteristics model of transient friction in pipes”, Journal of Hydraulic Engineering, 29(5), 669–684. [95] Prandtl, L., 1925, Z. angew. Math. Mech., 5(1), 136–139. [96] Nikuradse J., 1932, “Laws of turbulent flow in smooth pipes”, NASA TT F-10, 359 (Translation of VDI-Forsch). [97] Nikuradse J., 1933, “Laws of flow in rough pipes”, NASA TM 1292, 359 (Translation of VDI-Forsch). [98] Jones W.P., Launder B.E., 1972, “The prediction of laminarization with a twoequation model of turbulence”, International Journal of Heat and Mass Transfer, 15, 301–314. [99] Eichinger P., Lein G., 1992, “The influence of friction on un steady pipe flow”, Unsteady flow and fluid transients, R. Bettes, J. Watts (eds), Balkema, Rotterdam, The Netherlands, 41–50. [100] Pezzinga G., Brunone B., 2006, “Turbulence, fiction and Energy dissipation in transient pipe flow”, Vorticity and turbulence effects in fluid structures interactions: an application to hydraulic structure design, M. Brocchini, F. Trivellato (eds), WIT Press, UK, 213–236.

186 REFERENCES [101] Zhao M., Ghidaui M.S., 2006, “Investigation of turbulence behaviour in pipe transients using a k–ε model”, Journal of Hydraulic Research, 44(5), 682–692. [102] Fan S., Lakshminarayana B., Barnett M., 1993, “Low-Reynolds number k–ε model for unsteady turbulent boundary-layer flows”, AIAA Journal, 31(10), 1777–1784. [103] Lam, C.K.G., Bremhorst, K.A., 1981, “Modified form of the k–ε model for predicting wall turbulence.” Journal of Fluids Engineering, 103(3), 456–460. [104] Wylie E.B., 1984, “Simulation of Vaporous and Gaseous Cavitation”, Journal of Fluids Engineering, ASME, 106, 307–311. [105] Cannizzaro D., Pezzinga G., 2002, “Influenza del rilascio di gas sulle dissipazioni in transitorio con cavitazione”, Proc. XXVIII Convegno di Idraulica e Costruzioni Idrauliche, Potenza, 1, 619–626. [106] Zielke W., Perko H.D., Keller A., 1990, “Gas Release in Transient Pipe Flow”, Proceedings of the 6th International Conference on Pressure Surge, BHRA, C1, Cambridge, 3–13. [107] Wylie E.B., 1992, “Low void fraction two-component two-phase flow”, in R. Bettes, J. Watts (eds), Unsteady flow and Fluid Transients, Balkema, Rotterdam, 3–9. [108] Scaccia M., Cannizzaro, D., Pezzinga, G., 2004, “Calibrazione di modelli di cavitazione gassosa per mezzo di micro-algoritmi genetici”, Proc. XXIX Convegno di Idraulica e Costruzioni Idrauliche, Trento, 299–305. [109] Ewing D. J. F., 1980, “Allowing for free air in water hammer analysis”, Proceedings of the 3rd International Conference on Pressure Surges, BHRA, Canterbury, UK, 127–146. [110] Pezzinga G., 2003, “Second viscosity in transient cavitating pipe flows”, Journal of Hydraulic Research, 41(6), 656–665. [111] Cannizzaro, D., Pezzinga, G., 2005, “Energy dissipations in transient gaseous cavitation,” Journal of Hydraulic Engineering, ASCE, 131, 724–732. [112] Ming Z., Ghidaoui M.S., 2003, “Efficient quasi-two-dimensional model for water hammer problems”, Journal of Hydraulic Engineering, 1007–1013. [113] Stover R.L., 2004, “Development of a fourth generation energy recovery device. A CTO’s Notebook”, Desalination, 165, 313–321. [114] Stover R.L., Martin J., Nelson M., 2007, “The 200,000 m3 /day Hamma seawater desalination plant largest single-train SWRO capacity in the world and alternative to pressure center design”, Proceedings IDA World Congress, Gran Canaria, Spain. [115] Mambretti S., Orsi E., Gagliardi S., Stover R., 2009, “The behaviour of energy recovery devices in unsteady flow conditions and application in the modelling of the Hamma desalination plant”, Desalination, 238, 233–245.

...for scientists by scientists

Stochastic Methods in Engineering I. DOLTSINIS, University of Stuttgart, Germany

The increasing industrial demand for reliable quantification and management of uncertainty in product performance forces engineers to employ probabilistic models in analysis and design, a fact that has occasioned considerable research and development activities in the field. Notes on Stochastics eventually address the topic of computational stochastic mechanics. The single volume uniquely presents tutorials on essential probabilistics and statistics, recent finite element methods for stochastic analysis by Taylor series expansion as well as Monte Carlo simulation techniques. Design improvement and robust optimisation represent key issues as does reliability assessment. The subject is developed for solids and structures of elastic and elasto-plastic material, large displacements and material deformation processes; principles are transferable to various disciplines. A chapter is devoted to the statistical comparison of systems exhibiting random scatter. Where appropriate examples illustrate the theory, problems to solve appear instructive; applications are presented with relevance to engineering practice. The book, emanating from a university course, includes research and development in the field of computational stochastic analysis and optimization. It is intended for advanced students in engineering and for professionals who wish to extend their knowledge and skills in computational mechanics to the domain of stochastics. Contents: Introduction, Randomness, Structural analysis by Taylor series expansion, Design optimization, Robustness, Monte Carlo techniques for system response and design improvement, Reliability, Time variant phenomena, Material deformation processes, Analysis and comparison of data sets, Probability distribution of test functions. ISBN: 978-1-84564-626-4 eISBN: 978-1-84564-627-1 Published 2012 / 378pp / £158.00

WIT Press is a major publisher of engineering research. The company prides itself on producing books by leading researchers and scientists at the cutting edge of their specialities, thus enabling readers to remain at the forefront of scientific developments. Our list presently includes monographs, edited volumes, books on disk, and software in areas such as: Acoustics, Advanced Computing, Architecture and Structures, Biomedicine, Boundary Elements, Earthquake Engineering, Environmental Engineering, Fluid Mechanics, Fracture Mechanics, Heat Transfer, Marine and Offshore Engineering and Transport Engineering.


Water Hammer Simulations by S. Mambretti

Recommend Documents