## Abstract

Life prediction in energy infrastructure such as gas pipelines is important to maintain the integrity of such systems. This paper explores a life prediction model for polyethylene materials in natural gas distribution pipelines under creep damage. The model uses a power law equation to describe the crack growth rate and an asymptotic solution for the stress intensity factor (SIF) calculation considering local geometry variations. The SIF solution considers the effect of stress concentration introduced by common damages in pipes such as rock impingement and slit. An effective initial crack size model is proposed for the life prediction of plastic pipes considering the intrinsic initial defect. Large loading-induced plastic deformation is included by a correction factor in the crack growth model. The model is calibrated and validated using experimental data on Aldyl-A pipes with different types of damage. Due to the stochastic nature of the crack growth process, uncertainty quantification is performed, and Monte Carlo (MC) simulation is used to estimate the failure probability. The predicted probabilistic life distributions under different loading conditions are compared with the experimental data. Some conclusions and future work are drawn based on the proposed study and experimental validation.

## 1 Introduction

Many failure modes can occur during the operation of gas pipeline depending on the material system (e.g., steel pipeline for transmission and plastic pipeline for distribution pipeline systems), environment and operation conditions (e.g., fatigue and corrosion), and installation and manufacturing process (e.g., welding for steel pipes and butt fusion for plastic pipes). The gas distribution pipelines are mainly made of polymer materials and they are working under constant pressure, which makes them under substantial creep damage under room temperature [1]. The safety issue arises in the aging pipeline systems due to several recent accidents [2,3]. This study focuses on the creep damage in plastic pipelines. The creep process is also referred to as slow crack growth (SCG) [4] and other types of damages (e.g., rock impingement, slit or squeeze-off) may accelerate the failure process [5]. These damages are often introduced by excavation or by surrounding soil movements. An accurate remaining life prediction method for these pipes is of critical importance for a good maintenance planning to reduce the risk of failure. The proposed research will develop a life prediction model for the plastic pipeline considering the effect of different types of damage as local stress risers.

Some earliest studies on creep crack growth could be found in the 1960s. The crack growth process is fundamentally different in polymeric materials than that in metallic materials. There are crazes forming at the crack tip that bridges the edges of the crack and stress in crazing regions has a complex distribution as the deformation changes. To model this behavior, Schapery [6] divided the crack into a process zone and a far field region and used viscoelastic fracture mechanics approach to determine the crack growth rate as a function of the stress intensity factor (SIF). A model based on linear elastic fracture mechanics was proposed by Williams and Marshall [7], which used a time-dependent relaxation modulus to replace the Young's modulus so that the stress at the craze decreases over time. A model based on cohesive zone was developed to describe the creep crack growth process [8]. The craze was modeled as a one-dimensional zone that grows and ruptures as a function of time. One common concept in the above research is that fracture mechanics and parameters such as SIF and energy release rate are used for the creep damage prediction in polymeric materials. Some recent studies [9,10] used extrapolation from the fatigue data to predict the creep failure in polyethylene. In Ref. [10], the author states that there exists a good correlation for the two failure modes (creep and fatigue) and they share some similar characteristics. A Paris' law equation was fitted to experimental data to predict fatigue failure in Ref. [9]. In previous studies [11,12], the fatigue damage accumulation was treated as an equivalent crack growth process and the method was used for fatigue life prediction in both smooth and notched metallic specimen. The notch effect is introduced as local stress concentration factor in the model. The equivalent crack growth model with fracture mechanics analysis can unify the crack initiation and propagation into a single framework. Thus, it might be feasible to use this concept and extend it for polymer material creep life prediction given many similarities found for polymer material creep and fatigue crack growth [9,10]. The motivation of this research is to extend the equivalent crack growth model for the life prediction of polymer materials under creep damage.

In this paper, the equivalent crack growth model for creep damage is developed. The equivalent crack growth rate is modeled using a power law equation. The crack is assumed as a semicircular surface crack in the longitudinal direction and the SIF solution is derived in Ref. [13]. An asymptotic solution for the SIF is adopted from Liu and Mahadevan [14] to consider the damages in pipes as local stress risers. Next, the model was calibrated and validated using experimental creep data on Aldyl-A pipe. The experiments were done on a total of five groups of pipes with different surface conditions. One group is used as the reference group to calibrate the model parameters and others are used to validate the model. Due to the stochastic nature of crack growth, uncertainty is introduced in the model by considering the model parameters as random variables. Monte Carlo (MC) simulation is used for probabilistic life prediction. Both mean value and probability distribution distributions are compared between model predictions and experimental observations.

## 2 Life Prediction Model Development

### 2.1 Equivalent Crack Growth Model.

*C*and

*m*in the equation are material constants, which will be calibrated using the experimental data.

*K*is the SIF. Assuming the crack to be a semicircular crack in the inner surface of the pipe along the direction of the pipe (Fig. 1), the solution of the SIF can be found in Ref. [13] as

*a*is the crack length,

*σ*is the hoop stress, and

*Q*and

*F*are shape factor and boundary correlation factor and can be calculated given the geometry of the crack. The effect of damage is regarded as local stress concentration at the crack tip. An asymptotic solution for SIF that considers the effect of stress concentration factor is defined in Ref. [14] as

where *K _{t}* is the stress concentration factor introduced by damage and

*d*is the measurement of the damage. Figure 1 showed an illustrative figure of the damage.

*Q*and boundary correction factor

*F*can be found in Ref. [13]. In our case,

*Q*=

*2.464 and*

*F*is a function related to the crack length

*a*and the thickness of the pipe

*t*as

*a*to critical crack length

_{i}*a*and get the failure time as a function of stress

_{c}*σ*

where $I=\u222baiacFa+d{1\u2212exp\u2009[\u2212ad(Kt2\u22121)]}da$. The initial crack length could be observed using the scanning electron microscopy (SEM) imaging or using the equivalent initial flaw size [11,15,16]. The critical crack length is usually in the range of mm and the crack growth rate near the critical length is very fast. A small variation of the critical crack length is not important for the life prediction and is fixed as a constant value in the later probabilistic analysis.

*σ*is the yield stress of the material under static loadings. In the proposed study, this value is estimated using the experimental data when the failure time is extrapolated to zero in a stress-life plot. The equivalent crack length

_{f}*a′*can then be calculated as

The equivalent crack length in Eq. (8) will be used in the integral in Eq. (6) for the life prediction model. The model parameters *C* and *m* can be calibrated using one set of the obtained stress-life data. In the proposed model, the initial crack size is not precisely known, and an equivalent initial size is determined to match the stress-life data from experiments. Thus, it is named “equivalent crack growth model.” The idea is to use an approximated model to match the macrolevel behavior of the SCG, i.e., the life of the material. It should be noted that the calibrated model may not be used to describe the actual physical crack growth rate and are only used for the life prediction.

### 2.2 Uncertainty Quantification and Probabilistic Life Prediction.

Due to the stochastic nature of SCG, the huge scatter in the experimental observation indicates that a deterministic model such as the one discussed above may not be accurate in predicting the life of the material. Comparing to the deterministic life prediction, a probabilistic approach would be more appropriate in describing the life distribution. To introduce uncertainty into the model, some parameters are regarded as random variables and the distribution of life prediction is approximated using MC simulation. Details are discussed below.

In the proposed model, some parameters are related to material properties, such as *C*, *m*, *σ _{f}*,

*a*, and

_{i}*a*while others related to geometry of the damage, such as

_{c}*d*and

*K*. In the uncertainty analysis study, only three parameters, namely,

_{t}*C*,

*σ*, and

_{f}*K*, will be regarded as random variables. These three parameters will each represent the variability from material property and the geometry of the damage and are assumed to be independent.

_{t}*m*is treated as a constant. It is a parameter related to the failure mechanism and should be the same for the same material. The observation for the initial and critical crack length

*a*

_{i}and

*a*is hard to measure and cannot be determined using one set of calibration. All the variabilities from the material property are “lumped” into the randomness of parameter

_{c}*C*and

*σ*. In Ref. [19], a similar approach was used in predicting the scatter bond in the stress-life curves. Although this is a highly simplified approach, it is convenient to calculate the statistics of the parameter from experimental data. The validation of the model showed good accuracy for the probabilistic life prediction comparing with the experiment observations. The variability of the parameter

_{f}*K*changes the local stress and directly affects the crack growth. The parameter

_{t}*d*is in the range of 10–20 mm. It is much larger than the crack length and its variability has minimal impact on the life prediction. Hence, the variability of

*d*is ignored here since the crack growth rate is mainly dependent on the local stress concentration. It should be noted that ignoring variability of

*d*may not be applicable to other materials and structural components since the initial crack size may be comparable to the notch size and its effect cannot be ignored. Furthermore,

*C*and

*σ*represents the material variability and

_{f}*K*represents the pipe geometry variability due to the manufacturing and/or environmental conditions. They are assumed to be independent. Statistical analysis is used to select the best probability distribution to describe the variability among several candidate models, such as Gaussian, lognormal, Weibull and logistics distributions. The probabilistic life prediction is based on the uncertainty quantification of the model parameters. MC simulation [9,20] is used to calculate the failure probability.

_{t}## 3 Model Validation With Experimental Data

### 3.1 Experimental Data Description.

The experimental datum is provided by Gas Technology Institute (GTI). It is divided into five groups by their fracture mode and the damage type, namely, ductile group where failure mode is dominated by ductile deformation, SCG group, low ductile inner wall (LDIW) group, LDIW with indentation, and LDIW with squeeze-off. There are 450 data points in total. The test used rate process method for accelerating the pressure test and the results are shifted to 23 °C (73.4 °F) using bidirectional shift factor [21]. The data contain the measured hoop stress and failure time, and shifted hoop stress and failure time, as well as the stress concentration factor. The experimental data are plotted in Fig. 2 in log–log scale. The trend of the data is showing that the life in each group of pipes is log-linearly related to the hoop stress. The transition from ductile failure to SCG-dominated failure is modeled in the equivalent crack length given in Eq. (8). The SCG group is used as reference data for calibrating the material parameter *C* and *m*. The ductile data are used to calculate the plasticity correction stress.

The geometry of the pipe is shown in Fig. 3. The outer diameter is 2.375 in. The diameter to thickness ratio (standard dimension ratio value) is 11.

It is observed that the average initial crack length is about 0.001in. A SEM image showing the initial crack can be found in Fig. 4. In this study, the initial crack length is set to be *a _{i}* = 10

^{−3 }in. and we assume the critical crack length is

*a*= 0.1 in.. The average value for

_{c}*d*is assumed to be

*d*=

*1 in. and regarded as a constant.*

### 3.2 Probabilistic Model Parameters Calibration and Characterization.

where *F _{X}*(

*x*) and

_{i}*S*(

_{n}*x*) are the theoretical and empirical CDF at

_{i}*x*, respectively.

_{i}*x*corresponds to the

_{i}*i*th data. The candidate distributions used for the experimental data includes: normal, lognormal, logistic and Weibull distribution. First, the data in each ground are fitted in each of the listed distribution type using maximum likelihood method. Then the empirical CDF is compared with the CDF from the fitted distribution. The test was done using the matlab built in function in the statistics and machine learning toolbox. Table 1 recorded the test result. The value represents the test decision for the proposed distribution. The higher the value indicates better the goodness-of-fit. As we can see that for SCG, LDIW and squeeze-off group, a logistic distribution can best model the

*K*distribution. Although for LDIW indent group, the best distribution is a lognormal, the score for logistic distribution is considerably well. Overall, a logistic distribution tends to match the data best.

_{t}Distribution type | SCG | LDIW | LDIW indent | Squeeze-off |
---|---|---|---|---|

Normal | 0.116 | 0.107 | 0.492 | 0.713 |

Lognormal | 0.489 | 0.0532 | 0.624 | 0.935 |

Weibull | 0.0140 | 0.0255 | 0.556 | 0.348 |

Exponential | 4.79 × 10^{−29} | 1.00 × 10^{−31} | 1.19 × 10^{−08} | 4.63 × 10^{−10} |

Logistic | 0.690 | 0.219 | 0.488 | 0.990 |

Distribution type | SCG | LDIW | LDIW indent | Squeeze-off |
---|---|---|---|---|

Normal | 0.116 | 0.107 | 0.492 | 0.713 |

Lognormal | 0.489 | 0.0532 | 0.624 | 0.935 |

Weibull | 0.0140 | 0.0255 | 0.556 | 0.348 |

Exponential | 4.79 × 10^{−29} | 1.00 × 10^{−31} | 1.19 × 10^{−08} | 4.63 × 10^{−10} |

Logistic | 0.690 | 0.219 | 0.488 | 0.990 |

Figure 5(a) showed the empirical CDF for *K _{t}* values for each group plotted against the CDF of the fitted logistic distribution. It can be concluded that a logistic distribution can be used to describe the randomness on parameter

*K*. Similar approach was done to parameter

_{t}*C*.

Equation (10) is the mean curve for the stress-life relation. Comparing with Eqs. (6) and (10) indicates that *m *=* *2.5632. The mean of the stress concentration factor from the fitted distribution is used as the value for $Kt$ to evaluate the integral $I$ in Eq. (6), with the parameter values characterized in Sec. 3.1. The mean value for *C* is calculated from the intersection with the vertical axis of the regressed function: log*C* = −13.3442. To generate the distribution for *C*, the data points in SCG group are substituted into the log-linear function in Eq. (10) with a fixed *m*. This gives an array of values for *C*. For the collected data, it is found that log*C* follows a normal distribution with mean −13.3442 and standard deviation 0.2960. A plot of empirical CDF for log*C* along the fitted normal distribution can be found in Fig. 5(b), the result is showing good fitness.

*σ*needs to be calibrated using the ductile data. Again, linear regression with mean square error method can give the log-linear relation for the ductile region and it is written as

_{f}Equation (11) is the mean curve for the stress-life relation in the ductile region. Similar to the process discussed with Eq. (10), the mean value for log(*σ _{f}*) can be calculated from the vertical intercept using the equation in Eq. (6). Due to the randomness in the ductile data, σ

_{f}is considered as a random variable. And the value of log(

*σ*) is fitted to have a Normal distribution with mean 3.33 and standard deviation 0.0264.

_{f}By MC simulation, the 90% confidence bounds for the stress-life curve can be plotted. The MC simulation is done by randomly sample from the fitted distribution for the random variables *C*, *σ _{f}*, and

*K*. For each set of random variables, the life is calculated at each stress level. The upper and lower 5% of the calculated life is chosen as the boundary. Figure 6 showed the curve for the average failure time against stress level in double log scale and its 90% confidence bounds along with the experimental data for SCG and ductile group. Since this is the calibration for the reference data, it is not surprising that the model prediction has agreement with the data.

_{t}Applying the equivalent crack length model in Eq. (8), the approximation for the transition from ductile to SCG can be modeled. Figure 7 showed the asymptotic stress-life curve. After tuning the parameter, it is found that when choosing log(*σ _{f}*) = 3.22 we have a better fit for the prediction curve.

To this point, the calibration for the model parameter is finished and the uncertainty quantification for the model parameter is concluded in Table 2.

Distribution parameter | |||
---|---|---|---|

Parameter name | Distribution | Mu | Sigma |

$ai$ | Constant | 0.001(in) | |

$ac$ | Constant | 0.1 (in) | |

$d$ | Constant | 1 (in) | |

$log\u2009C$ | Normal | −13.3442 | 0.2960 |

$m$ | Constant | 2.5632 | |

$Kt$(SCG) | Logistic | 2.6192 | 0.2445 |

$Kt$(LDIW indent) | Logistic | 8.2307 | 0.9203 |

$Kt$(LDIW) | Logistic | 5.5326 | 0.4906 |

$Kt$(squeeze-off) | Logistic | 6.4095 | 0.4745 |

$Kt$(ductile) | Logistic | 2.6578 | 0.0766 |

log(σ)_{f} | Normal | 3.2 (psi) | 0.0264 |

Distribution parameter | |||
---|---|---|---|

Parameter name | Distribution | Mu | Sigma |

$ai$ | Constant | 0.001(in) | |

$ac$ | Constant | 0.1 (in) | |

$d$ | Constant | 1 (in) | |

$log\u2009C$ | Normal | −13.3442 | 0.2960 |

$m$ | Constant | 2.5632 | |

$Kt$(SCG) | Logistic | 2.6192 | 0.2445 |

$Kt$(LDIW indent) | Logistic | 8.2307 | 0.9203 |

$Kt$(LDIW) | Logistic | 5.5326 | 0.4906 |

$Kt$(squeeze-off) | Logistic | 6.4095 | 0.4745 |

$Kt$(ductile) | Logistic | 2.6578 | 0.0766 |

log(σ)_{f} | Normal | 3.2 (psi) | 0.0264 |

### 3.3 Sensitivity Analysis for the Random Variables.

Sensitivity analysis is performed to study the effect of each random variable on the predicted life. In this study, three random variables are being considered, namely, *C*, *K _{t,}* and

*σ*. To evaluate the effect of the randomness in each parameter, the parameters are perturbed by ±10% from the mean value of their fitted distribution. In each case, only one parameter is perturbed while the other two parameters are fixed at their mean value. Three stress-life curves can be generated with the sensitivity analysis. Figure 8 shows the mean and perturbed stress-life prediction to illustrate the impact of each random variable on the predicted life-stress curve.

_{f}From the plots in Fig. 8, it can be concluded that the three variables *C*, *K _{t}*, and

*σ*will affect the stress-life prediction uncertainty. The effect of

_{f}*C*,

*K*is closely related to the behavior in the creep damage region while the effect of

_{t}*σ*is limited to the ductile region. It also appears that the parameter

_{f}*C*has a relatively larger impact on the life prediction in the long-life regime (e.g., creep damage region) compared to the parameter

*K*.

_{t}### 3.4 Comparison of Model Predictions and Experimental Observations.

By changing parameter *K _{t}* to the values in other groups of pipes with the calibrated model parameters, the stress-life for other types of damage can be predicted. The prediction results are plotted with the experimental data for LDIW, indentation, and squeeze-off group in Fig. 9. Although some discrepancy can be observed, the proposed model prediction showed good agreement with the experimental data since most data points fall in the confidence bounds. The results indicate that the model is capable of predicting the life of the pipe with different types of damage.

With the equivalent crack length concept, the model can approximate the transition from ductile to SCG process. By applying Eq. (8) to the linear model in Eq. (6), the asymptotic prediction curve can be plotted (Fig. 10). The ductile data are plotted to show the transition from ductile to creep. The curve captured the behavior of these two failure regions. The two regions are governed by the same α-relaxation process, only the degree of constraint at the damage tip is different. In the ductile region, the entire pipe wall is subject to the applied stress while in the SCG region, only a small local volume. A stress concentration introduced in a polyethylene pipe would create a complex situation at the crack tip. SCG is still ductile creep but only occurs in a small volume at the crack tip. The visual comparison with the experimental observation showed that the model has the ability to predict the life of the pipe material with different types of damage. Discrepancy can be observed for the ductile data as a large portion of data points falls out of the prediction bound at 10^{3.2} to 10^{3.4 }psi stress level.

The above validation focuses on the general visual comparison of the data and prediction. A more rigorous model validation is performed by the comparison of the empirical CDF between the experimental data scatter and the model predicted scatter. The error analysis was done by taking the difference in failure time of the experimental data against the predicted mean. For each data point, the error is calculated as the actual life from test minus the predicted mean failure time at the corresponding stress. The errors were collected and plotted as an empirical CDF in Fig. 11 compared with the MC sample errors. Figure 11(a) represents the calibration quality and it can be seen that the model predicted scatter is very close to the experimental observed scatter. Good agreement can be observed in Fig. 11(b), which represents the predicted scatter with experimental measured scatter for LDIW group. For indentation and squeeze-off condition (Figs. 11(c) and 11(d)), the experimental observed scatter is slightly less than the model predicted scatter and the experimental median value is also less than the experimental value. Future study is needed to explain this discrepancy. It is possible that the indentation and squeeze-off alter the material properties, which is weaker than the virgin materials. Thus, the median life will be shorter than the model predicted using the virgin material properties. The difference is not huge, but more detailed experimental observation and theoretical study is required.

## 4 Conclusion

The proposed study used an equivalent crack growth model to predict the creep life for polymer pipe materials with different types of damage. Overall, the prediction showed good agreement with experimental data. Conclusions can be drawn from the presented work:

The power law SCG equation can predict the creep life of plastic pipes. This is similar to the classical approach in fatigue crack growth prediction model.

The applied equivalent crack length model can approximate the transition from SCG to ductile behavior.

The uncertainty quantification of the model parameters has captured the variability in the experimental observations. The stress concentration factor

*K*follows a logistic distribution, the logarithm of material parameter_{t}*C*and frictional stress*σ*, log_{f}*C*and log(*σ*), each can be modeled as a normal distribution._{f}The proposed model can generally achieve good prediction with the calibrated model parameters. But the predicted life seems to be longer than the average of experimental data in LDIW indent and squeeze-off group. This may be due to that indentation and squeeze-off damage effect the material property by the large inelastic strain history.

Future studies are needed based on this research. More experimental data would greatly help the further validation of the proposed model. A direct measurement of the creep crack growth rate of the material may be useful to justify the underlying physics of the equivalent crack growth. The model should also include the loading history for prediction in practical applications. The effects of large inelastic damage (indentation and squeeze-off) need further experimental studies (e.g., microstructural imaging) to estimate the change of material properties. The discrepancy of the model prediction and the experimental data for ductile group at higher stress level indicates additional investigation is needed to better explain the observations.

## Acknowledgment

The authors thank the technical discussions and financial sponsorship from DOT-PHMSA.

## Funding Data

DOT-PHMSA (Contract Number: DTPH5615T00007; Funder ID: 10.13039/100006284).

## Nomenclature

*a*=crack length

*a*=_{i}, a_{c}initial, critical crack length

*C, m*=material constant

*d*=notch size

*F*=boundary correction factor

*K*=SIF

*K*=_{t}stress concentration factor

*Q*=shape factor

*t*=thickness

*t*=time

*T*=_{f}failure time

*σ*=hoop stress