# Transient Simulation of Wind Turbine Towers under Lightning Stroke

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

---

## Abstract

A simulation algorithm is proposed in this paper for lightning transient analysis of the wind turbine (WT) towers. In the proposed algorithm, the tower body is first subdivided into a discrete multiconductor system. A set of formulas are given to calculate the electrical parameters of the branches in the multiconductor system. By means of the electrical parameters, each branch unit in the multiconductor system is replaced as a coupledπ-type circuit and the multiconductor system is converted into a circuit model. Then, the lightning transient responses can be obtained in different parts on the tower body by solving the circuit equations of the equivalent discretization network. The laboratory measurement is also made by a reduced-scale tower for checking the validity of the proposed algorithm.


---

## Body

## 1. Introduction

Global warming effect accelerates the utilization of wind energy. As a clean energy source, wind energy can be used to generate electric power without emission of carbon dioxide into the atmospheric environment. In consequence of the rapid growth in the utilization of wind energy for electric power supply, wind turbines (WTs) have increased constantly in size and rated power during recent decades. WTs are particularly vulnerable to lightning strokes due to their great height, distinctive shape, and rather exposed position. The lightning stroke effect on WTs has become a major concern as the number of the installed WTs continues to increase. Therefore, lightning protection of WTs is crucially important for the operational reliability of large wind power generation systems. The lightning protection design needs to obtain the lightning transient responses on WT towers since the tower body is the main conducting path of lightning current. The simulation algorithms were presented in the literature [1–3], in which the tower body was simply represented by a transmission line or a chain capacitance circuit. However, these existing algorithms are difficult to use for calculating the transient responses in different parts on the tower body due to the fact that they ignore the structural feature of the tower body. To overcome the shortcoming in these existing algorithms, a novel simulation algorithm is proposed in this paper. The proposed algorithm subdivides the continuous conducting shell of the tower body into a discrete multiconductor system and can give due consideration for the structural feature of the tower body. In the multiconductor system, the electrical parameters of the branches are calculated by an efficient procedure, and the branch units are represented by the coupled π-type circuits. A circuit model, constituted by a series of the coupled π-type circuits, is built for the multiconductor system. On the basis of the circuit model, the equivalent discretization network is further formed. The transient calculation is performed for the equivalent discretization network, and then the lightning transient responses can be obtained in different parts on the tower body. The measurement of the transient response is also taken on a reduced-scale WT tower, and the validity of the proposed algorithm is examined by comparing simulated and measured results.


## 2. Calculation of Electrical Parameters

### 2.1. Discretization Description of the WT Tower

After a WT is struck by lightning, lightning current usually passes from the blade root to the tower and then dissipates in the ground through the earth-termination system. Since most manufactures have used many kinds of  brushes and sliding contact systems to divert the lightning current from the main shaft, only a much smaller part of lightning current flows through the gearbox and the generator. For the WT structure, the tower is the longest conducting path of lightning current. When the tower conducts lightning current, serious transient potential rise appears on the tower body and may result in flashback to the electrical and electronic components installed inside the tower. Much attention to the potential rise has been paid in the lightning protection design of WTs. For the sake of the transient analysis of the potential distribution on the tower, the blade and the conducting path in the nacelle are left out of consideration. Instead, lightning current is injected into the tower body from its top, as shown in Figure1 [4]. In order to build the circuit model of the tower, the continuous conducting shell of the tower body needs to be subdivided into a discrete multiconductor system formed by longitudinal and transverse branches, as illustrated in Figure 2 [5]. The WT earth-termination system is modeled as the grounding resistances (Rg). The electrical parameters of the branches in the multiconductor system are represented by resistances, inductances, and capacitances. For simplifying the parameter calculation, the segmental arc of each transverse branch is replaced by its chord and all the branches are taken as the cylindrical conductors. The conductor radii are estimated from the average cross-sections of the respective branches.Figure 1
WT tower.Figure 2
Multiconductor system.


### 2.2. Inductance and Resistance

The inductance parameters can be determined by Neumann’s integral method [6, 7]. The typical branches in the multiconductor system are illustrated in Figure 3, where the presence of the ground is considered by symmetrically installing the image branches (depicted in dotted lines) below the earth surface [7]. According to Neumann’s integral method, the mutual inductance between the branches j and k is calculated by
(1)Ljk=μ04π[∫sj∫skdsj·dskDjk+∫sj′∫skdsj′·dskDjk′].Branches in the multiconductor system.
(a)
(b)The dot products of the vector differential segments become(2)dsj·dsk=cosφdsjdsk,dsj′·dsk=cosφ′dsj′ds,
where φ is the direction angle between the branchesj and k, and φ′ is that between the branches j′ and k. Substituting (2) into (1) gives
(3)Ljk=μ04π[cosφN(j,k)+cosφ′N(j′,k)],
where N(j,k) and N(j′,k) are two double line integrals:
(4)N(j,k)=∫sj∫skdsjdskDjk,N(j′,k)=∫sj′∫skdsj′dskDjk′.
Equations (3) and (4) indicate that the calculation of the mutual inductance depends on the double line integrals N(j,k) and N(j′,k). From the viewpoint of integral operation, the way to evaluate N(j,k) is the same as that to evaluate N(j′,k). For this reason, the solution of N(j,k) is only presented for the typical space positions in the multiconductor system, which holds for the evaluation of N(j′,k).The integralN(j,k) is first evaluated in the case of the coplanar branch pairs. For the two coplanar longitudinal branches, as shown in Figure 4(a), the integral N(j,k) is given by
(5)N(j,k)=(lj+sj)lnD2+D4+skD2+D4-sk-ljlnD1+D3+skD1+D3-sk+(lk+sk)lnD3+D4+sjD3+D4-sj-lklnD1+D2+sjD1+D2-sj.
If the two longitudinal branches are flush with each other, the integral N(j,k) is evaluated by putting sj=sk and D2=D3 into (5). Under the flush condition, the integral N(j,j) is evaluated by again putting D1=D4=r0 into (5). This can result in determining the self-inductance Ljj of the longitudinal branch j by (3).Coplanar branch pairs, (a) longitudinal pair, and (b) transverse pair.
(a)
(b)For the two coplanar transverse branches (φ=0), as shown in Figure 4(b), the integral N(j,k) is given by
(6)N(j,k)=-(sk-sj)sinh-1sk-sj2h+(sj+sk)sinh-1sj+sk2h+2(sk-sj2)2+h2-2(sj+sk2)2+h2.
In (6), letting sj=sk and h=r0 gives the integral N(j,j). Then, N(j,j) can be used to further determine the self-inductance Ljj of the transverse branchj by (3).The noncoplanar branch pairs are shown in Figure5, where line OG is the common perpendicular to the two branches and line OF lying on the plane containing the branch j is parallel to the branch k. The integral N(j,k) for both the longitudinal branch pair (Figure 5(a)) and the transverse branch pair (Figure 5(b)) is given in the following complicated expression
(7)N(j,k)=2[(lj+sj)tanh-1skD1+D2-ljtanh-1skD3+D4+(lk+sk)tanh-1sjD1+D4-lktanh-1sjD2+D3]-Φhsinφ,
where Φ is the solid angle
(8)Φ=tanh-1[h2cosφ+(lj+sj)(lk+sk)sin2φhD1sinφ]-tanh-1[h2cosφ+lk(lj+sj)sin2φhD2sinφ]+tanh-1[h2cosφ+ljlksin2φhD3sinφ]-tanh-1[h2cosφ+lj(lk+sk)sin2φhD4sinφ].Noncoplanar branch pairs, (a) longitudinal pair, and (b) transverse pair.
(a)
(b)On the basis of the formulas given previously, the self and mutual inductances can be calculated for a branch unit withM coupled branches in the multiconductor system, thus forming the inductance matrix
(9)L={Ljk}M×M.The branch resistance per unit length is approximately calculated by [8]
(10)R=ρπσ[(1-exp(-r0/σ))][2r0-σ(1-exp(-r0/σ))],
where r0 is the radius of the branch, ρ the material resistivity, and σ the skin depth:
(11)σ=1πfμγ,
where f is the maximum frequency likely to affect the system transient, μ the material permeability, and γ=1/ρ. f may be roughly evaluated by the waveform parameter of the injected lightning current [9].


### 2.3. Capacitance

The mutual potential coefficient between the branchesj and k, as shown in Figure 3, can be calculated by the average potential method [10]:
(12)pjk=14πε0sjsk[∫sj∫skdsjdskDjk-∫sj′∫skdsj′dskDjk′]=14πε0sjsk[N(j,k)-N(j′,k)].
From (12), the self-potential coefficient pjj of the branch j can also be obtained in a similar way to the self-inductance Ljj. Consequently, the potential coefficient matrix of a branch unit with M coupled branches in the multiconductor system is formed by
(13)P={pjk}M×M.
Calculating the inverse of P gives
(14)P-1={qjk}M×M.
From the inverse matrix P-1, the capacitance matrix can be obtained as [10, 11]
(15)C={Cjk}M×M,
where the diagonal and off-diagonal elements are determined as
(16)Cjj=∑ζ=1Mqjζ(j=1,2,…,M),Cjk=-qjk(j,k=1,2,…,M,j≠k).


## 2.1. Discretization Description of the WT Tower

After a WT is struck by lightning, lightning current usually passes from the blade root to the tower and then dissipates in the ground through the earth-termination system. Since most manufactures have used many kinds of  brushes and sliding contact systems to divert the lightning current from the main shaft, only a much smaller part of lightning current flows through the gearbox and the generator. For the WT structure, the tower is the longest conducting path of lightning current. When the tower conducts lightning current, serious transient potential rise appears on the tower body and may result in flashback to the electrical and electronic components installed inside the tower. Much attention to the potential rise has been paid in the lightning protection design of WTs. For the sake of the transient analysis of the potential distribution on the tower, the blade and the conducting path in the nacelle are left out of consideration. Instead, lightning current is injected into the tower body from its top, as shown in Figure1 [4]. In order to build the circuit model of the tower, the continuous conducting shell of the tower body needs to be subdivided into a discrete multiconductor system formed by longitudinal and transverse branches, as illustrated in Figure 2 [5]. The WT earth-termination system is modeled as the grounding resistances (Rg). The electrical parameters of the branches in the multiconductor system are represented by resistances, inductances, and capacitances. For simplifying the parameter calculation, the segmental arc of each transverse branch is replaced by its chord and all the branches are taken as the cylindrical conductors. The conductor radii are estimated from the average cross-sections of the respective branches.Figure 1
WT tower.Figure 2
Multiconductor system.


## 2.2. Inductance and Resistance

The inductance parameters can be determined by Neumann’s integral method [6, 7]. The typical branches in the multiconductor system are illustrated in Figure 3, where the presence of the ground is considered by symmetrically installing the image branches (depicted in dotted lines) below the earth surface [7]. According to Neumann’s integral method, the mutual inductance between the branches j and k is calculated by
(1)Ljk=μ04π[∫sj∫skdsj·dskDjk+∫sj′∫skdsj′·dskDjk′].Branches in the multiconductor system.
(a)
(b)The dot products of the vector differential segments become(2)dsj·dsk=cosφdsjdsk,dsj′·dsk=cosφ′dsj′ds,
where φ is the direction angle between the branchesj and k, and φ′ is that between the branches j′ and k. Substituting (2) into (1) gives
(3)Ljk=μ04π[cosφN(j,k)+cosφ′N(j′,k)],
where N(j,k) and N(j′,k) are two double line integrals:
(4)N(j,k)=∫sj∫skdsjdskDjk,N(j′,k)=∫sj′∫skdsj′dskDjk′.
Equations (3) and (4) indicate that the calculation of the mutual inductance depends on the double line integrals N(j,k) and N(j′,k). From the viewpoint of integral operation, the way to evaluate N(j,k) is the same as that to evaluate N(j′,k). For this reason, the solution of N(j,k) is only presented for the typical space positions in the multiconductor system, which holds for the evaluation of N(j′,k).The integralN(j,k) is first evaluated in the case of the coplanar branch pairs. For the two coplanar longitudinal branches, as shown in Figure 4(a), the integral N(j,k) is given by
(5)N(j,k)=(lj+sj)lnD2+D4+skD2+D4-sk-ljlnD1+D3+skD1+D3-sk+(lk+sk)lnD3+D4+sjD3+D4-sj-lklnD1+D2+sjD1+D2-sj.
If the two longitudinal branches are flush with each other, the integral N(j,k) is evaluated by putting sj=sk and D2=D3 into (5). Under the flush condition, the integral N(j,j) is evaluated by again putting D1=D4=r0 into (5). This can result in determining the self-inductance Ljj of the longitudinal branch j by (3).Coplanar branch pairs, (a) longitudinal pair, and (b) transverse pair.
(a)
(b)For the two coplanar transverse branches (φ=0), as shown in Figure 4(b), the integral N(j,k) is given by
(6)N(j,k)=-(sk-sj)sinh-1sk-sj2h+(sj+sk)sinh-1sj+sk2h+2(sk-sj2)2+h2-2(sj+sk2)2+h2.
In (6), letting sj=sk and h=r0 gives the integral N(j,j). Then, N(j,j) can be used to further determine the self-inductance Ljj of the transverse branchj by (3).The noncoplanar branch pairs are shown in Figure5, where line OG is the common perpendicular to the two branches and line OF lying on the plane containing the branch j is parallel to the branch k. The integral N(j,k) for both the longitudinal branch pair (Figure 5(a)) and the transverse branch pair (Figure 5(b)) is given in the following complicated expression
(7)N(j,k)=2[(lj+sj)tanh-1skD1+D2-ljtanh-1skD3+D4+(lk+sk)tanh-1sjD1+D4-lktanh-1sjD2+D3]-Φhsinφ,
where Φ is the solid angle
(8)Φ=tanh-1[h2cosφ+(lj+sj)(lk+sk)sin2φhD1sinφ]-tanh-1[h2cosφ+lk(lj+sj)sin2φhD2sinφ]+tanh-1[h2cosφ+ljlksin2φhD3sinφ]-tanh-1[h2cosφ+lj(lk+sk)sin2φhD4sinφ].Noncoplanar branch pairs, (a) longitudinal pair, and (b) transverse pair.
(a)
(b)On the basis of the formulas given previously, the self and mutual inductances can be calculated for a branch unit withM coupled branches in the multiconductor system, thus forming the inductance matrix
(9)L={Ljk}M×M.The branch resistance per unit length is approximately calculated by [8]
(10)R=ρπσ[(1-exp(-r0/σ))][2r0-σ(1-exp(-r0/σ))],
where r0 is the radius of the branch, ρ the material resistivity, and σ the skin depth:
(11)σ=1πfμγ,
where f is the maximum frequency likely to affect the system transient, μ the material permeability, and γ=1/ρ. f may be roughly evaluated by the waveform parameter of the injected lightning current [9].


## 2.3. Capacitance

The mutual potential coefficient between the branchesj and k, as shown in Figure 3, can be calculated by the average potential method [10]:
(12)pjk=14πε0sjsk[∫sj∫skdsjdskDjk-∫sj′∫skdsj′dskDjk′]=14πε0sjsk[N(j,k)-N(j′,k)].
From (12), the self-potential coefficient pjj of the branch j can also be obtained in a similar way to the self-inductance Ljj. Consequently, the potential coefficient matrix of a branch unit with M coupled branches in the multiconductor system is formed by
(13)P={pjk}M×M.
Calculating the inverse of P gives
(14)P-1={qjk}M×M.
From the inverse matrix P-1, the capacitance matrix can be obtained as [10, 11]
(15)C={Cjk}M×M,
where the diagonal and off-diagonal elements are determined as
(16)Cjj=∑ζ=1Mqjζ(j=1,2,…,M),Cjk=-qjk(j,k=1,2,…,M,j≠k).


## 3. Simulation of Lightning Transients

By means of the electrical parameters obtained previously, a branch unit withM coupled branches is represented by a coupled π-type circuit formed by resistances, inductances, and capacitances. A coupled π-type circuit representing three coupled branches (M=3) is shown in Figures 6(a) and 6(b). Expression of the electrical parameters in matrix form gives a simplified circuit, as shown in Figure 6(c). With a time discretization treatment for the inductance and capacitance elements, these two kinds of circuit elements are replaced by current sources and parallel equivalent resistances [12, 13]. Thus, Figure 6(c) is further converted into a discretization π-type circuit, as shown in Figure 6(d), where Δt is the time step size, IRL(t-Δt) and IC(t-Δt) are the current source vectors determined by the current and voltage values from previous time steps, and RRL and RC are the equivalent resistance matrices determined by the electrical parameters and the time step size. The formulas for calculating IRL(t-Δt), IC(t-Δt), RRL, and RC have been given in [13, 14]. For M (M>3) coupled branches, the discretization π-type circuit is similar to that in Figure 6(d). After all the branch units in the multiconductor system are replaced by the discretization π-type circuits, the tower body can be converted into an equivalent discretization network. The lightning current source is injected to the top nodes of the equivalent discretization network. According to the topology connection mode of the equivalent resistances and current sources, the node voltage equations are formed and then solved numerically at each time step. The solving process has been described in detail in [13, 14]. As a result, lightning transient responses can be obtained in different parts on the tower body.Coupledπ-type circuit of three coupled branches: (a) three coupled branches, (b) coupled π-type circuit, (c) simplified π-type circuit expressed in matrix, and (d) discretization π-type circuit.
(a)
(b)
(c)
(d)If the blade and the conducting path in the nacelle are taken into account, the former can be converted into a series circuit unit consisting of a fewπ-type circuits and the latter into a parallel resistance-capacitance unit. In the parallel unit, the resistance is the contact one of the brush and the sliding contact system, while the capacitance represents the main shaft bearings [2, 15]. After the two circuit units are serially connected to the equivalent discretization network of the tower, the complete circuit model can be built for the WT. However, the two circuit units in the complete circuit model do not change the lightning current and transient potential distributions on the tower body due to their series connection to the tower. Thus, their removal from the complete circuit model will not make a significant influence on calculation of the potential responses on the tower body.


## 4. Verification of the Proposed Algorithm with Laboratory Measurement

An experimental setup was built in the laboratory space, as shown in Figure7. The reduced-scale tower is made of the iron sheet, whose dimensions are taken as r1=0.0283 m, r2=0.067 m, h=2 m, and δ=0.002 m (see Figure 1). The potential measurement wire and the current lead wire are orthogonal to weaken the electromagnetic induction between them.Figure 7
Experimental setup.The potential measurement wire is grounded, that is, connected to the steel plate, at a point 9 m apart from the tower. Five resistances of 5Ω, representing the grounding resistances, are connected between the tower bottom and steel plate. The impulse current with fast wavefront is injected to the top of the tower. Under the excitation of the injected current, the transient potential response is measured at the top of the tower. The frequency bandwidth of the potential probe is 0~500 MHz and that of the current probe 20 kHz~500 MHz, which is competent to measure the fast transient responses. The measured waveforms of the injected current and the transient potential are shown in Figure 8, where the corresponding waveform simulated by the proposed algorithm is given simultaneously for comparison. It can be found that a better agreement appears between measured and simulated waveforms. This agreement confirms the validity of the proposed algorithm.Measured and simulated waveforms: (a) injected current waveform and (b) potential waveform at the top of the tower.
(a)
(b)


## 5. Application Example

This example takes into account an actual WT with a rated power of 2.5 MW. The dimensions of the WT tower arer1=1.35 m, r2=2.17 m, h=82 m, and δ=0.025 m, and the grounding resistance Rg is 4 Ω (see Figure 1). The parameter of the injected lightning current is taken as 10/350 μs, 100 kA according to the technical specification of lightning protection [16]. By using the proposed algorithm to perform transient simulation, the transient potential responses are obtained on the tower body. The potential waveforms at the top and the bottom of the tower are shown in Figure 9, and the peak potential distribution along the height of the tower is also plotted in Figure 10. As seen in Figures 9 and 10, the transient potential rise on the tower body is extremely serious and may do damage to the facilities installed inside the tower under lightning stroke.Transient potential waveforms: (a) potential waveforms at the top of the tower and (b) potential waveforms at the bottom of the tower.
(a)
(b)Figure 10
Peak potential distribution along the height of the tower.In an actual wind farm, individual WT earth-termination systems are usually connected by the metallic armor of the power cable running between the WTs. An interconnected grounding system consisting of 5 WTs and a substation is shown in Figure11. The distance between two successive WTs is 380 m. The substation (SS) is located 390 m away from WT5. Inner and outer radii of the power cable are 0.0135 m and 0.016 m, respectively. The burial depth of the power cable is 1 m, and the grounding grid of the substation has a grounding resistance of 10 Ω. The power cable between two successive WTs or SS and WT5 is divided into 5 segments. Each segment is represented by a π-type circuit [17]. The peak earth potential distribution in the interconnected grounding system is shown in Figure 12, where the peak earth potential at each WT or SS is expressed by the percentage of the maximum earth potential. As expected, the maximum earth potential always appears at WT3, where the lightning strikes. The earth potential rise at WT2 and WT4 is more serious than that at WT1 and WT5, for WT2 and WT4 are adjacent to WT3. This transferred overvoltage is very harmful to the neighboring WT transformers and the power cables.Figure 11
Layout of the interconnected grounding system.Figure 12
Earth potential distribution in the interconnected grounding system.In view of the results obtained above, the efficient protective measures against the transient potential rise on the tower body and the transferred overvoltage in the interconnected grounding system need to be taken for the multimegawatt WTs [18].


## 6. Conclusions

The lightning transient analysis has been carried out in this paper for the WT towers. The subdividing treatment of the large-sized continuous conducting shell allows a WT tower to be converted into a discrete multiconductor system. For calculating the electrical parameters of the branches in the multiconductor system, a set of analytical formulas have been given. These formulas have the capability of considering the electromagnetic coupling between the branches in different space positions. By means of the electrical parameters, the branch units are represented by the coupledπ-type circuits and the equivalent discretization network is further built for the multiconductor system. The solution of lightning transients is then obtained from the equivalent discretization network. The proposed algorithm can take into account the structural feature of actual large-sized WT towers and give the lightning transient responses in different parts on the tower body. The laboratory measurement has been made by a reduced-scale tower. A better agreement appears between simulated and measured potential waveforms. The applicability of the proposed algorithm to the lightning transient analysis has also been examined by an actual example of a 2.5 MW WT tower and its interconnected grounding system.



---

*Source: 142765-2013-08-19.xml*