Transient Simulation of Wind Turbine Towers under Lightning Stroke

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.

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 Figure 1 [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 (𝑅_𝑔). 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 𝑗 and 𝑘 is calculated by 𝐿_𝑗𝑘=𝜇₀4𝜋_{𝑠𝑗𝑠𝑘}𝑑𝐬_𝑗⋅𝑑𝐬_𝑘𝐃_𝑗𝑘+_{𝑠′𝑗𝑠𝑘}𝑑𝐬^′_𝑗⋅𝑑𝐬_𝑘𝐃^′_𝑗𝑘. (1)

(a)

(b)

Figure 3 

Branches in the multiconductor system.

The dot products of the vector differential segments become 𝑑𝐬_𝑗⋅𝑑𝐬_𝑘=cos𝜑𝑑𝐬_𝑗𝑑𝐬_𝑘,𝑑𝐬^′_𝑗⋅𝑑𝐬_𝑘=cos𝜑^′𝑑𝑠^′_𝑗𝑑𝑠, (2) where 𝜑 is the direction angle between the branches𝑗 and 𝑘, and 𝜑^′ is that between the branches 𝑗^′ and 𝑘. Substituting (2) into (1) gives 𝐿_𝑗𝑘=𝜇₀4𝜋cos𝜑𝑁𝑗,𝑘+cos𝜑^′𝑁𝑗^′,𝑘, (3) where 𝑁(𝑗,𝑘) and 𝑁(𝑗^′,𝑘) are two double line integrals: 𝑁𝑗,𝑘=_{𝑠𝑗𝑠𝑘}𝑑𝑠_𝑗𝑑𝑠_𝑘𝐷_𝑗𝑘,𝑁𝑗^′,𝑘=_{𝑠′𝑗𝑠𝑘}𝑑𝑠^′_𝑗𝑑𝑠_𝑘𝐷^′_𝑗𝑘. (4) Equations (3) and (4) indicate that the calculation of the mutual inductance depends on the double line integrals 𝑁(𝑗,𝑘) and 𝑁(𝑗^′,𝑘). From the viewpoint of integral operation, the way to evaluate 𝑁(𝑗,𝑘) is the same as that to evaluate 𝑁(𝑗^′,𝑘). For this reason, the solution of 𝑁(𝑗,𝑘) is only presented for the typical space positions in the multiconductor system, which holds for the evaluation of 𝑁(𝑗^′,𝑘).

The integral 𝑁(𝑗,𝑘) 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 𝑁(𝑗,𝑘) is given by 𝑁𝑗,𝑘=𝑙_𝑗+𝑠_𝑗𝐷ln²+𝐷₄+𝑠_𝑘𝐷₂+𝐷₄−𝑠_𝑘−𝑙_𝑗𝐷ln¹+𝐷₃+𝑠_𝑘𝐷₁+𝐷₃−𝑠_𝑘+𝑙_𝑘+𝑠_𝑘𝐷ln³+𝐷₄+𝑠_𝑗𝐷₃+𝐷₄−𝑠_𝑗−𝑙_𝑘𝐷ln¹+𝐷₂+𝑠_𝑗𝐷₁+𝐷₂−𝑠_𝑗. (5) If the two longitudinal branches are flush with each other, the integral 𝑁(𝑗,𝑘) is evaluated by putting 𝑠_𝑗=𝑠_𝑘 and 𝐷₂=𝐷₃ into (5). Under the flush condition, the integral 𝑁(𝑗,𝑗) is evaluated by again putting 𝐷₁=𝐷₄=𝑟₀ into (5). This can result in determining the self-inductance 𝐿_𝑗𝑗 of the longitudinal branch 𝑗 by (3).

(a)

(b)

Figure 4 

Coplanar branch pairs, (a) longitudinal pair, and (b) transverse pair.

For the two coplanar transverse branches (𝜑=0), as shown in Figure 4(b), the integral 𝑁(𝑗,𝑘) is given by 𝑁𝑠𝑗,𝑘=−_𝑘−𝑠_𝑗sinh⁻¹𝑠_𝑘−𝑠_𝑗+𝑠2ℎ^𝑗+𝑠_𝑘sinh⁻¹𝑠_𝑗+𝑠_𝑘2ℎ+2𝑠_𝑘−𝑠_𝑗2²+ℎ²−2𝑠_𝑗+𝑠_𝑘2²+ℎ². (6) In (6), letting 𝑠_𝑗=𝑠_𝑘 and ℎ=𝑟₀ gives the integral 𝑁(𝑗,𝑗). Then, 𝑁(𝑗,𝑗) can be used to further determine the self-inductance 𝐿_𝑗𝑗 of the transverse branch𝑗 by (3).

The noncoplanar branch pairs are shown in Figure 5, where line OG is the common perpendicular to the two branches and line OF lying on the plane containing the branch 𝑗 is parallel to the branch 𝑘. The integral 𝑁(𝑗,𝑘) for both the longitudinal branch pair (Figure 5(a)) and the transverse branch pair (Figure 5(b)) is given in the following complicated expression 𝑁𝑙𝑗,𝑘=2_𝑗+𝑠_𝑗tanh⁻¹𝑠_𝑘𝐷₁+𝐷₂−𝑙_𝑗tanh⁻¹𝑠_𝑘𝐷₃+𝐷₄+𝑙_𝑘+𝑠_𝑘tanh⁻¹𝑠_𝑗𝐷₁+𝐷₄−𝑙_𝑘tanh⁻¹𝑠_𝑗𝐷₂+𝐷₃−Φℎ,sin𝜑 (7) where Φ is the solid angle Φ=tanh⁻¹ℎ²cos𝜑+𝑙_𝑗+𝑠_𝑗𝑙_𝑘+𝑠_𝑘sin²𝜑ℎ𝐷₁sin𝜑−tanh⁻¹ℎ²cos𝜑+𝑙_𝑘𝑙_𝑗+𝑠_𝑗sin²𝜑ℎ𝐷₂sin𝜑+tanh⁻¹ℎ²cos𝜑+𝑙_𝑗𝑙_𝑘sin²𝜑ℎ𝐷₃sin𝜑−tanh⁻¹ℎ²cos𝜑+𝑙_𝑗𝑙_𝑘+𝑠_𝑘sin²𝜑ℎ𝐷₄.sin𝜑 (8)

(a)

(b)

Figure 5 

Noncoplanar branch pairs, (a) longitudinal pair, and (b) transverse pair.

On the basis of the formulas given previously, the self and mutual inductances can be calculated for a branch unit with 𝑀 coupled branches in the multiconductor system, thus forming the inductance matrix 𝐋=𝐿_{𝑗𝑘𝑀×𝑀}. (9)

The branch resistance per unit length is approximately calculated by [8] 𝜌𝑅=𝜋𝜎1−exp−𝑟₀/𝜎2𝑟₀−𝜎1−exp−𝑟₀,/𝜎 (10) where 𝑟₀ is the radius of the branch, 𝜌 the material resistivity, and 𝜎 the skin depth: 1𝜎=,𝜋𝑓𝜇𝛾 (11) where 𝑓 is the maximum frequency likely to affect the system transient, 𝜇 the material permeability, and 𝛾=1/𝜌. 𝑓 may be roughly evaluated by the waveform parameter of the injected lightning current [9].

2.3. Capacitance

The mutual potential coefficient between the branches 𝑗 and 𝑘, as shown in Figure 3, can be calculated by the average potential method [10]: 𝑝_𝑗𝑘=14𝜋𝜀₀𝑠_𝑗𝑠_{𝑘𝑠𝑗𝑠𝑘}𝑑𝑠_𝑗𝑑𝑠_𝑘𝐷_𝑗𝑘−_{𝑠′𝑗𝑠𝑘}𝑑𝑠^′_𝑗𝑑𝑠_𝑘𝐷^′_𝑗𝑘=14𝜋𝜀₀𝑠_𝑗𝑠_𝑘𝑁𝑗,𝑘−𝑁𝑗^′,𝑘. (12) From (12), the self-potential coefficient 𝑝_𝑗𝑗 of the branch 𝑗 can also be obtained in a similar way to the self-inductance 𝐿_𝑗𝑗. Consequently, the potential coefficient matrix of a branch unit with 𝑀 coupled branches in the multiconductor system is formed by 𝐏=𝑝_{𝑗𝑘𝑀×𝑀}. (13) Calculating the inverse of 𝐏 gives 𝐏⁻¹=𝑞_{𝑗𝑘𝑀×𝑀}. (14) From the inverse matrix 𝐏⁻¹, the capacitance matrix can be obtained as [10, 11] 𝐂=𝐶_{𝑗𝑘𝑀×𝑀}, (15) where the diagonal and off-diagonal elements are determined as 𝐶_𝑗𝑗=^𝑀_𝜁=1𝑞_𝑗𝜁,𝐶𝑗=1,2,…,𝑀_𝑗𝑘=−𝑞_𝑗𝑘𝑗,𝑘=1,2,…,𝑀,𝑗≠𝑘. (16)

3. Simulation of Lightning Transients

By means of the electrical parameters obtained previously, a branch unit with 𝑀 coupled branches is represented by a coupled 𝜋-type circuit formed by resistances, inductances, and capacitances. A coupled 𝜋-type circuit representing three coupled branches (𝑀=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 Δ𝑡 is the time step size, 𝐈_𝑅𝐿(𝑡−Δ𝑡) and 𝐈_𝐶(𝑡−Δ𝑡) are the current source vectors determined by the current and voltage values from previous time steps, and 𝐑_𝑅𝐿 and 𝐑_𝐶 are the equivalent resistance matrices determined by the electrical parameters and the time step size. The formulas for calculating 𝐈_𝑅𝐿(𝑡−Δ𝑡), 𝐈_𝐶(𝑡−Δ𝑡), 𝐑_𝑅𝐿, and 𝐑_𝐶 have been given in [13, 14]. For 𝑀 (𝑀>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.

(a)

(b)

(c)

(d)

Figure 6 

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.

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 Figure 7. The reduced-scale tower is made of the iron sheet, whose dimensions are taken as 𝑟₁=0.0283 m, 𝑟₂=0.067 m, ℎ=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.

(a)

(b)

Figure 8 

Measured and simulated waveforms: (a) injected current waveform and (b) potential waveform at the top of the tower.

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 are 𝑟₁=1.35 m, 𝑟₂=2.17 m, ℎ=82 m, and 𝛿=0.025 m, and the grounding resistance 𝑅_𝑔 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.

(a)

(b)

Figure 9 

Transient potential waveforms: (a) potential waveforms at the top of the tower and (b) potential waveforms at the bottom of the tower.

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 Figure 11. 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.

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China under Award no. 509770926 and the Fundamental Research Funds for the Central Universities under Award no. 2012JBZ006. The author expresses his thanks to the foundation committees.