Samuel Gaucher

Member, IEEE Christophe Guiffaut

Member, IEEE Alain Reineix

Olivier Cessenat

Geneviève Mazé-Merceur

Wideband simulations of periodic structures by the Hybrid Spectral FDTD/TD-VFz Method

Keywords: periodic boundary condition, finite-difference time-domain, oblique incidence, Spectral FDTD, time-domain vector fitting

come

I. INTRODUCTION

M ETAMATERIALS are used in many applications such as antennas, cloaking devices or radar cross section (RCS) reduction. Some of them are composed of periodic patterns. Thus, periodic Floquet boundary conditions (PBC) are applied in periodic directions such that only one unit cell needs to be analyzed instead of the entire structure. Particularly, the finite-difference time-domain (FDTD) is a good candidate to compute the unit cell because of its simplicity, efficiency and wideband capabilities. However when the incident angle is fixed, time-advanced solutions are required on boundaries. Various techniques exist to deal with this difficulty and [START_REF] Taflove | Computational electrodynamics: the finitedifference time-domain method[END_REF] distinguishes them into two classes: 'direct field methods' and 'field-transformation methods'.

Spectral FDTD (SFDTD) [START_REF] Yang | A simple and efficient FDTD/PBC algorithm for scattering analysis of periodic structures[END_REF] [3] is a direct field method which keeps the usual electric field E and magnetic field H unknowns. Then, all existing FDTD models and processings can straightforwardly be adapted to the SFDTD. In this approach, the horizontal wavenumber is fixed rather than the incident angle. Then the implementation of PBC becomes simple because no time-advanced data are required. The best advantage of this technique is the ability to use the standard Yee scheme without restriction on the Courant-Friedrich-Levy (CFL) condition, which is not angle-dependent [START_REF] Yang | A simple and efficient FDTD/PBC algorithm for scattering analysis of periodic structures[END_REF]. On the other hand, when the horizontal wavenumber is fixed, the SFDTD scheme simultaneously excites propagating wave and evanescent wave regions [START_REF] Yang | A simple and efficient FDTD/PBC algorithm for scattering analysis of periodic structures[END_REF] which are separated by the cut-off frequency. As a result, a resonant behavior occurs and yields undamped oscillations in the time signal. Thus, its discrete S. Gaucher and O. Cessenat are with Commissariat à l'Energie Atomique et aux Energies Alternatives (CEA-CESTA), Le Barp 33116, France, (e-mail: samuel.gaucher@cea.fr; olivier.cessenat@cea.fr).

S. Gaucher, C. Guiffaut and A. Reineix are with CNRS Limoges, XLIM, France, Limoges 87000, (email: samuel.gaucher@xlim.fr, christophe.guiffaut@xlim.fr, alain.reineix@xlim.fr).

This work is funded by the DGA AID and the CEA CESTA. Manuscript received month day, year; revised month day, year.

Fourier transform (DFT) provides an inaccurate frequency response. To remedy this, the SFDTD method with the ARMA post-processing has been developed in [START_REF] Yang | A hybrid spectral-FDTD/ARMA method for periodic structure analysis[END_REF]. Starting from the same idea, we propose to combine the SFDTD scheme with the discrete-time time-domain vector fitting (TD-VFz) [START_REF] Lei | Efficient linear macromodeling via discretetime time-domain vector fitting[END_REF] extrapolation technique. This approach is very robust, guaranteed stable poles and is a much lower model order than ARMA [START_REF] Ubolli | Comparison of methods for rational approximation of simulated time-domain responses: ARMA, ZD-VF, and TD-VF[END_REF]. The organization of this paper is as follows. The SFDTD method and the resonance problem are explained in Section II. In Section III, the resonance problem is solved and the eigenfrequency is identified when consider a homogeneous layer. A comparative study of the time signal extrapolation techniques TD-VFz, ARMA and matrix-pencil (MP) [START_REF] Sarkar | Using the matrix pencil method to estimate the parameters of a sum of complex exponentials[END_REF] was performed to deduce the good efficiency of the TD-VFz method. In Section IV, the hybrid SFDTD/TD-VFz scheme is presented. First, a simple methodology to construct the spectral response for an arbitrary incident angle is proposed: The mapping between the wavenumber to angle response is directly obtained by the TD-VFz rational transfer function. Secondly, a convolution perfectly matched layers (CPML) setting is given to avoid late-time divergence and to improve the accuracy near the cut-off frequency. In section V, two numerical examples are implemented to demonstrate the hybrid method performance. The hybrid method is compared to the Material Independent (MI) scheme [START_REF] Liang | Wideband analysis of periodic structures at oblique incidence by material independent FDTD algorithm[END_REF] for verification. MI and SFDTD simulations are performed with the XLIM laboratory homemade Temsifd solver [START_REF]Time ElectroMagnetic Simulator-Finite Difference Software[END_REF] and SFDTD temporal data are extrapolated with TD-VFz implemented in Matlab scripting.

II. SPECTRAL FDTD REVIEW

Let us assume a (l x , l y ) periodicity in the (x, y) directions with a z-propagation of an incident wave. Like [START_REF] Zhou | Efficient simulations of periodic structures with oblique incidence using direct spectral FDTD method[END_REF], an arbitrary azimuth angle is chosen to ϕ = ϕ 0 . Then, the horizontal wavenumber depends on the angular frequency ω = 2πf and the elevation angle θ as k h = ω sin θ/c 0 . c 0 is the speed of light in free space. As ϕ, k h is also fixed for one simulation of the SFDTD method. The constant components k x and k y of the plane wave vector k are deduced as a function of the azimuth angle and the horizontal wavenumber

k x = k h cos ϕ 0 , k y = k h sin ϕ 0 . (1)
By fixing them, simple time-domain PBC are found for

Ψ = E or H Ψ(x, y, z, t) = Ψ (x + ml x , y + nl y , z, t) e j(kxmlx+kynly) , (2)
where m and n are any integers. Then, SFDTD scheme is exactly the same as the FDTD one with complex E-H fields unknowns instead of real for the latter. We deduce from the expression θ(ω) = arcsin k h k 0 that different frequencies correspond to different incident angles, and for a positive real incidence, the condition k 0 ≥ k h must be verified. In other words, this is referred to the propagating wave region if f ≥ f c else to the evanescent wave region, where

f c = c 0 k h 2π (3)
is the cut-off frequency. Because both the propagating wave region and the evanescent wave region are excited by the SFDTD method, surface waves are guided along horizontal directions at some frequencies called 'eigen-frequencies'. Thus, the energy of a horizontal wave continuously re-enters the computational domain through the PBC instead of being absorbed by the top and bottom CPML [START_REF] Yang | A simple and efficient FDTD/PBC algorithm for scattering analysis of periodic structures[END_REF]. As a result, a resonance behavior occurs: temporal data do not converge and an accurate spectral response cannot be obtained by DFT.

To overcome this difficulty, [START_REF] Yang | A simple and efficient FDTD/PBC algorithm for scattering analysis of periodic structures[END_REF] adopts a modulated Gaussian waveform to excite only the propagating wave region

F inc (t) = exp - (t -t 0 ) 2 T 2 exp (j2πf 0 t), (4)
where t 0 is the time delay, T the pulse width and the signal magnitude at f c should be at least 40 dB lower than the signal level at the center frequency f 0 . This simple approach improves the spectral response accuracy by attenuating the horizontal resonance. However, especially for very wideband excitation, there are still some artifacts near the cutoff frequency corresponding to high incident angles in the spectral domain. Moreover, only the propagating wave region is determined with this technique and not the surface wave region.

III. EIGEN-FREQUENCIES IDENTIFICATION WITH TD-VFZ

Let us consider a homogeneous layer of thickness 9.375 mm and dielectric constant r = 4. FDTD simulation runs with k h = 50 rad/m, ϕ = 90°and a Gaussian pulse excitation. The time domain reflected field of the fundamental mode on the interface between the air and the layer is shown in Fig. 1. Temporal data do not converge and oscillations are observed.

To obtain a correct spectral reflection coefficient R, the ratio between the output (real part of the reflected field) and input (excitation) response is decomposed with TD-VFz [START_REF] Lei | Efficient linear macromodeling via discretetime time-domain vector fitting[END_REF], as a rational transfer function

R(z = exp [jω∆ t ]) = d + M n=1 r n 1 -z -1 s n , ( 5
)
where d is a constant, r n and s n are the complex residues and poles respectively and ∆ t is the time sample. Fig. 2 illustrates the reflection coefficient over a wideband from 0 + GHz to 25 GHz obtained by the SFDTD scheme with the modulated Gaussian excitation (4) and a strong attenuation Att f max = 60 dB at f c , and by the SFDTD scheme combined with the extrapolation technique ARMA, MP or TD-VFz for the same M = 18 poles number (p = q = 18 for , TE mode Analytic SFDTD, Att_fmax = 60 dB SFDTD+MP, M " 18 SFDTD+ARMA, p " q " 18 SFDTD+ARMA, p " q " 114 SFDTD+TDVFz, M " 18 Fig. 2: Reflection coefficient for the homogeneous layer.

ARMA [START_REF] Yang | A hybrid spectral-FDTD/ARMA method for periodic structure analysis[END_REF]) and 5 iterations for TD-VFz poles convergence.

The FDTD spatial steps are all equal to ∆ = 0.1875 mm and CFL = 0.99. Excellent agreement is observed between the hybrid SFDTD/TD-VFz method and analytical reference, especially near the cut-off frequency, on the contrary to the alone Spectral FDTD scheme where the response is noisy.

The MP technique gives a good result but with an accuracy loss for high frequencies. In addition, the number of poles is insufficient for the ARMA method which does not give the expected result. An interesting point consists in the eigenfrequencies identification. When [START_REF] Lei | Efficient linear macromodeling via discretetime time-domain vector fitting[END_REF] gives poles on the unit circle boundary, resonances occur and eigen-frequencies are identified. For the sample k h = 50, there is only one reported in the Table I for several extrapolation techniques and poles number choices. For the order M = 18, the best matching with theory is obtained with the SFDTD/TD-VFz scheme. We note that p = q = 114 for ARMA are required to obtain the same result as TD-VFz with M = 18. Furthermore, M = 24 gives the same excellent result for MP and TD-VFz but M < 24 for MP gives an inaccurate response at high frequencies. Regarding the CPU time, ARMA (p = q = 114) has a slightly faster calculation time than TD-VFz (M = 18) and MP is a slower method. Also the memory usage is 30 MB for ARMA and TD-VFz, and 100 MB for MP (M = 18). IV. HYBRID SPECTRAL-FDTD/TD-VFZ

A. Mapping R(0 < k h , f ) to R(0 < θ, f )
The normal incidence case θ = 0°is trivial and corresponds to k h = 0. To deal with the oblique incident angles for the frequency band [0 < f min , f max ] and incident angle range [0 < θ min , θ max ] of interest, we first calculate the extremum horizontal wavenumber as

k hmax = 2πf max c 0 sin θ max , k hmin = 2πf min c 0 sin θ min . (6)
To compute the reflection coefficient as a function of the incident angle, we run the FDTD simulation for the N samples of k i h with the fixed azimuth angle ϕ = ϕ 0

k i h = k hmin + i k hmax -k hmin N -1 , i = 0, . . . , N -1, (7)
by using (1) to determine the phase values in the PBC formulation [START_REF] Yang | A simple and efficient FDTD/PBC algorithm for scattering analysis of periodic structures[END_REF]. Then, each response with sample k i h is extrapolated by the M-order TD-VFz algorithm to obtain the M residues r i n and poles s i n and the constant d i . After that, an incident angle-reconstruction is applied. For a given angle θ min ≤ θ j ≤ θ max , we compute the N frequencies

f i = ck i h 2π sin θ j , i = 0, . . . , N -1. (8)
The reflection coefficient R evaluated at angle θ j is calculated with the transfer function [START_REF] Lei | Efficient linear macromodeling via discretetime time-domain vector fitting[END_REF] for all frequencies

f i R(f i , θ j ) = d i + M n=1 r i n 1 -exp [-j2πf i ∆ t ]s i n . ( 9
)
As frequencies are imposed by the formula (8), the Lagrange interpolation is used to obtain the reflection coefficient at the desired frequencies. Assuming a desired frequency f i ≤ f ≤ f i+1 , the reflection coefficient for the angle θ j and frequency f is given by

R( f , θ j ) = f i+1 - f f i+1 -f i R(f i , θ j ) + f -f i f i+1 -f i R(f i+1 , θ j ). (10)

B. CPML setting

The stretched-coordinate metrics s z of the complex frequency shifted CPML (CFS-CPML) approach [START_REF] Roden | Convolution PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media[END_REF] is defined in the z-direction as

s z = κ z + σ z α z + jω 0 , ( 11
)
where σ z is the conductivity depending of the z-depth layer of the CPML to provide absorption for propagating waves. Parameters κ z and α z better attenuate the evanescent waves at low frequencies, by stretching the FDTD spatial step with s z as can be seen by taking the limit ω → 0. In this paper, the parameter α z is constant with the PML cut-off frequency defined by

f P M L c = α z 2π 0 . ( 12
)
The condition f P M L c < f c is a good compromise for absorbing propagative waves without neglecting the absorption of evanescent waves. Then, we propose

f P M L c = νf c , (13)
where ν should be a real between 0.15 and 0.75 according to our numerical experiences. The constant α z profile in CPML is found by using ( 3)

α z = ν k h Z 0 , (14)
where Z 0 ≈ 376.730 Ω is the impedance in free space. It is important to have α z > 0 to avoid late-time divergence. The real part of the time domain reflected field for the same case as section III and during a long time simulation of 20 ns is shown in Fig. 3 for two α z values. Time domain signal diverges when α z = 0 but the setting ( 14) with ν = 0.75 efficiently corrects the problem. On the other hand, the conductivity grows geometrically in CPML because it is the best profile for lowering the reflection of the evanescent waves [START_REF] Bérenger | An optimized CFS-PML for wave-structure interaction problems[END_REF]. Moreover, it is very interesting to enhance the absorption of the waves at grazing incidence since they are close to the cut-off frequency and therefore disturbed by the resonance frequencies. Consequently for a CPML of 12 cells in thickness, we adopt the optimum geometrical expansion g = 1.9 with the theoretical normal incidence reflection coefficient expected at the interface between air and CPML equals to R(0) = 10 -14 . Theses setting are obtained from the PML theoretical reflection of [START_REF] Bérenger | Evanescent waves in PML's: Origin of the numerical reflection in wave-structure interaction problems[END_REF]. Therefore PML outperforms for grazing incidence wave absorption. The counterpart can be a lesser absorption for normal incidence or weak incidence angle waves but the reflection coefficient remains below -50 dB. Besides this behaviour is stable with the frequency.

V. NUMERICAL RESULTS

For the two numerical examples of this section, periodic solvers run with a 2 ns simulation time and TD-VFz with M = 40 poles and 15 iterations for the poles convergence. The setting [START_REF] Demarest | FDTD modeling of scatterers in stratified media[END_REF] with ν = 0.35 is chosen for all simulations.

A. Pattern composed of two homogeneous lossy layers

Let us consider a stratified media composed of two arbitrary homogeneous layers. The geometry, dimensions and parameters of the layers are plotted in Fig. 4 (a). The model is simulated by using the proposed method for N = 100 horizontal wavenumber samples with f max = 10 GHz, θ max = 80°, ∆ = 0.2 mm and CFL = 0.99. The numerical reflection coefficient and phase magnitude for the TE mode with a constant azimuth angle ϕ = 90°is compared to the theory [START_REF] Demarest | FDTD modeling of scatterers in stratified media[END_REF] as shown in Fig. 5 and Fig. 6 respectively. A very good agreement with theoretical data is observed for many incident angles. Note that the accuracy is less good at high incidence angles 60°and 80°for low frequencies (8) below 0.5 GHz approximatively due to the strong proximity of them with the cut-off frequency.

B. Jerusalem Cross FSS with anisotropic and lossy media

The Jerusalem cross [START_REF] Liang | Wideband analysis of periodic structures at oblique incidence by material independent FDTD algorithm[END_REF] is a perfect electromagnetic conductor (PEC) between two anisotropic and lossy FSS media with relative permittivity and conductivity diagonal tensors given by r = diag (2.2, 1.1, 1.5) and σ = diag (0.05, 0.03, 0.01). Pattern geometry is plotted in Fig. 4 (b). The model is simulated by using the proposed method for N = 100 horizontal wavenumber samples over a wideband from 5 GHz to 25 GHz with θ max = 80°, ∆ = 0.125 mm and CFL = 0.99. Considering the TM mode and a constant azimuth angle ϕ = 90°, very good agreement between the hybrid and MI scheme is observed in Fig. 7 for various incident angles. Note that the CFL condition of the MI field transformation method is c 0 √ 3∆ t ≤ ∆(1-sin θ). Thus, the simulation time becomes prohibitively long at high incidence. The best advantage of the hybrid method is that the CFL condition is the same for each horizontal wavenumver sample k h as c 0 √ 3∆ t ≤ ∆ [START_REF] Yang | A simple and efficient FDTD/PBC algorithm for scattering analysis of periodic structures[END_REF].

VI. CONCLUSION

The Spectral FDTD scheme combined with the TD-VFz decomposition is an efficient hybrid approach to analyze periodic structures excited by oblique incident plane waves without reducing the CFL stability criterion of the FDTD method. The VF-poles perfectly fit the response resonant modes allowing the hybrid method to deal with high incident angles. CFS-PML is successfully used to avoid late time instabilities with a simple setting of the low-frequency stretching coefficient α from the cut-off frequency of the PBC. In addition, the PML conductivity profile is optimized for grazing incident angle increasing the accuracy computation of the reflection coefficient.

Fig. 1 :

 1 Fig. 1: Reflected wave in the time domain. Resonance behavior is observed.

0 Fig. 3 :

 03 Fig. 3: Time domain signal for two values of α (k h = 50).

Fig. 4 :

 4 Fig. 4: Unit cell geometry. (a) two homogeneous layers, (yOz) plane. (b) JCFSS, left: (xOy) plane, right: (yOz) plane.

Fig. 5 :

 5 Fig.5: Reflection coefficient magnitude for the two layers.

Fig. 6 :

 6 Fig. 6: Reflection coefficient phase for the two layers.

Fig. 7 :

 7 Fig. 7: Reflection coefficient magnitude for the JCFSS.

TABLE I :

 I Eigen-frequency prediction. The theory gives f r = 2.0751 GHz.

		Prediction (GHz)	Error (%)	CPU time (s)
	ARMA (p=q=18)	2.2614	8.978	0.06
	MP (M=18)	2.0758	0.034	3.3
	TD-VFz (M=18)	2.0746	0.024	1.38
	MP (M=24)	2.0752	0.005	3.11
	TD-VFz (M=24)	2.0752	0.005	1.64
	ARMA (p=q=114)	2.0746	0.024	1.2


Samuel Gaucher
email: samuel.gaucher@xlim.fr

Alain Reineix
email: alain.reineix@xlim.fr

Olivier Cessenat
email: olivier.cessenat@cea.fr

Hybrid Spectral FDTD/TD-VFz method for analysis of periodic absorbers 1 st

Keywords: periodic bounday condition, oblique incidence, Spectral FDTD, time-domain vector fitting, ultra-broadband absorbers, polarization-insensitive

The Spectral FDTD scheme is combined with the discrete-time time-domain vector fitting decomposition technique to analyze wide-angle ultra-broadband absorbers with polarization-insensitive characteristics.

I. INTRODUCTION

Periodic absorbers are often required for several applications. Their characterization versus frequency often requires numerical modeling. Thus, periodic boundary conditions (PBC) are applied to simulate only one pattern rather than the entire structure. For this purpose, the efficiency, wideband and simple finite-difference time-domain (FDTD) method is used to compute the electromagnetic (EM) fields. However, PBCs for oblique plane wave illumination cannot be performed with the standard FDTD scheme. A modified FDTD scheme such as Spectral FDTD (SFDTD) [START_REF] Yang | A simple and efficient fdtd/pbc algorithm for scattering analysis of periodic structures[END_REF] [2] method is a good candidate because few changes are required to adapt a FDTD solver and its numerical models. Moreover no additional constraint is needed on the CFL criterion. However the horizontal wavenumber is fixed with this approach and not the incident angle. As a result, unwanted evanescent and horizontal modes [START_REF] Yang | A simple and efficient fdtd/pbc algorithm for scattering analysis of periodic structures[END_REF] are excited leading to a low time convergence and then inaccurate spectral responses. To overcome this problem, we propose to combine the SFDTD method with the discrete-time time-domain vector fitting (TD-VFz) [START_REF] Lei | Efficient linear macromodeling via discretetime time-domain vector fitting[END_REF] to identify the eigen-frequencies of the periodic structure. In addition to efficiency solve the resonance problem, the frequency versus wavenumber mapping is straightforwardly transformed to frequency versus incident angle one by using the TD-Vfz rational transfer function. To validate the proposed method, angular-stable ultra-broadband absorbers from literature with polarization-insensitive characteristics are computed with the XLIM laboratory homemade Temsi-fd solver [START_REF]Time ElectroMagnetic Simulator-Finite Difference Software[END_REF]. For structures composed of a dielectric substrate with constant complex permittivity, a dispersive Debye model is designed by optimization process to approach it over a desired wideband.

II. HYBRID SPECTRAL-FDTD/TD-VFZ METHOD

A. Spectral-FDTD review

Let us consider a (l x , l y ) periodicity in the (x, y) directions with a z-propagation of a Gaussian incident wave. First, an This work is funded by the DGA AID and the CEA CESTA. arbitrary azimuth angle is chosen to ϕ = ϕ 0 . Then, the horizontal wavenumber depends on the angular frequency ω = 2πf and the elevation angle θ as k h = ω sin θ/c 0 where c 0 is the speed of light in free space. The term k h is fixed for one simulation of the SFDTD method. The constant components k x and k y of the plane wave vector k are deduced as a function of the azimuth angle and the constant horizontal wavenumber

k x = k h cos ϕ 0 , k y = k h sin ϕ 0 . (1)
By fixing them, simple time-domain PBC are found for Ψ = E or H Ψ(x, y, z, t) = Ψ (x + ml x , y + nl y , z, t) e j(kxmlx+kynly) ,

where m and n are any integers. Then, SFDTD scheme is exactly the same as the FDTD one with complex E-H fields unknowns instead of real for the latter. We deduce by the expression θ(ω) = arcsin k h k 0 that different frequencies correspond to different incident angles, and for a positive real incidence, the condition k 0 ≥ k h must be verified. In other words, this is referred to the plane wave region if f ≥ f c else to the evanescent wave region, where

f c = c 0 k h 2π (3)
is the cut-off frequency. When both the plane and evanescent wave region are excited simultaneously by the SFDTD method, surface waves are guided along horizontal directions at some frequencies. Thus, the energy of a horizontal wave continuously re-enters the computational domain through the PBC instead of being absorbed by the top and bottom convolution perfectly matched layers (CPML) [START_REF] Yang | A simple and efficient fdtd/pbc algorithm for scattering analysis of periodic structures[END_REF]. As a result, a resonance behavior occurs: temporal data do not converge and an accurate spectral response cannot be obtained by discrete Fourier transform (DFT).

B. TD-VFz method in Matlab post-processing

To obtain a correct spectral reflection coefficient R, the ratio between the input (excitation) and output (real part of the reflected field) response is decomposed with TD-VFz [START_REF] Lei | Efficient linear macromodeling via discretetime time-domain vector fitting[END_REF], as a rational transfer function

R(z = exp [jω∆ t ]) = d + M n=1 r n 1 -z -1 s n , ( 4
)
where d is a constant, r n and s n are the complex residues and poles respectively and ∆ t is the time sample. The TD-VFz extrapolation algorithm is an iterative finding poles method: conjugated complex poles are initialized inside the unit circle (linearly distributed near the unit circle boundary) and the new poles are deduced from an eigenvalue problem so-called 'poles relocation'. When poles have converged, the residues are obtained as the solution in the least squares sense of an oversized linear system. An interesting point consists in the eigen-frequencies identification. When (4) gives poles on the unit circle boundary, resonances occur and eigen-frequencies are identified.

C. Mapping R(0 < k h , f ) to R(0 < θ, f )
The normal incidence case θ = 0°is trivial and corresponds to k h = 0. To deal with oblique incident angles for the frequency band [0 < f min , f max ] and incident angle band [0 < θ min , θ max ] of interest, we first calculate the extremum horizontal wavenumber as

k hmax = 2πf max c 0 sin θ max , k hmin = 2πf min c 0 sin θ min .
(5) To compute the reflection coefficient as a function of the incident angle, we run the FDTD simulation for the N samples of k i h with the fixed azimuth angle ϕ = ϕ 0

k i h = k hmin + i k hmax -k hmin N -1 , i = 0, . . . , N -1, (6)
by using (1) to determine the phase values in the PBC formulation (2). Then, each k i h sample is extrapolated by the M-order TD-VFz algorithm to obtain the M residues r i n , the M poles s i n and the constant d i . After that, an incident angle-reconstruction is applied. For a given angle θ min < θ j ≤ θ max , we compute the N frequencies

f i = ck i h 2π sin θ j , i = 0, . . . , N -1. (7)
The reflection coefficient R evaluated at angle θ j is calculated with the transfer function (4) for all frequencies

f i R(f i , θ j ) = d i + M n=1 r i n 1 -exp [-j2πf i ∆ t ]s i n . (8)
As frequencies are imposed by the formula [START_REF] Chen | Wide-angle ultra-broadband metamaterial absorber with polarization-insensitive characteristics[END_REF], the Lagrange interpolation is used to obtain the reflection coefficient at the desired frequencies. Assuming a desired frequency f i ≤ f ≤ f i+1 , the reflection coefficient for the angle θ j and frequency f is given by

R( f , θ j ) = f i+1 - f f i+1 -f i R(f i , θ j ) + f -f i f i+1 -f i R(f i+1 , θ j ). (9)

III. DEBYE MODEL TO APPROXIMATE A CONSTANT COMPLEX PERMITTIVITY

A. One pole Debye model

Several metamaterials are specified by a constant complex permittivity = 0 ( r -j r ) , r , r > 0.

(10) [START_REF] Luebbers | Lossy dielectrics in fdtd[END_REF] proposed a dispersive one pole Debye model with a conductivity parameter to approach it over a desired wideband. Then, the Debye relative complex permittivity reads

r = ∞ + s -∞ 1 + jωτ 0 - jσ ω 0 = r debye -j r debye , ( 11
)
where ∞ is the infinite frequency permittivity, s is the static permittivity et τ 0 is the relaxation time. An efficient setting of (11) over a desired wideband f min to f max is obtained by adjusting the imaginary and real parts of (11) with (10) at the frequencies f min and f max . For example for f min

K r = ∞ + s -∞ 1 + (2πf min τ 0 ) 2 , (12)
K r = ( s -∞ ) 2πf min τ 0 1 + (2πf min τ 0 ) 2 + σ 2πf min 0 , (13)
where the approximation factor K > 1 should be close to one. Two other equations follow by replacing f min by f max and K by K 2 = 1/K in (12) and (13). We obtain a nonlinear system with four unknows ∞ , s , τ 0 and σ which can be solved with the 'fminsearch' Matlab function. To determine the optimum approximation factor K, 'fminsearch' is applied for several K between 1.001 and 1.2 with a step ∆ K = 0.001. K is determined when the error function E is minimum i.e

E( K) = min K E(K), (14)
with

E(K) = max (||A(K)|| L ∞ , ||B(K)|| L ∞ ), (15)
A(K) = r debye (K) -r , (16)
B(K) = r debye (K) -r (17)
where the norm L ∞ in ( 15) is taken on a 100 samples discrete frequency interval of [f min , f max ]. Finally, the dispersive one pole Debye model now obtained, the JE convolution (JEC) FDTD method [START_REF] Chen | An fdtd formulation for dispersive media using a current density[END_REF] is used to compute the EM fields in the unit cell.

B. Examples 1) Flame Resistant 4 (FR4):

The FR4 relative permittivity is given by r = 4.3 -j0.025.

(18)

For the FR4 case, we have r = 4.3 and r = 0.025. We want to find the four parameters ∞ , s , τ 0 and σ of (11) to approach (18) over the frequency band f min = 0.5 GHz to f max = 20 GHz. The minimum of the error function E is attained for K = 1.006 i.e our optimization process has given this optimum approximation factor where the Matlab function 'fminsearch' returns

F R4 ∞ = 4.26291, (19)

3) Polyethylene glycol terephthalate (PET):

The PET relative permittivity is given by r = 3.0(1 -j0.06). For the desired frequency band f min = 3 GHz to f max = 40 GHz, our optimization process has given the optimum approximation factor K = 1.045 where the Matlab function 'fminsearch' returns P ET ∞ = 2.774121, P ET s = 3.140665, τ P ET 0 = 6.647.10 -12 s and σ P ET = 0.02384704 S/m.

IV. NUMERICAL RESULTS

For the two numerical examples of this section, SFDTD solver runs with a 2 ns simulation time and TD-VFz with M = 20 poles and 10 iterations for the poles convergence.

A. PMI and PET absorber with resistive films

As shown in Fig. 3, the studied absorber [START_REF] Chen | Wide-angle ultra-broadband metamaterial absorber with polarization-insensitive characteristics[END_REF] is composed of three PMI substrates all separated by a PET film of thickness t p = 0.175 mm covered by a resistive film. Note also a metal ground simulated with perfect electric conductor (PEC) to prevent the transmission of the EM waves. The PMI and PET constant complex permittivities are approached with a one pole Debye model over a wideband from 3 GHz to 40 GHz (the parameters are listed in section III-B). The angular responses are deduced by the methodology of section II-C for N = 100 samples of horizontal wavenumber

B. FR4 absorber with ITO films

The studied absorber [START_REF] Wang | High-performance and ultra-broadband metamaterial absorber based on mixed absorption mechanisms[END_REF] is composed of a FR4 substrate, covered by thin indium tin oxyde (ITO) films. The unit cell geometry is plotted in Fig. 6. The pattern is a Huygens metasurface on three-layers slab impedance metasurface with a metal ground sheet (PEC) to prevent the transmission of EM waves. Dimensions and ITO film conductivities are optimized to promote a high absorption rate. The FR4 constant complex permittivity is approached with a one pole Debye model over a wideband from 0.5 GHz to 20 GHz (the parameters are listed in section III-B1). The angular responses are deduced by the methodology of section II-C for N = 100 samples of horizontal wavenumber k h with θ max = 60°. The FDTD spatial steps are all equal to ∆ = 0.25 mm. Note that to simulate the d = 50 nm thin conductive ITO films, an equivalent conductivity is applied on one FDTD cell i.e σ cell i = σ i .d/∆ z . Fig. 8 shows that this absorber can achieve an absorption rate better than 90 % over an ultra-broadband 1-18 GHz for the TM mode and for many incident angles. However for the TE mode, the absorption rate is a bit deteriorated at high incidence but remains superior to approximately 75 % over 1-18 GHz as shown in Fig. 7.

V. CONCLUSION

The Spectral FDTD scheme combined with the TD-VFz decomposition is an efficient hybrid approach to analyze periodic structures excited by oblique incident plane waves without reducing the CFL stability criterion of the FDTD method. The dispersive Debye model is successfully designed with optimization process to simulate two efficient absorbers.

Fig. 1 and

 1 Fig.1and Fig.2respectively show the real part and imaginary part of the Debye relative permittivity obtained with the parameters (19)-(22). Very good agreement with theory is observed over the desired frequency band.2) Polymethacrylimide (PMI): The PMI relative permittivity is given by r = 1.06(1 -j0.005). For the desired frequency band f min = 3 GHz to f max = 40 GHz, our optimization process has given the optimum approximation factor K = 1.004 where the Matlab function 'fminsearch' returns P M I ∞

Fig. 1 .Fig. 2 .

 12 Fig. 1. Real part of the FR4 Debye relative permittivity compared with theory.

Fig. 3 .- 1 .

 31 Fig. 3. Pattern geometry of the PMI (white) and PET (red) absorber with resistive films (blue) [7]. t 1 = 2 mm, t 2 = 4 mm, t 3 = 3 mm, tp = 0.175 mm, P = a 1 = a 3 = 14 mm, a 2 = 12.5 mm.

  0 θ " 40 0 θ " 60 0

Fig. 4 .

 4 Fig. 4. TE mode absorption curve for the PMI/ PET absorber.

Fig. 5 .

 5 Fig. 5. TM mode absorption curve for the PMI/PET absorber.

Fig. 6 .

 6 Fig. 6. Pattern geometry of the FR4 (blue) absorber with ITO films (yellow) [8]. h 1 = 21 mm, g = h 2 = 2.5 mm, w = 20 mm, h 3 = 3 mm, l 1 = 8 mm, l 2 = 35 mm, l 3 = 39 mm, σ 2 = 6.7 × 10 4 S/m, σ 1 = σ 3 = σ 4 = 2.0 × 10 5 S/m.

Fig. 7 .

 7 Fig. 7. TE mode absorption curve for the FR4/ITO absorber.

Fig. 8 .

 8 Fig. 8. TM mode absorption curve for the FR4/ITO absorber.


Susana Naranjo
email: susana.naranjo-villamil@edf.fr

Eric Piedallu
email: eric.piedallu@edf.fr

Julien Gazave
email: julien.gazave@edf.fr

Christophe Guiffaut
email: christophe.guiffaut@xlim.fr

Alain Reineix
email: alain.reineix@xlim.fr

Experimental Characterization of Reinforced Concrete Cable Ducts

Keywords: lightning, reinforced concrete, cable duct, impedance, current distribution

When a building is struck by lightning in an industrial facility, a partial lightning current can flow through the reinforcement of cable ducts, causing interference to the cables routing inside. The distribution of the lightning current among the paths leading away from the building depends on their impedance. In this paper, we present the impedance measured in two different cable ducts: one duct with a wire-tied reinforcement and the other with a welded reinforcement. The measurements are compared to the results of numerical simulations carried out using the FDTD method.

I. INTRODUCTION

Lightning direct strikes can have severe consequences in industrial facilities. To prevent potential differences and reduce the magnetic field, all conducting components in the facility are typically interconnected. Thus, when one of the buildings is struck by lightning, only a part of the current flows to the earth-termination system. The other part is distributed among the low-impedance paths leading away from the building. Moreover, interference to the cables routing inside reinforced concrete ducts can be caused by coupling to the reinforcement of the duct and by a partial lightning current entering the cables via the bonding network.

The portion of the lightning current flowing along the reinforcement of the duct depends on the installation conditions, such as the dimensions of the structures and the grounding arrangement. In general, calculating the current distribution among the paths can become a very complex problem, especially in large facilities. The international standard IEC 62305-1 [START_REF]Protection Against Lightning-Part 1: General principles[END_REF] assumes that half of the current goes straight to the building's earth-termination system and proposes a simple model to distribute the remaining half. The German standard KTA 2206 [START_REF]Design of Nuclear Power Plants Against Damaging Effects from Lightning[END_REF] specifies weighting factors to calculate the partial lightning current flowing through the reinforced concrete ducts and the soil-contacting conductors of the building, assuming that one-third of the current flows into earth. Also, it proposes a formula to calculate the axial voltage induced in cables routing inside the ducts, from the partial lightning current flowing through their reinforcement. The formula considers an inductive transfer mode and uses the coupling impedance obtained experimentally in [START_REF] Zischank | Shielding effectiveness of reinforced concrete cable duct carrying partial lightning currents[END_REF]. An example of the application of the KTA approach is given in [START_REF] Zischank | Assessment of the lightning transient coupling to control cables interconnecting structures in large industrial facilities and power plants[END_REF].

Since the weighting factor of reinforced concrete ducts is based on their inductance, the uncertainties of the implemented models with respect to actual constructions could have a significant influence on the results. Aiming to quantify the uncertainties, in this paper, we compare the impedance measured in various sections of two different reinforced concrete ducts and the results of numerical simulations.

II. II. EXPERIMENTAL SETUP

As a first step towards validating the numerical model of reinforced concrete ducts, we measured cable ducts at the new lightning testing facility of EDF Power Networks Lab.

A. Cable Ducts

The testing facility has two parallel reinforced concrete ducts, as shown in Fig. 1. The steel plates covering the ducts were removed during the measurements. The ducts are both 60-meters-long with a cross-sectional area of 1 m × 0.8 m. They are divided into three electrically isolated sections of 20 m, named S1 to S3. Their reinforcement is made up of a double-layered reinforcing grid embedded in 20 cm of concrete. The rebars have a radius of 5 mm, and they form squared meshes of 20 cm x 20 cm. In the first duct (C1), the rebars are wire-tied. In the second (C2), they are welded. The layers of the reinforcing grids are separated by 15 cm, and they are interconnected, as shown in Fig. 2. Note that all these values are approximated; they are more likely to correspond to an average of the real dimensions.

The expansion joints between the sections of the ducts are filled with foam. Each section has 16 fixed grounding terminals: eight at the beginning and eight at the end. The reinforcement of adjacent sections can be interconnected using bridging braids, as shown in Fig. 3. To minimize unwanted effects at the boundaries and aiming to inject a current uniformly into the reinforcement, the ducts are terminated by copper plates, as shown in Fig. 4.

B. Impedance Measurements

The impedance of a structure can provide useful information regarding its electromagnetic behavior. Using a vector network analyzer (VNA), we measured the impedance of each section of the ducts, then the impedance of two interconnected sections, and finally the impedance of three interconnected sections. Since the magnitude of the impedance was expected to be low, two different methods were used: One-Port reflection, also known as Shunt measurement, and Shunt-Thru. The connection setups are shown in Fig. 5.

Regardless of the method and the section of the duct being characterized, the length of the wires used to inject the current and create a return path was constant. Their impedance was measured independently to remove their contribution from the results. Also, all the measurements were repeated using One-Port reflection after performing a new open/short/load calibration in which the wires were taken into account. To differentiate the results, we named the method One-Port*.

Certainly, removing the covering plates could have a significant effect on the results, especially if they were electrically connected to the reinforcement. The values measured are not to be scaled-up and compared directly to the models of enclosed cable ducts one can find in the literature [START_REF] Zischank | Shielding effectiveness of reinforced concrete cable duct carrying partial lightning currents[END_REF], [START_REF] Heidler | Induced overvoltages in cable ducts taking into account the current flow into earth[END_REF]. The plates were removed to focus exclusively on the characteristics of the reinforcement and its representation in numerical simulations. Frequently, the rebars forming the reinforcement of cable ducts cannot be welded consistently at every intersection. One then wonders if it is necessary to consider a reinforcing grid with bigger meshes in the numerical model.

III. NUMERICAL MODEL

Considering a reinforced concrete duct with the dimensions given in section II.A, different simulations are carried out using TEMSI-FD [START_REF] Temsi-Fd | Time ElectroMagnetic Simulator -Finite Difference Time Domain[END_REF], a full-wave solver based on the finite-difference time-domain (FDTD) method. All the conductors are modeled as thin wires [START_REF] Guiffaut | New oblique thin wire formalism in the FDTD method with multiwire junctions[END_REF] with a conductivity of 59.6×10 6 S/m. Except for the rebars forming the reinforcement, their conductivity is set to 8.33×10 6 S/m. The permeability of the steel is neglected. The concrete is modeled as a lossy dielectric material with a conductivity of 0.0052 S/m and a relative permittivity of 8.6. These values were obtained by fitting a Debye model to the data of a concrete sample with a moisture content of 5.5% [START_REF] Sandrolini | Modelling the electrical properties of concrete for shielding effectiveness prediction[END_REF]. The soil is considered homogeneous with a relative permittivity of 10 and a resistivity of 100 Ωm.

The impedance of the duct is calculated by dividing the voltage of a perfect voltage source over the current injected and then removing the effect of the return path. However, we observed that even if we try to remove the effect, the imaginary part of the impedance varies slightly with the position of the wires. Probably, it also happens when making measurements.

The voltage waveform is defined as a Gaussian function covering a frequency range from 0 Hz to 10 MHz. The maximum amplitude of the pulse is set to 100 V. All the external surfaces of the analysis space are defined as PMLs [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] to assume an open space. The computation model is shown in Fig. 6.

IV. RESULTS

The measurements and simulations in Fig. 7 show that the reactance is predominantly inductive, and as expected, it does not vary significantly. All curves are well fitted with an inductance of approximately 35 µH. On the other hand, there is a considerable difference in the resistance. The results of the simulations are closer to the impedance measured in C2. The latter could be explained by the contact resistance of the wiretied rebars in C1. Nevertheless, it is clear from the results in Fig. 7(b) that there is an error at low frequency. The magnitude of the resistance measured in C2 may already be within the range of uncertainty of the measurement method. From the impedance measured on each section of the ducts, the resistance and inductance were calculated. The average of the values obtained with the different methods is presented in Table I. Even though the values differ between two sections of the same duct, the sum of the resistances and the inductances are not far from the values calculated with the sections interconnected. Moreover, the parameters of a transmission line model of the cable duct, as proposed in [START_REF] Heidler | Induced overvoltages in cable ducts taking into account the current flow into earth[END_REF], could easily be calculated.

Differences between the two sections of the same duct are expected since, in reality, the ducts are not divided into three equally long sections. However, it is still interesting to observe the significant variations of the resistances calculated for the sections of C1. When wire-tying the rebars, it is complex to control the quality of the contact, and therefore, to define a load to adjust the values in numerical simulations of industrial facilities.

Nevertheless, adjusting the values may not be necessary because the reactance of the duct quickly becomes higher than its resistance. Thus, with a transmission line model of the ducts and the specifications of the grounding system and the bonding network, the partial current that may flow through the reinforcement could be approximated. The question that remains is if, from that partial current, one can accurately calculate the surges in the cables routing inside.

V. CONCLUSION

The fact that the rebars of reinforced concrete ducts may not always form perfectly symmetrical grids or be welded, has a negligible influence on the reactance of the ducts. Thus, if the distribution of the lighting current between the ducts depends mostly on their inductance, the uncertainties of the numerical models with respect to the constructions could be neglected.

The inductance of a reinforced concrete duct can be well approximated considering the average dimensions of the meshes of the reinforcing grid in full-wave simulations. However, including the resistance of the duct in the simulations could be tricky if the rebars are wire-tied. We observed that in the 60-meters-long ducts at the testing facility, the resistance was 20 times higher when the rebars were wire-tied. Besides, poor contact between the rebars could create non-linearities.

Fig. 1 .

 1 Fig. 1. Reinforced concrete ducts in the testing facility.

Fig. 2 .

 2 Fig. 2. Front view of the reinforcement of the ducts.

Fig. 3 .

 3 Fig. 3. Interconnection of two sections of the concrete duct.

Fig. 4 .

 4 Fig. 4. Copper termination plates.

Fig. 5 .

 5 Fig. 5. Measurement setups.

Fig. 6 .

 6 Fig. 6. Computation model of the reinforced concrete duct.

Fig. 7 .

 7 Fig. 7. Impedance of the reinforced concrete duct. (a) Measured in C1. (b) Measured in C2. (c) Calculated using TEMSI-FD.

TABLE I .

 I RESISTANCES AND INDUCTANCES CALCULATED FROM THE IMPEDANCE MESURED ON DIFFERENT SECTIONS OF THE DUCTS

	Section	R (mΩ)	C1	L (µH)	R (mΩ)	C2	L (µH)
	S1	59.652		8.962	8.175		8.189
	S2	217.564		19.663	19.501		13.790
	S3	87.399		13.599	53.291		11.938
	S1 S2	258.329		25.832	12.474		16.934
	S2 S3	249.848		27.808	14.999		28.742
	S1 S2 S3	302.034		34.833	15.851		35.027


Susana Naranjo-Villamil
email: susana.naranjo-villamil@edf.fr

Christophe Guiffaut
email: christophe.guiffaut@xlim.fr

Julien Gazave
email: julien.gazave@edf.fr

Alain Reineix
email: alain.reineix@xlim.fr

On the Calculation of Electrical Surges in Underground Cables due to a Direct Lightning Strike

Keywords: lightning, lightning protection, surges, industrial facility, finite-difference time-domain

When a building is struck by lightning, the lightning current travels through its reinforcement and along the lightning channel, generating a transient electromagnetic field. A part of the current reaches the ground termination system and the other part is distributed among the cable ducts and soil-containing conductors leading away from the building. Electrical surges can be induced by the transient fields, the coupling between the structures, and a partial lightning current entering the cables via its grounding. Since the Lightning Protection System (LPS) is generally designed based on the worst-case scenario, a parametric study is conducted using the FDTD method to identify the configurations in which the surges in the cables interconnecting two buildings are maximized. The results are compared to the case in which the cables are grounded at the entrance to the buildings. The direct connection to the reinforcement increases the currents and shifts the resonances towards higher frequencies.

I. INTRODUCTION

Electrical surges are considered as one of the most harmful indirect effects of lightning. They can cause permanent physical damage and the malfunctioning of sensitive devices. These surges are generally induced in the cables nearby by the transient electromagnetic field and the ground potential rise generated by a lightning strike. During a direct strike, the current travels along the lightning channel and through the struck object; yet, only a part of the current goes straight to the earth termination system. The other part is distributed among the cable ducts and soil-contacting conductors leading away from the striking point. Therefore, electrical surges can also be induced by the coupling of the cables to the struck object itself and the ducts, or by a partial lightning current entering the wires via their grounding. An interesting example of the lightning current distribution in the grounding system and the cable networks of a radio base station is presented in [START_REF] Grcev | Lightning Current Distribution to Ground at a Power Line Tower Carrying a Radio Base Station[END_REF].

To appropriately design the Lightning Protection System (LPS) and the protection measures against the Lightning Electromagnetic Pulse (LEMP) in industrial facilities, electrical surges need to be estimated. Numerous studies have been carried out to deepen the understanding of the phenomenon [START_REF] Heidler | Induced overvoltages in cable ducts taking into account the current flow into earth[END_REF]- [START_REF] Tatematsu | Application of the FDTD-based simulation code VSTL REV to the lightning surge analysis of a nuclear power plant[END_REF], and one can find simplified methods to assess overvoltages in international standards, see e.g. [START_REF]Protection against lightning-part 4: Electrical and electronic systems within structures[END_REF]. However, experience has shown that the methods are not always adapted to the characteristics of French industrial facilities [START_REF] Duqueroy | Lightning risk assessment evaluation on French nuclear power plants[END_REF]. Since modeling all the elements of a full-scale industrial facility at once is almost impossible, this paper aims to identify the parameters that could be neglected by considering a worst-case scenario. We focus on parameters that are commonly assumed to be invariable, and may still have an important influence on the results. The study is conducted on a simplified industrial facility: two reinforced concrete buildings interconnected by a reinforced concrete cable duct, where two cables are laid. First, the cables are connected to an internal grounding system and then to the reinforcement at the entrance to the buildings. The results are compared to two worst-case scenarios.

Even though the Transmission Line (TL) theory is still widely used to calculate surges due to the simplicity of circuit-based simulations, the complexity of the installations and the need for more accurate results have made the Full Wave methods gain popularity. Hence, the simulations of the industrial facility are carried out using TEMSI-FD (Transient ElectroMagnetic Simulator -Finite Difference) [START_REF]TEMSI-FD: Time Electromagnetic Simulator -Finite Difference Time Domain[END_REF], a solver based on the finite-difference time-domain (FDTD) method. A non-uniform grid is implemented with cell sizes varying from 25 cm to 2.5 m, and all the conducting wires, including the cables, the rebars, and the channel, are modeled as thin wires [START_REF] Guiffaut | New oblique thin wire formalism in the FDTD method with multiwire junctions[END_REF].

II. REFERENCE CASE

To study the surges resulting from a direct lightning strike in an industrial facility, a simplified case is defined: two 20meters-tall reinforced concrete buildings facing each other, with a horizontal cross-sectional area of 20 m × 50 m. The external walls, the roof, and the foundation of the buildings are made up of a single-layered reinforcing grid, given that a multi-layered reinforcing grid can be replaced by an equivalent single-layered grid [START_REF] Naranjo-Villamil | Simplified 3-D modeling of reinforced concrete for the calculation of transient electromagnetic fields inside a building struck by lightning[END_REF]. The grids have a squared mesh size of 50 cm × 50 cm and are embedded in 25 cm of concrete, which is modeled as a lossy dielectric material with a conductivity of 0.0052 S/m and a relative permittivity of 8.6. The rebars forming the grid have a radius of 1 cm, and their conductivity is set to 8.33×10 6 S/m. The foundations are buried five meters below the surface. No internal walls, columns, or beams are considered. A 100-meters reinforced concrete duct interconnects the buildings in the middle, as shown in Fig. 1. The duct has a cross-sectional area of 3 m × 2.5 m and is positioned 50 cm above the foundations of the buildings. The reinforcing grid of the duct has a mesh size of 25 cm × 25 cm and is also embedded in 25 cm of concrete. Its rebars have the same conductivity but a radius of 8 mm. The external grounding system consists of two grounding rings, one around each building, connected to a grounding grid with squared meshes of 20 m × 20 m. The conductors are buried 87.5 cm below the surface, have a radius of 7.5 mm, and are made of copper. Accordingly, their conductivity is set to 59.6×10 6 S/m. The soil is considered homogeneous with a relative permittivity of 10 and a conductivity of 0.01 S/m. The buildings are connected to the grounding rings at each corner, as shown in Fig. 2. The internal grounding system of each building consists of collective ground conductors forming five rings. The grounding rings are interconnected in the corners, as shown in Fig. 3. The first ring is also connected in the corners to the reinforcement of the building and thus directly to the external grounding system.

There are two wires in parallel in the duct: an insulated wire loaded with 50 Ω at the beginning and the end, and an auxiliary copper wire. The wires follow the same path inside the buildings but are positioned differently inside the duct (see Fig. 4). They are connected to the fifth grounding ring inside the buildings. The distances from the wires to the walls are 5 m on the x-axis and 75 cm on the y-axis. The cable shelves and the support rails inside the duct are not modeled. Each cable is divided into nine segments named S1 to S9, respectively. Segments S1 to S4 are the ones inside the first building, segment S5 is inside the duct, and segments S6 to S9 are inside the second building. A detailed view of the path followed inside the first building is also shown in Fig. 3. The currents are computed in the middle of each segment.

The lightning channel is represented as a monopole antenna by a 200-meters vertical lossy wire, excited at its base by a lumped current source and connected at the top end to a perfectly matched layer (PML) [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. Aiming to conduct the study in the frequency domain, the current waveform is defined as a Gaussian function of 1 A, covering a frequency range from 0 Hz to 1 MHz. Its magnitude is attenuated by 10 at 1 MHz and an attenuation of 10 6 is set at the time instant t = 0 s to avoid any noise in the response. All the results in the frequency domain are normalized to this source, i.e., they are divided by the spectrum of the Gaussian pulse. Thus, they are independent of the current waveform.

A. Extrapolation with the Matrix Pencil Method

Even though the Gaussian pulse has a duration of just a few microseconds, the energy had not been completely dissipated after 50 µs due to the presence of resonances. The resonances are most likely associated with the total length of the path to the ground termination system. Since a single simulation up to 50 µs, running in parallel in 10 CPUs was taking 22 hours, instead of increasing the simulation time, we decided to extrapolate all the results applying the Matrix Pencil method [START_REF] Sarkar | Using the matrix pencil method to estimate the parameters of a sum of complex exponentials[END_REF]. The method consists of modeling the late-time response of a system as a sum of complex exponentials and estimating the parameters of the model using a linear technique. From the poles and the residues obtained after solving the optimization problem, the response can be well approximated in both time and frequency domains. Fig. 5 shows the extrapolation of the current in segment S1 of the auxiliary copper wire in the time domain. The normalized currents in all the segments in the frequency domain are shown in Fig. 6.

To estimate electrical surges, the cables are often considered to be grounded at the beginning and the end of the duct. It can be observed in Fig. 7 that if both wires in the reference case are connected to the reinforcement at the entrance to the buildings, instead of being connected to the internal grounding system, the currents rise significantly and the resonances are no longer present below 1 MHz. If only the auxiliary wire is grounded at the entrance to the buildings, the current in the insulated wire increases around 2 dB at low frequency and decreases at the frequency of the first resonance.

III. PARAMETRIC STUDY

Due to the complexity of industrial facilities, the effects of a direct lightning strike are influenced by numerous parameters. In this study, we decided to focus on parameters that are rarely mentioned in lightning protection standards. The objective is to define a worst-case scenario from which, eventually, a more generalized parametric study could be conducted. The parameters are varied independently, and the results are analyzed at low frequency and the frequency of the first resonance. Even though the effects of the resonances are reduced, if not completely erased, when using the waveforms of the first positive and the first negative return strokes, they

TABLE I CASES IN WHICH THE HIGHEST NORMALIZED CURRENTS ARE OBSERVED FOR EACH PARAMETER

No. Case

1
One connection from the internal to the external grounding system in corner C2.

2

One connection between the internal grounding rings in the second building.

3 Cables connected to the first grounding ring in the first building.

4

Cables connected to the first grounding ring in the second building.

5

One connection between the first building and the external grounding grid, in corner C1.

6

The external grounding system consists of interconnected rings. 7

Strike in corner C1.

8

Strike in corner C2 and duct positioned 1 m away from the edge.

may create a major difference when using the waveform of the subsequent stroke. Table I summarizes the cases from which the highest normalized currents are obtained for each parameter. The corresponding frequency responses in segment S1 are shown in Table II and Fig. 8.

A. Internal Grounding System

Significant variations in the currents are not observed when modifying the configuration of the internal grounding system. Nor are there significant variations when reducing the number of connections from the internal to the external grounding system in the second building. In the struck building, the currents are slightly higher if there is a reduction in the number of connections to the external grounding and the reinforcement. A small rise is also obtained when the internal grounding rings are connected once, rather than at each corner. Since the currents collected will converge to one point, the coupling between the rings and the wires may be stronger. In both cases, the length of the path to the ground is modified, and consequently, the frequency of resonance is shifted.

B. Connection of the Wires to the Internal Grounding System

At low frequency, as long as the wires are relatively far from the walls and they are connected to the same internal grounding ring, the changes are negligible. As the connection is moved away from the roof, the currents start decreasing and then increase when it approaches the foundation. The highest currents are observed when the wires are connected to the first grounding ring inside the first building. Slightly higher currents are also found when the length of the cables and thus their impedance is reduced. On the contrary, at the frequency of the first resonance, the effect of the distance to the walls is more pronounced, and the highest currents are observed when the wires are connected to the first grounding ring inside the second building.

C. External Grounding System

The effect of modifying the configuration of the external grounding system, or varying the number of connections to the buildings, is only noticeable at low frequency. The portion of the current flowing to the grounding system increases when its impedance decreases, resulting in a lower current flowing through the duct and consequently a decrease in the currents induced in the cables. Hence, the highest currents are observed when the grounding system consists only of two interconnected grounding rings, one around each building. Naturally, the currents also increase if the buildings are not connected to the external grounding system. Fortunately, there is always at least one connection. The analysis is then similar to the one made for the interconnection of the internal grounding rings: the coupling may be stronger if the number of connections is reduced and placed close to the wires' path.

D. Striking Point

Instead of striking the building in the center of the roof, the lightning channel is attached to the corners. When the wires are connected close to the striking point, the highest currents are observed at low frequency and the lowest at the frequency of the first resonance. The increase is comparable to the effect obtained when the building is struck in corner C2, probably because there is a higher current running through the duct and, therefore, the coupling is stronger. The highest currents at the frequency of the first resonance are observed in the reference case, suggesting that the symmetry reinforces the phenomenon.

E. Position of the Duct

We consider the duct to be positioned at a distance to the edges varying between 1 m and 24 m. The lightning channel is attached first to the center of the roof and then to the corner C2. As expected, the closer the duct is to the edges and to the striking point, the bigger the portion of the current that flows through it. Additionally, changing the position of the duct also changes the path the wires follow inside the building; thus, at low frequency, the highest currents are obtained when the building is struck in the corner because of the proximity of the down conductor in the edge to segment S4. At the frequency of the first resonance, the values are still lower than in the reference case.

IV. WORST-CASE SCENARIO

From the parameters studied, the one that increases the surges the most is the configuration of the external grounding grid, followed by the position of the duct. Yet, the place where the different elements are connected can play an important role. Accordingly, a first scenario is defined, including the following differences to the reference case:

• The building is struck in corner C2.

• The duct is placed one meter away from the down conductor in corner C2. • The external grounding system consists of two interconnected grounding rings, one around each building. • The first building and its internal grounding system are connected to the external grounding system only once, passing by the corner C1. • The rings of the internal grounding system of the first building are also connected once. The connection is close to the striking point. • The wires are connected to the first grounding ring in each building, one meter away from the down conductors in corners C1 and C1b. The connection in the second building is close to C1b to reduce the cable length. In the first building, it is close to C1 to strengthen the coupling between the wires and the other conductors. Given that the highest currents at the frequency of the first resonance were observed when the configuration was symmetrical, a second scenario is defined with the following differences:

• The external grounding system consists of two interconnected grounding rings, one around each building. • The first building and its internal grounding system are connected to the external grounding system only once, passing by the corner C2. • The rings of the internal grounding system of the first building are also connected once. The connection is close to C2. • The wires are connected to the first grounding ring in each building, one meter away from the down conductors in corners C2 and C1b. An upper view of the configurations is presented in Fig.

The results in frequency domain are shown in Fig. 10. From these curves, the currents in the time domain, for any slow waveform can be obtained, by multiplying the spectrum of the desired waveform and computing the inverse Fourier transform. They are transfer functions. As an example, Fig. 11 shows the currents computed for the waveform of a first positive stroke and a subsequent negative stroke as defined in [START_REF]Protection against lightning-part 4: Electrical and electronic systems within structures[END_REF], protection level I. Since the waveform of the subsequent stroke has a short rising time, to compute the currents from the transfer functions, the simulations and the extrapolation were repeated using a shorter Gaussian pulse covering a 10 MHz bandwidth.

In the first scenario, the currents increase by 13.13 dB in the insulated wire and 17.03 dB in the auxiliary copper wire at 10 kHz. At the frequency of the first resonance, they increase by 7.26 dB and 8.25 dB, respectively. The currents are lower at 10 kHz in the second scenario; however, they are around 10 dB higher at the frequency of the first resonance. Compared to this increase, the rise in the currents observed in Fig. 7, when connecting both wires to the reinforcement in the reference case, no longer appears to be a significant overestimation. Fig. 10 also shows that if the wires in the worst-case scenarios are grounded at the entrance to the buildings, the currents at low frequency reach similar values. Hence, as long as the external grounding system is modeled appropriately one could estimate the currents in the worst-case scenarios by connecting the wires to the reinforcement. Nevertheless, Fig. 12 shows that it would reduce the resonant behavior when computing the time response for a fast waveform. Also, the currents may still be overestimated for slow waveforms if the wires are not loaded.

V. CONCLUSION

Lightning-induced surges in industrial facilities are strongly dependent on the configuration of the external grounding system, the position of the ducts with respect to the striking point, and the path the wires follow inside the buildings. To design the LPS, it would be preferable to consider at least two worst-case scenarios, one for slow excitations and another in which the resonances are intensified. Under the simplifications and the conditions defined in the study, the currents calculated when the cables are connected to the reinforcement at both ends of the ducts can give a good idea of the currents expected in the worst-case scenarios. Moreover, they could be used to estimate the surges. Further work is required to establish if the approach can be generalized. By grounding the cables at the entrance of the buildings, all the parameters associated with the internal grounding system could be disregarded in the model.

Fig. 1 .

 1 Fig. 1. Calculation model of the reference case.

Fig. 2 .

 2 Fig. 2. Upper view of the first building in the reference case.

Fig. 3 .

 3 Fig. 3. Detailed view of the first building in the reference case.

Fig. 4 .

 4 Fig. 4. Position of the wires inside the reinforced concrete duct.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Extrapolation of the current in segment S1 of the auxiliary wire.

Fig. 7 .

 7 Fig. 7. Normalized currents in segment S1 when the wires are grounded at the entrance to the buildings. (a) Insulated wire. (b) Auxiliary wire.

Fig. 8 .

 8 Fig. 8. Cases in which the highest normalized currents are observed in segment S1. (a) Insulated wire, cases 1 to 4. (b) Insulated wire, cases 5 to 8. (c) Auxiliary wire, cases 1 to 4. (d) Auxiliary wire, cases 5 to 8.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Upper view of the calculation models. (a) First scenario. (b) Second scenario.

Fig. 11 .

 11 Fig. 11. Currents in segment S1 in the worst-case scenarios. (a) Insulated wire. (b) Auxiliary wire.