Transient behavior of wind towers grounding systems under lightning strikes

Wind turbines, because of their height and localization, are frequently subject to direct lightning strikes. Thus, the investigation on the performance of their grounding system is of paramount interest for the prediction of the potential threats either for people working in (or animals passing through) the wind farm area and for the power and control units installed in close proximity of the turbine towers. In this paper,we perform a comprensive study of the transient behavior of the earthing systems of wind turbines. The analysis is conducted in the frequency domain and an hybrid approach, based on circuit theory and Method of Moments, is adopted to fully account for resistive, inductive and capacitive couplings between elements of the ground system. The actual transient behavior is obtained by means of an Inverse Fourier Transform. The results, computed considering a typical wind turbine grounding system configuration, provide more insight on the nature of the early-time transient response of grounding systems and allow to draw up useful design guidelines.

Abstract Wind turbines, because of their height and localization, are frequently subject to direct lightning strikes. Thus, the investigation on the performance of their grounding system is of paramount interest for the prediction of the potential threats either for people working in (or animals passing through) the wind farm area and for the power and control units installed in close proximity of the turbine towers. In this paper,we perform a comprensive study of the transient behavior of the earthing systems of wind turbines. The analysis is conducted in the frequency domain and an hybrid approach, based on circuit theory and Method of Moments, is adopted to fully account for resistive, inductive and capacitive couplings between elements of the ground system. The actual transient behavior is obtained by means of an Inverse Fourier Transform. The results, computed considering a typical wind turbine grounding system configuration, provide more insight on the nature of the early-time transient response of grounding systems and allow to draw up useful design guidelines.

Introduction
Wind farms are an excellent solution for the generation of clean electric power because of their reduced costs for investments and maintenance and their diffusion has been observed worldwide [1][2][3]. Wind turbines energy have been widely considered as one of the best choices in terms of payback time and energy conversion; thus, their installations have recently increased, interesting also those regions where lightning activity is significant and consequent damages may very likely occur [4,5]. Wind turbines are candidate victims for cloud-to-ground lightning mainly due to their special shape and complexity of apparatus: moreover, they are often located in isolated locations, at hill or mountain altitudes, where the high values of the keraunic level increase the probability of damage. Relevant statistics indicates that fatal damages in blades, generators and especially in control circuits are caused either by direct lightning strokes to wind turbines or transferred overvoltages from nearby fault locations [6]. Grounding systems of wind turbines are designed to prevent overvoltages and potential gradients that could damage parts of the wind turbine and to reduce human hazards. However, the design of the grounding system calls for the thorough analysis of the following peculiarities: -the grounding system of a wind turbine is generally much smaller than those of buildings with equal height; -grounding systems of turbines in wind farms are electrically interconnected; -the lightning protection level for a wind turbine is much higher than that of a normal building having an equivalent foundation; -generally poor grounding conditions are observed, i.e., high values of soil resistivity, ranging between 100 and 2000 X m [7,8].
The hazard and the risk are strictly dependent upon either values and duration of the transient phenomenon; thus, the grounding systems behavior of actual wind turbines has to be analyzed in the time domain. In addition, the transient analysis [9,10] is important because the injection of high impulse currents into the grounding system leads to an increase of the grounding potential during the transient state, which can occur in a danger to humans, animals, installations and equipment or in a loss of data transmitted through power or telecommunication cables [11]. The analysis is performed in the frequency domain by means of a recently proposed hybrid approach [12,13] based on circuit theory and Method of Moments, capable to accurately account for resistive, inductive and capacitive couplings. However, in our formulation no approximations are introduced to model the interactions between elements and this improvement yields more accurate results. Then, the corresponding transient response is obtained by means of the inverse Fourier transform (IFT). Other efficient approaches have been proposed in the past, based either on analytical approximations or on numerical methods, like FD-TD, e.g. [14][15][16][17]. In the proposed method, solving accurately the Sommerfeld's integrals occurring in the Green's functions of a stratified medium, appears to be very accurate and reliable especially for the early-time response, which is responsible of the most severe solicitation for electrical and electronic apparatus and systems installed in close proximity of the wind towers.

Grounding system analysis
The system configuration is shown in Fig. 1: a wind turbine grounding system is buried and subjected to a transient current generated by the lighting strike at a certain point of blades or nacelle (the influence of the turbine tower is here neglected, since the body is in series with the considered current source). Generally, grounding systems configurations are prescribed or recommended by manufacturers: typical grounding systems are arranged in a ring shape around the tower's base and connected with the tower itself through its irons in the foundation. Low-impedance grounding system is a major prerequisite for an effective protection of wind turbines from lightning strikes. In the following, a time-harmonic variation of angular frequency x is considered and all the quantities, unless otherwise indicated, are represented in terms of their complex phasors. Relevant time-domain features of the grounding system are obtained by means of a standard IFT.
The goal of the proposed hybrid method [18][19][20][21] is to obtain a frequency-dependent equivalent electric circuit, as done in other fields [22][23][24][25] for the grounding system where all the inductive, capacitive and conductive couplings among different conductor elements are accounted for and their analysis can be carried out in the frequency domain.
It should be highlighted that no low-frequency approximations are introduced, neither intrinsically approximate techniques are adopted, like, e.g., [26][27][28] based on complex images, to overcome the direct computation of Sommerfeld's integrals arising in the Green's functions. The choice of a direct approach is motivated by the need to analyzing grounding systems even during the early time of the lightning strike when the extremely large dimensions reduce dramatically the accuracy of the currently available approximate approaches.

Mathematical model
We assume that the conductors of the grounding network are completely buried in the earth that, without loss of generality, is considered stratified in N S planar layers. In the frequency domain, each earth layer is characterized by a complex conductivity r i ¼ r i þ jxe 0 e r;i (with i ¼ 1; . . .; N S ), where r i is the electric conductivity, e r;i is the relative dielectric permittivity, e 0 is the vacuum dielectric permittivity, and x is the angular frequency. For all the layers the vacuum permeability, l 0 , is assumed. Frequency dependence of ground permittivity is neglected, because of its little influence on the relevant performance in the considered frequency range of interest.
The grounding system is partitioned into N B segments that, in the framework of the networks theory, can be studied as elemental oriented branches connecting a pair of nodes (the terminals of the branches). The discrete grounding system is assumed to have N N nodes. Each kth conductor carries a uniform longitudinal current I LO;k flowing between the two terminal nodes whose scalar electric potentials (SEP) V i are defined with respect to the infinitely remote earth chosen as the zero potential reference. The introduction of the nodal voltages is necessary because at high frequencies the voltage drops along the branches are not negligible due to the electromagnetic coupling and to the increase of the internal impedance. Additionally, each branch drains a radial leakage current I LE;k to the surrounding earth. Such a radial current is assumed to be uniform, provided that the discretization is sufficiently fine. The grounding system is energized by the injection of single-frequency sinusoidal source currents F m at M nodes (Fig. 2).
Due to these assumptions, the SEP along each branch is not constant and varies in someway between the two voltage values at the terminal nodes. However, since the segment length is much shorter than the minimum wavelength, we can assume that the potential U k over the kbranch is constant and equal to the average of the voltages V i and V j at the two terminal nodes [29] Consequently, by applying (1) at each branch, the relationship between the N B Â 1 column vector U of the average SEP of segments and the N N Â 1 column vector V n of the SEP of nodes can be expressed as where K is a matrix of dimensions N B Â N N whose elements are Looking at the leakage currents and the average SEP of all branches, we found that there is a matrix relationship ensuing from the application of the Galerkin's MoM as [18] where I LE is the N B Â 1 column vector of the leakage currents and Y GC is a N B Â N B square admittance matrix accounting for all the conductive and capacitive couplings among the discretized segments, as shown in the next subsection. Subdividing each leakage current I LE;k into two equal contributions I LE;k =2 that are concentrated in the terminal nodes [29] (the approximation is reasonable under the assumption that the segmentation is fine with respect to the wavelength), it is possible to obtain the following relationship between the concentrated currents J LE at the nodes and the leakage currents I LE assumed uniformly distributed over each short branch where J LE is a N N Â 1 column vector. Under the above assumptions, the equivalent circuit can be studied making use of the conventional nodal analysis, resulting in the following equation [12,13] where F is the vector of the known external current sources, Y RL is the N B Â N B square admittance matrix accounting for the internal impedance of the branches and the magnetic coupling among them, obtained through the Galerkin's MoM as described in the next paragraph, and A is the N N Â N B branch-to-node incidence matrix whose elements are defined as Substituting Eqs. (2), (4), (5) into (6), the final equation is readily obtained as Once the nodal voltages are obtained, all the other relevant quantities, i.e., branch currents, leakage currents and potential on the ground, can be obtained straightforwardly.

Computation of impedances
Each branch is characterized by a lumped self-impedance.
Inductive and capacitive couplings with other branches are accounted for through mutual impedances. The admittance (a) (b) Fig. 2 Typical wind turbine grounding system arrangement (a) and section of the discretized model (b). In a I F is the fault current flowing into the earthing system and J F is the fault current density injected into the ground. The ground is assumed stratified and each layer is characterized by the resistivity q i and the permittivity e i . In b J LE;j is the leakage current at node j, I LE;k is the leakage current dispersed by the element k, I BR;k is the current flowing through the element k, and V i is the potential at the node i matrix Y RL can be computed as the inverse of the impedance matrix Z RL which account for the resistive-inductive couplings, whose elements are obtained as [12,13] In (9), Z ij is the internal impedance of the cylindrical conductor with radius r 0 at frequency x that is given by where k r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Àjxl c r c p , r c and l c are, respectively, the conductivity and the absolute permeability of the conductors. J 0 and J 1 are the Bessel function of the first kind and order 0 and 1, respectively. The mutual induction between conductors i and j is computed as where G A r; r 0 ð Þ is the magnetic vector potential dyadic Green function in the layered earth, l i and l j are the lengths of the ith observation and jth source segments with tangent unit vectorst andt 0 , respectively.
As concerns the admittance matrix Y GC , it can be computed as the inverse of the impedance matrix Z GC which account for the conductive-capacitive couplings among elements [12,13]. By expanding the leakage currents in uniform zero-order basis functions N i defined over the ith branch as N i ¼ 1=l i and by applying a standard Galerkin's MoM [30][31][32], the element Z CG;ij is where G u r; r 0 ð Þ is the scalar potential Green's function in the stratified earth.

Green's functions
The dyadic magnetic Green's function in the stratified earth can be computed as follows [33]: with The scalar potential Green's function is The spectral domain expressions for the elements of the Green's functions can be recast by using the transmissionline (TL) analogy described in [34,35]. More specifically, according to the TL notation (Fig. 3), we obtain where l 0 and e 0 are the free space permeability and permittivity, respectively, k 0 is the wavenumber of the free In (16) Fig. 3: they consist of a cascade connection of uniform trasmission line sections where each section represents a layer whose boundaries are placed at z n and z nþ1 and is characterized by a propagation wavenumber k p zn and characteristic impedance Z p n , with p ¼ e or h. The computation of currents and voltages can be performed efficiently through the recursive computation of the generalized voltage reflection coefficients C p n looking to the left (towards the depth of the ground) at the n-interface between medium n À 1 and medium n and C ! p n (with p ¼ e or h) looking to the right (towards the air) at the n-interface between medium n and medium n þ 1. Numerical details are reported in [34].
The computation of the Sommerfeld integrals is performed with the double exponential-type quadrature formulas recently proposed in [36]. Since the ground is lossy, the poles of the Green's functions, corresponding to guided propagation modes in the structure, move from the real axis into the forth quadrant; in addition, the Green's functions in intermediate layers, where the earthing systems are placed in practice, do not present branch points; consequently the integration can be computed straightforwardly on the real axis without deforming the Sommerfeld integral path into the first quadrant as it is usually done in order to avoid poles or branch points.

Singular term
The evaluation of the admittance and impedance matrices requires the computation of each linear wire segment comprising the earthing system. Since the basis functions for currents and charges are constant and the vector bases for current have constant direction on each linear segment, when the observation j and source i segments coincide, both the integrals in (11) and (12) reduce to a scalar integral of the form where c a is an appropriate coefficient and K r; r 0 ð Þ is the corresponding kernel (e.g., the vector or scalar Green's function in the space domain). As it is well known, the Green's function in the spectral domain inside the n ground layer is always composed by a first term that represents the direct ray between the source and the observation point and a second term that represents the superposition of the rays that undergo partial reflections at the upper and lower slab boundaries before reaching the field point. The first term, when is transformed back from the spectral to the spatial domain, contains a singular term that is of the form 1 / R and that must be accurately treated when the MoM analysis is applied to wire structures in order to have nonessential (i.e., integrable) singularities.
Following the procedure presented in [37,38], since the sources and the potentials about a linear tubular section present a rotational symmetry, assuming, without loss of generality, sources distributed uniformly on a cylindrical tube of constant radius r 0 centered along the z-axis and an observation point in cylindrical coordinates r 0 ; 0; z ð Þ, the wire kernel K r; r 0 ð Þ is defined as where R ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi z À z 0 ð Þ 2 þ2r 2 0 À 2qr 0 cos / 0 q . With simple manipulations that are not here reported for the sake of conciseness, introducing the substitution n ¼ p À / 0 ð Þ=2, the distance R can be rewritten as , and the kernel integral (18) can be expressed as where F b ð Þ is the complete elliptic integral of the first kind: The first term in (19) contains a logarithmic singularity as b ! 1, i.e. z ! z 0 , while the second term is smooth and can be easily evaluated with standard Gaussian quadrature formulas. When the first term in (19) is included in the source linear integral in (17), the logarithmic singularity must be accurately integrated on the linear element, and this job can be carried out by using the modified quadrature rules proposed by Wandzura in [39].

Lightning current parameters
Lightning depends on a number of factors like type (i.e., positive or negative discharges), first or subsequent strokes, direction, channel shape and so forth. Each feature is characterized by its uncertainty; the phenomenon is typically analyzed on a statistical basis, corroborated by means of extensive measurements collected in instrumented towers [40][41][42][43] or in artificially-initiated lightning channels [44,45]. Several expressions have been considered in the past for the LEMP current waveform: all of them satisfy lightning current waveform parameters which are generally agreed upon [46,47]. In particular, the empirical equation that is most widely used to reproduce the channel-base currents is the Heidler's function [48,49], which reads where I 01 is the peak current, s 11 and s 21 are time constants determining current rise-time and decay-time, respectively, n is the current steepness factor, and g 1 is the correction factor of the peak current, given by Distinguishing between first and subsequent strokes, in this paper we used a single Heidler's function for the first stroke and a sum of two functions like (21) for the subsequent strokes; parameters are given in Table 1.
Other expressions have been also adopted [51], in order to use a unique formula lasting for the whole time period of interest, avoiding the possible time-derivatives discontinuities appearing in the switch between first and subsequent stroke expressions.

Testing and validation of the method
The accuracy and efficiency of the proposed hybrid method is tested at first comparing its results with the predictions of some analytical formulas for the calculation of resistances to ground valid at DC proposed by the IEEE Standard 142 [52]. Six relevant configurations and their analytical formulas are reported in Table 2. The following assumptions  have been considered: ground resistivity q ¼ 100 X Â m, length of arm L ¼ 10 m, depth s ¼ 0:8 m, wire radius a ¼ 3:5 mm. The comparison among the results of the hybrid method and the predictions of the analytical formulas is reported in Table 3, considering a unit injected current. Moreover, the 3D maps of the potential on the ground are reported in Fig. 4. Excellent agrement between the numerical results and the predictions through analytical formulas is shown in each of the examined case.
Successively, the proposed approach has been applied to classical transient test cases to demonstrate the effectiveness of the procedure and to simultaneously asses the influence of the dominant parameters on the dynamic behavior of grounding electrodes under lightning strokes.
The considered test cases are two and are taken from [53]; they are based on the following configurations: -horizontal copper-made electrode with length 30 m (denoted as long electrode) and diameter 1.4 cm, buried at a depth of 0.8 m in a resistive soil with e r;g ¼ 10 and q g ¼ 300 Xm, whose results are reported in Fig. 5; -horizontal copper-made electrode with length 3 m (denoted as short electrode) and diameter 1.4 cm, buried  at a depth of 0.8 m in a conductive soil with e r;g ¼ 10 and q g ¼ 300 Xm, whose results are reported in Fig. 6.
In each figure the transient potential v t ð Þ is reported, computed by means of the proposed approach under the current pulse i t ð Þ of a first or subsequent return stroke whose parameters are reported in Table 1. The transient voltage is subdivided into two terms as proposed in [53] v where R is the DC grounding resistance. The first term in (23), named resistive component of the voltage drop, approximates the linear resistive voltage drop in the grounding system, while the second term, named reactive component of the voltage drop, accounts for the frequencydependent phenomena, i.e., the inductive and capacitive couplings between the elements of the grounding system.
The comparison between the rise voltage predictions of the two proposed methods for all the considered cases shows that the hybrid approach tends to estimate lower picks and higher rise times than the Grcev [53] results, based on a electromagnetic full-wave approach: the hybrid approach tends to estimate higher capacitive couplings between the elements of the grounding system, flattening the surge front and reducing the voltage picks.

Results
Even if most of the wind turbines manufacturers specify the grounding system configurations for their own towers, a typical wind turbine earthing system arrangement is shown in Fig. 7 [5]. It mainly consists of a square of galvanized steel flanges (gray line) buried at 2 m depth, two copper ring wires (black line) at two different levels, the smaller one at 5 cm depth while the larger one at 55 cm depth, four wires connecting the inner circle with the outer one, and additional eight copper wires connecting the outer circle and the square with the tower. Additional horizontal and vertical electrodes, respectively of length L h and L v , may be added to the ground system in order to obtain prescribed values of the ground resistance.
The grounding system is placed in a homogenous soil of low conductivity, q g ¼ 1200 Xm, which is a realistic scenario for typical wind farms. The relative permittivity is assumed to be e r ¼ 9. Figure 8 show the transient potential of the basic grounding system without any additional vertical or horizontal electrode, under a first or subsequent stroke. The inset of Fig. 8a shows the potential distribution on the ground at DC when the injected current is equal to 1 A. It can be seen in Fig. 8a that the ground system shows a resistive behavior under the first stroke, but when a subsequent stroke is considered, which is characterized by a fast rising front with higher frequency content (up to 10 MHz), it can be observed in Fig. 8b that the reactive component of the voltage drop plays a paramount role in the initial transient response, giving rise to a voltage peak that is higher than the peak predicted by the resistive response only.
Successively, we investigated the influence of additional vertical electrodes on the transient behavior of the grounding system. The basic earthing system has been upgraded with four vertical electrodes of variable length L v as shown in Fig. 7. Results reported in Fig. 9 show that additional short electrodes (3 or 5 m long) do not significantly decrease the transient response of the grounding system, under either a first or a subsequent stroke, with respect to the basic case. In order to clearly decrease the transient voltage it is necessary to resort to relatively long electrodes, e.g., with lengths ranging between 10 and 15 m.
Then we assessed the effect of horizontal electrodes of various length L h ; the grounding system has been upgraded with four horizontal electrodes connected at the four vertices and buried at 2 m depth as shown in Fig. 7. Figure 10 show the transient voltage raise of the grounding system at (a) (b) Fig. 7 Typical wind turbine grounding system arrangement: a prospective view; b front view (a) (b) Fig. 8 Computed results for the wind grounding system without additional horizontal or vertical electrodes. The inset in a shows the DC potential map on the ground when the injected current is 1 A. The DC resistance is R ¼ 40:04 X. a First stroke, b subsequent stroke the feeding point. Compared to the case of vertical electrodes, these results clearly show that there is no appreciable difference between the vertical and horizontal case, leading to the conclusion that the horizontal case is preferable because placing horizontal electrodes is usually cheaper especially in rocky terrains. It should also be observed that the computed DC resistances for the case of horizontal electrodes are slightly lower than the corresponding resistances in the case of vertical electrodes; this observation corroborates that horizontal electrodes are preferable.
To conclude the discussion of the results, a last remark is necessary about the protection level ensured by the earthing system according to international standards such as the IEC 62305 [54]. The earthing system usually comprises a single integrated structure that is suitable for several purposes, i.e., lightning protection, power and telecommunication systems earthing. The standard recommends that in dealing with the dispersion of the lightning current into the ground, the shape and dimensions of the earthing system are the important criteria to minimize any potentially dangerous overvoltages. In general, a low earthing resistance at industrial frequency is recommended, usually below 10 X. The standard, in the part number three [54,Section 5.4], suggests operative formulas for two basic types of earth electrode arrangements: type-A that comprises horizontal or vertical earth electrodes installed outside the structure and type-B that comprises either a ring conductor external to the structure to be protected, in contact with the soil for at least 80 % of its total length, or a foundation earth electrode forming a closed loop. The Type-B clearly apply to the present case under analysis. In this case the standard prescribes that for the ring earth electrode (or foundation earth electrode), the mean radius r e of the area enclosed by the ring earth/fundation electrode shall be not less than a prescribed value l 1 which depends on the lightning protection system (LPS) class and on the resistivity of the soil q g . In our case, with q g ¼ 1200 Xm, the prescribed values l 1 are 27 m for Class I, 14 m for Class II, 5 m for Class III and IV. In our case, considering only the outer ring of diameter 13.6 m (r e ¼ 6.8 m, the Class should be III. Anyway, the standard adds that when the required value of l 1 is larger than the convenient value of r e , additional horizontal or vertical electrodes shall be added with individual lengths l h ¼ l 1 À r e and l v ¼ l 1 À r e ð Þ=2. The standard recommends that the number of electrodes shall be not less than the number of the down-conductors, with a minimum of two. Hence, adding horizontal electrodes of length 15 m, we could obtain the protection Class II while adding vertical electrodes of the same length, we could reach Class I according to the standard. It is interesting to note that the overperformance claimed by the standard of the vertical electrodes on the horizontal ones is not observed on the computed results. In order to better investigate this issue, we have repeated the simulations changing the value of the grounding resistivity q g between 30 and 3000 Xm, and we have always observed that the performance of vertical and horizontal electrodes of the same length is almost the same.

Conclusions
The suitability of a new hybrid approach for the computation of the behavior of ground systems in the frequency domain has been investigated. The formulation is based on circuit theory and method of moments, to fully account for resistive, inductive and capacitive couplings, and the analysis is carried out in the frequency domain. When needed, the corresponding transient response is obtained by means of the inverse Fourier transform. The method has been used to characterize the transient behavior of wind turbine grounding systems in several configurations, investigating the relevant aspects described in International Standards.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.