Accelerating the numerical computation of indirect lightning effects by means of vector fitting

In the context of numerical computation of indirect lightning effects it is customary to use volumediscretizing methods in time domain, such as the Finite Difference Time Domain (FDTD) method, the Finite Integration Technique (FIT), or the Transmission Line Matrix (TLM) method. If standard lightning electromagnetic pulses (LEMPs) of tenths of microseconds duration are used as excitations, these methods require long computation times, as implied by the Courant criterion. It is proposed to use shorter pulse forms and to compare the transfer functions obtained by different pulse durations by means of macromodels that are obtained from the vector fitting method. Numerical computation of lightning related transfer functions of a canonical structure indicate that the duration of the exciting pulse can typically be shortened by at least one order of magnitude if compared to a standard pulse.


Introduction
Modern computer-based computation techniques make it possible to simulate lightning related effects on aircraft structures (Apra et al., 2008).These techniques are defined either in time or frequency domain.In case of time domain computations long computation times usually are required.This is a consequence of the Courant criterion which relates the minimum discretisation length of the geometric model to the maximum allowed time step (Hoffman, 2001).Typically, if geometric details in the range of millimeters need to be resolved, the maximum allowed time step will be of the order of picoseconds.Then the simulation of the lightning strike over its whole duration of several tenths of microseconds usually is a time consuming task.
Correspondence to: F. Gronwald (gronwald@tu-harburg. de) In this contribution an adaptive macromodelling technique is introduced which is based on the idea to determine the LEMP transfer function of a system by the use of excitation pulses that are shorter than standard lightning pulses.This method is based on an adaptive stopping criterion which has been used in the context of time domain characterization of microwave components (Deschrijver et al., 2009).In the present case of lightning analysis, several simulations with comparatively short excitation pulses have to be performed.The resulting transfer functions are characterized by a small number of parameters which are obtained from the vector fitting procedure and represent macromodels of the system (Gustavsen and Semlyen, 1999;Gustavsen, 2006).These macromodels can be used to quantitatively compare the transfer functions that are obtained by means of the different excitation pulses.
In principle, transfer functions calculated from different excitations should be the same if they relate to the same system.However, due to numerical inaccuracies of the time domain data and rounding errors the use of short pulses is limited.On the other hand, the use of comparable long excitation pulses will result in long computation times.Therefore a compromise between these two constraints needs to be found.To this end, comparatively short pulses are used at the beginning of an adaptive procedure.The duration of the first pulse can be chosen to be about a hundredth of the duration of a regular LEMP pulse, yielding a first transfer function with an associated macromodel.Then the duration of the excitation pulses is successively increased, yielding further transfer functions with associated macromodels.The differences between the parameters of different macromodels are observed until they fall below a predefined limit and an approximate convergence has been achieved, thus finally yielding an acceptable transfer function.
In the following Sect.2, the methodology of the proposed method is described in more detail.In Sect.3, a numerical example illustrates the method by means of a canonical tank model and Sect. 4 provides the conclusion.
Published by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V.
In lightning analysis, LEMP transfer functions are the main quantities of interest.They relate a lightning current i in (t) as relevant input variable to an observable which often is given by an induced voltage v out (t).In case of numerical lightning analyis in time domain an observable, such as v out (t), is calculated from the excitation i in (t) for a number of discrete points in time.Then a discrete transfer function H D (k) is obtained by Fourier or Laplace transformation according to where k and m denote discrete indices in frequency and time domain, respectively.The vector fitting method allows to approximately fit a rational function to the discrete transfer function (Gustavsen and Semlyen, 1999;Gustavsen, 2006).This rational function can be written as a residue-pole expansion of the form The advantage of H R (s) is that, as a result, in the context of lightning analysis LEMP transfer functions usually are characterized by only a small number of parameters r j , p j , d, and e.As will be seen below, realistic LEMP transfer functions can already be represented by a first order approximation with N = 1.The proposed method can now be illustrated by the block diagram which is shown in Fig. 1.Starting from a comparatively short excitation pulse, the desired observable is obtained from a numerical time-domain calculation.This allows to calculate the desired transfer function H D (k) and a corresponding macromodel which is expressed by H R (s).A problem that is associated to the use of short pulses is that the discrete Fourier transform of such a signal will tend to have a sparse set of samples in the low frequency range.Due to dominant quasistationary contributions of typical lightning currents, the main frequency range of interest is actually located within this sparsely populated range of low frequencies.However, the sparsity makes the determination of a related transfer function susceptible to numerical errors.Therefore it is necessary to check the stability of results by using longer excitation pulses.It is then possible to compare the transfer functions H R (s) that are obtained by excitation pulses of different lengths.If an increase of pulse duration does not lead to a significant change of H R (s) it is assumed that the desired transfer function is obtained with sufficient accuracy.
To make the notion of a short excitation pulse more precise, in Fig. 2 a standard double exponential lightning pulse of the form 2 J. Anatzki and F. Gronwald: Accelerating the Numerical Computation of Indirect Lightning Effects by means of Vector Fitting

Methodology
In lightning analysis, LEMP transfer functions are the main quantities of interest.They relate a lightning current i in (t) as relevant input variable to an observable which often is given by an induced voltage v out (t).In case of numerical lightning analyis in time domain an observable, such as v out (t), is calculated from the excitation i in (t) for a number of discrete points in time.Then a discrete transfer function H D (k) is obtained by Fourier or Laplace transformation according to where k and m denote discrete indices in frequency and time domain, respectively.The vector fitting method allows to approximately fit a rational function to the discrete transfer function (Gustavsen and Semlyen, 1999;Gustavsen, 2006).This rational function can be written as a residue-pole expansion of the form The advantage of H R (s) is that, as a result, in the context of lightning analysis LEMP transfer functions usually are characterized by only a small number of parameters r j , p j , d, and e.As will be seen below, realistic LEMP transfer functions can already be represented by a first order approximation with N = 1.
The proposed method can now be illustrated by the block diagram which is shown in Fig. 1.Starting from a comparatively short excitation pulse, the desired observable is obtained from a numerical time-domain calculation.This allows to calculate the desired transfer function H D (k) and a corresponding macromodel which is expressed by H R (s).A problem that is associated to the use of short pulses is that the discrete Fourier transform of such a signal will tend to have a sparse set of samples in the low frequency range.Due to dominant quasistationary contributions of typical lightning currents, the main frequency range of interest is actually located within this sparsely populated range of low frequencies.However, the sparsity makes the determination of a related transfer function susceptible to numerical errors.Therefore it is necessary to check the stability of results by using longer excitation pulses.It is then possible to compare the transfer functions H R (s) that are obtained by excitation pulses of different lengths.If an increase of pulse duration does not lead to a significant change of H R (s) it is assumed that the desired transfer function is obtained with sufficient accuracy.
To make the notion of a short excitation pulse more precise, in Fig. 2 a standard double exponential lightning pulse of the form FDTD Simulation using a Gaussian pulse with defined duration t Determine the deviation between old and new model  5412A, 2005), where the time from beginning to peak value 200 kA is 6.4 us, the time from beginning to decay to 50 percent peak value, that is 100 kA, is 69 us.
is contrasted to a unipolar Gaussian pulse of the form The parameters of the pulses are adjusted such that they have a similar rise time but the Gaussian pulse decays much faster if compared to the double exponential pulse.The amplitude spectrum of both pulses is shown in Fig. 3.For the double exponential pulse, the amplitude spectrum is well-known and characterized by dominant low-frequency contributions (Lee, 1995).The amplitude spectrum of the unipolar Gaussian pulse is characterized by dominant low frequency contributions as well and thus the Gaussian pulse represents a suitable excitation to obtain LEMP transfer functions from time-domain analysis.
Adv   5412A, 2005), where the time from beginning to peak value 200 kA is 6.4 us, the time from beginning to decay to 50 percent peak value, that is 100 kA, is 69 us.  is contrasted to a unipolar Gaussian pulse of the form The parameters of the pulses are adjusted such that they have a similar rise time but the Gaussian pulse decays much faster if compared to the double exponential pulse.The amplitude spectrum of both pulses is shown in Fig. 3.For the double exponential pulse, the amplitude spectrum is well-known and characterized by dominant low-frequency contributions (Lee, 1995).The amplitude spectrum of the unipolar Gaussian pulse is characterized by dominant low frequency contributions as well and thus the Gaussian pulse represents a suitable excitation to obtain LEMP transfer functions from time-domain analysis.The duration of a Gaussian excitation is determined by the parameter T that is introduced in Eq. (4).It is inversely proportional to the frequency F max which is defined as the frequency where the amplitude spectrum of the Gaussian pulse has decayed by -20 dB if compared to its maximum value, Therefore a longer pulse duration implies a smaller value for F max and vice versa.This circumstance is illustrated in Fig. 4 where Gaussian pulses corresponding to F max = Gaussian Pulse Double−Exp.Pulse 3. Amplitude spectrum of both double exponential pulse and ssian pulse, as displayed in Fig. 2.
ntrasted to a unipolar Gaussian pulse of the form parameters of the pulses are adjusted such that they have ilar rise time but the Gaussian pulse decays much faster mpared to the double exponential pulse.The amplitude trum of both pulses is shown in Fig. 3.For the douxponential pulse, the amplitude spectrum is well-known characterized by dominant low-frequency contributions , 1995).The amplitude spectrum of the unipolar Gauspulse is characterized by dominant low frequency contions as well and thus the Gaussian pulse represents a ble excitation to obtain LEMP transfer functions from -domain analysis.he duration of a Gaussian excitation is determined by the meter T that is introduced in Eq. ( 4).It is inversely proional to the frequency F max which is defined as the frecy where the amplitude spectrum of the Gaussian pulse decayed by -20 dB if compared to its maximum value,

Numerical Example: LEMP Transfer Function of a Canonical Tank Model
As a specific example a canonical tank model is considered.
The tank model consists of a rectangular box which contains a fuel pipe.The fuel pipe is electrically bonded to one side of the interior of the tank and leaves the tank through a circular hole as shown in Fig. 6.The tank is subject to a lightning current source i in (t) and the induced voltage v out (t) between the fuel pipe and the boundary of the circular hole is the observable of interest, compare Fig. 7 .Then a LEMP transfer function is defined according to Eq. ( 1).
The tank model has been created within the CST Mi- is contrasted to a unipolar Gaussian pulse of the form The parameters of the pulses are adjusted such that they have a similar rise time but the Gaussian pulse decays much faster if compared to the double exponential pulse.The amplitude spectrum of both pulses is shown in Fig. 3.For the double exponential pulse, the amplitude spectrum is well-known and characterized by dominant low-frequency contributions (Lee, 1995).The amplitude spectrum of the unipolar Gaussian pulse is characterized by dominant low frequency contributions as well and thus the Gaussian pulse represents a suitable excitation to obtain LEMP transfer functions from time-domain analysis.
The duration of a Gaussian excitation is determined by the parameter T that is introduced in Eq. ( 4).It is inversely proportional to the frequency F max which is defined as the frequency where the amplitude spectrum of the Gaussian pulse has decayed by -20 dB if compared to its maximum value, Therefore a longer pulse duration implies a smaller value for F max and vice versa.This circumstance is illustrated in Fig. 4 where Gaussian pulses corresponding to F max = 1 MHz, 3 MHz, and 10 MHz are displayed.Additionally, the amplitude spectrum of the Gaussian pulse with F max = 3 MHz is shown in Fig. 5.It is suggested to use Gaussian excitations of different parameters F max , corresponding to different time durations, within the macromodelling routine of Fig. 1.In the following section this routine will be illustrated by a numerical example.

Numerical Example: LEMP Transfer Function of a Canonical Tank Model
As a specific example a canonical tank model is considered.
The tank model consists of a rectangular box which contains a fuel pipe.The fuel pipe is electrically bonded to one side of the interior of the tank and leaves the tank through a circular hole as shown in Fig. 6.The tank is subject to a lightning current source i in (t) and the induced voltage v out (t) between the fuel pipe and the boundary of the circular hole is the observable of interest, compare Fig. 7 .Then a LEMP transfer function is defined according to Eq. ( 1).
The tank model has been created within the CST Microwave Studio software which also is used to run the simulation model and to generate the time domain data v out (t) 1 MHz, 3 MHz, and 10 MHz are displayed.Additionally, the amplitude spectrum of the Gaussian pulse with F max = 3 MHz is shown in Fig. 5.
It is suggested to use Gaussian excitations of different parameters F max , corresponding to different time durations, within the macromodelling routine of Fig. 1.In the following section this routine will be illustrated by a numerical example.
www.adv-radio-sci.net/9/323/2011/Adv.Radio Sci., 9, 323-328, 2011  from the input variable i in (t) by means of the Finite Integration Technique (Clemens and Weiland, 2001).The result of a sample calculation is shown in Figs. 8 and 9.A mainly inductive response is observed, where voltage and current ap-  proximately are related by u(t) = L di(t) dt .This observation has been discussed in more detail in (Gronwald, 2010).As a general rule, LEMP transfer functions of aircraft structures are mainly characterized by inductive and resistive contributions.This is due to the quasistationary nature of the lightning current which makes radiation contributions negligible.Also capacitive effects can be neglected as long as the metallic parts of aircraft structures are properly bonded (Fisher et al., 1990).Therefore, in the present case, it is meaningful to consider the transfer admittance ) as an appropriate first order macromodel for the (inverse) transfer function of the canonical tank model.For details on macromodels of elementary networks it is referred to (Antonini, 2003).
To calculate the parameters L and R which, according to (6), determine the LEMP transfer admittance Y R (s), it is first necessary to perform a time domain calculation for a given excitation i in (t) to compute v out (t).Then the corresponding discrete transfer function H D (k) can be obtained from (1).Next, the inverse of this transfer function is fitted by means of the vector fitting method to the first order macromodel Y R (s) which is given by ( 6).This yields the pole p and residue r which both determine the parameters L and R.
Following the procedure described in the previous paragraph, macromodels for the LEMP transfer function of the canonical tank model have been obtained for nine different Gaussian excitations.The results are given in Tab. 1.The chosen pulse durations vary by a factor of 80, that is, the computation times to compute v out vary by a factor of about 80 as well.The computations were performed on a single PC   from the input variable i in (t) by means of the Finite Integration Technique (Clemens and Weiland, 2001).The result of a sample calculation is shown in Figs. 8 and 9.A mainly inductive response is observed, where voltage and current ap-  proximately are related by u(t) = L di (t)  dt .This observation has been discussed in more detail in (Gronwald, 2010).As a general rule, LEMP transfer functions of aircraft structures are mainly characterized by inductive and resistive contributions.This is due to the quasistationary nature of the lightning current which makes radiation contributions negligible.Also capacitive effects can be neglected as long as the metallic parts of aircraft structures are properly bonded (Fisher et al., 1990).Therefore, in the present case, it is meaningful to consider the transfer admittance ) as an appropriate first order macromodel for the (inverse) transfer function of the canonical tank model.For details on macromodels of elementary networks it is referred to (Antonini, 2003).
To calculate the parameters L and R which, according to (6), determine the LEMP transfer admittance Y R (s), it is first necessary to perform a time domain calculation for a given excitation i in (t) to compute v out (t).Then the corresponding discrete transfer function H D (k) can be obtained from (1).Next, the inverse of this transfer function is fitted by means of the vector fitting method to the first order macromodel Y R (s) which is given by ( 6).This yields the pole p and residue r which both determine the parameters L and R.
Following the procedure described in the previous paragraph, macromodels for the LEMP transfer function of the canonical tank model have been obtained for nine different Gaussian excitations.The results are given in Tab. 1.The chosen pulse durations vary by a factor of 80, that is, the computation times to compute v out vary by a factor of about 80 as well.The computations were performed on a single PC The tank model consists of a rectangular box which contains a fuel pipe.The fuel pipe is electrically bonded to one side of the interior of the tank and leaves the tank through a circular hole as shown in Fig. 6.The tank is subject to a lightning current source i in (t) and the induced voltage v out (t) between the fuel pipe and the boundary of the circular hole is the observable of interest, compare Fig. 7 .Then a LEMP transfer function is defined according to Eq. ( 1).
The tank model has been created within the CST Microwave Studio software which also is used to run the simulation model and to generate the time domain data v out (t) from the input variable i in (t) by means of the Finite Integration Technique (Clemens and Weiland, 2001).The result of a sample calculation is shown in Figs. 8 and 9.A mainly inductive response is observed, where voltage and current approximately are related by u(t) = L di (t)  dt .This observation has been discussed in more detail in (Gronwald, 2010).As inductive response is observed, where voltage and curren    proximately are related by u(t) = L di (t)  dt .This observatio has been discussed in more detail in (Gronwald, 2010).A a general rule, LEMP transfer functions of aircraft structure are mainly characterized by inductive and resistive contribu tions.This is due to the quasistationary nature of the light ning current which makes radiation contributions negligible Also capacitive effects can be neglected as long as the metal a general rule, LEMP transfer functions of aircraft structures are mainly characterized by inductive and resistive contributions.This is due to the quasistationary nature of the lightning current which makes radiation contributions negligible.Also capacitive effects can be neglected as long as the metallic parts of aircraft structures are properly bonded (Fisher et al., 1990).Therefore, in the present case, it is meaningful to consider the transfer admittance as an appropriate first order macromodel for the (inverse) transfer function of the canonical tank model.For details on macromodels of elementary networks it is referred to (Antonini, 2003).
To calculate the parameters L and R which, according to (6), determine the LEMP transfer admittance Y R (s), it is first necessary to perform a time domain calculation for a given excitation i in (t) to compute v out (t).Then the corresponding discrete transfer function H D (k) can be obtained from (1).Next, the inverse of this transfer function is fitted by means of the vector fitting method to the first order macromodel Y R (s) which is given by ( 6).This yields the pole p and residue r which both determine the parameters L and R.
Following the procedure described in the previous paragraph, macromodels for the LEMP transfer function of the canonical tank model have been obtained for nine different Gaussian excitations.The results are given in Tab. 1.The chosen pulse durations vary by a factor of 80, that is, the computation times to compute v out vary by a factor of about 80 as well.The computations were performed on a single PC with 3.0 GHz CPU, 4 GB RAM, using the CST Microwave Studio 2010 software on a 32-bit Windows platform.
In order to judge the resulting parameters L and R one should note that the material of the tank model is modelled as highly conducting aluminium (σ = 3.72 • 10 7 S/m).Therefore the voltage response is dominantly inductive, such that the resistive response is somewhat hidden and does not significantly contribute to the voltage response.To judge the quality of the fitted macromodels the relative root mean square error is included in Tab. 1 which quantifies the deviation between the numerically calculated transfer function H −1 D (k) and the fitted first order LEMP transfer admittance Y R (s).It is seen that this deviation is small, thus justifying the assumption of a simple first order macromodel.Towards shorter pulse durations the fit error increases.This indicates that the possibility to shorten the pulse duration is limited due to numerical inaccuracies that accompany the computation of the Fourier transformation in order to obtain H −1 D (k).However, the fit error increases towards longer pulse durations as well.To explain this observation it is noted that for long computation times numerical instabilities tend to build up, at least for the present numerical model which only contains minor losses.The smallest fit error occurs for the Gaussian pulse of time duration 1.83 us.
For further validation, it is of interest to compare the results that were obtained from the Gaussian excitations to an analogous calculation using a standard double exponential pulse as displayed in Fig. 2.However, as it turns out, the long pulse duration leads to strong oscillations that are interpreted as numerical instabilities and make it not feasible to extract a meaningful LEMP transfer function.As an alternative, an independent calculation was performed using the Method of Moment code CONCEPT-II which is defined in frequency domain (Brüns et al., 2011).To this end, a surface discretized model of the canonical tank was created and simulations at 161 discrete frequency samples between 1 kHz and 8 MHz were performed.From the resulting data a corresponding first order macromodel was obtained with parameters L = 0.5192 nH, R = 2.8887 u , and a relative root mean square error RMSE= 0.0023, thus supporting the results of Tab. 1.

Conclusions
Standard lightning pulses are characterized by durations of tenths of microseconds, thus leading to long computation times for time domain calculations of lightning current induced quantities.However, in order to obtain LEMP transfer functions it is also possible to use shorter excitation pulses, involving less computation time.Typically, LEMP transfer functions are characterized by only a small number of parameters.Therefore the vector fitting procedure is useful to obtain these parameters which, in turn, can be used to quantitatively compare transfer functions that are obtained from different excitation pulses.This makes it possible to quantify whether a short excitation pulse is sufficient to determine a LEMP transfer function with sufficient accuracy.

Fig. 2 .
Fig.2.Graphical comparison between a double exponential pulse and a Gaussian pulse.Both pulses have the same maximum value and are characterized by a similar rise rime.The double exponential pulse is a standard pulse, defined within the standard(SAE  5412A, 2005), where the time from beginning to peak value 200 kA is 6.4 us, the time from beginning to decay to 50 percent peak value, that is 100 kA, is 69 us.J. Anatzki and F. Gronwald: Accelerating the Numerical Computation of Indirect Lightning Effects by means of Vector Fitting 3

Fig. 3 .
Fig. 3. Amplitude spectrum of both double exponential pulse and Gaussian pulse, as displayed in Fig. 2.
natzki and F. Gronwald: Accelerating the Numerical Computation of Indirect Lightning Effects by means of Vector Fitting

Fig. 4 .
Fig. 4. Gaussian pulses corresponding to F max = 1 MHz, 3 MHz, and 10 MHz.A larger value of F max implies a shorter pulse duration.

Fig. 5 .
Fig. 5. Amplitude spectrum of a Gaussian pulse with F max = 3 MHz.At 3 MHz its value has decayed by 20 dB if compared to the maximum value.

Fig. 5 .
Fig. 5. Amplitude spectrum of a Gaussian pulse with F max = 3 MHz.At 3 MHz its value has decayed by 20 dB if compared to the maximum value.

Fig. 6 .
Fig. 6.Canonical tank geometry with interior tank pipe.For better visibility the front side of the tank is removed.The outer metallic frame serves as a return conductor.

Fig. 7 .
Fig. 7. Enlarged view of tank pipe and circular exit hole.The arrows between tank pipe and boundary of the circular hole indicate voltage monitors which record the voltage v out (t) during the numerical simulation.

Fig. 9 .
Fig. 9. Induced voltage v out (t) which is due to the current excitation of Fig. 8.

Fig. 6 .
Fig. 6.Canonical tank geometry with interior tank pipe.For better visibility the front side of the tank is removed.The outer metallic frame serves as a return conductor.

Fig. 6 .
Fig. 6.Canonical tank geometry with interior tank pipe.For better visibility the front side of the tank is removed.The outer metallic frame serves as a return conductor.

Fig. 7 .
Fig. 7. Enlarged view of tank pipe and circular exit hole.The arrows between tank pipe and boundary of the circular hole indicate voltage monitors which record the voltage v out (t) during the numerical simulation.

Fig. 9 .
Fig. 9. Induced voltage v out (t) which is due to the current excitation of Fig. 8.

Fig. 7 .
Fig. 7. Enlarged view of tank pipe and circular exit hole.The arrows between tank pipe and boundary of the circular hole indicate voltage monitors which record the voltage v out (t) during the numerical simulation.
ronwald: Accelerating the Numerical Computation of Indirect Lightning Effects by means of Vector Fitting nk geometry with interior tank pipe.For better ide of the tank is removed.The outer metallic turn conductor.ew of tank pipe and circular exit hole.The ar-

Fig. 9 .
Fig. 9. Induced voltage v out (t) which is due to the current excita tion of Fig. 8.

Fig. 9 .
Fig. 9. Induced voltage v out (t) which is due to the current excitation of Fig. 8.

Table 1 .
Resulting macromodels for various Gaussian excitation pulses of different durations.