# Calculation of Lightning Transient Responses on Wind Turbine Towers

Authors: Xiaoqing Zhang; Yongzheng Zhang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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 transform with 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.


---

## Body

## 1. 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 frequency-dependent 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.


## 2. Discretization Treatment of WT Tower

An actual WT tower is usually a hollow circular truncated cone, as shown in Figure1(a). Since h≫r1 and r2, and r2 is not much larger than r1, 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 l has to respect the following condition:
(1)l<ξmin10,
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.Geometric simplification of WT tower: (a) hollow circular truncated cone; (b) cylindrical shell.
(a)
(b)Multiconductor grid: (a) cylinder; (b) prism.
(a)
(b)


## 3. Derivation of Impedance and Capacitance Formulas

In view of the electromagnetic couplings between the branches, the impedances and capacitances are expressed as the matricesZ 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 Figure2(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 d is introduced into the image installation [10], as shown in Figure 3. The complex depth is defined by
(2)d=ρjωμ0,
where j=-1, ω is the angular frequency, ρ the earth resistivity, and μ0 the space permeability (4π×10-7 H/m). The real branch j and its image branch j′ are considered to be symmetric about a plane which is located at a distance d below the earth surface. In terms of the Neumann’s integral formula [10, 11], the mutual impedance between the branches j and k is expressed by
(3)zjk(ω)=Rjk(ω)+jωLjk(ω)=jωμ04π[∫x1x1+lh∫x1x1+lhdxjdxksjk-∫x1x1+lh∫x1x1+lhdxj′dxksjk′],
where
(4)sjk=(xj-xk)2+a2,sjk′=(xk-xj′)2+(2y1+a+2d)2.
Evaluating the double integrals in (3) leads to
(5)zjk(ω)=jωμ02π×[lhsinh-1lha-a2+lh2z+a-lhsinh-1lh2y1+a+2d+lh2+(2y1+a+2d)2-(2y1+a+2d){lhsinh-1lha-a2+lh2z}].
The self-impedance Zjj(ω) of the transverse branch j is obtained by substituting a=r0 into (5). In a similar way, the mutual impedance between two longitudinal branches, as shown in Figure 4, is given by
(6)Zmn(ω)=jωμ02π×[(2y1+lv+2d)sinh-12y1+lv+2db+(2y1+lv+2d)2+b2-(y1+lv+d)×sinh-12(y1+lv+d)b-124(y1+lv+d)2+b2-(y1+d)sinh-12(y1+d)b+124(y1+d)2+b2+lvsinh-1lvb-b2+lv2+b{(2y1+lv+2d)sinh-12y1+lv+2db}].
In (6), letting b=r0 gives the self-impedance Zmm(ω) of the longitudinal branch m.Figure 3
Two parallel transverse branches.Figure 4
Two parallel longitudinal branches.After the self and mutual impedances are calculated for a coupled branch unit withN transverse or longitudinal branches by using (5) or (6), the impedance matrix Z can be formed as
(7)Z=[Rjk(ω)]N×N+jω[Ljk(ω)]N×N=[zjk(ω)]N×N
or
(8)Z=[Rmn(ω)]N×N+jω[Lmn(ω)]N×N=[zmn(ω)]N×N.
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 d. If the earth is assumed to be perfectly conducting (ρ=0), the complex depth d is zero according to (2). Substitution of d=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
(9)Z0=jωL0,
where L0 is the constant inductance matrix which is independent on the angular frequency ω.


### 3.2. 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 Z0=jωL0 and the admittance matrix Y=jωC of the coupled branches becomes diagonal with negative elements determined from space permeability μ0 and permittivity ε0 [10] as
(10)Z0·Y=(jωL0)·(jωC)=-(ω2ε0μ0)E,
where C is capacitance matrix of the coupled branches and E is the unit matrix. Consequently, the capacitance matrix can be obtained as
(11)C=ε0μ0Z0.


## 3.1. Impedance

Consider two transverse branches in the multiconductor grid shown in Figure2(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 d is introduced into the image installation [10], as shown in Figure 3. The complex depth is defined by
(2)d=ρjωμ0,
where j=-1, ω is the angular frequency, ρ the earth resistivity, and μ0 the space permeability (4π×10-7 H/m). The real branch j and its image branch j′ are considered to be symmetric about a plane which is located at a distance d below the earth surface. In terms of the Neumann’s integral formula [10, 11], the mutual impedance between the branches j and k is expressed by
(3)zjk(ω)=Rjk(ω)+jωLjk(ω)=jωμ04π[∫x1x1+lh∫x1x1+lhdxjdxksjk-∫x1x1+lh∫x1x1+lhdxj′dxksjk′],
where
(4)sjk=(xj-xk)2+a2,sjk′=(xk-xj′)2+(2y1+a+2d)2.
Evaluating the double integrals in (3) leads to
(5)zjk(ω)=jωμ02π×[lhsinh-1lha-a2+lh2z+a-lhsinh-1lh2y1+a+2d+lh2+(2y1+a+2d)2-(2y1+a+2d){lhsinh-1lha-a2+lh2z}].
The self-impedance Zjj(ω) of the transverse branch j is obtained by substituting a=r0 into (5). In a similar way, the mutual impedance between two longitudinal branches, as shown in Figure 4, is given by
(6)Zmn(ω)=jωμ02π×[(2y1+lv+2d)sinh-12y1+lv+2db+(2y1+lv+2d)2+b2-(y1+lv+d)×sinh-12(y1+lv+d)b-124(y1+lv+d)2+b2-(y1+d)sinh-12(y1+d)b+124(y1+d)2+b2+lvsinh-1lvb-b2+lv2+b{(2y1+lv+2d)sinh-12y1+lv+2db}].
In (6), letting b=r0 gives the self-impedance Zmm(ω) of the longitudinal branch m.Figure 3
Two parallel transverse branches.Figure 4
Two parallel longitudinal branches.After the self and mutual impedances are calculated for a coupled branch unit withN transverse or longitudinal branches by using (5) or (6), the impedance matrix Z can be formed as
(7)Z=[Rjk(ω)]N×N+jω[Ljk(ω)]N×N=[zjk(ω)]N×N
or
(8)Z=[Rmn(ω)]N×N+jω[Lmn(ω)]N×N=[zmn(ω)]N×N.
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 d. If the earth is assumed to be perfectly conducting (ρ=0), the complex depth d is zero according to (2). Substitution of d=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
(9)Z0=jωL0,
where L0 is the constant inductance matrix which is independent on the angular frequency ω.


## 3.2. 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 Z0=jωL0 and the admittance matrix Y=jωC of the coupled branches becomes diagonal with negative elements determined from space permeability μ0 and permittivity ε0 [10] as
(10)Z0·Y=(jωL0)·(jωC)=-(ω2ε0μ0)E,
where C is capacitance matrix of the coupled branches and E is the unit matrix. Consequently, the capacitance matrix can be obtained as
(11)C=ε0μ0Z0.


## 4. 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 Figure5, where the earthing system is simply represented as the earthing resistances (Rg). 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) and 6(b) show the π-circuit representing three coupled branches. If the circuit parameters are expressed in matrix form, Figure 6(b) can be further depicted as a compact π-circuit, as shown in Figure 6(c).Figure 5
Injection of lightning current to tower body.π-circuit of three branches: (a) three coupled branches, (b) coupled π-circuit, and (c) compact π-circuit expressed in matrix.
(a)
(b)
(c)With all the coupled branch units in the multiconductor grid replaced by theπ-circuits, Figure 5 is 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 circuit. The lightning transient responses in different parts on the tower body can be obtained by performing a transient calculation for the equivalent circuit.


## 5. 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(12)A=[AC,Aim,AS],
and the relevant current vector is
(13)I(ω)=[IC(ω),Iim(ω),IS(ω)]T,
where the subscripts C, im, and S 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
(14)AI(ω)=ACIC(ω)+AimIim(ω)+ASIS(ω)=0.
Substituting ACYCACTUn(ω) for ACIC(ω) in (14) gives
(15)ACYCACTUn(ω)+AimIim(ω)+ASIS(ω)=0,
where YC is the shunt capacitive admittance matrix. Let Un(ω) and Zim be the node voltage vector and series impedance matrix, respectively. The branch equation of the series impedance is given by
(16)AimTUn(ω)-ZimIim(ω)=0.
Hence, the hybrid equation of the equivalent circuit can be built by merging (15) and (16) as
(17)[ACYCACTAimAimT-Zim][Un(ω)Iim(ω)]=[-ASIS(ω)0].
On the right hand of (17), lightning current source IS(ω) is the Fourier transform of the injected lightning current i (see Figure 5). i is usually expressed by the double exponential function [14]
(18)i=I(e-αt-e-βt).
In accordance with the Chinese national standard [2], the waveform parameter of i is 2.6/50 μs. This gives α=1.5×10-2 μs−1and β=1.86 μs−1 by using the least square fitting method [15]. The peak current I takes the value of 100 kA in light of the severe lightning stroke condition. Thus, finding the Fourier transform of (18) gives
(19)IS(ω)=(β-α)I(α+jω)(β+jω).
By taking the IS(ω) as the excitation, the frequency responses Un(ω) and Iim(ω) 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
(20)ΛX(ω)=b,
where
(21)Λ=[Λjk]M×M=[ACYCACTAimAimT-Zim]M×M,X(ω)=[Xj(ω)]M×1=[Un(ω)Iim(ω)]M×1,b=[bj]M×1=[-ASIS(ω)0]M×1,
where M 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 M-1 step elimination operations, the original matrix Λ in (20) is converted into an upper triangular matrix(22)[Λ11(1)Λ12(1)⋯Λ1k(1)⋯Λ1M(1)0Λ22(2)⋯Λ2k(2)⋯Λ2M(2)⋮⋮⋯⋮⋯⋮00⋯0Λkk(k)⋯ΛkM(k)⋮⋮⋯⋮⋯⋮00⋯0⋯0ΛMM(M)][X1(ω)X2(ω)⋮Xk(ω)⋮XM(ω)]=[b1(1)b2(2)⋮bk(k)⋮bM(M)].The detailed elimination procedure has been stated in [16, 17]. Subsequently, the solution of vector X(ω) is obtained by the back substitution operations
(23)XM(ω)=bM(M)ΛMM(M),Xk(ω)=(bk(k)-∑j=k+1MΛjk(k))Λkk(k)(k=M-1,M-2,…,2,1).
An element of vector X(ω), Xj(ω)  (j=1,2,…,M), is the frequency response and needs to be transformed into time-domain to give the corresponding time response. The inverse 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,
(24)γΔω,γ2Δω,γ3Δω,…,γjΔω,…,γλmΔω,
where Δω(=2πΔf) is the angular frequency step, λm 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 fm(=γλmΔω/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 λm is evaluated by
(25)λm=ln(fm/Δf)lnγ.
Therefore, for a frequency response Xj(ω), its time response xj(t) can be obtained by
(26)xj(t)=1π∑k=1λmXj[(γ2k-1-1)Δω]×exp[j(γ2k-1-1)Δωt]σγ2(k-1)(γ2-1)Δω,
where σ is the standard smoothing factor
(27)σ=sin[(λ2k-1-1)π/(γ2λm-1)](γ2k-1-1)π/(γ2λm-1).
Owing to a high accuracy and long observation time with a small number of sampling steps, (26) is suitable for the lightning transient calculation.


## 6. 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 andh=82m, r1=2.7 m, r2=4.3 m, and w=0.025 m (see Figure 1(a)). The relevant discrete multiconductor grid is shown in Figure 7, where r0=3.5 m. The earthing resistance Rg is 3.5Ω. The transient responses on the tower body are calculated by using the proposed method. The waveforms of current on branch P2,1-P3,1 and potential at node P1,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.Figure 7
Multiconductor grid of a 1.5 MW WT tower.Current and potential waveforms: (a) current waveform on branchP2,1-P3,1 and (b) potential waveform of node P1,1.
(a)
(b)Figure 9
Peak potentials at different nodes on the tower body.Research Centre for electromagnetic transient simulation and widely applied in electric power systems. It is clear from Figure9 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.


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



---

*Source: 757656-2013-04-09.xml*