Derivation and Validation of a Leakage Model for Longitudinal Slits in Polyethylene Pipes

: Although we do not understand the complex behavior of leaks in thick-walled plastic pipes, the ability to model them is essential to the assessment and mitigation of the real losses from drinking water – distribution systems. A methodology was followed whereby numerical and physical experimental data were used to develop a generalized model that quantified both the structural and leak dynamics of longitudinal slits in thick-walled viscoelastic pipes. The dimensionally homogeneous time- and pressure-dependent leakage model was validated using experimental data independent from the physical observations used for the viscoelastic calibration. This was shown to be an effective tool to predict the distinct short- and long-term responses. Such accurate models, which describe the characteristic response of complex leaks, are crucial in improving the understanding of leakage


Introduction
Real losses from water-distribution systems (WDS) that result from background leakage and bursts have a significant effect on the overall sustainability of this vital infrastructure, both economically and environmentally. Leakage statistics from 2013 showed that the level of leakage in the UK stood at 172 ML=day (Ofwat 2013). This value has remained approximately static for over a decade. The ability to assess, analyze, and mitigate the impact of such losses is dependent on the fundamental understanding of the behavior of leaks and the capability to predict their behavior. Studies such as Greyvenstein and van Zyl (2006) have emphasized the need to consider the sensitivity of leakage to pressure due to the dynamic nature of the leak areas.
Plastic pipes are a popular choice for hydraulic pipelines due to their cost effectiveness. Leaks occurring in this type of pipe, particularly the longitudinal cracks/slits that are a dominant failure mode in plastic pipes, result in a complex leakage behavior due to the inherent material rheology. Quantifying the pressure-dependent structural response of the leaks in these viscoelastic pipes and the interdependence with the leak hydraulics is crucial in predicting the time-and pressure-dependent leakage and allowing accurate leak management strategies.

Background
Leaks occur in a diverse range of shapes and sizes, depending on factors that include the pipe material, loading, ground conditions, age, and manufacturing process. Idealized leakage flow rates can be estimated using the orifice equation, which assumes a constant leak area and a square-root relationship between the pressure and the flow rate. In May (1994) it was surmised that leaks could be more accurately modeled by assuming that the leak area varied under different internal pressures and that the leak flow rate could be described by a fixed and variable area discharge (FAVAD) model where Q = leakage flow rate; h = difference in the pressure head between the inside and outside of the pipe (the driving head across the leak orifice); C d = coefficient of discharge; A 0 is the initial leak area and dAðhÞ is a function that determines the rate of change of area with internal head; and g = gravitational acceleration. Experimental studies have subsequently confirmed that the leak area is the primary causative factor of this observed behavior [e.g., Cassa and van Zyl (2008); Ferrante (2012)]. The importance of understanding the dynamic nature of leaks has resulted in a number of studies investigating the function that determines the rate of change of the area with the internal pressure head. In particular, studies have explored how the area of the leak is affected by the geometry of the system, the material properties of the pipe, the loading conditions, and time dAðhÞ¼KðGeometry; Material; Loading; TimeÞð 2Þ where dAðhÞ should be considered to be an unknown function, K, of the geometry, material, loading, and time. Much of this work has focused on linear elastic materials, considering longitudinal cracks (Grebner and Strathmeier 1984;Bhandari and Leroux 1993;Ávila Rangel and Gonzalez Barreto 2006;Al-khomairi 2005;De Miranda et al. 2012). For clarity, herein cracks refers to the naturally occurring leak apertures in pipes and slits refers to the artificially manufactured failure apertures.

Numerical Studies
Numerical simulations offer a powerful tool to explore the structural behavior of leaks. Cassa and van Zyl (2011) used finiteelement analysis (FEA) to investigate the significance of various parameters in the leak behavior of longitudinal, circumferential, and spiral slits. The results compiled were used to statistically derive an equation defining the pressure-dependent change of area for each leak type and were subsequently input to a reformulation of the FAVAD model defined by May (1994). The finalized leakage model presented by van Zyl and Cassa (2014) provides an effective and simple methodology to quantify the pressure-dependent leak area and hence to predict an explicit leakage exponent value. However, the dimensionally inconsistent nature of the derived expression means that little physical meaning can be attached to the components of the predictive model.
Accurate definitions of the boundary conditions are critical to the validity of the numerical simulations, including both the fixed boundary conditions defining the degree of freedom of a given model and the variable boundary conditions (e.g., the applied loadings). Methods to simulate the pressure loading due to hydraulic conditions are well practiced (Rahman et al. 1998;Cassa and van Zyl 2008;De Miranda et al. 2012). However, there is still uncertainty regarding the magnitude and hence the significance of the slit face loading for a given leak. Such localized pressures are commonly excluded from analyses because of the uncertainty in the definition of the true pressure distribution (Takahashi 2002).

Empirical Investigations
Experimental data from real and artificial leaks have been used to verify and validate the application of leakage models (De Miranda et al. 2012;Franchini and Lanza 2014). Buckley (2007) explored the dynamic leak area of a range of different failure types including circular orifices and longitudinal slits, using a hydraulic pressurized bladder. The work provided an insight into the influence of different parameters on the leak area, in particular, the significance of the slit length for the scale of deformation. However, the work fitted linear pressure-area relationships and therefore ignored the viscoelastic behavior of PVC pipes, observable in the results.

Polyethylene Pipes
Polyethylene (PE) pipes are frequently used in the water industry due to the cost benefits offered by the inherent durability and flexibility of the material (GPS PE Pipe Systems 2014). In contrast to the typical steel, PVC, and some cast-iron pipes where the diameter-to-wall thickness ratio (D=s) is commonly greater than 20, classifying such pipes as thin-walled, PE pipes within the diameter range 20-1,000 mm are classed as thick-walled (BSI 2011). The spatially dependent cross-sectional stress distributions and the viscoelastic nature of PE result in the complex behavior of failures, primarily due to the time, temperature, and pressure dependence of the structural response. Additionally, the extrusion-based manufacturing process tends to result in the formation of residual stresses within the material due to the differential cooling of the internal and external faces (Hutar et al. 2012). The solidified external surface and the slowly cooling interior create a through-wall stress distribution. The average measured residual stresses may be simply approximated as a tensile stress on the external face and a compressive stress on the internal face in the range 2.0-5.0 MPa, with a nonlinear distribution across the wall thickness (Guan and Boot 2004;Frank et al. 2009).
The residual stresses are important in quantifying the material behavior under hydraulic loading conditions and, more critically, in quantifying the structural behavior of pipe failures such as longitudinal slits. Massari et al. (2012) presented the results from a series of physical observations of the behavior of longitudinal slits in highdensity polyethylene (HDPE), emphasizing the time and pressure dependence of the localized strain around the leak. The article also demonstrated that the longitudinal slits, a common failure mode in the direction of the pipe extrusion (O'Connor 2011), demonstrate the complex leakage behavior (Massari et al. 2012). Ssozi et al. (2015) produced FEA models of the viscoelastic response of longitudinal slits and demonstrated the main impacts on slit openings in plastic pipes, namely, the dependence on the time after loading to determine the instantaneous leak area. However, there remains a need to determine all the key parameters controlling the time-dependent leakage behavior of such longitudinal slits in viscoelastic pipes, equivalent to the studies of linear elastic pipes.

Research Aim
The aim of this research was to derive a simple, dimensionally homogeneous model to quantify the dynamic leak area of stable longitudinal slits in pressurized medium-density polyethylene (MDPE) pipes. The research sought to investigate the influence of the geometrical characteristics, material properties, and internal pressure loading conditions, along with additional, more complex loading conditions. Ultimately, the research intended to validate an analytical model, integrating this pressure-dependent leak-area expression into an expression that describes the leakage behavior of this important leak type.

Research Methodology
The research methodology was to develop a generalized analytical model to describe the leakage behavior of stable longitudinal slits, built on a combination of numerical simulations and experimental data. The data presented here also appeared as part of the data used and analyzed in Fox (2015). The steps taken in this study were 1. numerical simulations of the pressure-dependent leak area in linear elastic pipes; 2. exploration of complex loading cases; 3. derivation of a simple linear elastic pressure-dependent leak area model from the simulation results; 4. calibration of the leak area model for linear viscoelastic pipe materials; 5. incorporation of the leak area model into the FAVAD equation; and 6. validation of the dynamic leakage model with independent flow measurements. Finite-element analyses were used to empirically derive an expression based on Eq. (2) defining the pressure-dependent variable leak area of the longitudinal slits in thick-walled linear elastic pipes. In addition to the standard internal pressurization previously studied, the numerical simulations explored the effect of the slit face loading, the residual stress due to the manufacturing process, and the longitudinal pipe stresses. A simple and dimensionally homogeneous form was targeted in the analysis to provide a functional and efficient model. The experimental data of the synchronous pressure head and leak area data for the longitudinal slits in medium-density polyethylene (MDPE) were subsequently used to calibrate the time-dependent viscoelastic modulus in the variable leak area model. The leak area model was then incorporated into Eq. (1), and it was validated against the further independent experimental leak flow data.

Numerical Simulations of the Pressure-Dependent Leak Area in Linear Elastic Pipes: Model and Boundary Conditions
In order to quantify the significance and the influence of a range of parameters on the structural behavior of longitudinal slits in a pressurized thick-walled pipe in an efficient and reliable manner, a three-dimensional (3D) finite-element model was created. The methodology explored a range of key parameters, which were categorized as (1) geometric configurations, (2) material properties, and (3) loading conditions, comparable to the work conducted by Cassa and van Zyl (2008) but regarding thick-walled pipes. Table 1 provides the details of the parameters investigated.
A 3D linear elastic FEA model was created, with a nominal internal pipe diameter of 50 mm, a wall thickness of 6.5 mm, and longitudinal slit dimensions of 60 mm by 1 mm (SDR11 PN16). This and all the derivatives of it were developed in the finiteelement software APDL Mechanical (ANSYS Inc. 2015).
The ends of the pipe were modeled as fully constrained and the pipe length modeled as long enough that these boundary conditions did not affect the leak behavior, which was effectively unconstrained. Simulations with different end conditions and pipe lengths were run to select the boundary conditions that had insignificant effects on the dynamic leak areas relative to the localized boundary conditions. The fixed boundary conditions were used to replicate the physical boundary conditions used in the Fox et al. (2016) laboratory investigations, in order to achieve an accurate comparison and subsequent calibration of the structural behavior between the two sets of empirical data. A plane of symmetry was used to reduce the model size and consequently the required simulation time, achieved by splitting the model in half along the axial length of the pipe (z-direction) through the center of the longitudinal slit. Full-and half-pipe models were run and compared in order to confirm that this had a negligible impact on the observed dynamic leak area. The model was then fixed in the x-direction along the line of symmetry, as shown in Fig. 1.
The internal pipe pressure was simulated by applying a nodal pressure loading. The range of simulated pressures investigated were representative of the typical operating pressures in waterdistribution pipes in the UK. Pressure loading steps equivalent to the pressure heads of 10, 20, 30, 40, 50, and 60 m were implemented.
The influence of the residual stresses inherent in polyethylene pipes, neglecting time-dependent effects, were also explored using the linear elastic FEA model. The approximate residual stress distributions quantified in existing studies (Guan and Boot 2004;Frank et al. 2009), where the average measured stresses were 4 MPa (tension) and −4 MPa (contraction) at the internal and external faces of the pipe, respectively, were adopted for the investigation.
A standardized meshing scheme was used in order to maintain a consistent level of accuracy and comparability among all the models developed, using 20-node 3D SOLID186 elements, which have the capability to accurately model large deflections and strains and effectively simulate the residual pipe stresses (ANSYS Inc. 2015). A refined 3D tetrahedral-shaped element mesh focused around the leak opening using line sizing, with a minimum resolution of 2 mm (100 nodes along the slit edge), was implemented. A gradient mesh was then created adjacent to the leak, increasing in coarseness toward the pipe ends. Fig. 2 shows the final standardized mesh for the platform model.
The mesh sizing provided an efficient and accurate model to simulate the structural behavior of the specific leak type. A mesh invariance analysis (Fig. 3) showed that increasing the fineness/ resolution of the developed mesh did not significantly alter the simulation solution, thus validating the developed platform model. The developed mesh size was comparable to those used in the studies presented by van Zyl (2008) andDe Miranda et al. (2012).
Each variable parameter (provided in Table 1) was then changed independently to explore its effect on the dynamic leak area; a total of 755 simulations were run. Fig. 1 shows the structure of the pipe model used in the simulations, including the applied fixed boundary conditions.

Exploration of Complex Loading Cases
In addition to the simple internal pressurization loading, this work analyzed three complex loading conditions. The first complex loading condition considered included the presence of a slit face loading, the second included the impact of residual manufacturing stresses in the pipe walls, and the third included an applied longitudinal stress.

Slit Face Loading
As the fluid exits the pipe through the slit, it applies pressure to the slit face; this additional load causes an increased opening of the slit compared with the internal pressure alone. Slit face loading varies as a function of the internal and external pressure and the external pipe burial conditions. The three load cases were simulated as potential representations of the real slit face loading based on , where s n is the discretized wall thickness over n steps). Fig. 4 shows an example from the results of the analysis, using a 20 × 1 mm slit, where the increases in the central deflection of Load Cases B and C relative to Load Case A were 12.25 and 7.28%, respectively. The results indicated that the slit face loading was significant with respect to the change of the leak area. However, the explicit relationship between this pressure distribution and the pipe geometry and the external boundary conditions were beyond the scope of the investigation. This parameter was therefore not included in any further analysis aimed at deriving an analytical leak area model, for which a quantification of the actual slit face pressure distribution was required.

Residual Stress
The presence of residual stresses in the pipe walls was simulated using the residual stress distribution previously described. The leak    areas were quantified and compared with models with no residual stress, for a range of pressures. Fig. 5 shows a selection of results from the analyses for the 40, 60, and 80 × 1 mm slit models. The figure shows that the residual stresses resulted in a relative offset (reduction) of the measured leak area, which was approximately constant for each case. In other words, the residual stress altered the initial area, A 0 , but the change of area, dA, remained constant for simulations with and without an applied residual stress. This result indicated that the effect of the residual stresses may be directly integrated in a leak area model by the inclusion of the initial area, which is a function of the inherent residual stress distribution. Consequently, the residual stress parameter was not included in the derivation of the leak area model.

Longitudinal Stress
As anticipated, the longitudinal stresses were shown to have a negligible influence on the dynamic leak area. This corresponded with the findings of Cassa and van Zyl (2008). The longitudinal stresses were therefore not included in the further analysis.

Derivation of a Simple Linear Elastic Pressure-Dependent Leak Area Model from Simulation Results
The objective of the analysis was to develop a dimensionally consistent analytical model, in a simple form, that included the parameters from the initial FEA that were observed to have an appreciable effect, in order to provide an efficient and practical tool for the assessment of the pressure-dependent area of the leaks, of the form of Eq. (2). In order to assess the impact of each parameter on the structural behavior, a power term-based analytical expression was evaluated using statistical multiple regression analysis. The power term formulation was used for the analysis, akin to the work conducted by Cassa and van Zyl (2011). An initial analysis of the individual parameters suggested that the behavior of the leak area could not be represented by a single input parameter but rather was the coupled effect of all the significant parameters. This was evidenced through the varying power relationships observed when individual component analyses were undertaken, as shown, for example, in Fig. 5, where the apparent relationship between P and the change in the leak area varied as the leak length changed. Table 2 provides the results of the multiple regression fitting process.
The regression fit produced Eq. (3), describing the change of the leak area with an associated R 2 value of 0.894; C 0 is a constant coefficient equal to 1.87 where the symbols are as defined in Table 2. A research aim was to achieve a dimensional homogeneous equation to quantify the physical significance of the independent and coupled parameters.
Despite the quality of the fit of the mathematical expression, Eq. (3) does not meet this requirement. To achieve dimensional homogeneity, the lower-order terms (the width and Poisson ratio) were neglected. The remaining terms were then refitted and the powers rounded to the nearest whole number. The simplest expression evaluated that captured the modeled behavior of the change of area is shown in Eq. (4), demonstrating that the slit length (L c ) and pipe wall thickness (s) are the most important geometrical parameters In this equation, the factor P=E represents the relative stiffness of the structure with respect to the loading. From the fundamental theory defining the hoop stresses in thick-walled cylinders, the linear inverse relationship between the pressure and the elastic modulus was noted, whereby a reduction in E was equivalent to an increase in P.
In Eq. (4), the term C 1 represents a dimensionless constant that was simultaneously evaluated with the parameter exponents; physically, it describes all the dynamics of the slit expansion not explained by the exponent terms. Further understanding of C 1 was gained by exploring its relationship to a range of dimensionless parameters, based on reasonable engineering judgment. The ratio of crack length (L c ) and pipe circumference (πD) is a dimensionless parameter commonly used in structural mechanicstodefinetherelativesizeofacrackinapipeorpressurized vessel, primarily for circumferential cracks. Employing this term and plotting the relationship with the coefficient term C 1 , as shown in Fig. 6, highlights the capacity of L c =πD as a predictor of C 1 .E q .( 5) was fitted to the data shown in Fig. 6 for integration into the final dimensionally homogeneous formulation   Fig. 6 shows that this relationship is representative of all pipes tested with a slit length-to-pipe circumference ratio less than 1. When the slit length increases greatly, the representation is less valid. In reality, such large-scale slits may result in the total failure of the structural integrity of the pipe, which means that the elastic deformation of the slit is inconsequential. When Eq. (6)w a s substituted into Eq. (4), it was shown to provide a very good fit to all the collated data with a mean Area model =Area measured ratio of 1.01 and a standard deviation of 0.12 (excluding the data for a slit length-to-pipe circumference ratio greater than 1). Compared with the regression fitted model with a constant scaling coefficient [Eq. (3)], the dynamic leak area model presented in Eq. (4), including the dimensionless coefficient C 1 , produced an R 2 value of 0.969, highlighting the improved accuracy of this simple model form.

Calibration of the Leak Area Model for Linear Viscoelastic Pipe Materials
Considering the fundamental structural behavior of polyethylene pipes, a constant elastic modulus may be replaced by an empirically calibrated time-dependent elastic modulus to capture the pressuredependent linear viscoelastic structural response. This corresponds with the constitutive linear viscoelastic equations that define the relationship between stress and strain using a Volterra integral equation, where time is the variable limit of integration (Ssozi et al. 2015).
Eq. (4) is a predictor of the change of area for longitudinal slits in linear elastic thick-walled pipes under applied hydraulic pressurized loading. Substituting E for Eðt; TÞ, a time (t) and temperature (T) dependent elastic modulus that is equal to the reciprocal of creep compliance, JðtÞ, enables the definition of the leak area of a longitudinal slit in a viscoelastic thick-walled pipe The time-and temperature-dependent elastic modulus is the multiplicative inverse of the generalized Kelvin-Voigt creep compliance model, with the instantaneous elastic modulus component (E inst ) accounted for by the empirically derived temperaturedependent formula presented by Bilgin et al. (2008). The chosen mathematical representation of the viscoelastic behavior has been previously demonstrated to be an effective model of polyethylene material rheology (Covas et al. 2004;Ferrante et al. 2010). The retardation time components, τ n , were predetermined before the fitting process to cover the relative short-and long-term behavior of the material; this approach was previously validated in Fox et al. (2016). The time-and temperature-dependent elastic modulus adopted for the linear viscoelastic calibration is Experimental data from a single MDPE test section, with a 50-mm internal diameter and 6.5-mm wall thickness (d=s ¼ 9.69)an dc o ntaining a 60 × 1 mm slit at three discrete pressure heads (10, 20, and 25 m), were used for the calibration of the viscoelastic constants, using the experimental setup fully described in Fox et al. (2016). A recirculating pipe loop was used, consisting of a 141-m length of pipe of the same specification as the test section, fed by a variable-speed pump. A removable section 62 m downstream of the pump allowed the installation of the test section in the main pipe loop. A single valve downstream of the test section was closed so that the total system flow rate was equal to the leak flow rate through the longitudinal slit. Synchronous measurements of the leak flow rate, pressure, and leak area were recorded under quasi-steady-state conditions (slowly changing) in the controlled laboratory environment. The leak area was measured using a nonintrusive imaging technique, in which images were recorded for the period of pressurization and depressurization and the leak area was calculated using a calibrated binary pixel count.
A nonlinear least-squares methodology was used, using the Levenberg-Marquardt algorithm and a function tolerance of 1 × 10 −12 , to fit the creep compliance terms (J n ), using the experimental data of the leak area, pressure head, test section parameters, mean daily temperature, and Eq. (6). As directly measured leak areas were available, the calibration of the viscoelastic model did not require fitting of the coefficient of discharge. Table 3 provides the results and standard errors of the regression for the calibration process. The increase in the standard error of regression with the increasing pressure was primarily a function of the increase in the uncertainty of the experimentally measured leak areas (Fox et al. 2016).
The generalized Kelvin-Voigt viscoelastic model [Eq. (7)] provided sufficient detail to account for the relative short-term and long-term structural responses using the predefined time components of τ ; it also accounted for the observed hysteresis, the difference between the creep and recovery phases. The experimental data indicated that while the 5-day material behavior did not form a constant hysteresis cycle, it was tending toward this pseudo-equilibrium state with time. That is, further repeat loading cycles beyond the 5-day limit presented would have resulted in a consistent daily material response. The results therefore indicated that the modeled leak area provided a suitable fit to the Fig. 6. Coefficient (C 1 ) analysis from finite-element data.
experimental data with greater error associated with the first and second pressurization phases.

Incorporation of the Leak Area Model into the FAVAD Equation
To validate the leak area model, including the creep compliance calibration terms, the model was integrated into the FAVAD model, as shown in Eq. (1), where the pressure term in Eq. (6), P, has been converted to the pressure head, hðtÞ where ρ = fluid density. This equation describes the leakage flow rate for longitudinal slits over time due to changes in the internal pressure head as a function of the initial leak area, slit length, pipe diameter, pipe wall thickness, coefficient of discharge, and viscoelastic properties.

Validation of the Dynamic Leakage Model with Independent Flow Measurements
The new relationship [Eq. (8)] was then validated using additional independent physical observations of the leak flow and pressure head in longitudinal slits in MDPE pipes. The leakage flow rates were recorded during 3-to 5-day quasi-steady-state tests, which involved repeated cycles of an 8-h pressurization phase followed by a 16-h depressurization phase. These tests were conducted using the experimental facility described in Fox et al. (2016). To model the leakage flow rate in Eq. (8), a constant discharge coefficient (C d ¼ 0.6) was assumed. In Fox et al. (2016), the C d values for the slits were shown to vary with different manufactured leaks.
To demonstrate the effect of uncertainty in C d , models were also produced with C d AE 5%. Fig. 7 shows the measured and modeled leak flow rates for the 60 × 1 mm slit.
The results generally showed a very good correlation with the experimental results. The results indicated that the correlation between the measured and modeled leak flow rates increased with repeated pressurizations. This was a direct result of the calibrated viscoelastic model, which as fitted here favors the long-term response of the leaks. In reality, leaks rarely exist in previously unloaded pipe sections; therefore the effectiveness of the model in predicting the leakage behavior under truly representative distribution-system conditions was not undermined.
Step changes in the measured leakage during each test phase (not including the initial pressurization) were surmised to be the result of the sensitivity of the slits to blockages from system debris. This was observed to be most significant for the 10-m pressure head tests where there was less driving force to cause the expulsion of any debris that reduced the cross-sectional leak area. It is also noteworthy that the pressure heads of 10 m or less were extremely conservative in comparison with real system pressures, due to the minimum service pressure requirements, outlined by Ofwat (2008), to ensure that a minimum pressure of 7.14 m (≈0.7 bar) is maintained at consumer taps.
To confirm the validity of the generalized leakage model, the leakage flow rates from two supplementary test sections were modeled. The test sections contained 20 × 1 mm and 40 × 1 mm artificially manufactured longitudinal slits in a pipe of the same specification as the 60 × 1 mm test section previously described. Figs. 8 and 9 show the results of the model validation, for which the tests assumed that the value of C d did not vary during the test; to assess the effect of uncertainty in C d , the results are also shown for C d ¼ 0.6 AE 5%.
Figs. 8 and 9 show the results, which correspond with those shown in Fig. 7, whereby the correlation of the shape of the curves between the measured and modeled flows increased with time. The offset visible was speculated to be a result of the underprediction of the initial leak area from the experimental data (measurement error AE3.82 mm 2 ).
It may be inferred from the observations that the relative change of the leak flow rate over time decreased as the initial slit length was reduced. This might be explained by the knowledge that the relative change of area over time decreases as the slit length decreases, as described in Eq. (6), due to the increased structural stiffness of the pipe. The results shown support confidence in the validity and effectiveness of the developed model in capturing the pressure-dependent leakage behavior of longitudinal slits in a thick-walled viscoelastic pipe.

Discussion
The methodology and results presented here aimed to develop a generalized analytical model to describe the leakage behavior of stable longitudinal slits in a thick-walled viscoelastic pipe, through the combination of numerical simulations and experimental data. Initially, a model form was sought that could accurately capture the pressure-dependent area behavior of longitudinal slits in a linear elastic pipe prior to the calibration of the viscoelastic components.
The assessment of the slit face loading highlighted the need to fully understand the influence of this loading and how it would change with the increasing area. An appreciation of the external ground conditions and the leak hydraulics is necessary to fully explore this phenomenon. The current leakage model [Eq. (8)] presents an idealized representation of the leakage behavior of fully submerged longitudinal slits in MDPE pipes, assuming negligible external loading, and is therefore not truly representative of leaks found in real water-distribution systems. Further work is required to fully understand the behavior of buried leaking pipes. Consideration of the structural interaction between an external medium and the pipe as well as the effects on the leakage hydraulics will further enhance the ability to accurately quantify the behavior of this particular failure type. This does prevent the dimensionally homogeneous model developed in this paper providing a unique insight into the fundamental structural behavior of this failure type. Furthermore, the residual manufacturing stresses were found predominantly to affect the initial leak size and had a negligible impact on how the slits changed shape with time.
The statistical evaluation of the FEA results indicated the most dominant parameters (the slit length, pipe diameter, wall thickness, pressure, and elastic modulus) in the observed structural behavior, allowing less significant parameters (the slit width and Poisson ratio) to be removed from further analysis. Exploring the dimensionless coefficient C 1 showed that a better model fit could be achieved if this were assumed to be a function of πD=L c , the relative size of the leak to the pipe diameter.
An interesting finding was the triviality of the slit width, which was shown to directly influence the initial area only. Similarly, a maximum relative slit length limit was highlighted in the analysis, whereby the test sections with a slit length-topipe circumference ratio greater than 1 deviated from the derived model predictions.
The fundamental difference between the leak area models for longitudinal slits in linear elastic pipes, as presented in Eq. (6), and the expression from van Zyl and Cassa (2014), is that the models were derived from thick-walled and thin-walled FEA, respectively. Comparison of the predictive capabilities of the two models ( Fig. 10) highlights the consequence of this, whereby the change of the leak area was significantly overpredicted by Eq. (6) for thin-walled pipes and vice versa for the van Zyl and Cassa (2014) model. Note that there is a much higher frequency for the thickwalled simulated models in comparison with the thin-walled models.
The disparity in the results for each predictive model presented is surmised to result principally from the inherent difference in the cross-sectional material stress distribution and hence the localized deformations between thick-and thin-walled pipes. Additionally, the boundary conditions used by Cassa and van Zyl (2008) predefined the mode of deformation, limiting the displacement in the direction of the leak and potentially reducing the change of the leak area. The net deformations of the pressurized leaks were judged to be dependent on two primary modes of deformation, material bulging and bending moment, recognized by De Miranda et al. (2012). The boundary conditions used by Cassa and van Zyl (2008) may, however, reflect the boundary conditions for buried pipes where the surrounding soil acts as an additional pipe restraint, limiting the deformation due to the applied moment and thereby increasing the significance of the localized pipe bulging. The boundary conditions presented here reflect the experimental setup, thereby providing a platform to accurately calibrate and validate the derived dynamic leakage model using a combination of numerical simulations and physical observations. Eq. (4) represents a simple and computationally efficient expression to describe the dependent dynamic leak area, based on the numerical modeling. This is particularly important in incorporating the time-dependent viscoelastic behavior, which would require extreme levels of data processing if other derived models were used (De Miranda et al. 2012).
There is no clear means of determining the most effective representation of viscoelasticity to use (Purkayastha and Peleg 1984). The decision criterion is therefore based on the most efficient model to describe the observed physical behavior. The generalized Kelvin-Voigt mathematical representation was demonstrated to be an effective means of predicting the structural response of such dynamic leaks. The calibrated creep compliances here are relevant only for the exact MDPE used in these tests; it is therefore not feasible to directly equate the components for the different viscoelastic materials, including the various grades of MDPE, high-density polyethylene (HDPE), and PVC.
The aim of the calibration process was to derive a dimensionally homogeneous model that captured the full characteristic response of the observed dynamic behavior, without directly apportioning this to specific molecular processes. An 11-component viscoelastic model was shown to accurately describe the short-term and longterm structural responses of the leak and the instantaneous elastic and retarded viscous components. If the long-term response only (for example, the leak area's opening after 24 h) were required, it would be possible to simplify the model by removing the lower-order retardation time components. However, this would, for example, greatly reduce the capacity of the model to accurately quantify the total daily leakage volume.
The presented leakage model [Eq. (8)] offers potential benefits to understanding leakage control through pressure management. Currently, the simple exponent approach does not accurately reflect the true complexity of the pressure-leakage relationship observed in viscoelastic pipes. The model proposed in this article will allow engineers to more accurately understand the implications of changing the pressures in their systems, and the resulting increase or decrease in the leakage rates.

Conclusion
A generalized analytical model to describe the structural response and subsequent leakage through a longitudinal slit in a thick-walled viscoelastic polyethylene pipe was developed, calibrated, and validated through the use of numerical simulations and experimental data. The derivation of a dimensionally homogeneous expression defined the structural dynamics of a leak in a linear elastic pipe material, with experimental data used to calibrate the timedependent elastic modulus in a PE pipe. Validation of the developed relationship was successfully demonstrated by the use of independent experimental results.
The proposed model provides a means to assess the timedependent leakage response for longitudinal slits for a range of parameters, including the leak geometry, pipe dimensions, and hydraulic loading conditions. This fundamental understanding of how leaks behave in water-distribution systems is crucial to improving the current leakage levels by enabling new and improved leakage management strategies.

Acknowledgments
The project reported in this paper was supported by the Engineering and Physical Sciences Research Council (EPSRC) Grant No. EP/G015546/1.

Notation
The following symbols are used in this paper: A 0 = constant initial leak area; C d = coefficient of discharge; C i = empirical constants of the derived equations; D = pipe diameter; dA = change of leak area; E = pipe material Young's modulus; E inst = instantaneous elastic modulus of viscoelastic material; g = gravitational acceleration; h = driving head across the leak orifice; JðtÞ = viscoelastic creep compliance function; J i = viscoelastic creep compliance coefficients; L c = manufactured slit length; m = constant that determines the rate that area changes with internal head; P slit = slit face pressure loading; P = internal pipe pressure; Q = leakage flow rate; s = pipe wall thickness; T = temperature of the pipe; t = time coordinate;  Table 1. W = manufactured slit width; ν = pipe material Poisson's ratio; ρ = fluid density; and τ i = viscoelastic retardation time components.