High Voltage Transmission Line Vibration: Using MATLAB to Implement the Finite Element Model of a Wind-Induced Power-Line Conductor Vibration

Abstract

Wind-induced vibration affects the performance and structural integrity of high voltage transmission lines. The finite element method (FEM) is employed to investigate wind-induced vibration in MATLAB. First, the FEM model was used to develop the equation of motion of the power line conductor. In addition, dampers, conditions for damping, free and forced vibrations of the overhead conductor were considered in the FEM model. Wind-induced experiments were conducted in the laboratory using an actual overhead power conductor. The developed FEM models were simulated in the MATLAB computing environment. The results from the MATLAB simulation, finite element and experimental recordings were compared in order to evaluate the efficacy of models simulated in MATLAB and developed using the FEM.

1. Introduction

The availability and use of electrical power in the society is crucial in the development and growth of the society. Power generated from power stations is transmitted using high voltage transmission lines. The transmission of power from the point of generation to the point of use requires complex network of high voltage lines, systems and components [1]. High voltage conductors are usually subjected to vibration and the vulnerabilities of the power lines to vibration can lead to fatigue failure. Thus, power loading determination and control on the power grid can influence the integrity of the transmission network. High voltage conductor vibration is very difficult to model due to the fact that the responses exhibit a non-linear behaviour. There has been concerted effort to try and predict the conductor response as a result of aeolian vibration. Evaluations of conductor vibration caused by aeolian forces have been investigated be several researchers [2, 3, 4, 5, 6].

Recent researches developed models to investigate wind-induced vibration using nonlinear time history, expert systems, the concept of principal modes, aero-elastic and bending stiffness [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. The investigation of vortex formation and the phenomenon of wind-induced vibration was done using the concept fluid–solid dynamic excitation [18, 19, 20]. The models were used to determine how wind loading influenced the oscillation of transmission lines. This form of investigation was done by experimental studies carried out in a wind tunnel [21, 22]. The outcomes of these experiments were used to determine conductor input loading. Several models developed by various researchers can be used to determine conductor damping and also the placement of vibration absorbers on the line conductors to curtail the effect of cable mechanical oscillation [23, 24, 25, 26, 27]. Based on the various models that have been developed by researchers as indicated in the first and second paragraphs, there is a need to further analyze wind-induced vibration using finite element method (FEM) in MATLAB.

The design, construction and maintenance of power transmission network requires adequate understanding of the system dynamics that occurs when subjected to vortex induced vibration [28]. Various analysis can be conducted using techniques that suit certain objectives. System integrity in high voltage transmission lines is of paramount importance. MATLAB is a multi-model simulation environment used for numerical computing. It can be integrated with physical hardware or systems in order to determine real-time performance, characteristics and behavior. MATLAB also provides a platform for special hardware in loop simulations [29]. These functionalities amongst others are vital in determining various characteristics and behaviors in high voltage transmission lines.

High voltage transmission lines and grid can experience vulnerabilities such as vibration, electromagnetic transients, fatigue, transmission loss, switching surges, conductor sag fluctuation [30, 31]. When the conductor experiences vibration, the transmission lines experience high amplitudes of vibrations from wind forces and can lead to fatigue of the transmission lines [1]. The use of systems simulation and analysis provides the platform to understand the response of the transmission conductor. The responses considered in the chapter include transmission line excitation through wind loading, conductor properties such as damping and damper placement used in mitigating the vibration.

The chapter discussed the development and implementation of a wind-induced high voltage transmission line vibration using finite element method (FEM) in MATLAB. The sections in the chapter discussed the development of transmission line equation of motion, the solution to the equation of motion, free and forced vibration of the transmission line, dampers and conductor self-damping, FEM MATLAB setup and implementation, simulation of FEM models. The chapter also discussed results from FEM models, simulation and experimental investigation. The chapter is focused towards the development of a finite element method and its implementation on the MATLAB software. The developed finite element method (FEM) was based on the concept of the simply supported beam model and it was used in modeling the transverse vibration of power line conductors. The results from the FEM were then compared with results from the analytical model and results obtained from experimental studies documented in [1]. The results from MATLAB simulations from the finite element models and experimental results were compared in order to determine the accuracy of the models. The developed FEM was then used as the means to verify the effect of varying the conductor axial tension on the natural frequencies of the conductors.

2. Transmission line equation of motion (EOM)

The transverse displacement of high voltage transmission line conductor is generally caused by wind loading. This form of vibration with small displacement is known as aeolian vibration and it is a source of concern to the power lines reliability. One the vulnerabilities is that it can cause fatigue failure of the transmission lines. Conductors are example of continuous or distributed systems and modeling its mechanical vibration can either be as a beam or taut string. In [18, 19], it was ascertained that modeling a conductor as a beam is more accurate than modeling it as a taut string due to the effect of the bending stiffness. Hence, in line with the above, the conductor transverse vibration was modeled as a beam, simply supported or pinned at both ends. The distributed loading on the conductor is replaced by effective point load that can effectively have the same resultant effect as that of the actual distributed load.

The high voltage transmission line equation of motion was formulated by assuming that power conductors can modeled as beams with fixed ends. The following assumptions were considered [1]:

• The power conductor is uniform along its length and it is slender

• The power conductor is a solid with cylindrical body having both linear and homogeneous physical properties throughout its cross-sectional area

• The power conductor has a symmetrical plane which acts as the plane of vibration such that there is the decoupling of translational and rotational motion.

The assumptions were based on beam theory. In considering the power conductor as a beam, sagged by a tensile force $S$ , being acted upon by a concentrated wind load $f\left(x,t\right)$ , with cross-sectional area $A$ , density $\rho$ , flexural rigidity $\mathit{EI}$ , displaced at a distance of $x$ after time. In Eq. (1), the high voltage transmission line equation of motion is expressed as:

$$f\left(x,t\right)=\mathit{EI}\frac{{\partial }^{4}y\left(x,t\right)}{\partial x^4}-S\frac{{\partial }^{2}y\left(x,t\right)}{\partial x^2}+\mathit{\rho A}\frac{{\partial }^{2}y\left(x,t\right)}{\partial t^2}$$ E1

For $x\in \left(0,1\right),t>0$ . The boundary conditions are expressed and indicated in Eqs. (2) and (3):

$$y\left(0,t\right)=\frac{{\partial }^2\left(0,t\right)}{\partial x^2}=0$$ E2

$$y\left(l,t\right)=\frac{{\partial }^2\left(l,t\right)}{\partial x^2}=0$$ E3

The initial conditions at $t=0$ are indicated in Eq. (4), Eq. (5) and expressed as:

$$y\left(x,0\right)=y_o\left(x\right)$$ E4

$$\dot{y}\left(x,0\right)=\dot{y}\left(x\right)$$ E5

Introducing the mass per unit length of the power conductor, the new equation of motion indicated in Eq. (6) is expressed as:

$$f\left(x,t\right)=\mathit{EI}\frac{{\partial }^{4}y\left(x,t\right)}{\partial x^4}-S\frac{{\partial }^{2}y\left(x,t\right)}{\partial x^2}+m\frac{{\partial }^{2}y\left(x,t\right)}{\partial t^2}$$ E6

In order to derive a possible solution, the model was simplified using dimensionless functions and Dirac delta functions. In Eqs. (7)–(12), the variables are expressed in dimensionless form and expressed as:

$$Y=\frac{y\left(x,t\right)}{D}$$ E7

$$X=\frac{x}{L}$$ E8

$$\tau =\frac{t}{f}$$ E9

$$I_p=\frac{{\mathit{Df}}^2}{g}$$ E10

$$S_p=\frac{SD}{\gamma L^2}$$ E11

$$M_p\frac{\mathit{EI}D}{\gamma L^4}$$ E12

Eq. (13) indicates the revised equation of motion and it is expressed as:

$$M_p.\frac{{\partial }^{4}Y}{\partial X^4}-S_p\frac{{\partial }^{2}Y}{\partial X^2}+I_p\frac{{\partial }^{2}Y}{\partial {\tau }^2}=\frac{1}{\gamma }\left[F\left(X,\tau \right)+\sum _n\delta \left(X-X_n\right)F_n\left(\tau \right)\right]$$ E13

Where $\gamma$ represents the power conductor weight per unit length and $g$ represents gravitational constant. $X_n\delta \left(X-X_n\right)$ represents the Dirac delta function, $F\left(X,\tau \right)$ denotes the net transverse force per unit length acting on the conductor and $F_n\left(\tau \right)$ denotes the nth concentrated force acting transversely on the conductor.

3. Solution to the EOM

The general solution to the high voltage transmission line equation of motion was derived using Euler-Bernoulli equation. The particular solution to the equation of motion was derived using a product of two functions. The two functions were first separated using the principle of variable separation as expressed in Eq. (14) [19]:

$$Y\left(x,t\right)=X\left(x\right)T\left(t\right)$$ E14

Where $X\left(x\right)$ is the normalized function representing the mode shape of the equation of motion. The normalized function ensures that orthogonality condition was satisfied in the derivation of the EOM model solution. Applying the normalized function in the EOM yields Eqs. (15) and (16):

$$\mathit{EI}\stackrel{////}{X}\left(x\right)-S\stackrel{//}{X}\left(x\right)-{\omega }^2\mathit{\rho AX}\left(x\right)=0$$ E15

$$\ddot{T}\left(t\right)+{\omega }^{2}T\left(t\right)=0$$ E16

Where $\stackrel{////}{X}\left(x\right)=\frac{d^{4}y}{{\mathit{dx}}^4}$ , $\stackrel{//}{X}\left(x\right)=\frac{d^{2}y}{{\mathit{dx}}^2}$ , $\ddot{T}\left(t\right)=\frac{d^{2}y}{{\mathit{dt}}^2}$ and ${\omega }^2$ is a constant that equates $x$ and $t$ . Assuming that $X\left(x\right)={\mathit{Ze}}^{\mathrm{\Psi }x}$ , the model is expressed in Eq. (17) as:

$${\mathit{Ze}}^{\mathrm{\Psi }x}\left(\mathit{EI}{\mathrm{\Psi }}^4-S{\mathrm{\Psi }}^2-\mathit{\rho A}{\omega }^2\right)=0$$ E17

Considering that ${\mathit{Ze}}^{\mathrm{\Psi }x}\ne 0$ , hence $\left(\mathit{EI}{\mathrm{\Psi }}^4-S{\mathrm{\Psi }}^2+\mathit{\rho A}{\omega }^2\right)=0$ . The general solution of the Euler-Bernoulli equation which represents the solution to the equation of the motion of the transmission line is expressed in Eqs. (18) and (19) as [9]:

$${\mathrm{\Omega }}^2,{\mathrm{\Psi }}^2=-\frac{\left(-S\right)\pm \sqrt{S^2-4\left(\mathit{EI}\right)\left(-\mathit{\rho A}{\omega }^2\right)}}{2\mathit{EI}}$$ E18

$$\mathrm{\Omega },\mathrm{\Psi }=\left(\pm \right)\sqrt{\frac{S\pm \sqrt{S^2+4\mathit{EI}(\mathit{\rho A}}{\omega }^2)}{2\mathit{EI}}}$$ E19

The values of $\mathrm{\Omega }$ and $\mathrm{\Psi }$ represents the general solution of the equation of motion. The practical implication of the derived solution is that it represents the transverse vibration of the high voltage transmission line. The derived solution has infinite number of solutions and the solution is indexed to accommodate all the possible solutions from the model. The indexed solution is expressed in Eqs. (20) and (21) as:

$${\mathrm{\Omega }}_n=\sqrt{\frac{S}{2\mathit{EI}}+\sqrt{\frac{S^2}{{\left(2\mathit{EI}\right)}^2}+m_L\frac{{\left(2\pi f_n\right)}^2}{\mathit{EI}}}}$$ E20

$${\mathrm{\Psi }}_n=\sqrt{-\frac{S}{2\mathit{EI}}+\sqrt{\frac{S^2}{{\left(2\mathit{EI}\right)}^2}+m_L\frac{{\left(2\pi f_n\right)}^2}{\mathit{EI}}}}$$ E21

Where ${\omega }_n=2\pi f_n$ and for $n=1,2,3,\dots$ .

In Eqs. (22)–(24), the infinite natural frequencies of the power conductor were derived while considering that the mode shape is the same as a pinned-pinned beam eigenfunction model with no external force. Hence,

$$Y_n\left(x,t\right)=\text{sin}\frac{\mathit{n\pi x}}{l}\text{cos}{\omega }_{n}t$$ E22

$$\mathit{EI}{\left(\frac{\mathit{n\pi }}{l}\right)}^4\text{sin}\frac{\mathit{n\pi x}}{l}\text{cos}{\omega }_{n}t-S{\left(-\frac{\mathit{n\pi }}{l}\right)}^2\text{sin}\frac{\mathit{n\pi x}}{l}\text{cos}{\omega }_{n}t+\mathit{\rho A}\left(-{\omega }_n\right)\text{sin}\frac{\mathit{n\pi x}}{l}\text{cos}{\omega }_{n}t=0$$ E23

$$\text{sin}\frac{\mathit{n\pi x}}{l}\text{cos}{\omega }_{n}t\left[\frac{\mathit{EI}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4+\frac{S}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^2-{{\omega }_n}^2\right]=0$$ E24

The natural frequency of the power conductor in rad/s is expressed in Eqs. (25) and (26) as:

$${\omega }_n^2={\left(\frac{\mathit{n\pi }}{l}\right)}^2\frac{S}{\mathit{A\rho }}+{\left(\frac{\mathit{n\pi }}{l}\right)}^4\frac{\mathit{EI}}{\mathit{A\rho }}$$ E25

$${\omega }_n=\sqrt{{\left(\frac{\mathit{n\pi }}{L}\right)}^2\frac{S}{m_L}\left[1+{\left(\frac{\mathit{n\pi }}{L}\right)}^2\frac{\mathit{EI}}{S}\right]}$$ E26

The natural frequency in Hz is expressed in Eq. (27) as:

$$F_n=\frac{1}{2\pi }\sqrt{{\left(\frac{\mathit{n\pi }}{L}\right)}^2\frac{S}{m_L}\left[1+{\left(\frac{\mathit{n\pi }}{L}\right)}^2\frac{\mathit{EI}}{S}\right]}$$ E27

4. Free vibration of power conductor

The self-damping model of the power conductor provided the basis to analyze free vibration experienced by the conductor. Free vibration occurs when the forcing function causing the power conductor to vibrate become zero. Hence the equation of motion is expressed in Eq. (28) as [19]:

$$\mathit{EI}\frac{{\partial }^{4}y\left(x,t\right)}{\partial x^4}-S\frac{{\partial }^{2}y\left(x,t\right)}{\partial x^2}+\mathit{\beta I}\frac{{\partial }^{5}y\left(x,t\right)}{\partial x^4\partial t}+C\frac{\partial y\left(x,t\right)}{\partial t}+\mathit{\rho A}\frac{{\partial }^{2}y\left(x,t\right)}{\partial t^2}=0$$ E28

Applying the principle of separation of variable to the equation of motion yields Eq. (29):

$$\mathit{EI}\stackrel{////}{X}\left(x\right)T\left(t\right)-S\stackrel{//}{X}\left(x\right)T\left(t\right)+\mathit{\beta I}\stackrel{///}{X}\left(x\right)\dot{T}\left(t\right)+\mathit{CX}\left(x\right)\dot{T}\left(t\right)+\mathit{\rho AT}\left(t\right)X\left(x\right)=0$$ E29

Integrating the eigenfunction $X_n\left(x\right)=\text{sin}\left(\frac{\mathit{n\pi x}}{l}\right)$ in the model yields Eqs. (30)–(33):

$$\begin{array}{l}\mathit{EI}\left[{\left(\frac{\mathit{n\pi }}{l}\right)}^4\text{sin}\left(\frac{\mathit{n\pi x}}{l}\right)\right]T\left(t\right)-S\left[\left(\frac{-\mathit{n\pi }}{l}\right)\text{sin}\left(\frac{\mathit{n\pi x}}{l}\right)\right]T\left(t\right)\\ +\mathit{\beta I}\left[{\left(\frac{\mathit{n\pi }}{l}\right)}^4\text{sin}\left(\frac{\mathit{n\pi x}}{l}\right)\right]\dot{T}\left(t\right)+C\left[\text{sin}\left(\frac{\mathit{n\pi x}}{l}\right)\right]\dot{T}\left(t\right)+\mathit{\rho A}\left[\text{sin}\left(\frac{\mathit{n\pi x}}{l}\right)\right]\ddot{T}\left(t\right)=0\end{array}$$ E30

$$\text{sin}\left(\frac{\mathit{n\pi x}}{l}\right)\left[\begin{array}{l}\mathit{EI}{\left(\frac{\mathit{n\pi }}{l}\right)}^{4}T\left(t\right)+S{\left(\frac{\mathit{n\pi }}{l}\right)}^2\dot{T}\left(t\right)\\ +\mathit{\beta I}{\left(\frac{\mathit{n\pi }}{l}\right)}^{4}T\left(t\right)+C\dot{T}\left(t\right)+\mathit{\rho A}\ddot{T}\left(t\right)\end{array}\right]=0$$ E31

$$\mathit{\rho A}\ddot{T}\left(t\right)+\left[\mathit{\beta I}{\left(\frac{\mathit{n\pi }}{l}\right)}^4+C\right]\dot{T}\left(t\right)+\left[S{\left(\frac{\mathit{n\pi }}{l}\right)}^2+\mathit{EI}\left(\frac{\mathit{n\pi }}{l}\right)\right]T\left(t\right)=0$$ E32

$$\ddot{T}\left(t\right)+\left[\frac{\mathit{\beta I}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4+\frac{C}{\mathit{\rho A}}\right]\dot{T}\left(t\right)+\left[\frac{S}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^2+\frac{\mathit{EI}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4\right]T\left(t\right)=0$$ E33

Considering that the vibration model represents a multi-degree vibration system. The natural frequency of the power conductor is determined and expressed in Eqs. (34) and (35):

$${\omega }_n^2=\frac{S}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^2+\frac{\mathit{EI}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4$$ E34

$$2{{\mathit{\xi \omega }}_n}^2=\left[\frac{\mathit{\beta I}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4+\frac{C}{\mathit{\rho A}}\right]$$ E35

The temporal solution to the free vibration model is expressed in Eq. (36) as:

$$T_n=A_1{e}^{-{\xi }_n{\omega }_{n}t}\text{sin}\left({\omega }_{d}t+\varphi \right)$$ E36

The solution can also be represented in Eq. (37) and expressed as:

$$T_n=e^{-{\xi }_n{\omega }_{n}t}\left(B_1\text{sin}{\omega }_{d}t+B_2\text{cos}{\omega }_{d}t\right)$$ E37

Where the damped frequency of the power conductor is expressed in Eq. (38) as:

$${\omega }_d={\omega }_n\sqrt{1-{\xi }^2}$$ E38

The system response is expressed in Eq. (39) as:

$$y\left(x,t\right)=\sum _{n=1}^{\infty }{A}_1{e}^{-{\xi }_n{\omega }_{n}t}\text{sin}\left({\omega }_{d}t+\phi \right)\text{sin}\frac{\mathit{n\pi x}}{l}$$ E39

The response can also be represented in Eq. (40) and expressed as:

$$y\left(x,t\right)=\sum _{n=1}^{\infty }\left[e^{-{\xi }_n{\omega }_{n}t}\left(B_1\text{sin}{\omega }_{d}t+B_2\text{cos}{\omega }_d\right)\right]\text{sin}\frac{\mathit{n\pi x}}{l}$$ E40

5. Forced vibration of power conductor

High voltage transmission lines are exposed to loading from the wind. The actual system representation through system simulation strategy considers a case of distributed load through the span of the conductor. In order to simplify simulations, the external force acting on the conductor is represented as a point load. In Eqs. (41)–(43), the equation of motion is solved with an excitation force in order to evaluate the actual response of high voltage transmission lines under aeolian vibration [1]. Hence,

$$\mathit{EI}\frac{{\partial }^{4}y\left(x,t\right)}{\partial x^4}-S\frac{{\partial }^{2}y\left(x,t\right)}{\partial x^2}+\mathit{\beta I}\frac{{\partial }^{5}y\left(x,t\right)}{\partial x^4\partial t}+C\frac{\partial y\left(x,t\right)}{\partial t}+\mathit{\rho A}\frac{{\partial }^{2}y\left(x,t\right)}{\partial t^2}=f\left(x,t\right)$$ E41

$$\text{sin}\left(\frac{\mathit{n\pi x}}{l}\right)\left[\begin{array}{l}\mathit{EI}{\left(\frac{\mathit{n\pi }}{l}\right)}^{4}T\left(t\right)+S{\left(\frac{\mathit{n\pi }}{l}\right)}^2\dot{T}\left(t\right)\\ +\mathit{\beta I}{\left(\frac{\mathit{n\pi }}{l}\right)}^{4}T\left(t\right)+C\dot{T}\left(t\right)+\mathit{\rho A}\ddot{T}\left(t\right)\end{array}\right]=F\text{sin}{\omega }_{\mathit{dr}}t$$ E42

$$\ddot{T}\left(t\right)+\left[\frac{\mathit{\beta I}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4+\frac{C}{\mathit{\rho A}}\right]\dot{T}\left(t\right)+\left[\frac{S}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^2+\frac{\mathit{EI}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4\right]T\left(t\right)=F\text{sin}{\omega }_{\mathit{dr}}t$$ E43

Expressing the model as a multi-degree system yields Eq. (44):

$$T\left(t\right)={\mathit{Ae}}^{-{\mathit{\zeta \omega }}_{n}t}\text{sin}\left({\omega }_{d}t+\varphi \right)+X\text{cos}\left(\mathit{\omega t}-\theta \right)$$ E44

The natural frequency of the power conductor under forced vibration is expressed in Eqs. (45) and (46) as:

$${\omega }_n^2=\frac{S}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^2+\frac{\mathit{EI}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4$$ E45

$$2{{\mathit{\xi \omega }}_n}^2=\left[\frac{\mathit{\beta I}}{\mathit{\rho A}}{\left(\frac{\mathit{n\pi }}{l}\right)}^4+\frac{C}{\mathit{\rho A}}\right]$$ E46

The solution to the equation of motion under forced vibration is expressed in Eq. (47) as:

$$y\left(x,t\right)=\text{sin}\frac{\mathit{n\pi x}}{l}\left[{\mathit{Ae}}^{-\mathit{\zeta \omega t}}\text{sin}\left(\text{sin}\mathit{\omega t}+\varphi \right)+X\text{cos}\left(\mathit{\omega t}-\theta \right)\right]$$ E47

6. Conductor self-damping and dampers

The influence of external and internal damping mechanisms was considered in the conductor vibration model. The factors considered included the following [1, 32]:

• The power conductor inter-strand motion and fluid damping. This is proportional to the conductor velocity and represented as viscous damping in the model.

• The rate of strain in the power conductor. This proportional to the internal damping of the power conductor.

The high voltage transmission line damped model is expressed in Eq. (48) as:

$$\mathit{EI}\frac{{\partial }^{4}y\left(x,t\right)}{\partial x^4}-S\frac{{\partial }^{2}y\left(x,t\right)}{\partial x^2}+\mathit{\beta I}\frac{{\partial }^{5}y\left(x,t\right)}{\partial x^4\partial t}+C\frac{\partial y\left(x,t\right)}{\partial t}+\mathit{\rho A}\frac{{\partial }^{2}y\left(x,t\right)}{\partial t^2}=f\left(x,t\right)$$ E48

Where $C$ and $\beta$ represent damping constants. In the presence of axial load, viscous air damping, strain rate damping or Kelvin-Voigt damping, high voltage transmission line integrity can be managed.

There are various types of dampers that can be used to reduce vibration. The dampers are excited by the vibration of the power conductor and the vibration of their masses connected by the massager cable help to damp out energy. Stockbridge dampers are commonly installed on high voltage transmission lines to reduce aeolian vibrations. Stockbridge dampers can be symmetrical or asymmetrical in their design. An example of dampers installed on high voltage transmission lines is shown in Figure 1 . The design of Stockbridge dampers follows the principle of cantilever beams with mass at the free ends. The contribution of dampers to power conductor vibration mitigation is to lower the severity of the vibration to a level that might prevent failure to the line.

Figure 1.

Asymmetrical damper.

7. FEM MATLAB model setup, formulation and implementation

In order to implement the conductor model in MATLAB environment, finite element analysis formulation was done as function of the physical state of power transmission line conductor. The models developed using finite element analysis can then be implemented in MATLAB. The FEM model enables the analysis of the dynamic behavior and response of power line conductor to the dynamic forces of wind [33]. Consider a power transmission line subjected to dynamic aeolian vibration as an assembly of thin strands having distributed mass and elasticity. The physical model can be represented used partial differential equations. Each strand in the transmission line experiences axial, bending and torsional loads from the wind [34, 35]. The accurate representation of each factor is critical in the determination of the dynamic behavior of power transmission lines [36, 37]. Euler-Bernoulli curved beam theory was used to formulate the finite element model of power transmission lines.

Consider a power transmission line experiencing a vertical force, curvature and an axial force has an axial displacement modeled in Eq. (49) as [21]:

$$u\left(x,y\right)=u_0\left(x\right)+\frac{v}{R}-y\left({\theta }_x+\frac{u\left(x\right)}{R}\right)$$ E49

Where ${\theta }_x$ represents the rotation of the power line due to flexural effect, $R$ represents the radius of rotation, $y$ represents the distance from the axis of rotation to the centroidal axis of the conductor or transverse displacement, $v$ represents tangential displacement and $u\left(x\right)$ represents the axial displacement of the power lines. The shape function for power transmission lines having rotation, bending and axial motion components is modeled using discretization techniques and represented in Eqs. (50)–(52) as:

$$u\left(s\right)=b_0+b_{1}s$$ E50

$$v\left(s\right)=c_0+c_{1}S+c_2{S}^2+c_3{S}^3$$ E51

$$\theta =\frac{dv\left(S\right)}{ds}=c_1+2c_{2}S+3c_3{S}^2$$ E52

The solution to the discretization of the models yields Eq. (53) to Eq. (59):

$$\left[\begin{array}{l}u\\ v\\ \theta \end{array}\right]=\left[\begin{array}{lll}{N}_1& 0& N_2\\ 0& N_3& 0\\ 0& 0& N_4\end{array}\begin{array}{lll}0& 0& 0\\ N_4& N_5& N_6\\ 0& N_5& N_6\end{array}\right]\left[\begin{array}{l}{u}_1\\ \begin{array}{l}\begin{array}{l}{v}_1\\ {\theta }_1\end{array}\\ u_2\\ v_2\end{array}\\ {\theta }_2\end{array}\right]$$ E53

Where

$$N_1=\frac{1}{2}\left(1-\zeta \right)$$ E54

$$N_2=\frac{1}{2}\left(1+\zeta \right)$$ E55

$$N_3=\frac{1}{4}\left(2-3\zeta +{\zeta }^3\right)$$ E56

$$N_4=\frac{1}{4}\left(1-\zeta -{\zeta }^2-{\zeta }^3\right)$$ E57

$$N_5=\frac{1}{4}\left(2+3\zeta -{\zeta }^3\right)$$ E58

$$N_6=\frac{1}{4}\left(1-\zeta +{\zeta }^2+{\zeta }^3\right)$$ E59

The power line matrix model contains the strand stiffness $K$ , mass matrix $M$ and the load vector $F$ . They are expressed in Eq. (60) as:

$$\left[K\right]=\frac{1}{2}\int {N_a}^T\left(\mathit{EA}\right)N_a\mathit{\delta \zeta }+\frac{1}{2}\int {N_B}^T\left(\mathit{EI}\right)N_B\mathit{\delta \zeta }+\frac{1}{2}\int {N_B}^T\left(T\right)N_{B,T}\mathit{\delta \zeta }$$ E60

Where $A$ represents the cross-sectional area of the power line, $E$ represents the young modulus of the power line material, $I$ represents polar moment of area, $T$ represents the kinetic energy of the system. The matrix is modeled in Eq. (61).

$$\left[M\right]=\frac{1}{2}\int {\dot{u}}^T\mathit{\rho A}\dot{u}+\frac{1}{2}\int {\dot{g}}^T\mathit{\rho A}\dot{v}$$ E61

Where $\rho$ represents the density of the power line material and Eq. (62) indicates external excitation.

$$\mathit{\delta W}=\frac{1}{2}\int F\mathit{\delta u}$$ E62

The power line conductor model is constructed using Euler-Bernoulli theories and summarized in Eq. (63) as:

$$\left[\begin{array}{ll}{M}_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right]\left(\begin{array}{l}\ddot{u}\\ \ddot{v}\end{array}\right)+\left[\begin{array}{ll}{K}_{11}& K_{12}\\ K_{21}& K_{22}\end{array}\right]\left(\begin{array}{l}u\\ v\end{array}\right)=\left(\begin{array}{l}{F}_1\\ F_2\end{array}\right)$$ E63

The finite element analysis follows a step by step numerical computation in the MATLAB environment as documented in [38, 39]. The dynamic response analysis assumes continuous displacement, velocity and acceleration [40, 41]. The numerical integration technique utilized was based on Newmark integration method. The compact form of the high voltage transmission line model is expressed in Eqs. (64)–(69) as [18]:

$$\left[M\right]\left\{\ddot{y}\right\}+\left[C\right]\left\{\dot{y}\right\}+\left[K\right]\left\{y\right\}=\left[F\right]$$ E64

$${\left[\stackrel{\wedge }{K}\right]}_{s+1}={\left[K\right]}_{s+1}+a_3{\left[m\right]}_{s+1}$$ E65

$${\left[\stackrel{\wedge }{F}\right]}_{s,s+1}={\left\{F\right\}}_{s+1}+{\left[m\right]}_{s+1}\left(a{\left\{y\right\}}_s+a{\left\{\dot{y}\right\}}_s+a{\left\{\ddot{y}\right\}}_s\right)$$ E66

$$a_3=\frac{2}{\gamma {\left(\mathrm{\Delta }t\right)}^2}$$ E67

$$a_4=\frac{2}{\gamma \mathrm{\Delta }t}$$ E68

$$a_5=\frac{1}{\gamma }-1$$ E69

The initial conditions are expressed in Eq. (70) as:

$${\left[\ddot{y}\right]}_0={\left[M\right]}^{-1}{\left[F\right]}_0-{\left[K\right]}^{-1}{\left[y\right]}_0$$ E70

The acceleration vector is expressed in Eqs. (71)–(74) as:

$${\left[\ddot{y}\right]}_{s+1}=a_3\left({\left\{y\right\}}_{s+1}-{\left\{y\right\}}_s\right)-a_4\left\{\dot{y}\right\}-a_5{\left\{\ddot{y}\right\}}_s$$ E71

$${\left[\ddot{y}\right]}_{s+1}={\left\{\ddot{y}\right\}}_s+a_2{\left\{\ddot{y}\right\}}_s=a_1{\left\{\ddot{y}\right\}}_{s+1}$$ E72

$$a_1=\alpha \mathrm{\Delta }t$$ E73

$$a_2=\left(1-\alpha \right)\mathrm{\Delta }t$$ E74

8. FEM MATLAB model implementation strategy

In order to test the validity of the models discussed earlier using MATLAB, an aluminum power conductor with a steel core having a total diameter of 35.56 mm and having an ultimate tensile strength of 51.51kN was used in setting up the MATLAB simulation. Further physical properties of the power cable are shown in Table 1 . The power conductor had a minimum bending stiffness ${\mathit{EI}}_{\min }$ of 8.66 Nm² and maximum bending stiffness ${\mathit{EI}}_{\max }$ of 433 Nm². The wholistic finite element models where implemented in MATLAB using strategy expressed in Figure 2 .

Strand layer	Strand material	Diameter (mm)	No. of strands	Pitch per length (cm)	Lay direction
Layer 0	Steel	2.25	1
Layer 1	Aluminum	3.38	6	16.1	Left hand lay
Layer 2	Aluminum	3.38	12	22.2	Right hand lay

Table 1.

Power transmission conductor physical properties.

Figure 2.

FEM MATLAB implementation strategy.

The inputs in the MATLAB algorithm were bending and axial loads, the cross-sectional area of the power conductor, strand radius and strand material type. The type of analysis which can be either static or dynamic was also specified as part initial and boundary conditions. Also included in the algorithm was to specify if the computation focuses on local vibration of the power conductor or the global vibration model.

9. Experimental investigation of conductor vibration

MATLAB code was written for the FEM and this was used to model the dynamic analysis of the problem of conductor vibration. To validate the FEM model an experimental study was conducted at the Vibration and Research Testing Centre (VRTC) situated at the University of KwaZulu-Natal which comprises of apparatus similar to that shown in Figure 3 . The sweep tests (resonance search) were carried out and the test results were used to obtain natural frequencies and the modes of vibration for a Pelican conductor. The frequency range for the Pelican conductor was between 5 and 50 Hz and testing was done for three axial tensions of 20, 25, 30 and 35% of its ultimate tensile strength (UTS). The experimental results obtained were used to validate the developed FEM model. The comparison between results from the experimental data, FEM and the theoretical model for the three different axial tensions for high voltage conductors are reported in the next section.

Figure 3.

Experimental test set-up [7] .

10. Simulation and experimental results

The results from the MATLAB simulations were compared with results from the finite element models (FEM) and experimental recordings. These are shown in Figures 4 – 7 . The results were compared in terms of the natural frequency of vibration or vertical displacement of the power conductor.

Figure 4.

Frequency of vibration at 20% UTS.

Figure 5.

Frequency of vibration at 25% UTS.

Figure 6.

Frequency of vibration at 30% of UTS.

Figure 7.

Frequency of vibration at 35% of UTS.

11. Conclusion

The results showed that the implementation of the derived models in MATLAB provided a reliable strategy in the determination of the wind-induced dynamic properties of high voltage transmission lines. The results from MATLAB simulation, finite element method and experimental recordings were similar in values and showed similar trend. MATLAB as an environment can be used as a reliable simulation tool to implement and analyze high voltage conductor dynamics. The parameters obtained from the results, to some degree of accuracy can be used to predict the response of conductors due to aeolian vibration caused by wind loading.