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

*e 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. Inmost 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. *e 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. *e results indicate that the displacement time-history root-mean-square error (RMSE) of the open-source transmission tower-line model is less than 1% and with a 66% computational time reduction compared with the ANSYS model. *e application of the open-source package transmission tower-line model on extreme wind speed considering the aerodynamic damping is further implemented.


Introduction
Transmission tower-line systems connect power plants to customers and are widely distributed throughout the country. e 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. e 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. e 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 crossbracing system. Xie and Sun [5] further investigated the tower failure mechanism under bending and flexure-rotation loading. e 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 userdefined material model in ABAQUS to simulate the behavior of the tower under various loading conditions. e 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. e transmission line in the transmission tower-line system contains a set of wires connected to the tower, including conductors and ground wires. ese 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. ey 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. ey 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][13][14][15][16][17][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. e 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 trialand-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. e 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. ai 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 highcomputation efficient transmission tower-line model. To achieve such an objective, the performance of an individual transmission line with different elements should be investigated first. en a transmission tower-line model opensource MATLAB software package with the transmission line utilizing elastic catenary cable element and nonlinear truss element is developed and compared. e initial strain condition of nonlinear truss element is modeled in two forms. e 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. e 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.
e rest of the article is organized as follows: in Section 2, 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. e open-source package will be distributed on GitHub at [28] after the article is published.

Methodology
e transmission tower-line system contains two separately designed components: the overhanging wires and the supporting structures. e 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. e 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 towerline structure is fixed at tower bottom position and the transmission lines hanging between adjacent towers.

Transmission-Line
Modeling. Transmission line under gravity forms a natural catenary shape [22]. us, the shape of the transmission line can be found by the elastic catenary element under loading without approximating the gravity. e 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. erefore, the calculation of the initial strain of the nonlinear truss element needs to be addressed. e initial strain of the nonlinear truss element can be calculated as 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 M L 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.

Nonlinear Truss Element.
e tension-only two-node three-dimensional truss element is employed both in ANSYS and MATLAB models. e 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): where dσ * is the stress change, dλ * is the strain change, A * is the cross-sectional area, l 0 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

Shock and Vibration 3
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. e 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. e logarithmic strain and cross section are defined as where l is the deformed length, l 0 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 crosssection area at current state. In equation (5d), the strain is positive definite because the transmission line is tensiononly structure. For dynamic analysis, the equation of motion under external loading can be written as 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 P(t) 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. erefore, to calculate the aerodynamic damping, the method proposed by Wang et al. [31] is adopted.
e 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) e aerodynamic damping matrix can be calculated as By considering the Rayleigh damping, the damping matrix can be represented as 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.

Elastic Catenary Cable Element.
e basic assumptions of the elastic catenary element are as follows: (1) the strain stress relationship obeys Hooke's law, and the crosssectional 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. e coordinates of nodes I and J in Cartesian coordinate are (0, 0, 0) and (l 1 , l 2 , l 3 ). In this figure, the external loads are thermal load ΔT and uniformly distributed loads w 1 , w 2 , and w 3 . In the Lagrangian coordinate, the deformed and undeformed coordinates of the cable are p and s, respectively. e equations of equilibrium are given by where w i is the uniformly distributed load in the direction x i and f i is the tension force in the direction x i . e cable tension at Lagrangian coordinate is e cable tension T is related to the strain by Hooke's law as

Input:
Unstrained element length: L 0 ; Poisson's ratio: ]; Initial element cross section area: A 0 ; Initial element strain: σ 0 ; Elastic modulus: E; Line initial coordinate coord 0 ; Mass matrix: M L ; Initial displacement, velocity, and acceleration matrix u 0 , _ u 0 , and € u 0 ; Damping matrix C; Convergence criteria: eps Output: displacement u i+1 ; velocity _ u i+1 ; acceleration € u i+1 Begin: 1.0 Initial parameters calculation 1.1 State determinations: the force (f s ) 0 and the initial tangent stiffness matrix ( i+1 < eps, skip the following steps and go to step 4.0; otherwise, implement the following steps: 3.8 Update the line tangent stiffness 3.8.1 From coord j+1 i+1 , calculate the element length L j+1 3.8.2 Calculate element length incremental: dL � L j+1 − L 0 3.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 area A * using equation (6) 3.8.7 Calculate the tangent stiffness K LT at iteration j using equation (2) 3.9 Replace j by j + 1 and repeat steps 3.1 to 3.8; after converge, denote final displacement value as u i+1 , and the coordinate of the line as coord i+1 4.0 Calculate the velocity and acceleration for time x 2 Figure 1: Single cable hanging at two nodes with applied loads.

Shock and Vibration
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. e relationship between the Cartesian and Lagrange coordinate is expressed as Substituting equations (9)-(11) into (12), the cable projected length x i can be expressed as a function of the undeformed Lagrangian coordinate s: e boundary conditions of the cable are where l 0 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 f i as Differentiating both sides of equation (15), the relation between projected length and the internal forces is In matrix form, equation (14b) can be expressed as where [F] is the symmetric flexibility matrix with terms given by e stiffness matrix k is the inverse of the flexibility matrix.
e global tangent stiffness and the internal force vector are determined by the six degrees of freedom matrix and vector as e nodal force f 4 , f 5 , and f 6 can be determined through the force equilibrium in x 1 , x 2 , and x 3 direction as 6 Shock and Vibration e stiffness matrix in equation (21) is a function of the internal nodal forces f 1 , f 2 , and f 3 , which is unknown. e internal nodal forces can be initialized based on the equations given by Jayaraman and Knudson [24] and Irvine [26] , which are as follows: in which In other cases, when the initial tension T 0 instead of the unstressed length l 0 is given, the stiffness matrix has four unknowns: f 1 , f 2 , f 3 , and l 0 . e unstressed length is initialized and the forces are modified by Equations (19), (27), and (28) can be solved using Newton-Raphson iteration. e Jacobian of the nonlinear systems of equations is computed from the following equations: e same procedures are employed to solve the nonlinear equation of motion under external loading using the elastic catenary finite element model. Algorithm 2 illustrates the steps to solve equation (7) using elastic catenary cable element.
To obtain the tangent stiffness of the transmission line in Algorithm 2, 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. e 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.

Transmission Tower Modeling.
e 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 timeconsuming. 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. e finite element model and the lumped mass model are shown in Figure 2. e 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 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. e 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 as m � m t /N. After the mass of each node is calculated, the assumed mass matrix is assembled to

Shock and Vibration
Here, K is the stiffness matrix extracted from ANSYS from the abovementioned method; N is the number of nodes; M A is the assumed tower mass matrix; and M T is the calculated tower mass matrix. e 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.

Transmission Tower-Line Model Development.
e 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

Input:
Unstrained element length: L 0 ; Elastic modulus: E; Linear thermal expansion coefficient: α Temperature change: ΔT; Distributed loads in each direction x i : w i , i � 1, 2, 3 Line initial coordinate coord 0 ; Mass matrix M L ; Damping matrix C; Initial stiffness matrix (K LC ) 0 Initial displacement, velocity and acceleration matrix u 0 , _ u 0 , and € u 0 ; Newmark-Beta convergence criteria: eps; Elastic catenary convergence criteria: ϵ Output: displacement u i+1 ; velocity _ u i+1 ; acceleration € u i+1 Begin: 1.0 Initial parameters calculation 1.1 State determinations: the force (f s ) 0 and the initial tangent stiffness matrix 3 Select Newmark-Beta parameters c and β, and time interval dt i+1 < eps, skip the following steps and go to step 4.0; otherwise, implement the following steps: Initialize the convergent parameter dl 3.8.7 While dl > ϵ 3.8.8 Calculate parameter c 1 using (16) 3.8.9 Calculate the force in x i at J using (24c) 3.8. 10 Calculate the end force at node I and J using (10) 3.8.11 Calculate the projection length L xi of the element at Cartesian coordinates using (15) 3.8.12 Calculate dl � L xi − L 0xi 3.8. 13 Calculate the flexibility matrix F e using (19) 3.8.14 Update the internal force incremental df � F e /dl 3.8. 15 Update the internal force vector at node I( Calculate the stiffness matrix k e using (21) 3.8.17 Assemble the element stiffness k E using (22) 3.8. 18 Assemble the global stiffness matrix K LC 3.9 Replace j by j + 1 and repeat steps 3.1 to 3.8; after converge, denote final displacement value as u i+1 , and the coordinate of the line as coord i+1 4.0 Calculate the velocity and acceleration for time 8 Shock and Vibration considered as pinned connection [18]. At the tower and line coupled nodes, the stiffness and mass in each direction is directly added as K C and M C . With the transmission towerline 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.
In Figure 3, 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.

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. Figure 4 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.  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.

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. e 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.

Initial Shape and Response
Assessment. e 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.

Shape of Cable Structure under Gravity and Static Wind Loading.
e initial shape of the ground wire and conductor under gravity loading is shown in Figures 5 and 6. In Table 2, the sag difference between SAP2000 and elastic catenary finite element model is 0.001%. e 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. e 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. en, 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. e 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: where v is the wind speed at height z; Ψ is the wind incident angle; c w is the load factor; Q is the numerical constant; K z is the velocity pressure exposure coefficient; K zt is the topographic factor; G is the tower gust response factor; A mt and A ml is the projection area in transverse and longitudinal direction; and C ft and C fl is the force coefficient in transverse and longitudinal direction. e wind speed at height z is governed by the power law as where v s is the basic wind speed at standard height z s and α is the power law exponent that represents the surface roughness. e basic wind speed defined in ASCE 7-10 [35] is a three-second gust speed at 10 m above the ground in Exposure C. us, the standard height z s is 10 m. 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], 54 m/s is the most dominant basic wind speed in Texas region. Accordingly, the basic wind speed for the transmission line static analysis takes 54 m/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 20 m 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. Figures 8 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 autogravimetric method.

Dynamic Response of Cable Structure with Earthquake and Dynamic Wind Loading.
e 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-meansquared-error (RMSE) is used as the error indicator.      (1) Earthquake loading: the El-Centro ground motion acceleration data is used as the benchmark to evaluate the performance of the models. e earthquake loading is applied in the transverse direction of the transmission line as the transmission line has smaller stiffness in that direction. Figure 10 shows the transmission line middle position displacement time-history in y and z directions. e 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.
ose large RMSEs demonstrate that the LS-DYNA model does not obtain the correct transmission line displacement time-history. e 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.
(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). e 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   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. H jm is the S(ω) Cholesky decomposed matrix as in 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 Figure 11, the nonlinear truss calibrated model and the elastic catenary model displacement time-history is overlapped, and the ANSYS auto-gravimetric method displacement timehistory 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.
In summary, from the dynamic wind loading cases, transmission line displacement response magnitude and trend in LS-DYNA is far off from the rest models. erefore, 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. e 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. e 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. e initial shape of the transmission line determines the coordinates of the nonlinear truss element. In Algorithm 1, 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.

Development of Transmission Tower-Line System Model
In Section 3, the single cable lumped mass models have been

Single Tower Lumped Mass Tower
Model. e transmission tower studied is a steel-made, L-shaped cross-sectional suspension tower with a height of 31.5 m. e 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. e 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. e tower stiffness and mass matrix are calculated based on the method in Section 2.2. e static wind loading calculation is based on Section 3.2.1. In the extreme wind map from ASCE Manual 74, 54 m/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 54 m/s. e wind angle is   chosen as 30°to verify the lumped mass tower model x and y direction stiffness formulation. e lumped mass transmission tower static displacement under static wind loading displacement difference with respect to that from the ANSYS detailed model is shown in Figure 12. 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. e dynamic wind force generation method is described in Section 2.3. e 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] In Figure 13, 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%. e 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.
To conclude, from the single tower static and dynamic analysis, the lumped mass can be treated as an alternative model to the ANSYS detailed model. e LS-DYNA model has relatively large displacement difference to the ANSYS model.

Transmission Line Element Size Reduction.
e span of the transmission line is 200m; 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. e 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.

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. e two transmission lines are ground wires, whose material properties are shown in Table 1. e span between two towers is 200 m. 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 10 m/s and wind angle 30°are chosen. e 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 16 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.
e 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.
In Tables 6 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%. e 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 Figure 16, the LS-DYNA models displacement time-history is different from other models. In Tables 6 and 7, the LS-DYNA maximum displacement timehistory 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.  e 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. e time interval for the lumped mass model is selected as 0.0025 s for optimal computational efficiency and accuracy trade-off.

Performance of Full Transmission Tower-Line System Model during Wind Loadings
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. Figure 18 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.
In Table 8, 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%. e 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. Table 9 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 autogravimetric model, at the same wind speed 10 m/s, the largest RMSE is 0.1% at wind angle 90 and at wind speed 20 m/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 autogravimetric model and the NLTc-ANSYS model is in the transmission line initial strain and initial shape determination. e small initial strain and initial shape difference     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. Table 10 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. erefore, 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.

Application of Open-Source Package Full Transmission Tower-Line Model on Extreme Wind Speed Scenario with Aerodynamic Damping Implementation
In Moreover, as mentioned in Section 4.1, the most dominant extreme wind speed in Texas region is 54 m/s. Hence, the basic wind speed utilizes 50 m/s to generate the dynamic wind speed and wind force for the full transmission tower-line system. Figure 19 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 outof-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]. e displacement time-history of the transmission tower-line system with and without aerodynamic damping in Figure 20 illustrates that the aerodynamic damping has suppressed the vibration of the system. However, in Figure 20, the effect of the aerodynamic damping is not

24
Shock and Vibration significant to the system. e reason is that, as shown in Figure 19, the aerodynamic damping ratios are small because of the taut transmission line configurations (small sag-tospan ratio) in the model.

Conclusion
In this paper, two different three-dimensional transmission tower-line models are investigated and developed in opensource MATLAB software package and compared with ANSYS commercial software model. e 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) e 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) e LS-DYNA default nonlinear truss element model uses uncalibrated initial strain; consequently, the LS- In summary, this paper develops high fidelity threedimensional 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. ose 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.