Kaddour Arzag*
| Zin-Eddine Azouz
© 2026 The authors. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
This paper introduces an investigation into the influence of a hill presence on the electromagnetic (EM) field generated by a lightning strike to a wind tower (WT) blade. The simulation method involves the three-dimensional finite-differences time-domain (3D-FDTD) considering a perfectly conducting ground. The first part compares the results computed in the presence and in the absence of the hill. After that, we examine the influence of the distance between the observation points and the hill on the EM field components. We then present a comparison between three cases of the WT location on the hill. At the end of this study, the obtained results show that the presence of the hill causes a reduction in both the vertical electric and azimuthal magnetic field magnitudes, and it increases the time of arrival (ToA) of these waves. However, the rise time is slightly affected (difference of 0.1 µs). Furthermore, the wave shape, rise time, and the ToA are affected by the change in distance. The change in the WT locations affects the initial peak of the EM field components. The 3D-FDTD approach yields efficient results of lightning EM fields. In future studies, we plan to consider more complex geometries that are closer to reality, particularly by taking into account soil conductivity.
3D finite difference time domain method, wind turbine tower, hill, lightning, electromagnetic fields, transmission line model
Wind turbine towers (WTT) are characterized by their high altitude and pointed blades (appearance of high electric fields at the pointed electrode); these two conditions make wind towers (WT) exposed to lightning strikes. Thus, the interaction between lightning and WTs results in the propagation of a current along the WT and along the lightning channel. Furthermore, the current propagation generates an electromagnetic (EM) field, which can be coupled to electric, electronic, and telecommunication systems. The coupling phenomenon can affect the proper functioning or cause the complete shutdown of these systems.
Thus, to quantify the lightning EM disturbances on the different systems, it is essential to adopt numerical methods such as the finite-difference time-domain (FDTD) method [1], the Method of Moments (MoM) [2], and the Finite-Element Method (FEM) [3]. In this context, several papers have been published, aiming to study lightning strikes on WTs and their effects on systems installed on these WTs or on those located at close or far distances from them. Among these works, Wang and Zhang [4] developed a numerical method for determining the transient electric field inside WTs struck by lightning. This method is based on the use of a coupled Resistor, Inductor, and Capacitor (RLC) network to model the tower and the computation of the network lightning transient responses by solving its equivalent circuit. The validity of the proposed method has been confirmed by comparing the obtained results with simulation results. Nagao et al. [5] used the FDTD method to calculate the transient current distributions within the grounding systems of a wind-turbine-generator (WTG) tower struck by lightning and its neighboring WTG tower. They then calculated, using the same method, the transient electric field inside these towers. Relevant conclusions were reached in this study, particularly regarding the reduction of the electric field. Yamamoto et al. [6] carried out analytical and experimental investigations of magnetic-field transitions in a WT nacelle. The authors used the FDTD method for the analytical studies and a reduced-scale WT generator model for the experimental investigation. Wang and Zhang [7] presented an efficient algorithm for computing the effects of magnetic transients in a WT tower struck by lightning. The obtained results made it possible to identify potential hazards and to design appropriate lightning protection systems. Ahmed and Ishii [8] analyzed the lightning surge response of interconnected WT grounding systems using both the Numerical Electromagnetic Code (NEC-4), based on the MoM, and a Virtual Surge Test Lab (VSTL) employing the FDTD method. The results obtained using the NEC code and the VSTL were then compared to assess the usefulness of these two methods for the transient analysis of large grounding systems. Rizk et al. [9] used the FDTD method to compute the earth potential rise (EPR) and the transient grounding resistance (TGR) in order to investigate the grounding system of the lightning protection system associated with wind power generation systems. Several grounding system configurations were investigated using the FDTD method, which is well suited for modeling the three-dimensional geometric complexity of the structure under study.
Rizk et al. [10] studied, using the 3D-FDTD method, the impact of ground conductivity, relative permittivity, the position of the struck WTT, and the lightning current rise time on the lightning EM fields impinging on the underground cable sheath. The results clearly demonstrated the influence of soil conductivity, ground permittivity, and the position of the struck WTT on the lightning current flowing through the cable sheath and its associated electric field. The study [11] used the same method to calculate, for different ground configurations, the lightning-induced voltages (LIVs) appearing on a three-phase overhead line due to strokes to a nearby grounded WTT. The influence of ground permittivity on the peak values of the LIVs was also examined in this study. Araneo and Celozzi [12] proposed a hybrid approach to analyze the transient behavior of WT earthing systems under lightning strikes. The hybrid approach combines circuit theory with the MoM in the frequency domain. The results obtained provided a better understanding of the nature of the early-time transient response of the grounding systems under consideration and enabled the development of useful design guidelines.
Pham [13] addressed the effects of lightning-induced overvoltages on WTs and mixed overhead-cable lines using the Electromagnetic Transients Program (EMTP). Sun et al. [14] presented a surge analysis of lightning strikes on overhead lines in a wind farm. In this study, mathematical models of the different system components were implemented in the EMTP environment to obtain relevant results and insights. Nam et al. [15] analyzed, using the Electromagnetic Transient Program/Alternative Transients Program (EMTP/ATP), several factors affecting lightning overvoltages when lightning strikes a WT blade in a wind farm. Based on this analysis, the authors established guidelines for selecting an effective WT protection system. Upward lightning (UL) at WTs was investigated by Stucke et al. [16] using direct UL measurements conducted at instrumented towers. Their study aimed to evaluate the effectiveness of upward lightning detection by lightning location systems (LLSs). Cai et al. [17] proposed introducing an impedance discontinuity within the turbine into the electric-field calculation model for WT return strokes. The results demonstrated that models assuming constant reflection coefficients are inaccurate, whereas adopting variable reflection coefficient values yields a more accurate agreement. Da Silva et al. [18] examined the performance of a wind farm turbine grounding system (WTGS) interconnected by an overhead ground wire (OHGW) and compared it with a traditional WTGS interconnected by underground bare cables for soils with low-frequency resistivities. The full-wave EM software XGSLab was used to model the wind farm, and the transient responses were calculated using the Partial Element Equivalent Circuit (PEEC) method. The study showed that the transient responses increase with increasing soil resistivity. Moreover, the WTGS interconnected by an OHGW generates higher transient currents and ground potential rise (GPR) at adjacent transition poles, whose function is to directly connect the power cables to the medium-voltage distribution line, than those obtained with the traditional WTGS interconnected by underground bare cables.
Indeed, all previous studies [4-18] were carried out assuming flat ground conditions. However, in reality, wind farms are often located in regions characterized by complex topography where WTs are installed on hills. The presence of such geographical reliefs can significantly affect the propagation of EM disturbances radiated by lightning strikes to WTs. This effect is related to the elevation, shape, and relative position of the WTs on the hills. Moreover, all studies employing the FDTD method for modeling purposes are based on the Yee algorithm [1]. When applied to three-dimensional problems, this algorithm requires the implementation of complex boundary conditions, such as perfectly matched layers (PMLs) [19].
Against this background, the novelty of the present work lies in studying, analyzing, and evaluating for the first time the influence of hills on the EM field components generated by lightning strikes to WT. To the best of the authors’ knowledge, this is the first study in which a three-dimensional FDTD method with uniaxial perfectly matched layer (UPML) boundary conditions, based on the Taflove formulation [20-22], is applied to model this problem geometry. The Taflove formulation enables the simultaneous computation of EM field components and the implementation of UPML boundary conditions, which results in a significant reduction in computational time.
In this investigation, three different WT positions and two observation points are considered. To achieve this objective, the transmission line (TL) model [23], extended to account for a tall strike object [24], is employed to describe the current distribution along the lightning channel and within the WT. A comparative analysis of the results obtained in the presence and absence of a hill is then performed to draw meaningful conclusions. The EM field components are computed at two different distances from the hill foot. Three parameters are considered in this comparative study: the magnitudes of the EM field components, the rise time, and the time of arrival (ToA). The numerical code developed for this purpose was validated in study [22].
2.1 Geometry
In order to achieve the main objective described previously in the Introduction section, the geometry illustrated in Figure 1 is adopted. As shown in Figure 1(a), a WTT is placed at the top of a hill. The WTT has a total height of 155 m, including three blades, each with a length of 60 m. As illustrated in Figures 1(b) and 1(c), the considered hill is characterized by a trapezoidal profile in the y-z plane. The top of the hill has a rectangular shape with a length of 800 m and a width of 96 m in the x-y plane. The base of the hill is also rectangular, with dimensions of 800 m × 440 m in the x-y plane. The height of the hill is 200 m, and its slope angle is 30°. Indeed, the ground depth beneath the base of the hill is 300 m. The WT is struck by lightning with a channel height of 850 m, measured between the tip of the WT blade and the cloud. The computational domain accommodating the previously described system has dimensions of 800 m × 1360 m × 1505 m along the x-, y-, and z-axes, respectively. The observation points, at which the EM field components are computed, are located at horizontal distances of 50 m and 500 m from the foot of the hill along the y-axis direction.
Figure 1. Geometry adopted for the computation of the lightning electromagnetic field
2.2 Modeling
The EM field components were calculated using the 3D-FDTD method based on the Taflove formulation. This method requires discretization of the entire computational domain into cubic or parallelepiped cells. In this study, a cell size of 5 m × 4.8 m × 5 m was adopted. The discretization step along the y-axis (4.8 m) was selected to ensure a hill slope angle of 30°. The time step was set to 3 ns in accordance with the Courant-Friedrichs-Lewy (CFL) stability condition [25], and the total observation time was 10 μs. UPML boundary conditions were applied at the six faces of the computational domain to eliminate unwanted reflections from the domain boundaries. Since the distance between the observation points and the WT does not exceed one kilometer, the soil was assumed to be perfectly conducting (i.e., having infinite conductivity). This assumption has been validated in several previous studies [26-28].
The adopted 3D-FDTD method based on the Taflove formulation is expressed by Eqs. (1(a-d)) as follows:
$\begin{gathered}D_x^{n+1}\left(i+\frac{1}{2}, j, k\right)=\left(\frac{2 \varepsilon K_y-\sigma_y \Delta t}{2 \varepsilon K_y+\sigma_y \Delta t}\right) \cdot D_x^n\left(i+\frac{1}{2}, j, k\right) \\ +\left(\frac{2 \varepsilon \Delta t}{2 \varepsilon K_y+\sigma_y \Delta t}\right) \cdot\left[\frac{H_z^{n+1 / 2}\left(i+\frac{1}{2}, j+\frac{1}{2}, k\right)-H_z^{n+1 / 2}\left(i+\frac{1}{2}, j-\frac{1}{2}, k\right)}{\Delta y}-\frac{H_y^{n+1 / 2}\left(i+\frac{1}{2}, j, k+\frac{1}{2}\right)-H_y^{n+1 / 2}\left(i+\frac{1}{2}, j, k-\frac{1}{2}\right)}{\Delta z}\right]\end{gathered}$ (1a)
$\begin{gathered}E_x^{n+1}\left(i+\frac{1}{2}, j, k\right)=\left(\frac{2 \varepsilon K_z-\sigma_z \Delta t}{2 \varepsilon K_z+\sigma_z \Delta t}\right) \cdot E_x^n\left(i+\frac{1}{2}, j, k\right) \\ +\left[\frac{1}{\left(2 \varepsilon K_z+\sigma_z \Delta t\right) \varepsilon}\right] \cdot\left[\left(2 \varepsilon K_x+\sigma_x \Delta t\right) \cdot D_x^{n+1}\left(i+\frac{1}{2}, j, k\right)-\left(2 \varepsilon K_x-\sigma_x \Delta t\right) \cdot D_x^n\left(i+\frac{1}{2}, j, k\right)\right]\end{gathered}$ (1b)
$\begin{gathered}B_x^{n+3 / 2}\left(i, j+\frac{1}{2}, k+\frac{1}{2}\right)=\left(\frac{2 \varepsilon K_y-\sigma_y \Delta t}{2 \varepsilon K_y+\sigma_y \Delta t}\right) \cdot B_x^{n+1 / 2}\left(i, j+\frac{1}{2}, k+\frac{1}{2}\right) \\ +\left(\frac{2 \varepsilon \Delta t}{2 \varepsilon K_y+\sigma_y \Delta t}\right) \cdot\left[\frac{E_z^{n+1}\left(i, j+1, k+\frac{1}{2}\right)-E_z^{n+1}\left(i, j, k+\frac{1}{2}\right)}{\Delta y}-\frac{E_y^{n+1}\left(i, j+\frac{1}{2}, k+1\right)-E_y^{n+1}\left(i, j+\frac{1}{2}, k\right)}{\Delta z}\right]\end{gathered}$ (1c)
$\begin{gathered}H_x^{n+3 / 2}\left(i, j+\frac{1}{2}, k+\frac{1}{2}\right)=\left(\frac{2 \varepsilon K_z-\sigma_z \Delta t}{2 \varepsilon K_z+\sigma_z \Delta t}\right) \cdot H_x^{n+1 / 2}\left(i, j+\frac{1}{2}, k+\frac{1}{2}\right) \\ +\left[\frac{1}{\left(2 \varepsilon K_z+\sigma_z \Delta t\right) \mu}\right] \cdot\left[\left(2 \varepsilon K_x+\sigma_x \Delta t\right) \cdot B_x^{n+3 / 2}\left(i, j+\frac{1}{2}, k+\frac{1}{2}\right)-\left(2 \varepsilon K_x-\sigma_x \Delta t\right) \cdot B_x^{n+1 / 2}\left(i, j+\frac{1}{2}, k+\frac{1}{2}\right)\right]\end{gathered}$ (1d)
Similar expressions can be derived for the remaining electric and magnetic field components along the y- and z-axes.
$\sigma_x(x)=\left(\frac{x}{d}\right)^m \cdot \sigma_{\max }$ (2)
$K_x(x)=1+\left(K_{x, \max }-1\right) \cdot\left(\frac{x}{d}\right)^m$ (3)
$\sigma_{\max }=-\frac{(m+1) \ln (R(0))}{2 \eta d}$ (4)
$R(\theta)=e^{-2 \eta \cos \theta \int_0^d \sigma(x) d x}$ (Reflection error) (5)
$\eta=\sqrt{\frac{\mu}{\varepsilon}}$ (6)
where,
d: The PML area thickness.
x: A positive integer corresponding to the layer number ($(0<x<d)$.
In Eqs. (1a) and (1b), $D_x^n$ and $E_x^n$ represent the electric density and the electric field intensity along the $x$- axis and at the step of time $n$, respectively. The parameters $B_x^{n+1 / 2}$ and $H_x^{n+1 / 2}$ represent the magnetic flux density and the magnetic field intensity, respectively, as given in Eqs. (1c) and (1d). These quantities are also defined along the x-axis but at the half time step, since each magnetic field component is located between two electric field components and vice versa. The indices i, j, and k denote the spatial position indices. By defining the multiplicative coefficients of the electric and magnetic field components and their corresponding flux densities, a unified treatment of these quantities can be achieved in both the interior computational domain and the PML region. Moreover, the parameters σ and K, associated with the computational domain, depend on the nature of the medium. For example, in free space, their values are σ = 0 and K = 1, whereas in the PML region these parameters are assumed to follow a polynomially graded profile, as given by Eqs. (2) and (3).
In this study, the source of EM disturbances is represented by the lightning current. The mathematical modeling of the latter can be carried out in two steps. The first step consists of specifying the injected current at the top of the WT, which corresponds to the attachment point between the lightning channel and the WT. The second step concerns the spatiotemporal distribution of the lightning current along the WT and the lightning channel.
The injected current at the top of the WT is modeled as the sum of two Heidler functions [29], representing a subsequent return stroke, as expressed in Eq. (7). All numerical parameters used in this simulation are listed in Table 1.
$\begin{gathered}i(0, t)=\frac{i_{01}}{\eta_1} \frac{\left(t / \tau_{11}\right)^{n_1}}{1+\left(t / \tau_{11}\right)^{n_1}} \exp \left(\frac{-t}{\tau_{12}}\right) +\frac{i_{02}}{\eta_2} \frac{\left(t / \tau_{21}\right)^{\mathrm{n}_2}}{1+\left(t / \tau_{21}\right)^{\mathrm{n}_2}} \exp \left(\frac{-\mathrm{t}}{\tau_{22}}\right)\end{gathered}$ (7)
where, i01 and i02 denote the current amplitudes, $\tau_{11}$ and $\tau_{12}$ are the front time constants, $\tau_{21}$ and $\tau_{22}$ are the decay time constants, and n1, and n2 are the exponents.
$\eta_1=\left[-\left(\frac{\tau_{11}}{\tau_{12}}\right)\left(n_1 \cdot \frac{\tau_{12}}{\tau_{11}}\right)^{\frac{1}{n_1}}\right], \eta_2=\left[-\left(\frac{\tau_{21}}{\tau_{22}}\right)\left(n_2 \cdot \frac{\tau_{22}}{\tau_{21}}\right)^{\frac{1}{n_2}}\right]$
The spatial and temporal distribution of the lightning current along the WT and the lightning channel is modeled using the TL model [23], extended to account for the presence of a tall object, as developed by Baba and Rakov [24]. The mathematical expressions of this model are given in Eqs. (8) and (9), where 
and 
denote the reflection coefficients of the current wave at the bottom and at the top of the WT, respectively, as defined in Eqs. (10) and (11).
Along the strike WT:
$\begin{gathered}0 \leq z^{\prime} \leq h \\ I\left(z^{\prime}, t\right)= \frac{1-\rho_{\text {top }}}{2} \sum_{n=0}^{\infty}\left[\begin{array}{c}\rho_{\text {bot }}^n \rho_{\text {top }}^n I_{s c}\left(h, t-\frac{h-z^{\prime}}{c}-\frac{2 n h}{c}\right) \\ +\rho_{\text {bot }}^{n+1} \rho_{\text {top }}^n I_{s c}\left(h, t-\frac{h+z^{\prime}}{c}-\frac{2 n h}{c}\right)\end{array}\right]\end{gathered}$ (8)
Along the lightning channel:
$z^{\prime} \geq h$$I\left(z^{\prime}, t\right)=\frac{1-\rho_{\text {top }}}{2} .\left[\begin{array}{c}I_{s c}\left(h, t-\frac{h-z^{\prime}}{c}\right) \\ +\sum_{n=1}^{\infty} \rho_{\text {bot}}^n \rho_{\text {top}}^{n-1}\left(1+\rho_{\text {top}}\right) I_{s c}\left(h, t-\frac{h+z^{\prime}}{v}-\frac{2 n h}{c}\right)\end{array}\right]$ (9)
$\rho_{\text {bot}}=\frac{Z_{\text {tow}}-Z_{g r}}{Z_{\text {tow}}+Z_{g r}}$ (10)
$\rho_{\text {top}}=\frac{Z_{\text {tow}}-Z_{\text {ch}}}{Z_{\text {tow}}+Z_{\text {ch}}}$ (11)
where, “n” is an index representing the number of reflections occurring between the top and the bottom of the tall object; “c” is the speed of light speed (3 × 108 m/s); “v” is the lightning current propagation speed along the lightning channel (v = 1.5 × 108 m/s); and Isc denotes the lightning short-circuit current. In addition, Zch is the characteristic impedance of the lightning channel, Zgr is the ground impedance, and Ztow is the tower impedance. The reflection coefficients (ρbot, ρtop) are set to ρbot = 1 and ρtop = −0.5; these values were also adopted by Baba and Rakov [24].
Table 1. Channel-base current parameters
|
i01 (kA) |
$\tau_{11}$(µs) |
$\tau_{12}$ (µs) |
i02 (kA) |
$\tau_{21}$ (µs) |
$\tau_{22}$ (µs) |
n1 |
n1 |
|
10.5 |
0.25 |
2.5 |
6.5 |
2.1 |
230 |
2 |
2 |
3.1 Lightning current
Figure 2 illustrates the injected current used as a source at the top of the WT blade, which represents a pointed location capable of attracting the lightning leader. The current waveform is computed using Eq. (7) and the parameters listed in Table 1. It can be observed that the resulting current exhibits a pulse-shaped waveform, with a peak magnitude of 12 kA reached at 0.7 µs, which corresponds to the current rise time.
The propagation of the current wave along the WT and the lightning channel has been simulated and analyzed. The spatiotemporal current distribution along the WT is computed using Eq. (8), whereas the current distribution along the lightning channel is obtained using Eq. (9).
Figure 2. Waveform of the uncontaminated current used in the computations
Figure 3 illustrates the space-time variations of the lightning current wave. According to these results, the maximum current value is observed at the bottom of the WT, as it corresponds to the superposition of the incident and reflected waves at ground level. The results also show that the current waveforms computed at the top of the WT and at heights of 250 m, 500 m, and 750 m above the turbine top exhibit identical shapes and amplitudes. This behavior is attributable to the use of the TL model, in which the current propagates along the lightning channel without attenuation or distortion. All current pulses exhibit multiple oscillations caused by multiple reflections of the current wave between the upper and lower points of the WT. However, these oscillations are damped as the simulation time progresses, since they are associated with transient phenomena. The rise time of all current pulse is 0.7 µs.
3.2 Electromagnetic field components
Three case studies related to the EM field computation are carried out. The first case investigates the influence of the presence of hills on the EM field components. The second case examines the effect of the WT location relative to the hill on the EM field waveforms and their magnitudes. The third case addresses the influence of the distance between the observation points and the hill on the EM field components.
3.2.1 Influence of hill presence on the electromagnetic field
Figure 4 presents the time variations of the vertical electric field “Ez”. The results compare two cases: the presence of the hill and the absence of the hill. This comparison is performed to assess the influence of the hill on the propagation of the EM field generated by lightning strikes to the WT. The EM field is evaluated at an observation point located 500 m from the bottom side of the hill, as shown in Figure 1. Table 2 summarizes the comparison between the two cases, considering three parameters: the peak magnitude, the rise time, and the ToA.
According to the results plotted in Figure 4 and the numerical values reported in Table 2, the initial peak of the vertical electric field Ez is significantly affected by the presence of the hill. Specifically, the initial peak decreases from 1.1 kV/m in the absence of the hill to 0.55 kV/m when the hill is present. In contrast, the rise time is only slightly influenced by the ground geometry in both cases (flat ground and ground with a hill). In the presence of the hill, Ez reaches its maximum value at 1.1 µs, whereas in the absence of the hill, the maximum occurs at 1.0 µs.
Table 2. Effect of the presence of a hill on the initial peak magnitude, rise time, and time of arrival (ToA) of Ez
|
Parameters |
Presence of Hill |
Absence of Hill |
|
Ez max (kV/m) |
0.55 |
1.1 |
|
Rise time (µs) |
1.1 |
1 |
|
ToA (µs) |
2.95 |
2.77 |
Furthermore, the summarized results in Table 2 indicate that the ToA is influenced by changes in the ground structure. Specifically, the ToA in the presence of the hill reaches 2.95 µs, which is longer than the value obtained in the absence of the hill (2.77 µs). In both computed electric field waveforms, multiple peaks are observed, resulting from multiple reflections of the lightning current wave between the upper point of the WT blade and the base of the WT. In addition, all computed electric field values are positive. Overall, the simulated vertical electric field is significantly affected by the presence of the hill in terms of peak magnitude and ToA, whereas the rise time is only slightly influenced by the change in ground geometry from flat terrain to hilly terrain.
The simulation results for the azimuthal magnetic field are presented in Figure 5. In this figure, the case with the hill present is represented by a solid blue line, whereas the case without the hill is depicted by a dashed red line. Table 3 summarizes the comparison of the parameters obtained for each case considered in this study.
The results presented in Figure 5 and Table 3 allow the following conclusions to be drawn:
•The peak magnitude of the azimuthal magnetic field is significantly affected by the ground structure. It reaches 1.65 A/m in the presence of the hill and 2.73 A/m in the absence of the hill.
•The rise time is only slightly influenced, being equal to 1.0 µs when the hill is present and 1.1 µs in the absence of the hill.
•A difference of 0.18 µs is observed in the ToA between the two cases. The ToA is 2.95 µs in the presence of the hill and 2.77 µs in its absence.
•Both azimuthal magnetic field waveforms exhibit multiple oscillations, caused by reflections of the lightning current wave between the top of the WT blade and the base of the WT.
Figure 5. Azimuthal magnetic field waveforms at a distance of 500 m from the hill
Note: The solid line represents the case with the hill present, whereas the dashed line corresponds to the case without the hill.
Table 3. Effect of the presence of a hill on the initial peak magnitude, rise time, and time of arrival (ToA) of Hx
|
Parameters |
Presence of Hill |
Absence of Hill |
|
Hx max (A/m) |
1.65 |
2.73 |
|
Rise time (µs) |
1.1 |
1 |
|
ToA (µs) |
2.95 |
2.77 |
The differences observed in the parameter values are mainly attributed to variations in the distance between the WT and the observation points. This distance is 720 m in the absence of the hill and 747 m when the hill is present. In addition, the differences are related to the elevation of the WT relative to the flat ground, which is 0 m in the absence of the hill and 200 m in the presence of the hill. These elevation differences modify the propagation paths of the EM field waves, thereby directly influencing the parameters considered in the comparative analysis.
3.2.2 Effect of the distance between the hill and the observation points
The vertical electric field and the azimuthal magnetic field are simulated at two horizontal distances from the base of the hill, namely 50 m and 500 m. Figure 6 shows the resulting vertical electric field waveforms, where the dashed red line corresponds to the results obtained at 50 m, while the solid blue line represents those computed at 500 m.
Figure 6. Vertical electric field Ez waveforms computed at two distances from the hill
Note: The solid line represents the case at 500 m, whereas the dashed line corresponds to the case at 50 m.
Table 4. Effect of the distance from the hill on the initial peak magnitude, rise time, and time of arrival (ToA) of Ez
|
Parameters |
50 m |
500 m |
|
Ez max (kV/m) |
0.73 |
0.55 |
|
Rise time (µs) |
1.48 |
1 |
|
ToA (µs) |
2 |
2.95 |
The numerical comparison of the parameters obtained at the two considered distances is presented in Table 4.
According to the results plotted in Figure 6 and the numerical values reported in Table 4, the vertical electric field waveform is significantly affected by the change in distance from 50 m to 500 m. The waveform computed at 50 m exhibits a pulse shape with a peak magnitude of 0.73 kV/m, reached at a rise time of 1.48 µs, and a ToA of 2.0 µs. This pulse is directly damped without the appearance of oscillations. In contrast, the electric field waveform, calculated at a distance of 500 m, clearly exhibits oscillations. In this case, the peak magnitude is 0.55 kV/m, occurring at a rise time of 1.0 µs, and the electric field reaches the observation point (at 500 m) at a ToA of 2.95 µs.
Figure 7. Azimuthal magnetic field waveforms computed at two distances from the hill
Note: The solid line represents the case at 500 m, whereas the dashed line corresponds to the case at 50 m.
Figure 7 presents the computed azimuthal magnetic field results at horizontal distances of 50 m and 500 m from the bottom of the hill. Table 5 provides a numerical comparison between the two cases considered in this part of the study.
For the case of the magnetic field computed at an observation point located 50 m from the hill base (represented by the dashed red line), the waveform exhibits a pulse shape with a peak magnitude of 1.78 A/m, a rise time of 1.44 µs, and a ToA of 2.0 µs. After reaching its peak, the pulse decays exponentially and exhibits small oscillations at approximately 4.5 µs and 6.5 µs. In contrast, for the second case corresponding to the magnetic field computed at an observation point located 500 m from the hill base, the waveform displays an oscillatory pulse shape. The first peak reaches a magnitude of 1.65 A/m, with a rise time of 1.0 µs, and the signal arrives at the observation point after 2.95 µs (ToA).
In conclusion, the EM field components, namely the vertical electric field and the azimuthal magnetic field, are significantly affected by the change in distance between the observation point and the bottom of the hill from 50 m to 500 m.
Table 5. Influence of the distance from the hill on the initial peak magnitude, rise time, and time of arrival (ToA) of Hx
|
Parameters |
50 m |
500 m |
|
Hx max (A/m) |
1.78 |
1.65 |
|
Rise time (µs) |
1.44 |
1 |
|
ToA (µs) |
2 |
2.95 |
The physical origin of the observed differences between the EM field components computed at the two considered distances from the hill (50 m and 500 m) is related to the relative positions of the observation points. The first observation point, located at a short distance of 50 m from the bottom of the hill, is weakly affected by multiple reflections of the lightning current wave between the top and the base of the WT. In contrast, the second observation point, situated at a distance of 500 m, is more exposed to interaction phenomena between the lightning channel and the WT. As a result, multiple oscillations appear in the corresponding field waveforms. These oscillations are caused by multiple reflections of the lightning current wave between the top of the WT and the bottom of the hill (Table 5).
3.2.3 Influence of the wind tower location on the hill
The investigation of the effect of the WT position on the EM field components constitutes the final part of this study. For this purpose, three cases are considered. In the first case, the WT is placed at the top of the hill. In the second case, the WT is located at the midpoint of the inclined surface of the hill.
Furthermore, although WTs are rarely installed at the base of hills, this configuration is considered here solely for comparative purposes. It should be noted that, in all three cases, the WT is located at a distance of 400 m from the front boundary of the computational domain. Moreover, two observation points are considered in this investigation: the first is located at a horizontal distance of 50 m from the bottom of the hill, and the second at a distance of 500 m from the hill bottom.
Figure 8. Vertical electric field Ez waveforms computed at 50 m from the hill for three wind turbine (WT) positions on the hill
Figure 8 shows the computed vertical electric field at the first observation point (50 m) for the three considered WT positions on the hill. The solid blue line represents the case where the WT is located at the top of the hill. The red dashed line corresponds to the WT positioned at the midpoint of the inclined surface, while the green dotted line represents the case where the WT is placed at the bottom of the hill.
Table 6 presents the numerical results of the vertical electric field component, including the maximum value of the first peak, the rise time, and the ToA.
Table 6. Influence of the wind turbine (WT) location on the initial peak magnitude, rise time, and time of arrival (ToA) of Ez
|
Parameters |
Top of Hill |
Middle of Hill |
Bottom of Hill |
|
Ez max (kV/m) |
0.73 |
1.32 |
−5.3 |
|
Rise time (µs) |
1.48 |
1.2 |
0.96 |
|
ToA (µs) |
2 |
1.05 |
0.6 |
The comparison between the plotted waveforms and the numerical results corresponding to the three investigated cases leads to the following observations:
•WT located at the top of the hill: The computed vertical electric field exhibits a pulse waveform with a maximum magnitude of 0.73 kV/m, a rise time of 1.48 µs, and a ToA of 2 µs. The pulse increases from zero to its maximum value and then decreases exponentially until it returns to zero.
•WT located at the middle of the inclined surface of the hill: In this case, the electric field waveform rises from zero to a maximum value of 1.32 kV/m, with a rise time of 1.2 µs and a ToA of 1.05 µs. After reaching its peak, the field decays exponentially toward zero.
•WT located at the bottom of the hill: The electric field waveform is characterized by a polarity reversal, changing from negative to positive values. This behavior is attributed to the reflection of the EM field wave at the hill wall, which is located at a very short distance from the observation point. The dominant peak occurs in the negative polarity, with a maximum magnitude of −5.3 kV/m and a rise time of 0.96 µs. Subsequently, the field increases to reach zero at 2.66 µs, then continues increasing to attain a positive value of 1.9 kV/m. Afterward, oscillations appear in the waveform. The ToA in this case is 0.6 µs.
Figure 9 shows the computed vertical electric field component at a horizontal distance of 500 m from the bottom of the hill. Table 7 presents a comparison of the numerical values of the vertical electric field for the different WT positions.
Figure 9. Vertical electric field Ez waveforms computed at a horizontal distance of 500m from the hill for the three considered wind turbine (WT) positions on the hill
It should be noted that, for all WT locations, the electric field remains positive, with no polarity reversal. The magnitude and rise time of the first peak obtained for the WT located at the top, middle, and bottom of the hill are, respectively, 0.55 kV/m at 1 µs, 0.6 kV/m at 1 µs, and 1.05 kV/m at 1.1 µs. Furthermore, the maximum electric field values observed for these three cases are 0.55 kV/m at 4 µs, 0.9 kV/m at 5.25 µs, and 1.16 kV/m at 6.8 µs, respectively. The corresponding ToA are 2.95 µs for the WT located at the hilltop, 2.1 µs for the WT positioned at the middle of the inclined surface, and 1.7 µs for the WT placed at the bottom of the hill.
Table 7. Effect of the wind turbine (WT) location on the initial peak magnitude, rise time, and time of arrival (ToA) of Ez
|
Parameters |
Top of Hill |
Middle of Hill |
Bottom of Hill |
|
Ez max (kV/m) |
0.55 |
0.6 |
1.05 |
|
Rise time (µs) |
1 |
1 |
1.1 |
|
ToA (µs) |
2.95 |
2.1 |
1.7 |
These results indicate that the maximum amplitude of the first peak decreases when moving from the case where the WT is located at the bottom of the hill to the case where it is placed at the top, with intermediate values for the WT positioned on the inclined surface. Conversely, the ToA increases from the bottom to the top position of the hill. The rise time is only slightly affected by the change in the WT position relative to the hill. The observed differences in the electric field parameters are mainly attributed to the variation in the direct distance between the observation point and the different WT locations on the hill.
The azimuthal magnetic field is also investigated for the three WT positions on the hill and for two horizontal distances from the bottom of the hill (50 m and 500 m). Figure 10 shows the temporal variations of the azimuthal magnetic field component Hx, while Table 8 summarizes the numerical values of the maximum magnitude of the initial peak, the rise time, and the ToA of Hx obtained for the three cases.
Figure 10. Azimuthal magnetic field Hx waveforms computed at 50 m from the hill for three wind turbine (WT) positions on the hill
Table 8. Influence of the wind turbine (WT) location on the maximum magnitude of the initial peak, the rise time, and the time of arrival (ToA) of Hx
|
Parameters |
Top of Hill |
Middle of Hill |
Bottom of Hill |
|
Hx max (A/m) |
1.78 |
3.27 |
32.9 |
|
Rise time (µs) |
1.44 |
1.67 |
0.97 |
|
ToA (µs) |
1.88 |
1 |
0.6 |
Thus, at the observation point located 50 m from the bottom of the hill, the first peak reaches a maximum magnitude of 1.78 A/m, with a rise time of 1.44 µs and a ToA of 1.88 µs, for the case where the WT is installed at the top of the hill. For the WT positioned at the middle of the inclined surface of the hill, the first peak is characterized by a maximum amplitude of 3.27 A/m, a rise time of 1.67 µs, and a ToA of 1 µs. When the WT is installed at the bottom of the hill, the first peak reaches 32.9 A/m, with a rise time of 0.97 µs and a ToA of 0.6 µs.
It should be noted that the highest amplitude of the azimuthal magnetic field is obtained for the WT located at the bottom of the hill, as this configuration corresponds to the shortest distance between the WT and the observation point. Moreover, the maximum amplitude of the first peak, the rise time, and the ToA are significantly affected by the WT position on the hill.
Figure 11 shows the simulation results of the azimuthal magnetic field Hx calculated at a horizontal distance of 500 m from the bottom of the hill. The corresponding numerical results are reported in Table 9. According to these results, it is clearly observed that the waveforms, the maximum magnitude of the initial peak, and the ToA are significantly affected by the change in the WT position on the hill. However, the rise time is not influenced by the WT location.
Figure 11. Azimuthal magnetic field Hx waveforms computed at 500 m from the hill for three wind turbine (WT) positions on the hill
Table 9. Influence of the WT location on the maximum magnitude of the initial peak, the rise time, and the time of arrival (ToA) of Hx
|
Parameters |
Top of Hill |
Middle of Hill |
Bottom of Hill |
|
Hx max (A/m) |
1.65 |
1.73 |
2.75 |
|
Rise time (µs) |
1 |
1 |
1 |
|
ToA (µs) |
2.95 |
2.1 |
1.7 |
For the first case, where the WT is located at the top of the hill, the maximum magnitude of the initial peak is 1.65 A/m, the rise time is 1 µs, and the ToA is 2.95 µs. When the WT is placed at the middle of the inclined surface of the hill, the initial peak reaches a maximum value of 1.73 A/m, with a rise time of 1 µs and a ToA of 2.1 µs. In the last case, corresponding to the WT installed at the bottom of the hill, the initial peak exhibits a maximum magnitude of 2.75 A/m, the rise time remains equal to 1 µs, and the ToA decreases to 1.7 µs.
Therefore, for an observation distance of 500 m from the bottom of the hill, the highest initial peak magnitude of Hx is obtained when the WT is located at the bottom of the hill, as this configuration corresponds to the shortest distance between the WT and the observation point.
In this paper, the 3D-FDTD method, combined with the TL model extended to include a tall object, has been implemented to investigate the EM fields generated by lightning strikes to a WT blade in the presence of a hill. The main objective was to extract relevant information regarding the key parameters of the resulting EM fields. To this end, several comparative analyses were conducted.
The obtained results demonstrate the effectiveness of the proposed approach based on the Taflove formulation of the 3D-FDTD method. This numerical method enables the investigation of the influence of geographical reliefs on the propagation of lightning-induced EM fields. The quantification of the EM field components facilitates the evaluation of their coupling with various systems, thereby contributing to the design of appropriate protection measures against EM disturbances.
Future work will focus on modeling and simulating more complex scenarios, including grounds with finite conductivity and the presence of water bodies such as rivers and lakes. In addition, the coupling between lightning-induced EM fields and overhead or buried systems connected to WTs will be considered as a perspective for further investigations.
|
Ex, Ey, Ez |
Electric field components |
|
Hx, Hy, Hz |
Magnetic field components |
|
i(0, t) |
Lightning current at the channel base |
|
i01, i02 |
Current amplitudes. |
|
Kx, Ky, Kz |
Real number |
|
n1, n2 |
Exponents |
|
R(θ) |
Reflection error |
|
t |
Time |
|
v |
Lightning return stroke velocity |
|
∆x, ∆y, ∆z |
Space increment |
|
Isc |
Lightning short-circuit current |
|
Zch |
Lightning channel characteristic impedance |
|
Zgr |
Ground impedance |
|
Ztow |
Tower impedance |
|
Greek symbols |
|
|
ε |
Absolute permittivity |
|
εr |
Relative permittivity |
|
η |
Characteristic wave impedance |
|
η1, η2 |
Numerical coefficients |
|
θ |
Incidence angle |
|
μ |
Absolute permeability |
|
σ |
Electric conductivity |
|
$\tau_{11}$, $\tau_{12}$ |
Front time constants |
|
$\tau_{21}$, $\tau_{22}$ |
Decay time constants |
|
$\rho_{bot}$, $\rho_{top}$ |
Reflection coefficients of the current wave at the bottom and at the top of the WT |
|
Subscripts |
|
|
i, j, k |
Position indices |
|
x, y, z |
Cartesian coordinates |
|
n |
Time step indices |
[1] Yee, K.S. (1966). Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Transactions on Antennas and Propagation, 14(3): 302-307. https://doi.org/10.1109/TAP.1966.1138693
[2] Harrington, R.F. (1968). Field Computation by Moment Methods. Macmillan.
[3] Sadiku, M.N.O. (1989). A simple introduction to finite element analysis of electromagnetic problems. IEEE Transactions on Education, 32(2): 85-93. https://doi.org/10.1109/13.28037
[4] Wang, X., Zhang, X. (2008). Electric field distribution inside wind turbine towers struck by lightning. In 2008 8th International Symposium on Antennas, Propagation and EM Theory, Kunming, China, pp. 496-499. https://doi.org/10.1109/ISAPE.2008.4735258
[5] Nagao, M., Nagaoka, N., Baba, Y., Ametani, A. (2009). FDTD electromagnetic analysis of a wind turbine generator tower struck by lightning. IEEJ Transactions on Power and Energy, 129(10): 1181-1186. https://doi.org/10.1541/ieejpes.129.1181
[6] Yamamoto, K., Chikara, T., Ametani, A. (2009). Numerical electromagnetic field analysis of transient magnetic fields in a nacelle caused by a lightning stroke to a wind turbine generator system. In International Conference on Power Systems Transients, Kyoto, Japan.
[7] Wang, X., Zhang, X. (2010). Calculation of electromagnetic induction inside a wind turbine tower struck by lightning. Wind Energy, 13(7): 615-625. https://doi.org/10.1002/we.382
[8] Ahmed, M.R., Ishii, M. (2012). Electromagnetic analysis of interconnected groundings at wind farm. IEEJ Transaction on Power and Energy, 132(7): 684-689. https://doi.org/10.1541/ieejpes.132.684
[9] Rizk, M.E.M., Badran, E.A., Abdel-Rahman, M.H. (2014). FDTD calculation of transient electromagnetic fields in the grounding system of wind towers due to direct lightning strikes. Journal of Lightning Research, 6: 1-13. https://doi.org/10.2174/1652803401306010001
[10] Rizk, M.E., Mahmood, F., Lehtonen, M., Badran, E.A., Abdel-Rahman, M.H. (2015). Investigation of lightning electromagnetic fields on underground cables in wind farms. IEEE Transactions on Electromagnetic Compatibility, 58(1): 143-152. https://doi.org/10.1109/TEMC.2015.2493206
[11] Rizk, M.E., Mahmood, F., Lehtonen, M., Badran, E.A., Abdel-Rahman, M.H. (2016). Induced voltages on overhead line by return strokes to grounded wind tower considering horizontally stratified ground. IEEE Transactions on Electromagnetic Compatibility, 58(6): 1728-1738. https://doi.org/10.1109/TEMC.2016.2588000
[12] Araneo, R., Celozzi, S. (2016). Transient behavior of wind towers grounding systems under lightning strikes. International Journal of Energy and Environmental Engineering, 7: 235-247. https://doi.org/10.1007/s40095-015-0196-7
[13] Pham, C.T. (2021). Evaluation of the influence of lightning overvoltage propagation on wind turbine and mix overhead-cable lines. Journal of Science and Technology, 31(2): 108-116. https://doi.org/10.51316/jst.152.ssad.2021.31.2.14
[14] Sun, Q., Yang, L., Li, Q., et al. (2021). Surge analysis for lightning strike on overhead lines of wind farm. Electric Power Systems Research, 194: 107066. https://doi.org/10.1016/j.epsr.2021.107066
[15] Nam, N.V., Top, T.V., Huy, N.Q. (2022). Factors affecting to transient overvoltage due to direct lightning strikes on wind turbine blade in wind farm. Journal of Sciences and Technology (Hanoi University of Industry), 58(4): 16-21. https://jst-haui.vn/media/30/uffile-upload-no-title30858.pdf.
[16] Stucke, I., Morgenstern, D., Diendorfer, G., et al. (2023). Upward lightning at wind turbines: Risk assessment from larger-scale meteorology. Journal of Geophysical Research: Atmospheres, 129(1): e2023JD039505. https://doi.org/10.1029/2023JD039505
[17] Cai, L., Fan, W., Yan, R., et al. (2024). Electromagnetic fields of lightning return stroke to wind turbines with discontinuous impedance model. Electric Power Systems Research, 233: 110515. https://doi.org/10.1016/j.epsr.2024.110515
[18] Da Silva, W.C., de Azevedo, W.L.M., D’Annibale, J.L.A., de Araújo, A.R.J., Filho, J.P. (2025). Lightning analysis of grounding system in wind farms: Comparison between interconnected overhead ground wire and underground bare cables. Electric Power Systems Research, 242: 111417. https://doi.org/10.1016/j.epsr.2025.111417
[19] Berenger, J.P. (1996). Three-dimensional perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 127: 363-379. https://doi.org/10.1006/jcph.1996.0181
[20] Taflove, A., Hagness, S.C. (2000). Computational Electrodynamics: The Finite-Difference Time-Domain Method (2nd ed.). Artech House.
[21] Arzag, K., Azzouz, Z., Ghemri, B. (2016). Lightning electromagnetic pulse simulation using 3D-FDTD method (comparison between PEC and UPML boundary conditions). In Proceedings of the International Conference on Lightning Protection, Estoril, Portugal, pp. 1-6. https://doi.org/10.1109/ICLP.2016.7791354
[22] Arzag, K., Azzouz, Z., Ghemri, B. (2016). 3D-FDTD computation of lightning return stroke current and associated electromagnetic field using electromagnetic models. International Review of Electrical Engineering, 11(5): 517-525. https://doi.org/10.15866/iree.v11i5.9697
[23] Uman, M.A., McLain, D.K., Krider, E.P. (1975). The electromagnetic radiation from a finite antenna. American Journal of Physics, 43: 33-38. https://doi.org/10.1119/1.10027
[24] Baba, Y., Rakov, V.A. (2005). Lightning electromagnetic environment in the presence of a tall grounded strike object. Journal of Geophysical Research, 110(D9): 1-18. https://doi.org/10.1029/2004JD005505
[25] Courant, R., Friedrichs, K., Lewy, H. (1967). On the partial difference equations of mathematical physics. IBM Journal of Research and Development, 11(2): 215-234. https://doi.org/10.1147/rd.112.0215
[26] Rachidi, F., Nucci, C.A., Ianoz, M., Mazzetti, C. (2002). Influence of a lossy ground on lightning-induced voltages on overhead lines. IEEE Transactions on Electromagnetic Compatibility, 38(3): 250-264. https://doi.org/10.1109/15.536054
[27] Rubinstein, M. (1996). An approximate formula for the calculation of the horizontal electric field from lightning at close, intermediate, and long range. IEEE Transactions on Electromagnetic Compatibility, 38(3): 531-535. https://doi.org/10.1109/15.536087
[28] Zeddam, A., Degauque, P. (2017). Current and voltage induced on telecommunication cables by a lightning stroke. In Lightning Electromagnetics, Routledge, pp. 377-400. https://doi.org/10.1080/02726348708908197
[29] Heidler, F. (1985). Travelling current source model for LEMP calculation. In 6th Symposium and Technical Exhibition on Electromagnetic Compatibility, Zurich, Switzerland, pp. 157-162. https://doi.org/10.23919/EMC.1985.10798870