Multiphysics Modeling of Power Transmission Line Failures Across Four US States

Abstract

The failure of overhead transmission lines in the United States can lead to significant economic losses and widespread blackouts, affecting the lives of millions. This study focuses on analyzing the failure of transmission lines, specifically considering the effects of wind, ambient temperature, and current demands, incorporating minimal and significant pre-existing damage. We propose a multiphysics framework to analyze the transmission line failures across sensitive and affected states of the United States, integrating historical data on wind and ambient temperature. By combining numerical simulation with historical data analysis, our research assesses the impact of varying environmental conditions on the reliability of transmission lines. Our methodology begins with a deterministic approach to model temperature and damage evolution, using phase-field modeling for fatigue and damage coupled with electrical and thermal models. Later, we adopt the probability collocation method to investigate the stochastic behavior of the system, enhancing our understanding of uncertainties in model parameters, conducting sensitivity analysis to identify the most significant model parameters, and estimating the probability of failures over time. This approach allows for a comprehensive analysis of factors affecting transmission line reliability, contributing valuable insights into improving power line’s resilience against environmental conditions.

1. Introduction

The power grid’s components are so interconnected that the failure of a single component can cause a widespread outage. Overhead transmission lines exposed to dynamic weather conditions and current loads are particularly vulnerable. Such vulnerabilities not only lead to frequent blackouts across the United States, affecting millions annually, but also disrupt daily life and impose significant economic burdens. A report by the Department of Energy estimated that weather-related power outages cost the US economy between USD 18 billion and USD 33 billion from 2003 to 2012 [1]. Some states such as California, Texas, and Michigan are particularly more vulnerable, with Florida having the highest total number of 25,348,824 affected by power outages [2]. Hence, understanding the reliability of transmission cables in such sensitive and affected regions is crucial.

Several studies have investigated the impact of environmental conditions on overhead power lines. The effect on the life of aluminum conductor steel-reinforced (ACSR) cables due to thermal stress was investigated by integrating a strength reliability analysis based on Monte Carlo simulation and fuzzy logic-based ductility analysis [3]. Additionally, the damage caused by wind loading has been extensively studied, affecting the transmission lines and the supporting towers [4]. Studies by [5,6,7] have extensively examined the dynamic responses and structural effects of wind on transmission line infrastructure. In [8], the overheating effect of wildfire due to the radiative heat transfer was studied. Similarly, the ice and wind combined effect was studied in [9]. Several studies have been made to understand the effect of weather conditions on the thermal rating of the conductor [10,11,12]. In addition to environmental factors, the choice of conductor material plays a crucial role in determining the reliability of transmission wires.

All these studies laid a solid basis for understanding the failure mechanisms of power transmission. While much of this research has focused on the thermal and electrical modeling of transmission lines, none considered the mechanical model that relates to fatigue and damage that influence the overall conductivity of the material. Transmission lines are subjected to dynamic wind, temperature, and electrical loads, which can induce fatigue and damage to the wires. This damage compromises the conductivity of the conductor, leading to increased heat generation. Therefore, a coupled model that integrates mechanical properties with electrical and thermal aspects can significantly improve failure prediction.

Recently, phase-field models have been integrated with electrical and thermal models to study various phenomena [13,14,15]. Over the years, these models have been crucial in simulating brittle [16,17,18,19,20,21], ductile [22,23,24,25,26], isothermal fatigue fracture [27], and non-isothermal fatigue fracture mechanisms [28], They accurately capture phenomena such as crack initiation, propagation, branching, and coalescence—events commonly observed in dynamic fractures. Based on these applications, incorporating the phase-field model with electrical and thermal models becomes crucial in studying transmission line reliability. Further incorporation of stochastic analysis can enhance the precision of lifespan predictions. This approach enables a more comprehensive understanding of the factors that affect transmission line durability in real-world scenarios.

Although the phase-field methods are extensively used for damage evolution, discrete dislocation dynamics (DDD) can provide valuable insights into the underlying microstructural mechanisms [29]. The interaction of dislocation can significantly affect the behavior resulting in failure [30,31]. The normal behavior of materials also changes under high operating temperatures, resulting in visco-elastic behavior. Integrating fractional visco-elasto-plastic models as in the studies [32,33,34,35] can significantly improve the reliability analysis. Even considering the uncertainty in the fractional order further improves capturing the behavior of the system [36].

The classical Monte Carlo (MC) method [37,38] is a standard benchmark for studying stochastic solutions that compute the quantity of interest (QoI) in a straightforward fashion. However, this method is characterized by a slow convergence rate and requires a large number of realizations. Some well-established methods, such as polynomial chaos [39,40] or its generalization via Galerkin projection [41,42,43], require modifications to the governing equations for stochastic analysis. This requirement makes these methods intrusive and potentially impractical for complex problems since such modifications may not be feasible or could overly complicate the analysis. To address this, non-intrusive techniques must be used; one such technique is the probability collocation method (PCM) [44,45]. PCM preserves the simplicity of solution structures and allows for independently sampled realizations, thereby achieving better convergence compared to the traditional Monte Carlo (MC) method. Although the curse of dimensionality challenges PCM due to tensorial products, this issue can be effectively addressed using techniques such as sparse grids [46] or active subspace methods [47,48,49].

In a recent study [50], the authors integrated the phase-field model with thermal and electrical modules to predict the lifespan of overhead transmission lines. The study initially examined a deterministic solution and later implemented the PCM for stochastic analysis. However, a more comprehensive analysis involving historical data across various specific states would provide a more realistic assessment. The reliability and durability of transmission lines vary under different weather conditions, highlighting the need to adjust the weather condition assessments to particular cases. Furthermore, analyzing the lifespan of transmission lines with insignificant initial damage could offer valuable insights into the extent to which pre-existing damage influences overall reliability.

In this study, we develop a coupled fatigue and damage-based phase-field model with thermal and electrical aspects to investigate the reliability of overhead transmission lines under environmental conditions and the presence of initial damage, including the uncertainties associated with these. We first investigate the reliability of Transmission Lines deterministically using the finite element method. To explore the uncertainties associated with environmental conditions and initial damage, we incorporate PCM for uncertainty quantification (UQ), sensitivity analysis (SA), and probability of failure using the limit state function. To understand the reliability of transmission lines more realistically, we incorporate several realistic scenarios for a detailed historical analysis of wind and temperature data rather than using parameterized data by focusing on the four most affected U.S. states: Texas (TX), California (CA), Michigan (MI), and Florida (FL) according to [2]. Additionally, we assess the system’s reliability under minimal damage to establish a baseline and examine how precursor damage affects the reliability of transmission lines. By analyzing these scenarios using our developed model, this research offers critical insights into the conditions-including pre-existing defects-that most significantly influence the durability and reliability of transmission lines. These insights can guide the development of more resilient power line infrastructures designed to withstand diverse environmental conditions and extend their lifespan. Ultimately, this work contributes to improving the reliability of the transmission line.

This paper is structured as follows: Section 2 presents the problem statement, including four representative scenarios. Section 3 presents a multiphysics model and discusses each model. Additionally, the section discusses the historical data analysis for wind and temperature of each specific scenario, including the reliability analysis procedure using a limit state function. Section 4 presents the deterministic approach and solution, including one-dimensional finite element discretization. Section 5 presents the stochastic analysis methods and solutions using PCM as the building block. Section 6 presents the verification and analysis of the model. Section 7 addresses the conclusion.

2. Problem Statement

The operating temperature is crucial for the reliable operation of overhead transmission lines. High operating temperature affects the wire’s life and reliability. Even the presence of damage in the wire further elevates the temperature These issues often remain undetected until a cable rupture occurs. Additionally, environmental conditions, including temperature, speed of wind, and current load, can adversely affect these parameters, pushing them beyond acceptable levels.

Therefore, operating temperature is considered the primary factor for our analysis. Along with physical effects, material parameters, and loading conditions, these elements determine the failure state of the transmission line over its lifespan. This study aims to understand the multiphysics effects on these critical factors, which are crucial for assessing the transmission line’s reliability.

Representative Scenarios

Previous research on the reliability of transmission lines often relied on parameterized values for wind and temperature, which lack precision. This study aims to improve this by incorporating real data from four representative scenarios: TX, CA, MI, and FL. We source wind and temperature data from [51,52]. The data for each state are presented in

Table 1. Wind speed (ft/s) [51] and temperature (K) [52] data by state.

TX	ft/s	11.59	12.32	12.76	13.05	12.17	10.41	9.09	8.65	9.24	10.12	10.85	11.15
	K	287.54	286.98	294.82	294.71	298.71	301.21	303.48	303.48	300.37	295.21	292.04	286.15
CA	ft/s	2.05	2.79	3.08	3.52	3.23	3.08	3.08	2.93	2.49	2.20	2.05	2.05
	K	286.65	288.04	287.76	290.09	293.37	294.21	295.71	297.54	297.15	295.15	289.76	287.65
MI	ft/s	16.13	15.40	15.25	15.25	13.35	12.32	11.59	10.71	11.29	13.49	15.11	15.40
	K	273.37	272.09	278.26	280.98	287.43	294.87	298.43	296.21	290.59	283.71	280.76	273.65
FL	ft/s	12.61	12.91	14.23	14.23	13.20	11.00	10.56	10.41	10.85	12.61	12.91	12.17
	K	295.09	295.98	298.26	300.87	299.48	302.26	303.09	302.76	301.93	301.04	298.43	293.93

.

3. Methodology

This study aims to analyze the failure mechanisms in overhead transmission lines by developing a coupled model that integrates thermal, mechanical, and electrical aspects. We systematically assess how mechanical properties, in terms of damage, affect the electrical and thermal aspects of the cable in the presence of environmental conditions contributing to line reliability. Figure 1 illustrates the transmission line mounted on the transmission tower, highlighting the cable failure.

We consider a cable to be a one-dimensional domain. The cable is supported by two towers to maintain the sag on it. We consider the projected span equal to its effective operational length for simplification, as shown in Figure 2. While sag can be an essential consideration, we only account for its influence on horizontal tension due to temperature variations. We do not explicitly model the sag; instead, we focus on understanding how the horizontal tension contributes to material damage and fatigue over time.

The operating temperature of the conductor affects the tension acting on it. High temperature reduces the tension, while low temperature increases it due to the effect of sag. The temperature of the cable is mainly driven by current passing through it and convective cooling from the wind acting on it. Additionally, any damage to the cable reduces its conductivity, resulting in further heat generation.

The electrical model here does not consider the electrical design of transmission lines. Rather, it focuses on the effect of current passing through the conductor independent of the voltage levels delivered to consumers.

In this study, we propose that cyclic loading accelerates the damage to the material and reduces the conductivity. This further increases the temperature due to resistivity effects and leads to greater heat generation due to Joule heating, ultimately leading to failure. The heat generation elevates the temperature of the material and accelerates aging.

This section introduces the mechanical, thermal, and electrical models that collectively simulate these processes. We detail the method for deriving cyclic loading using DFT and assessing transmission line reliability under these conditions.

3.1. Mechanical Model

For the mechanical model, we followed the principles in [28], which provides the governing equations for displacement $u_d$, damage ${\phi }_d$, and fatigue ${\mathcal{F}}_f$. ${\phi }_d=0$ represents virgin material, and ${\phi }_d=1$ indicates complete damage, while values $0<{\phi }_d<1$ represent intermediate levels of damage. ${\mathcal{F}}_f$ is considered an internal variable. Specific material evolution behaviors can be obtained from the governing equations by choosing appropriate free-energy potentials. Alternatively, one can approach this by considering the free-energy function, as follows:

$$\Psi (\nabla u_d,{\phi }_d,\nabla {\phi }_d,{\mathcal{F}}_f)=d\left({\phi }_d\right)Y{\left(\nabla u_d\right)}^2+g_c{\displaystyle \frac{\gamma }{2}}{\left(\nabla {\phi }_d\right)}^2+\mathcal{K}({\phi }_d,{\mathcal{F}}_f),$$

(1)

where Y (GPa) denotes the Young’s modulus, $g_c$ (kN/m) denotes the fracture energy release rate, and $\gamma >0$ denotes the phase-field layer width parameter (m). The degradation function is defined as $d\left({\phi }_d\right)={(1-{\phi }_d)}^2$, affecting both elastic response and the electrical conductivity. Additionally, $\mathcal{K}({\phi }_d,{\mathcal{F}}_f)$ represents the coupling of damage and fatigue, describing how damage evolves due to fatigue over time.

The original model described in [28] is time-dependent. However, we considered the quasi-static form for our problem, which allowed us to simplify the governing equations for $u_d$ and ${\phi }_d$. Still, the evolution of ${\mathcal{F}}_f$, which represents long-term aging, is modeled as a time-dependent ordinary differential equation (ODE), as follows:

$$\nabla ·\left({(1-{\phi }_d)}^{2}Y\nabla u_d\right)-\gamma g_c\nabla ·(\nabla {\phi }_d\otimes \nabla {\phi }_d)+f=0,$$

(2)

$$\gamma g_c\Delta {\phi }_d+(1-{\phi }_d){(\nabla u_d)}^{T}Y(\nabla u_d)-{\displaystyle \frac{1}{\gamma }}[g_c{\mathcal{H}}^{\prime }\left({\phi }_d\right)+{\mathcal{F}}_f{\mathcal{H}}_f^{\prime }\left({\phi }_d\right)]=0,$$

(3)

$$\dot{{\mathcal{F}}_f}=-{\displaystyle \frac{{\widehat{F}}_f}{\gamma }}{\mathcal{H}}_f\left({\phi }_d\right),$$

(4)

The potentials $\mathcal{H}\left({\phi }_d\right)$ and ${\mathcal{H}}_f\left({\phi }_d\right)$ describe the damage evolution from 0 to 1 as fatigue changes from zero to $g_c$. ${\mathcal{H}}^{\prime }\left({\phi }_d\right)$ and ${\mathcal{H}}_f^{\prime }\left({\phi }_d\right)$ are the derivatives of these potentials with respect to ${\phi }_d$. The appropriate choices of these potentials are as follows:

$$\mathcal{H}\left({\phi }_d\right)=\left\{\begin{array}{cc}0.5{\phi }_d^2\hfill & \mathrm{for}0\le {\phi }_d\le 1,\hfill \\ 0.5+\delta ({\phi }_d-1)\hfill & \mathrm{for}{\phi }_d>1,\hfill \\ -\delta {\phi }_d\hfill & \mathrm{for}{\phi }_d<0.\hfill \end{array}\right.$$

(5)

$${\mathcal{H}}_f\left({\phi }_d\right)=\left\{\begin{array}{cc}-\phi \hfill & \mathrm{for}0\le {\phi }_d\le 1,\hfill \\ -1\hfill & \mathrm{for}{\phi }_d>1,\hfill \\ 0\hfill & \mathrm{for}{\phi }_d<0.\hfill \end{array}\right.$$

(6)

We describe the evolution of ${\mathcal{F}}_f$ through ${\widehat{F}}_f$, representing the formation and growth of micro-cracks under cyclic loading conditions and the influence of temperature. ${\widehat{F}}_f$ has a linear relationship with the level of stress associated with the following virgin material:

$${\widehat{F}}_f={\rho }_{m}a\left({\displaystyle \frac{{\theta }_c}{{\theta }_0}}\right)(1-{\phi }_d)\left|Y\nabla u_d\right|,$$

(7)

where the parameter a represents the aging rate, ${\theta }_c$ (K) represents the conductor temperature, ${\theta }_0$ (K) denotes the reference temperature, and ${\rho }_m\left({\mathrm{kg}/\mathrm{m}}^3\right)$ denotes the density of the material.

The mechanical model allows damage healing when tensile stress decreases. To simulate an irreversible damage process and avoid healing mechanisms, we adopt an approach similar to that in [19]. We define ${\mathbb{H}}_e$ as the local maximum strain energy history:

$${\mathbb{H}}_e(x,t)=max({(\nabla u_d(x,t))}^{T}Y(\nabla u_d(x,t)),\mathcal{H}(x,t)).$$

(8)

We integrate the local maximum strain energy history variable ${\mathbb{H}}_e$ into the damage equation and obtain a modified equation for damage as follows:

$$\gamma g_c\Delta {\phi }_d+(1-{\phi }_d){\mathbb{H}}_e-{\displaystyle \frac{1}{\gamma }}[g_c{\mathcal{H}}^{\prime }\left({\phi }_d\right)+{\mathcal{F}}_f{\mathcal{H}}_f^{\prime }\left({\phi }_d\right)]=0.$$

(9)

3.2. Thermal Model

The original model [28] relates fatigue with the increase in temperature due to repetitive and fast loading; however, here, we focus on the long-term damage rather than the short-term increase in temperature. Although the short-term effects are crucial, the long-term reliability of the transmission lines is mainly affected by static loads, environmental conditions, and gradual material degradation. Therefore, we assume a quasi-static regime to simplify our analysis by focusing on steady-state conditions that significantly impact long-term performance. The steady-state heat equation is given by the following:

$$\nabla ·({\kappa }_m\nabla {\theta }_c)+q=0.$$

(10)

Equation (10) describes the heat balance within a system, where ${\kappa }_m$ (W/(m K)) represents the thermal conductivity, and $q\left({\mathrm{J}/\mathrm{m}}^2\right)$ denotes the net heat exchange. Joule heating, $q_j\left({\mathrm{J}/\mathrm{m}}^2\right)$, resulting from the current passing through the conductor, serves as a heat source, while convective cooling, $q_c\left({\mathrm{J}/\mathrm{m}}^2\right)$, driven by the wind, acts as a heat sink. Thus, the overall heat exchange can be represented as follows:

$$q=q_j-q_c,$$

(11)

We consider a convective heat transfer using the following relation:

$$q_c=h({\theta }_c-{\theta }_a),$$

(12)

where $h\left({\mathrm{W}/(\mathrm{m}}^2\mathrm{K})\right)$ is the convective heat transfer coefficient. We consider forced convection due to cross-flow over the cylinder using the relation given by [53], where h can be calculated using the following relation:

$$N{u}_D={\displaystyle \frac{hD}{{\kappa }_{air}}}=C_{e}R{e}_D^{m_e}P{r}^{{\displaystyle \frac{1}{3}}}$$

(13)

where $N{u}_D$ represents the Nusselt number and $Pr$ represents the Prandtl number. We first determine the Reynolds number $R{e}_D$ using a relation given by the following:

$$R{e}_D={\displaystyle \frac{v_{air}D}{{\nu }_{air}}},$$

(14)

where ${\nu }_{air}({\mathrm{m}}^2/\mathrm{s})$ denotes the kinematic viscosity of air, $v_{air}\left(\mathrm{m}/\mathrm{s}\right)$ denotes the velocity of air, and D (m) denotes the diameter of the conductor.

$C_e$ and $m_e$ are experimentally determined, and depend upon the value of $R{e}_D$.

Table 2. Experimental values of $C_e$ and $m_e$ for different $R{e}_D$ ranges.

ReD Range	Ce	me
0.4–4	0.989	0.330
4–40	0.911	0.385
40–4000	0.683	0.466
4000–40,000	0.193	0.618
40,000–400,000	0.027	0.805

shows the values of $C_e$ and $m_e$ for different ranges of $R{e}_D$.

3.3. Electrical Model

In the transmission line design, various factors such as frequency, inductance, reactance, and electromagnetic interactions with the environment and nearby conductors are considered. However, we are not considering the overall design of the transmission lines; rather, we focus on understanding how current-driven Joule heating affects the reliability of cables with existing damage. Additionally, we do not consider the transient effects of AC currents since they average out over time. Therefore, we use a DC-equivalent mean current, simplifying our approach to measure the heat source term efficiently. Further, to reduce complexity, we parameterize the base value of the current $I_{base}$, corresponding to the allowable ampacity for all-aluminum conductors with a diameter of approximately 40 mm. While the current data can be obtained through sophisticated methods, we do not focus on a specific transmission line. Instead, our emphasis is on the voltage drop due to temperature and damage-induced resistivity.

From the above consideration, we focus on solving the conservation of current through the following set of equations:

$$\nabla ·J_d=0,$$

(15)

$$J_d={\sigma }_{E}E,$$

(16)

$$E=-\nabla V,$$

(17)

where $J_d\left({\mathrm{A}/\mathrm{m}}^2\right)$ denotes current density, E (V) denotes the electric field generated by the voltage V, and ${\sigma }_E$ represents the electric conductivity at the operating temperature, which depends upon the degradation function, $d\left({\phi }_d\right)$, and non-degraded conductivity of damage, ${\sigma }_{E,T}$:

$${\sigma }_E={(1-{\phi }_d)}^2{\sigma }_{E,T}.$$

(18)

${\sigma }_{E,T}$ at the operating temperature, is related to the conductivity ${\sigma }_{E,0}$ at a reference temperature and can be obtained by the following:

$${\sigma }_{E,T}={\displaystyle \frac{{\sigma }_{E,0}}{1+{\alpha }_c({\theta }_c-{\theta }_0)}},$$

(19)

where ${\alpha }_c$ is the coefficient of resistivity of the conductor. Combining the above equations, the equation for the voltage field is as follows:

$$\nabla ·(-{\sigma }_E\nabla V)=0,$$

(20)

$q_j$ can be defined as follows:

$$q_j=J_d·E.$$

(21)

Overall, this model links damage to voltage drop, increasing the power loss due to increased resistance, thus aggravating the thermal load due to Joule heating.

3.4. Sag Consideration

The temperature of the cable influences the horizontal tension acting on it. This section focuses on obtaining the mechanical load acting on a cable fixed at two ends, forming a catenary curve. While several studies have been conducted to understand the mechanical behavior of such cables [6,54], we simplify our approach using a one-dimensional damage phase-field model. We assume the cable’s length L (m) is approximately equal to the span length $S_l$ (m). This assumption simplifies the model but does not ignore sag S; rather, it highlights the effect of sag on horizontal tension H that results in fatigue and damage.

The sag in the cable is directly related to the temperature variation affecting the horizontal tension. Higher temperatures increase the length of the cable due to thermal elongation, which increases the sag and reduces tension. Conversely, lower temperatures contract the cable, reducing the length and ultimately reducing the sag and increasing the tension. Therefore, the horizontal tension at either end of the cable results from an initial pre-tension combined with tension adjustments due to temperature variations. Based on these parameters, this model determines the appropriate mechanical loading conditions in terms of horizontal tension.

We follow the methodology presented in [55] to calculate the horizontal tension in the cable. We consider $W_b$ (kg/m) as the cable weight per unit length and $H_0$ as the initial tension due to its load, which is considered $20\%$ of the aluminum’s ultimate strength. The initial slag $S_0$ is calculated as follows:

$$S_0={\displaystyle \frac{W_b{S}_l^2}{8H_0}},$$

(22)

Although this is a simplification, the theoretical length $L_0$ needs to be considered, which is given by the following relation:

$$L_0=S_l+{\displaystyle \frac{8S_0^2}{3S_l}}.$$

(23)

The length of the cable changes with the change in the temperature according to the following classical formula:

$$L=L_0(1+{\alpha }_L\Delta {\theta }_c),$$

(24)

where ${\alpha }_L$ represents the thermal expansion coefficient. The sag due to the change in length is as follows:

$$S=\sqrt{{\displaystyle \frac{3S_l(L-S_l)}{8}}},$$

(25)

which then leads to a new calculation for the horizontal tension as follows:

$$H={\displaystyle \frac{W{S}_l^2}{8S}},$$

(26)

In this model, W denotes the total weight that could account for the weight of additional factors such as ice and wind. For this study, we consider only the wind component and calculate W as follows:

$$W=\sqrt{W_b^2+W_w^2}.$$

(27)

We follow a similar approach as in [56,57] to calculate the wind component. The wind-induced component $W_w$ (N/m) is obtained from the wind pressure $P_w\left({\mathrm{N}/\mathrm{m}}^2\right)$ and related to the wind velocity $v_{air}$ as follows:

$$P_w={\displaystyle \frac{1}{2}}{\rho }_{air}{v}_{air}^2,$$

(28)

$$W_w=P_w{C}_{D}D{sin}^2\left({\theta }_w\right){\alpha }_s,$$

(29)

where ${\rho }_{air}$ denotes the density of air, D denotes the diameter exposed to wind, ${\theta }_w$ denotes the angle between the transmission line and wind flow, ${\alpha }_s$ denotes the span factor, and $C_D$ denotes the drag coefficient; the values are shown in Figure 3 for different $R{e}_D$.

3.5. Discrete Fourier Transform

The available data for wind and temperature are discrete. We implement the discrete Fourier transform (DFT) and Fourier series to obtain the loading conditions. The following expression defines the DFT:

$$X_k=\sum _{z=0}^{Z-1}{a}_z{e}^{-i2\pi {\displaystyle \frac{k}{Z}}z}$$

(30)

where $a_z$ denotes the zth sample of a, and Z denotes the total number of the sample. We obtain the mean, $A_0$, and frequency coefficients, $A_z$ and $B_z$, for each sample from DFT. The mean and frequency coefficients are used in the Fourier series to obtain the cyclic loading equation using the following:

$$f\left(t\right)=A_0+\sum _{z=1}^{Z/2}[A_{z}cos\left(2\pi z{\displaystyle \frac{t}{T}}\right)+B_{z}sin\left(2\pi z{\displaystyle \frac{t}{T}}\right)]$$

(31)

where t is time and T represents the time period. The discrete (original data) and continuous data (obtained using the Fourier series analysis) plots of each state are shown in Figure 4, Figure 5, Figure 6 and Figure 7:

3.6. Reliability of Transmission Lines

In this study, we define a limit state function $g(R,S;t)$ for failure analysis. R denotes the threshold temperature ${\theta }_{\lim }$ beyond which failure occurs, S represents the maximum temperature ${\theta }_{\max }$ on the transmission line, and t indicates time. We specifically focus on the ${\theta }_{\max }$ that the line can withstand before the failure. This makes our function as follows:

$$g({\theta }_{\lim },{\theta }_{\max };t)={\theta }_{\lim }-{\theta }_{\max }\left(t\right)$$

(32)

While both R and S are treated as random variables in general cases, R is held constant in this study. Only ${\theta }_{\max }$ is considered a random variable. Therefore, the probability of failure $P_f\left(t\right)$ is as follows:

$$P_f\left(t\right)=P\{g({\theta }_{\lim },{\theta }_{\max };t)<0\}.$$

(33)

3.7. Multiphysics Framework

In Figure 8, we present a schematic representation outlining our approach to modeling transmission line failures. The model consists of the main governing equations, including Equations (2), (4), (9), (10), and (20), which describe the interconnection between the thermal properties, electrical properties, and mechanical properties of the system. Additionally, we introduce an environmental module for initial and boundary conditions, as follows:

4. Deterministic Solution

In this section, we describe the finite element discretization, considering a one-dimensional domain, and outline the deterministic solution process. The deterministic solution has two primary objectives: it acts as a black box for stochastic analysis for UQ and provides an interpretable basis for evaluating uncertainty propagation in a higher dimension.

4.1. Finite-Element Discretization

We consider the domain length, L = 200 m, and number of elements, N = 1000, for our one-dimensional problem. By multiplying the governing Equations (2), (4), (9), (10), and (20) with a test function w, integrating by parts, and considering the volume differential expressed as the cross-sectional area $A\left(x\right)$, we derive their corresponding weak forms as follows:

$${\int }_0^L-{(1-{\phi }_d)}^{2}YA\left(x\right){\displaystyle \frac{d{u}_d}{dx}}{\displaystyle \frac{dw}{dx}}dx+{\int }_0^L\gamma g_{c}A\left(x\right){\left({\displaystyle \frac{d{\phi }_d}{dx}}\right)}^2{\displaystyle \frac{dw}{dx}}dx+{\int }_0^{L}fA\left(x\right)wdx=0,$$

(34)

$$\begin{array}{c}\hfill {\int }_0^L-\gamma g_{c}A\left(x\right){\displaystyle \frac{d{\phi }_d}{dx}}{\displaystyle \frac{dw}{dx}}dx+{\int }_0^{L}A\left(x\right)\mathbb{H}wdx-{\int }_0^{L}A\mathbb{H}{\phi }_{d}wdx\\ \hfill -{\int }_0^L{\displaystyle \frac{g_{c}A\left(x\right)}{\gamma }}{\phi }_{d}wdx+{\int }_0^L{\displaystyle \frac{A\left(x\right)}{\gamma }}{\mathcal{F}}_{f}wdx=0,\end{array}$$

(35)

$${\int }_0^L{\dot{\mathcal{F}}}_{f}wA\left(x\right)dx={\int }_0^L{\displaystyle \frac{-\rho a(1-{\phi }_d)Y\left|\frac{d{u}_d}{dx}\right|(-{\phi }_d)}{\gamma }}{\displaystyle \frac{{\theta }_c}{{\theta }_0}}wA\left(x\right)dx,$$

(36)

$$\begin{array}{c}\hfill {\int }_0^L-{\kappa }_{m}A\left(x\right){\displaystyle \frac{d{\theta }_c}{dx}}{\displaystyle \frac{dw}{dx}}dx+{\int }_0^L{\sigma }_{E}A\left(x\right){\left({\displaystyle \frac{dV}{dx}}\right)}^{2}wdx\\ \hfill -{\int }_0^{L}h{\theta }_c{A}_s\left(x\right)wdx+{\int }_0^{L}h{\theta }_a{A}_s\left(x\right)wdx=0,\end{array}$$

(37)

$${\int }_0^L{\sigma }_E{\displaystyle \frac{dV}{dx}}{\displaystyle \frac{dw}{dx}}A\left(x\right)wdx=0,$$

(38)

For each element k, we use a linear approximation, as follows:

$$\begin{array}{cc}\hfill u_d^k& =N{\widehat{u}}_d^k,\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\phi }_d^k& =N{\widehat{\phi }}_d^k,\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {{\mathcal{F}}_f}^k& =N{\widehat{{\mathcal{F}}_f}}^k,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill {\theta }_c^k& =N{\widehat{\theta }}_c^k,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill V^k& =N{\widehat{V}}^k.\hfill \end{array}$$

(43)

We calculate the finite-element interpolation of spatial derivatives using linear combinations of the derivatives of shape functions:

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{d{u}_d}{dx}}\right)}^k& =B{\widehat{u}}_d^k,\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{d{\phi }_d}{dx}}\right)}^k& =B{\widehat{\phi }}_d^k,\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{d{\theta }_c}{dx}}\right)}^k& =B{\widehat{\theta }}_c^k,\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{dV}{dx}}\right)}^k& =B{\widehat{V}}^k,\hfill \end{array}$$

(47)

where N, B, ${\widehat{u}}_d^k$, ${\widehat{\phi }}_d^k$, ${\widehat{{\mathcal{F}}_f}}^k$, ${\widehat{\theta }}_c^k$, and ${\widehat{V}}^k$ are defined as follows:

$$\begin{array}{cc}\hfill N& =\left[\begin{array}{cc}{N}_1& N_2\end{array}\right],\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill B& =\left[\begin{array}{cc}{N}_{1,x}& N_{2,x}\end{array}\right],\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill {\widehat{u}}_d^k& =\left[\begin{array}{cc}{u}_1^k& u_2^k\end{array}\right],\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill {\widehat{\phi }}_d^k& =\left[\begin{array}{cc}{\phi }_1^k& {\phi }_2^k\end{array}\right],\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill {\widehat{{\mathcal{F}}_f}}^k& =\left[\begin{array}{cc}{\mathcal{F}}_1^k& {\mathcal{F}}_2^k\end{array}\right],\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill {\widehat{\theta }}_c^k& =\left[\begin{array}{cc}{\theta }_1^k& {\theta }_2^k\end{array}\right],\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill {\widehat{V}}^k& =\left[\begin{array}{cc}{V}_1^k& V_2^k,\end{array}\right].\hfill \end{array}$$

(54)

where $N_1$ and $N_2$ are linear interpolation functions.

We use a forward Euler method to evolve ${\mathcal{F}}_f$, resulting in the discretization of kth element as follows:

$$\begin{array}{cc}\hfill K_u{\widehat{u}}_d^k& =w_u+M{\widehat{f}}^k,\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill K_{\phi }{\widehat{\phi }}_d^k& =w_{\phi },\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill M{\widehat{\mathcal{F}}}_f^{{n+1}^k}& =M{\widehat{\mathcal{F}}}_f^{n^k}+\Delta t{w}_{\mathcal{F}},\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill K_{\theta }{\widehat{\theta }}_c^k& =w_{\theta },\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill K_V{\widehat{V}}^k& =0.\hfill \end{array}$$

(59)

where the discrete forms are as follows:

$$\begin{array}{cc}\hfill K_u& ={\int }_k{(1-N{\widehat{\phi }}_d^k)}^{2}YA\left(x\right)B^{T}Bdx,\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill w_u& ={\int }_k\gamma g_{c}A\left(x\right){\left(B{\widehat{\phi }}_d^k\right)}^{2}Bdx,\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill M& ={\int }_{k}A\left(x\right)N^{T}Ndx,\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill K_{\phi }& ={\int }_k\gamma g_{c}A\left(x\right)B^{T}Bdx+{\int }_k\mathbb{H}A\left(x\right)N^{T}Ndx+{\int }_k{\displaystyle \frac{g_{c}A\left(x\right)}{\gamma }}{N}^{T}Ndx,\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill w_{\phi }& ={\int }_k\mathbb{H}A\left(x\right)Ndx+{\int }_k{\displaystyle \frac{A\left(x\right)}{\gamma }}{N}^T{\widehat{\mathcal{F}}}_f^{n^k}Ndx,\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill K_{\theta }& ={\int }_k{\kappa }_{m}A\left(x\right)B^{T}Bdx+{\int }_{k}h{A}_s\left(x\right)N^{T}Ndx,\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill w_{\theta }& ={\int }_k{\sigma }_{E}A\left(x\right){\left(B{\widehat{V}}^k\right)}^{2}Ndx+{\int }_{k}h{A}_s\left(x\right){\theta }_{a}Ndx,\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill K_v& ={\int }_k{(1-N{\widehat{\phi }}_d^k)}^2{\sigma }_{E,T}A\left(x\right)B^{T}Bdx.\hfill \end{array}$$

(67)

We derive the global forms of these matrices and vectors using standard finite-element assembly procedures.

$$\begin{array}{cc}\hfill K_u{\widehat{u}}_d& =w_u+M\widehat{f},\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill K_{\phi }{\widehat{\phi }}_d& =w_{\phi },\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill M{\widehat{\mathcal{F}}}_f^{n+1}& =M{\widehat{\mathcal{F}}}_f^n+\Delta t{w}_{\mathcal{F}},\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill K_{\theta }{\widehat{\theta }}_c& =w_{\theta },\hfill \end{array}$$

(71)

$$\begin{array}{cc}\hfill K_V\widehat{V}& =0.\hfill \end{array}$$

(72)

The formulation supports a staggered solution scheme at each time step, implemented according to the following Algorithm 1:

Algorithm 1 Multiphysics framework solution.

Set initial conditions and boundary conditions. for Each $\Delta t$ do     Calculate the horizontal tensile load.     Solve the governing equation for displacements.     Update strain energy history.     Solve the equation for the damage field.     Update fatigue.     Solve the heat equation to update the temperature field.     Solve for the voltage field.

4.2. Results and Discussion

We consider an all-aluminum conductor (AAC) subjected to cyclic loading conditions under wind, temperature, and current. The details on wind and temperature loading conditions are provided in Section 2. The current loading is parameterized as follows:

$$I\left(t\right)=-I_{base}-I_{amp}(sin4\pi t),$$

(73)

where $I_{base}$ is the base current, set at 1500 A, and $I_{amp}$ is the amplitude, set at 100 A, corresponding to the allowable ampacity for a 40 mm aluminum conductor.

We represent the boundary conditions for displacement setting $u=0$ and $u=H$ at $x=0$ and $x=L$, respectively. At both ends, we specify ${\displaystyle \frac{d{\phi }_d}{dx}}=0$. Similarly, we apply boundary conditions for the current conservation equation setting $V=0$ and $V=J_d$ at $x=0$ and $x=L$, respectively.

In practice, all materials feature inherent imperfections, so we consider variable cross-section areas in the center of the domain to represent the cumulative effects of these multiple imperfections. We also consider a base scenario where the materials are assumed to have insignificant imperfections, serving as a comparison point with the significant initial damage cases. We define the cross-section area as follows:

$$A\left(x\right)=A_0\left(1-{\displaystyle \frac{1}{A_{\sigma }\sqrt{2\pi }}}exp\left({\displaystyle \frac{-{(x-L/2)}^2}{2A_{\sigma }^2}}\right)\right),$$

(74)

where $A_0$ denotes the cross-section area without any damage, while $A_{\sigma }$ indicates the spread-to-depth ratio at the center of the line, representing different levels of damage. Figure 9 illustrates various area profiles based on $A_{\sigma }$, as follows:

The material parameter values are given in

Table 3. Air, Material, and Geometry parameters.

Parameter	Value	Unit
Domain length L	200	m
Number of elements N	1000
Initial Tension $H_0$	40	kN
Diameter D	0.04	m
Young modulus Y	69	GPa
Damage layer width $\gamma$	0.02	m
Fracture energy $g_c$	10	kN/m
Density ${\rho }_m$	2700	kg/m3
Aging coefficient a	$1\times {10}^{-10}$	m5/(y kg)
Thermal conductivity ${\kappa }_m$	237	W/(m K)
Electrical conductivity ${\sigma }_{E,0}$	$3.77\times {10}^7$	S/m
Temperature coefficient ${\alpha }_l$	$3.9\times {10}^{-3}$	K−1
Density of air ${\rho }_{air}$	1.225	kg/m3
Kinematic viscosity of air ${\nu }_{air}$	$15\times {10}^{-6}$	m2/s
Thermal conductivity of air ${\kappa }_{air}$	0.0295	W/(m K)
Prandtl $Pr$	0.71

. The simulations are considered for 60 years with a 6000 run period for which the time-step, $\Delta t$ = 0.01. Aluminum begins to anneal at temperatures exceeding 366 K [3,58], and rupture occurs when the temperature surpasses 373 K [59]. Consequently, we establish the maximum temperature limit as ${\theta }_{\lim }=373$ K (100 °C). Once this limit is reached, all simulations are stopped.

We begin our analysis by examining the evolution of field quantities in the Texas scenario, using the parameters in Table 3. Every five years, we plot the evolution of field variables along the length of the transmission line as shown in Figure 10. We note that damage typically begins and builds up in areas with smaller cross-sections, causing increased temperatures and notable disturbances in the voltage fields. As damage and temperature rise, the voltage drop along the line also increases over time due to higher electrical resistance in the line.

Next, we study the effect of $A_{\sigma }$ on the maximum values of field variables. Figure 11 shows that the presence of damage affects the evolution of field variables at the center of the line.

Moreover, we compare the life span of transmission lines for each specific scenario varying the cross-sectional area of damage $A_{\sigma }$ under respective loading conditions in Figure 12. In the Texas scenario, we observe the life span of the transmission line reduces from 51 years to 44 years with moderate damage. However, with severe damage, the life span reduces to 28 years, indicating how detrimental initial damage could be to the longevity of the material. Similarly, in the California scenario, the lifespan of the transmission line reduces from 52 years to 26 years under severe damage. This reduction is due to the state’s high temperatures and low wind speeds, significantly lowering the convective cooling effect. Consequently, in the absence of the cooling effect, the accumulation of damage raises the material temperature, leading to overheating and early failure. In the case of the Michigan scenario, the state shows the most fluctuating temperature and wind over a year, which progresses the damage, raising the temperature of the material and leading to early failure under insignificant damage conditions among the four states. However, due to the high wind speed and low temperature enhancing convective cooling during some seasons, the life span under severe damage is higher than the California scenario. The Florida scenario shows high temperatures over a year, similar to Texas and California. However, the high wind speed in the state enhances convective cooling, reducing the early failure of the transmission line. Overall, every scenario shows a significant reduction in the life span of the transmission line under severe damage.

5. Stochastic Solution

In this section, we address the challenge of performing UQ, SA, and the probability of failure ($P_f$) on the model treated as a black box. We focus on the conductor temperature as our quantity of interest (QoI). The non-intrusive nature of our approach is particularly beneficial, allowing us to apply the same methods of deterministic solutions to stochastic problems without altering the model’s governing equations.

We employ the PCM for three main purposes. First, we use it to compute the QoI. Second, we conduct global sensitivity analysis by calculating Sobol sensitivity indices ($S_i$). These indices highlight the relative importance of each model parameter by quantifying the contribution to the variance of the QoI. Lastly, we use PCM to facilitate the computation of the probability of failure in a straightforward fashion by calculating the expectation of the Bernoulli random variable.

5.1. Uncertainty Quantification

We followed the approach outlined in [27] to conduct UQ using PCM. We assume $({\Omega }_s,\mathcal{G},\mathbb{P})$ as a complete probability space, where ${\Omega }_s$ denotes the space of outcomes $\omega$, $\mathcal{G}$ is the $\sigma -$ algebra, and $\mathbb{P}$ is a probability measure mapping $\mathcal{G}$ to the interval [0, 1]. In our model, we treat the material and loading parameters as random variables with two different sets, ${\xi }_m\left(\omega \right)$ and ${\xi }_l\left(\omega \right)$. We model our QoI temperature as random variables and denoted random parameters as $\xi =\xi \left(\omega \right)$.

Our QoI, denoted as $Q_{oI}$, involves calculating the expectation $\mathbb{E}\left[Q_{oI}(x,t;\xi )\right]$ as

$$\mathbb{E}\left[Q_{oI}(x,t;\xi )\right]={\int }_a^b{Q}_{oI}(x,t;\xi )\rho \left(\xi \right)d\xi ,$$

(75)

where the probability density function (PDF) of $\xi$, denoted as $\rho \left(\xi \right)$, is evaluated by mapping the physical parametric space to the standard space $[-1,1]$, followed by integration using the Gauss quadrature. This transformation allows the integral to be expressed within this standard interval as follows:

$$\mathbb{E}\left[Q_{oI}(x,t;\xi )\right]={\int }_{-1}^1{Q}_{oI}(x,t;\xi \left(\eta \right))\rho \left(\xi \left(\eta \right)\right)Jd\xi \left(\eta \right),$$

(76)

where $J=d\xi /d\eta$ represents the Jacobian of the transformation. To approximate the expectation, we use polynomial interpolation in the stochastic space, expressed as ${\widehat{Q}}_{oI}(x,t;\xi )$, which approximates the exact solution as follows:

$$\mathbb{E}\left[Q_{oI}(x,t;\xi )\right]\approx {\int }_{-1}^1{\widehat{Q}}_{oI}(x,t;\xi \left(\eta \right))\rho \left(\xi \left(\eta \right)\right)Jd\xi \left(\eta \right).$$

(77)

We use Lagrange polynomials $L_i\left(\xi \right)$ to interpolate the solution in the stochastic space, as follows:

$${\widehat{Q}}_{oI}(x,t;\xi )=\sum _{i=1}^I{Q}_{oI}(x,t;{\xi }_i)L_i\left(\xi \right),$$

(78)

which satisfies the Kronecker delta property at the interpolation points as follows:

$$L_i\left({\xi }_j\right)={\delta }_{ij}.$$

(79)

We substitute the polynomial approximation from Equation (78) into Equation (77) and use the quadrature rule to approximate the integral and compute the expectation as follows:

$$\mathbb{E}\left[Q_{oI}(x,t;\xi )\right]\approx \sum _{m=1}^n{w}_m\rho \left(\xi \left(\eta \right)\right)J\sum _{i=1}^I{Q}_{oI}(x,t;\xi \left(\eta \right))L_i\left(\xi \left(\eta \right)\right),$$

(80)

where ${\eta }_m$ and $w_m$ are coordinates and weights, respectively, for each integration point $m=1,2,\cdots ,n$. We choose the same collocation and integration points using Equation (79) and simplify Equation (80) as follows:

$$\mathbb{E}\left[Q_{oI}(x,t;\xi )\right]=\sum _{m=1}^n{w}_m\rho \left({\xi }_m\left({\eta }_m\right)\right)J{Q}_{oI}(x,t;{\xi }_m\left({\eta }_m\right)).$$

(81)

We use ${\xi }_m\left({\eta }_m\right)=a+{\displaystyle \frac{(b-a)}{2}}({\eta }_m+1)$ as linear affine mapping that calculates the Jacobian as $J=(b-a)/2$. Finally, we approximate the integration and rewrite as a summation over the collocation points, assuming a uniform distribution for the parameters throughout the interval $[a,b]$, with $\rho \left(\xi \right)=1/(b-a)$. The expectation is as follows:

$$\mathbb{E}\left[Q_{oI}(x,t;\xi )\right]={\displaystyle \frac{1}{2}}\sum _{m=1}^n{w}_m{Q}_{oI}(x,t;{\xi }_m).$$

(82)

Similar to the MC method, the standard deviation is computed as follows:

$$\sigma \left[Q_{oI}(x,t;\xi )\right]=\sqrt{{\displaystyle \frac{1}{2}}\sum _{m=1}^n{w}_m{\left(Q_{oI}(x,t;{\xi }_m)-\mathbb{E}\left[Q_{oI}(x,t;\xi )\right]\right)}^2}.$$

(83)

For the higher dimensions, Equation (75) includes additional integration that reduces to the following:

$$\begin{array}{cc}\hfill \mathbb{E}\left[Q_{oI}(x,t;{\xi }^1,\cdots ,{\xi }^k)\right]& ={\mathbb{E}}_{PCM}\left[Q_{oI}(x,t;{\xi }^1,\dots ,{\xi }^k)\right]\hfill \\ & \approx \sum _{m=1}^n\dots \sum _{l=1}^L{w}_m\dots w_l\rho \left({\xi }_m\right)\dots \rho \left({\xi }_l\right)J_m\dots J_l{Q}_{oI}(x,t;{\xi }_m^1,\dots ,{\xi }_l^k)\hfill \end{array}$$

(84)

In the expanded version of PCM for higher dimensions, we have k summations corresponding to each dimension in the random space. For simplicity, we represent the expectation as $\mathbb{E}\left[Q_{oI}(x,t;{\xi }^1,,\dots ,,{\xi }^k)\right]=\mathbb{E}\left[Q_{oI}\right]$. We then formulate the standard deviation as follows:

$$\begin{array}{c}\sigma \left[Q_{oI}(x,t;{\xi }^1,\dots ,{\xi }^k)\right]={\sigma }_{PCM}\left[Q_{oI}(x,t;{\xi }^1,\dots ,{\xi }^k)\right]\hfill \\ \approx \sqrt{\sum _{m=1}^n\dots \sum _{l=1}^L{w}_m\dots w_l\rho \left({\xi }_m\right)\dots \rho \left({\xi }_l\right)J_m\dots J_l{\left(Q_{oI}(x,t;{\xi }_m^1,\dots ,{\xi }_l^k)-\mathbb{E}\left[Q_{oI}\right]\right)}^2}.\hfill \end{array}$$

(85)

Here, we assume discretization in the parametric space is isotropic, and the random variables are mutually independent. Additionally, PCM is sufficient for the model in this study, as it involves fewer than six dimensions. However, this approach becomes inefficient for models with more than six dimensions due to the exponential increase in the number of simulations required by the tensor product. The problem of dimensionality can be reduced using Smolyak sparse grids [46], Principal Component Analysis [60], active subspaces methods [47], and low-rank approximations [61].

5.2. Sensitivity Analysis

We examine the sensitivity analysis using Sobol indices [62]. For details on the derivation, we refer to the work by Saltelli et al. [63]. In our global sensitivity analysis, each parameter is denoted as ${\xi }^u$, where $u=1,2,...,v$. We evaluate the influence of ${\xi }^u$ on the variance of the QoI as follows:

$$V_{{\xi }^u}\left({\mathbb{E}}_{{\tilde{\xi }}^u}\left(Q_{oI}\right|{\xi }^u)\right)$$

(86)

where ${\tilde{\xi }}^u$ represents all possible values of the random parameters except for ${\xi }^u$, which is held fixed. Equation (86) describes the process of computing the expected value of $Q_{oI}$ with ${\xi }^u$ fixed. Based on the law of total variance, we have the following:

$$V_{{\xi }^u}\left({\mathbb{E}}_{{\tilde{\xi }}^u}\left(Q_{oI}\right|{\xi }^u)\right)+{\mathbb{E}}_{{\tilde{\xi }}^u}\left(V_{{\tilde{\xi }}^u}\left(Q_{oI}\right|{\xi }^u)\right)=V\left(Q_{oI}\right)$$

(87)

We normalize Equation (87) to obtain the first-order $S_i$. This index quantifies the impact of the random variable ${\xi }^u$ on the total variance and is calculated as follows:

$$S_i={\displaystyle \frac{V_{{\xi }^u}\left({\mathbb{E}}_{{\tilde{\xi }}^u}\left(U\right|{\xi }^u)\right)}{V\left(U\right)}}$$

(88)

$S_i$ only quantifies the first-order effects of the variable ${\xi }^u$ on the variance. Following normalization, the sum of all $S_i$ values is less than 1, with the remaining portion representing interactions among the parameters, which are not considered in this paper.

5.3. Probability of Failure

In this section, we calculate the $P_f$ over time. Traditionally, MC methods are employed in reliability analysis to compute the probability of failure [64]. Stochastic collocation methods are used for similar analysis; however, the method primarily yields the moments of a limit state function $g(R,S)$. To compute $P_f$, these moments must be transformed into a PDF, which can be achieved using various techniques such as the method of moments [65,66], polynomial chaos [67,68], Gaussian transformations [69], or entropy optimization methods [70]. However, these processes add a layer of complexity.

In this study, we present an alternative way that utilizes the efficiency of the PCM to compute the probability of failure ($P_f$). Rather than calculating and approximating PDF to determine $P(g<0)$, we simplify the process by transforming g into a Bernoulli random variable $h_B$ with coefficient $p_h$.

The transformation of g to $h_B$ is defined as follows:

$$h_B=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}g\ge 0,\hfill \\ 1,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(89)

In single realization, we obtain $h_B$, which remains zero until the point where ${\theta }_{\max }$ exceeds the threshold ${\theta }_{\lim }$, at which $h_B$ becomes one and remains at this value. This represents a step function for each realization.

At a fixed time-step, considering the expectation of $h_B$, due to the smoothness of the QoI, $h_B$ can be assumed as a real value between 0 and 1 which reflects the probability of the maximum temperature exceeding the limit up to that time.

5.4. Results and Discussion

In the deterministic case, we understood how varying initial damage under cyclic loading conditions influences the lifespan of the transmission line. Now, we shift our focus to understanding the impact of uncertainty in the parameters on maximum temperature. For this analysis, we model the parametric uncertainty with a uniform distribution, where each parameter varies by $10\%$ around its mean value.

5.4.1. Texas Scenario

We first consider the Texas scenario as our baseline for conducting preliminary analyses on the uncertainties associated with parametric material properties and input loading conditions. Our primary goal is to analyze UQ and perform SA on factors influencing the maximum temperature. We first focus on material parameters, represented by the set ${\xi }_m\left(\omega \right)=\{A_{\sigma }\left(\omega \right),\gamma \left(\omega \right),g_c\left(\omega \right),a\left(\omega \right)\}$. These parameters are selected due to their inherent measurement inaccuracies or their assumptions within our modeling. We assume the rest of the material parameters to be deterministic. Following the analysis of material parameters, we then examine the uncertainty in loading conditions, denoted by the set ${\xi }_l\left(\omega \right)={\theta }_{base}\left(\omega \right),w_{base}\left(\omega \right),I_{base}\left(\omega \right),I_{amp}\left(\omega \right)$. From both analyses using 5 PCM points per dimension, we identify two material parameters and three loading parameters that are most influential according to global SA results.

We begin by examining the uncertainty in the material parameters from the set ${\xi }_m\left(\omega \right)$. We analyze the expectation and standard deviation of the temperature field over time, as presented in Figure 13. Similarly to the deterministic case, the location at the conductor with initial damage has the maximum temperature and standard deviation. Considering this, we plot the evolution of the expected maximum temperature and the standard deviation at the center of the transmission line, where both the maximum temperature and standard deviation increased over time, as shown in Figure 14. We plot the evolution of the expected maximum temperature and the standard deviation similarly for loading parametric set ${\xi }_l\left(\omega \right)$ in Figure 15.

We use Equation (88) to calculate the Sobol indices $S_i$ and shown in Figure 16 for material parameters set ${\xi }_m\left(\omega \right)$ and loading parameters set ${\xi }_l\left(\omega \right)$. Initially, $A_{\sigma }$ has a significant impact, but it is later surpassed by $g_c$ and a. The parameter $I_{base}$, among the loading parameters, sets ${\xi }_l\left(\omega \right)$ from the initial state and emerges as the most crucial factor influencing the uncertainty of ${\theta }_{max}$ due to its direct impact on Joule heating. Meanwhile, the wind base parameter, $w_{base}$, and the temperature base parameter also play significant roles due to the interplay between Joule heating and convective cooling:

We combine the two most influential material parameters with the three most influential loading parameters to form a new set, ${\xi }_1\left(\omega \right)=g_c\left(\omega \right),a\left(\omega \right),{\theta }_{base}\left(\omega \right),w_{base}\left(\omega \right),I_{base}\left(\omega \right)$, and conduct a sensitivity analysis for Texas. Figure 17 shows that initially, the loading conditions have a more substantial impact than material parameters, but over time, as aging effects become more significant, the relative importance of material parameters increases.

5.4.2. California, Michigan, and Florida

In California, Michigan, and Florida scenarios, we combine the two most influential material parameters and the three most influential loading parameters to obtain a new parametric set, $\xi \left(\omega \right)=\{g_c\left(\omega \right),a\left(\omega \right),{\theta }_{base}\left(\omega \right),w_{base}\left(\omega \right),I_{base}\left(\omega \right)\}$, to perform the SA.

In the California scenario, when the time series for maximum temperature is truncated at the earliest occurrence of failure, the analysis reveals that failures happen within a span shorter than five years. However, the pattern in Figure 18 aligns closely with observations from the Texas Scenario. The consistency in the pattern highlights the dominant influence of the current base parameter, $I_{base}$, which is intensified by low wind speeds that reduce convective cooling. As a result, the conductor temperature rapidly increases to critical thresholds, driven by the Joules heating in the absence of effective cooling:

In the Michigan scenario, shown in Figure 19, $I_{base}$ initially has the most significant impact on the failure of the transmission line. Over time, however, $g_c$ becomes increasingly influential, eventually surpassing the impact of $I_{base}$. Historical data from Michigan shows that the state experiences higher wind speeds compared to the other states analyzed. Although the wind has a convective effect, the impact of increased mechanical load on the cable is more significant, accelerating the transmission line’s aging:

In the Florida scenario, as illustrated in Figure 20, $I_{base}$ initially shows the most significant influence on the transmission line’s failure. However, over time, the fracture energy parameter, $g_c$, becomes more dominant. Historical data considered in this study show that Florida experiences high temperatures throughout the year. Despite higher wind speeds, the effectiveness of convective cooling is diminished due to the high ambient temperatures.

5.4.3. Probability of Failure

In this section, we first calculate the expected value of h, with the application of PCM considering $n=5$. The results are displayed in Figure 21 for the reference mean parameter values from the parameter sets ${\xi }_1\left(\omega \right)$, ${\xi }_2\left(\omega \right)$, ${\xi }_3\left(\omega \right)$, and ${\xi }_4\left(\omega \right)$. Under the consideration of four specific states, the probability of failure curve for California initially shows the early higher possibility of occurrence of failure. In contrast, the Florida scenario shows the failure curve on the right shift, indicating a lower possibility of failure among the four states.

Finally, we analyze the impact of varying mean values of five influential parameters ${\xi }_1\left(\omega \right)=g_c\left(\omega \right),a\left(\omega \right),{\theta }_{base}\left(\omega \right),w_{base}\left(\omega \right),I_{base}\left(\omega \right)$ under three different levels of $A_{\sigma }$ in the time-series probability of failure ($P_f$) for each scenario. The comparison includes a baseline scenario representing a transmission line with minimal initial damage. Unlike previous sensitivity analyses of parameters, which provided insights into levels of importance relative to their uncertainty, this section focuses on the effects of changing the baseline values of influential parameters on the $P_f$ curves under minimal, moderate, and severe initial damage.

In each scenario, we observe three distinct curve shapes representing different levels of damage, as depicted in Figure 22. The analysis shows that in the presence of moderate initial damage, the chance of failure increases from $20\%$ to $60\%$ after 45 years in the Texas, California, and Michigan scenarios. In Florida, the probability of failure increases from $20\%$ to $60\%$ after 45 years. In severe initial damage, the probability of failure reaches $100\%$ well before the failure initiates in a minimal damage case. All the states show the life span of the transmission lines to be around 35 to 40 years under severe damage. This indicates how detrimental the initial damage could be to the life span of the transmission line. Among all the four scenarios, the chance of failure was initiated early in the case of California, which has low wind and high ambient temperature, reducing the impact of convective cooling. Also, the plots reveal that the probability of failure is around $90\%$ in 60 years, indicating a $10\%$ chance the transmission line will last more than 60 years. However, in the case of Michigan, as it has higher fluctuations in wind and temperature, the chance of failure is maximum at 60 years compared to the three states.

6. Verification and Analysis of the Model

As the model involves the coupling of several governing equations, obtaining an analytical solution is not feasible. Therefore, we perform convergence analysis using a refined PCM solution as our reference:

$$ϵ={\displaystyle \frac{\parallel \theta -{\theta }_{ref}{\parallel }_2}{\parallel {\theta }_{ref}{\parallel }_2}}.$$

(90)

In our convergence study, we defined the problem as a 1D problem, considering the most influential parameter $I_{base}$ as the only model parameter. We considered 100 collocation points as the reference solution and calculated the relative errors until the earliest failure, considering the Texas Scenario as a reference. The comparison is graphically represented in Figure 23 to highlight the accuracy of PCM relative to Monte Carlo methods:

The probabilistic collocation method achieved significantly lower errors, approximately two orders of magnitude less, compared to 10,000 Monte Carlo simulations using just 5 collocation points.

7. Conclusions

We developed an integrated model to study the reliability of transmission lines across four specific states in the US: Texas, California, Michigan, and Florida. We considered the historical wind and temperature data of each state. The model integrates the fatigue and damage phase-field model with a thermal model that accounts for Joule heating and convective cooling, as well as an electrical model that addresses how accumulated damage and temperature-induced resistance lead to voltage drops along the transmission line.

We used the finite element method to solve all the governing equations. We studied four different states of the US, Texas, California, Michigan, and Florida, to understand the long-term behavior of the overhead power lines in the presence of minimal, moderate, and severe damage. We used discrete Fourier analysis to obtain the continuous loading condition from the discrete historical data of wind and temperature. To reduce the complexity, we parameterized the current loading. Subsequently, we implemented the PCM for UQ, SA, and probability of failure assessments.

The deterministic solution revealed how the temperature in the conductor under cyclic wind, air temperature, and current loading in the presence of initial damage affects the overall longevity of the transmission lines. Even the moderate initial damage significantly reduced the material’s life expectancy. In the global sensitivity analysis using PCM, we identified that $I_{base}$ initially had the most significant effect, which was later surpassed by $g_c$ over time. Additionally, exploring the further application of PCM, we studied the probability of failure for each specific scenario, varying the most influential parameters for three damage cases. The analysis revealed that the chance of failure was significantly increased even with moderate initial damage.

Overall, this research provides critical insights into how initial damage significantly reduces the lifespan of transmission lines—an aspect often overlooked in previous studies. By highlighting key parameters that most influence transmission line reliability, this predictive model offers direct applications for improving transmission line design and enhancing resilience. Furthermore, by incorporating realistic environmental data, the model can be adapted to various regions, making it a flexible and practical tool for assessing transmission line reliability across diverse climates and environmental conditions.

However, the model incorporates simplifying assumptions such as static conditions and the one-dimensional domain of the cable, which could be further explored in the future. Even the integration of real current data for specific states, along with the effects of solar radiation and solar heat gain, could enhance the predictive accuracy of failure probabilities. All these modifications shift the dependency from abstract stochastic process modeling to a more robust physics-based approach, significantly improving the reliability assessments of the transmission line.