# Investigation and Development of a Three-Dimensional Transmission Tower-Line System Model Using Nonlinear Truss and Elastic Catenary Elements for Wind Loading Dynamic Simulation

Authors: Xiao Zhu; Ge Ou
Journal: Shock and Vibration
 (2021)

---

## Abstract

The accuracy of transmission tower-line system simulation is highly impacted by the transmission line model and its coupling with the tower. Owing to the high geometry nonlinearity of the transmission line and the complexity of the wind loading, such analysis is often conducted in the commercial software. In most commercial software packages, nonlinear truss element is used for cable modeling, whereas the initial strain condition of the nonlinear truss under gravity loading is not directly available. Elastic catenary element establishes an analytical formulation for cable structure under distributed loading; however, the nonlinear iteration to reach convergence can be computational expensive. To derive an optimal transmission tower-line model solution with high fidelity and computational efficiency, an open-source three-dimensional model is developed. Nonlinear truss element and elastic catenary element are considered in the model development. The results of the study imply that both elements are suitable for the transmission line model; nevertheless, the initial strain in nonlinear truss element largely impacts the model accuracy and should be calibrated from the elastic catenary model. To cross-validate the developed models on the coupled transmission tower and line, a one-span eight-line system is modeled with different elements and compared with several state-of-the-art commercial packages. The results indicate that the displacement time-history root-mean-square error (RMSE) of the open-source transmission tower-line model is less than1% and with a 66% computational time reduction compared with the ANSYS model. The application of the open-source package transmission tower-line model on extreme wind speed considering the aerodynamic damping is further implemented.


---

## Body

## 1. Introduction

Transmission tower-line systems connect power plants to customers and are widely distributed throughout the country. The failure of the transmission systems can result in tremendous economic and social life loss. In 2018, Hurricane Michael caused a widespread power outage that affected 1.7 million customers in six states in the US [1]. In the report, 116 transmission lines were damaged and led to a blackout with 16 deaths and $25 million in economic loss. Because of the high-rise of the transmission tower and the long span of the transmission line, the transmission system is sensitive to environmental wind loadings [2]. Consequently, the structural response of the transmission tower-line system under wind loading has drawn lots of researchers’ attention. The development of a high fidelity and computationally efficient model which can estimate the dynamic response of transmission tower-line system during extreme wind condition is prominent.The transmission tower-line system contains two major components: the transmission towers and transmission lines. Computational models are used to evaluate the linear and nonlinear time-history of the transmission tower systems. Most of the previous studies have been made to investigate the stand-alone transmission tower performance. Kemp and Behncke [3] and Alam and Santhakumar [4] conducted full-scale experiments to indicate that large bending moments existed in the tower legs and cross-bracing system. Xie and Sun [5] further investigated the tower failure mechanism under bending and flexure-rotation loading. The results embodied that, at bending load, tower failure was dominated by the leg buckling and at flexure-rotation loading, the diagonal bracing member buckling was the main reason. Tian et al. [6] proposed a beam finite element tower model associated with a user-defined material model in ABAQUS to simulate the behavior of the tower under various loading conditions. The full-scale test obtained tower load-bearing capacity and failure mode verified the accuracy of the proposed finite element model. Battista et al. [7] summarized that the bending stress caused by the tower member connections was important to evaluate the ultimate strength of the tower, and consequently, spatial frame element should be used to build the tower finite element model. From the literature, the failure mechanism of the transmission tower is mainly due to the buckling of the elements.The transmission line in the transmission tower-line system contains a set of wires connected to the tower, including conductors and ground wires. These wires are considered as cable structures that are highly geometrically nonlinear. Hence, flexibility and large deformation must be taken into consideration in the transmission line model. Owing to the computational complexity in modeling the dynamics of nonlinear transmission line under seismic or wind loading and the coupling between transmission tower and line, researchers started to explore the impact of transmission line dynamic to the transmission tower response. Momomura et al. [8] conducted experiments in mountain areas and concluded that the presence of the conductor significantly affected the vibration of the tower. Deng et al. [9, 10] implemented a series of wind tunnel tests to investigate the influence of the conductor on the transmission tower-line system from different wind angles. They concluded that when the wind direction was perpendicular to the transmission tower-line system, the vibration of the tower increased compared with a single tower, and the conductors enhanced the damping of the tower. Xie et al. [11] performed a wind tunnel test to study the displacement of the tower with and without conductors. They found that 70–90% of the tower displacement was induced by the conductors.Because of the low flexural stiffness of the transmission line and free of rotational degree of freedom, nonlinear truss element is conventionally used for transmission lines in the transmission system modeling [12–18]. However, the gravity loading cannot be automatically accounted for in the nonlinear truss element; accordingly, calibration of its initial strain is required to be able to well represent the internal stress of the cable as well as the accurate cable shape. McClure et al. [19, 20] proposed that, for the transmission line modeling using the nonlinear truss element, approximated initial strain should be applied at first to avoid singularities in the initial stiffness matrix formulation. The initial strain was calculated from the catenary equation. Zhang et al. [21] stated that the self-weight of the transmission line introduced initial strains in the transmission line, and to find the initial shape and initial strain, a trial-and-error auto-gravimetric analysis in ANSYS should be performed. Although Keyhan et al. [19, 20] and Zhang et al. [21] mentioned how to set the initial strain of the nonlinear truss element, very few research documents have discussed the impact of the importance of calibrated initial strain on static performance of cable and its impacts on the global performance of coupled transmission tower-line systems.Another candidate for transmission line modeling is the cable element. The elastic catenary cable element takes the self-weight of the transmission line into the internal force vector and stiffness matrix formulation directly without any approximations [22]. Although the elastic catenary equation was first obtained by Leibniz in 1691 [23], the explicit tangent stiffness matrix and the internal force vector have only recently been well studied. Jayaraman and Knudson [24] firstly developed the stiffness matrix and internal force vector of the cable structure in two-dimensional space. Thai and Kim [22] extended the elastic catenary cable elements in three-dimensional space. Salehi et al. [25] further extended the elastic catenary cable element by considering the uniformly distributed load in all directions. Experiments and numerical comparison showed that the elastic catenary element could calculate the displacement of the cable accurately and efficiently [26, 27]. Eventually, the elastic catenary cable element can be treated as a reference to verify the accuracy of the nonlinear truss element formulation under static and dynamic loading. However, the elastic catenary cable element is not supported by most of the popular commercial software packages (ANSYS and ABAQUS). SAP2000 commercial software has built-in elastic catenary cable element, but user-defined dynamic loading is not supported. Besides, the stiffness matrix of the elastic catenary element is calculated by inversing the flexibility matrix, in which the computational complexity increases rapidly with the increasing number of elements.In summary, the accuracy of transmission line structure will be compromised by using the nonlinear truss element if its initial strain is not carefully calibrated. Elastic catenary element offers analytical formulation of cable structures. However, the computational efficiency of a multiple-line structure can be jeopardized by large matrix inversion operation, specifically during a dynamic simulation. Owing to such stated limitations and drawbacks in the state-of-the-art option, we propose to identify an optimal candidate to represent transmission lines in a high-fidelity and high-computation efficient transmission tower-line model. To achieve such an objective, the performance of an individual transmission line with different elements should be investigated first. Then a transmission tower-line model open-source MATLAB software package with the transmission line utilizing elastic catenary cable element and nonlinear truss element is developed and compared. The initial strain condition of nonlinear truss element is modeled in two forms. The first one is using an uncalibrated strain which is the default of most commercial software, and the initial strain of the other nonlinear truss element is calibrated to the elastic catenary element under gravity loading. The developed models are compared and cross validated with ANSYS transient, SAP2000, and ANSYS LS-DYNA.To further investigate the accuracy and the impact of using different models on modeling the large nonlinear transmission line and the coupling between line and tower, the dynamic displacement response of a one-span eight-line transmission tower-line model under different wind speeds and wind angles is compared. In the one-span eight-line transmission system case, four transmission tower-line models are compared: (1) lumped mass transmission tower coupled with elastic catenary cable element line in MATLAB (ECE-MATLAB); (2) lumped mass transmission tower coupled with calibrated nonlinear truss element line in MATLAB (NLTc-MATLAB); (3) transmission tower coupled with calibrated nonlinear truss element line in ANSYS (NLTc-ANSYS); (4) ANSYS auto-gravimetric model (ANSYS Auto-Gravimetric), which is developed by Zhang et al. [21] and discussed in Section 3.2.The rest of the article is organized as follows: in Section2, the methodology of developing the transmission tower, transmission line, and transmission tower-coupled finite element models are described. In Section 3, the transmission line initial shape and initial strain finding algorithms are described in detail. In Section 4, the single transmission tower, single transmission line, and a one-span two-line transmission system numerical model comparisons are presented. In Section 5, the numerical model in full transmission tower-line model setting up will be implemented and compared. In Section 6, the effects of the aerodynamic damping on the transmission tower-line system are investigated. Section 7 will conclude and summarize the whole study. The open-source package will be distributed on GitHub at [28] after the article is published.


## 2. Methodology

The transmission tower-line system contains two separately designed components: the overhanging wires and the supporting structures. The wire system comprises the ground wires and the conductors. Because the functionality of the ground wire and conductor is different, the material and geometric properties are generally different too. The supporting structure consists of transmission towers and foundations. In engineering, the transmission towers are considered fixed at the foundations and provide constraints to the wire system. Consequently, the transmission tower-line structure is fixed at tower bottom position and the transmission lines hanging between adjacent towers.


### 2.1. Transmission-Line Modeling

Transmission line under gravity forms a natural catenary shape [22]. Thus, the shape of the transmission line can be found by the elastic catenary element under loading without approximating the gravity. The two-node nonlinear truss element not having the rotational degree of freedom is widely utilized for transmission line modeling. Nevertheless, the initial strain of the nonlinear truss element has a significant influence on the structural behavior of the transmission line. Therefore, the calculation of the initial strain of the nonlinear truss element needs to be addressed. The initial strain of the nonlinear truss element can be calculated as(1)ε0=FEA0,where ε0 is the initial strain of the nonlinear truss element, F is the transmission line end force, E and A are Young’s modulus and cross-sectional area, and 0 indicates the initial state. Yet, in most cases, instead of the cable end force, the sag, the initial length of the cable is provided. Hence, the method to find the initial shape and initial strain of the cable under gravity loading needs to be explored. For the static analysis, the mass of each node in the transmission line is not used, whereas in the dynamic analysis, the mass matrix ML is used to solve the nonlinear equation of motions. For the node in between two elements, the mass of that node sums up half of the two adjacent element mass, whereas for the boundary nodes, the mass takes half of the associated element.


#### 2.1.1. Nonlinear Truss Element

The tension-only two-node three-dimensional truss element is employed both in ANSYS and MATLAB models. The tangent stiffness matrix and internal force vector equations of the nonlinear truss element are derived based on Jürgen and Bathe [29]. With the load applied, the Newton–Raphson iteration is needed to iteratively calculate the current state (∗) tangent stiffness and the internal force vector, which is shown in equations (2) and (3):(2)k∗=dσ∗dλ∗A∗1l0M∗+σ∗A∗1l∗N∗,(3)fi∗=σ∗A∗V∗,where dσ∗ is the stress change, dλ∗ is the strain change, A∗ is the cross-sectional area, l0 is the initial truss length, l∗ is the deformed truss length, and M∗ and N∗ are orthogonal coordinate transformation matrices. In this paper, only the geometric nonlinearity is considered; thus, dσ∗/dλ∗=E. For two-dimensional truss element, M∗, N∗, and V∗ are given by(4a)M∗=c2cs−c2−cscss2−cs−s2−c2−csc2cs−cs−s2css2∗,(4b)N∗=s2−cs−s2cscss2−cs−s2−c2−csc2cs−cs−s2css2∗,(4c)V∗=−c,−s,c,s∗T,where c is cosθ and s is sinθ.When the difference between the internal force vector and the applied load vector is smaller than the defined convergence criteria, the iteration will stop. The tangent stiffness matrix, internal force vector, and the displacement can then be obtained. In equations (2) and (3), the nonlinear strain has different definitions in different material laws. In ANSYS, the strain for large elongation is defined as a logarithmic strain. In this paper, the strain adopts logarithmic strain because the developed open-source transmission system finite element model is compared with ANSYS models. The logarithmic strain and cross section are defined as(5a)λ=ll0,(5b)ε=lnλ,(5c)dσ∗=Eε,(5d)σ∗=dσ∗+σ0k,ifdσ∗+σ0≥0,0,ifdσ∗+σ0<0,(6)A∗=A0λ−2ν,where l is the deformed length, l0 is the original length, λ is the axial elongation factor, ε is the logarithmic strain, dσ∗ is the stress change, σ∗ is the current state stress, σ0 is the initial stress, E is the material elastic modulus, and A∗ is the cross-section area at current state. In equation (5d), the strain is positive definite because the transmission line is tension-only structure.For dynamic analysis, the equation of motion under external loading can be written as(7)MU¨+CU˙+K∗U=Pt,where M, C, and K∗ are the mass, damping, and stiffness matrices; U¨, U˙, and U are the vectors of the acceleration, velocity, and displacement; and Pt is the external loading vector.For the transmission line, the aerodynamic damping greatly influences the dynamic displacement response under large wind speed [30]. Consequently, the damping matrix in equation (7) should take the aerodynamic damping into consideration. Wang et al. [31] proposed a closed-form formulation of the aerodynamic damping ratios for the transmission lines. Therefore, to calculate the aerodynamic damping, the method proposed by Wang et al. [31] is adopted. The detailed derivation to calculate the aerodynamic damping can be found in [31]. However, in [31], instead of the explicitly calculating the damping matrix, the aerodynamic damping modal ratio is calculated. To calculate the damping matrix, the following procedures are applied:(1)
From [31], the in-plane and out-of-plane frequencies Ω and the corresponding mode shapes V are calculated and assembled the mass matrix as M(2)
From [31], the in-plane and out-of-plane direction modal aerodynamic damping ratios are calculated as ξ(3)
The aerodynamic damping matrix can be calculated asCA=MVM−12ξMΩMV′MBy considering the Rayleigh damping, the damping matrix can be represented as(8a)CR=αMM+αKK,(8b)CA=MVM−12ξMΩMV′M,(8c)C=CR+CA,where αM and αK are the mass and stiffness proportional damping coefficients, respectively.To solve the equation of motion, nonlinear Newmark-Beta method is employed [32]. In Newmark-Beta method, at each time instant i, the displacement is updated at each iteration j until the residual force in the system is less than the convergence criteria. For the material nonlinear model, at each time instant, the stiffness of the system remains constant. However, for the geometric nonlinearity line model, on each time instant i at each iteration j, the stiffness is updated according to the node position and element strain. Consequently, on each time instant i and at each iteration j, the stiffness of the system will be updated. Algorithm 1 shows the steps to solve equation (7) using nonlinear truss element.Algorithm 1: Nonlinear Newmark-Beta method for nonlinear truss element formulation.
Input:Unstrained element length:L0; Poisson’s ratio: ν; Initial element cross section area: A0;Initial element strain:σ0; Elastic modulus: E; Line initial coordinate coord0;Mass matrix:ML; Initial displacement, velocity, and acceleration matrix u0, u˙0, and u¨0;Damping matrixC; Convergence criteria: epsOutput:displacementui+1; velocity u˙i+1; acceleration u¨i+1Begin:1.0 Initial parameters calculation1.1 State determinations: the forcefs0 and the initial tangent stiffness matrix KLT01.2 SolveMLu¨0=P0−Cu˙0−KLT0u0 ⟶ u¨01.3 Select Newmark-Beta parametersγ and β, and time interval dt1.4 Calculatea1=1/βdt2ML+γ/βdtC; a2=1/βdtML+γ/β−1C2.0 Calculations for each time instant,i=0,1,2,…2.1 Initializej=1, ui+1j=ui, fsi+1j=fsi and KLTi+1j=KLTi2.1P^i+1=Pi+1+a1ui+a2u˙i+MLu¨i3.0 For each iteration,j=1,2,3,…3.1 The residual forceR^i+1j=P^i+1−fsi+1j−a1ui+1j3.2 Check convergence; ifR^i+1j<eps, skip the following steps and go to step 4.0; otherwise, implement the following steps:3.3K^LTi+1j=KLTi+1j+a13.4 SolveK^LTi+1jΔui+1j=R^i+1j ⟶ Δui+1j3.5ui+1j+1=ui+1j+Δui+1j3.6 Update the line coordinatecoordi+1j+1=coordi+1j+Δui+1j3.7 Updatefsi+1j+1=fsi+1j+KLTi+1jΔui+1j3.8 Update the line tangent stiffness3.8.1 Fromcoordi+1j+1, calculate the element length Lj+13.8.2 Calculate element length incremental:dL=Lj+1−L03.8.3 Computing axial elongation factorλ using equation (5a)3.8.4 Computing logarithmic strainϵ using equation (5b)3.8.5 Computing the current state stressσ∗ using equations (5c) and (5d)3.8.6 Computing the cross section areaA∗ using equation (6)3.8.7 Calculate the tangent stiffnessKLT at iteration j using equation (2)3.9 Replacej by j+1 and repeat steps 3.1 to 3.8; after converge, denote final displacement value as ui+1, and the coordinate of the line as coordi+14.0 Calculate the velocity and acceleration for time instanti+1u˙i+1=γ/βdtui+1−ui+1−γ/βu˙i+dt1−γ/βu¨i,u¨i+1=1/βdt2ui+1−ui−γ/βdtu˙i+1/2β−1u¨i5.0 Replacei by i+1 and implement steps 2.0 to 4.0 for time instant i+1∗KLT: the tangent stiffness of the transmission line using nonlinear truss element; L refers to transmission line; T refers to truss element


#### 2.1.2. Elastic Catenary Cable Element

The basic assumptions of the elastic catenary element are as follows: (1) the strain stress relationship obeys Hooke’s law, and the cross-sectional area remains unchanged; (2) the geometric nonlinearity considered is a small strain with large deformations; and (3) the cable is perfectly flexible [25]. Figure 1 shows a single elastic catenary element hanging at two nodes I and J. The coordinates of nodes I and J in Cartesian coordinate are 0,0,0 and l1,l2,l3. In this figure, the external loads are thermal load ΔT and uniformly distributed loads w1, w2, and w3. In the Lagrangian coordinate, the deformed and undeformed coordinates of the cable are p and s, respectively.Figure 1
Single cable hanging at two nodes with applied loads.The equations of equilibrium are given by(9)Tdxidp=−wis+fi,i=1,2,3,where wi is the uniformly distributed load in the direction xi and fi is the tension force in the direction xi.The cable tension at Lagrangian coordinate is(10)Ts=∑i=13wis+fi2.The cable tensionT is related to the strain by Hooke’s law as(11)T=EAε−εt=EAdp−dsds−αΔT=EAdpds−1−αΔT,where E is the linear elastic modulus, A is the cross-sectional area, α is the linear thermal expansion coefficient, and ΔT is the temperature change.The relationship between the Cartesian and Lagrange coordinate is expressed as(12)xis=∫0sdxi=∫0sdxidpdpdsds,i=1,2,3.Substituting equations (9)–(11) into (12), the cable projected length xi can be expressed as a function of the undeformed Lagrangian coordinate s:(13)xis=∫0s−wis+fi∑i=13wis+fi2∑i=13wis+fi2EA+1+αΔTds,i=1,2,3.The boundary conditions of the cable are(14a)xilo=li,i=1,2,3,(14b)xi0=0,i=1,2,3,where l0 is the undeformed cable length.Solving equation (13) and applying the boundary conditions in equations (14a) and (14b), the cable length in Cartesian coordinates can be expressed as a function of the internal force fi as(15)li=−l0fiEA−L02wi2EA+1+αΔTw3wwiT1−T2+w2fi−a1wilna1w+T1−lnlow+a1w+T2,(16)w=∑j=13wj2,c1=∑j=13fjwj,T1=T0,T2=Tl0,j=1,2,3.Differentiating both sides of equation (15), the relation between projected length and the internal forces is(17)dli=∑j=13∂li∂fidfi,i=1,2,3.In matrix form, equation (14b) can be expressed as(18)dl1dl2dl3=f11f12f13f21f22f23f31f32f33dF1dF2dF3=FdF1dF2dF3,where F is the symmetric flexibility matrix with terms given by(19)fij=b0i,j−1+αΔTw3b1i,j+b2i,jlna1w+T1−lnlow+a1w+T2,(20a)b0i,j=−lEA,i=j,0,i≠j,(20b)b1i,j=−wwifj+3T2+FjT1+w2fi−a1wiwfj+wjl0w+T2T2l0w2+a1+wT2−wfj+wjT1T1a1+wT1,(20c)b2i,j=−lEA,i=j,0,i≠j,The stiffness matrixk is the inverse of the flexibility matrix.(21)k=F−1.The global tangent stiffness and the internal force vector are determined by the six degrees of freedom matrix and vector as(22)KLC=−kkk−k,(23)Fint=f1,f2,f3,f4,f5,f6T.The nodal forcef4,f5,andf6 can be determined through the force equilibrium in x1,x2,andx3 direction as(24a)f4=−w1l0+f1,(24b)f5=−w2l0+f2,(24c)f6=−w3l0+f3.The stiffness matrix in equation (21) is a function of the internal nodal forces f1,f2,andf3, which is unknown. The internal nodal forces can be initialized based on the equations given by Jayaraman and Knudson [24] and Irvine [26] , which are as follows:(25a)f1=−wl12λ0,(25b)f2=−wl22λ0,(25c)f3=w2−l3coshλ0sinhλ0+l0,in which(26)λ0=106,ifl12+l22=0,0.2,ifl02≤l12+l22+l32,3l02−l32l12+l22−1,ifl02>l12+l22+l32.In other cases, when the initial tensionT0 instead of the unstressed length l0 is given, the stiffness matrix has four unknowns: f1,f2,f3,andl0. The unstressed length is initialized and the forces are modified by(27)l0=l12+l22+l32,(28)fi=−lil0T0,i=1,2,3.Equations (19), (27), and (28) can be solved using Newton–Raphson iteration. The Jacobian of the nonlinear systems of equations is computed from the following equations:(29)H=F∂l∂l0T∂T∂fiT0,(30)∂T∂fiT=f1T,f2T,f3TT,(31)∂l∂l0T=∂l1∂l0,∂l2∂l0,∂l3∂l0T,(32)∂li∂l0=−1+αΔTT2+l0EAfi+3,i=1,2,3.The same procedures are employed to solve the nonlinear equation of motion under external loading using the elastic catenary finite element model. Algorithm2 illustrates the steps to solve equation (7) using elastic catenary cable element.Algorithm 2: Nonlinear Newmark-Beta method for elastic catenary element formulation.
Input:Unstrained element length:L0; Elastic modulus: E; Linear thermal expansion coefficient: αTemperature change:ΔT; Distributed loads in each direction xi: wi,i=1,2,3Line initial coordinatecoord0; Mass matrix ML; Damping matrix C; Initial stiffness matrix KLC0Initial displacement, velocity and acceleration matrixu0, u˙0, and u¨0;Newmark-Beta convergence criteria:eps; Elastic catenary convergence criteria: ϵOutput:displacementui+1; velocity u˙i+1; acceleration u¨i+1Begin:1.0 Initial parameters calculation1.1 State determinations: the forcefs0 and the initial tangent stiffness matrix KLC01.2 SolveMLu¨0=P0−Cu˙0−KLC0u0 ⟶ u01.3 Select Newmark-Beta parametersγ and β, and time interval dt1.4 Calculatea1=1/βdt2ML+γ/βdtCn!/r!n−r!; a2=1/βdtML+γ/β−1C2.0 Calculations for each time instant,i=0,1,2,…2.1 Initializej=1, ui+1j=ui, fsi+1j=fsi and KLCi+1j=KLCi2.1P^i+1=Pi+1+a1ui+a2u˙i+MLu¨i3.0 For each iteration,j=1,2,3,…3.1 The residual forceR^i+1j=P^i+1−fsi+1j−a1ui+1j3.2 Check convergence; ifR^i+1j<eps, skip the following steps and go to step 4.0; otherwise, implement the following steps:3.3K^LTi+1j=KLCi+1j+a13.4 SolveK^LCi+1jΔui+1j=R^i+1j ⟶ Δui+1j3.5ui+1j+1=ui+1j+Δui+1j3.6 Update the line coordinatecoordi+1j+1=coordi+1j+Δui+1j3.7 Updatefsi+1j+1=fsi+1j+KLCi+1jΔui+1j3.8 Update the line tangent stiffness3.8.1 For elementel=1:N3.8.2  Extract the two element ends coordinate fromcoordi+1j+1 as Ix1,x2,x3 and Jx1,x2,x33.8.3  From (16), calculate w3.8.4  From (25c), calculate the internal nodal force vector fi at If1,f2,f33.8.5  Compute target length in each directionL0xi=Jx1,x2,x3−Ix1,x2,x33.8.6  Initialize the convergent parameterdl3.8.7  Whiledl>ϵ3.8.8   Calculate parameterc1 using (16)3.8.9   Calculate the force inxi at J using (24c)3.8.10   Calculate the end force at nodeI and J using (10)3.8.11    Calculate the projection lengthLxi of the element at Cartesian coordinates using (15)3.8.12   Calculatedl=Lxi−L0xi3.8.13   Calculate the flexibility matrixFe using (19)3.8.14   Update the internal force incrementaldf=Fe/dl3.8.15   Update the internal force vector at nodeIx1,x2,x3 as fi+df3.8.16   Calculate the stiffness matrixke using (21)3.8.17   Assemble the element stiffnesskE using (22)3.8.18  Assemble the global stiffness matrixKLC3.9 Replacej by j+1 and repeat steps 3.1 to 3.8; after converge, denote final displacement value as ui+1, and the coordinate of the line as coordi+14.0 Calculate the velocity and acceleration for time instanti+1u˙i+1=γ/βdtui+1−ui+1−γ/βu˙i+dt1−γ/βu¨i,u¨i+1=1/βdt2ui+1−ui−γ/βdtu˙i+1/2β−1u¨i5.0 Replacei by i+1 and implement steps 2.0 to 4.0 for time instant i+1∗KLC: the tangent stiffness of the transmission line using nonlinear truss element; L refers to transmission line; C refers to catenary elementTo obtain the tangent stiffness of the transmission line in Algorithm2, during each Newmark-Beta iteration, the elastic catenary finite element method needs to iteratively calculate the flexibility matrix first and then take the inverse to obtain the tangent stiffness matrix for each element and then assemble the global stiffness matrix. The additional iterations and the matrix inverse operations will slow down the computational speed when the number of elements becomes large. However, in Algorithm 1 for the nonlinear truss element formulation, the stiffness matrix is directly calculated based on the current state transmission line coordinate and strains.


### 2.2. Transmission Tower Modeling

The transmission tower is a three-dimensional structure with hundreds of elements. In ANSYS, the tower is modeled based on its actual material and geometry properties. However, the implementation of detailed transmission tower model in MATLAB is time-consuming. Li et al. [33] proposed that the linear lumped mass system can be employed to model the transmission towers; therefore, the complex finite element transmission tower model can be reduced to the lumped mass model with mass concentrated at N critical deformation points. The finite element model and the lumped mass model are shown in Figure 2. The stiffness of the tower is extracted from the ANSYS; first, the flexibility matrix F is computed, and then the stiffness matrix is obtained by the inverse flexibility matrix KT=F−1.Figure 2
Single transmission tower model. (a) ANSYS/LS-DYNA model. (b) MATLAB lumped mass model.
(a)(b)In the ANSYS model, the mass of the transmission tower is distributed in each element, while in the lumped mass model, the mass is concentrated at each node. To find the concentrated mass on each node, eigen analysis is conducted. The procedures to find the mass of each node are as follows:(1)
From ANSYS modal analysis, the first 3N frequencies as diagonal element ofΩ are extracted.(2)
In the lumped mass model, assuming for each node, the mass is the total mass divided by the number of nodes asm=mt/N. After the mass of each node is calculated, the assumed mass matrix is assembled to MA.(3)
The eigenvectors of the mass and stiffness matrixMA,K are calculated as V.(4)
The mass matrix is updated asMT=KVΩ−1V−1.Here,K is the stiffness matrix extracted from ANSYS from the abovementioned method; N is the number of nodes; MA is the assumed tower mass matrix; and MT is the calculated tower mass matrix. The fidelity of the tower lumped mass model is validated through static analysis and dynamic response under static and dynamic wind loading in Section 4.2.


### 2.3. Transmission Tower-Line Model Development

The coupling effects between the transmission tower-line system cannot be ignored during the static and dynamic analyses [18]. For this reason, the coupling terms must be taken into consideration during the formulating of the tower-line system global stiffness and mass matrix. In the numerical and ANSYS models, the transmission tower and line is considered as pinned connection [18]. At the tower and line coupled nodes, the stiffness and mass in each direction is directly added as KC and MC. With the transmission tower-line system assembled, the system will deform to balance the transmission line end forces under gravity loading. Hence, the static analysis of the transmission tower-line system under gravity loading should be conducted first to establish the system initial state before the static and dynamic analyses [15]. Figure 3 describes the flowchart of the development of the transmission tower-line numerical model procedure.Figure 3
Flowchart of developing coupled transmission tower-line system model.In Figure3, the cable model has two candidates, in which the formulation of the two candidates is derived in Section 2.1. From the 3D nonlinear truss element stiffness formulation derivation, an initial strain should be applied to make the stiffness matrix stable. However, the default commercial software assigned initial strain to the nonlinear truss element model does not reflect the actual strain state of the cable under gravity loading. As a result, an initial strain calibrated nonlinear truss element model is developed by adopting the initial strain and initial shape calculated from elastic catenary finite element model. To validate the accuracy and fast convergence of the calibrated nonlinear truss element model, the ANSYS auto-gravimetric initial strain and shape finding algorithm in ANSYS by Zhang et al. [21] is also implemented.


## 2.1. Transmission-Line Modeling

Transmission line under gravity forms a natural catenary shape [22]. Thus, the shape of the transmission line can be found by the elastic catenary element under loading without approximating the gravity. The two-node nonlinear truss element not having the rotational degree of freedom is widely utilized for transmission line modeling. Nevertheless, the initial strain of the nonlinear truss element has a significant influence on the structural behavior of the transmission line. Therefore, the calculation of the initial strain of the nonlinear truss element needs to be addressed. The initial strain of the nonlinear truss element can be calculated as(1)ε0=FEA0,where ε0 is the initial strain of the nonlinear truss element, F is the transmission line end force, E and A are Young’s modulus and cross-sectional area, and 0 indicates the initial state. Yet, in most cases, instead of the cable end force, the sag, the initial length of the cable is provided. Hence, the method to find the initial shape and initial strain of the cable under gravity loading needs to be explored. For the static analysis, the mass of each node in the transmission line is not used, whereas in the dynamic analysis, the mass matrix ML is used to solve the nonlinear equation of motions. For the node in between two elements, the mass of that node sums up half of the two adjacent element mass, whereas for the boundary nodes, the mass takes half of the associated element.


### 2.1.1. Nonlinear Truss Element

The tension-only two-node three-dimensional truss element is employed both in ANSYS and MATLAB models. The tangent stiffness matrix and internal force vector equations of the nonlinear truss element are derived based on Jürgen and Bathe [29]. With the load applied, the Newton–Raphson iteration is needed to iteratively calculate the current state (∗) tangent stiffness and the internal force vector, which is shown in equations (2) and (3):(2)k∗=dσ∗dλ∗A∗1l0M∗+σ∗A∗1l∗N∗,(3)fi∗=σ∗A∗V∗,where dσ∗ is the stress change, dλ∗ is the strain change, A∗ is the cross-sectional area, l0 is the initial truss length, l∗ is the deformed truss length, and M∗ and N∗ are orthogonal coordinate transformation matrices. In this paper, only the geometric nonlinearity is considered; thus, dσ∗/dλ∗=E. For two-dimensional truss element, M∗, N∗, and V∗ are given by(4a)M∗=c2cs−c2−cscss2−cs−s2−c2−csc2cs−cs−s2css2∗,(4b)N∗=s2−cs−s2cscss2−cs−s2−c2−csc2cs−cs−s2css2∗,(4c)V∗=−c,−s,c,s∗T,where c is cosθ and s is sinθ.When the difference between the internal force vector and the applied load vector is smaller than the defined convergence criteria, the iteration will stop. The tangent stiffness matrix, internal force vector, and the displacement can then be obtained. In equations (2) and (3), the nonlinear strain has different definitions in different material laws. In ANSYS, the strain for large elongation is defined as a logarithmic strain. In this paper, the strain adopts logarithmic strain because the developed open-source transmission system finite element model is compared with ANSYS models. The logarithmic strain and cross section are defined as(5a)λ=ll0,(5b)ε=lnλ,(5c)dσ∗=Eε,(5d)σ∗=dσ∗+σ0k,ifdσ∗+σ0≥0,0,ifdσ∗+σ0<0,(6)A∗=A0λ−2ν,where l is the deformed length, l0 is the original length, λ is the axial elongation factor, ε is the logarithmic strain, dσ∗ is the stress change, σ∗ is the current state stress, σ0 is the initial stress, E is the material elastic modulus, and A∗ is the cross-section area at current state. In equation (5d), the strain is positive definite because the transmission line is tension-only structure.For dynamic analysis, the equation of motion under external loading can be written as(7)MU¨+CU˙+K∗U=Pt,where M, C, and K∗ are the mass, damping, and stiffness matrices; U¨, U˙, and U are the vectors of the acceleration, velocity, and displacement; and Pt is the external loading vector.For the transmission line, the aerodynamic damping greatly influences the dynamic displacement response under large wind speed [30]. Consequently, the damping matrix in equation (7) should take the aerodynamic damping into consideration. Wang et al. [31] proposed a closed-form formulation of the aerodynamic damping ratios for the transmission lines. Therefore, to calculate the aerodynamic damping, the method proposed by Wang et al. [31] is adopted. The detailed derivation to calculate the aerodynamic damping can be found in [31]. However, in [31], instead of the explicitly calculating the damping matrix, the aerodynamic damping modal ratio is calculated. To calculate the damping matrix, the following procedures are applied:(1)
From [31], the in-plane and out-of-plane frequencies Ω and the corresponding mode shapes V are calculated and assembled the mass matrix as M(2)
From [31], the in-plane and out-of-plane direction modal aerodynamic damping ratios are calculated as ξ(3)
The aerodynamic damping matrix can be calculated asCA=MVM−12ξMΩMV′MBy considering the Rayleigh damping, the damping matrix can be represented as(8a)CR=αMM+αKK,(8b)CA=MVM−12ξMΩMV′M,(8c)C=CR+CA,where αM and αK are the mass and stiffness proportional damping coefficients, respectively.To solve the equation of motion, nonlinear Newmark-Beta method is employed [32]. In Newmark-Beta method, at each time instant i, the displacement is updated at each iteration j until the residual force in the system is less than the convergence criteria. For the material nonlinear model, at each time instant, the stiffness of the system remains constant. However, for the geometric nonlinearity line model, on each time instant i at each iteration j, the stiffness is updated according to the node position and element strain. Consequently, on each time instant i and at each iteration j, the stiffness of the system will be updated. Algorithm 1 shows the steps to solve equation (7) using nonlinear truss element.Algorithm 1: Nonlinear Newmark-Beta method for nonlinear truss element formulation.
Input:Unstrained element length:L0; Poisson’s ratio: ν; Initial element cross section area: A0;Initial element strain:σ0; Elastic modulus: E; Line initial coordinate coord0;Mass matrix:ML; Initial displacement, velocity, and acceleration matrix u0, u˙0, and u¨0;Damping matrixC; Convergence criteria: epsOutput:displacementui+1; velocity u˙i+1; acceleration u¨i+1Begin:1.0 Initial parameters calculation1.1 State determinations: the forcefs0 and the initial tangent stiffness matrix KLT01.2 SolveMLu¨0=P0−Cu˙0−KLT0u0 ⟶ u¨01.3 Select Newmark-Beta parametersγ and β, and time interval dt1.4 Calculatea1=1/βdt2ML+γ/βdtC; a2=1/βdtML+γ/β−1C2.0 Calculations for each time instant,i=0,1,2,…2.1 Initializej=1, ui+1j=ui, fsi+1j=fsi and KLTi+1j=KLTi2.1P^i+1=Pi+1+a1ui+a2u˙i+MLu¨i3.0 For each iteration,j=1,2,3,…3.1 The residual forceR^i+1j=P^i+1−fsi+1j−a1ui+1j3.2 Check convergence; ifR^i+1j<eps, skip the following steps and go to step 4.0; otherwise, implement the following steps:3.3K^LTi+1j=KLTi+1j+a13.4 SolveK^LTi+1jΔui+1j=R^i+1j ⟶ Δui+1j3.5ui+1j+1=ui+1j+Δui+1j3.6 Update the line coordinatecoordi+1j+1=coordi+1j+Δui+1j3.7 Updatefsi+1j+1=fsi+1j+KLTi+1jΔui+1j3.8 Update the line tangent stiffness3.8.1 Fromcoordi+1j+1, calculate the element length Lj+13.8.2 Calculate element length incremental:dL=Lj+1−L03.8.3 Computing axial elongation factorλ using equation (5a)3.8.4 Computing logarithmic strainϵ using equation (5b)3.8.5 Computing the current state stressσ∗ using equations (5c) and (5d)3.8.6 Computing the cross section areaA∗ using equation (6)3.8.7 Calculate the tangent stiffnessKLT at iteration j using equation (2)3.9 Replacej by j+1 and repeat steps 3.1 to 3.8; after converge, denote final displacement value as ui+1, and the coordinate of the line as coordi+14.0 Calculate the velocity and acceleration for time instanti+1u˙i+1=γ/βdtui+1−ui+1−γ/βu˙i+dt1−γ/βu¨i,u¨i+1=1/βdt2ui+1−ui−γ/βdtu˙i+1/2β−1u¨i5.0 Replacei by i+1 and implement steps 2.0 to 4.0 for time instant i+1∗KLT: the tangent stiffness of the transmission line using nonlinear truss element; L refers to transmission line; T refers to truss element


### 2.1.2. Elastic Catenary Cable Element

The basic assumptions of the elastic catenary element are as follows: (1) the strain stress relationship obeys Hooke’s law, and the cross-sectional area remains unchanged; (2) the geometric nonlinearity considered is a small strain with large deformations; and (3) the cable is perfectly flexible [25]. Figure 1 shows a single elastic catenary element hanging at two nodes I and J. The coordinates of nodes I and J in Cartesian coordinate are 0,0,0 and l1,l2,l3. In this figure, the external loads are thermal load ΔT and uniformly distributed loads w1, w2, and w3. In the Lagrangian coordinate, the deformed and undeformed coordinates of the cable are p and s, respectively.Figure 1
Single cable hanging at two nodes with applied loads.The equations of equilibrium are given by(9)Tdxidp=−wis+fi,i=1,2,3,where wi is the uniformly distributed load in the direction xi and fi is the tension force in the direction xi.The cable tension at Lagrangian coordinate is(10)Ts=∑i=13wis+fi2.The cable tensionT is related to the strain by Hooke’s law as(11)T=EAε−εt=EAdp−dsds−αΔT=EAdpds−1−αΔT,where E is the linear elastic modulus, A is the cross-sectional area, α is the linear thermal expansion coefficient, and ΔT is the temperature change.The relationship between the Cartesian and Lagrange coordinate is expressed as(12)xis=∫0sdxi=∫0sdxidpdpdsds,i=1,2,3.Substituting equations (9)–(11) into (12), the cable projected length xi can be expressed as a function of the undeformed Lagrangian coordinate s:(13)xis=∫0s−wis+fi∑i=13wis+fi2∑i=13wis+fi2EA+1+αΔTds,i=1,2,3.The boundary conditions of the cable are(14a)xilo=li,i=1,2,3,(14b)xi0=0,i=1,2,3,where l0 is the undeformed cable length.Solving equation (13) and applying the boundary conditions in equations (14a) and (14b), the cable length in Cartesian coordinates can be expressed as a function of the internal force fi as(15)li=−l0fiEA−L02wi2EA+1+αΔTw3wwiT1−T2+w2fi−a1wilna1w+T1−lnlow+a1w+T2,(16)w=∑j=13wj2,c1=∑j=13fjwj,T1=T0,T2=Tl0,j=1,2,3.Differentiating both sides of equation (15), the relation between projected length and the internal forces is(17)dli=∑j=13∂li∂fidfi,i=1,2,3.In matrix form, equation (14b) can be expressed as(18)dl1dl2dl3=f11f12f13f21f22f23f31f32f33dF1dF2dF3=FdF1dF2dF3,where F is the symmetric flexibility matrix with terms given by(19)fij=b0i,j−1+αΔTw3b1i,j+b2i,jlna1w+T1−lnlow+a1w+T2,(20a)b0i,j=−lEA,i=j,0,i≠j,(20b)b1i,j=−wwifj+3T2+FjT1+w2fi−a1wiwfj+wjl0w+T2T2l0w2+a1+wT2−wfj+wjT1T1a1+wT1,(20c)b2i,j=−lEA,i=j,0,i≠j,The stiffness matrixk is the inverse of the flexibility matrix.(21)k=F−1.The global tangent stiffness and the internal force vector are determined by the six degrees of freedom matrix and vector as(22)KLC=−kkk−k,(23)Fint=f1,f2,f3,f4,f5,f6T.The nodal forcef4,f5,andf6 can be determined through the force equilibrium in x1,x2,andx3 direction as(24a)f4=−w1l0+f1,(24b)f5=−w2l0+f2,(24c)f6=−w3l0+f3.The stiffness matrix in equation (21) is a function of the internal nodal forces f1,f2,andf3, which is unknown. The internal nodal forces can be initialized based on the equations given by Jayaraman and Knudson [24] and Irvine [26] , which are as follows:(25a)f1=−wl12λ0,(25b)f2=−wl22λ0,(25c)f3=w2−l3coshλ0sinhλ0+l0,in which(26)λ0=106,ifl12+l22=0,0.2,ifl02≤l12+l22+l32,3l02−l32l12+l22−1,ifl02>l12+l22+l32.In other cases, when the initial tensionT0 instead of the unstressed length l0 is given, the stiffness matrix has four unknowns: f1,f2,f3,andl0. The unstressed length is initialized and the forces are modified by(27)l0=l12+l22+l32,(28)fi=−lil0T0,i=1,2,3.Equations (19), (27), and (28) can be solved using Newton–Raphson iteration. The Jacobian of the nonlinear systems of equations is computed from the following equations:(29)H=F∂l∂l0T∂T∂fiT0,(30)∂T∂fiT=f1T,f2T,f3TT,(31)∂l∂l0T=∂l1∂l0,∂l2∂l0,∂l3∂l0T,(32)∂li∂l0=−1+αΔTT2+l0EAfi+3,i=1,2,3.The same procedures are employed to solve the nonlinear equation of motion under external loading using the elastic catenary finite element model. Algorithm2 illustrates the steps to solve equation (7) using elastic catenary cable element.Algorithm 2: Nonlinear Newmark-Beta method for elastic catenary element formulation.
Input:Unstrained element length:L0; Elastic modulus: E; Linear thermal expansion coefficient: αTemperature change:ΔT; Distributed loads in each direction xi: wi,i=1,2,3Line initial coordinatecoord0; Mass matrix ML; Damping matrix C; Initial stiffness matrix KLC0Initial displacement, velocity and acceleration matrixu0, u˙0, and u¨0;Newmark-Beta convergence criteria:eps; Elastic catenary convergence criteria: ϵOutput:displacementui+1; velocity u˙i+1; acceleration u¨i+1Begin:1.0 Initial parameters calculation1.1 State determinations: the forcefs0 and the initial tangent stiffness matrix KLC01.2 SolveMLu¨0=P0−Cu˙0−KLC0u0 ⟶ u01.3 Select Newmark-Beta parametersγ and β, and time interval dt1.4 Calculatea1=1/βdt2ML+γ/βdtCn!/r!n−r!; a2=1/βdtML+γ/β−1C2.0 Calculations for each time instant,i=0,1,2,…2.1 Initializej=1, ui+1j=ui, fsi+1j=fsi and KLCi+1j=KLCi2.1P^i+1=Pi+1+a1ui+a2u˙i+MLu¨i3.0 For each iteration,j=1,2,3,…3.1 The residual forceR^i+1j=P^i+1−fsi+1j−a1ui+1j3.2 Check convergence; ifR^i+1j<eps, skip the following steps and go to step 4.0; otherwise, implement the following steps:3.3K^LTi+1j=KLCi+1j+a13.4 SolveK^LCi+1jΔui+1j=R^i+1j ⟶ Δui+1j3.5ui+1j+1=ui+1j+Δui+1j3.6 Update the line coordinatecoordi+1j+1=coordi+1j+Δui+1j3.7 Updatefsi+1j+1=fsi+1j+KLCi+1jΔui+1j3.8 Update the line tangent stiffness3.8.1 For elementel=1:N3.8.2  Extract the two element ends coordinate fromcoordi+1j+1 as Ix1,x2,x3 and Jx1,x2,x33.8.3  From (16), calculate w3.8.4  From (25c), calculate the internal nodal force vector fi at If1,f2,f33.8.5  Compute target length in each directionL0xi=Jx1,x2,x3−Ix1,x2,x33.8.6  Initialize the convergent parameterdl3.8.7  Whiledl>ϵ3.8.8   Calculate parameterc1 using (16)3.8.9   Calculate the force inxi at J using (24c)3.8.10   Calculate the end force at nodeI and J using (10)3.8.11    Calculate the projection lengthLxi of the element at Cartesian coordinates using (15)3.8.12   Calculatedl=Lxi−L0xi3.8.13   Calculate the flexibility matrixFe using (19)3.8.14   Update the internal force incrementaldf=Fe/dl3.8.15   Update the internal force vector at nodeIx1,x2,x3 as fi+df3.8.16   Calculate the stiffness matrixke using (21)3.8.17   Assemble the element stiffnesskE using (22)3.8.18  Assemble the global stiffness matrixKLC3.9 Replacej by j+1 and repeat steps 3.1 to 3.8; after converge, denote final displacement value as ui+1, and the coordinate of the line as coordi+14.0 Calculate the velocity and acceleration for time instanti+1u˙i+1=γ/βdtui+1−ui+1−γ/βu˙i+dt1−γ/βu¨i,u¨i+1=1/βdt2ui+1−ui−γ/βdtu˙i+1/2β−1u¨i5.0 Replacei by i+1 and implement steps 2.0 to 4.0 for time instant i+1∗KLC: the tangent stiffness of the transmission line using nonlinear truss element; L refers to transmission line; C refers to catenary elementTo obtain the tangent stiffness of the transmission line in Algorithm2, during each Newmark-Beta iteration, the elastic catenary finite element method needs to iteratively calculate the flexibility matrix first and then take the inverse to obtain the tangent stiffness matrix for each element and then assemble the global stiffness matrix. The additional iterations and the matrix inverse operations will slow down the computational speed when the number of elements becomes large. However, in Algorithm 1 for the nonlinear truss element formulation, the stiffness matrix is directly calculated based on the current state transmission line coordinate and strains.


## 2.1.1. Nonlinear Truss Element

The tension-only two-node three-dimensional truss element is employed both in ANSYS and MATLAB models. The tangent stiffness matrix and internal force vector equations of the nonlinear truss element are derived based on Jürgen and Bathe [29]. With the load applied, the Newton–Raphson iteration is needed to iteratively calculate the current state (∗) tangent stiffness and the internal force vector, which is shown in equations (2) and (3):(2)k∗=dσ∗dλ∗A∗1l0M∗+σ∗A∗1l∗N∗,(3)fi∗=σ∗A∗V∗,where dσ∗ is the stress change, dλ∗ is the strain change, A∗ is the cross-sectional area, l0 is the initial truss length, l∗ is the deformed truss length, and M∗ and N∗ are orthogonal coordinate transformation matrices. In this paper, only the geometric nonlinearity is considered; thus, dσ∗/dλ∗=E. For two-dimensional truss element, M∗, N∗, and V∗ are given by(4a)M∗=c2cs−c2−cscss2−cs−s2−c2−csc2cs−cs−s2css2∗,(4b)N∗=s2−cs−s2cscss2−cs−s2−c2−csc2cs−cs−s2css2∗,(4c)V∗=−c,−s,c,s∗T,where c is cosθ and s is sinθ.When the difference between the internal force vector and the applied load vector is smaller than the defined convergence criteria, the iteration will stop. The tangent stiffness matrix, internal force vector, and the displacement can then be obtained. In equations (2) and (3), the nonlinear strain has different definitions in different material laws. In ANSYS, the strain for large elongation is defined as a logarithmic strain. In this paper, the strain adopts logarithmic strain because the developed open-source transmission system finite element model is compared with ANSYS models. The logarithmic strain and cross section are defined as(5a)λ=ll0,(5b)ε=lnλ,(5c)dσ∗=Eε,(5d)σ∗=dσ∗+σ0k,ifdσ∗+σ0≥0,0,ifdσ∗+σ0<0,(6)A∗=A0λ−2ν,where l is the deformed length, l0 is the original length, λ is the axial elongation factor, ε is the logarithmic strain, dσ∗ is the stress change, σ∗ is the current state stress, σ0 is the initial stress, E is the material elastic modulus, and A∗ is the cross-section area at current state. In equation (5d), the strain is positive definite because the transmission line is tension-only structure.For dynamic analysis, the equation of motion under external loading can be written as(7)MU¨+CU˙+K∗U=Pt,where M, C, and K∗ are the mass, damping, and stiffness matrices; U¨, U˙, and U are the vectors of the acceleration, velocity, and displacement; and Pt is the external loading vector.For the transmission line, the aerodynamic damping greatly influences the dynamic displacement response under large wind speed [30]. Consequently, the damping matrix in equation (7) should take the aerodynamic damping into consideration. Wang et al. [31] proposed a closed-form formulation of the aerodynamic damping ratios for the transmission lines. Therefore, to calculate the aerodynamic damping, the method proposed by Wang et al. [31] is adopted. The detailed derivation to calculate the aerodynamic damping can be found in [31]. However, in [31], instead of the explicitly calculating the damping matrix, the aerodynamic damping modal ratio is calculated. To calculate the damping matrix, the following procedures are applied:(1)
From [31], the in-plane and out-of-plane frequencies Ω and the corresponding mode shapes V are calculated and assembled the mass matrix as M(2)
From [31], the in-plane and out-of-plane direction modal aerodynamic damping ratios are calculated as ξ(3)
The aerodynamic damping matrix can be calculated asCA=MVM−12ξMΩMV′MBy considering the Rayleigh damping, the damping matrix can be represented as(8a)CR=αMM+αKK,(8b)CA=MVM−12ξMΩMV′M,(8c)C=CR+CA,where αM and αK are the mass and stiffness proportional damping coefficients, respectively.To solve the equation of motion, nonlinear Newmark-Beta method is employed [32]. In Newmark-Beta method, at each time instant i, the displacement is updated at each iteration j until the residual force in the system is less than the convergence criteria. For the material nonlinear model, at each time instant, the stiffness of the system remains constant. However, for the geometric nonlinearity line model, on each time instant i at each iteration j, the stiffness is updated according to the node position and element strain. Consequently, on each time instant i and at each iteration j, the stiffness of the system will be updated. Algorithm 1 shows the steps to solve equation (7) using nonlinear truss element.Algorithm 1: Nonlinear Newmark-Beta method for nonlinear truss element formulation.
Input:Unstrained element length:L0; Poisson’s ratio: ν; Initial element cross section area: A0;Initial element strain:σ0; Elastic modulus: E; Line initial coordinate coord0;Mass matrix:ML; Initial displacement, velocity, and acceleration matrix u0, u˙0, and u¨0;Damping matrixC; Convergence criteria: epsOutput:displacementui+1; velocity u˙i+1; acceleration u¨i+1Begin:1.0 Initial parameters calculation1.1 State determinations: the forcefs0 and the initial tangent stiffness matrix KLT01.2 SolveMLu¨0=P0−Cu˙0−KLT0u0 ⟶ u¨01.3 Select Newmark-Beta parametersγ and β, and time interval dt1.4 Calculatea1=1/βdt2ML+γ/βdtC; a2=1/βdtML+γ/β−1C2.0 Calculations for each time instant,i=0,1,2,…2.1 Initializej=1, ui+1j=ui, fsi+1j=fsi and KLTi+1j=KLTi2.1P^i+1=Pi+1+a1ui+a2u˙i+MLu¨i3.0 For each iteration,j=1,2,3,…3.1 The residual forceR^i+1j=P^i+1−fsi+1j−a1ui+1j3.2 Check convergence; ifR^i+1j<eps, skip the following steps and go to step 4.0; otherwise, implement the following steps:3.3K^LTi+1j=KLTi+1j+a13.4 SolveK^LTi+1jΔui+1j=R^i+1j ⟶ Δui+1j3.5ui+1j+1=ui+1j+Δui+1j3.6 Update the line coordinatecoordi+1j+1=coordi+1j+Δui+1j3.7 Updatefsi+1j+1=fsi+1j+KLTi+1jΔui+1j3.8 Update the line tangent stiffness3.8.1 Fromcoordi+1j+1, calculate the element length Lj+13.8.2 Calculate element length incremental:dL=Lj+1−L03.8.3 Computing axial elongation factorλ using equation (5a)3.8.4 Computing logarithmic strainϵ using equation (5b)3.8.5 Computing the current state stressσ∗ using equations (5c) and (5d)3.8.6 Computing the cross section areaA∗ using equation (6)3.8.7 Calculate the tangent stiffnessKLT at iteration j using equation (2)3.9 Replacej by j+1 and repeat steps 3.1 to 3.8; after converge, denote final displacement value as ui+1, and the coordinate of the line as coordi+14.0 Calculate the velocity and acceleration for time instanti+1u˙i+1=γ/βdtui+1−ui+1−γ/βu˙i+dt1−γ/βu¨i,u¨i+1=1/βdt2ui+1−ui−γ/βdtu˙i+1/2β−1u¨i5.0 Replacei by i+1 and implement steps 2.0 to 4.0 for time instant i+1∗KLT: the tangent stiffness of the transmission line using nonlinear truss element; L refers to transmission line; T refers to truss element


## 2.1.2. Elastic Catenary Cable Element

The basic assumptions of the elastic catenary element are as follows: (1) the strain stress relationship obeys Hooke’s law, and the cross-sectional area remains unchanged; (2) the geometric nonlinearity considered is a small strain with large deformations; and (3) the cable is perfectly flexible [25]. Figure 1 shows a single elastic catenary element hanging at two nodes I and J. The coordinates of nodes I and J in Cartesian coordinate are 0,0,0 and l1,l2,l3. In this figure, the external loads are thermal load ΔT and uniformly distributed loads w1, w2, and w3. In the Lagrangian coordinate, the deformed and undeformed coordinates of the cable are p and s, respectively.Figure 1
Single cable hanging at two nodes with applied loads.The equations of equilibrium are given by(9)Tdxidp=−wis+fi,i=1,2,3,where wi is the uniformly distributed load in the direction xi and fi is the tension force in the direction xi.The cable tension at Lagrangian coordinate is(10)Ts=∑i=13wis+fi2.The cable tensionT is related to the strain by Hooke’s law as(11)T=EAε−εt=EAdp−dsds−αΔT=EAdpds−1−αΔT,where E is the linear elastic modulus, A is the cross-sectional area, α is the linear thermal expansion coefficient, and ΔT is the temperature change.The relationship between the Cartesian and Lagrange coordinate is expressed as(12)xis=∫0sdxi=∫0sdxidpdpdsds,i=1,2,3.Substituting equations (9)–(11) into (12), the cable projected length xi can be expressed as a function of the undeformed Lagrangian coordinate s:(13)xis=∫0s−wis+fi∑i=13wis+fi2∑i=13wis+fi2EA+1+αΔTds,i=1,2,3.The boundary conditions of the cable are(14a)xilo=li,i=1,2,3,(14b)xi0=0,i=1,2,3,where l0 is the undeformed cable length.Solving equation (13) and applying the boundary conditions in equations (14a) and (14b), the cable length in Cartesian coordinates can be expressed as a function of the internal force fi as(15)li=−l0fiEA−L02wi2EA+1+αΔTw3wwiT1−T2+w2fi−a1wilna1w+T1−lnlow+a1w+T2,(16)w=∑j=13wj2,c1=∑j=13fjwj,T1=T0,T2=Tl0,j=1,2,3.Differentiating both sides of equation (15), the relation between projected length and the internal forces is(17)dli=∑j=13∂li∂fidfi,i=1,2,3.In matrix form, equation (14b) can be expressed as(18)dl1dl2dl3=f11f12f13f21f22f23f31f32f33dF1dF2dF3=FdF1dF2dF3,where F is the symmetric flexibility matrix with terms given by(19)fij=b0i,j−1+αΔTw3b1i,j+b2i,jlna1w+T1−lnlow+a1w+T2,(20a)b0i,j=−lEA,i=j,0,i≠j,(20b)b1i,j=−wwifj+3T2+FjT1+w2fi−a1wiwfj+wjl0w+T2T2l0w2+a1+wT2−wfj+wjT1T1a1+wT1,(20c)b2i,j=−lEA,i=j,0,i≠j,The stiffness matrixk is the inverse of the flexibility matrix.(21)k=F−1.The global tangent stiffness and the internal force vector are determined by the six degrees of freedom matrix and vector as(22)KLC=−kkk−k,(23)Fint=f1,f2,f3,f4,f5,f6T.The nodal forcef4,f5,andf6 can be determined through the force equilibrium in x1,x2,andx3 direction as(24a)f4=−w1l0+f1,(24b)f5=−w2l0+f2,(24c)f6=−w3l0+f3.The stiffness matrix in equation (21) is a function of the internal nodal forces f1,f2,andf3, which is unknown. The internal nodal forces can be initialized based on the equations given by Jayaraman and Knudson [24] and Irvine [26] , which are as follows:(25a)f1=−wl12λ0,(25b)f2=−wl22λ0,(25c)f3=w2−l3coshλ0sinhλ0+l0,in which(26)λ0=106,ifl12+l22=0,0.2,ifl02≤l12+l22+l32,3l02−l32l12+l22−1,ifl02>l12+l22+l32.In other cases, when the initial tensionT0 instead of the unstressed length l0 is given, the stiffness matrix has four unknowns: f1,f2,f3,andl0. The unstressed length is initialized and the forces are modified by(27)l0=l12+l22+l32,(28)fi=−lil0T0,i=1,2,3.Equations (19), (27), and (28) can be solved using Newton–Raphson iteration. The Jacobian of the nonlinear systems of equations is computed from the following equations:(29)H=F∂l∂l0T∂T∂fiT0,(30)∂T∂fiT=f1T,f2T,f3TT,(31)∂l∂l0T=∂l1∂l0,∂l2∂l0,∂l3∂l0T,(32)∂li∂l0=−1+αΔTT2+l0EAfi+3,i=1,2,3.The same procedures are employed to solve the nonlinear equation of motion under external loading using the elastic catenary finite element model. Algorithm2 illustrates the steps to solve equation (7) using elastic catenary cable element.Algorithm 2: Nonlinear Newmark-Beta method for elastic catenary element formulation.
Input:Unstrained element length:L0; Elastic modulus: E; Linear thermal expansion coefficient: αTemperature change:ΔT; Distributed loads in each direction xi: wi,i=1,2,3Line initial coordinatecoord0; Mass matrix ML; Damping matrix C; Initial stiffness matrix KLC0Initial displacement, velocity and acceleration matrixu0, u˙0, and u¨0;Newmark-Beta convergence criteria:eps; Elastic catenary convergence criteria: ϵOutput:displacementui+1; velocity u˙i+1; acceleration u¨i+1Begin:1.0 Initial parameters calculation1.1 State determinations: the forcefs0 and the initial tangent stiffness matrix KLC01.2 SolveMLu¨0=P0−Cu˙0−KLC0u0 ⟶ u01.3 Select Newmark-Beta parametersγ and β, and time interval dt1.4 Calculatea1=1/βdt2ML+γ/βdtCn!/r!n−r!; a2=1/βdtML+γ/β−1C2.0 Calculations for each time instant,i=0,1,2,…2.1 Initializej=1, ui+1j=ui, fsi+1j=fsi and KLCi+1j=KLCi2.1P^i+1=Pi+1+a1ui+a2u˙i+MLu¨i3.0 For each iteration,j=1,2,3,…3.1 The residual forceR^i+1j=P^i+1−fsi+1j−a1ui+1j3.2 Check convergence; ifR^i+1j<eps, skip the following steps and go to step 4.0; otherwise, implement the following steps:3.3K^LTi+1j=KLCi+1j+a13.4 SolveK^LCi+1jΔui+1j=R^i+1j ⟶ Δui+1j3.5ui+1j+1=ui+1j+Δui+1j3.6 Update the line coordinatecoordi+1j+1=coordi+1j+Δui+1j3.7 Updatefsi+1j+1=fsi+1j+KLCi+1jΔui+1j3.8 Update the line tangent stiffness3.8.1 For elementel=1:N3.8.2  Extract the two element ends coordinate fromcoordi+1j+1 as Ix1,x2,x3 and Jx1,x2,x33.8.3  From (16), calculate w3.8.4  From (25c), calculate the internal nodal force vector fi at If1,f2,f33.8.5  Compute target length in each directionL0xi=Jx1,x2,x3−Ix1,x2,x33.8.6  Initialize the convergent parameterdl3.8.7  Whiledl>ϵ3.8.8   Calculate parameterc1 using (16)3.8.9   Calculate the force inxi at J using (24c)3.8.10   Calculate the end force at nodeI and J using (10)3.8.11    Calculate the projection lengthLxi of the element at Cartesian coordinates using (15)3.8.12   Calculatedl=Lxi−L0xi3.8.13   Calculate the flexibility matrixFe using (19)3.8.14   Update the internal force incrementaldf=Fe/dl3.8.15   Update the internal force vector at nodeIx1,x2,x3 as fi+df3.8.16   Calculate the stiffness matrixke using (21)3.8.17   Assemble the element stiffnesskE using (22)3.8.18  Assemble the global stiffness matrixKLC3.9 Replacej by j+1 and repeat steps 3.1 to 3.8; after converge, denote final displacement value as ui+1, and the coordinate of the line as coordi+14.0 Calculate the velocity and acceleration for time instanti+1u˙i+1=γ/βdtui+1−ui+1−γ/βu˙i+dt1−γ/βu¨i,u¨i+1=1/βdt2ui+1−ui−γ/βdtu˙i+1/2β−1u¨i5.0 Replacei by i+1 and implement steps 2.0 to 4.0 for time instant i+1∗KLC: the tangent stiffness of the transmission line using nonlinear truss element; L refers to transmission line; C refers to catenary elementTo obtain the tangent stiffness of the transmission line in Algorithm2, during each Newmark-Beta iteration, the elastic catenary finite element method needs to iteratively calculate the flexibility matrix first and then take the inverse to obtain the tangent stiffness matrix for each element and then assemble the global stiffness matrix. The additional iterations and the matrix inverse operations will slow down the computational speed when the number of elements becomes large. However, in Algorithm 1 for the nonlinear truss element formulation, the stiffness matrix is directly calculated based on the current state transmission line coordinate and strains.


## 2.2. Transmission Tower Modeling

The transmission tower is a three-dimensional structure with hundreds of elements. In ANSYS, the tower is modeled based on its actual material and geometry properties. However, the implementation of detailed transmission tower model in MATLAB is time-consuming. Li et al. [33] proposed that the linear lumped mass system can be employed to model the transmission towers; therefore, the complex finite element transmission tower model can be reduced to the lumped mass model with mass concentrated at N critical deformation points. The finite element model and the lumped mass model are shown in Figure 2. The stiffness of the tower is extracted from the ANSYS; first, the flexibility matrix F is computed, and then the stiffness matrix is obtained by the inverse flexibility matrix KT=F−1.Figure 2
Single transmission tower model. (a) ANSYS/LS-DYNA model. (b) MATLAB lumped mass model.
(a)(b)In the ANSYS model, the mass of the transmission tower is distributed in each element, while in the lumped mass model, the mass is concentrated at each node. To find the concentrated mass on each node, eigen analysis is conducted. The procedures to find the mass of each node are as follows:(1)
From ANSYS modal analysis, the first 3N frequencies as diagonal element ofΩ are extracted.(2)
In the lumped mass model, assuming for each node, the mass is the total mass divided by the number of nodes asm=mt/N. After the mass of each node is calculated, the assumed mass matrix is assembled to MA.(3)
The eigenvectors of the mass and stiffness matrixMA,K are calculated as V.(4)
The mass matrix is updated asMT=KVΩ−1V−1.Here,K is the stiffness matrix extracted from ANSYS from the abovementioned method; N is the number of nodes; MA is the assumed tower mass matrix; and MT is the calculated tower mass matrix. The fidelity of the tower lumped mass model is validated through static analysis and dynamic response under static and dynamic wind loading in Section 4.2.


## 2.3. Transmission Tower-Line Model Development

The coupling effects between the transmission tower-line system cannot be ignored during the static and dynamic analyses [18]. For this reason, the coupling terms must be taken into consideration during the formulating of the tower-line system global stiffness and mass matrix. In the numerical and ANSYS models, the transmission tower and line is considered as pinned connection [18]. At the tower and line coupled nodes, the stiffness and mass in each direction is directly added as KC and MC. With the transmission tower-line system assembled, the system will deform to balance the transmission line end forces under gravity loading. Hence, the static analysis of the transmission tower-line system under gravity loading should be conducted first to establish the system initial state before the static and dynamic analyses [15]. Figure 3 describes the flowchart of the development of the transmission tower-line numerical model procedure.Figure 3
Flowchart of developing coupled transmission tower-line system model.In Figure3, the cable model has two candidates, in which the formulation of the two candidates is derived in Section 2.1. From the 3D nonlinear truss element stiffness formulation derivation, an initial strain should be applied to make the stiffness matrix stable. However, the default commercial software assigned initial strain to the nonlinear truss element model does not reflect the actual strain state of the cable under gravity loading. As a result, an initial strain calibrated nonlinear truss element model is developed by adopting the initial strain and initial shape calculated from elastic catenary finite element model. To validate the accuracy and fast convergence of the calibrated nonlinear truss element model, the ANSYS auto-gravimetric initial strain and shape finding algorithm in ANSYS by Zhang et al. [21] is also implemented.


## 3. Investigation of Single Transmission Line Model Performance

To first investigate the performance of transmission lines with different element models, both nonlinear truss element model and elastic catenary element model are developed. In the transmission tower-line system, the transmission line contains ground wires and conductors installed in different locations of the tower. Figure4 illustrates the layout and coordinate system of the transmission lines in the transmission tower-line system. In this figure, two groups of transmission lines are referred, namely, the ground wire and conductor. In each transmission line group, except the location differences, the material and span are identical. Consequently, to find the initial shape of the transmission lines, only one transmission line in each group is implemented. For the ground wire, the line is fixed at the coordinate of 0,−3,31.5 meter and 200,−3,31.5 meter is utilized to conduct the static and dynamic analysis, whereas for the conductor, the line placed at coordinate 0,−4.44,26.9 meter and 200,−4.44,26.9 meter is used. Table 1 shows the material properties and the geometry of the ground wire and conductor of transmission lines on account of gravity loading. Here, the sag of the ground wire and conductor is different due to the material property differences. Hence, the initial shape and initial strain of ground wire and conductor need to be found separately.Figure 4
Layout and coordinate system of the transmission lines in the transmission tower-line system.Table 1
Transmission line properties.
Elastic modulus (GPa)Density (kg/m3)Cross-sectional area (m2)Maximum sag (m)Unstrained length (m)Ground wire7846023.29E − 042.9199.913Conductor6717806.05E − 043.7200.033


### 3.1. Nonlinear Truss Element Model Shape Calibration and Initial Strain Calibration for Transmission Line

At initial state, the transmission line deforms to a shape resulted to its self-weight. From cable structure analysis, the deformed shape of the transmission line under gravity loading is catenary with strains caused by the deformation in the model. As Zhang et al. [21] mentioned that the initial strain and initial shape was critical to establish the transmission line finite element model, it is necessary to accurately find the initial shape and the corresponding initial strain of the transmission line under self-weight accurately. Zhang et al. [21] proposed an initial shape and initial strain finding method using nonlinear truss in ANSYS software. The key concepts of the proposed method are as follows: first, a broad grid range of the transmission line initial strains is assumed and the initial strain in ANSYS transmission line model is applied to do autogravimetric analysis; then, the search grid is gradually narrowed down so that the calculated value (for example, the sag of the cable) and the target value difference are less than the predefined convergence criteria (for example, the relative error is less than 1%); and finally, the shape and initial strain of the transmission line are determined. From the key concepts, it is clear that lots of iterations are needed to achieve the predefined convergent criteria. In this paper, elastic catenary finite element model is utilized to find the initial shape and strain of the transmission line under gravity loading. Although the catenary equation or parabolic equation will give similar initial shape and initial strain of the transmission line, the elastic catenary finite element model can further be applied to transmission line and transmission tower-line system static and dynamic analysis.


### 3.2. Initial Shape and Response Assessment

The assessment of single transmission line performance is developed and compared with different models and commercial software packages: (1) elastic catenary cable element in MATLAB; (2) calibrated nonlinear truss element in ANSYS; (3) ANSYS autogravimetric initial shape and strain finding; (4) SAP2000, and (5) ANSYS LS-DYNA model. In this section, the initial shape of the transmission line and their static responses will be discussed first, and later a comparison of transmission dynamic response due to earthquake and wind loading presented.


#### 3.2.1. Shape of Cable Structure under Gravity and Static Wind Loading

The initial shape of the ground wire and conductor under gravity loading is shown in Figures5 and 6. In Table 2, the sag difference between SAP2000 and elastic catenary finite element model is 0.001%. The small sag difference between the two models is made because SAP2000 makes use of the same elastic catenary finite element for the transmission line modeling. Using calibrated nonlinear truss element in ANSYS, the ground wire and conductor difference with respect to elastic catenary finite element is less than 0.02%. However, ANSYS autogravimetric method calculated sag difference corresponding to elastic catenary finite element is 0.50% and 0.71% for the ground wire and conductor, respectively. Compared with the ANSYS autogravimetric method, the calibrated nonlinear truss element gives more accurate results. Besides, to find the initial shape and initial strain of the transmission lines, the ANSYS autogravimetric method tried and compared several initial strains to find the optimal one. Figure 7 illustrates the ANSYS auto-gravimetric method initial strain calculation and initial shape determination procedure. At first round grid search, the ANSYS auto-gravimetric method implements a broad range of initial strains, which is 1.0E−7 to 1.0E−3. The sag of each strain is calculated and then compared with the target sag to determine whether the applied initial strain is small or large. After several strains are tried, a smaller range of initial strains can be determined, which is 1.0E−4 to 1.0E−3. Then, the process is repeated to gradually narrow down the initial strain search range until the desired initial strain is obtained. With the optimal initial strain obtained, the initial shape of the transmission line can be found by using the auto-gravimetric analysis in ANSYS.Figure 5
Ground wire initial shape.Figure 6
Conductor initial shape.Table 2
Transmission line initial shape finding method calculated sag difference with respect to elastic catenary finite element model.
Calculated line sag difference (%)ANSYS auto-gravimetricSAP 2000Calibrated nonlinear truss (ANSYS)Ground wire0.5040.0010.000Conductor0.7080.0010.014Figure 7
ANSYS auto-gravimetric method initial strain finding iterations for the conductor.The input static wind force is computed from the wind speed. Based on ASCE manual 74 [34] and ASCE 7–10 [35], the static wind force in the transverse and longitudinal direction is as follows:(33)Ft=γwQKzKztv2GtcosΨCftAmt,(34)Fl=γwQKzKztv2GtsinΨCflAml,where v is the wind speed at height z; Ψ is the wind incident angle; γw is the load factor; Q is the numerical constant; Kz is the velocity pressure exposure coefficient; Kzt is the topographic factor; G is the tower gust response factor; Amt and Aml is the projection area in transverse and longitudinal direction; and Cft and Cfl is the force coefficient in transverse and longitudinal direction.The wind speed at heightz is governed by the power law as(35)vv¯s=zzsα,where v¯s is the basic wind speed at standard height zs and α is the power law exponent that represents the surface roughness. The basic wind speed defined in ASCE 7-10 [35] is a three-second gust speed at 10m above the ground in Exposure C. Thus, the standard height zs is 10m.In equations (33) and (34), the basic wind speed and wind angle should be determined. In the extreme wind map from ASCE Manual 74 [34], 54m/s is the most dominant basic wind speed in Texas region. Accordingly, the basic wind speed for the transmission line static analysis takes 54m/s and the wind angle is perpendicular to the span of the line, which is the y direction in Figure 4. Since there are 100 elements in the transmission line finite element model with the maximum height difference to be the sag value, which is small comparing with the transmission line length, only several representative nodes are necessary to generate the wind speed and apply the wind force. For the ground wire and conductor, the node at every 20m in x direction is chosen as the representative node; as a result, there are nine representative nodes in the transmission line, which are shown as the dots in Figure 4 for the ground wire and conductor.Figures8 and 9 illustrate the static deformation using different transmission line models and their differences with respect to elastic catenary finite element model. From the results, the SAP2000 model and the nonlinear truss calibrated model differences in all directions are less than 0.1%, whereas the ANSYS auto-gravimetric method maximum difference is more than 4%. While in the initial strain finding, the calculated sag difference between the elastic catenary element and the ANSYS auto-gravimetric method is less than 1%, and the small calculated sag difference enlarges the error in the static analysis using ANSYS auto-gravimetric method.Figure 8
Ground wire static response.Figure 9
Conductor static response.


#### 3.2.2. Dynamic Response of Cable Structure with Earthquake and Dynamic Wind Loading

The displacement time-history analysis of the system under dynamic loading is studied to verify the dynamic properties of the transmission line model. Because the transmission lines cross different terrains, which may be intense wind flow and earthquake zone, it is necessary to compare the different transmission line model displacement responses under seismic and dynamic wind loading. In the static analysis, the ground wire and conductor are compared separately but with the same methodology. For dynamic analysis, only the conductor is taken to investigate the performance of different finite element models. To quantitatively compare the displacement response difference, root-mean-squared-error (RMSE) is used as the error indicator.(36)RMSE=∑n=0N−1xn−yn2N×100%.(1) Earthquake loading: the El-Centro ground motion acceleration data is used as the benchmark to evaluate the performance of the models. The earthquake loading is applied in the transverse direction of the transmission line as the transmission line has smaller stiffness in that direction.Figure10 shows the transmission line middle position displacement time-history in y and z directions. The x direction displacement is not shown because the displacement at that direction is zero due to structural symmetry. In this figure, it is clear that the LS-DYNA model displacement is different from other models. Table 3 quantitatively illustrates that the RMSE between the LS-DYNA model and the elastic catenary element in the y and z is 67.54% and 78.78%, respectively. Those large RMSEs demonstrate that the LS-DYNA model does not obtain the correct transmission line displacement time-history. The reason behind those large RMSEs is that, in the LS-DYNA model, the initial strain of the transmission line is implicitly implemented to the stiffness formulation. From the LS-DYNA displacement time-history, this implicitly applied initial strain is a small value that only aimed at stabilizing the nonlinear truss element stiffness formulation. However, in Table 3, SAP2000 model and the calibrated nonlinear truss model maximum RMSE is less than 4%, demonstrating that those two models are comparable with the elastic catenary finite element model. Additionally, the ANSYS auto-gravimetric method gives acceptable results, but it is clear that the displacement at transverse direction has larger difference, which is 10.89%, compared with SAP2000 and calibrated nonlinear truss model.Figure 10
Conductor middle position displacement time-history under El-Centro earthquake.Table 3
Conductor middle position displacement RMSE between each model and the elastic catenary model under El-Centro earthquake loading.
Transmission modelSAP2000ANSYS auto-gravimetricCalibrated nonlinear trussLS-DYNAY direction (%)0.2510.893.5067.54Z direction (%)0.131.640.1778.78(2) Dynamic wind loading: the dynamic wind speed consists of basic and fluctuating wind speed. According to Davenport [36], the basic wind speed profile over height is governed by the power law equation (35). The fluctuating wind speed is computed from the basic wind speed, in which the spatial and temporal correlation is taken into consideration. Based on the Shinozuka theory [37], the fluctuating wind speed v1t,v2t,…,vnpt at time t is(37)vkyk,zk,t=2Δω∑m=1k∑l=1NHjmωmlcosωmlt−θkmωml+Φml,k=1,2…,np,(38)ωml=l−1Δω+mNΔω,l=1,2,…,N,where N is an arbitrarily large positive number; Δω=ωup/N, the frequency increment; ωup is the cut-off frequency; that is, when ω>ωup, Sω=0; Φml is uniformly distributed random phase angle in 0,2π. θkm is the Hω phase angle. In engineering, Sω and Hω are real matrices; thus, θkm=0. Hjm is the Sω Cholesky decomposed matrix as in (39)Sω=HωH∗ωT.In equation (39), the fluctuating wind spectral density matrix Sω is calculated from Davenport auto-correlation spectrum and cross-correlation spectrum [36].Due to the limitations in SAP2000 software, the dynamic wind loading cannot be applied. For this reason, the elastic catenary, ANSYS auto-gravimetric method, nonlinear truss calibrated model, and LS-DYNA model are compared.In Figure11, the nonlinear truss calibrated model and the elastic catenary model displacement time-history is overlapped, and the ANSYS auto-gravimetric method displacement time-history deviates after some time instant, whereas the LS-DYNA model displacement time-history far off the rest models. In Table 4, the RMSE at the transverse direction for the nonlinear truss calibrated model, the ANSYS auto-gravimetric model, and the LS-DYNA model is 1.14%, 14.93%, and 133.02%, which quantitatively illustrates the accuracy of the models.Figure 11
Conductor middle position displacement time-history under wind loading.Table 4
Conductor middle position displacement RMSE between each model and the elastic catenary model under dynamic wind loading.
Transmission line modelANSYS auto-gravimetric methodNonlinear truss calibratedLS-DYNAx direction (%)0.030.010.19y direction (%)14.931.14133.02z direction (%)8.510.28121.56In summary, from the dynamic wind loading cases, transmission line displacement response magnitude and trend in LS-DYNA is far off from the rest models. Therefore, with the default software parameter settings, LS-DYNA is not suitable for the transmission line and transmission tower-line coupling modeling unless the initial strain of the transmission line can be accurately applied to the model. The ANSYS auto-gravimetric transmission line model gives relative accurate static and dynamic response compared to the elastic catenary finite element model, whereas the nonlinear truss calibrated model structural static and dynamic properties are close to elastic catenary finite element model. The differences between the nonlinear truss calibrated model and ANSYS auto-gravimetric model, both being implemented in ANSYS software, come from the initial strain and the initial shape calculation differences under gravity loading. The initial shape of the transmission line determines the coordinates of the nonlinear truss element. In Algorithm1, the nonlinear truss element stiffness formulation takes the coordinate and the strain of the truss element into consideration. Eventually, the difference in the initial strain and initial shape brings about different nonlinear truss element stiffness.


## 3.1. Nonlinear Truss Element Model Shape Calibration and Initial Strain Calibration for Transmission Line

At initial state, the transmission line deforms to a shape resulted to its self-weight. From cable structure analysis, the deformed shape of the transmission line under gravity loading is catenary with strains caused by the deformation in the model. As Zhang et al. [21] mentioned that the initial strain and initial shape was critical to establish the transmission line finite element model, it is necessary to accurately find the initial shape and the corresponding initial strain of the transmission line under self-weight accurately. Zhang et al. [21] proposed an initial shape and initial strain finding method using nonlinear truss in ANSYS software. The key concepts of the proposed method are as follows: first, a broad grid range of the transmission line initial strains is assumed and the initial strain in ANSYS transmission line model is applied to do autogravimetric analysis; then, the search grid is gradually narrowed down so that the calculated value (for example, the sag of the cable) and the target value difference are less than the predefined convergence criteria (for example, the relative error is less than 1%); and finally, the shape and initial strain of the transmission line are determined. From the key concepts, it is clear that lots of iterations are needed to achieve the predefined convergent criteria. In this paper, elastic catenary finite element model is utilized to find the initial shape and strain of the transmission line under gravity loading. Although the catenary equation or parabolic equation will give similar initial shape and initial strain of the transmission line, the elastic catenary finite element model can further be applied to transmission line and transmission tower-line system static and dynamic analysis.


## 3.2. Initial Shape and Response Assessment

The assessment of single transmission line performance is developed and compared with different models and commercial software packages: (1) elastic catenary cable element in MATLAB; (2) calibrated nonlinear truss element in ANSYS; (3) ANSYS autogravimetric initial shape and strain finding; (4) SAP2000, and (5) ANSYS LS-DYNA model. In this section, the initial shape of the transmission line and their static responses will be discussed first, and later a comparison of transmission dynamic response due to earthquake and wind loading presented.


### 3.2.1. Shape of Cable Structure under Gravity and Static Wind Loading

The initial shape of the ground wire and conductor under gravity loading is shown in Figures5 and 6. In Table 2, the sag difference between SAP2000 and elastic catenary finite element model is 0.001%. The small sag difference between the two models is made because SAP2000 makes use of the same elastic catenary finite element for the transmission line modeling. Using calibrated nonlinear truss element in ANSYS, the ground wire and conductor difference with respect to elastic catenary finite element is less than 0.02%. However, ANSYS autogravimetric method calculated sag difference corresponding to elastic catenary finite element is 0.50% and 0.71% for the ground wire and conductor, respectively. Compared with the ANSYS autogravimetric method, the calibrated nonlinear truss element gives more accurate results. Besides, to find the initial shape and initial strain of the transmission lines, the ANSYS autogravimetric method tried and compared several initial strains to find the optimal one. Figure 7 illustrates the ANSYS auto-gravimetric method initial strain calculation and initial shape determination procedure. At first round grid search, the ANSYS auto-gravimetric method implements a broad range of initial strains, which is 1.0E−7 to 1.0E−3. The sag of each strain is calculated and then compared with the target sag to determine whether the applied initial strain is small or large. After several strains are tried, a smaller range of initial strains can be determined, which is 1.0E−4 to 1.0E−3. Then, the process is repeated to gradually narrow down the initial strain search range until the desired initial strain is obtained. With the optimal initial strain obtained, the initial shape of the transmission line can be found by using the auto-gravimetric analysis in ANSYS.Figure 5
Ground wire initial shape.Figure 6
Conductor initial shape.Table 2
Transmission line initial shape finding method calculated sag difference with respect to elastic catenary finite element model.
Calculated line sag difference (%)ANSYS auto-gravimetricSAP 2000Calibrated nonlinear truss (ANSYS)Ground wire0.5040.0010.000Conductor0.7080.0010.014Figure 7
ANSYS auto-gravimetric method initial strain finding iterations for the conductor.The input static wind force is computed from the wind speed. Based on ASCE manual 74 [34] and ASCE 7–10 [35], the static wind force in the transverse and longitudinal direction is as follows:(33)Ft=γwQKzKztv2GtcosΨCftAmt,(34)Fl=γwQKzKztv2GtsinΨCflAml,where v is the wind speed at height z; Ψ is the wind incident angle; γw is the load factor; Q is the numerical constant; Kz is the velocity pressure exposure coefficient; Kzt is the topographic factor; G is the tower gust response factor; Amt and Aml is the projection area in transverse and longitudinal direction; and Cft and Cfl is the force coefficient in transverse and longitudinal direction.The wind speed at heightz is governed by the power law as(35)vv¯s=zzsα,where v¯s is the basic wind speed at standard height zs and α is the power law exponent that represents the surface roughness. The basic wind speed defined in ASCE 7-10 [35] is a three-second gust speed at 10m above the ground in Exposure C. Thus, the standard height zs is 10m.In equations (33) and (34), the basic wind speed and wind angle should be determined. In the extreme wind map from ASCE Manual 74 [34], 54m/s is the most dominant basic wind speed in Texas region. Accordingly, the basic wind speed for the transmission line static analysis takes 54m/s and the wind angle is perpendicular to the span of the line, which is the y direction in Figure 4. Since there are 100 elements in the transmission line finite element model with the maximum height difference to be the sag value, which is small comparing with the transmission line length, only several representative nodes are necessary to generate the wind speed and apply the wind force. For the ground wire and conductor, the node at every 20m in x direction is chosen as the representative node; as a result, there are nine representative nodes in the transmission line, which are shown as the dots in Figure 4 for the ground wire and conductor.Figures8 and 9 illustrate the static deformation using different transmission line models and their differences with respect to elastic catenary finite element model. From the results, the SAP2000 model and the nonlinear truss calibrated model differences in all directions are less than 0.1%, whereas the ANSYS auto-gravimetric method maximum difference is more than 4%. While in the initial strain finding, the calculated sag difference between the elastic catenary element and the ANSYS auto-gravimetric method is less than 1%, and the small calculated sag difference enlarges the error in the static analysis using ANSYS auto-gravimetric method.Figure 8
Ground wire static response.Figure 9
Conductor static response.


### 3.2.2. Dynamic Response of Cable Structure with Earthquake and Dynamic Wind Loading

The displacement time-history analysis of the system under dynamic loading is studied to verify the dynamic properties of the transmission line model. Because the transmission lines cross different terrains, which may be intense wind flow and earthquake zone, it is necessary to compare the different transmission line model displacement responses under seismic and dynamic wind loading. In the static analysis, the ground wire and conductor are compared separately but with the same methodology. For dynamic analysis, only the conductor is taken to investigate the performance of different finite element models. To quantitatively compare the displacement response difference, root-mean-squared-error (RMSE) is used as the error indicator.(36)RMSE=∑n=0N−1xn−yn2N×100%.(1) Earthquake loading: the El-Centro ground motion acceleration data is used as the benchmark to evaluate the performance of the models. The earthquake loading is applied in the transverse direction of the transmission line as the transmission line has smaller stiffness in that direction.Figure10 shows the transmission line middle position displacement time-history in y and z directions. The x direction displacement is not shown because the displacement at that direction is zero due to structural symmetry. In this figure, it is clear that the LS-DYNA model displacement is different from other models. Table 3 quantitatively illustrates that the RMSE between the LS-DYNA model and the elastic catenary element in the y and z is 67.54% and 78.78%, respectively. Those large RMSEs demonstrate that the LS-DYNA model does not obtain the correct transmission line displacement time-history. The reason behind those large RMSEs is that, in the LS-DYNA model, the initial strain of the transmission line is implicitly implemented to the stiffness formulation. From the LS-DYNA displacement time-history, this implicitly applied initial strain is a small value that only aimed at stabilizing the nonlinear truss element stiffness formulation. However, in Table 3, SAP2000 model and the calibrated nonlinear truss model maximum RMSE is less than 4%, demonstrating that those two models are comparable with the elastic catenary finite element model. Additionally, the ANSYS auto-gravimetric method gives acceptable results, but it is clear that the displacement at transverse direction has larger difference, which is 10.89%, compared with SAP2000 and calibrated nonlinear truss model.Figure 10
Conductor middle position displacement time-history under El-Centro earthquake.Table 3
Conductor middle position displacement RMSE between each model and the elastic catenary model under El-Centro earthquake loading.
Transmission modelSAP2000ANSYS auto-gravimetricCalibrated nonlinear trussLS-DYNAY direction (%)0.2510.893.5067.54Z direction (%)0.131.640.1778.78(2) Dynamic wind loading: the dynamic wind speed consists of basic and fluctuating wind speed. According to Davenport [36], the basic wind speed profile over height is governed by the power law equation (35). The fluctuating wind speed is computed from the basic wind speed, in which the spatial and temporal correlation is taken into consideration. Based on the Shinozuka theory [37], the fluctuating wind speed v1t,v2t,…,vnpt at time t is(37)vkyk,zk,t=2Δω∑m=1k∑l=1NHjmωmlcosωmlt−θkmωml+Φml,k=1,2…,np,(38)ωml=l−1Δω+mNΔω,l=1,2,…,N,where N is an arbitrarily large positive number; Δω=ωup/N, the frequency increment; ωup is the cut-off frequency; that is, when ω>ωup, Sω=0; Φml is uniformly distributed random phase angle in 0,2π. θkm is the Hω phase angle. In engineering, Sω and Hω are real matrices; thus, θkm=0. Hjm is the Sω Cholesky decomposed matrix as in (39)Sω=HωH∗ωT.In equation (39), the fluctuating wind spectral density matrix Sω is calculated from Davenport auto-correlation spectrum and cross-correlation spectrum [36].Due to the limitations in SAP2000 software, the dynamic wind loading cannot be applied. For this reason, the elastic catenary, ANSYS auto-gravimetric method, nonlinear truss calibrated model, and LS-DYNA model are compared.In Figure11, the nonlinear truss calibrated model and the elastic catenary model displacement time-history is overlapped, and the ANSYS auto-gravimetric method displacement time-history deviates after some time instant, whereas the LS-DYNA model displacement time-history far off the rest models. In Table 4, the RMSE at the transverse direction for the nonlinear truss calibrated model, the ANSYS auto-gravimetric model, and the LS-DYNA model is 1.14%, 14.93%, and 133.02%, which quantitatively illustrates the accuracy of the models.Figure 11
Conductor middle position displacement time-history under wind loading.Table 4
Conductor middle position displacement RMSE between each model and the elastic catenary model under dynamic wind loading.
Transmission line modelANSYS auto-gravimetric methodNonlinear truss calibratedLS-DYNAx direction (%)0.030.010.19y direction (%)14.931.14133.02z direction (%)8.510.28121.56In summary, from the dynamic wind loading cases, transmission line displacement response magnitude and trend in LS-DYNA is far off from the rest models. Therefore, with the default software parameter settings, LS-DYNA is not suitable for the transmission line and transmission tower-line coupling modeling unless the initial strain of the transmission line can be accurately applied to the model. The ANSYS auto-gravimetric transmission line model gives relative accurate static and dynamic response compared to the elastic catenary finite element model, whereas the nonlinear truss calibrated model structural static and dynamic properties are close to elastic catenary finite element model. The differences between the nonlinear truss calibrated model and ANSYS auto-gravimetric model, both being implemented in ANSYS software, come from the initial strain and the initial shape calculation differences under gravity loading. The initial shape of the transmission line determines the coordinates of the nonlinear truss element. In Algorithm1, the nonlinear truss element stiffness formulation takes the coordinate and the strain of the truss element into consideration. Eventually, the difference in the initial strain and initial shape brings about different nonlinear truss element stiffness.


## 3.2.1. Shape of Cable Structure under Gravity and Static Wind Loading

The initial shape of the ground wire and conductor under gravity loading is shown in Figures5 and 6. In Table 2, the sag difference between SAP2000 and elastic catenary finite element model is 0.001%. The small sag difference between the two models is made because SAP2000 makes use of the same elastic catenary finite element for the transmission line modeling. Using calibrated nonlinear truss element in ANSYS, the ground wire and conductor difference with respect to elastic catenary finite element is less than 0.02%. However, ANSYS autogravimetric method calculated sag difference corresponding to elastic catenary finite element is 0.50% and 0.71% for the ground wire and conductor, respectively. Compared with the ANSYS autogravimetric method, the calibrated nonlinear truss element gives more accurate results. Besides, to find the initial shape and initial strain of the transmission lines, the ANSYS autogravimetric method tried and compared several initial strains to find the optimal one. Figure 7 illustrates the ANSYS auto-gravimetric method initial strain calculation and initial shape determination procedure. At first round grid search, the ANSYS auto-gravimetric method implements a broad range of initial strains, which is 1.0E−7 to 1.0E−3. The sag of each strain is calculated and then compared with the target sag to determine whether the applied initial strain is small or large. After several strains are tried, a smaller range of initial strains can be determined, which is 1.0E−4 to 1.0E−3. Then, the process is repeated to gradually narrow down the initial strain search range until the desired initial strain is obtained. With the optimal initial strain obtained, the initial shape of the transmission line can be found by using the auto-gravimetric analysis in ANSYS.Figure 5
Ground wire initial shape.Figure 6
Conductor initial shape.Table 2
Transmission line initial shape finding method calculated sag difference with respect to elastic catenary finite element model.
Calculated line sag difference (%)ANSYS auto-gravimetricSAP 2000Calibrated nonlinear truss (ANSYS)Ground wire0.5040.0010.000Conductor0.7080.0010.014Figure 7
ANSYS auto-gravimetric method initial strain finding iterations for the conductor.The input static wind force is computed from the wind speed. Based on ASCE manual 74 [34] and ASCE 7–10 [35], the static wind force in the transverse and longitudinal direction is as follows:(33)Ft=γwQKzKztv2GtcosΨCftAmt,(34)Fl=γwQKzKztv2GtsinΨCflAml,where v is the wind speed at height z; Ψ is the wind incident angle; γw is the load factor; Q is the numerical constant; Kz is the velocity pressure exposure coefficient; Kzt is the topographic factor; G is the tower gust response factor; Amt and Aml is the projection area in transverse and longitudinal direction; and Cft and Cfl is the force coefficient in transverse and longitudinal direction.The wind speed at heightz is governed by the power law as(35)vv¯s=zzsα,where v¯s is the basic wind speed at standard height zs and α is the power law exponent that represents the surface roughness. The basic wind speed defined in ASCE 7-10 [35] is a three-second gust speed at 10m above the ground in Exposure C. Thus, the standard height zs is 10m.In equations (33) and (34), the basic wind speed and wind angle should be determined. In the extreme wind map from ASCE Manual 74 [34], 54m/s is the most dominant basic wind speed in Texas region. Accordingly, the basic wind speed for the transmission line static analysis takes 54m/s and the wind angle is perpendicular to the span of the line, which is the y direction in Figure 4. Since there are 100 elements in the transmission line finite element model with the maximum height difference to be the sag value, which is small comparing with the transmission line length, only several representative nodes are necessary to generate the wind speed and apply the wind force. For the ground wire and conductor, the node at every 20m in x direction is chosen as the representative node; as a result, there are nine representative nodes in the transmission line, which are shown as the dots in Figure 4 for the ground wire and conductor.Figures8 and 9 illustrate the static deformation using different transmission line models and their differences with respect to elastic catenary finite element model. From the results, the SAP2000 model and the nonlinear truss calibrated model differences in all directions are less than 0.1%, whereas the ANSYS auto-gravimetric method maximum difference is more than 4%. While in the initial strain finding, the calculated sag difference between the elastic catenary element and the ANSYS auto-gravimetric method is less than 1%, and the small calculated sag difference enlarges the error in the static analysis using ANSYS auto-gravimetric method.Figure 8
Ground wire static response.Figure 9
Conductor static response.


## 3.2.2. Dynamic Response of Cable Structure with Earthquake and Dynamic Wind Loading

The displacement time-history analysis of the system under dynamic loading is studied to verify the dynamic properties of the transmission line model. Because the transmission lines cross different terrains, which may be intense wind flow and earthquake zone, it is necessary to compare the different transmission line model displacement responses under seismic and dynamic wind loading. In the static analysis, the ground wire and conductor are compared separately but with the same methodology. For dynamic analysis, only the conductor is taken to investigate the performance of different finite element models. To quantitatively compare the displacement response difference, root-mean-squared-error (RMSE) is used as the error indicator.(36)RMSE=∑n=0N−1xn−yn2N×100%.(1) Earthquake loading: the El-Centro ground motion acceleration data is used as the benchmark to evaluate the performance of the models. The earthquake loading is applied in the transverse direction of the transmission line as the transmission line has smaller stiffness in that direction.Figure10 shows the transmission line middle position displacement time-history in y and z directions. The x direction displacement is not shown because the displacement at that direction is zero due to structural symmetry. In this figure, it is clear that the LS-DYNA model displacement is different from other models. Table 3 quantitatively illustrates that the RMSE between the LS-DYNA model and the elastic catenary element in the y and z is 67.54% and 78.78%, respectively. Those large RMSEs demonstrate that the LS-DYNA model does not obtain the correct transmission line displacement time-history. The reason behind those large RMSEs is that, in the LS-DYNA model, the initial strain of the transmission line is implicitly implemented to the stiffness formulation. From the LS-DYNA displacement time-history, this implicitly applied initial strain is a small value that only aimed at stabilizing the nonlinear truss element stiffness formulation. However, in Table 3, SAP2000 model and the calibrated nonlinear truss model maximum RMSE is less than 4%, demonstrating that those two models are comparable with the elastic catenary finite element model. Additionally, the ANSYS auto-gravimetric method gives acceptable results, but it is clear that the displacement at transverse direction has larger difference, which is 10.89%, compared with SAP2000 and calibrated nonlinear truss model.Figure 10
Conductor middle position displacement time-history under El-Centro earthquake.Table 3
Conductor middle position displacement RMSE between each model and the elastic catenary model under El-Centro earthquake loading.
Transmission modelSAP2000ANSYS auto-gravimetricCalibrated nonlinear trussLS-DYNAY direction (%)0.2510.893.5067.54Z direction (%)0.131.640.1778.78(2) Dynamic wind loading: the dynamic wind speed consists of basic and fluctuating wind speed. According to Davenport [36], the basic wind speed profile over height is governed by the power law equation (35). The fluctuating wind speed is computed from the basic wind speed, in which the spatial and temporal correlation is taken into consideration. Based on the Shinozuka theory [37], the fluctuating wind speed v1t,v2t,…,vnpt at time t is(37)vkyk,zk,t=2Δω∑m=1k∑l=1NHjmωmlcosωmlt−θkmωml+Φml,k=1,2…,np,(38)ωml=l−1Δω+mNΔω,l=1,2,…,N,where N is an arbitrarily large positive number; Δω=ωup/N, the frequency increment; ωup is the cut-off frequency; that is, when ω>ωup, Sω=0; Φml is uniformly distributed random phase angle in 0,2π. θkm is the Hω phase angle. In engineering, Sω and Hω are real matrices; thus, θkm=0. Hjm is the Sω Cholesky decomposed matrix as in (39)Sω=HωH∗ωT.In equation (39), the fluctuating wind spectral density matrix Sω is calculated from Davenport auto-correlation spectrum and cross-correlation spectrum [36].Due to the limitations in SAP2000 software, the dynamic wind loading cannot be applied. For this reason, the elastic catenary, ANSYS auto-gravimetric method, nonlinear truss calibrated model, and LS-DYNA model are compared.In Figure11, the nonlinear truss calibrated model and the elastic catenary model displacement time-history is overlapped, and the ANSYS auto-gravimetric method displacement time-history deviates after some time instant, whereas the LS-DYNA model displacement time-history far off the rest models. In Table 4, the RMSE at the transverse direction for the nonlinear truss calibrated model, the ANSYS auto-gravimetric model, and the LS-DYNA model is 1.14%, 14.93%, and 133.02%, which quantitatively illustrates the accuracy of the models.Figure 11
Conductor middle position displacement time-history under wind loading.Table 4
Conductor middle position displacement RMSE between each model and the elastic catenary model under dynamic wind loading.
Transmission line modelANSYS auto-gravimetric methodNonlinear truss calibratedLS-DYNAx direction (%)0.030.010.19y direction (%)14.931.14133.02z direction (%)8.510.28121.56In summary, from the dynamic wind loading cases, transmission line displacement response magnitude and trend in LS-DYNA is far off from the rest models. Therefore, with the default software parameter settings, LS-DYNA is not suitable for the transmission line and transmission tower-line coupling modeling unless the initial strain of the transmission line can be accurately applied to the model. The ANSYS auto-gravimetric transmission line model gives relative accurate static and dynamic response compared to the elastic catenary finite element model, whereas the nonlinear truss calibrated model structural static and dynamic properties are close to elastic catenary finite element model. The differences between the nonlinear truss calibrated model and ANSYS auto-gravimetric model, both being implemented in ANSYS software, come from the initial strain and the initial shape calculation differences under gravity loading. The initial shape of the transmission line determines the coordinates of the nonlinear truss element. In Algorithm1, the nonlinear truss element stiffness formulation takes the coordinate and the strain of the truss element into consideration. Eventually, the difference in the initial strain and initial shape brings about different nonlinear truss element stiffness.


## 4. Development of Transmission Tower-Line System Model

In Section3, the single cable lumped mass models have been compared with the detailed line model in ANSYS Transient and LS-DYNA. To develop an accurate transmission tower-line model, the transmission tower dynamic model and the coupling between tower and line requires careful validation. To achieve computational efficiency, the optimal number of elements in both nonlinear truss element cable model and elastic catenary element cable model is studied. To verify the coupling in the developed transmission tower-line system are captured, five models are cross compared: (1) lumped mass transmission tower coupled with elastic catenary cable element line in MATLAB (ECE-MATLAB); (2) lumped mass transmission tower coupled with calibrated nonlinear truss element line in MATLAB (NLTc-MATLAB); (3) transmission tower coupled with calibrated nonlinear truss element line in ANSYS (NLTc-ANSYS); (4) ANSYS auto-gravimetric model (ANSYS Auto-Gravimetric); and (5) ANSYS LS-DYNA model (LS-DYNA). Moreover, as Fei et al. [38] showed that the tower top displacement is a good indicator to reflect the dynamic response of the transmission tower-line system; eventually, the tower top maximum displacement and the displacement time-history root difference indicator are compared with respect to the ANSYS model to verify the accuracy of the lumped mass model. In the transmission tower-line model, the insulators are not considered. Moreover, from the state-of-the-art transmission tower-line models, each transmission line is aggregated to single transmission line instead of split transmission line.


### 4.1. Single Tower Lumped Mass Tower Model

The transmission tower studied is a steel-made, L-shaped cross-sectional suspension tower with a height of 31.5 m. The transmission tower is a design for the Texas region from Tort et al. [39], which is based on the ASCE Manual 74 [34] using the commercial software PLS-TOWER. The properties of the lines are in Table 1, which is provided by Xue et al. [18]. As shown in Figure 1, the lumped mass tower model has 13 nodes. The tower stiffness and mass matrix are calculated based on the method in Section 2.2.The static wind loading calculation is based on Section3.2.1. In the extreme wind map from ASCE Manual 74, 54m/s basic wind speed is the most dominant wind speed in Texas region. Hence, the basic wind speed for the transmission tower static analysis takes 54m/s. The wind angle is chosen as 30° to verify the lumped mass tower model x and y direction stiffness formulation.The lumped mass transmission tower static displacement under static wind loading displacement difference with respect to that from the ANSYS detailed model is shown in Figure12. From this figure, it is clear that, in all directions, the maximum displacement difference is less than 0.15%, indicating that the stiffness matrix in the lumped mass model is with high accuracy.Figure 12
The lumped mass tower model static displacement difference corresponding to ANSYS (%).The dynamic wind force generation method is described in Section2.3. The displacement time-history responses are compared at two selected nodes: one is on the tower main body (Node 3) and the other is at the tower top (Node 12). To quantitatively measure the difference between the transmission tower lumped mass and ANSYS detailed model displacement time-history without the effects of the displacement magnitude, the difference indicator is considered as [40](40)De=∑n=0N−1xn−yn2/N∑n=0N−1xn2/N.In Figure13, the lumped mass model and the ANSYS model node 3 and node 12 displacement time-history in all directions are close to each other. In Table 5, the displacement difference indicator between the lumped mass model and ANSYS model is less than 4%. The small discrepancy implies that the lumped mass model captured the dynamic properties of the ANSYS model. However, in Table 5, the LS-DYNA model displacement time-history minimum difference with respect to the ANSYS model reaches 17.75%, which is relatively large.Figure 13
Single tower dynamic displacement at nodes 3 and 12 under wind loading.Table 5
Single tower dynamic displacement time-historyDe indicator between each model and the ANSYS model.
x direction (%)y direction (%)z direction (%)Node 3Lumped mass model2.951.963.78ANSYS LS-DYNA model29.0426.2745.86Node 12Lumped mass model1.612.033.05ANSYS LS-DYNA model17.7519.2222.23To conclude, from the single tower static and dynamic analysis, the lumped mass can be treated as an alternative model to the ANSYS detailed model. The LS-DYNA model has relatively large displacement difference to the ANSYS model.


### 4.2. Transmission Line Element Size Reduction

The span of the transmission line is200m; therefore, to accurately simulate the transmission line, large number of finite elements will be used. However, the number of elements of the transmission line has a great impact on the computational time and accuracy of the model. Hence, the frequency domain analysis is performed to find the optimal number of elements. Figure 14 shows the frequency components of the ground wire and the conductor in different transmission line models. The transmission line has been divided into 10, 20, 40, and 100 elements. From the frequency differences plot, using 40 elements for all the transmission line models, the first ten frequencies difference is less than 1%. As a result, 40 element numbers is chosen as the optimal element number for the transmission line.Figure 14
Transmission line size reduction.


### 4.3. Tower-Line Coupling Validation on a One-Span Two-Line System

To validate the accurate representation of transmission tower-line coupling in the proposed model, responses of the one-span two-line transmission tower-line system are compared. The two transmission lines are ground wires, whose material properties are shown in Table1. The span between two towers is 200m. Figure 15 illustrates the two different directions of the wind speed. When the wind speed direction is along the longitudinal direction, the wind angle is annotated 0°, and when the wind speed direction is along the transverse direction, the wind angle is annotated 90°. To cross validate the transmission tower-line models, coupling effects are successfully captured and the wind speed 10m/s and wind angle 30° are chosen. The displacement of the left tower and one ground wire middle position displacement time history is compared because those two points are good indicators to reflect the tower and cable dynamic properties.Figure 15
The illustration of the dynamic wind applied force points and the direction.Figure16 shows the structure tower top and ground wire middle position displacement time-history between the different models under wind loading. In this figure, the NLTc-ANSYS and ANSYS auto-gravimetric model displacement time-history are overlapped. Tables 6 and 7 quantitatively measure the displacement time-history difference indicator using equation (40) between the ANSYS auto-gravimetric model and the NLTc-ANSYS model at tower top and line middle position. In Table 6, the tower top displacement difference between the ANSYS auto-gravimetric model and the NLTc-ANSYS model is less than 1%. In Table 7, the ground wire displacement difference between those two models is 4.14% and 8.23% for y and z directions, which are larger than the tower displacement time-history difference. The reason is that the initial strain calculated from the ANSYS auto-gravimetric model is different from the NLTc-ANSYS model, so that the ground wire end force transferred to the tower is different. However, since the stiffness of the tower is much larger than the ground wire, the displacement difference is smaller at the tower top.Figure 16
One-span two-line transmission tower-line model displacement at tower top and cable middle position.Table 6
One-span two-line tower top dynamic displacement time-historyDe difference indicator between each model and the ANSYS model.
x direction (%)y direction (%)z direction (%)ANSYS auto-gravimetric0.73070.71590.9544ECE-MATLAB0.87620.28809.2953NLTc-MATLAB0.87220.32769.2916LS-DYNA12.671816.130018.1636Table 7
One-span two-line line middle dynamic displacement time-historyDe difference indicator between each model and the ANSYS model.
x direction (%)y direction (%)z direction (%)ANSYS auto-gravimetric0.09804.13738.2330ECE-MATLAB0.94750.23651.6704NLTc-MATLAB0.99880.89101.6977LS-DYNA11.781154.350284.4399In Tables6 and 7, the ECE-MATLAB and NLTc-MATLAB displacement difference with respect to each other is small. In Table 6, the tower top displacement difference in the x and y directions are smaller than 1%. However, in the z direction, the displacement difference is around 10%. The reason behind the relatively large RMSE is that, in ANSYS, the geometric nonlinearity is also applied to the transmission tower. With the applied gravity loading in the system, the geometric nonlinearity will stiffen the tower, which leads to the smaller displacement comparing with MATLAB lumped mass model.Additionally, in Figure16, the LS-DYNA models displacement time-history is different from other models. In Tables 6 and 7, the LS-DYNA maximum displacement time-history difference is 18.16% and 84.44% for the tower top and ground wire position, which manifests that the LS-DYNA model does not simulate the structure dynamic behavior with high accuracy.In general, the ECE-MATLAB, NLTc-MATLAB, NLTc-ANSYS, and ANSYS auto-gravimetric model are comparable with each other and with high fidelity, while the errors from transmission line modeling in LS-DYNA have further propagated and impacted largely on the transmission tower top response.


## 4.1. Single Tower Lumped Mass Tower Model

The transmission tower studied is a steel-made, L-shaped cross-sectional suspension tower with a height of 31.5 m. The transmission tower is a design for the Texas region from Tort et al. [39], which is based on the ASCE Manual 74 [34] using the commercial software PLS-TOWER. The properties of the lines are in Table 1, which is provided by Xue et al. [18]. As shown in Figure 1, the lumped mass tower model has 13 nodes. The tower stiffness and mass matrix are calculated based on the method in Section 2.2.The static wind loading calculation is based on Section3.2.1. In the extreme wind map from ASCE Manual 74, 54m/s basic wind speed is the most dominant wind speed in Texas region. Hence, the basic wind speed for the transmission tower static analysis takes 54m/s. The wind angle is chosen as 30° to verify the lumped mass tower model x and y direction stiffness formulation.The lumped mass transmission tower static displacement under static wind loading displacement difference with respect to that from the ANSYS detailed model is shown in Figure12. From this figure, it is clear that, in all directions, the maximum displacement difference is less than 0.15%, indicating that the stiffness matrix in the lumped mass model is with high accuracy.Figure 12
The lumped mass tower model static displacement difference corresponding to ANSYS (%).The dynamic wind force generation method is described in Section2.3. The displacement time-history responses are compared at two selected nodes: one is on the tower main body (Node 3) and the other is at the tower top (Node 12). To quantitatively measure the difference between the transmission tower lumped mass and ANSYS detailed model displacement time-history without the effects of the displacement magnitude, the difference indicator is considered as [40](40)De=∑n=0N−1xn−yn2/N∑n=0N−1xn2/N.In Figure13, the lumped mass model and the ANSYS model node 3 and node 12 displacement time-history in all directions are close to each other. In Table 5, the displacement difference indicator between the lumped mass model and ANSYS model is less than 4%. The small discrepancy implies that the lumped mass model captured the dynamic properties of the ANSYS model. However, in Table 5, the LS-DYNA model displacement time-history minimum difference with respect to the ANSYS model reaches 17.75%, which is relatively large.Figure 13
Single tower dynamic displacement at nodes 3 and 12 under wind loading.Table 5
Single tower dynamic displacement time-historyDe indicator between each model and the ANSYS model.
x direction (%)y direction (%)z direction (%)Node 3Lumped mass model2.951.963.78ANSYS LS-DYNA model29.0426.2745.86Node 12Lumped mass model1.612.033.05ANSYS LS-DYNA model17.7519.2222.23To conclude, from the single tower static and dynamic analysis, the lumped mass can be treated as an alternative model to the ANSYS detailed model. The LS-DYNA model has relatively large displacement difference to the ANSYS model.


## 4.2. Transmission Line Element Size Reduction

The span of the transmission line is200m; therefore, to accurately simulate the transmission line, large number of finite elements will be used. However, the number of elements of the transmission line has a great impact on the computational time and accuracy of the model. Hence, the frequency domain analysis is performed to find the optimal number of elements. Figure 14 shows the frequency components of the ground wire and the conductor in different transmission line models. The transmission line has been divided into 10, 20, 40, and 100 elements. From the frequency differences plot, using 40 elements for all the transmission line models, the first ten frequencies difference is less than 1%. As a result, 40 element numbers is chosen as the optimal element number for the transmission line.Figure 14
Transmission line size reduction.


## 4.3. Tower-Line Coupling Validation on a One-Span Two-Line System

To validate the accurate representation of transmission tower-line coupling in the proposed model, responses of the one-span two-line transmission tower-line system are compared. The two transmission lines are ground wires, whose material properties are shown in Table1. The span between two towers is 200m. Figure 15 illustrates the two different directions of the wind speed. When the wind speed direction is along the longitudinal direction, the wind angle is annotated 0°, and when the wind speed direction is along the transverse direction, the wind angle is annotated 90°. To cross validate the transmission tower-line models, coupling effects are successfully captured and the wind speed 10m/s and wind angle 30° are chosen. The displacement of the left tower and one ground wire middle position displacement time history is compared because those two points are good indicators to reflect the tower and cable dynamic properties.Figure 15
The illustration of the dynamic wind applied force points and the direction.Figure16 shows the structure tower top and ground wire middle position displacement time-history between the different models under wind loading. In this figure, the NLTc-ANSYS and ANSYS auto-gravimetric model displacement time-history are overlapped. Tables 6 and 7 quantitatively measure the displacement time-history difference indicator using equation (40) between the ANSYS auto-gravimetric model and the NLTc-ANSYS model at tower top and line middle position. In Table 6, the tower top displacement difference between the ANSYS auto-gravimetric model and the NLTc-ANSYS model is less than 1%. In Table 7, the ground wire displacement difference between those two models is 4.14% and 8.23% for y and z directions, which are larger than the tower displacement time-history difference. The reason is that the initial strain calculated from the ANSYS auto-gravimetric model is different from the NLTc-ANSYS model, so that the ground wire end force transferred to the tower is different. However, since the stiffness of the tower is much larger than the ground wire, the displacement difference is smaller at the tower top.Figure 16
One-span two-line transmission tower-line model displacement at tower top and cable middle position.Table 6
One-span two-line tower top dynamic displacement time-historyDe difference indicator between each model and the ANSYS model.
x direction (%)y direction (%)z direction (%)ANSYS auto-gravimetric0.73070.71590.9544ECE-MATLAB0.87620.28809.2953NLTc-MATLAB0.87220.32769.2916LS-DYNA12.671816.130018.1636Table 7
One-span two-line line middle dynamic displacement time-historyDe difference indicator between each model and the ANSYS model.
x direction (%)y direction (%)z direction (%)ANSYS auto-gravimetric0.09804.13738.2330ECE-MATLAB0.94750.23651.6704NLTc-MATLAB0.99880.89101.6977LS-DYNA11.781154.350284.4399In Tables6 and 7, the ECE-MATLAB and NLTc-MATLAB displacement difference with respect to each other is small. In Table 6, the tower top displacement difference in the x and y directions are smaller than 1%. However, in the z direction, the displacement difference is around 10%. The reason behind the relatively large RMSE is that, in ANSYS, the geometric nonlinearity is also applied to the transmission tower. With the applied gravity loading in the system, the geometric nonlinearity will stiffen the tower, which leads to the smaller displacement comparing with MATLAB lumped mass model.Additionally, in Figure16, the LS-DYNA models displacement time-history is different from other models. In Tables 6 and 7, the LS-DYNA maximum displacement time-history difference is 18.16% and 84.44% for the tower top and ground wire position, which manifests that the LS-DYNA model does not simulate the structure dynamic behavior with high accuracy.In general, the ECE-MATLAB, NLTc-MATLAB, NLTc-ANSYS, and ANSYS auto-gravimetric model are comparable with each other and with high fidelity, while the errors from transmission line modeling in LS-DYNA have further propagated and impacted largely on the transmission tower top response.


## 5. Performance of Full Transmission Tower-Line System Model during Wind Loadings

Finally, the developed transmission tower-line model is applied to the full setup. The transmission tower is designed with eight transmission lines attached: six conductors and two ground wires. As is illustrated in the one-span two-line model, the LS-DYNA model displacement time-history response deviates from other models. In addition, the performance of transmission tower-line model with calibrated nonlinear truss element line in ANSYS (NLTc-ANSYS) delivers the most accurate result, which is considered as the reference model. Hence, in the full model, the rest of the three different models are compared: (1) lumped mass transmission tower coupled with elastic catenary cable element line in MATLAB (ECE-MATLAB); (2) lumped mass transmission tower coupled with calibrated nonlinear truss element line in MATLAB (NLTc-MATLAB); and (3) ANSYS auto-gravimetric model.The total GPU time in each model has been recorded and compared to illustrate the efficiency of the models under different wind loadings. There are totally six dynamic wind loadings applied to the transmission tower-line models. The wind speeds are10m/s and 20m/s. For each wind speed, three wind angles 0°, 30°, and 90° are used to generate the dynamic wind loading. For the total displacement time-history comparison, only wind speed 20m/s at angle 30° is chosen to illustrate the difference of the proposed models. At each wind speed and wind angle, two nodes in the tower and the cable middle position at each line are chosen to illustrate the displacement time-history. The selected points are illustrated in Figure 17. For dynamic wind loading simulation of the transmission tower-line model, the time interval in ANSYS model is determined at 0.01 s. ANSYS software also automatically interpolates the input force when the convergence is not meet. The time interval for the lumped mass model is selected as 0.0025 s for optimal computational efficiency and accuracy trade-off.Figure 17
The displacement time-history extracting nodes.To quantitatively measure the response of the three models, the tower top maximum displacement and the RMSE between the three different models and the NLTc-ANSYS model are chosen as difference indicators. Figure18 shows the displacement time-history of the full transmission system using different models, given the same dynamic wind loading with wind speed 20m/s at angle 30°. It also demonstrates similar results between four models. To quantify the result accuracy of different models, the tower top displacement and the RMSE between the models under different wind speeds and wind angles are listed in Tables 8 and 9.Figure 18
One-span eight-line nodes 8 and 13 total displacement time-history comparison at wind speed20m/s and angle 30°.Table 8
One-span eight-line tower top maximum displacement difference corresponding to NLTc-ANSYS model at different wind speeds and wind angles.
Maximum displacement (m)AbsoluteDifference (%)NLTc-ANSYSNLTc-MATLABECE-MATLABANSYS auto-gravimetricS10 m/s A000.01320.02980.02930.0023S10 m/s A300.01120.00670.00530.0072S10 m/s A900.04530.09310.08140.0427S20 m/s A000.04550.05340.05280.0024S20 m/s A300.07090.00790.00900.0063S20 m/s A900.16100.08130.09010.0126S: wind speed; A: wind angle.Table 9
One-span eight-line tower top displacement time-history RMSE corresponding to NLTc-ANSYS model at different wind speeds and wind angles.
Tower top displacement RMSE (%)NLTc-MATLABECE-MATLABANSYS auto-gravimetricS10 m/s A000.00760.00720.0053S10 m/s A300.01300.01230.0048S10 m/s A900.17900.15890.1007S20 m/s A000.03910.04010.0122S20 m/s A300.13570.10970.0945S20 m/s A900.68720.79700.1690S: wind speed; A: wind angle.In Table8, the ANSYS auto-gravimetric model maximum tower top displacement difference corresponding to the NLTc-ANSYS model at all tested wind speeds and wind angles is less than 0.05%. Additionally, both the NLTc-MATLAB and ECE-MALTAB model maximum tower top displacement differences are less than 0.1%. The maximum tower top displacement is a single value that is used to determine the accuracy of the models. Consequently, to check the differences between the models at all-time instant, the RMSE is employed.Table9 summarizes the RMSE between the three different models and the NLTc-ANSYS model at different wind speeds and wind angles. In Table 9, for the ANSYS auto-gravimetric model, at the same wind speed 10m/s, the largest RMSE is 0.1% at wind angle 90 and at wind speed 20m/s; as the wind angle increases from 0° to 90°, the RMSE increases from 0.012% to 0.17%. Comparing between different models, the only difference between the ANSYS auto-gravimetric model and the NLTc-ANSYS model is in the transmission line initial strain and initial shape determination. The small initial strain and initial shape difference adds to the differences in the dynamic analysis enlarging the RMSE. However, such an impact can be negligible. For the two proposed open-source models, at wind speed 10m/s, the maximum RMSE errors of NLTc-MATLAB and ECE-MATLAB are 0.18% and 0.16% respectively, which are similar to the performance of the ANSYS auto-gravimetric model. At wind speed 20m/s, the largest RMSE for the NLTc-MATLAB and ECE-MATLAB model is 0.68% and 0.80%, which are 4 and 4.7 times larger than the ANSYS auto-gravimetric model. The relatively large RMSE at a larger wind speed is due to the plastic deformation of the tower in the ANSYS model. It is suggested that nonlinearity should be considered in the lumped mass model, to encounter the plastic deformation, local buckling, and the failure mechanism of the transmission tower in the future model development.Another important aspect to evaluate the transmission tower-line model is the computational efficiency. Table10 summarizes the computational time in each wind speed at each angle and the averaged computational time for different models. In Table 10, ANSYS-based models cost similar computational time around 500 sec, whereas the NLTc-MATLAB and ECE-MATLAB model take 183 and 338 seconds in average to run the simulation, respectively. Therefore, the benefit of using open-source models, especially the ones with nonlinear truss element, is clearly demonstrated. Such a benefit will be further significant when longer simulation time history is needed, or more transmission towers are included. Comparing between the two transmission lines element models, the elastic catenary element based model consumes additional 90% computational time as compared to the nonlinear truss element based model. For calculating the dynamic response of elastic catenary cable element, as stated in Algorithm 2, the determination of the transmission line stiffness matrix in each sub-iteration at each single time step requires the flexibility matrix to be first computed, and then the stiffness matrix is obtained by the inverse calculation of the flexibility matrix. Consequently, the computational burden increased significantly due to the complexity of inverting the flexibility matrix with the increased number of nodes. By contrast, calculating the transmission line dynamic response with the nonlinear truss element, as introduced in Algorithm 1, the stiffness matrix is directly obtained by implementing the coordinates and strains of the element to the stiffness matrix.Table 10
One-span eight-line tower computational time for difference models at different wind speeds and wind angles.
Computation GPU time (s)NLTc-ANSYSNLTc-MATLABANSYS auto-gravimetricECE-MATLABS10 m/s A00532191514354S10 m/s A30512188494352S10 m/s A90470195481357S20 m/s A00520175518325S20 m/s A30515175493322S20 m/s A90502172491323Averaged508.5182.67498.5338.83S: wind speed; A: wind angle.


## 6. Application of Open-Source Package Full Transmission Tower-Line Model on Extreme Wind Speed Scenario with Aerodynamic Damping Implementation

In Section5, the comparison between the ANSYS models and the open-source MATLAB models showed that the open-source MATLAB model consumes less computational time with the same level of accuracy. Among the open-source MATLAB models, the transmission line with initial strain calibrated nonlinear truss element (NTLc) is computationally superior to the transmission line model with elastic catenary element formulation (ECE). Therefore, to investigate the effects of the aerodynamic damping on the transmission line system, the open-source NLTc-MATLAB model is implemented. Moreover, as mentioned in Section 4.1, the most dominant extreme wind speed in Texas region is 54m/s. Hence, the basic wind speed utilizes 50m/s to generate the dynamic wind speed and wind force for the full transmission tower-line system.Figure19 shows the correlation between the aerodynamic damping ratios and the basic wind speed for the conductors and the ground wires. In Figure 19, the aerodynamic damping ratios are positively correlated with the wind speed. For the first out-of-plane symmetric mode, the aerodynamic damping ratios are not monotonically increased with the wind speed due to the influence of the static position of the transmission line [31].Figure 19
The correlation of aerodynamic damping ratios and wind speeds.The displacement time-history of the transmission tower-line system with and without aerodynamic damping in Figure20 illustrates that the aerodynamic damping has suppressed the vibration of the system. However, in Figure 20, the effect of the aerodynamic damping is not significant to the system. The reason is that, as shown in Figure 19, the aerodynamic damping ratios are small because of the taut transmission line configurations (small sag-to-span ratio) in the model.Figure 20
One-span eight-line model dynamic displacement response comparison with and without aerodynamic damping.
(a)(b)


## 7. Conclusion

In this paper, two different three-dimensional transmission tower-line models are investigated and developed in open-source MATLAB software package and compared with ANSYS commercial software model. The paper compares the performance of transmission tower-line system model with using nonlinear truss element or elastic catenary element for transmission line modeling. From the results of the study, the following conclusions can be drawn and suggestions made:(i)
The initial strain of the nonlinear truss element has a great impact on the shape finding of the transmission line and consequently the small discrepancy of the initial strain will propagate to the transmission line static and dynamic analysis, which will enlarge the errors in the transmission tower-line models. Calibrating the nonlinear truss model initial strain and initial shape with elastic catenary finite element model will improve the dynamic response accuracy.(ii)
The LS-DYNA default nonlinear truss element model uses uncalibrated initial strain; consequently, the LS-DYNA computed responses for both the single line and tower-line coupled system simulation deviate largely from other model outputs.(iii)
Two open-source MATLAB transmission tower-line full models are successfully developed and implemented. The developed models demonstrate accurate representation of transmission tower-line coupling phenomenon.(iv)
The computational time for running a two-tower eight-line model is reduced by 66% with using the open-source model as compare to the commercial software.(v)
Both the elastic catenary finite element and calibrated nonlinear truss element model give accurate static and dynamic responses as the transmission line model candidates. The calibrated nonlinear truss transmission line model computational efficiency is better than the elastic catenary transmission line model.(vi)
A linear lumped mass tower model can reproduce the static and dynamic responses of the transmission tower model in ANSYS in low to middle wind speed, where plastic deformation in the transmission tower does not dominate the dynamic response.(vii)
The open-source transmission tower-line package can integrate the aerodynamic damping into consideration. For the specific transmission tower-line model, the effect of the aerodynamic is not significant due to the taut transmission line configurations.In summary, this paper develops high fidelity three-dimensional transmission tower-line system models in open-source MATLAB package by using lumped tower model with two different transmission line models, namely, the elastic catenary finite element line model and the nonlinear truss calibrated line model. Those models can be implemented to stand-alone transmission line and the coupled transmission tower-line system static and dynamic analysis. However, if the transmission tower elastic-plastic deformation is considered, the lumped mass models in the open-source MATLAB package results will be different from the model in ANSYS because the lumped mass tower model is assumed to be elastic. In the future, transmission tower material nonlinearity should also be considered.



---

*Source: 1959683-2021-12-20.xml*