Modeling and Analysis of Progressive Ice Shedding along a Transmission Line during Thermal De-Icing

Abstract

Progressive ice shedding (PIS) along transmission lines is a common type of ice shedding during thermal de-icing that requires investigation to ensure the security of transmission lines. In current research, PIS is commonly analyzed using a constant speed for ice detaching from the conductor, which is not accurate for PIS simulation. Therefore, a mechanical model of PIS is established in this study to analyze PIS during thermal de-icing. First, an ice detachment model during thermal de-icing is built to determine the detachment times of the initial ice and remaining ice. Then, a two-node isoparametric truss element is employed to derive the static and dynamic equilibrium equations of an iced conductor to simulate the dynamic response of PIS. Relative to commercial software, these equations can easily accommodate the changing mass of ice with the flow of melted water. The dynamic equilibrium equations are then solved using the ice detachment model to obtain the dynamic response of PIS. Finally, small-scale and full-scale experimental results are employed to verify the proposed method. The simulation results show that the results of the proposed method are more consistent with the experimental results than are the results of existing methods that assume a constant propagation speed. The proposed method can be further applied to optimize transmission line designs and evaluate the application of thermal de-icing devices.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Ice shedding is one of the major sources of iced transmission line faults in cold regions. The statistics of historical ice disasters, such as the ice storm in North America in 1998 and the ice disaster in southern China in 2008, show that conductor vibration caused by ice shedding can result in flashover, fire burn, burnout and other electrical accidents, insulator rupture, cable breakage and tower deformation, and collapse [1]. The electrical and mechanical faults of transmission lines threaten the reliability of power system operation and the security of transmission lines. To enhance the design of transmission lines and ensure the security of power systems, analysis of the dynamic response of ice shedding is necessary.

Current research on the dynamic response of ice shedding from transmission lines generally focuses on ice shedding scenarios with a fixed amount of ice shedding. Various ice shedding scenarios have been analyzed, including ice shedding from a single-span transmission line [2, 3], a continuous-span line [4], an overhead ground wire [5], a tower-line system [6], a transmission line with bundled conductors and spacers [3, 7], and a high-voltage overhead transmission line [8, 9].

This research assumes that, after the initial ice shedding, the remaining ice does not detach from the vibrating conductor. However, the remaining ice may fracture into small fragments from conductor vibration due to its plastic strain [10] and result in progressive ice shedding (PIS) along the transmission line [5]. When transmission lines are exposed to thermal de-icing [11], the maximum allowed plastic strain of the ice decreases as the inner ice melts during thermal de-icing, and the ice detaches from the conductor due to the propagation of transverse waves along the transmission line. This phenomenon was confirmed through experiments on thermal de-icing [12]. Consequently, PIS is a common type of ice shedding during thermal de-icing. Therefore, it is significant to analyze the dynamic response of PIS during thermal de-icing for the security of transmission lines and the application of thermal de-icing.

Analyzing PIS during thermal de-icing involves two challenges: the criterion of ice detachment from the conductor and the dynamic model for analyzing the PIS dynamic response during thermal de-icing. For the ice detachment criterion (IDC), in current research, the time of ice detachment from the conductor is determined by a constant propagation speed of the transverse wave along the transmission line [13]. Nevertheless, the constant propagation speed is variable for various ice and conductor parameters, and the ice may not break apart with the propagation of transverse waves. Consequently, the IDC for a constant ice shedding speed is inaccurate for analyzing PIS. The ice detaches from the conductor due to interaction between the inner and external forces in the ice [5, 14, 15]. The maximum bending stress and maximum effective plastic strain of ice were employed as the IDC in [5, 14], and the adhesive force and cohesive force on the ice were adopted to calculate the IDC in [15]. These IDCs are effective for solid ice, but the ice in thermal de-icing is hollow ice with an air gap that involves a different IDC than that of solid ice. Although the model to calculate the air gap by thermal de-icing was proposed in [16], the IDC in thermal de-icing must be further researched.

The current research considers relatively simple ice shedding scenarios, such as a fixed amount of ice shedding from conductor. It is not easy for existing methods to simulate progressive ice melting and the resultant complex dynamics, such as the mass variation and dynamic force variation on the remaining ice caused by conductor fluctuation [2–9]. Moreover, the ability to simultaneously analyze many ice shedding scenarios with various transmission line parameters, ice thicknesses, and ice shedding positions also needs to be improved.

To consider the complicated PIS during thermal de-icing, this paper proposes a new modeling method, which is able to deal with the progressive dynamics of PIS. This method is based on the finite element mathematical model which is widely used in the research of suspended cable structures in mechanics [17]. The finite element method can consider the mass variation of ice and simulate the dynamic force variation in the ice for PIS. The proposed method first analyzes the characteristics of PIS during thermal de-icing. Compared with existing work, the detachment time of the initial ice shedding and the dynamic detaching model of the remaining ice are described in detail. Then, these characteristics are integrated into a finite element model to develop a more practical and exact PIS model. The proposed method has four advantages: (1) the mass variation of ice with thermal de-icing and the fluctuating conductor is considered; (2) the dynamic force on the ice is calculated dynamically; (3) the time of ice detachment from the conductor is more accurate because the IDC is obtained from the dynamic force of the ice; and (4) various PIS scenarios can be built flexibly because of the transparent modeling and solution process.

The paper is organized as follows. Section 2 presents the ice detachment model for thermal de-icing. Section 3 describes the static and dynamic equilibrium equations of a conductor based on a finite-element mathematical model. Section 4 explains the solution of the dynamic equilibrium equations considering the integration of ice detachment and conductor vibration. Case studies are presented and discussed in Section 5, followed by conclusions.

2. Characteristic Analysis of PIS during Thermal De-Icing

In this section, two types of ice detachment are introduced by analyzing the dynamic process of PIS. Then, two ice detachment models for the initial ice and remaining ice are proposed based on the ice melting model and the mechanical model of ice.

2.1. Ice Detachment Process in Thermal De-Icing

Due to the torque of the conductor and the gravity of the ice and wind, accreted ice usually has an irregular shape, which complicates the analysis of ice shedding. To analyze the PIS, we assume that the cross section of the ice is eccentric and round and that the conductor is farther from the center of the round from the middle of the conductor to the suspension point [12, 15].

When the de-icing current flows through the conductor, the ice at different positions on the conductor melts simultaneously. However, because the conductor near the suspension point is closest to the upper surface of the ice, the ice near the suspension point detaches from the conductor first and produces the initial ice shedding and conductor vibration. At this time, the ice remaining on the conductor vibrates with the conductor. When the external force on the remaining ice is greater than the inner force of the ice, the ice detaches from the conductor, resulting in PIS. Therefore, the detachment models of the ice near the suspension point and the ice remaining on the conductor are the foundation for modeling PIS.

2.2. Classic Thermal De-Icing Model

Although the sun and wind can melt ice, thermal de-icing is the preferred solution in severe weather. Additionally, the duration of thermal de-icing is between 0.5 and 3 hours [18]. Relative to the energy produced by thermal de-icing, ice melting by severe weather can be ignored. Therefore, the following assumptions are adopted in the study: (1) external conditions will not result in ice melting or ice accumulation and (2) the conductor resistance is constant, and the ice melting current evenly heats the conductor.

Based on the above assumptions, the heat conduction process of thermal de-icing can be described as follows. An electric current flowing through the transmission conductor produces Joule heat. Part of the heat warms the ice layer, conductor, and air gap, and part of the heat (namely, the latent heat) melts the ice, and the remaining heat reaches the outer surface of the ice layer and is dissipated through convection and radiation [11, 18]. This heat conductor process can be formulated as $$\begin{array}{}\text{(1)}& Q_{\mathrm{J}}=R_{\mathrm{0}}{I}^{\mathrm{2}}t=Q_{\mathrm{1}}+Q_{\mathrm{2}}+Q_{\mathrm{3}}\end{array}$$ where $Q_{\mathrm{J}}$ is to the Joule heat produced by the ice melting current per unit length, J/m; R ₀ is the conductor resistance per unit length, Ω/m; I is the ice melting current, A; t is the duration of thermal de-icing, s; Q ₁ is the heat dissipated through convection and radiation; Q ₂ is the latent heat for melting ice; and Q ₃ is the heat that warms the ice layer (the heat to warm the conductor and air gap is ignored) [11]. $$\begin{array}{}\text{(2)}& Q_{\mathrm{1}}=\mathrm{2}\pi r_{\mathrm{i}\mathrm{c}\mathrm{e}}h\left(T_{\mathrm{o}\mathrm{i}}-T_{\mathrm{a}\mathrm{i}\mathrm{r}}\right)t\text{(3)}& Q_{\mathrm{2}}={\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{L}_{\mathrm{f}}{V}_{\mathrm{m}}\text{(4)}& Q_{\mathrm{3}}={\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{V}_{\mathrm{i}\mathrm{c}\mathrm{e}}{C}_{\mathrm{i}\mathrm{c}\mathrm{e}}{T}_{\mathrm{i}\mathrm{c}\mathrm{e}}\end{array}$$ Here, $r_{\mathrm{i}\mathrm{c}\mathrm{e}}$ is the radius of the iced conductor, m; h is the heat-exchange coefficient [18], W/(m²·K); $T_{\mathrm{o}\mathrm{i}}$ is the outer surface temperature of the ice layer [18], °C; $T_{\mathrm{a}\mathrm{i}\mathrm{r}}$ is the ambient temperature, °C; ${\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}$ is the density of ice, kg/m³; $L_{\mathrm{f}}$ is the latent heat for melting ice ( $L_{\mathrm{f}}$ = 335,000 J/kg); $V_{\mathrm{m}}$ is the volume of melted ice per unit length, m²; $V_{\mathrm{i}\mathrm{c}\mathrm{e}}$ is the volume of the ice layer per unit length, m²; $C_{\mathrm{i}\mathrm{c}\mathrm{e}}$ is the specific heat of ice, J/(kg·°C); and $T_{\mathrm{i}\mathrm{c}\mathrm{e}}$ is the temperature of the ice layer, which can be simplified as ( $T_{\mathrm{o}\mathrm{i}}$ /2- $T_{\mathrm{a}\mathrm{i}\mathrm{r}}$ ) [18], °C.

The entire equation can be expressed as $$\begin{array}{}\text{(5)}& \left[R_{\mathrm{0}}{I}^{\mathrm{2}}-\mathrm{2}\pi r_{\mathrm{i}\mathrm{c}\mathrm{e}}h\left(T_i-T_a\right)\right]t={\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{L}_{\mathrm{f}}{V}_{\mathrm{m}}+{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{V}_{\mathrm{i}\mathrm{c}\mathrm{e}}{C}_{\mathrm{i}\mathrm{c}\mathrm{e}}\left(\frac{T_{\mathrm{o}\mathrm{i}}}{\mathrm{2}}-T_{\mathrm{a}\mathrm{i}\mathrm{r}}\right)\end{array}$$ where R ₀ and $r_{\mathrm{i}\mathrm{c}\mathrm{e}}$ are the conductor parameters; $d_{\mathrm{i}\mathrm{c}\mathrm{e}}$ , $L_{\mathrm{f}}$ , $C_{\mathrm{i}\mathrm{c}\mathrm{e}}$ , $T_{\mathrm{o}\mathrm{i}}$ , and $V_{\mathrm{i}\mathrm{c}\mathrm{e}}$ are the ice parameters; h and $T_{\mathrm{a}\mathrm{i}\mathrm{r}}$ are environmental parameters; and I, t, and $V_{\mathrm{m}}$ are the parameters associated with thermal de-icing. In the process of thermal de-icing, I and $V_{\mathrm{m}}$ are determined through an initial calculation. In addition, the duration t of thermal de-icing is determined to calculate the mechanical parameters of the ice remaining on the conductor. Therefore, based on formula (5), we can build the detachment models of the initial ice shedding and the ice remaining on the conductor.

2.3. The Detaching Time of Initial Ice Shedding

The cross section of the eccentric iced conductor before and after ice melting is illustrated in Figure 1, in which the shadow area denotes melted ice (the air gap).

[figures omitted; refer to PDF]

When the elliptical air gap is tangent to the outer surface of both the conductor and the ice layer, the ice hangs from the conductor [18], which is the detachment criterion of the initial ice shedding in the static state. Therefore, at the moment of detachment, the volume of the air gap per length [18] can be expressed as $$\begin{array}{}\text{(6)}& V_m=\mathrm{0.5}\pi {\left(r_{\mathrm{c}}+\frac{d_{\mathrm{u}}}{\mathrm{2}}\right)}^{\mathrm{3}/\mathrm{2}}\left({r_{\mathrm{c}}}^{\mathrm{1}/\mathrm{2}}+{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{1}/\mathrm{2}}\right)-\pi {r_{\mathrm{c}}}^{\mathrm{2}}\end{array}$$ where $r_{\mathrm{c}}$ is the radius of the conductor, m, and $d_{\mathrm{u}}$ denotes the minimum distance between the upper surface of the conductor and the ice layer, m. For eccentric iced conductor, this distance can be calculated as $$\begin{array}{}\text{(7)}& d_{\mathrm{u}}=r_{\mathrm{i}\mathrm{c}\mathrm{e}}-r_{\mathrm{c}}-d_{\mathrm{e}}\end{array}$$ where $d_{\mathrm{e}}$ is the distance between the center of the iced conductor and the center of the cable before thermal de-icing, m.

By substituting formula (6) into formula (5), the duration of thermal de-icing of the iced layer is obtained. $$\begin{array}{c}\text{(8)}& t_{\mathrm{d}}=\frac{\left[{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{L}_{\mathrm{f}}{V}_{\mathrm{m}}+{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{V}_{\mathrm{i}\mathrm{c}\mathrm{e}}{C}_{\mathrm{i}\mathrm{c}\mathrm{e}}\left(T_{\mathrm{o}\mathrm{i}}/\mathrm{2}-T_{\mathrm{a}\mathrm{i}\mathrm{r}}\right)\right]}{R_{\mathrm{0}}{I}^{\mathrm{2}}-\mathrm{2}\pi r_{\mathrm{i}\mathrm{c}\mathrm{e}}h\left(T_{\mathrm{o}\mathrm{i}}-T_{\mathrm{a}\mathrm{i}\mathrm{r}}\right)}{V}_{\mathrm{m}}=\mathrm{0.5}\pi {\left(r_{\mathrm{c}}+\frac{d_{\mathrm{u}}}{\mathrm{2}}\right)}^{\mathrm{3}/\mathrm{2}}\left({r_{\mathrm{c}}}^{\mathrm{1}/\mathrm{2}}+{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{1}/\mathrm{2}}\right)-\pi {r_{\mathrm{c}}}^{\mathrm{2}}\\ V_i=\pi {r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}-\pi r_{\mathrm{c}}^{\mathrm{2}}\end{array}$$

2.4. Dynamic Detachments of Remaining Ice

When the initial ice detaches from the conductor, the remaining ice hangs from the conductor. Because the remaining ice melts simultaneously along with the initially detached ice, the remaining ice will move down with the cross section shown in Figure 2, in which O₁ is the center of the circle ice layer, O₂ and O₃ are the center of the conductor before and after ice melting, separately, $d_{\mathrm{m}}$ is the rising distance of the conductor after ice melting, and the shadowed area is the air gap. Due to the flowing melted water, the air gap of the ice near the middle point is filled with water, which increases the entire mass of the remaining ice. The air gap of the remaining ice can be expressed as $$\begin{array}{}\text{(9)}& V_{\mathrm{m}\mathrm{e}\mathrm{l}\mathrm{t}}=\frac{\left[R_{\mathrm{0}}{I}^{\mathrm{2}}-\mathrm{2}\pi r_{\mathrm{i}\mathrm{c}\mathrm{e}}h\left(T_{\mathrm{o}\mathrm{i}}-T_{\mathrm{a}\mathrm{i}\mathrm{r}}\right)\right]t_{\mathrm{d}}-{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{V}_{\mathrm{i}\mathrm{c}\mathrm{e}}{C}_{\mathrm{i}\mathrm{c}\mathrm{e}}\left(T_{\mathrm{o}\mathrm{i}}/\mathrm{2}-T_{\mathrm{a}\mathrm{i}\mathrm{r}}\right)}{{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{L}_{\mathrm{f}}}\end{array}$$ The mass of the remaining ice is $$\begin{array}{}\text{(10)}& m_{\mathrm{r}}={\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}\pi \left[{\left(r_{\mathrm{c}}+d_{\mathrm{i}\mathrm{c}\mathrm{e}}\right)}^{\mathrm{2}}-r_{\mathrm{c}}^{\mathrm{2}}\right]-{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}{V}_{\mathrm{m}\mathrm{e}\mathrm{l}\mathrm{t}}\end{array}$$ The rising distance of the conductor after ice melting [10, 15] is $$\begin{array}{}\text{(11)}& d_{\mathrm{m}}=\mathrm{2}{\left[\frac{{\left(r_{\mathrm{c}}^{\mathrm{2}}+V_{\mathrm{m}\mathrm{e}\mathrm{l}\mathrm{t}}/\pi \right)}^{\mathrm{2}}}{r_{\mathrm{c}}}\right]}^{\mathrm{1}/\mathrm{3}}-\mathrm{2}{r}_{\mathrm{c}}\end{array}$$ The initial ice detachment from the conductor creates conductor vibration that may break the remaining ice. According to the detachment criterion of solid ice under instantaneous ice shedding in [15], the ice breaks when the inertia force of the ice is greater than the composition of the adhesive force on the ice and the cohesive force in the ice, and the ice is assumed to be split into two parts along line A-B-C-D as shown in Figure 2.

[figures omitted; refer to PDF]

The detachment of solid ice in [15] differs in two ways from the detachment of the remaining ice in this study. First, the adhesive force between the ice and conductor is small and can be ignored in this study [19]. Second, the mass of the remaining ice changes with the flow of melted water. Therefore, when the resultant force of the inertia force and gravity on the ice layer is greater than the cohesive force in the ice, the remaining ice detaches from the conductor. A force analysis of the remaining ice is shown in Figure 3. The detachment criterion can be written as $$\begin{array}{}\text{(12)}& F_{\mathrm{v}\mathrm{i}}+G_{\mathrm{r}}\ge F_{\mathrm{c}\mathrm{o}}\end{array}$$ where $F_{\mathrm{v}\mathrm{i}}$ is the inertia force, $G_{\mathrm{r}}$ is gravity, and $F_{\mathrm{a}\mathrm{d}}$ is the cohesive force of the ice.

[figure omitted; refer to PDF]

The forces in (12) can be expressed as $$\begin{array}{c}\text{(13)}& F_{\mathrm{v}\mathrm{i}}=m_{\mathrm{r}\mathrm{l}}{a}_{\mathrm{r}}{G}_{\mathrm{r}}=m_{\mathrm{r}\mathrm{l}}g\\ F_{\mathrm{c}\mathrm{o}}=\mathrm{2}\left(\sqrt{{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}-{d_{\mathrm{m}}}^{\mathrm{2}}}-r_{\mathrm{c}}\right){\tau }_{\mathrm{c}\mathrm{o}}\end{array}$$ where $m_{\mathrm{r}\mathrm{l}}$ is the mass of the ice under line A-B-C-D; $a_{\mathrm{r}}$ is the vertical acceleration of the remaining ice; g is the acceleration of gravity; and ${\tau }_{\mathrm{c}\mathrm{o}}$ is the cohesive strength of the ice [15].

Due to the flowing of melted water, there is no water in the air gap of the remaining ice near the suspension point, and the air gap of the remaining ice near the middle point is filled with water. With melted water in the air gap, the mass of the remaining ice under line A-B-C-D can be expressed as $$\begin{array}{}\text{(14)}& m_{\mathrm{r}\mathrm{l}\_\mathrm{1}}={\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}\left[\pi {r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}-{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}\arccos \frac{\mathrm{2}{d}_{\mathrm{m}}}{r_{\mathrm{c}}+d_{\mathrm{i}\mathrm{c}\mathrm{e}}}+d_{\mathrm{m}}\sqrt{\frac{{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}}{\mathrm{4}}-d_{\mathrm{m}}^{\mathrm{2}}}-\frac{\pi {r_{\mathrm{c}}}^{\mathrm{2}}}{\mathrm{2}}+\frac{\left({\rho }_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}-{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}\right)V_{\mathrm{m}\mathrm{e}\mathrm{l}\mathrm{t}}}{{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}}\right]\end{array}$$

From formulas (12)-(14), the critical detachment acceleration of the remaining ice with melted water can be calculated as follows: $$\begin{array}{}\text{(15)}& a_{\mathrm{c}\_\mathrm{1}}=\frac{\mathrm{2}\left(\sqrt{{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}-\mathrm{4}{d}_{\mathrm{m}}^{\mathrm{2}}}-r_{\mathrm{c}}\right)\ast {\tau }_{\mathrm{c}}}{m_{\mathrm{r}\mathrm{l}\_\mathrm{1}}}-\mathrm{g}\end{array}$$

Without melted water in the air gap, the mass of the remaining ice under line A-B-C-D and the critical detachment acceleration of the remaining ice can be expressed as $$\begin{array}{}\text{(16)}& m_{\mathrm{r}\mathrm{l}\_\mathrm{2}}={\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}\left[\pi {r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}-{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}\arccos \frac{\mathrm{2}{d}_{\mathrm{m}}}{r_{\mathrm{c}}+d_{\mathrm{i}\mathrm{c}\mathrm{e}}}+d_{\mathrm{m}}\sqrt{\frac{{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}}{\mathrm{4}}-d_{\mathrm{m}}^{\mathrm{2}}}-\frac{\pi {r_{\mathrm{c}}}^{\mathrm{2}}}{\mathrm{2}}-V_{\mathrm{m}\mathrm{e}\mathrm{l}\mathrm{t}}\right]\text{(17)}& a_{c\_\mathrm{2}}=\frac{\mathrm{2}\left(\sqrt{{r_{\mathrm{i}\mathrm{c}\mathrm{e}}}^{\mathrm{2}}-\mathrm{4}{d}_{\mathrm{m}}^{\mathrm{2}}}-r_{\mathrm{c}}\right)\ast {\tau }_{\mathrm{c}}}{m_{\mathrm{r}\mathrm{l}\_\mathrm{2}}}-\mathrm{g}\end{array}$$

3. Dynamic PIS Model Based on Finite Element Mathematical Model

Since the volume of the air gap, the mass of the remaining ice, and the cohesive force of the ice all vary with conductor vibration, it is challenging to use commercial software to model the dynamic process of PIS during thermal de-icing. In this section, a finite element mathematical model is built to model the dynamic process of PIS.

3.1. Static Model of an Iced Transmission Line

Building a finite element mathematical model of a complex scenario is a common method by which to characterize the dynamic processes of complex scenarios in mechanics research. First, the conductor can be uniformly divided into multiple segments. Then, the mathematical model of each segment, including static and dynamic equilibrium equations, is built using element models. The equations of every segment can be combined into the conductor model, which can be solved to obtain the dynamic response of the conductor.

For a suspended cable structure, the most common element models are the two-node isoparametric truss element (TNITE) [20], two-node parabolic element (TNPE) [21], two-node catenary element (TNCE) [22], and multinode isoparametric curve element (MNICE) [23]. The TNPE model and TNCE model are suitable for static analysis but are inaccurate for dynamic analysis of conductor vibration because of the deformation in the vibration. The MNICE models require more calculation resources. The TNITE model has similar accuracy but requires less calculation resources, which is more suitable for the dynamic analysis of PIS.

According to Hooke’s law, the relation between stress and strain in the element is [22] $$\begin{array}{}\text{(18)}& \left\{\sigma \right\}=E\left\{\epsilon \right\}+\left\{{\sigma }_{\mathrm{0}}\right\}=\frac{\mathrm{4}E}{L^{\mathrm{2}}}\left({\left\{X_e\right\}}^T\left[C\right]\left\{u_e\right\}+\frac{\mathrm{1}}{\mathrm{2}}{\left\{\delta \right\}}^T\left[C\right]\left\{u_e\right\}\right)+{\sigma }_{\mathrm{0}}\end{array}$$ where σ , ε , and σ ₀ are the axis stress, axis strain, and initial axis stress of the element, respectively; the axial direction means the direction perpendicular to the cross section of the conductor; E is the elastic modulus of the cable; X _e is the integral coordinate matrix of two nodes in an element before cable deformation; and H is the correlation matrix of the shape function, which can be expressed as $$\begin{array}{}\text{(19)}& \left[H\right]=\frac{d{\left[N\right]}^T}{d\xi }\frac{d\left[N\right]}{d\xi }\end{array}$$ where N is a shape function matrix defined as $$\begin{array}{}\text{(20)}& N=\left[\begin{array}{cccccc}{N}_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& N_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& N_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& N_{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& N_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& N_{\mathrm{2}}\end{array}\right]\end{array}$$ where $N_{\mathrm{1}}=\mathrm{1}/\mathrm{2}-\left(\mathrm{1}/\mathrm{2}\right)\xi$ , $N_{\mathrm{2}}=\mathrm{1}/\mathrm{2}+\left(\mathrm{1}/\mathrm{2}\right)\xi$ , and $\xi =\mathrm{2}\mathrm{s}/L$ ; ξ is a relative coordinate; s is the length from a point to the midpoint of the element; and L is the element length of the power line.

According to the virtual work principle ( δU– δW = 0), the virtual work by external forces equals that by internal forces, which can be written as $$\begin{array}{}\text{(21)}& \delta {\left\{\epsilon \right\}}^T\left\{\sigma \right\}dV-\delta {\left\{u_e\right\}}^T\left\{R_e\right\}=\mathrm{0}\end{array}$$ where $\{R_e\}$ is the matrix of the equivalent nodal load, which can be defined as $$\begin{array}{}\text{(22)}& \left\{R_e\right\}=\frac{L}{\mathrm{2}}{{\displaystyle \int }}_{-\mathrm{1}}^{\mathrm{1}}{\left[N\right]}^T\left\{q\left(s\right)\right\}d\xi \end{array}$$ where q(s) is the uniformly distributed load on the element, which is the integration of the self-weight load and the ice load before ice shedding and of only the self-weight load after ice shedding. The equivalent form of formula (22) is R _e = Kδ , which is the static equilibrium equation of the element.

By integrating formulas (18)-(22), the element equilibrium equation can be written as $$\begin{array}{}\text{(23)}& \frac{\mathrm{2}A}{L}{{\displaystyle \int }}_{-\mathrm{1}}^{\mathrm{1}}{\left[C\right]}^T\left(\left\{X_e\right\}+\left\{\delta \right\}\right)\left(\frac{\mathrm{4}E}{L^{\mathrm{2}}}\left({\left\{X_e\right\}}^T\left[C\right]\left\{\delta \right\}\right.\right.+\left.\left.\frac{\mathrm{1}}{\mathrm{2}}{\left\{\delta \right\}}^T\left[C\right]\left\{\delta \right\}\right)+{\sigma }_{\mathrm{0}}\right)d\xi -\left\{R_e\right\}=\mathrm{0}\end{array}$$

3.2. Dynamic Model of PIS

In the vibration process of the conductor, the equivalent nodal load R _e consists of the self-weight load and the ice load, in addition to the inertial force and damping force. The dynamic equilibrium equation can be expressed as $$\begin{array}{c}\text{(24)}& R_e=M\stackrel{••}{\delta }+C\stackrel{•}{\delta }+K\delta M=\frac{L}{\mathrm{2}}\ast {{\displaystyle \int }}_{-\mathrm{1}}^{\mathrm{1}}\rho {\left[N\right]}^{T}Nd\xi \\ C=\frac{L}{\mathrm{2}}\ast {{\displaystyle \int }}_{-\mathrm{1}}^{\mathrm{1}}u{\left[N\right]}^{T}Nd\xi \end{array}$$ where M, C, and K are the mass, damping, and stiffness matrices, respectively, and $\stackrel{••}{\mathit{\delta }}$ , $\stackrel{•}{\mathit{\delta }}$ , and $\mathit{\delta }$ are the nodal acceleration, velocity, and displacement, respectively.

The damping matrix C is defined in terms of Rayleigh damping: $$\begin{array}{}\text{(25)}& C=\alpha M+\beta K\text{(26)}& \alpha =\frac{\mathrm{2}\left({\xi }_i{\omega }_j-{\xi }_j{\omega }_i\right)}{{\omega }_j^{\mathrm{2}}-{\omega }_i^{\mathrm{2}}}{\omega }_i{\omega }_j\text{(27)}& \beta =\frac{\mathrm{2}\left({\xi }_j{\omega }_j-{\xi }_i{\omega }_i\right)}{{\omega }_j^{\mathrm{2}}-{\omega }_i^{\mathrm{2}}}\end{array}$$ where α and β are the mass and stiffness ratio coefficients, ω _i and ω _j are the cut-off frequencies of the lower and upper bounds in the frequency domain of interest, and ξ _i and ξ _j are the corresponding damping ratios.

4. Dynamic Model Solution

In the PIS during thermal de-icing, the flow of melted water and the detachment of ice affect the dynamic response of PIS and result in perturbation of the mass matrix, which must be corrected with the conductor vibration. In this section, the correction method of the mass matrix is presented, and the solution method for the dynamic equilibrium equations is described.

4.1. Mass Correction for the Analysis of PIS

From formulas (12)~(13), we can observe that the mass of the ice can impact the force in the ice and further impact the time of ice detachment and the accuracy of the dynamic response solution. Therefore, it is necessary to calculate the mass mutation in every element during solving of the dynamic model of PIS. The mechanism of mass mutation for an element is different before and after the initial ice detachment. Before the initial ice detachment, the iced conductor is in a static state, and the mass mutation is caused by the melted water, which flows along the air gap to the bottom of the ice layer. Nevertheless, after the initial ice detachment, the mass mutation results from the PIS and the melted water that spills from the element when the original low point of the two end nodes is higher than the original high point. Therefore, the following two mass mutation methods are discussed.

4.1.1. Mass Correction before the Initial Ice Shedding

The melted ice results in an air gap in the ice layer, as shown in Figure 1. For one element, the gap volume, which is the volume of melted ice, is greater than the volume of melted water because the density of ice is less than the density of water. The melted water flows from the element near the suspension point to the element near the lowest point.

The number of elements with water can be obtained from the volume of the melted water. Taking a conductor with an equal altitude at two end nodes for example, the volume of melted ice in the ice layer can be expressed as $$\begin{array}{}\text{(28)}& V_{\mathrm{i}\mathrm{c}\mathrm{e}}={\displaystyle \sum }_{j=\mathrm{1}}^k{V}_{\mathrm{m}\_j}{L}_{\mathrm{m}\_j}\end{array}$$ where $V_{\mathrm{m}\_j}$ refers to the volume of melted ice of element j; $L_{\mathrm{m}\_j}$ refers to the length of element j; and k refers to the number of elements.

The volume of the melted water in the ice layer can be expressed as $$\begin{array}{}\text{(29)}& V_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}=\frac{V_{\mathrm{i}\mathrm{c}\mathrm{e}}{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}}{{\rho }_{\mathrm{w}}}\end{array}$$ where ${\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}$ and ${\rho }_{\mathrm{w}}$ refer to the densities of ice and water, respectively.

The elements without water can be written as follows: $$\begin{array}{}\text{(30)}& {\displaystyle \sum }_{j=\mathrm{1}}^n{V}_{\mathrm{m}\_j}{L}_{\mathrm{m}\_j}\ge \frac{\left(V_{\mathrm{i}\mathrm{c}\mathrm{e}}-V_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}\right)}{\mathrm{2}}\end{array}$$ where n refers to the number of elements without water counted from the suspension point.

After the number of elements without water is obtained, the mass of the elements without water can be calculated using formula (14), and the mass of the other elements can be computed using formula (15).

4.1.2. Mass Correction after Initial Ice Shedding

After the initial ice shedding, the element mass can mutate with conductor vibration and the melted flowing water. For the dropped ice, the element mass shifts to the self-weight by modifying the parameters in formula (24). While the ice does not fall, the element mass is determined by the flow of melted water. When the altitude of the original low point $P_{\mathrm{m}\mathrm{1}}$ at the two end nodes is greater than the altitude of the original high point $P_{\mathrm{m}\mathrm{2}}$ , the element mass can be expressed as formula (14); otherwise, the mass can be expressed as formula (16). Therefore, the element mass with conductor vibration can be obtained as follows: $$\begin{array}{}\text{(31)}& m_{\mathrm{r}\mathrm{l}}=\left\{\begin{array}{cc}{m}_{\mathrm{r}\mathrm{l}\_\mathrm{1}}\hfill & if\hspace{0.17em}\hspace{0.17em}{P}_{\mathrm{m}\mathrm{1}}\ge P_{\mathrm{m}\mathrm{2}}\hfill \\ m_{\mathrm{r}\mathrm{l}\_\mathrm{2}}\hfill & if\hspace{0.17em}\hspace{0.17em}{P}_{\mathrm{m}\mathrm{1}}<P_{\mathrm{m}\mathrm{2}}\hfill \end{array}\right.\end{array}$$

The altitude of the two end nodes must be updated in each iteration.

4.2. Solution of the Dynamic Equilibrium Equations

The dynamic equilibrium equations for analyzing the PIS are higher-order differential-algebra equations, which makes obtaining an analytic solution challenging. The method integrating the Newmark-β and Newton-Raphson methods is applied to the dynamic equilibrium equations in this study.

The dynamic equilibrium equation (24) at time $t+\Delta t$ can be written as $$\begin{array}{}\text{(32)}& M{\stackrel{••}{\delta }}_{t+\Delta t}+C{\stackrel{•}{\delta }}_{t+\Delta t}+K{\delta }_{t+\Delta t}=R\left(e\right)\end{array}$$

The solution procedure for PIS induced by thermal de-icing is shown in Figure 4 and consists of the following steps:

(1)

Initialize. Set the transmission line parameters (such as the span length and conductor type), icing parameters (such as the shape, density, and thickness), and ice melting parameters (such as the ice melting current and ambient temperature).

(2)

Calculate the melting time of the initial ice shedding from the ice detachment model of initial ice shedding in Section 2 and the mass of melted water by the mass correction before initial ice shedding in Section 4.

(3)

Employ the above parameters to solve the static equilibrium equations before the initial ice shedding.

(4)

Calculate the critical vertical acceleration of the remaining ice on each element according to the ice detachment model of the remaining ice in Section 2.

(5)

Set the initial displacement ${\delta }_{t+\Delta t}$ equal to zero.

(6)

Modify the mass matrices in Section 4.

(7)

Update the matrix values corresponding to ${\delta }_{t+\Delta t}$ , including the value of the equivalent nodal load matrix, internal force matrix, tangent stiffness matrix, and strain matrix at $t+\Delta t$ .

(8)

Modify the tension and tangent stiffness matrix while ε ≤ 0.

(9)

The values of the displacement, velocity, and accelerated speed at t clock and the values of the displacement, internal force matrix, and mass matrix at $t+\Delta t$ are substituted into (32) to calculate the value of deviation $\Psi$ . By substituting $\Psi$ and the tangent stiffness matrix at time $t+\Delta t$ into formula (33), the new displacement ${\delta }_{t+\Delta t}$ can be obtained. $$\begin{array}{}\text{(33)}& {\delta }_{t+\Delta t}={\delta }_t-{\left({K_T}^{t+\Delta t}+\frac{\mathrm{1}}{\alpha \Delta t^{\mathrm{2}}}M+\frac{\delta }{\alpha \Delta t}C\right)}^{-\mathrm{1}}\psi \end{array}$$

(10)

Go back to step (9) if the difference in internal force between two continuous iterations is greater than the liminal value. Otherwise, continue.

(11)

Calculate the values of velocity and acceleration. If the stop time is not reached, go back to step (6) to compute the parameters in the next time step. Otherwise, end the calculation.

[figure omitted; refer to PDF]

5. Validation of the Proposed Method

In this section, the results of the proposed method are compared with the results of a small-scale experiment to confirm the effectiveness of the TITNE method. Then, the results of a full-scale experiment are employed to verify the effectiveness of the proposed ice detachment model by comparison with the results of PIS with constant speed.

5.1. Effectiveness of the Proposed TNITE Model

5.1.1. Experimental Configuration

Because the uniform accreted ice along the conductor is challenging to simulate in the warm regions and warm seasons, it is common to simulate the icing on conductors as lumped loads [8, 24]. The experimental results of a single span in [8] are employed to validate the effectiveness of the proposed TNITE model.

The accreted ice in the experiment is simulated as 10 lumped loads that are fixed along the 235 m-long conductor with remote-controlled cutters. The cutters can release the lumped loads partially or simultaneously. In this paper, the scenario in which the lumped loads are released simultaneously is employed to verify the effectiveness of proposed method. The parameters of the conductor (LGJ 630/45) are shown in Table 1. The equivalent ice thickness simulated by the loads is 15 mm. The lumped loads are released simultaneously in the experiment.

Table 1

Mechanical parameters of LGJ 630/45.

Parameter	Unit	LGJ 630/45
Wind span	m	235
Cross-sectional area	mm2	666.55
Young’s modulus	MPa	63,000
Weight per unit length	kg/m	2.06
Diameter	mm	33.6
Rated tensile strength	N	148,700
Initial horizontal stress	N/mm2	39.6069

To obtain the modal damping ratios, a spectral analysis of experimental data is conducted, and the results are shown in Figure 5. Two significant frequencies, ω ₁ = 0.398 Hz and ω ₂ = 0.596 Hz, are observed. Meanwhile, the damping ratios calculated by Half power bandwidth method are about ξ ₁ = 0.1055 and ξ ₂ = 0.0713. Consequently, the mass and stiffness ratio coefficients are calculated by (25)~(27) as α = 0.0832 and $\beta$ =0.0051.

[figure omitted; refer to PDF]

5.1.2. Results Comparison

The vertical displacement of the experiment, ANSYS, and the proposed method are shown in Figure 6. In the simulation using ANSYS, Link10 is used to simulate the transmission conductor. The transmission conductor is divided into 100 elements. The simulation time is 10 s, and the simulation step size is 0.05 s for both ANSYS and the method developed in this study. The amplitude and fluctuation trends of the displacement are similar for the three methods, and the amplitude decreases under the action of damping with increasing vibration time. The amplitude and fluctuation trends of tension are very similar for the three methods. Some displacement errors appear between the proposed method and the experiment, possibly due to the damping coefficients in the simulation not accurately representing the experimental value and the fact that certain devices used in the experiment, such as the tension sensor and insulators, were not modeled in the simulation. However, we can conclude that the conductor model TNITE has a similar precision to that of the experiment and commercial software and has sufficient precision to simulate PIS by thermal de-icing.

[figures omitted; refer to PDF]

Because the tension curve is very smooth, we need to confirm whether the FE model properly captures the shock wave induced in the conductor by the weight loads dropping. The effectiveness is confirmed in terms of mesh size and time steps. The vertical displacement and tension results of the proposed method with different mesh sizes and time steps are illustrated in Figure 7. When the mesh sizes are 200 and 300 and the time steps are 0.01 s and 0.005 s, the vertical displacement and tension result are almost the same as the result with 100 elements and a 0.05 s time step. This result means that the model in this paper can capture the shock wave properly.

[figures omitted; refer to PDF]

5.2. Effectiveness of the Proposed PIS Model

PIS is the cracking process of ice on a conductor. Its cracking speed is difficult to determine. Therefore, PIS is a discrete dynamic process of some small segments. When the segments are small enough, the dynamic response of the small segments is equivalent to the PIS. In this subsection, this idea is verified by comparing the result of proposed method with the result of thermal de-icing, as well as by comparison with other PIS models.

5.2.1. Experimental Configuration

A full-scale experiment of thermal de-icing was conducted in the Xuefeng Mountain Natural Icing Station established by Chongqing University [18]. There are two towers in the station. The distance between the towers is 80 m, between which various types of conductors are installed [25]. The transmission conductors, type LGJ-300, are powered by DC current for thermal de-icing. The parameters of the ice shedding scenario are illustrated in Table 2. The PIS during thermal ice shedding is recorded and described in [12], which is employed to validate the proposed method.

Table 2

Mechanical parameters of the full scale experiment.

Parameter	Unit	Value
$${\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}$$	kg/m3	670
$$T_{\mathrm{o}\mathrm{i}}$$	°C	-1.093
$$T_{\mathrm{a}\mathrm{i}\mathrm{r}}$$	°C	-4
R 0	Ω	0.09614×10−3
I	A	600
h	W/(m2·K)	4.7835
$$r_{\mathrm{i}\mathrm{c}\mathrm{e}}$$	M	12.1×10−3
d u	m	19.26×10−3
$$r_{\mathrm{c}}$$	m	36.2×10−3

5.2.2. Simulation Configuration

The experiment in [12] was simulated using uniform loads fixed along the 80 m-long conductor. The shape of the ice coating and the conductor in the experiment is equivalent to the eccentric circle illustrated in Figure 8. The eccentricity of the ice coating gradually increases from the middle to the suspension point, which can be measured by the distance between the center of the iced conductor and the center of the cable before thermal de-icing. The ice coating diameter is 24.1 mm, and the distance between the center of the conductor and the center of the ice coating is 12.0 mm and 7.2 cm for ices 1 and 2, respectively. The influence of the insulator string and vibration coupling effect of the transmission line and tower is ignored in this isolated-span power line model. The damping coefficients cannot be obtained by the method in section A, because the ice is shed with the fluctuation of the conductor. The value selection of damping coefficients needs to depend on experience. The damping of cable is modeled as equivalent viscous damping based on a lumped parameter model in many papers [4–6]. In this paper, the method used to select damping refers to the method in [8], which compares the results of numerical computations and physical tests. A series of damping parameters is set, and the damping coefficients which have the best fit for the results of numerical computation and physical tests are adopted to be the damping coefficients for the numerical computations. The simulation results of different mass ratio coefficients are shown in Figure 9. When the mass ratio coefficient is 0.17 and the stiffness ratio coefficient is 0, the simulation result best fits the experimental result. The temperature is -5°C, and the ice melting current is 1000 A.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

5.2.3. Results Comparison

Based on the experimental results, the results of the proposed method are compared with the results of PIS with a constant ice shedding speed. The vertical displacement of the middle point was illustrated in [18] with details and is employed as a benchmark to compare the different methods. The results of the proposed method and the PIS with a constant ice shedding speed are demonstrated in Figure 10.

[figure omitted; refer to PDF]

(1) Accuracy of the Proposed Method. Before PIS, the conductor is in a static state, and the sag of the iced conductor is 2.39 m. When the initial ice sheds from the conductor, the middle point jumps with the fluctuation of the conductor. The jump height of the middle point in the first peak is 4.05 m as determined by the proposed method, while it is 4.06 m in the experiment. The time to reach the peak is 1.8 s for the proposed method and the experiment. The vertical displacement and peak times of the proposed method and the experiment are similar. The jump height of the middle point in the first valley is 0.93 m in the experiment and 0.75 m as determined by the proposed method. The time to reach the valley is 2.55 s and 2.8 s for the experiment and the proposed method, respectively. The peak time, the maximum vertical displacement, and the vibration period of the experiment and proposed method are similar. The vertical displacements of the proposed method and the experiment contain some error because the ice was not completely shed in the experiment. Therefore, the proposed method is accurate for simulating PIS.

(2) Accuracy of PIS with Constant Speed. PIS with constant speed is simulated using the same simulation scenario as in the proposed method. For PIS with a constant ice shedding speed, the dynamic response is related to the ice shedding speed. When the ice shedding speed is 60 m/s, the vertical displacement of the middle point in the first peak is 3.52 m, the time to reach the peak is 1.5 s, the position of the middle point in the valley is -1.52 m, and the time to reach the valley is 2.6 s. When the ice shedding speed is 84.7 m/s, the vertical displacement of the middle point in the first peak is 3.84 m, the time to reach the peak is 1.4 s, the position of the middle point in the valley is -1.74 m, and the time to reach the valley is 2.4 s. When the ice shedding speed is 100 m/s, the vertical displacement of the middle point in the first peak is 4.59 m, the time to reach the peak is 1.4 s, the position of the middle point in the valley is -1.60 m, and the time to reach the valley is 2.4 s. All the simulation results for various ice shedding speeds have some errors relative to the experimental results. The times to reach the first peak for all ice shedding speeds are less than those in the experiment. For the PIS with a speed of 60 m/s, the vibration period is similar to that in the experiment, but the vertical displacement in the peak and valley is less than that in the experiment. For the PIS with a speed of 84.7 m/s, the vertical displacement in the first peak is similar to that in the experiment, but the position of the valley and the time to reach the peak and valley are less than those in the experiment. For the PIS with a speed of 100 m/s, the vertical displacement in the peak is larger than that in the experiment, while the position of the valley and the time to reach the peak and valley are less than those in the experiment. Therefore, PIS with a constant ice shedding speed is less accurate than the proposed method.

In [8], the PIS speed is set to be less than 50 m/s. The results of this study are inaccurate for analyzing PIS from the results of constant speed in this study.

Ice shedding can result in a transversal wave along the span. The transversal wave was assumed to generate ice shedding in [15]. The theoretical wave speed can be obtained using the following equation: $$\begin{array}{}\text{(34)}& v_{\mathrm{w}\mathrm{a}\mathrm{v}\mathrm{e}}=\sqrt{\frac{T}{m}}\end{array}$$ where T is the conductor tension and m is the mass per unit length of the string.

The transverse wave speed of the experiment is 84.7 m/s according to the theoretical model (34). The speed is consistent with the observation in the experiment [12]. However, the dynamic response of PIS with a speed of 84.7 m/s is inconsistent with the results of the experiment illustrated in Figure 10, indicating that the simulation results of ice shedding with a constant ice shedding speed are inaccurate, even though the ice shedding speed is accurate.

(3) Reason for Inaccuracies in PIS with a Constant Speed. The mechanism of the dynamic response difference with different shedding speeds is that the ice shedding speed can impact the mass and vertical speed of the segments. The theoretical transverse wave speed caused by initial ice shedding can be employed as a reference. If the ice shedding speed is faster than the transverse wave speed, the ice is detached from the segment when the wave arrives at the segment. On one hand, the mass of the segment is smaller than that of the iced segment, which results in a larger vertical displacement of the conductor by the energy of the transverse wave. On the other hand, the detachment of ice will generate vertical movement of the conductor. The vertical movement with the addition of the theoretical wave will result in a higher vertical displacement. If the ice shedding speed is slower than the transverse wave speed, the energy of the transverse wave needs to afford the movement of the iced conductor, which will result in a smaller vertical displacement of the segment. Because the actual ice shedding time is different from the theoretical transverse wave speed, the dynamic response of PIS with the theoretical transverse wave speed is inconsistent with the results of the experiment.

The difference between the PIS with theoretical transverse wave and the proposed method is the IDC. The IDC of the proposed method is that the inertia force of ice is greater than the composition of the adhesive force on the ice and the cohesive force in the ice, while the IDC of PIS with theoretical transverse wave is that the detachment time and the time the wave arrives are according to the speed. The ice shedding time of the proposed method is longer than the arriving time of the transversal wave because the inertia force for the ice shedding is calculated from the vertical acceleration. The vertical acceleration of the middle point is shown in Figure 11. The time at which the transverse wave arrives is 0.47 s, and the acceleration begins to increase with the wave propagation. The time of maximum acceleration is approximately 0.6 s, which is approximately 0.13 s later than when the transverse wave arrives. The delay time of the ice detachment from the conductor is the main reason for the inaccuracy of PIS with constant speed.

[figure omitted; refer to PDF]

6. Conclusion

To study the dynamic response of PIS, a mechanical model of PIS during thermal de-icing is established in this study. An ice detachment model for the initial ice and the remaining ice is built based on an ice melting model and an ice mechanical model. Then, a finite-element mathematical model is proposed to simulate the dynamic process of PIS. Finally, the proposed model is verified with experimental results. Comparison with a small-scale experiment shows that the dynamic simulation result is as accurate as that of the commercial finite element software ANSYS. Based on the results of a full-scale experiment, a comparison of the proposed method and PIS with a constant ice shedding speed shows that the proposed method is more accurate for simulating PIS. Because of the transparent modeling and solving process, the proposed method can be used to analyze complex ice shedding scenarios and can serve as a reference tool for system operators to optimize transmission line designs and evaluate the application of thermal de-icing devices.

Because the TNITE model cannot bend, mechanical ice fracturing cannot be considered in this model. Neglecting mechanical ice fracturing will result in errors in the ice detachment time. This may be one of the reasons that the vertical displacement of the proposed method has some error relative to the experiment. In future work, an element model which can consider the bend and compression should be employed to simulate the dynamic process of PIS.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (51507080, 61673213) and the Fundamental Research Funds for the Central Universities (30918011330).

References

[1] X. Xu, D. Niu, P. Wang, Y. Lu, H. Xia, "The weighted support vector machine based on hybrid swarm intelligence optimization for icing prediction of transmission line," Mathematical Problems in Engineering, vol. 2015, 2015.

[2] L. E. Kollár, M. Farzaneh, "Vibration of bundled conductors following ice shedding," IEEE Transactions on Power Delivery, vol. 23 no. 2, pp. 1097-1104, DOI: 10.1109/TPWRD.2007.915876, 2008.

[3] L. E. Kollar, M. Farzaneh, "Modeling sudden ice shedding from conductor bundles," IEEE Transactions on Power Delivery, vol. 28 no. 2, pp. 604-611, DOI: 10.1109/TPWRD.2012.2227281, 2013.

[4] M. R. Fekr, G. McClure, "Numerical modelling of the dynamic response of ice-shedding on electrical transmission lines," Atmospheric Research, vol. 46 no. 1-2,DOI: 10.1016/s0169-8095(97)00046-x, 1998.

[5] T. Kálmán, M. Farzaneh, G. McClure, "Numerical analysis of the dynamic effects of shock-load-induced ice shedding on overhead ground wires," Computers & Structures, vol. 85 no. 7, pp. 375-384, DOI: 10.1016/j.compstruc.2006.11.026, 2007.

[6] G. McClure, M. Lapointe, "Modeling the structural dynamic response of overhead transmission lines," Computers & Structures, vol. 81 no. 8, pp. 825-834, DOI: 10.1016/S0045-7949(02)00472-8, 2003.

[7] L. E. Kollár, M. Farzaneh, "Modeling the dynamic effects of ice shedding on spacer dampers," Cold Regions Science and Technology, vol. 57 no. 2, pp. 91-98, DOI: 10.1016/j.coldregions.2009.03.004, 2009.

[8] X. Meng, L. Wang, L. Hou, G. Fu, B. Sun, M. MacAlpine, W. Hu, Y. Chen, "Dynamic characteristic of ice-shedding on UHV overhead transmission lines," Cold Regions Science and Technology, vol. 66 no. 1, pp. 44-52, DOI: 10.1016/j.coldregions.2011.01.005, 2011.

[9] F. Yang, J. Yang, Z. Zhang, "Unbalanced tension analysis for UHV transmission towers in heavy icing areas," Cold Regions Science and Technology, vol. 70, pp. 132-140, DOI: 10.1016/j.coldregions.2011.07.008, 2012.

[10] K. Ji, X. Rui, L. Li, A. Leblond, G. McClure, "A novel ice-shedding model for overhead power line conductors with the consideration of adhesive/cohesive forces," Computers & Structures, vol. 157, pp. 153-164, DOI: 10.1016/j.compstruc.2015.05.014, 2015.

[11] X. L. Jiang, S. H. Fan, "Simulation and experimental investigation of dc ice-melting process on an iced conductor," IEEE Transactions on Power Delivery, vol. 25 no. 2, pp. 919-929, DOI: 10.1109/TPWRD.2009.2037632, 2010.

[12] X. Jiang, M. Bi, Z. Li, Z. Xiang, B. Dong, S. Zhao, "Study on DC ice melting and ice shedding process under natural condition," Dianwang Jishu/Power System Technology, vol. 37 no. 9, pp. 2626-2631, 2013.

[13] X. Meng, L. Hou, L. Wang, M. MacAlpine, G. Fu, B. Sun, Z. Guan, W. Hu, Y. Chen, "Oscillation of conductors following ice-shedding on UHV transmission lines," Mechanical Systems and Signal Processing, vol. 30, pp. 393-406, DOI: 10.1016/j.ymssp.2011.10.020, 2012.

[14] F. Mirshafiei, G. McClure, M. Farzaneh, "Modelling the dynamic response of iced transmission lines subjected to cable rupture and ice shedding," IEEE Transactions on Power Delivery, vol. 28 no. 2, pp. 948-954, DOI: 10.1109/TPWRD.2012.2233221, 2013.

[15] K. P. Ji, X. M. Rui, L. Li, "Dynamic response of overhead transmission lines with eccentric ice deposits following shock loads," IEEE Transactions on Power Delivery, vol. 32 no. 3, pp. 1287-1294, DOI: 10.1109/TPWRD.2015.2501029, 2017.

[16] W. Yaoxuan, J. Xingliang, F. Songhai, M. Zhigao, "Asynchronism of ice shedding from the de-iced conductor based on heat transfer," IET Science, Measurement & Technology, vol. 10 no. 4, pp. 389-395, DOI: 10.1049/iet-smt.2015.0258, 2016.

[17] H.-T. Thai, S.-E. Kim, "Nonlinear static and dynamic analysis of cable structures," Finite Elements in Analysis and Design, vol. 47 no. 3, pp. 237-246, DOI: 10.1016/j.finel.2010.10.005, 2011.

[18] X. Jiang, Y. Wang, L. Shu, Z. Zhang, Q. Hu, Q. Wang, "Control scheme of the de-icing method by the transferred current of bundled conductors and its key parameters," IET Generation, Transmission & Distribution, vol. 9 no. 15, pp. 2198-2205, DOI: 10.1049/iet-gtd.2015.0196, 2015.

[19] S. Y. Sadov, P. N. Shivakumar, D. Firsov, S. H. Lui, R. Thulasiram, "Mathematical model of ice melting on transmission lines," Journal of Mathematical Modelling and Algorithms, vol. 6 no. 2, pp. 273-286, DOI: 10.1007/s10852-006-9043-4, 2007.

[20] A. Tiar, W. Zouari, H. Kebir, R. Ayad, "A nonlinear finite element formulation for large deflection analysis of 2D composite structures," Composite Structures, vol. 153, pp. 262-270, DOI: 10.1016/j.compstruct.2016.05.102, 2016.

[21] W.-X. Ren, M.-G. Huang, W.-H. Hu, "A parabolic cable element for static analysis of cable structures," Engineering Computations (Swansea, Wales), vol. 25 no. 4, pp. 366-384, DOI: 10.1108/02644400810874967, 2008.

[22] M.-G. Yang, Z.-Q. Chen, X.-G. Hua, "A new two-node catenary cable element for the geometrically non-linear analysis of cable-supported structures," Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 224 no. 6, pp. 1173-1183, DOI: 10.1243/09544062JMES1816, 2010.

[23] Y. Wang, S. R. Zuo, C. Wu, "A finite element method with six-node isoparametric element for nonlinear analysis of cable structures," Applied Mechanics and Materials, vol. 275, pp. 1132-1135, DOI: 10.4028/www.scientific.net/AMM.275-277.1132, 2013.

[24] L. E. Kollár, M. Farzaneh, P. van Dyke, "Modeling ice shedding propagation on transmission lines with or without interphase spacers," IEEE Transactions on Power Delivery, vol. 28 no. 1, pp. 261-267, DOI: 10.1109/TPWRD.2012.2212918, 2013.

[25] S. H. Fan, X. L. Jiang, "DC ice-melting model for elliptic glaze iced conductor," IEEE Transactions on Power Delivery, vol. 26 no. 4, pp. 2697-2704, DOI: 10.1109/TPWRD.2011.2156433, 2011.