# Hydraulic Simulation of Leakage in Post-Earthquake Water Supply Pipeline Network and Its Algorithm Improvement

Authors: Tianyang Yu
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

For the post-earthquake water supply pipeline network, the method of generating random numbers to determine the pipelines with leakage damage according to the earthquake damage rate is adopted. Select the most suitable seismic damage rate calculation model for this study. The proportion of different leakage points of pipelines with different pipe materials and interface forms summarized by Shi and O’Rourke from the actual earthquake damage is used to randomly sample the pipelines that cause earthquake damage, determine the leakage loss form of each leakage point, and use the corresponding leakage flow model and the leakage point area model for hydraulic analysis of the leakage of the post-earthquake pipeline network. The pressure-dependent nodal demand (PDD) model is added to the hydraulic analysis to avoid the problem of node negative pressure in the hydraulic analysis of low-pressure pipe networks. The linear search and backtracking algorithm was used to control the iterative step size of the nonlinear equations of the pipe network nodes, which ensures the global convergence of the low-pressure hydraulic analysis of the pipe network.


---

## Body

## 1. Introduction

It is difficult to locate and simulate the leakage point and leakage form of the water supply pipeline network damaged by the earthquake. In order to meet the actual needs of the hydraulic function analysis of the post-earthquake water supply network, a more reasonable leakage point location model, leakage point area calculation model, leakage flow model, and hydraulic function analysis model must be used to analyze the pipeline network. In the calculation of the leakage area of the pipeline, Xing Yan calculated the displacement of the pipeline with the socket-type interface under the action of earthquake [1–6]. This method is only suitable for calculating the pipeline with the socket-type interface, such as when the allowable displacement is known, the ductile iron pipe and the cast iron pipe, but in the actual pipe network composition, there are pipes with different pipe materials and interface forms, such as the leakage analysis of welded steel pipes. For the consideration of the leakage area of the pipeline, it cannot be considered only from the point of the pull-off of the interface. According to the actual earthquake damage, Shi and O’Rourke considered the leakage forms of pipelines as leakage points with different shapes and gave the proportion of different leakage point forms in pipelines with different pipe materials and joint forms, as well as the different leakage point forms calculation method [7]. In terms of the point leakage model for calculating leakage flow, China Urban Water Supply Association gave the C model in 1988 [8]; the S model is a similar model to the C model. The form of the leakage point has not been considered in the C model, and from the value of its coefficient, the calculation result of the leakage flow is bound to be small. The S model considers that all leak points are equivalent round holes and experiments show that the shape of the outlet will directly affect the relationship between flow and pressure [9]. The shape of pipeline leaks under earthquake damage is irregular in most cases, so the calculation results of the S model will obviously be quite different from the actual situation. Under the earthquake damage, the water leakage of the water supply network increases, which leads to a decrease in the water supply pressure of the pipeline network, and the user’s water supply is insufficient. At this time, the user’s water distribution can be simulated by the pressure-dependent demand (PDD) relationship. In terms of low-pressure hydraulic analysis of the post-earthquake pipeline network, GIRAFFE and Gao Huiying both use the nodal fixed water demand model [10, 11]. When there is negative pressure at the node in the calculation result, the negative pressure node is processed and then the hydraulic analysis of the pipeline network is performed again until there is no negative pressure node. For a single failure condition, this analysis method requires multiple trial calculations to obtain reasonable results. Chen et al., Liu et al. used the PDD model for hydraulic calculation, but they did not point out the convergence problem in the hydraulic solution of the pipeline network when the PDD model was used [9, 12]. Most of the existing PDD models are piecewise functions, which will reduce the convergence performance of the iterative solution method for nonlinear hydraulic equations in the pipeline network [13, 14]. When the point leakage model and PDD model proposed in this paper are used to calculate the water leakage of the pipe section, the water demand function of the node will also be complicated. Therefore, it is necessary to apply a global convergence algorithm to solve the hydraulic equations of the pipe network.The innovation of this paper is by generating random numbers, the pipelines with earthquake damage are determined according to the earthquake damage rate. According to the proportion of different leakage point forms of pipelines with different pipe materials and interface forms, random sampling is used to determine the leakage point forms. The leakage area calculation model and the corresponding point leakage model summarized by Shi and O’Rourke based on the actual earthquake damage are used to analyze the leakage hydraulics of the pipeline network after the earthquake. The user node adopts the pressure-based water distribution (PDD)model simulation. The linear search and backtracking method is applied to the iterative solution of the nonlinear equations of the pipeline network node flow to ensure the global convergence of the solution of the low-pressure hydraulic analysis equations of the pipeline network. More than 5000 hydraulic simulations were carried out on the water supply network with different pipe materials under the earthquake intensity of 6 to 10 degrees, the probability mass function, and cumulative distribution function of the water supply rate were drawn, and the general law of the water supply rate changing with the seismic intensity was summarized.


## 2. Determination of Pipeline Leaks and Forms of Leaks

The probability of earthquake damage to the pipeline is(1)PL=1−e−λL.In the formula,L is the length of the pipeline, km, and λ is the seismic damage rate of the pipeline, points/km.Using the empirical formula of earthquake damage rate provided by the Japan Waterworks Association [15, 16],(2)λ=Cp×Cd×Cl×Sdp,Sdp=6.33×10−5×v2.10.In the formula,CP, Cd, and Cl are the influence coefficients of pipe material, pipe diameter, and liquefaction, respectively. v is the peak seismic velocity, cm/s; and Sdp is the standard seismic damage rate. Refer [15] for the value of the influence coefficient.The seismic damage probabilityPL of each pipeline is obtained according to the parameters of each pipeline, and then a random number equal to the number of pipelines in the pipeline network is generated, and a random number r is assigned to each pipeline. When 0<r<PL, the pipeline is damaged. In the case of a short pipeline (not more than 1 km), it is assumed that the damage point is at the midpoint of the length of the pipeline.Shi and O’Rourke divided leakage into five forms [7, 17] and gave the calculation method of leakage area for each leakage form, which is suitable for all kinds of socket-type interface pipes and continuous non-interface pipes, as shown in Table 1. The proportion of different leakage forms in pipelines with different materials and interface forms is given. In the pipelines that have been damaged in the pipeline network, the leakage forms of leakage points are determined by random sampling.Table 1
Pipe leakage forms and their calculation formula (according to document [7]).
Leak formCircumferential looseningLateral cracksLongitudinal cracksDamaged pipe wallWall tearLeakage area calculation formulaAL=t⋅k⋅D⋅πAL=0.5⋅θ1⋅D2⋅πAL=L⋅D⋅θ2AL=k1⋅k2⋅D2⋅πAL=k⋅D⋅W⋅πParameter descriptionD is the diameter of the pipe (mm)θ1 is the opening angle of the circular crack (°), which is taken as 0.5L is the longitudinal crack length (mm), take 13k is a dimensionless constant, taking 0.3The opening angle ofθ2 longitudinal crack (°), take 0.1W is the tearing width of the pipe wall (mm), take 12t is the thickness of annular separation (mm), take 10k1, k2 are dimensionless constants, take 0.05k is a dimensionless constant, taking 0.3Schematic diagram


## 3. Leakage Point Form and Its Corresponding Point Leakage Model

In the post-earthquake hydraulic analysis of the pipeline network, the seepage node where the pipeline is damaged is usually assumed to be a new user node, and the seepage flow is calculated using a point leakage model:(3)QL=C⋅μ⋅AL⋅H,QL is leakage flow m3/s, AL is leakage area m2, and H is leakage point water pressure (m).Four leakage models corresponding to the seepage forms listed in Table1 are selected for the post-earthquake hydraulic analysis of the pipeline network in the example in this paper, as shown in Table 2.Table 2
Different damage forms and their corresponding leakage model.
Model name/sourceLeakage modelPipe damage formThe corresponding leakage point form in this exampleC ModelQL=0.421⋅AL⋅HNo special restrictionsCircumferential looseningTabesh et al. [18]QL=2.418⋅AL⋅HCrackLateral cracksAWWA(1999) [19]QL=0.658⋅AL⋅HCracks or interface failureLongitudinal cracksTabesh et al. [18]QL=3.486⋅AL⋅HOrificeDamaged pipe wall


## 4. Pressure-Dependent Nodal Demand (PDD) Model

In the traditional hydraulic analysis of the water distribution network, it is assumed that the pressure of the pipe network meets the user’s requirements; that is, the user’s water demand is fully met, and the water distribution of the pipe network nodes is fixed [20–25]. Under the action of earthquake, the pressure of the pipeline network will be too low to provide enough water demand due to the leakage of the pipe section and the failure of the pumping station. The above assumptions do not hold. The nodal water distribution at this time is related to the nodal pressure. If the fixed water distribution is still used instead of the actual water distribution for the hydraulic calculation of the pipe network, the iterative process of the hydraulic equations is based on the assumption of the continuity of the node flow. At this time, solving the result that satisfies the fixed water demand of the node will reduce the total water head of the node. At the same time, the node negative pressure that does not conform to the actual situation occurs.Since the nodal flow model (PDD) based on water pressure is adopted, the water demandQi of node i is determined by the total water head Hi at this point (QiHi), and its general expression is as follows [26, 27]:(4)Qiavl=0,Hiavl<Himin,Qireq⋅DSRHiavl,Himin≤Hiavl<Hireq,Qireq,Hireq≤Hiavl.In the formula,Qiavl is the actual water distribution of the node, m3/s; Qireq is the node water demand when the pressure is sufficient, m3/s; DSRHiavl is the node water supply satisfaction rate; Hireq is the total water head required to meet the water volume of the node, m; Hiavl is the actual total water head of the node, m; and Himin is the limit total head of the node with no water and partial water supply, m.


## 5. Low-Pressure Hydraulic Simulation Algorithm of Post-Earthquake Pipeline Network Based on Nodal Flow Continuity Equation

The nodal flow continuity equations of the pipeline network are expressed as the function of nodal water pressure:(5)FH=A⋅Qp−QN=0,where H is the nodal water pressure vector (m); H=H1,H2,…,Hn1 (n1 is the nodes’ quantities of original pipeline network). Use Newton–Raphson iteration to carry out the hydraulic analysis of the pipeline network; the hydraulic equation of the pipeline network is a nonlinear equation system about the water pressure {H}. A is the correlation matrix of the pipeline network structure. QPH is the pipeline flow vector (L/s): QP=QP1,QP2,…,QPm (m is the quantities of the pipes). QNH is the nodal water demand vector(L/s). FH=f1,f2,…,fn1 is the hydraulic imbalance vector of the nodes.Due to the virtual leakage nodesQL1,QL2,…,QLn2 (n2 is the quantities of leakage nodes) in the leakage state, the original balanced hydraulic equation appears unbalanced, that is, FH0≠0. The hydraulic equation of the leaking network is rebalanced.(6)FH=A⋅QP−QN−QL=0.The Newton–Raphson iteration can make equationF(H) approach 0. Set the accuracy and calculate the solution to the hydraulic equation of the leaking network.(1)
Given the initial nodal water pressureH0=H10,H20,…,Hn0, calculate FH0. If each element in FH0 is less than the calculation accuracy, then H0 are the desired water pressure; otherwise, set k=0 and go to step (2).(2)
Calculate the value of∂F/∂Hk at H=Hk to get the Jacobian matrix Jk.(3)
Solve for the corrected water headΔHk=−Jk−1⋅FHk for each step in the Newton–Raphson iteration and compute the water pressure approximation values Hk+1 for the next step.(4)
CalculateFHk+1, if each element in FHk+1 is less than the calculation accuracy; then, Hk+1 are the required water pressure; otherwise, set k=k+1 and go to step (2).Adding the pressure-determined flow relationship to the Jacobian matrixJ increases the nonlinearity of equation (6), resulting in a decrease in its solution convergence performance. The Newton iteration step size ΔHk directly updates the independent variable Hk+1=Hk+ΔHk. If the initial value is not close enough to the actual equation root, the iteration will cause the result to deviate to an irregular distance, resulting in nonconvergence. A strategy for determining a reasonable iterative step size is that the step size reduces F2=FTF. Assuming f=1/2FT⋅F, the above problem is transformed into a f-minimization problem. Note that the Newton step size ΔH is the descending direction of f [28]:(7)∇f⋅δH=F⋅J⋅−J−1⋅F=−FT⋅F<0.This ensures that a step size adjustment factorε (0 <ε ≤ 1.0) makes fHk+1 drop sufficiently, where Hk+1=Hk+ε⋅δH.(8)fHk+1≤fHk+α⋅∇f⋅Hk+1−Hk.To speed up the linear search, another linear search and backtracking method is used here:(9)gε=fHk+ε⋅δH.In order to solve the minimum value ofgε (that is, fHk+1), it is necessary to solve the value of ε such that g′ε=∇f⋅δH=0.(10)gε≈g1−g0−g′0ε2+g′0ε+g0,(11)ε2=−g′02g1−g0−g′0,(12)gε≈aε3+bε2+g′0ε+g0,(13)ab=1ε1−ε21ε12−1ε22−ε2ε12ε1ε22⋅gε1−g′0ε1−g0gε2−g′0ε2−g0,(14)ε3=−b+b2−3ag′03a.In thei-th backtracking process after the second time, εi, εi−1, and εi−2 can be used to replace ε3, ε2, and ε1 in Equation (13) and (14), respectively. Programming the improved hydraulic analysis algorithm with MATLAB software and the algorithm flow is shown in Figure 1.Figure 1
Flow chart of the modified post-earthquake hydraulic analysis simulation algorithm.The post-earthquake water supply rateCs can be expressed as follows:(15)Cs=∑i=1niQi∑i=1ni0Qi.In the above formula,Cs is the ratio of the total water supply immediately after the earthquake to the total water supply before the earthquake; Qi is the water demand of node i; and ni0 and ni are the number of pre- and post-earthquake nodes to meet water supply demand.Assuming that the water flow demand of nodei is Qi, if node i finds m times of normal water supply in n times of Monte Carlo simulation, the reliability of node i can be expressed as the probability as follows:(16)PQi=mn.


## 6. Examples

This paper takes the water supply pipeline network (Apulian) of the industrial area of Apulia in southern Italy as the research object. For the basic data of the pipe network, see reference [29], and the schematic diagram of the model is shown in Figure 1. All algorithms used in this study were compiled by the authors using MATLAB software.Using the leakage model and the corresponding point leakage model introduced in Parts 2 and 3 and the pressure-dependent nodal demand model and nodal flow continuity equation using linear search and backtracking introduced in Parts 4 and 5, the leakage hydraulic simulation is carried out on the water supply network shown in Figure2. The Monte Carlo simulation program was used to simulate the two kinds of pipes, cast iron pipe, and steel pipe, respectively. The calculation results are shown in Table 3.Figure 2
Schematic diagram of Apulian water supply network (spatial distribution).Table 3
Evaluation results of water supply rate of water supply network using different pipes.
IntensityLiquefactionMaterialSamples numberCs meanCs standard deviationRisk valueVIINoCast iron pipe66280.9850.0590.015VIIINoCast iron pipe55480.9610.1110.039IXintermediateCast iron pipe55970.8060.0940.194XSeverityCast iron pipe64830.2960.0860.704IXintermediateSteel pipe53140.9070.0920.093From Figures3(a) to 7(a), it can be seen that the peak frequency distribution of the post-earthquake water supply rate Cs of the Apulian water supply network moves to the direction of water supply rate decreasing with the increase in the intensity. At VII degree and VIII degree, the peak frequency appears above the water supply rate of 0.9 (that is, the section with a large probability mass function), the peak frequency occurs at the water supply rate of 0.8 at IX degree, and the section where the peak frequency of the water supply rate of X degree occurred is 0.2-0.3. And with the increase in intensity, the peak frequency decreases, and the distribution range of water supply rate increases. This reflects the distribution of the simulated water supply rates in the sample. As the intensity increases, the value of the distribution range of the water supply rate decreases, and the frequency of lower water supply rates increases significantly, which will lead to a decrease in the water supply rate of the network as the intensity increases [30, 31]. At the same time, the change of the slope of the curves in Figures 3(b) to 7(b) and the distribution of the straight line segments also illustrate this trend. For example, when the intensity is lower, the slope of the curve in the range of the water supply rate less than 0.8 is smaller, indicating that the density of the samples distributed in the smaller water supply rate segment is lower and the slope of the curve in the range of larger water supply rate (above 0.8) increases sharply, indicating that the sample density distributed in the higher water supply rate section is high, the water supply rate of the network is at a high level, and the earthquake disaster risk is low. When the intensity is high, the slope of the curve increases sharply in the section where the water supply rate is small, indicating that the sample of water supply rate is mainly distributed in a lower value range, water supply rate of the water supply network is low, and the earthquake disaster risk is high.Figure 3
Post-earthquake water supply rate (cast iron pipe) with no liquefaction of intensity VII. (a) Probability mass function. (b) Cumulative distribution function.
(a)(b)Figure 4
Post-earthquake water supply rate (cast iron pipe) with no liquefaction of intensity VIII. (a) Probability mass function. (b) Cumulative distribution function.
(a)(b)Figure 5
Post-earthquake water supply rate (cast iron pipe) with intermediate liquefaction of intensity IX. (a) Probability mass function. (b) Cumulative distribution function.
(a)(b)Figure 6
Post-earthquake water supply rate (cast iron pipe) with severity liquefaction of intensity X. (a) Probability mass function. (b) Cumulative distribution function.
(a)(b)Figure 7
Post-earthquake water supply rate (steel pipe) with intermediate liquefaction of intensity IX.
(a)(b)In this paper, the improved hydraulic analysis algorithm, the original algorithm (the method that hardly uses any improved low pressure hydraulic calculation method), and “two-step iteration” low pressure hydraulic analysis method [32] are used to perform hydraulic analysis on the post-earthquake seepage pipe network. In this paper, the Wagner model is selected as the PDD model of the user node water distribution. It is generally believed that the continuity of the pressure-based water distribution function curve will affect the convergence performance of the iterative calculation. The results of the three analysis algorithms are shown in Figure 8. It can be seen that both the original method and “two-step iteration” method decline rapidly in the early stage, and the decline ability degenerates when it is gradually approaching the true solution in step 7. Under the same initial value of iteration and allowable convergence error, the method in this paper iterates for 11 steps, while the original method iterates for 29 steps and “two-step iteration” method iterates for 20 steps. Also, the method in this paper takes less time. The linear search and backtracking method proposed in this paper has the characteristics of small computation and fast convergence.Figure 8
(a) Comparison between the improved hydraulic analysis algorithm in this paper and the original method. (b) Comparison between the improved hydraulic analysis method in this paper and “two-step iteration” low pressure hydraulic analysis method.
(a)(b)Based on Monte Carlo simulation, the spatial distribution ofPQi of 24 user nodes in the Apulian water supply network under the action of degree IX earthquake is shown in Figure 9.Figure 9
Spatial distribution of node reliability of Apulian water supply pipe network (cast iron pipe) after moderate liquefaction earthquake with IX intensity.The above reliability analysis results show that when the Apulian water supply pipe network adopts cast iron pipes and suffers from a degree IX earthquake, the number of nodes with a reliability of 85% of the user nodes accounts for 70% of the summary points of the pipe network. The water supply network still has water supply capacity after the earthquake.


## 7. Conclusion

(1)
On the basis of applying Shi and O’Rourke’s leakage point model and its corresponding point leakage model for the hydraulic analysis of the post-earthquake pipeline network, a global convergence algorithm is proposed to simulate the low pressure hydraulic of the water supply network. The convergence problem of Newton’s iterative algorithm in network leakage analysis is improved, so that the calculation results of hydraulic leakage rate of pipe network and nodal pressure are more accurate and reliable.(2)
Taking a water supply network (Apulian) in Italy as an example, the Monte Carlo simulation results of different pipeline materials encountering different intensities are compared, and the reliability results of the user nodes of the water supply network are given.



---

*Source: 8915561-2022-06-17.xml*