Calculation of Lightning Transient Responses on Wind Turbine Towers

An efficient method is proposed in this paper for calculating lightning transient responses on wind turbine towers. In the proposed method, the actual tower body is simplified as a multiconductor grid in the shape of cylinder. A set of formulas are given for evaluating the circuit parameters of the branches in the multiconductor grid. On the basis of the circuit parameters, the multiconductor grid is further converted into an equivalent circuit. The circuit equation is built in frequency-domain to take into account the effect of the frequency-dependent characteristic of the resistances and inductances on lightning transients. The lightning transient responses can be obtained by using the discrete Fourier transformwith exponential sampling to take the inverse transform of the frequency-domain solution of the circuit equation. A numerical example has been given for examining the applicability of the proposed method.


Introduction
Lightning is a complex atmospheric discharge phenomenon.The recorded lightning current has a highest peak value of around 300 kA and a duration of a few hundred microseconds [1].In fact, this peak value rarely occurs in China.The median peak value is 26 kA occurring with 50% probability according to Chinese measured data [2].With respect to wind turbines (WTs), the lightning strokes can be classified into two main types, downward flash and upward flash.The WTs with the height exceeding 100 m are mostly struck by upward flash.The impacts of lightning on WTs range from disturbances on control electronics, damages to single components, such as generator and sensor, to fires resulting in a complete loss of the installation.The control system frequently suffering from lightning transient overvoltages is typically the electrical control cabinet in the bottom of the WT tower.Its damage accounts for about 40% of all lightning damage events of WTs in the west of China.Furthermore, offshore WTs are often equipped with relatively advanced equipment such as communication antenna, transponder, GPS receiver, sea marking light, air traffic warning light, fog-horns, and meteorological instrument, all of which are especially susceptible to lightning damage.Since accessing WTs offshore is not as easy as on land, lightning damage to WTs offshore can result in significantly higher costs of repair and maintenance.In view of the seriousness of lightning stroke to WTs, the lightning protection design has been regarded with more and more attention.In the lightning protection design of WTs, the need exists for determining the lightning transient responses on the WT tower since the tower body is the longest conducting path of lightning current in the WT structure.The previous methods usually modeled the WT tower as a transmission line [3,4] or a capacitance chain [5].Although the previous methods are simple, they cannot give the lightning transient responses in different parts on the tower body owing to their neglecting the structural feature of the tower body.For an improvement in the lightning transient calculation of WTs, a novel method is proposed in this paper.The proposed method represents the actual WT tower in the shape of the circular truncated cone as an equivalent cylindrical shell.The cylindrical shell is further subdivided into a discrete multiconductor grid.A set of analytic formulas are deduced for evaluating the circuit parameters of the branches in the multiconductor grid.Based on the circuit parameters, the multiconductor grid is converted into an equivalent circuit.The hybrid equation is built in the frequency-domain for the equivalent circuit, and then the discrete Fourier transform with exponential sampling is used to obtain the lightning transient responses in different parts on the tower body.As compared to the traditional time-domain method [6], the proposed method can give due consideration to the frequencydependent characteristic of the branch impedance directly in the frequency-domain and present better applicability to the lightning transient calculation.A numerical example has also been given for checking the validity of the proposed method.

Discretization Treatment of WT Tower
An actual WT tower is usually a hollow circular truncated cone, as shown in Figure 1(a).Since ℎ ≫  1 and  2 , and  2 is not much larger than  1 , the actual tower body may be simplified as a cylindrical shell, as shown in Figure 1(b).For the purpose of transient calculation, the continuous cylindrical shell can be subdivided into a multiconductor grid constituted by longitudinal and transverse branches [6,7], as shown in Figure 2(a).As the lightning current flowing through the WT tower presents significant harmonic components not exceeding a few MHz, the lightning transient responses should be calculated by the distribution parameter circuit model.However, if the branch length is taken as significantly small, an approximate calculation can be conducted referring to lumped parameter circuit model.For this reason, the branch length  has to respect the following condition: where  min is the minimum wavelength associated to the maximum frequency of the spectrum [8].The circuit parameters of the branches in the multiconductor grid are represented by the impedances and capacitances.For convenience of the parameter calculation, each transverse arc is approximately replaced by its corresponding chord, as shown in Figure 2(b), and all the branches are taken as cylindrical conductors whose radii are estimated from their respective average cross-sections.

Derivation of Impedance and Capacitance Formulas
In view of the electromagnetic couplings between the branches, the impedances and capacitances are expressed as the matrices Z and C, respectively.Z and C are symmetric matrices according to the reciprocity principle [9].The formulas for evaluating Z and C are derived later.
3.1.Impedance.Consider two transverse branches in the multiconductor grid shown in Figure 2(b).The earth is not perfectly conducting in reality, and its resistivity  has to be taken into account in the impedance calculation.For this purpose, a complex depth  is introduced into the image . . .installation [10], as shown in Figure 3.The complex depth is defined by where j = √ −1,  is the angular frequency,  the earth resistivity, and  0 the space permeability (4 × 10 −7 H/m).
The real branch  and its image branch   are considered to be symmetric about a plane which is located at a distance  below the earth surface.In terms of the Neumann's integral formula [10,11], the mutual impedance between the branches  and  is expressed by where Evaluating the double integrals in (3) leads to The self-impedance   () of the transverse branch  is obtained by substituting  =  0 into (5).In a similar way, the mutual impedance between two longitudinal branches, as shown in Figure 4, is given by In (6), letting  =  0 gives the self-impedance   () of the longitudinal branch .
After the self and mutual impedances are calculated for a coupled branch unit with  transverse or longitudinal branches by using ( 5) or ( 6), the impedance matrix Z can be formed as (7) or It is clear from (7) and ( 8) that the resistance and inductance in the impedance exhibit a pronounced frequency-dependent characteristic owing to introducing the complex depth .If the earth is assumed to be perfectly conducting ( = 0), the complex depth  is zero according to (2).Substitution of  = 0 into (5) or (6) shows that the real part of the impedance is also zero.As a result, the impedance matrix Z is further simplified as where L 0 is the constant inductance matrix which is independent on the angular frequency .

Capacitance.
In the calculation of capacitances, the frequency effect on the charge density distribution can be significantly neglected since the charge attenuation constant  =  ( is the permittivity) is rather small in the frequency range of lightning transient for both the branch conductor and the earth [12].Therefore, the capacitances of the branches can still be calculated on the assumption that the earth is perfectly conducting.On the basis of the electromagnetic analogy [9], the product of the impedance matrix Z 0 = jL 0 and the admittance matrix Y = jC of the coupled branches becomes diagonal with negative elements determined from space permeability  0 and permittivity  0 [10] as where C is capacitance matrix of the coupled branches and E is the unit matrix.Consequently, the capacitance matrix can be obtained as

Equivalent Circuit Model
During a lightning stroke, the lightning current is usually injected to the tower body from its top and then dissipated in the earth by the earthing system.This process can be illustrated in Figure 5, where the earthing system is simply represented as the earthing resistances (  ).Once the impedance and capacitance parameters are obtained using the formulas given earlier, each coupled branch unit in the multiconductor grid can equivalently be represented by a -circuit consisting of the series impedance and shunt capacitance [6,13].Figures 6(a

Circuit Equation and Its Numerical Solution
In order to take into account the frequency-dependent characteristic of the resistances and inductances in series impedances on lightning transients, the circuit equation is built in frequency-domain.The circuit components of the equivalent circuit are numbered in the sequence of shunt capacitance, series impedance, and lightning current source.
In accordance with the component number, the incidence matrix of the equivalent circuit is written as and the relevant current vector is where the subscripts , im, and  denote shunt capacitance, series impedance, and lightning current source, respectively.From Kirchhoff 's current law, the node current equation is derived from ( 12) and (13) as Substituting   Y     U  () for   I  () in (14) gives where Y  is the shunt capacitive admittance matrix.Let U  () and Z im be the node voltage vector and series impedance matrix, respectively.The branch equation of the series impedance is given by Hence, the hybrid equation of the equivalent circuit can be built by merging ( 15) and ( 16) as On the right hand of ( 17), lightning current source   () is the Fourier transform of the injected lightning current  (see Figure 5). is usually expressed by the double exponential function [14]  =  ( − −  − ) .
In accordance with the Chinese national standard [2], the waveform parameter of  is 2.6/50 s.This gives  = 1.5 × 10 −2 s −1 and  = 1.86 s −1 by using the least square fitting method [15].The peak current  takes the value of 100 kA in light of the severe lightning stroke condition.Thus, finding the Fourier transform of (18) gives By taking the   () as the excitation, the frequency responses U  () and I im () can be obtained by solving (17).For each given angular frequency , (17) is a system of complex linear algebraic equations and can be generally written as where where  is the order of the coefficient matrix in (17).
The Gauss elimination with column pivoting is employed    to solve (20).It consists of two main phases: the forward elimination and the back substitution.After performing −1 step elimination operations, the original matrix Λ in ( 20) is converted into an upper triangular matrix (1)  11 . . .
The detailed elimination procedure has been stated in [16,17].Subsequently, the solution of vector X() is obtained by the back substitution operations  Fourier transform with exponential sampling is chosen for this purpose [18,19].The inverse transform method takes the sampling frequency in an exponential manner, that is, Δ,  2 Δ,  3 Δ, . . .,   Δ, . . .,    Δ, where Δ (= 2Δ) is the angular frequency step,   the total number of frequency samples, and  the frequency ratio between two neighbor sampling points.The value of  is experientially chosen in range of 1.03 to 1.05 for solution of the fast electrical transients [20].In (24), the maximum sampling frequency   (=   Δ/2) is determined from the frequency spectrum of lightning current and may be of the order (0.5∼5) MHz [19], depending on the waveform of lightning current.The total number of frequency samples   is evaluated by Therefore, for a frequency response   (), its time response   () can be obtained by where  is the standard smoothing factor Owing to a high accuracy and long observation time with a small number of sampling steps, (26) is suitable for the lightning transient calculation.

A Numerical Example
A typical Chinese-built WT with 1.5 MW is considered here for practical application purpose.The dimensions of the WT tower are and ℎ = 82m,  1 = 2.7m,  2 = 4.3 m, and  = 0.025 m (see Figure 1(a)).The relevant discrete multiconductor grid is shown in Figure 7, where  0 = 3.5 m.The earthing resistance   is 3.5Ω.The transient responses on the tower body are calculated by using the proposed method.The waveforms of current on branch  2,1 −  3,1 and potential at node  1,1 are shown in Figure 8, where the corresponding waveforms excluding the frequency-dependent characteristic of the resistances and inductances (see (9)) are also given for comparison.As seen from Figure 8 that the waveforms excluding the frequency-dependency contain the obvious numerical oscillations, while the waveforms including frequency-dependency are relatively smooth.This indicates that inclusion of the frequency-dependency can give a higher accuracy for the lightning transient calculation.In order to check the validity of the proposed method, the peak potentials at different nodes on the tower body are calculated by the PSCAD/EMTDC software, the time-domain method [6], and the proposed method, respectively.The calculated results are shown in Figure 9.The PSCAD/EMTDC software used here has been specifically developed by Manitoboa HVDC.
Research Centre for electromagnetic transient simulation and widely applied in electric power systems.It is clear from Mathematical Problems in Engineering  Figure 9 that a close agreement appears between the potential values obtained from the proposed method and those from the PSCAD/EMTDC software and the time-domain method.

Conclusions
The lightning transient calculation has been performed in this paper for obtaining the transient responses on WT towers.Discretization representation of the actual tower body in the shape of circular truncated cone as a cylindrical multiconductor grid makes a significant simplification to the transient calculation.A set of analytical formulas have been given for evaluation of the circuit parameters of the coupled branches in the multiconductor grid, and the discrete Fourier transform with exponential sampling has been employed to calculate the transient responses on the WT tower.This calculation procedure can take into account the effect of frequency-dependent characteristic of resistances and inductances on lightning transients.The practical applicability of the proposed method has been examined by a numerical example of 1.5 MW WT tower, which shows that proposed method is useful in lightning transient calculation of WT towers.
) and 6(b) show the -circuit representing three coupled branches.If the circuit parameters are expressed in matrix form, Figure6(b) can be further depicted as a compact -circuit, as shown in Figure6(c).With all the coupled branch units in the multiconductor grid replaced by the -circuits, Figure5is converted into an equivalent circuit consisting of resistances, inductances, and capacitances.The lightning current is modeled as a current source and applied to the top node of the equivalent

Figure 5 :
Figure 5: Injection of lightning current to tower body.

Figure 9 :
Figure 9: Peak potentials at different nodes on the tower body.