An algorithm for evaluation of lightning electromagnetic fields at difference distances with respect to lightning channel

The analytical field expressions are proposed to estimate the electromagnetic fields associated with vertical lightning channel at an observation point directly in the time domain using current and geometrical parameters. The proposed expressions embrace a large number of channel base current functions and widely used engineering current models. Moreover, they can be joined with different coupling models for evaluation of lightning induced voltage when they provide different required field components directly in the time domain. The data from measured fields and simulated fields from previous works were employed to validate the proposed field expressions. The simulated results were compatible with both computed field values from previous studies and the measured fields.


Introduction
Lightning is a natural phenomenon that can give effect on the power lines in two ways, that is, direct and indirect.In the direct effects, lightning strikes to power lines or towers or substations directly.On the other hand, the indirect effect considers electromagnetic fields associated with lightning channel and induced voltage due to coupling between generated electromagnetic fields and power system when lightning strikes to the ground or any objective near to power lines or substations [1,2].This study considers the evaluation of electromagnetic fields due to lightning channel.Several methods consider the calculation of electromagnetic fields due to lightning channel while they are validated using measured fields at one or a number of observation points [1] while the most common calculation methods can be concluded by (i) the MONOPOLE method [3,4], (ii) the DIPOLE method [1,3], (iii) the FDTD method [5], (iv) the FDTD-HYBRID methods [6,7].
It should be mentioned that the electromagnetic fields associated with the lightning channel rely basically on some geometrical, channel, and ground parameters.The basic input parameters in almost all calculation methods in perfect ground conductivity case are (i) the channel base current that is usually prepared with direct measurement [2] from the triggered lightning technique or obtained from inverse procedure algorithms from measured electromagnetic fields [8,9], (ii) the current model for prediction of return stroke current wave shape along lightning channel, (iii) the radial distance from lightning channel, (iv) the observation point height with respect to ground surface, (v) the return stroke current velocity along lightning channel.Two first calculation methods are expressed in the cylindrical [1] and spherical [3,4] domains using current and charge density [4] parameters along the lightning channel while they have complexity in the calculation of induction partial when some realistic channel base current functions are used (as Heidler function) and they do not have an analytical solution of integral to time [7,10].On the other hand, the FDTD method just is validated for closed distances from lightning channel [5].Moreover, the FDTD HYBRID methods need to estimate magnetic fields at six points around observation point for calculation of electric fields at an observation point [6,7,11].In order to provide the general analytical electromagnetic field expressions in the time domain, the widely used channel base current functions are introduced and categorized.Therefore, after generalization of current functions based on proposed expressions, the general electromagnetic field expressions due to vertical lightning channel on the prefect ground are proposed.These equations can be applied to a large number of leading current functions and their combinations and can support different generalized engineering current models directly in the time domain.The basic assumptions of this study are the following.
(i) The ground conductivity is infinite.
(ii) The ground surface is flat.
(iii) The lightning channel is perpendicular to the ground and has no branches.
(iv) The corona effect is negligible.
In the next section, the widely used current functions and current models are reviewed.Therefore, the proposed electromagnetic field expressions are presented.Finally, a justification of the proposed method and its validation by measuring fields and simulated fields from previous studies is given.

Channel Base Current Functions
Previous studies presented several channel base current functions to simulate return stroke current wave shape at channel base on the ground surface where the fixed coefficients of the current functions were determined by the data from measured current or measured fields using inverse procedure algorithms.The significant channel base current functions are tabulated in Table 1 where ,   2) by changing the constant coefficients or the combination of (1) and (2).Hence, this study adapts these two basic equations and considers the widely used current functions using a combination of them: where ] . (3)

Return Stroke Current Models
One of the effective factors on electromagnetic field calculations is the behavior of return stroke currents at different heights along the lightning channel that is modeled via current models.Generally, the return stroke current models can be classified into four groups as follows [12,13]: (i) the gas-dynamic or physical models, (ii) the current distributed models, (iii) the electromagnetic models, (iv) the engineering models.
To develop general expressions for electromagnetic fields due to lightning channel, the engineering models are applied in this study.In the engineering models, the current at different heights along the lightning channel can be expressed by a special function based on the channel base current and an attenuation height dependent factor.Most of the engineering models can be presented by a general equation such as (4) [12,13].The internal parameters of (4) for different engineering current models are listed in Table 2 [1,12] where  is considered as the speed of light in free space,  as a decay constant (1∼2 km), and   as the cloud height with respect to the ground surface [1]: where   is the temporary charge height along lightning channel, (  , ) is current distribution along the lightning channel at any height   and any time , (0, ) is channel base current, (  ) is the attenuation height dependent factor, V is the current-wave propagation velocity, V  is the upward propagating front velocity, and  is the Heaviside function defended as Based on (4), return stroke current at the height of   along the lightning channel can be estimated by using channel base current and the height dependent attenuation factor.Real channel

Image channel
The perfect ground Note that, in recent years, several models based on (4) and using measured electromagnetic fields and inverse procedure algorithms were presented [14][15][16].Similarly, different current models can be examined by the comparison between simulated electromagnetic fields using current model and measured fields at different observation points [1].In this model the assumed return stroke velocity is set at a constant value that is typically between /2 and 2/3 [17].

Electromagnetic Fields Associated with Lightning Channel
In this section, electromagnetic field expressions due to two expressed functions in ( 1) and ( 2) at an observation point above the ground surface are estimated whereas the current model is set on (4).According to the geometry of the situation (see Figure 1), the potential vector ( ⃗ ) at the observation point above ground surface can be obtained from ( 6) when Maxwell's equation [3,18] and Lorentz gauge [19] are applied: where ⃗  is the current density,  0 = 4 × 10 −7 V⋅s/(A⋅m), and  0 = 8.85 × 10 −12 F/m.Moreover, Figure 1 demonstrates that the solution to (6) can be obtained from (7) whereas the infinitesimal current source is located along -axis at space in position  →  2 from the origin.Also, the observation point is located at space and  →  1 from the origin point; therefore, ⃗  can be estimated by where ⃗  is line current density and   is unit vector along lightning channel.
Likewise, by replacing return stroke current and   with ⃗  and unit vector, respectively, in (7) can convert it to (8) where Further, the magnetic flux density and electric fields can be obtained from potential vector ⃗  as expressed in ( 10) and (11), respectively [3,18]: where ⃗  is the magnetic flux density and ⃗  is the electric field vector.
Therefore, based on the general current model in (4), the field components for first current function as presented in (1) can be obtained from (12), when (8), (10), (11), the trapezoid, and FDTD methods are applied [20][21][22].Note that the internal parameters are collected from Table 3.Also, ⃗   RS1 (, , ,   ) = 0: where ⃗   RS1 is the magnetic flux density at -axis due to return stroke current function from (1), ⃗   RS1 is the magnetic flux density at -axis due to return stroke current function from (1), ⃗   RS1 is the magnetic flux density at -axis due to return stroke current function from (1), ⃗   RS1 is the electric field at -axis due to return stroke current function from (1), ⃗   RS1 is the electric field at -axis due to return stroke current function from (1), and ⃗   RS1 is the electric field at -axis due to return stroke current function from (1).Δ is the time step,  is the number of time steps, Δℎ   for others is division factor (≥ 2), In addition, by using the same methodology that was used for the first current function, the components of electromagnetic field for the second current function as expressed in (2) are proposed by (15).Note that ⃗   RS2 (, , ,   ) = 0: where ⃗   RS2 is the magnetic flux density at -axis due to return stroke current function from (2), ⃗   RS2 is the magnetic flux density at -axis due to return stroke current function from (2), ⃗   RS2 is the magnetic flux density at -axis due to return stroke current function from (2), ⃗   RS2 is the electric field at -axis due to return stroke current function from (2), ⃗   RS2 is the electric field at -axis due to return stroke current function from (2), and ⃗   RS2 is the electric field at -axis due to return stroke current function from (2).
Consequently, due to the whole return stroke current functions from Table 1, the electromagnetic fields can be easily estimated via the proposed field expressions whereas the total fields due to two part functions are equal to the sum of the electromagnetic fields due to each part independently.On top of that the first part current functions in Table 3 should be converted to either (1) or (2) format before getting into proposed expressions.In other words, the proposed algorithm can support the whole current models related to (4).In the same way, the components of the electromagnetic fields in the Cartesian coordinates can be converted to   Mathematical Problems in Engineering 7 cylindrical domain using ( 16) and ( 17) for the magnetic flux density and the electric field components, respectively: where  = arc cos (/√ 2 +  2 ), ⃗   (, , , ) / ⃗   (, , , ) is the magnetic flux density/electric field at -direction (horizontal), ⃗   (, , , ) / ⃗   (, , , ) is the magnetic flux density/electric field at -direction, and ⃗   (, , , ) / ⃗   (, , , ) is the magnetic flux density/electric field at -direction (vertical).

Result and Discussion
This section deals with the validation of the electromagnetic fields due to lightning channel using some channel base current functions.According to the proposed method, the widely used current functions can be converted into sum of individual parts using (1) and (2).Therefore, the electromagnetic fields due to each part can be calculated by ( 12) and (15).Then, the values of total electromagnetic field are equal to the sum of individual components of electromagnetic field.When the MTLE model is applied as current model, the simulated magnetic flux density and the vertical electric field caused by a sum of two Heidler function as a channel base current function at an observation point with 4.6 km distance from lightning channel are depicted in Figures 2  and 3, respectively.The measured magnetic flux density and vertical electric field with similar primary conditions as the ones in Figures 2 and 3 are demonstrated in Figures 4 and 5, respectively.Therefore, comparisons between the measured and simulated fields for magnetic flux density and vertical eclectic field are tabulated in Tables 4 and 5, respectively.It should be noted that the current parameters are  01 = 17 kA,  02 = 8 kA, Γ 11 = 0.4 × 10 −6 , Γ 12 = 4 × 10 −6 , Γ 21 = 4 × 10 −6 , Γ 22 = 50 × 10 −6 ,  1 = 2,  2 = 2, V = 1 × 10 8 , and  = 1500 [23].Additionally, the sum of two Heidler functions as a channel base current function can be expressed by the sum of two general current functions based on (2).Therefore, the total electromagnetic fields are estimated by using the proposed method and considering each current part separately.The comparison between simulated and measured fields shows the overall good agreement with experimental measurement for proposed field expressions.Similarly, the simulated fields   [23,24].are in accord with other previous method of calculation in [23,24].
The proposed method is also validated by Nucci current function in Table 1 which is defined by the combination of ( 1) and ( 2) when the MTLL model is applied as a current model.The simulated magnetic flux density and vertical electric field are formed as in Figures 6 and 7, respectively.By the same token, Figure 8 depicts the measured vertical electric field.
In addition, the numerical comparison between the measured fields and simulated vertical fields is listed in Table 6 suggesting that the simulated field is in agreement with the corresponding measured field.Moreover, magnetic flux densities simulated using the proposed field expressions are compatible with the simulated field calculated using other methods given in the reference [5].Note that the current parameters are  01 = 3.25 kA,  02 = 8.95 kA, Γ 11 = 0.072 × 10 −6 , Γ 12 = 16.67 × 10 −6 , Γ 21 = 100 × 10 −6 , Γ 22 = 0.5 × 10 −6 ,  1 = 2, and V = 1.5 × 10 8 m/s [5].The proposed expressions can directly be applied to widely used current functions in the time domain just by substituting the current and geometrical parameters in the proposed field expressions without needing to apply any extra algorithms and extra conversions (such as Fourier transform) when the realistic current functions are applied [1].Besides, the proposed algorithm expressions can   support different engineering current models, provided that they follow (4).The proposed method can corroborate most of the coupling models that compute lightning induced overvoltage (LIOV), directly in the time domain, because not only are all of the electromagnetic field components accessible but also the fundamental equations are in the Cartesian coordinates.This adds a new dimension to the development of LIOV coupling models which are restricted, in the present context, to apply a fast and direct computational block in the time domain for estimation of electromagnetic fields generated by lightning channel without needing to use extra computational stage as Fourier transform for a large number of channel base current functions and engineering current models.Likewise, the proposed field expressions can be used in the inverse procedure algorithms to evaluate lightning return stroke current using measured electromagnetic fields directly in the time domain [8,9].Moreover, the proposed field expressions can be developed for the case of nonperfect ground conductivity using Cooray-Rubinstein equation [25].
Consequently, the evaluation of electromagnetic fields at different points along the power line can easily be computed with different values of  or  [25][26][27][28].The accuracy of the results can be increased by changing the values of  and Δ.However, the processing time will increase with such increasing of .Further, the proposed field expressions can be used for close and medium distances from the lightning channel; however, some of the calculation methods in the time domain are merely confined to close distances from the lightning channel.

Conclusion
In this study, the general electromagnetic field expressions due to vertical lightning channel at an observation point are proposed directly in the time domain just by substituting of current and geometrical parameters in the proposed field expressions without needing to apply any algorithms and extra conversions.The proposed field expressions can support the widely used current functions and current models.Besides, the proposed method was validated with two return stroke current samples.The results showed that there is a good agreement between the electromagnetic fields computed by proposed method and both the measured fields and those calculated by previous simulation methods.In addition, the proposed method is applicable for both close and intermediate distances from the lightning channel and it is completely in accordance with many coupling models in evaluating the induced lightning voltage and also inverse procedure algorithms in the time domain.This would be very useful especially for the utility in assessing the lightning performance of an overhead line, particularly the distribution line.

Figure 1 :
Figure 1: Geometry of the problem.

Table 3 :
The internal parameters of proposed electromagnetic fields expressions.

Table 4 :
Numerical comparison between the simulated magnetic flux density and the corresponding measured fields from Figures2 and 4after subtracting initial decay time.

Table 5 :
Numerical comparison between the simulated vertical electric field and the corresponding measured fields from Figures3 and 5after subtracting initial decay time.

Table 6 :
Numerical comparison between the simulated magnetic flux density and the corresponding measured fields from Figures7 and 8after subtracting initial decay time.