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(\mathbf{x}, y) = u\_0(\mathbf{x}) + \frac{v}{R} - y\left(\theta\_\mathbf{x} + \frac{u(\mathbf{x})}{R}\right) \tag{49}$$

Where θ<sup>x</sup> 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 xð Þ 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(s) = b\_0 + b\_1 s \tag{50}$$

$$
\sigma(\mathbf{s}) = \mathbf{c}\_0 + \mathbf{c}\_1 \mathbf{S} + \mathbf{c}\_2 \mathbf{S}^2 + \mathbf{c}\_3 \mathbf{S}^3 \tag{51}
$$

$$\theta = \frac{dv(S)}{ds} = c\_1 + 2c\_2S + 3c\_3S^2 \tag{52}$$

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

$$
\begin{bmatrix} u \\ v \\ \theta \end{bmatrix} = \begin{bmatrix} N\_1 & 0 & N\_2 & 0 & 0 & 0 \\ 0 & N\_3 & 0 & N\_4 & N\_5 & N\_6 \\ 0 & 0 & N\_4 & 0 & N\_5 & N\_6 \end{bmatrix} \begin{bmatrix} u\_1 \\ v\_1 \\ \theta\_1 \\ u\_2 \\ u\_3 \\ v\_2 \\ \theta\_2 \end{bmatrix} \tag{53}
$$

Where

$$N\_1 = \frac{1}{2}(1 - \zeta) \tag{54}$$

$$N\_2 = \frac{1}{2}(1+\zeta) \tag{55}$$

$$N\_3 = \frac{1}{4} \left( 2 - 3\zeta + \zeta^3 \right) \tag{56}$$

$$N\_4 = \frac{1}{4} \left( 1 - \zeta - \zeta^2 - \zeta^3 \right) \tag{57}$$

$$N\_5 = \frac{1}{4} \left( 2 + 3\zeta - \zeta^3 \right) \tag{58}$$

$$N\_{\theta} = \frac{1}{4} \left( 1 - \zeta + \zeta^2 + \zeta^3 \right) \tag{59}$$

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:

High Voltage Transmission Line Vibration: Using MATLAB to Implement the Finite Element Model of a Wind… 51 http://dx.doi.org/10.5772/intechopen.75186

$$\mathbf{N}[K] = \frac{1}{2} \left[ \mathbf{N}\_a{}^T (EA) \mathbf{N}\_a \ \delta \zeta + \frac{1}{2} \right] \mathbf{N}\_B{}^T (EI) \mathbf{N}\_B \ \delta \zeta + \frac{1}{2} \left[ \mathbf{N}\_B{}^T (T) \mathbf{N}\_{B,T} \ \delta \zeta \tag{60}$$

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

$$[M] = \frac{1}{2} \left[ \dot{\boldsymbol{\mu}}^T \rho A \dot{\boldsymbol{\mu}} + \frac{1}{2} \right] \dot{\boldsymbol{g}}^T \rho A \dot{\boldsymbol{\psi}} \tag{61}$$

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

$$
\delta W = \frac{1}{2} \int F \cdot \delta u \tag{62}
$$

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

$$
\begin{bmatrix} M\_{11} & M\_{12} \\ M\_{21} & M\_{22} \end{bmatrix} \begin{pmatrix} \ddot{u} \\ \ddot{v} \end{pmatrix} + \begin{bmatrix} K\_{11} & K\_{12} \\ K\_{21} & K\_{22} \end{bmatrix} \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} F\_1 \\ F\_2 \end{pmatrix} \tag{63}
$$

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]:

$$[M]\{\ddot{y}\} + [\mathbb{C}]\{\dot{y}\} + [\mathbb{K}]\{y\} = [F] \tag{64}$$

$$\mathbb{E}\left[\overset{\wedge}{K}\right]\_{s+1} = [K]\_{s+1} + a\_3[m]\_{s+1} \tag{65}$$

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

$$a\_3 = \frac{2}{\chi \left(\Delta t\right)^2} \tag{67}$$

$$a\_4 = \frac{2}{\gamma \Delta t} \tag{68}$$

$$a\_5 = \frac{1}{\mathcal{V}} - 1\tag{69}$$

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

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

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

$$\left[\ddot{y}\right]\_{s+1} = a\_3 \left( \{y\}\_{s+1} - \{y\}\_s \right) - a\_4 \{\dot{y}\} - a\_5 \{\ddot{y}\}\_s \tag{71}$$

$$\left[\ddot{y}\right]\_{s+1} = \left\{\ddot{y}\right\}\_{s} + a\_{2}\{\ddot{y}\}\_{s} = a\_{1}\{\ddot{y}\}\_{s+1} \tag{72}$$

$$a\_1 = a \Delta t \tag{73}$$

$$a\_2 = (1 - \alpha)\Delta t\tag{74}$$
