Instruction FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 30, N o 2, June 2017, pp. 245 - 256 DOI: 10.2298/FUEE1702245L NOVEL APPROACH TO MODELLING OF LIGHTNING CURRENT DERIVATIVE  Karl Lundengård 1 , Milica Rančić 1 , Vesna Javor 2 , Sergei Silvestrov 1 1 Mälardalen University, UKK, Division of Applied Mathematics, Västerås, Sweden 2 University of Niš, Faculty of Electronic Eng., Dept. of Power Engineering, Niš, Serbia Abstract. A new approach to mathematical modelling of lightning current derivative is proposed in this paper. It builds on the methodology, previously developed by the authors, for representing lightning currents and electrostatic discharge (ESD) currents waveshapes. It considers usage of a multi-peaked form of the analytically extended function (AEF) for approximation of current derivative waveshapes. The AEF function parameters are estimated using the Marquardt least-squares method (MLSM), and the framework for fitting the multi- peaked AEF to a waveshape with an arbitrary number of peaks is briefly described. This procedure is validated performing a few numerical experiments, including fitting the AEF to single- and multi-peaked waveshapes corresponding to measured current derivatives. Key words: analytically extended function, lightning current derivative, lightning current function, lightning stroke, Marquardt least-squares method 1. INTRODUCTION Besides different parameters of lightning electromagnetic field and lightning discharge currents, which greatly endanger the functionality of power systems, electrical equipment and electronic devices, lightning current derivative signal is often measured at tall instrumented towers, towers at elevated terrain and at rocket-triggered stations, [1]-[8]. Current derivatives approximation is important for calculation of lightning induced overvoltages and for further improvements of lightning discharge models [1], [5], [9]. Generalizing the function for representing lightning currents from [10]-[12], the proposed multi-peaked analytically extended function (AEF) has been applied by the authors to modelling of different lightning currents, including those defined in the IEC Standard 62305-1 [13], slow and fast-decaying ones, as well as measured ones, see e.g. [14]-[16]. Furthermore, it has been recently used in [9] and [17] for representation of the electrostatic discharge (ESD) current corresponding to the IEC Standard 61000-4-2 waveshape as given Received August 29, 2016; received in revised form November 26, 2016 Corresponding author: Vesna Javor University of Niš, Faculty of Electronic Engineering, Aleksandra Medvedeva 14, 18000 Niš, Serbia (E-mail: vesna.javor@elfak.ni.ac.rs) 246 K. LUNDENGÅRD, M. RANĈIĆ, V. JAVOR, S. SILVESTROV in [18]-[19]. The AEF’s parameters were fitted to the desired current waveshapes using the Marquardt least-squares method (MLSM), [20]. In this paper we explore the possibility of reproducing the waveshape of the lightning current first derivative using the AEF and adjusting its non-linear parameters employing the MLSM. The validity of approximation and this methodology is tested by performing a few numerical experiments related to modelling of lightning current derivative signals measured at the CN Tower [5]. Since installation, simultaneous measurements of currents and current derivatives by Rogowski coils, corresponding electromagnetic field values detected by sensors and high-speed cameras at a few km distance from the tower have been providing useful data for analysis, [1]-[5]. Reflection coefficients are estimated for the CN Tower and employed for magnetic field calculation in [5]. Reflections occur from the tip of this tower, top and bottom of its restaurant and from the ground, so as at the upward-propagating lightning return- stroke channel front, and produce peaks in the current derivative waveshape. In this paper, lightning current derivative approximation is done taking into account the initial peak and subsequent peaks in the derivative waveshape, regardless of their cause. The same procedure may be used in the case when measured current derivatives have multi- peaked waveshapes for other reasons, e.g. due to various configurations of the terrain and some tall structures, or due to lightning current channel discontinuities and branching. 2. MODELLING OF THE LIGHTING CURRENT DERIVATIVE 2.1. Analytically extended function (AEF) and some of its properties The basic building block of the multi-peaked AEF is, as referred to in [18], the power exponential function (PEF) given by 1 ( ; ) ( ) , 0 , t x t te t      (1) where the β-parameter determines the steepness of both its rising and its decaying part. The AEF is constructed as a function consisting of piecewise linear combinations of PEFs that have been scaled and translated to ensure that the resulting function is continuous. In [18], it is defined as 1 1 1 , , 1 1 , , 1 1 ( ), ,1 , d ( ) d ( ), , 1, q k q q q p k p nq Dm Dm q k q k m m k k np Dm q k q k m k k I I x t t t t q p i t t I x t t t q p                            (2) where:  1 2 , ,..., pDm Dm Dm I I I - the difference in height between each pair of peaks,  1 2 , ,..., pm m m t t t - the times corresponding to these peaks,  0 q n  - the number of terms in each time interval,  ,q k  - real values so that ,1 1 qn q kk    , and  , ( ) q k x t - PEFs defined by ,q k  parameters in the following way: Novel Approach to Modelling of Lighting Current Derivative 247 2 , 1 2 , 1 , exp 1 , ( ) exp 1 1, q k q q q q q k q q m m m m q k m m t t t t q p t t x t t t q p t t                                      (3) where 1q q qm m m t t t     . Expression (2) can be written more compactly as     1 1 T 1 T 1 , , 1 , d ( ) d , , 1, k q q q k p q Dm Dm q q m m k p D m q q m k I I t t t t q p i t t I t t t q p                     x x   (4) after introducing T ,1 ,2 , [ ] , qq q q q n     ,1 ,2 , ( ) [ ( ) ( ) ( )]. qq q q q n t x t x t x tx The first derivative of the multi-peaked AEF corresponds to the second derivative of the lightning discharge current, i(t), and can be easily found since the AEF consists of elementary functions. Compact form is given by     1 1 T 2 2 T ( ) , , 1 , d ( ) d ( ) , , 1, q q q q q q q q p q q m q Dm q q q m m m m m q Dm q q q m m m t t x t I t t t t q p t t t i t t t t x t I t t t q p t t                  η B x η B x (5) where q B are diagonal matrices: 2 2 2 ,1 ,2 , 2 2 2 ,1 ,2 , diag( 1, 1, , 1), 1 , diag( , , , ), 1. q q q q q n q q q q n q p q p                    B Based on this expression, it is easy to see that the current’s second derivative is also continuous since it will be zero at each qm t . The integral of the AEF corresponds to the lighting discharge current i(t) and is also relatively straightforward to find, since the integral of the PEF can be written using the lower incomplete gamma function ([21]) i.e. 248 K. LUNDENGÅRD, M. RANĈIĆ, V. JAVOR, S. SILVESTROV 1 0 1 0 1 01 1 ( ; ) d ( ( 1, ) ( 1, )) ( , , ) t t e e x t t t t t t                    , (6) where 1 0 ( , ) d t t e        is the incomplete gamma function. Combining (6) and (2) we obtain the integral of the rising part of the AEF 1 1 , 1 1 11 2 , , 1 1 1 1 , 1 1 d ( ) d ( ) ( , ) d ˆ ( 1) ( ) ( , ) , 0 , a b a k a a a q q k q q b k b b p na t m a Dm Dm a k a m a t k k nqb m Dm Dm q k q k q a k k nb b m Dm Dm b k b b m a b m m k k i t t t t I I g t t t t I I g t t I I g t t t t t t                                                          1 1 , , a a b ba m m b m t t t t t       (7) with 2 , 1 1 2 , 1 1 02 1 0 ,22 , ( , ) 1, , ( 1) q k q q q k q q m m q q k m mq k t t t te g t t t t                     and 1 ˆ ( ) ( 1, ) e g         . The integration formula corresponding to the decaying part is 1 1 0 1, 1 1 0 0 1 1 1 d ( ) d ( , ), d p k p np t Dm p k p mt k k i t t I g t t t t t t             , i.e. 1 2 1, 1, 1 1 d ( ) d ( ) d p k m p np Dm p k p kt k k i t t I g t           , (8) with 1 ( ) ( ( 1) ( 1, )) e g            and 1 0 ( ) de         , the Gamma function, [21]. 2.2. Marquardt least-squares method (MLSM) Detailed explanation of the MLSM algorithm is given in [15], [16], here we just go over the parts specific for the multi-peaked AEF. The MLSM is used for estimating β-parameters, and from these, the corresponding η– parameters are calculated. In each iteration step, η–parameters are obtained using the regular least-square method since for fixed β-parameters the AEF is linear in η. Based on these η–parameters, a new set of β-parameters is found. The MLSM uses a Jacobian matrix, denoted by J, containing partial derivatives of the residuals. The least square fitting of the multi-peaked AEF to a set of data points can be done separately between each peak (and after the final one), and the corresponding J matrix is ,1 ,1 ,2 ,1 , ,1 ,1 ,2 ,2 ,2 , ,2 ,1 , ,2 , , , ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) q q q q q q q q q q q n q q q q q q n q q q k q q k q n q k p t p t p t p t p t p t p t p t p t                J , (9) Novel Approach to Modelling of Lighting Current Derivative 249 where q k is the number of data points between the (q-1)th and qth peak, and ,q r t is the time corresponding to these data points, and 2 , , , , , 2 , , , , , d ( ) 2 ( ) ( 1), 1 , d ( ) ( ) ( ), 1, q q Dm q k q k q q k q r q r q r Dm q k q k q q k t t q r i t I h t x q p t p t I h t x q p                          with 1 1 ln 1, 1 , ( ) ln 1, 1. q q q q q q m m m m q m m t t t t q p t t h t t t q p t t                                  3. AEF REPRESENTING MEASURED LIGHTING CURRENT DERIVATIVE - EXAMPLES In this section we validate our model by attempting to represent measured lightning current derivatives data obtained at the 553m CN Tower, Toronto, Canada, [5]. Time and current values corresponding to AEF peaks were chosen manually and the rest of the AEF parameters were obtained using the framework briefly described in Section 2.2. The number of time intervals and terms in each of them vary from example to example. General notation, AEFp(n1, …, np) for nq, q=1, …, p, is used to denote an AEF with p peaks and chosen number of terms nq in each time interval q. 3.1. Single-peaked waveshape The first example illustrates the application of a single-peaked AEF to representation of the measured initial current derivative impulse occurring in the first 0.5 s given in [5, Fig. 4]. The best fitting was obtained choosing two terms in each of the two time intervals: 0-tm and tm-0.5 s (the moment tm corresponds to the maximum current derivative). Current derivative value at t=0 is treated as the first point of approximation, so there are 4 terms in total, for these 2 intervals. Obtained AEF2(2,2) model is illustrated in Fig. 1a along with the measured data, data points used for the MLSM fitting, and the locations of peaks observed in this waveshape. Using the expressions (7) and (8) we also obtained the AEF’s integral, i.e. the lighting discharge current. It can be observed in Fig. 1b along with the numerically integrated measured data. 250 K. LUNDENGÅRD, M. RANĈIĆ, V. JAVOR, S. SILVESTROV a) b) Fig. 1 a) AEF2(2,2) representing measured lightning current derivative from [5, Fig. 4], and b) the corresponding lightning current 3.2. Multi-peaked waveshapes In this part we attempt modeling of the measured current derivatives data that include the initial and a number of subsequent impulses. The recoded waveshapes have great number of peaks and therefore is harder to model them using standard functions, but these are more suitable for modelling by the multi-peaked AEF. The first example corresponds to an event of lightning discharge measured at the CN Tower, using the Rogowski coil positioned at 474 m, illustrated in [5, Fig. 2] in 10 s. Such current derivative waveshape corresponds to typical fast-rising negative lightning discharge which occurs in about 80% of the registered cases (in 126 flashes out of 160 [5]). The complexity of the AEF used for modelling of such multi-peaked waveshapes depends on the desired level of accuracy of the data representation. Novel Approach to Modelling of Lighting Current Derivative 251 In Fig. 2 are presented two AEFs with different number of peaks, including the starting current derivative value at t = 0 and other peaks which are chosen such that they correspond to local maxima only: a) AEF6(1,2,2,2,2,2) with 6 intervals and 11 terms in total, and b) AEF8(1,2,2,2,2,2,2,2) with 8 intervals and 15 terms in total. The increased number of time intervals fixes representation of the waveshape part corresponding to the period between the fourth and fifth peak of AEF6, and also after its sixth peak, so that the total number of intervals in AEF8 is increased by 2, whereas the number of terms by 4. a) b) Fig. 2 Multi-peaked AEFs (using starting point and maxima only) representing measured lightning current derivative from [5, Fig. 2]: a) AEF6(1,2,2,2,2,2) with 6 peaks, b) AEF8(1,2,2,2,2,2,2,2) with 8 peaks 252 K. LUNDENGÅRD, M. RANĈIĆ, V. JAVOR, S. SILVESTROV Additional improvement is needed and could be achieved by further segmentation and including also local minima, so as by increasing the number of terms. Two such AEF models are illustrated in Figs. 3a and 3b, both with thirteen peaks, but for different number of terms chosen to represent some of its intervals. Thirteen peaks in AEF13 include 4 minima added to AEF8 and also one more maximum at its ending part, so that the number of peaks is increased from 8 (in Fig.2b) to 13 (in Figs. 3a and 3b). These two AEFs are denoted by a) AEF13a(1,1,1,1,1,2,1,1,1,1,2,2,1) with 13 intervals and 16 terms in total, and b) AEF13b(1,1,2,1,2,2,1,2,1,2,2,2,2) with 13 intervals and 21 terms in total, where the bold numbers in brackets point out to the changed number of terms, in some intervals increased from 1 to 2. a) b) Fig. 3 Multi-peaked AEFs with 13 peaks (using starting point, 8 maxima and 4 minima) representing measured lightning current derivative from [5, Fig. 2]: a) AEF13a(1,1,1,1,1,2,1,1,1,1,2,2,1), b) AEF13b (1,1,2,1,2,2,1,2,1,2,2,2,2) Novel Approach to Modelling of Lighting Current Derivative 253 Results for the same lightning current derivative measured at CN Tower are given in first 7s in Figs. 4a and 4b for fitting by AEFs corresponding to data from [5, Fig. 6]. Model AEF7(1,2,2,2,2,2,2) with 7 peaks (starting point and maxima only) and 13 terms is presented in Fig. 4a, able to capture the initial impulse and subsequent peaks due to reflections at the tower discontinuities. AEF7 has one more peak added at the end of AEF6, and 2 more terms in total. In Fig. 4b, AEF13c(1,2,2,2,2,2,2,2,2,1,1,2,2) model is presented with the total of 13 peaks (the starting point, 8 maxima and 4 minima), which almost perfectly models measured set of data. It has 23 terms in total, 4 terms added and 2 excluded compared to AEF13b. The difference between those two is that 1 peak was added for AEF13c between tenth and eleventh peak of AEF13b, which improved significantly the approximation, and the thirteenth peak was excluded from the end of AEF13b. a) b) Fig. 4 Multi-peaked AEF representing measured current derivative from [5, Fig. 6]: a) AEF7(1,2,2,2,2,2,2) with 7 peaks (using starting point and maxima only), b) AEF13c(1,2,2,2,2,2,2,2,2,1,1,2,2) with 13 peaks (starting point, 8 maxima & 4 minima) 254 K. LUNDENGÅRD, M. RANĈIĆ, V. JAVOR, S. SILVESTROV Figure 5 illustrates lighting discharge currents corresponding to above modelled multi- peaked current derivative waveshapes. Again, expressions (7) and (8) are employed to calculate them, and the numerically integrated measured data is also given for comparison. Fig. 5a corresponds to AEF13b(1,1,2,1,2,2,1,2,1,2,2,2,2) model shown in Fig. 3b, while Fig. 5b relates to model AEF13c(1,2,2,2,2,2,2,2,2,1,1,2,2) from Fig. 4b. a) b) Fig. 5 Lightning currents corresponding to derivatives modelled by AEFs: a) AEF13b (1,1,2,1,2,2,1,2,1,2,2,2,2) from Fig. 3b, b) AEF13c(1,2,2,2,2,2,2,2,2,1,1,2,2) from Fig. 4b Novel Approach to Modelling of Lighting Current Derivative 255 4. CONCLUSIONS Approximation of lightning current derivatives is needed for calculation of lightning induced effects and for improvements of lightning discharge models. Suitability of the multi- peaked AEF to represent lightning current derivatives is presented in this paper through a few examples. AEF’s non-linear parameters are calculated using Marquardt least-squares method (MLSM), so that the measured current derivatives signals [5] are well approximated. The approximation by AEFs in this paper is done for single- and multi-peaked current derivative waveshapes. Increasing the number of maxima and minima, so as the number of terms in total, improves the approximation of the current derivative by AEF. The lightning current waveshape is obtained with great accuracy as analytically integrated AEF representation of the measured derivative. Multi-peaked lightning current derivatives are characteristic for lightning discharges to tall towers and high structures at elevated terrain, but also for subsequent lightning strokes and lightning current channels with discontinuities and branching. Further work should be aimed at including such current and its derivative function into lightning stroke models in order to obtain measured lightning electromagnetic field at certain distances. REFERENCES [1] K. Elrodesly and A. M. Hussein, "CN Tower Lightning Return-Stroke Current Simulation”, Journal of Lightning Research, vol. 4, Suppl 2: M3, pp. 60-70, 2012. [2] A. M. Hussein, M. Milewski, and W. Janischewskyj, "Correlating the Characteristics of the CN Tower, Lightning Return-Stroke Current with those of its Generated Electromagnetic Pulse”, IEEE Transactions on Electromagnetic Compatibility, vol. 50, no. 3, pp. 642-650, Aug. 2008. [3] A. M. Hussein, M. Milewski, E. Burnazovic and W. Janischewskyj, "Current Waveform Characteristics of CN Tower Negative and Positive Lightning", In Proceedings of the X International Symposium on Lightning Protection, Curitiba, Brazil, 2009, pp. 451-456. [4] B. Kordi, R. Moini, W. Janischewskyj, A. M. Hussein, V. O. Shostak and V. A. Rakov, " Application of the antenna theory model to a tall tower struck by lightning”, Journal of Geophysical Research, vol. 108, no. D17, 4542, doi: 10.1029/2003JD003398, 2003. [5] M. Milewski and A. M. Hussein, "Tall-Structure Lightning Return-Stroke Modelling", In Proceedings of the 14th International Middle East Power Systems Conference (MEPCON’10), Cairo University, Egypt, 2010, Paper ID 313, pp. 947-952. [6] F. Rachidi, W. Janischewskyj, A. M. Hussein, C. A. Nucci, S. Guerrieri, B. Kordi and J-S. Chang, "Current and Electromagnetic Field Associated with Lightning–Return Strokes to Tall Towers", IEEE Transactions on Electromagnetic Compatibility, vol. 43, no. 3, pp. 356-367, Aug. 2001. [7] V. A. Rakov, “Transient Response of a Tall Object to Lightning”, IEEE Transactions on Electromagnetic Compatibility, vol. 43, no. 4, pp. 654-661, 2001. [8] M. A. Uman, J. Schoene, V. A. Rakov, K. J. Rambo and G. H. Schnetzer, “Correlated Time Derivatives of Current, Electric Field Intensity, and Magnetic Flux Density for Triggered Lightning at 15 m”, Journal of Geophysical Research, vol. 107, no. D13, doi: 10.1029/2000JD000249, 2002. [9] V. Javor, "An Analytically Extended Function for Representing the Lightning Current First Derivative", In Proceedings of the Int. Colloquium on Lightning and Power Systems, Bologna, Italy, 2016, P13_S3.2, pp. 1-8. [10] V. Javor and P. D. Rancic, “A Channel-Base Current Function for Lightning Return-Stroke Modeling”, IEEE Transactions on Electromagnetic Compatibility, vol. 53, no. 1, pp. 245-249, Feb. 2011. [11] V. Javor, "Multi-Peaked Functions for Representation of Lightning Channel-Base Currents", In Proceedings of 2012 International conference on lightning protection - ICLP, Vienna, Austria, 2012, pp. 1–4. 256 K. LUNDENGÅRD, M. RANĈIĆ, V. JAVOR, S. SILVESTROV [12] V. Javor, "New Function for Representing IEC 61000-4-2 Standard Electrostatic Discharge Current", Facta Universitatis, Series: Electronics and Energetics, vol. 27(4), pp. 509-520, 2014. [13] IEC 62305-1, Protection against lightning - Part I: General Principles Ed. 2.0, 2010-12. [14] K. Lundengård, M. Ranĉić, V. Javor and S. Silvestrov, "Application of the Multi-Peaked Analytically Extended Function to Representation of Some Measured Lightning Currents", Serbian Journal of Electrical Engineering, vol. 13(2), pp. 1-11, 2016. [15] K. Lundengård, M. Ranĉić, V. Javor and S. Silvestrov, "Estimation of Parameters for the Multi-Peaked AEF Current Functions", Methodology and Computing in Applied Probability, Springer, pp. 1-15, 2016. DOI: 10.1007/s11009-016-9501-z [16] K. Lundengård, M. Ranĉić, V. Javor and S. Silvestrov, “On Some Properties of the Multi-Peaked Analytically Extended Function for Approximation of Lightning Discharge Currents”, Engineering Mathematics I: Electromagnetics, Fluid Mechanics, Material Physics and Financial Engineering, Series: Springer Proceedings in Mathematics & Statistics, Vol. 178, Eds. S. Silvestrov and M. Ranĉić, Springer, Heidelberg, 2016, pp. 151-172, eBook ISBN 978-3-319-42082-0; Hardcover ISBN 978-3-319-42081-3; DOI 10.1007/978-3-319-42082-0 [17] K. Lundengård, M. Ranĉić, V. Javor and S. Silvestrov, "Multi-Peaked Analytically Extended Function Representing Electrostatic Discharge (ESD) Currents", in AIP Conference Proceedings of ICNPAA 2016, La Rochelle, France, 2016, pp. 1-10. [18] EMC - Part 4-2: Testing and Measurement Techniques - Electrostatic Discharge Immunity Test. IEC International Standard 61000-4-2, basic EMC publication, 1995+A1:1998+A2:2000. [19] EMC - Part 4-2: Testing and Measurement Techniques - Electrostatic Discharge Immunity Test. IEC International Standard 61000-4-2, basic EMC publication, Ed. 2, 2009. [20] D. M. Marquardt, "An Algorithm for Least-Squares Estimation of Nonlinear Parameters", Journal of the Society for Industrial and Applied Mathematics, vol. 11(2), pp. 431-441, 1963. [21] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. 1964, Dover, New York.