Numerical computation of lightning transfer functions for layered, anisotropically conducting shielding structures by the method of moments

A formalism for the computation of lightning transfer functions by the method of moments, which involves shielding structures that may consist of layered, anisotropically conducting composite materials, is presented in this contribution. The composite materials, being of a type that is widely used in spaceand aircraft design, are electrically characterized by an equivalent conductivity. As basis for the quantitative analysis the method of moments is used where shielding surfaces can be treated by a thin layer technique which utilizes analytical solutions inside the layer. Also the effect of an extended lightning channel can be taken into account. The method is applied to geometries that resemble an actual airplane fuselage.


Introduction
Despite the constantly ongoing progress in the development and application of numerical methods in electromagnetics it turns out that the computation of lightning-related effects in the framework of Electromagnetic Compatibility (EMC) still constitutes a highly challenging task.This is due to a number of difficulties that can be characterized as follows (Apra et al., 2008;Anatzki and Gronwald, 2012;Prost et al., 2013): First, the modeling of an actual lightning channel as electromagnetic source requires to turn a physically complicated and geometrically extended excitation into a numerical model.Second, the rather long duration of a lightning electromagnetic pulse (LEMP), together with its associated low-frequency spectrum, requires both in time and frequency domain stable and efficient numerical algorithms.Third, in actual applications it is often necessary to calculate LEMP transfer functions of complex systems, such as aircraft, for example.This, in turn, often involves to model advanced materials and to deal with a high degree of complexity.
In this contribution it is outlined how to numerically calculate LEMP transfer functions by the method of moments (MoM), taking into account the main difficulties mentioned above.To this end, in Sect. 2 it is outlined, based on early work (Brüns, 1985), how to take into account the influence of the lightning channel in a realistic way.The modeling of composite materials as an example for an advanced material is explained in Sect.3. It is also mentioned that lowfrequency shielding can be modeled by a thin layer technique (Bürger et al., 1995).Then, in Sect.4, it is summarized how to process frequency domain results obtained by the MoM in order to obtain time domain results that are useful for lightning analysis.Simulation examples that implement these methods are given in Sect.5, followed by a short summary in Sect.6.

Numerical modeling of the lightning channel
There are several engineering approaches for the modeling of the electric current of a lightning channel (Rakov and Uman, 2003).They relate the current at the channel base to the current at any position on the channel, where the channel itself can be modeled as a wire of length of several kilometers.Most models are categorized to be either of the Transmission Line (TL) or of the Traveling Current Source (TCS) type.The TL type is used here, therefore it shortly is described in the following subsection.More details on its implementation within the MoM can be found in (Brüns, 1985).
Published by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V.

The transmission line model
In the TL type model (Uman and McLain, 1969) a current wave is traveling upwards with propagation speed v f , starting at the channel base, as depicted in Fig. 1.At positions above z = v f • t the current is zero.The current at any position z along the channel and at time t is given by where i 0 (t) is the current pulse at the channel base and v f is the propagation speed, which typically is chosen to be one third of the vacuum speed of light v f = 1 3 • c.

The current pulse
For the current i 0 (t) the current pulse according to (Heidler, 1985), which can be found in the VG 95371-10 standard, is used.The slope of this curve is continuous, which has proven to be consistent with the results of lightning measurements.
The entire pulse can be divided into three components that are multiplied by each other and determine the amplitude, the smooth rising edge, and the exponential decay, respectively.The corresponding equation is given by where i max is the maximum value of the current, η is the amplitude correction factor, which ensures that the pulse reaches the value i max , T rise is the rise time, and T hold the decay time of the pulse.For the following investigations the parameters i max = 100 kA, η = 0.986, T rise = 1.82 µs, and T hold = 285 µs are used.This choice refers to a negative first stroke of threat level "severe".A pulse with these parameter values is shown in Fig. 2.

Modeling of composite material
Composite materials are widely used in aircraft design for weight reduction and improvement of mechanical strength.In this publication a type of carbon fiber composite is considered which consists of several layers of carbon fibers with different orientations, enclosed in resin.

Equivalent conductivity
Instead of modeling single fibers one may assume an anisotropic conductivity for each layer (Holloway et al., 2005;Happ et al., 2013), where the direction of the fibers in the layer corresponds to the direction of the highest conductivity in the conductivity tensor.
A further simplification can be made by replacing the various anisotropically conducting layers by a single layer with the same overall thickness and an equivalent isotropic conductivity σ eq .An analytical investigation of the shielding effectiveness of a plane shield with several anisotropically conducting layers has shown that this simplification is valid up to several tens of MHz (Happ et al., 2015).Since lightning is a low frequency phenomenon, the equivalent conductivity can be used to model the shielding properties of the composite material with respect to lightning effects.As a first step to calculate σ eq the diagonal entries of the average conductivity tensor are calculated by where L is the number of layers.Here it is assumed that each layer of fibers has the conductivity σ ξ in fiber direction and σ ψ in cross-fiber direction.The orientation of the fibers is specified by the angle α i .The thickness of layer i is denoted by d i and the total thickness of the multilayer material is denoted by d.Then the equivalent conductivity is chosen to be the smaller of the two conductivity values, to give a valid value for a worst case approximation.In the following, the conductivities are chosen to be σ ξ = 10 kS m −1 and σ ξ = 100 S m −1 and each layer has thickness For a symmetric layer pattern, e.g., a multiple of four layers with a relative rotation angle of 45 • between the fibers of adjacent layers, the equivalent isotropic conductivity is given by σ eq = 5050 S m −1 .

Thin layer technique
There are several methods to numerically model thin layers.
The surface impedance boundary condition method (Harrington and Mautz, 1975;Chiang and Chew, 2006), for example, is applicable to scattering problems, but it is not suitable for the calculation of shielding effectiveness, since the coupling between inside and outside region of the shielding geometry has to accurately be taken into account.It might also come to mind to apply a Green's function of layered media (Michalski and Mosig, 1997;Ginste et al., 2010), but this typically requires some kind of canonical symmetry of the geometry considered.Therefore it is not immediate to apply this method to arbitrarily shaped three dimensional structures.The thin layer technique that is applied here to efficiently model thin layers of finite size in conjunction with the MoM (Bürger et al., 1995) is based on an analytically formulated coupling matrix.In this case, the layer has to fulfill the requirement to be thin compared to the overall dimension of the body to be modeled and it has to provide a sufficiently high conductivity.If the conductivity is large enough then the wave propagation inside the layer is perpendicular to its surface and can be described by an analytical solution in the form of a coupling matrix which also correctly incorporates all effects related to the wave propagation in a lossy medium.As a consequence, two regions that are separated by a two-dimensional layer, compare Fig. 4, can be treated by the MoM, where the coupling through the layer is taken into account by a coupling matrix which relates the tangential fields at both sides of the layer to each other.This hybrid-technique, which combines the MoM with an analytical solution, has proven to be stable down to frequencies in the kHz range (Happ et al., 2014).

Time domain simulations with a frequency domain solver
A priori, the lightning current is formulated in time domain while the MoM and the layer technique are formulated in frequency domain.In this section it is explained how to relate both formulations to calculate the impulse response of a lightning current.
The spectrum I 0 (ω) of the current i 0 (t) at the channel base can be determined analytically, as described in (Andreotti et al., 2005).The spatial distribution of the current at any position z of the channel is a time shifted version of i 0 (t), therefore the shift property of the Fourier transformation can be used to write the corresponding spectrum as which clearly states that the magnitude of the spectrum does not depend on the position and only the phase is influenced by z .Now the system response G(ω) to a unit excitation, in this case an impressed current with a constant amplitude with respect to frequency, which flows spatially distributed on the channel, can be calculated.Then both spectra can be multiplied in frequency domain to find the spectrum of the pulse response R(ω), which could, e.g., be a voltage or a field value.The function R(ω) can be transformed to time domain via an inverse Fourier transformation to yield the pulse response r(t).
The calculated spectrum consists of values that refer to discrete frequencies.This leads to a periodic extension of the pulse response in time domain.The time t ex , when the periodic extension starts, is related to the frequency difference f between the discrete frequencies of the spectrum by To ensure meaningful simulation results, f has to be small enough such that the impulse response has decayed before the periodic extension starts.This implies that a long excitation pulse requires a small frequency step width, which sets up a limit for this method, due to the fact that the MoM becomes unstable at very low frequencies.Hence, excitation signals with a duration of several milliseconds cannot accurately be modeled by this method.For the pulse in Fig. 2 a frequency step width of f = 600 Hz, which corresponds to the time t ex = 1.666 ms, is used for the simulations in the next section.Another important point is the maximum frequency that has to be considered to reproduce the steep rise of the pulse.In case of a lightning pulse with a rise time of a few microseconds a maximum frequency in the range of a few MHz is sufficient.
Finally, it should be mentioned that this is a linear formulation, hence non-linear effects, which could result from matter interaction at very high field magnitudes, are not taken into account.

Simulation examples
In this section the proposed formalism is illustrated by means of several examples.In all cases the pulse introduced in Sect.2.2 is used as excitation.In the case of a nearby lightning the channel starts at (x, y, z) = (0, 0, 0) on a perfectly conducting ground and ends at (0, 0, l), where in our case the channel length l is chosen to be of length 3 km.As depicted in Fig. 1, the current pulse starts at the channel base an moves upwards with speed v f .
The first configuration to be considered is a single lightning channel without a neighboring structure.Then a transmission line structure, which is loaded by 50 resistors, is placed close to the lightning channel.In a next step, this transmission line structure is located inside a cylindrical cavity.The geometrical details of the resulting setup are shown in Fig. 5.The cylinder is chosen either as a closed one with conductivity σ = 5 kS m −1 , which is a realistic equivalent conductivity for a carbon fiber composite, or a perfectly electrically conducting (PEC) cylinder with 12 apertures on both sides, representing a fuselage of a passenger aircraft.As excitation a lightning channel close to the structure or a lightning channel directly attached to the structure is assumed.The resulting four different cases are summarized in Fig. 6, where only short sections of the lightning channels are shown.In the case of a direct strike the impressed current flows from the top side of the structure upwards along the lightning channel.To model the discharge along the fuselage and towards the ground a second wire connects the bottom side of the cylinder to the PEC ground.The current on this second wire is not an impressed one, as the one on the long wire that models the lightning channel, it rather results from the MoM simulation.

Field of the lightning channel
As a prerequisite, in this subsection the electric and magnetic fields of the lightning channel without a neighboring structure are investigated and compared to the semi-analytical formulas derived in (Uman et al., 1975).With these formulas the fields at any position on the ground at z = 0 can be calculated as a function of the horizontal distance to the channel. .Electric field of the channel at different distances from the channel at z = 0, calculated by the numerical method presented in this paper and a semi-analytical formula taken from (Uman et al., 1975). .Magnetic field of the channel at different distances from the channel at z = 0, calculated by the numerical method presented in this paper and a semi-analytical formula taken from (Uman et al., 1975).
The corresponding curves for the electric and magnetic fields are shown in Figs.7 and 8, respectively.The numerical results obtained from the MoM formalism are in excellent agreement with the semi-analytical results, as exemplified by three observation points at distances 25, 50 and 100 m.For these distances the effect of the finite length of the lightning channel in the used model is negligible.
In Fig. 9 the same curves for the electric field as in Fig.

Structure close to the lightning channel
The electric field at the center of the transmission line structure, which corresponds to the point (25, 0, 100) m, with and without the finitely conducting cylinder as shield is illustrated in Fig. 10.In both cases the y component of the field is zero.As expected, the amplitudes of the electric field in the presence of the shield are much smaller if compared to the situation without shield.
The magnetic fields for these two cases are plotted in Fig. 11.The rise time of the magnetic field inside the finitely conducting cylinder is much larger if compared to the rise time of the magnetic field without the cylinder.This is due to the low pass behavior of the shield with respect to the magnetic field, i.e., the high frequency components of the field are attenuated and therefore the rise time is increased.
The related transfer function G(ω), as defined by the ratio between the voltage U (ω) at the termination resistance and an impressed current of constant amplitude of 1A per frequency, is shown in Fig. 12.This transfer function includes the following effects: the creation of the fields of the lightning channel, the coupling into the cylindrical cavity by diffusion or aperture coupling, field distribution inside the cavity, and coupling into the transmission line structure.At low frequencies the curves for the case of the transmission line alone and the transmission line inside the finitely conducting cylinder are on top of each other.This shows that in this frequency region the presence of the conducting cylinder is hardly relevant.The magnitude of the transfer function for the case of the PEC cylinder with apertures is lower compared to the other two cases because the magnetic field may only couple through the apertures into the fuselage.This effect is small for the considered wavelengths which are much larger than the aperture dimensions.Finally, the time-domain voltage at the terminating resistor is shown in Fig. 13 as a response to the current pulse according to Fig. 2. A delay time of 1 µs can be observed, which corresponds to the time the pulse needs to travel from the ground to the height of 100 m where the structure is located.The maximum voltage for the cases of no cylinder, finitely con- ducting cylinder, and PEC cylinder with apertures are 17 kV, 800 V, and 1.5 V, respectively.The main contribution to the voltage is due to the magnetic field.Due to the low pass characteristic of the finitely conducting cylinder the voltage has a decreased rise time, similar to the magnetic field itself.

Direct lightning strike
In this subsection cases are considered where the lightning channel is directly attached to the structure, compare Fig. 6b  and d, such that the lightning current flows directly on the surface of the cylindrical structure.
In Fig. 14 the spectrum of the transfer functions for both the finitely conducting cylinder and the PEC cylinder with apertures is shown.Two maximum values can be identified in the considered frequency range.The first maximum occurs at a frequency of 372 kHz.This frequency is too low for being a resonance frequency of the cylinder.Analysis shows that this frequency is related to the wire that connects the cylinder to the ground.More precisely, the frequency turns out to be the antiresonance frequency of the wire (Schelkunoff and Friis, 1952), i.e., the imaginary part of the impedance of the wire as seen from the current source is zero and the real part of the impedance exhibits a maximum.The second maximum of the transfer function occurs at the frequency 1.6 MHz, where the length of the wire from the cylinder to the ground equals one half of a wavelength.
The time domain response of the voltages are shown in Fig. 15.Oscillations are clearly visible, where the time period of the oscillation corresponds to the frequency where the first maximum of the transfer function occurs.For the case of the finitely conducting cylinder the voltage is considerably higher if compared to the case of the PEC cylinder.Therefore, in this example, the diffusion coupling is larger if compared to the aperture coupling, at least in the considered frequency range which is relevant for lightning analysis.

Conclusions
A formalism for the calculation of LEMP transfer functions by means of the MoM has been proposed.Composite materials are modeled by an equivalent conductivity which is applicable for the frequency range of the lightning spectrum.
In the examples considered it turned out that diffusion coupling has a larger influence on the transfer function than the aperture coupling.Therefore, the conductivity and thickness of the shell of the fuselage are very important parameters and it clearly is not sufficient to approximate composite materials by PEC material for lightning analysis.

FFigure 1 .
Figure 1.Illustration of the upward moving current front in the lightning channel.

Figure 3 .
Figure 3. Local coordinate system for one layer of carbon fibers, where the ξ -direction corresponds to the orientation of the fibers.

Figure 4 .
Figure 4. Two MoM regions that are separated by a thin finitely conducting layer.Its electromagnetic properties can be described by an analytical formulation.

Figure 5 .Figure 6 .
Figure 5. Dimensions and positions of the transmission line structure and its surrounding cylinder, representing a simple model of an aircraft fuselage.
Figure7.Electric field of the channel at different distances from the channel at z = 0, calculated by the numerical method presented in this paper and a semi-analytical formula taken from(Uman et al., 1975).
Figure8.Magnetic field of the channel at different distances from the channel at z = 0, calculated by the numerical method presented in this paper and a semi-analytical formula taken from(Uman et al., 1975).

Figure 9 .Figure 10 .Figure 11 .Figure 12 .
Figure 9. Electric field of channel at different distances from the channel at z = 0, calculated for a longer time interval.

Figure 13 .
Figure 13.Voltage at the termination resistor of the transmission line with and without conducting cylinder (a) and with PEC cylinder with apertures (b) for the structures being close to the lightning channel.

Figure 14 .Figure 15 .
Figure14.Frequency spectrum of the system response for the cases of a closed conductive cylinder and a PEC cylinder with apertures if a direct lightning strike is applied.