Numerical simulation of corona-induced vibration of high voltage conductor

Abstract

When it rains, electric power transmission lines start vibrating due to corona effect. This type of vibration is known as “corona-induced vibration”. The aim of this paper is to elaborate a mathematical model for numerical simulation of the corona-induced vibration, with consideration of the influence of the magnitude and the polarity of the electric field on the conductor surface. Finite element method was employed to develop the numerical model, and the finite difference method was used for the time discretisation. The moment of application of the corona-induced force is evaluated using the resultant vertical force applied to a water drop, suspended under a high voltage conductor. Some experimental results of other authors are exploited to evaluate the precision of the simulation and the validation of numerical results.

Keywords

corona-induced vibration / corona wind / finite element method

Introduction

One of the consequences of high voltage electric power systems is the corona effect. This phenomenon is the source of electromagnetic interference, audible noises, important energy losses and mechanical vibrations. This latter consequence, called “corona-induced vibration”, can lead to the fatigue of overhead conductors and supporting elements [1]. It has been established that the intermittent presence of corona space charge and the ionic wind are the main causes of this phenomenon. Research in this field began in 1970 by an analytical study focusing on the determination of vibration amplitudes. Following studies, this time with experimental detail, realized in Canada, led to a significant result and especially interest in the mechanism of vibrations [2,3]. Then, in 1986 a precise mechanism was proposed and accepted by the scientific community [4]. During these years many researchers have studied the different aspects of this subject. Diverse experimental models and laboratory mechanisms were used to simulate this phenomenon. However, most of the researchers accomplished the successive results but few of these results were based on a numerical model and numerical simulation. Therefore, the present work is based on the numerical simulations of the corona-induced vibration [5]. Two simulation techniques are used: the modal superposition for the discretisation of the movement, and the central difference method for the discretisation of the time. Some experimental results of other authors are exploited to evaluate the accuracy of the numerical simulation.

Description of vibration mechanism

The vibration mechanism can be described by the following steps (Fig. 1) [6-11]:

1) The conductor is attracted to the ground surface, due to the electric image force.

2) Under wet conditions and in the presence of electric field, suspended drops are formed at the lower surface of the conductor.

3) The suspended drops at the bottom of the conductor surface take on a conical shape. The formation of cones results from the interaction between the forces due to the electrostatic field on the surface of the conductor, surface tension and gravity.

4) Due to the field intensification at the tip of the cones, corona discharge increases the space charge around the suspended water drops.

5) The increase of the space charge around water drops produces a partial shielding effect between the conductor and the ground.

6) The electric image force is eliminated and thus the conductor moves upward.

7) The suspended drops reach a critical size, instability occurs and a certain quantity of the water drop is ejected from the conductor.

8) Corona discharge is attenuated and therefore the current decreases to a small value. The water drops remaining on the conductor surface do not have the conical shape, consequently the corona discharge is very weak and the space charge becomes small.

9) Therefore, there is no shielding produced by the space charge and thus the attractive force between the conductor and the ground becomes important again; the conductor moves downward.

10) As the rain continues, more water flows, and the conical shape of the drops produces again the space charge, thus the process is repeated (continuous vibration).

Development of basic differential equation

We first determine the basic differential equation that describes the tense-conductor vertical displacements submitted to a distributed external force.

First, let us consider the governing differential equations of the suspended cable:

(1)\rho (x){\displaystyle \frac{{\partial }^{\mathrm{2}}U(x,t)}{\partial t^{\mathrm{2}}}}+\mu (x){\displaystyle \frac{\partial U(x,t)}{\partial t}}-{\displaystyle \frac{\partial }{\partial x}}\left[\alpha (x){\displaystyle \frac{\partial U(x,t)}{\partial x}}\right]=f(x,t),

where \rho (x){\displaystyle \frac{{\partial }^{\mathrm{2}}U(x,t)}{\partial t^{\mathrm{2}}}} is inertia force term, \mu (x){\displaystyle \frac{\partial U(x,t)}{\partial t}} is damping force term, {\displaystyle \frac{\partial }{\partial x}}\left[\alpha (x){\displaystyle \frac{\partial U(x,t)}{\partial x}}\right] is conductor tension force term, f(x,t) is external forces, U(x,t) is vertical conductor displacement.

Equation (1) is solved using the finite element method; this method is a numerical analysis technique for obtaining approximate solutions to a wide variety of engineering problems. The GALERKIN technique leads to the decrease in the integration order [12,13].

The typical pondered residual equation can be written as follows:

(2){\displaystyle \underset{}{\overset{e}{\int }}R(x,t;a){\mathrm{\Phi }}_i(x)\mathrm{d}x=\mathrm{0}},

where

R(x,t;a)=\rho (x){\displaystyle \frac{{\partial }^{\mathrm{2}}U(x,t)}{\partial t^{\mathrm{2}}}}+\mu (x){\displaystyle \frac{\partial U(x,t)}{\partial t}}-{\displaystyle \frac{\partial }{\partial x}}\left[\alpha (x){\displaystyle \frac{\partial U(x,t)}{\partial x}}\right]-f(x,t),

and {\mathrm{\Phi }}_i(x) is a polynomial interpolation function.

When developing the integral part, we obtain

(3){\displaystyle \underset{}{\overset{e}{\int }}{\mathrm{\Phi }}_i^e(x)\rho (x){\displaystyle \frac{{\partial }^{\mathrm{2}}{\tilde{U}}^e(x,t)}{\partial t^{\mathrm{2}}}}\mathrm{d}x}+{\displaystyle \underset{}{\overset{e}{\int }}{\mathrm{\Phi }}_i^e(x)\mu (x){\displaystyle \frac{\partial {\tilde{U}}^e(x,t)}{\partial t}}\mathrm{d}x}+{\displaystyle \underset{}{\overset{e}{\int }}{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_i^e(x)}{\mathrm{d}x}}\alpha (x){\displaystyle \frac{\partial {\tilde{U}}^e(x,t)}{\partial x}}\mathrm{d}x}={\displaystyle \underset{}{\overset{e}{\int }}f(x,t){\mathrm{\Phi }}_i^e(x)\mathrm{d}x}-{\left[\left(-\alpha (x){\displaystyle \frac{\partial {\tilde{U}}^e(x,t)}{\partial x}}\right){\mathrm{\Phi }}_i^e(x)\right]}_{x_{\mathrm{1}}}^{x_n}.

We adopt now an approximate solution to the problem:

(4){\tilde{U}}^e(x,t;a)={\displaystyle \sum _{j=\mathrm{1}}^n{a}_j(t){\mathrm{\Phi }}_j^e(x)},

where a_j(t) represents the values of the function \tilde{U} at the different nodes, and n is the liberty degrees.

When substituting the approximate solution and its derivative:

\begin{array}{l}{\displaystyle \frac{\partial {\tilde{U}}^e(x,t)}{\partial x}}={\displaystyle \sum _{j=\mathrm{1}}^n{a}_j(t)}{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_j(x)}{\mathrm{d}x}},\hfill \\ {\displaystyle \frac{\partial {\tilde{U}}^e(x,t)}{\partial t}}={\displaystyle \sum _{j=\mathrm{1}}^n{\displaystyle \frac{\mathrm{d}{a}_j(t)}{\mathrm{d}t}}}{\mathrm{\Phi }}_j(x),\hfill \\ {\displaystyle \frac{{\partial }^{\mathrm{2}}{\tilde{U}}^e(x,t)}{\partial t^{\mathrm{2}}}}={\displaystyle \sum _{j=\mathrm{1}}^n{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{a}_j(t)}{\mathrm{d}{t}^{\mathrm{2}}}}}{\mathrm{\Phi }}_j(x)\hfill \end{array}

into Eq. (3), we obtain:

(5)M^e{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}a(t)}{\mathrm{d}{t}^{\mathrm{2}}}}+C^e{\displaystyle \frac{\mathrm{d}a(t)}{\mathrm{d}t}}+K^{e}a(t)=F^e(t),

where

\begin{array}{l}{M}_{ij}^e={\displaystyle \underset{}{\overset{e}{\int }}{\mathrm{\Phi }}_i^e(x)\rho (x){\mathrm{\Phi }}_j^e(x)\mathrm{d}x},\hfill \\ C_{ij}^e={\displaystyle \underset{}{\overset{e}{\int }}{\mathrm{\Phi }}_i^e(x)\mu (x){\mathrm{\Phi }}_j^e(x)\mathrm{d}x},\hfill \\ K_{ij}^e={\displaystyle \underset{}{\overset{e}{\int }}{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_i^e(x)}{\mathrm{d}x}}\alpha (x){\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_j^e(x)}{\mathrm{d}x}}\mathrm{d}x},\hfill \\ F_i^e(t)={\displaystyle \underset{}{\overset{e}{\int }}f(x,t){\mathrm{\Phi }}_i^e(x)\mathrm{d}x}-{[{\tau }^e(x,t){\mathrm{\Phi }}_i^e(x)]}_{x_{\mathrm{1}}}^{x_n},\hfill \\ {\tau }^e(x,t)=-\alpha (x){\displaystyle \frac{\partial {\tilde{U}}^e(x,t)}{\partial x}}.\hfill \end{array}

We develop a specific expression of {\mathrm{\Phi }}_i(x) corresponding to a linear element (Fig. 2) [7]:

(6)\begin{array}{l}{\mathrm{\Phi }}_{\mathrm{1}}^e(x)={\displaystyle \frac{x_{\mathrm{2}}-x}{x_{\mathrm{2}}-x_{\mathrm{1}}}},{\mathrm{\Phi }}_{\mathrm{2}}^e(x)={\displaystyle \frac{x-x_{\mathrm{1}}}{x_{\mathrm{2}}-x_{\mathrm{1}}}},\hfill \\ {\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{1}}^e(x)}{\mathrm{d}x}}=-{\displaystyle \frac{\mathrm{1}}{x_{\mathrm{2}}-x_{\mathrm{1}}}},{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_2^e(x)}{\mathrm{d}x}}={\displaystyle \frac{\mathrm{1}}{x_{\mathrm{2}}-x_{\mathrm{1}}}}.\hfill \end{array}

Substituting Eq. (6) into Eq. (5), we obtain

\left[\begin{array}{ll}{\displaystyle \frac{{\rho }^{e}L}{\mathrm{3}}}\hfill & {\displaystyle \frac{{\rho }^{e}L}{\mathrm{6}}}\hfill \\ {\displaystyle \frac{{\rho }^{e}L}{\mathrm{6}}}\hfill & {\displaystyle \frac{{\rho }^{e}L}{\mathrm{3}}}\hfill \end{array}\right]\left\{\begin{array}{l}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{a}_{\mathrm{1}}(t)}{\mathrm{d}{t}^{\mathrm{2}}}}\hfill \\ {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{a}_{\mathrm{2}}(t)}{\mathrm{d}{t}^{\mathrm{2}}}}\hfill \end{array}\right\}+\left[\begin{array}{ll}{\displaystyle \frac{{\rho }^{e}L}{\mathrm{3}}}\hfill & {\displaystyle \frac{{\rho }^{e}L}{\mathrm{6}}}\hfill \\ {\displaystyle \frac{{\rho }^{e}L}{\mathrm{6}}}\hfill & {\displaystyle \frac{{\rho }^{e}L}{\mathrm{3}}}\hfill \end{array}\right]\left\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{a}_{\mathrm{1}}(t)}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}{a}_{\mathrm{2}}(t)}{\mathrm{d}t}}\end{array}\right\}+\left[\begin{array}{ll}{\displaystyle \frac{{\alpha }^e}{L}}\hfill & {\displaystyle \frac{-{\alpha }^e}{L}}\hfill \\ {\displaystyle \frac{-{\alpha }^e}{L}}\hfill & {\displaystyle \frac{{\alpha }^e}{L}}\hfill \end{array}\right]\left\{\begin{array}{l}{a}_{\mathrm{1}}(t)\\ a_{\mathrm{2}}(t)\end{array}\right\}=\left\{\begin{array}{l}{\displaystyle \frac{f^{e}L}{\mathrm{2}}}+{\tau }^e(x,t)\\ {\displaystyle \frac{f^{e}L}{\mathrm{2}}}-{\tau }^e(x,t)\end{array}\right\},

and assembling of all elements

(7)\mathit{M}\left\{\stackrel{••}{a}\right\}+\mathit{C}\left\{\stackrel{•}{a}\right\}+\mathit{K}\{a\}=\{F\},

where a_i(t) represents the displacements of each node.

At this step, we integrate Eq. (7) in order to determine the amplitude of the vibrations. To solve this equation, modal superposition should be a good choice and there will be no need to use a direct integration.

Modal superposition method

This method allows transposing basic differential equations to a system of independent equations, while each equation represents one mode of vibrations [7,13].

First, let us calculate the eigenvalues (frequencies) and eigenvectors (modes) of the system.

The undamped modes in Eq. (7) are time-harmonic (sinusoidal) solutions when there are no loads ({F} = {0}) and damping is neglected (BoldItalic= 0); that is

(8)\mathit{M}\left\{\stackrel{••}{a}\right\}+\mathit{K}\{a\}=\{\mathrm{0}\}.

Let

(9)\{a\}=\{v\}{\mathrm{e}}^{\mathrm{i}\omega t},

substituting Eq. (9) into Eq. (8) yields

(10)\mathit{K}\{a\}-\lambda \mathit{M}\{a\}=\{\mathrm{0}\}.

This latter equation represents the generalized algebraic eigen problem where the eigenvalue \lambda is the square of the circular frequency, that is,

{\omega }_i^{\mathrm{2}}={\lambda }_i.

The eigenvectors are orthogonal with respect to BoldItalicand BoldItalic, and they are orthonormal with respect to BoldItalic, that is,

{\{v\}}_i^{\mathrm{T}}\mathit{K}{\{v\}}_j=\{\begin{array}{ll}{\omega }_i^{\mathrm{2}},\hfill & i=j,\hfill \\ \mathrm{0},\hfill & i\ne j,\hfill \end{array}

{\{v\}}_i^{\mathrm{T}}\mathit{M}{\{v\}}_j=\{\begin{array}{ll}\mathrm{1},\hfill & i=j,\hfill \\ \mathrm{0},\hfill & i\ne j.\hfill \end{array}

It is convenient to define an N-by-N square matrix BoldItalic that contains the N eigenvectors as columns:

\mathit{V}=\{v_{\mathrm{1}},v_{\mathrm{2}},\mathrm{...},v_n\},

and an N-by-N diagonal matrix {\mathbf{\Omega }}^{\mathrm{2}} that contains the N eigenvalues on the diagonal.

The global solution to Eq. (7) can be written as a linear superposition of the N modes:

(11)\{a(t)\}={\displaystyle \sum _{j=\mathrm{1}}^n{A}_j(t)\{v\}}=\mathit{V}\{A(t)\}.

Substituting Eq. (11) into Eq. (7), pre-multiplying the latter by {\mathit{V}}^{\mathrm{T}}, and using the orthogonal relationships developed earlier, we can write

(12)\left\{\stackrel{••}{A}(t)\right\}+{\mathit{V}}^{\mathrm{T}}\mathit{C}\mathit{V}\left\{\stackrel{•}{A}(t)\right\}+{\mathbf{\Omega }}^{\mathrm{2}}\left\{\stackrel{}{A(t)}\right\}={\mathit{V}}^T\{F\}.

Equation (12) can be decoupled, combining the damping factors on the diagonal:

(13){\mathit{V}}^{\mathrm{T}}\mathit{C}\mathit{V}=\left[\begin{array}{lllll}\mathrm{2}{\omega }_{\mathrm{1}}^{\mathrm{2}}\xi & \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{0}& \mathrm{2}{\omega }_{\mathrm{1}}^{\mathrm{2}}\xi & \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{2}{\omega }_{\mathrm{1}}^{\mathrm{2}}\xi & \cdots & \mathrm{0}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \cdots & \mathrm{2}{\omega }_{\mathrm{1}}^{\mathrm{2}}\xi \end{array}\right].

Substituting Eq. (13) into Eq. (12), enabling them to be written as N separate equations:

(14)\begin{array}{ll}\stackrel{••}{A}(t)+\mathrm{2}{\omega }_i^{\mathrm{2}}{\xi }_i\stackrel{•}{A}(t)+{\omega }_i^{\mathrm{2}}A(t)=f_i(t),& i=\mathrm{1},\mathrm{2},\cdots ,m,\end{array}

where {\xi }_i are the modal damping ratios, and m is the necessary mode number to really represent the system, and

f_i(t)={\{v\}}_i^{\mathrm{T}}\{F(t)\}.

The system of Eq. (14) can be numerically solved by using one of the time-stepping methods.

Central difference method

One can use the central difference method to solve the system of Eq. (14). This method requires three times either: t_{n-\mathrm{1}}, t_n or t_{n-\mathrm{2}}, the system is evaluated at the central time [1,6,7]:

(15)\stackrel{••}{A}{(t)}_{n-\mathrm{1}}+{\mathit{C}}_D\stackrel{•}{A}{(t)}_{n-\mathrm{1}}+{\mathbf{\Omega }}^{\mathrm{2}}A{(t)}_{n-\mathrm{1}}={\{f\}}_{n-\mathrm{1}}.

The two derivatives are approximated by central differences:

(16){\left\{\stackrel{•}{A}\right\}}_{n-\mathrm{1}}={\displaystyle \frac{{\{A\}}_n-{\{A\}}_{n-\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }t}},

(17){\left\{\stackrel{••}{A}\right\}}_{n-\mathrm{1}}={\displaystyle \frac{{\{A\}}_n-\mathrm{2}{\{A\}}_{n-\mathrm{1}}+{\{A\}}_{n-\mathrm{2}}}{\mathrm{\Delta }{t}^{\mathrm{2}}}}.

Substituting Eqs. (16) and (17) into Eq. (15) yields:

(18){\{A\}}_n={\displaystyle \frac{\mathrm{2}-{\omega }_i^{\mathrm{2}}\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{1}+{\omega }_i{\xi }_i\mathrm{\Delta }t}}{\{A\}}_{n-\mathrm{1}}+{\displaystyle \frac{{\omega }_i\xi \mathrm{\Delta }t-\mathrm{1}}{\mathrm{1}+{\omega }_i{\xi }_i\mathrm{\Delta }t}}{\{A\}}_{n-\mathrm{2}}+{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{1}+{\omega }_i{\xi }_i\mathrm{\Delta }t}}{\{f\}}_{n-\mathrm{1}}.

We can now solve the system of Eq. (18) to find the modal amplitudes, and then each A_i(t) summed according to Eq. (11) to produce \{a(t)\}.

Corona-induced force calculation

To solve the system of Eq. (18), it is necessary to evaluate the variation of corona-induced force as a function of time. According to Farzaneh [8], the value of this force is greater just before the ejection, when the drop reaches its maximal length, and smaller for all other vibration cycles. Furthermore, it has been observed that there is synchronization between the drops ejection and the conductor movement [7,8,11,12]. The ejection of drops always occurs at the lowest position of the conductor. For this reason, we can represent the corona-induced force by an impulse force with conserving the same quantity of energy transmitted to the conductor while using a sinusoidal force that has been evaluated by Farzaneh in the laboratory (as a function of the electric field intensity and polarity), as

(19)F_{\mathrm{i}\mathrm{m}\mathrm{p}}=\mathrm{1.84}{F}_{\sin }.

The period of application of the impulse force has been estimated to 20 ms, and the phase shift between the moment when the conductor arrives at its lower position and the moment of application of the force has been estimated to be 10 ms [7].

Moment of application of corona force

The moment of application of the corona force is evaluated, while comparing the vertical force balance applied to a water drop suspended under a high voltage (HV) conductor in movement [6,7,11]. The forces applied to the water drop are:

1) Gravity force

(20)F_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}=DV{a}_{\mathrm{G}},

where D is the density of water, V is the volume of the drop, a_G is the acceleration of gravity.

The density of water and the acceleration of gravity are constant, but the volume of the drop is variable.

2) Inertia force

(21)F_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{a}}=DV{a}_{\mathrm{M}},

where a_M is the acceleration of movement.

Acceleration is calculated by using the second derivative of displacement; this is given in Eq. (5).

3) Electrostatic force

(22)F_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{o}}=\mathrm{58.06}\times {\mathrm{10}}^{-\mathrm{18}}{E}^{\mathrm{2}}{r}^{\mathrm{2}},

where r is the average radius of the suspended drop, E is the value of the electrical field on the conductor surface.

4) Surface tension

(23)F_{\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}}=\mathrm{2}\mathrm{\pi }r\gamma ,

where \gamma =\mathrm{1.28}{\mathrm{e}}^{-\mathrm{2}}\mathrm{N}/\mathrm{m} at 20ºC.

5) Induced corona force

This force is due to the charge carriers moving from a high voltage conductor under rain to the ground. Farzaneh has measured the value of this force at the ejection time for each drop [8] and he reported this value as

(24)F_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{o}\mathrm{n}\mathrm{a}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{N}/\mathrm{d}\mathrm{r}\mathrm{o}\mathrm{p}(\mathrm{o}\mathrm{r}{\mathrm{10}}^{-\mathrm{3}}\mathrm{N}/\mathrm{t}\mathrm{w}\mathrm{o}\mathrm{d}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{s}).

At each iteration, when the sum of the three forces (inertia, gravity, and electrostatic) becomes larger than the sum of the other two forces (corona force and tension force), the ejection of the water drop occurs.

Drop’s volume

The volume of a water drop suspended under a high voltage conductor is a time variable. When a drop enlarges sufficiently, depending on the electrical field strength on the surface of the conductor and the rain precipitation, it will be ejected at a certain level of elongation. The ejection cannot detach the whole volume of the drop from the conductor. The remaining volume of the drop after ejection is called residual. Parameter K is defined as the ratio of drop volume before and after ejection. For a specific electrical field and rain precipitation, ratio K and the maximum volume of water drop just before ejection are inserted into the numerical model as input values. Based on this data, the program calculates the residual volume of the water after ejection and uses that result for the calculation of the applied force to a suspended water drop.

The calculation of the volume water supplied to the conductor at each time step is realized as follows [7]:

{\displaystyle \frac{p(\mathrm{m}\mathrm{m}/\mathrm{h})\times l(\mathrm{m}\mathrm{m})\times d(\mathrm{m}\mathrm{m})}{\mathrm{3600}(\mathrm{s}/\mathrm{h})}}=\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}(\mathrm{m}{\mathrm{m}}^{\mathrm{3}}/s),

where p is intensity of rain precipitation, l is total length of the conductor, d is diameter of the conductor.

Numerical simulation

A calculation program has been developed by using MATLAB software, in order to evaluate the amplitude of corona-induced vibration. The modal superposition method used in this study, requires first the calculation of eigenvalues of the system and the corresponding eigenvectors. The conductor’s initial position is calculated with the following equation [6]:

y={\displaystyle \frac{s}{w}}\times \left(\cosh {\displaystyle \frac{wl}{\mathrm{2.0}s}}-\mathrm{1}\right),

where s is mechanical tension applied to the extremities of the conductor, w is the weight of the conductor by unit of length.

Then, at each time step, we calculate the conductor’s displacement by using Eq. (19). The resolution of the uncoupled equation system is achieved with a time step equal to 0.001 s.

To validate the numerical simulations presented in this paper, we decided to compare the numerical results with the experiment obtained for the same conditions and situations. Since our numerical model is based on the physical model used by Farzaneh, we opted to make a comparison of the numerical results with his results.

The model is an aluminum conductor steel reinforced (ACSR) conductor, 3.58 m in length and 3.05 cm in diameter, the mass of the conductor is 5.92 kg. The conductor is subjected to an artificial rain and is placed along the axis of a cylindrical metallic cage. The cage has an interior diameter of 1 m and the cylindrical cage has a length of 2.0 m. Two insulators of 1.0 m height support the two extremities of the conductor [8,10].

Results and discussion

Figures 3-5 illustrate the variation of the central node versus electrical field for the two polarities, under a rain precipitation rate above 25 mm/h. The curves are shown for both the numerical and the experimental simulation.

One can observe that for the same value of electrical field supplied to the ACSR conductor, negative direct current (DC) has the most amplitude of vibration while the alternative has less. These results can be explained by the fact that at the same applied voltage, corona activities in the air are usually stronger under negative fields than under positive fields, because of the lower negative corona-onset voltage. The relatively low amplitude of the vibration under alternating current (AC) fields, compared to that under DC, may be attributed to the residual corona-space charge from the previous half cycle of each period of instantaneous voltage.

Note that experimental and simulation curves are concordant; there is a little divergence between the two curves but have nearly the same behavior. The differences between the numerical and experimental results could be explained as follows:

First, laboratory experimentations may include some errors, due to the imprecision of measurement devices or the wrong readings, especially for very small or very high fields.

Second, the real number of water drops suspended to the conductor can vary from a cycle of vibrations to others.

Third, application moments and the corona force can cause the divergences between the numerical and experimental results.

Figure 6 shows the variation of amplitude vibrations of the central node during 20 seconds, and Fig. 7 during the first 4 seconds, for an electric field of 13.9 kV/cm. It can be seen that there are fluctuations of the vibration amplitude at the beginning, and then the vibration becomes stable. Indeed, it has been observed experimentally that at the beginning, ejected drops and discharge current are distributed randomly; and at steady state vibration,le ejection of drops is synchronized with the movement of the conductor [12].

Furthermore, the value of vibration frequency is of the same order as the natural frequency of the conductor (5.6 Hz). This result is in agreement with the experimental results.

Conclusions

The results obtained from the research work presented in this paper, which are related to the corona-induced vibration, leads to the following conclusions:

1) The finite element method has been permitted to simulate the corona-induced vibration.

2) The modal superposition method that is employed in this paper allows making use of smaller time steps and consequently increases the simulation accuracy.

3) The corona-induced force has an impulse shape. It gives more accurate results than a sinusoidal shape as used in the earlier studies.

4) The variation of conductor’s surface electric field or its polarity causes the variation of the vibration amplitude.

5) For the same value of conductor’s surface electrical field, negative DC has the most amplitude of vibration and the alternative has less.

6) At steady state vibration, ejection drop is synchronized with the movement of the conductor and the vibration frequency is of the same order of magnitude as the natural frequency of the conductor.

7) The numerical model designed for the laboratory model, described in this paper, is so close to the reality that it can be a good base for the numerical simulation of high voltage transmission lines.

References

[1] Burnett D S. Finite Element Analysis. Addison-Wesley Publishing Company, 1988

[2] Kawamoto H, Umezu S. Electrohydrodynamic deformation of water surface in a metal pin to water plate corona discharge system. Journal of Physics D: Applied Physics, 2005, 38(6): 887-894

[3] Kawamoto H, Umezu S. Force at spark discharge in pin-to-plate system. Journal of Electrostatics, 2007, 65(2): 75-81

[4] Kollar L E, Farzaneh M, Karev A R. Modeling droplet collision and coalescence in an icing wind tunnel and the influence of these processes on droplet size distribution. International Journal of Multiphase Flow, 2005, 31(1): 69-92

[5] Kollar L E, Farzaneh M. Vibration of bundled conductors following ice shedding. IEEE Transactions on Power Delivery, 2008, 23(2): 1097-1104

[6] Derakhshanin M. Simulations numériques des vibrations induites par effet de couronne sur un court conducteur soumis à une pluie artificielle. Dissertation for the Master’s Degree. Chicoutimi: Université of Chicoutimi, 2001

[7] Demers P. Simulations numériques des vibrations induites par effet de couronne sur les conducteurs à haute tension. Dissertation for the Master’s Degree. Chicoutimi: Université of Chicoutimi, 1994

[8] Farzaneh M. Contribution à l’étude des mécanismes de vibrations induites par effet de couronne. Dissertation for the Doctoral Degree. Toulouse: Université of Paul Sabatier, 1986

[9] Farzaneh M, Teisseyre Y. Mechanical vibration of HV conductors induced by corona: roles of the space charge and ionic wind. IEEE Transactions on Power Delivery, 1988, 3(3): 1122-1130

[10] Farzaneh M. Effects of the intensity of precipitation and transverse wind on the corona-induced vibration of HV conductors. IEEE Transactions on Power Delivery, 1992, 7(2): 674-680

[11] Hamel M. Influence de la variation de la température ambiante sur les vibrations induites par effet de couronne. Dissertation for the Master’s Degree. Chicoutimi: université of Chicoutimi, 1991

[12] Phan L C, Adachi T, Allaire M C. Experimental investigations of corona-induced vibration on high voltage conductors with different types of supports. IEEE Transactions on Power Apparatus and Systems, 1981, PAS-100(4): 1975-1984

[13] Staub C. Modélisation dynamique de procédés de forgeage. Dissertation for the Doctoral Degree. Lyon: The National Applied Science Institute, 1998

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg