I would like to express my deepest appreciation to the examination committee for taking the time to be part of this journey. Special thanks to professor Farhad Rachidi and professor Marc Hélier, who accepted to report on my thesis. Many thanks to professor Françoise Paladian for her efficiency as the chair of the committee, to Dr.

Jean-Philippe Parmantier for his valuable and constructive comments, and to Dr.

Alexandre Laisné for engaging with my thesis enthusiastically. Thanks also to Charlotte Trouilloud, who played a decisive role in defining the scope of this work.

I cannot begin to express my gratitude to my supervisors professor Alain Reineix, Dr. Christophe Guiffaut, and Dr. Julien Gazave. I will never be able to make up for my terrible speech when the jury awarded me the Ph.D. degree, but please know that I meant it when I said that I couldn't have done it without you. None of this would have been possible without the guidance and the support of every single one of you.

Alain, since I met you I was amazed by how someone so bright and knowledgeable could be so humble. Talking to you felt effortless. You got me, you understood my thinking process, and you always offered helpful advice. Christophe, you contributed the most to my success. I cannot think of a single time you didn't make yourself available to discuss my ideas or answer my questions. Thank you for your outstanding patience. Julien, I'm indebted to you for believing in me from the start. You were always there, cheering for me, even during tough times. I would also like to extend my sincere thanks to Dr. Olivier Tantot, Dr. Claire Dalmay, and all the professors who allowed me to be part of their team. Teaching has always been one of my greatest passions and it was such a pleasure to work with you. Olivier, I will never forget that day you came into my office, it was past 7 pm, and I had had a particularly bad day. I don't remember exactly what you said but I do remember it made me feel exceptional and worthy of being there. It wasn't the only time your kindness made my day. Thank you. I am also grateful to Nicolas Bui and Eric Piedallu, who made my time much more enjoyable. Nico, I was able to settle in quickly thanks to you. I very much appreciate all the help and practical suggestions you provided. Eric, your positive attitude is infectious. Thanks for always being willing to lend a hand. Not only are you a valuable teammate, but you are an amazing friend.

I would like to recognize the assistance that I received from Nadine Aymard and Claire Lemarchand. Thanks should also go to Elie Durcik for his vote of confidence, to the National University of Colombia for giving me the opportunity to come to France in the first place, and to all those who made the university one of my favorite places in the world. List of Tables

Preface

From mechanical to electrical effects, a direct lightning strike to an industrial facility can cause severe damage. Knowing what to expect is not only fundamental to design the lightning protection system of the facility, but also necessary to achieve reliability and safety.

When a reinforced concrete building in the facility is struck by lightning, the lightning current flows along the lightning channel and through the reinforcement of the building, generating a transient electromagnetic field. Only a part of that current goes straight to the earth-termination system, the other part is distributed among the ducts and the conductors leading away from the building. The transient electromagnetic field could lead to an upset of electric installations and sensitive electronic devices inside the building. In addition, the electrical surges it induces in the cables nearby may result in physical damage. Electrical surges can also be induced by the coupling of the cables to the reinforcement of the building, the reinforcement of the cable duct, or by a partial lightning current entering the conductors via their grounding.

One can find valuable principles in the existing literature to protect structures and equipment against the adverse effects of lightning. However, experience has shown that in general, the methods suggested by current standards are not adapted to the characteristics of French industrial facilities. Also, the principles on which they are based cannot be used in the safety demonstration. There is hence a growing need for accurate and validated tools to estimate the effects of lightning strikes and assess the risk associated.

In this thesis, we have decided to focus on the indirect effects of direct lightning strikes on large industrial facilities. More specifically, on estimating the peak-values of the transient magnetic field generated by a direct strike inside industrial buildings;

and the surges induced in underground cables when one of the buildings in the facility is struck by lightning. We carry out simulations using TEMSI-FD, a full-wave solver based on the finite-difference time-domain (FDTD) method. Then, using the data obtained from the simulations, we build different surrogate models to create a set of engineering tools to calculate the peak-values. Full-wave simulations can be computationally expensive and time-consuming. Sometimes all one needs is a quick yet reasonable estimation of the values.

This thesis is organized as follows:

Chapter 1 gives a brief overview of lightning and lightning protection principles and introduces the industrial context, which defines the main objectives of the work.

Chapter 2 presents the numerical approach and techniques that will be used to reach the objectives. It includes a theoretical background on full-wave simulations and surrogate modeling. Chapter 4 is dedicated to estimating the peak-values of the electrical surges in underground cables. Similar to chapter 3, a parametric study is conducted to identify the components of the electromagnetic environment that must be considered in numerical models. Then, a dataset is built from a worst-case scenario, and twelve surrogate models to calculate the peak-values of the surges are proposed.

Finally, the conclusions of the thesis and the perspectives for further research are presented in chapter 5.

Chapter 1 Introduction

Lightning

Lightning is a transient phenomenon that has always fascinated mankind; yet, it has only started to be understood since Franklin's kite experiments. We now know, partly because of Franklin, that clouds are electrified and lightning is nothing but an electrical discharge. Although, he was not the first person to see the similarities between lightning and sparks [1].

Many features of lightning could be derived from the studies of long sparks in high-voltage laboratories [2]. However, it is worth noting that there are significant differences between what happens in nature and what we can reproduce in a laboratory. From all the definitions one can find in the literature, the one given by Komarek in [3] stands out to show us how much there is left to learn: "Now lightning is but an expression, spectacular to be sure but still only an expression, of a vital force or process in nature".

The primary source of lightning is thunderstorms, in which the downdrafts and updrafts cause the separation of electrical charges in clouds [3]. As shown in figure 1-1, the idealized charge structure of thunderclouds consists of three charge regions:

one main positive (P), one main negative (N), and a small positive at the base of the cloud (LP) [4]. In general, the amount and the density of lightning increase with the size of the thundercloud [2]. A lightning discharge can be classified into two categories: ground flashes and cloud flashes, depending on whether a charge is brought down to earth or not. Ground flashes constitute about 25% of global lightning activity [4] and can be classified into four categories [START_REF] Cooray | Lightning parameters of engineering interest[END_REF]: downward negative ground flashes, downward positive ground flashes, upward positive ground flashes, and upward negative ground flashes. Upward and downward, indicate the direction of propagation of visible leader channels, and positive and negative indicate the polarity of the charges carried by those leaders [2].

About 90% of ground flashes are downward negative flashes [4], i.e., with a leader traveling from cloud to ground and bringing down a negative charge.

The leader travels in a stepped manner, creating a path from the negative charge region in the thundercloud to the ground. When it approaches the ground, a connecting leader may rise from a grounded structure, bridging the gap to the ground.

As soon as the stepped leader and the connecting leader meet, the potential of the tip of the stepped leader changes from cloud potential to ground potential and the first return-stroke process begins [START_REF] Cooray | Lightning parameters of engineering interest[END_REF].

The return-stroke is a wave of near-ground potential traveling upwards along the channel to neutralize the leader charge. Each flash typically contains 3 to 5 strokes [4]. All the strokes that propagate along the formed channel after the first returnstroke are called subsequent return-strokes. Since the current of the return-stroke is thought to produce most of the damage associated with a lightning strike [START_REF] Rakov | Fundamentals of Lightning[END_REF], lightning protection systems are generally designed based on its parameters.

Lightning protection principles

Lightning interacts with humans and objects in different ways, depending on whether we consider a direct strike or a nearby lightning strike. A direct strike to a human, as well as a direct strike to an object in the vicinity, can be lethal. Also, a direct strike to an object could have electrical and mechanical effects. For example, a direct strike to a tree can lead to an explosive expansion that may shatter the tree, but it can also create a potential difference between the trunk of the tree and a person standing nearby. If the potential difference is large enough, a side flash may occur, i.e., a portion of the lightning current may flow to the ground through the person.

One can also receive side flashes from other structures which are not equipped with a proper lightning protection system (LPS) [START_REF] Cooray | An Introduction to Lightning[END_REF].

Similarly, during a lightning strike to a structure, the potential difference between a conductor through which lightning current is flowing and a person at ground potential touching the conductor would result in a current flowing through the person.

The potential difference is called a touch voltage. In addition, the flow of lightning current through the conductors could melt them or cause mechanical stress in the structure [START_REF] Cooray | An Introduction to Lightning[END_REF].

Another effect of direct lightning strikes is electrical surges. 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. A portion of the lightning current may therefore enter the electrical system, creating sparks and a dangerous voltage rise. 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].

Electrical surges can also be induced by coupling to the struck object, and by the transient electromagnetic field and the ground potential rise generated by both direct and nearby strikes. Also, direct and nearby strikes can cause injuries in humans due to step voltages. The flow of the current in the ground creates a voltage that decreases with the distance to the striking point. Thus, if a person is standing with one foot directed towards the striking point and the other away from it, the potential difference between the feet could cause a current to flow through the person [START_REF] Cooray | An Introduction to Lightning[END_REF].

To minimize the unwanted effects of a lightning strike, there is a large variety of measures that could be taken,including a large variety of equipment that could be installed. Different standards (see e.g. [10,[START_REF]IEC 62561: Lightning protection system components (lpsc)[END_REF][START_REF]KTA 2206: Design of Nuclear Power Plants Against Damaging Effects from Lightning[END_REF][START_REF]NFPA 780: Standard for the installation of lightning protection systems[END_REF][START_REF]Recommended practice for equipment grounding and bonding in industrial and commercial power systems[END_REF]) suggest different approaches.

Yet, most of them agree on the necessity to have an external lightning protection system to safely conduct the current to the ground, and an internal lightning protection system to limit the surges. Generally, when it comes to protecting a structure, the external lightning protection system consists of an air-termination system, down conductors, and an earth-termination system (see figure 12). The air-termination system intercepts lightning flashes, the down-conductors transport the lightning current to the ground, and the earth-termination system helps to dissipate the current transported by the downconductors into the ground.

External lightning protection system

As in figure 1-2, an air-termination system can simply consist of the standard Franklin's lightning rods, also known as air terminals or strike-termination devices.

Or, it could also be, among other arrangements, a metal mesh covering the upper surface of the structure, as shown in figure 123. The latter is called the mesh method.

There are several methods to determine the arrangement and the position of the airtermination system. The rolling sphere method proposed by the IEC 62305-3 standard [START_REF]3: Physical damage to structures and life hazard[END_REF] is the universal method of design. Regardless of the method, air-termination components should preferably be mounted at corners and edges. To create parallel current paths, there must be several equally spaced downconductors of minimum length. The number of down-conductors, their characteristics, and their position, depend on the structure. When possible, they should be kept away from internal circuits and metallic parts, and there should be one downconductor installed at each corner of the structure [16]. The earth-termination system consists of buried interconnected conductors. The lower the impedance to earth, the more effective the earth-termination system would be. The impedance decreases with the number of parallel earth electrodes, and the lowest is obtained when the largest surface of the electrodes is in contact with soil.

According to the IEC 62305-3 standard [START_REF]3: Physical damage to structures and life hazard[END_REF], there are two basic earthing arrangements:

Type A: Earth electrodes connected to the down conductors, without forming a loop.

Type B: Earth electrodes forming a closed ring around the structure, or a foundation earth electrode.

The type B arrangement is recommended in large facilities. It is worth noting that since the soil is an imperfect conductor, the earth ring electrodes of different buildings in a facility must be interconnected to avoid potential differences. One way of interconnecting the rings is with a meshed grounding network. This is called a meshed earth-termination system. The smaller the mesh size of the grounding network, the lower the potential differences between the buildings. An example of a meshed earth-termination system of an industrial facility is shown in figure 1234.

Internal lightning protection system

The basis of the internal lightning protection system is surge suppression and potential equalization. Potential equalization is necessary to prevent flashovers inside the structure, hazardous touch voltages, and the accumulation of a static charge.

In general, all the conductors entering the structure are connected to a protective equipotential bonding system. The bonding system is also connected to the external lightning protection system and to other conductive components of the structure, such as water pipes. Since some cables cannot be directly bonded to the earth-termination system, surge protective devices (SPDs) are required [START_REF] Cooray | An Introduction to Lightning[END_REF]. An SPD is a device that changes the impedance between its terminals to divert surge currents or unwanted signals to the ground.

Structures in industrial facilities typically have a coordinated SPD system, with SPDs installed at different levels. For example, at the entrance to the structure, at the entrance to the power system, and to connect sensitive electronic devices. According to the lightning protection zone (LPZ) concept [START_REF]Electrical and electronic systems within structures[END_REF], SPDs should be installed at the boundary of two zones. The zones are theoretically assigned parts of space with different expositions to the effects of the lightning electromagnetic pulse (LEMP).

Hence, they are defined by the protection measures employed. Basic protection measures include structural components which can be used as spatial shields, such as the reinforcement of a building.

Industrial context and background

The design and the operation of Nuclear Power Plants (NPPs) in France are based on the installation nucléaire de base (INB) regulation of February 2012 [18], wherein lightning is considered a natural external event that has to be taken into account in the safety demonstration. Even though the regulation does not impose a procedure to assess the risk associated with a lightning strike, it demands a deterministic safety demonstration. Therefore, all the methods and tools used in the analysis must be validated and their limits and accuracy well estimated.

In the European Union, all nations are encouraged to transfer the IEC 62305 standards [10] to national standards [START_REF] Heidler | Lightning standards[END_REF]. They provide valuable principles to design the LPS and the protection measures against the LEMP for equipment and structures.

The principles could be used to protect industrial facilities; however, the standard is not dedicated to NPPs and does not meet the requirements to be used in the safety demonstration. There are other standards, such as IEC 62003 [20] and KTA 2206 [START_REF]KTA 2206: Design of Nuclear Power Plants Against Damaging Effects from Lightning[END_REF], that were specifically created for NPPs. The first provides the principles for electromagnetic compatibility testing of instrumentation and control equipment, and the second specifies the safety precautions against adverse effects of lightning. Yet, experience has shown that in general, the current standards are not adapted to the characteristics of French NPPs [START_REF] Duqueroy | Lightning risk assessment evaluation on French nuclear power plants[END_REF][START_REF] Duquerroy | Evolution of lightning protection of nuclear power plants: An overview of edf's experience[END_REF].

On average, Météorage detects 453 000 ground flashes in metropolitan France every year [START_REF]Foudre, éclair et tonnerre[END_REF]. Each flash is a natural source of electromagnetic interference, which may cause physical damage and lead to an upset of electric installations and sensitive electronic devices. In industrial facilities, even a momentary malfunction can have severe consequences.

Let us consider a direct lightning strike to a reinforced concrete building in the facility. During the strike, the lightning current flows along the channel and the reinforcement of the building, generating a transient electromagnetic field. Now, if we consider the case shown in figure 12345, a part of the lightning current would dissipate into earth through the foundation of the building, another part would flow through the meshed grounding network, and the rest would be distributed among the ducts and the conductors leading away from the building. Moreover, for the safety demonstration, it would be essential to estimate, among others, the electromagnetic field inside the building and the surges in the cables running in the ducts. Because of the complexity of industrial facilities, the estimation is far from being straightforward.

Numerous studies have been carried out to deepen the understanding of the phenomenon, see e.g. [START_REF] Heidler | Induced overvoltages in cable ducts taking into account the current flow into earth[END_REF][START_REF] Zischank | Shielding effectiveness of reinforced concrete cable duct carrying partial lightning currents[END_REF][START_REF] Zischank | Assessment of the lightning transient coupling to control cables interconnecting structures in large industrial facilities and power plants[END_REF][START_REF] Diaz | Lightning transient voltages in cables of a large industrial site using a FDTD thin wire model[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]. In general, even calculating the distribution of the lightning current among the conductors can become a complicated problem. The international standard IEC 62305-1 [29] assumes that half of the current goes straight to the earth-termination system of the building and proposes a simple model to distribute the remaining half. The German standard KTA 2206 [START_REF]KTA 2206: Design of Nuclear Power Plants Against Damaging Effects from Lightning[END_REF] specifies weighting factors to calculate the portion of the current flowing through the ducts and the soilcontacting conductors, assuming that one-third of the current flows into the ground.

Also, it proposes a formula to calculate the axial voltage induced in cables running in the ducts from the partial lightning current flowing through the reinforcement of the ducts. 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].

Regarding the electromagnetic field, a very good estimation could be obtained using full-wave methods. A full-wave numerical approach allows one to consider the majority of the characteristics and components of the electromagnetic environment, see e.g. [START_REF] Tatematsu | Application of the FDTD-based simulation code VSTL REV to the lightning surge analysis of a nuclear power plant[END_REF][START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF]. The problem is that the simulations can quickly become computationally expensive. Alternatively, the IEC 62305-4 standard [START_REF]Electrical and electronic systems within structures[END_REF] suggests a simple formula to calculate the maximum value of the magnetic field strength in single-layer grid-like shields. Yet, as all formulas, it comes with a certain degree of uncertainty.

Motivation and objectives

The accuracy of the results obtained from full-wave simulations depends on the representativeness of the computation model, and when it comes to large industrial facilities, making the right choices can be challenging. It is virtually impossible to model all the elements of a full-scale industrial facility and consider all the components of the electromagnetic environment. The amount of computational resources required to carry out the simulations can increase exponentially with the level of detail. Hence, carrying out full-wave simulations to evaluate protection measures is not always feasible nor practical daily.

From all the effects a direct lightning strike to a building could have on large industrial facilities, we decided to focus on the transient magnetic field inside reinforced concrete buildings and the surges in underground cables. Nearby strikes are not considered because lightning is more likely to strike high grounded vertical objects [START_REF] Rakov | Fundamentals of Lightning[END_REF].

Nor are the surges induced by lightning strikes on overhead lines because they have been studied extensively (see e.g. [START_REF] Baba | Voltages Induced on an Overhead Wire by Lightning Strikes to a Nearby Tall Grounded Object[END_REF][START_REF] Fernando H Silveira | Voltages Induced in Single-Phase Overhead Lines by First and Subsequent Negative Lightning Strokes: Influence of the Periodically Grounded Neutral Conductor and the Ground Resistivity[END_REF][START_REF] Rachidi | A Review of Field-to-Transmission Line Coupling Models With Special Emphasis to Lightning-Induced Voltages on Overhead Lines[END_REF][START_REF] Fernando | Lightning-Induced Voltages Over Lossy Ground: The Effect of Frequency Dependence of Electrical Parameters of Soil[END_REF][START_REF] Tatematsu | Three-Dimensional FDTD Calculation of Lightning-Induced Voltages on a Multiphase Distribution Line With the Lightning Arresters and an Overhead Shielding Wire[END_REF]). Moreover, this work aims to propose reliable methods adapted to large industrial facilities, as an alternative to full-wave simulations, to estimate the peak-values of both the magnetic field and the surges. Three objectives are defined consequently:

• Conceive a traceable electromagnetic model to study the phenomena using a full-wave numerical approach.

• Identify the most influential parameters to reduce the complexity of the model.

Note that this would result in two different 3-D models because some of the parameters that influence the magnetic field inside the building can have a negligible effect on the surges.

• Build different surrogate models and efficient algorithms to create a set of tools one could use to calculate the peak-values of the magnetic field and the surges.

Research with similar objectives was conducted from 2013 to 2016 [START_REF] Luis | Lightning Induced Voltages in Cables of Power Production Centers[END_REF]. It showed the pertinence of full-wave methods to study the phenomena and that it is possible to model the primary structures of industrial facilities using the thin-wire formalism of Holland [START_REF] Holland | Finite-Difference Analysis of EMP Coupling to Thin Struts and Wires[END_REF]. The research was concluded with six polynomial models of the peak-values of lightning-induced voltages (LIVs) in instrumentation and measurement cables. Different models were necessary due to the dependence of the effect of different parameters on the type of excitation. In this work, besides studying the lightninginduced magnetic fields, we include additional components and characteristics of large industrial facilities that were not considered to compute LIVs.

Chapter 2

Towards a set of tools to calculate the peak-values of the magnetic field and the surges

As an alternative to analytical approaches, several computation methods have been developed over the years to assess electromagnetic problems. Due to the simplicity of circuit-based simulations, the Transmission Line (TL) theory is still widely used to calculate the electromagnetic fields and surges generated by a lightning strike.

Nevertheless, the need for more accurate results and the development of affordable high-performance computers have made the full-wave methods gain popularity.

Full-wave methods

Full-wave methods can be subdivided into several categories based on the formulation of Maxwell's equations and the numerical approximation used to solve them.

Typically, the lightning community uses software based on the finite-difference timedomain (FDTD) method and the method of moments (MoM) due to their capabilities to handle three-dimensional structures, as well as cables and wiring structures (see e.g. [START_REF] Tatematsu | A Technique for Calculating Voltages Induced on Twisted-Wire Pairs Using the FDTD Method[END_REF][START_REF] Baba | Electric and Magnetic Fields Predicted by Different Electromagnetic Models of the Lightning Return Stroke Versus Measured Fields[END_REF][START_REF] Fridolin | Share of the Lightning Current on Electric Installation Lines Within a Building Considering Equipotential Bonding at Roof Level[END_REF][START_REF] Metwally | Computation of Transient Overvoltages in Low-Voltage Installations During Direct Strikes to Different Lightning Protection Systems[END_REF][START_REF] Tatematsu | Lightning Current Simulation of 66-kV Substation With Power Cables Using the Three-Dimensional FDTD Method[END_REF]). Both methods have advantages and drawbacks. The choice of the right technique depends on the specificities of the environment and the com-putational resources available.

The method of moments (MoM)

The method of moments (MoM), formalized by Harrington in 1968 [START_REF] Roger F Harrington | Field computation by moment methods[END_REF], is a frequencydomain technique used to solve Maxwell's equations in the integral form. The technique projects the Maxwell equations into a system of linear equations that can be solved numerically, by enforcing boundary conditions. As a frequency-based technique, it does not suffer from numerical dispersion.

The unknowns of the system are the current and the charge densities induced on the surface of the structures. Since only the surfaces are discretized, the dimensions of the problem are reduced by one. That is the main advantage of the MoM over other full-wave methods. However, the number of unknowns can quickly grow into the thousands for a complex geometry, which can be problematic given that the iterative matrix formulation is used to solve the system. The bigger the number of unknowns, the higher the memory requirements to store the matrix, and the longer the time required to find a solution. Nevertheless, it is worth noting that methods such as the fast multipole method (FMM) [START_REF] Coifman | The fast multipole method for the wave equation: a pedestrian prescription[END_REF] can be applied to accelerate the iterative solver. In addition, the computation and storage cost can be reduced by implementing matrix compression methods, such as hierarchal matrices (H-matrix) [START_REF] Negi | Memory reduced half hierarchal matrix (h-matrix) for electrodynamic electric field integral equation[END_REF].

Also, the MoM, like all integral methods, has difficulties handling inhomogeneous materials, and the electric field integral equation (EFIE) it typically uses is not adapted for low-frequency applications. The method is more suitable for narrowband problems and frequency sweep responses than for problems that are of transient nature. Since lightning processes are time-dependent and radiate over a broad range of frequencies, a finite-difference (FD) method solving Maxwell's equations in the time domain could be more appropriate to study the phenomena.

The finite-difference time-domain (FDTD) method

The finite-difference time-domain (FDTD) method, introduced by Yee in 1966 [START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF], is one of the most commonly used time-domain methods. It solves Maxwell's equations in the partial differential form, by discretizing the spatial domain and the time interval. The FDTD method is highly versatile and accurate because the equations are solved without making any physical approximation or assumption about the system.

However, like for all systems of partial differential equations, the initial conditions and the boundary conditions of the system must be specified. Thus, to assume an open space, absorbing boundaries must be defined. To do so, one can implement absorbing boundary conditions (ABCs) or perfectly matched layers (PMLs).

There are several formulations of ABCs depending on the materials at the boundaries and the angle of incidence of radiating waves (see e.g. [START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF][START_REF] Lindman | Free-space" boundary conditions for the time dependent wave equation[END_REF][START_REF] Mur | Absorbing Boundary Conditions for the Finite-Difference Approximation of the Time-Domain Electromagnetic-Field Equations[END_REF]). Yet, existing

ABCs restrict themselves to absorbing waves at certain angles [START_REF] Steven | Notes on Perfectly Matched Layers (PMLs)[END_REF]. As a solution,

Berenger introduced the concept of matched absorbing layers in 1994 [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. An absorbing layer is an absorbing material that is placed adjacent to the edges of the computational space to absorb outgoing waves [START_REF] Steven | Notes on Perfectly Matched Layers (PMLs)[END_REF]. Moreover, a PML is an absorbing medium for radiating waves, theoretically independent of the incident angle and effective over a large frequency band. It is constructed to prevent reflections at the interface between materials. Implementing PMLs in FDTD simulations is straightforward, but it comes with an additional computational cost.

Regular PMLs have difficulties absorbing evanescent waves, which can be dominant below a certain cutoff frequency. The cutoff frequency depends on the dimensions of the structures considered in the model. To overcome this limitation, a complex frequency shifted (CFS) tensor coefficient can be included in the formulation of the PMLs. This new formulation is known as CFS-PML (CMPL) and was introduced by Kuzuoglu and Mittra in 1996 [START_REF] Kuzuoglu | Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers[END_REF].

If to the additional computational costs of PMLs, we add the fact that in FDTD simulations the electromagnetic fields are computed at every point in space, it is clear that the FDTD method might not be the most efficient. Nevertheless, it is advantageous when it comes to modeling industrial facilities, mainly because of the large number of elements that must be taken into account. In FDTD simulations, as long as the spatial domain remains constant, the number of elements in the model does not have an important effect on the calculation time. In addition, subcellular models, such as thin-wires [START_REF] Holland | Finite-difference analysis of emp coupling to thin struts and wires[END_REF][START_REF] Edelvik | A new technique for accurate and stable modeling of arbitrarily oriented thin wires in the fdtd method[END_REF] and thin-slots [START_REF] Edelvik | Stable modeling of arbitrarily oriented thin slots in the fdtd method[END_REF], can be implemented to optimize the space discretization when handling multiscale elements. For example, the effect of a wire the transversal area of which is small compared to the other elements of the model can be taken into account without reducing the size of the FDTD cell, by hybridizing Maxwell's equations and 1-D transmission line equations.

Another important advantage of the FDTD method is that broadband results can be obtained with a single simulation, and also, the incorporation of nonlinear materials is fairly simple. Hence, we have decided to carry out the simulations using transient electromagnetic simulator -finite difference (TEMSI-FD) [START_REF]TEMSI-FD: Time Electromagnetic Simulator -Finite Difference Time Domain[END_REF], an FDTDbased solver developed by the XLIM institute.

TEMSI-FD

The XLIM Institute has been developing and updating TEMSI-FD since 2002. It includes various techniques to represent thin-wires with oblique trajectories [START_REF] Guiffaut | Cartesian shift thin wire formalism in the fdtd method with multiwire junctions[END_REF][START_REF] Guiffaut | New oblique thin wire formalism in the FDTD method with multiwire junctions[END_REF],

insulated wires [START_REF] Guiffaut | Insulated oblique thin wire formalism in the FDTD method[END_REF], and coaxial cables [START_REF] Guiffaut | Des fils obliques pour une modélisation conforme et sans maillage des câbles dans la méthode FDTD. Bilan et extensions[END_REF][START_REF] Diaz | Lightning transient voltages in cables of a large industrial site using a FDTD thin wire model[END_REF], and can be used in parallel computations to reduce the calculation time. In addition, since TEMSI-FD is developed internally, it comes with additional benefits: fixed and clear assumptions behind the calculations, the possibility to add new features, and full technical support.

In the last five years, TEMSI-FD has been used to calculate, among others, LIVs in cables [START_REF] Diaz | Lightning transient voltages in cables of a large industrial site using a FDTD thin wire model[END_REF][START_REF] Luis | Lightning Induced Voltages in Cables of Power Production Centers[END_REF][START_REF] Bdour | Emulation of an Expensive FDTD Code With Mixed Quantitative and Qualitative Factors for Analysis of Lightning-Induced Transient Responses[END_REF], lightning effects on aircraft [START_REF] Monferran | Modèles numérique et stochastique des fixations pour la contrainte foudre des lignes d'assemblages sur avion[END_REF][START_REF] Monferran | Lightning currents on fastening assemblies of an aircraft fuel tank, part ii: Fdtd modeling merged with a circuit model supplemented by a statistical model[END_REF][START_REF] Bastard | Modélisation électromagnétique unifiée pour les contraintes foudre et champ fort sur aéronefs[END_REF], and electromagnetic fields inside reinforced concrete buildings [START_REF] Laisné | Eléments de modélisation de la foudre sur les bâtimentspartie 2 : simulation numérique 3d d'un bâtiment soumis à une onde plane[END_REF][START_REF] Laisné | Eeléments de modélisation de la foudre sur les bâtimentspartie 3 : simulation numérique 3d d'un bâtiment soumis à un impact direct et à une boucle d'injection[END_REF].

Computation models

Once again, the accuracy of the results obtained using any full-wave method depends on the representativeness of the computation model, and unfortunately, there will always be uncertainties. First, modeling every single detail of the electromagnetic environment of an industrial facility is unfeasible. Second, our knowledge of physical phenomena, such as lightning, is limited.

To calculate the electromagnetic fields and surges generated by a lightning strike in an industrial facility, we focus on modeling three main components of the electromagnetic environment: the lightning return stroke, the soil, and reinforced concrete structures.

Lightning return stroke

The electromagnetic field generated by the return stroke depends on the variation of the current along the lightning channel. Until now, this variation has been inferred from the observation of the optical radiation of the channel [START_REF] Cooray | An Introduction to Lightning[END_REF], which in addition to the current and the velocity of the return stroke measured at the channel base, has resulted in different models.

The return stroke models have been categorized into four classes, based on their governing equations [START_REF] Rakov | Review and evaluation of lightning return stroke models including some aspects of their application[END_REF]: gas dynamic models, electromagnetic (EM) or antennatheory models, distributed circuit models, also known as transmission line models, and engineering models. The latter are the most popular in engineering applications.

They can easily be implemented as a set of equations relating the temporal-spatial variation of the current along the channel to the channel base current and the return stroke speed. To this class belong the well-known transmission line (TL) model introduced by Uman and McLain [START_REF] Uman | Magnetic field of lightning return stroke[END_REF], not to be confused with the class mentioned before, the modified transmission line model with linear decay with height (MTLL) [START_REF] Rakov | A modified transmission line model for lightning return stroke field calculations[END_REF], and the modified transmission line model with exponential decay with height (MTLE) [START_REF] Nucci | On lightning return stroke models for LEMP calculations[END_REF].

However, since we chose a full-wave method to carry out the simulations, we will implement an EM model, which is also based on Maxwell's equations. Some advantages of EM models include rigor because a transverse electromagnetic (TEM) mode is not assumed, and the capability of reproducing both the features of the lightning current distribution and the resulting electromagnetic fields [START_REF] Baba | Electromagnetic models of the lightning return stroke[END_REF]. The lightning channel will therefore be represented as a monopole antenna by a vertical thin-wire [START_REF] Guiffaut | New oblique thin wire formalism in the FDTD method with multiwire junctions[END_REF].

The wire is excited at its base by a lumped current source and connected at the top end to a PML.

Following the recommendation given in [START_REF] Karami | On Practical Implementation of Electromagnetic Models of Lightning Return-Strokes[END_REF], instead of embedding the wire in a fictitious half-space dielectric medium or a dielectric coating, the wire is loaded by distributed inductance and resistance to reproduce the features observed using optical techniques. We consider a distributed series resistance of 1 Ω/m to control the current attenuation and match the impedance of the channel, which is expected to be between 0.6 and 2.5 kΩ [START_REF] Gorin | Measurements of lightning currents at the Ostankino tower (in Russian)[END_REF]. In addition, we consider distributed series inductance of 4 𝜇H/m to adjust the propagation speed of the current along the channel. Assuming that the speed of the current is equal to the optical speed measured at the channel base, it is expected to be within one-third and one-half of the speed of light. The radius of the lightning channel is set to 1 cm as in the parametric study presented in [START_REF] Luis | Lightning Induced Voltages in Cables of Power Production Centers[END_REF].

Current parameters

The effects of lightning depend on the current parameters. Generally, lightning protection standards consider four parameters: the current peak, the charge, the specific energy, and the current steepness. The charge and the specific energy are responsible for mechanical and thermal effects. The current peak and the current steepness are responsible for induced effects, such as voltage drops and surges [START_REF] Heidler | Parameters of lightning current given in iec 62305 -background, experience and outlook[END_REF]. The first stroke has the highest current peak, and the subsequent strokes have the highest current steepness. The current of the first stroke typically rises to a peak-value of 30 kA in a few microseconds, whereas the current of the subsequent strokes rises to a peak-value of 10 kA in less than a microsecond [4]. The current waveform can be analyzed considering the commonly used doubleexponential function given by the following formula:

𝑖 = 𝑖 𝑚𝑎𝑥 𝑘 ( 𝑒 𝑡/𝑡 1 -𝑒 𝑡/𝑡 2 ) , (2.1)
where 𝑖 𝑚𝑎𝑥 is the maximum value of the current, 𝑘 a correction factor, and 𝑡 1 and 𝑡 2 are once again the front time and the decay time, respectively. However, the waveform of the double exponential function cannot reproduce the features of subsequent strokes.

With the parameters defined for subsequent strokes, the current would unrealistically start with the maximum current steepness at 𝑡 = 0 s [START_REF] Heidler | Parameters of lightning current given in iec 62305 -background, experience and outlook[END_REF].

Hence, to adjust the values the IEC 62305 standard suggests the Heidler function given by the following formula:

𝑖 = 𝑖 𝑚𝑎𝑥 𝑘 × (𝑡/𝑡 1 ) 10 1 + (𝑡/𝑡 1 ) 10 × 𝑒 -𝑡/𝑡 2 , (2.2)
where 𝑖 𝑚𝑎𝑥 is still the maximum value of the current, 𝑘 a correction factor, 𝑡 1 the front time, and 𝑡 2 the decay time. This is the function considered in this work.

In addition, to reduce the calculation time when analyzing the results of the simulations in the frequency domain, we consider a Gaussian function

𝑔(𝑡) = 𝐴 0 𝑒 ( - (𝑡-𝑡 0 ) 2 𝑇 2 ) FT -→ 𝑒(𝑓 ) = 𝐴 0 𝑇 √ 𝜋𝑒 -(𝜋𝑓 𝑇 ) 2 𝑒 -𝑗2𝜋𝑓 𝑡 0 , (2.3)
where 𝐴 0 is the maximum amplitude of the pulse, 𝑒 is the exponential function, 𝑡 0 is the position of the center of the pulse, and 𝑇 is the width of the pulse at 77 % of the maximum. The frequency range covered by the pulse depends on the latter.

The width of the pulse at 77 % of the maximum (𝑇 ) is defined by

𝑇 = √ 𝑙𝑛(𝐴𝑡𝑡 𝑓𝑚𝑎𝑥 ) 𝜋𝑓 𝑚𝑎𝑥 , (2.4)
where 𝑓 𝑚𝑎𝑥 is the upper limit of the frequency rage covered by the pulse, 𝐴𝑡𝑡 𝑓𝑚𝑎𝑥 is the attenuation at 𝑓 𝑚𝑎𝑥 with respect to the magnitude spectrum at 𝑓 = 0 Hz, and 𝑙𝑛 is the natural logarithm.

The position of the center of the pulse (𝑡 0 ) is defined by

𝑡 0 = 𝑇 √ 𝑙𝑛(𝐴𝑡𝑡 0 ), (2.5)
where 𝐴𝑡𝑡 0 is the attenuation of the signal at 𝑡 = 0 s, and 𝑙𝑛 is once again the natural logarithm.

To cover a frequency range from 0 Hz to 10 MHz, 𝑡 0 and 𝑇 are set to 0.179 𝜇𝑠 and 48.3 𝑛𝑠, respectively. In addition, 𝐴𝑡𝑡 0 is set to 10 6 to avoid any noise in the time-domain response.

When considering a Gaussian function, we apply a Fourier transform (FT) to the results and normalize them to the source, i.e., we divide the results in the frequency domain by the spectrum of the Gaussian pulse. Thus, the results are independent of the current waveform; they become transfer functions. Moreover, instead of carrying out a simulation for each of the different strokes, we can carry out a single simulation using a Gaussian pulse and then compute the results for any other current waveform.

It is as simple as multiplying the spectrum of the waveform of interest by the transfer function and then computing the inverse Fourier transform. Note, however, that the complete convergence of the response is necessary to apply a Fourier transform, and the Gaussian pulse used in the simulation must cover the frequency range of the waveform of interest.

Matrix Pencil method MPM

Gaussian pulses covering the frequency ranges of the lightning return stroke have a duration of just a few microseconds. Nevertheless, it could happen that even after 50 𝜇𝑠, the response has not yet converged. This could be explained by limited lowfrequency spectral energy. If in addition to the inductive effect, there are low resistive losses, it would be difficult to achieve convergence. Also, the presence of resonances can delay energy dissipation.

Due to the stability criteria of the FDTD method, the FDTD time step oversamples the wavelength of the source. Thus, both waiting for the response to converge and considering the Heidler function lead to expensive simulations. Instead of increasing the simulation time to achieve convergence, one can extrapolate the results.

Several extrapolation methods are available in the literature (see e.g. [START_REF] Sidi | Practical Extrapolation Methods: Theory and Applications[END_REF]). Among the most widely used methods in electromagnetic problems, one can find the vectorfitting (VF) method [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF] and the matrix pencil method (MPM) [START_REF] Sarkar | Using the matrix pencil method to estimate the parameters of a sum of complex exponentials[END_REF]. The MPM consists in 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. As opposed to the VF method, the MPM does not require starting poles. The method has already been used for high-intensity radiated field (HIRF) applications [START_REF] Bastard | Modélisation électromagnétique unifiée pour les contraintes foudre et champ fort sur aéronefs[END_REF] to improve the efficiency of the FDTD method at low frequencies [START_REF] Perrin | Using transfer function calculation and extrapolation to improve the efficiency of the finite-difference time-domain method at low frequencies[END_REF]. Also, it has been applied to infer a high-frequency model of a power transformer and fit the frequency responses of a cable system [START_REF] Sheshyekani | Application of the Matrix Pencil Method to Rational Fitting of Frequency-Domain Responses[END_REF].

The main advantage of the MPM is the reduction in the computation time. Yet, even with a complete time convergence of the response, it could facilitate signal processing. The method gives the option to set the frequency step, which can be useful when multiplying the normalized results by the spectrum of the waveform of interest and applying the inverse Fourier transform. It is worth noting that one can also adjust the frequency step of the results using a moving average filter or applying the VF method.

Soil

The soil constitutes an important part of the electromagnetic environment, and it is one of the components that come with most of the uncertainty. The electrical parameters of the soil can experience variations of several orders of magnitude with small changes in the temperature, the moisture content, the salt concentration, and so forth.

Typically, the soil is assumed homogeneous and characterized by a dielectric constant 𝜖 and a low-frequency resistivity 𝜌, which depends on the type of soil. However, these parameters are known to be frequency-dependent [START_REF] Fernando | Lightning-Induced Voltages Over Lossy Ground: The Effect of Frequency Dependence of Electrical Parameters of Soil[END_REF]. In addition, the lightning radiated electromagnetic field cannot be fully reproduced under the assumption of a homogeneous ground [START_REF] Li | Importance of Taking Into Account the Soil Stratification in Reproducing the Late-Time Features of Distant Fields Radiated by Lightning[END_REF]. It has been shown that the frequency-dependent parameters reduce the amplitude of induced voltages in overhead lines [START_REF] Fernando | Lightning-Induced Voltages Over Lossy Ground: The Effect of Frequency Dependence of Electrical Parameters of Soil[END_REF], influence the performance of the grounding grid [START_REF] Visacro | Frequency Dependence of Soil Parameters: Experimental Results, Predicting Formula and Influence on the Lightning Response of Grounding Electrodes[END_REF], and affects the early time response of the electromagnetic field, especially in a poorly-conducting ground [START_REF] Li | On the influence of the soil stratification and frequencydependent parameters on lightning electromagnetic fields[END_REF]. Nevertheless, the soil parameters are considered frequency-independent in this work due to their minor influence on the variables of interest compared to the influence of other parameters.

Disregarding this dependence is similar to neglecting the soil ionization, which is not considered either. The soil ionization increases the conductivity, and therefore, it also reduces the adverse effects of lightning. Hence, in both cases, we remain on the safe side.

Regarding the electromagnetic fields, even when considering the soil as a perfectly conducting plane, the vertical electric field and the azimuthal magnetic field at close range can be well approximated. Yet, a horizontally stratified model is necessary to reproduce the horizontal electric field [START_REF] Li | On the influence of the soil stratification and frequencydependent parameters on lightning electromagnetic fields[END_REF]. While the early response of the field can be reproduced by adjusting the resistivity of a homogeneous ground, the late-time response of the field is very sensitive to ground stratification. A good agreement can be found using a two-layer soil [START_REF] Li | Importance of Taking Into Account the Soil Stratification in Reproducing the Late-Time Features of Distant Fields Radiated by Lightning[END_REF]. However, once again because of the little influence on the variables of interest, we consider a homogeneous soil in this work.

Reinforced concrete structures

In industrial facilities, the buildings and the cable ducts are generally made of reinforced concrete, one of the most widely used materials in modern constructions. As its name indicates, it is made of concrete and reinforcement. Concrete is a composite material that is strong under compression but less so in tension, which makes it on its own unsuitable for structures exposed to tensile stress. To increase the concrete's tensile strength, reinforcement is therefore included.

As figure 2-2 illustrates, the reinforcement usually consists of at least two interconnected layers of reinforcing grids, where cylindrical reinforcing steel bars (rebars) are arranged periodically to form square or rectangular meshes. Rebars are typically welded or wire-tied to ensure that they deliver the expected tensile strength. It has been demonstrated that the grid-like shield formed by the rebars increases the through-wall attenuation of signals [START_REF] Casey | Electromagnetic shielding behavior of wire-mesh screens[END_REF][START_REF] Dalke | Effects of reinforced concrete structures on RF communications[END_REF][START_REF] Richalot | Electromagnetic propagation into reinforced-concrete walls[END_REF][START_REF] Pena | Measurement and modeling of propagation losses in brick and concrete walls for the 900-MHz band[END_REF] and enhances the shielding effectiveness of the structure against both high-altitude electromagnetic pulse [START_REF] Hyun | Analysis of shielding effectiveness of reinforced concrete against high-altitude electromagnetic pulse[END_REF] and lightning electromagnetic fields [START_REF] Zischank | Laboratory simulation of direct lightning strokes to a modeled building: measurement of magnetic fields and induced voltages[END_REF][START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF][START_REF] Kern | Magnetic fields and induced voltages in case of a direct strike-comparison of results obtained from measurements at a scaled building to those of IEC 62305-4[END_REF][START_REF] Metwally | Reduction of lightning-induced magnetic fields and voltages inside struck double-layer grid-like shields[END_REF][START_REF] Mardiguian | Lightning-generated fields in reinforced concrete buildings[END_REF][START_REF] Maksimowicz | Investigation of models of grid-like shields subjected to lightning electromagnetic field: experiments in the frequency domain[END_REF][START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF].

The layers are interconnected by vertical and horizontal rebars known as hoops.

In this work, we consider the hoops to be positioned at the intersection of the rebars, their radius to be equal to the radius of the rebars, and the meshes of the layers to be aligned. We decided to assume that all the rebars forming the reinforcement of the same structure are identical and there is a perfect electrical connection between them at every intersection. Also, that the properties of the materials are not affected by external factors. Moreover, we are disregarding, among others, the contact resistance, thermal effects, and non-linearities. Further studies are required to understand the implications of this decision and calculate the uncertainty associated.

Since it is necessary to represent the reinforcement as realistic as possible to obtain accurate results, the rebars are modeled as thin-wires [START_REF] Guiffaut | New oblique thin wire formalism in the FDTD method with multiwire junctions[END_REF]. We thought about using thin slots; however, we realized that it could lead to an important error, especially when calculating the electromagnetic field inside the building. There are two main reasons for that: first, the lightning current does not distribute uniformly in the reinforcement, and therefore, a compensation of the field may not occur exactly in the same areas. Second, because the electromagnetic behavior of reinforced concrete changes when the reinforcement is characterized in transmission or conduction [START_REF] Mardiguian | Lightning-generated fields in reinforced concrete buildings[END_REF][START_REF] Maksimowicz | Investigation of models of grid-like shields subjected to lightning electromagnetic field: experiments in the frequency domain[END_REF][START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF]. Clearly, a thin-slot model would not be able to reproduce these features.

If the conductivity of the slot is adjusted to match the impedance of the reinforcing grid, it would not match the shielding effectiveness. Even if the average does, it would not be able to reproduce the electromagnetic field locally or close to the walls.

Surrogate modeling

The first objective of this work is to conceive a 3-D model of industrial facilities for both understanding the effects of a direct lightning strike to a building and evaluating the magnetic field and the surges. By modeling an industrial facility, we are considering it as a set of components forming a system, which is delimited by spatial and temporal boundaries and influenced by the environment. Based on the initial conditions and the inputs, the model reflects the behavior of the original system with a certain degree of uncertainty. Since carrying out simulations using a full-wave approach can be computationally expensive, we will consider the response of the 3-D model, i.e. the simulation results obtained using TEMSI-FD, as the response of the original system and use them to build more efficient computational models.

A computational model maps the inputs of the system into the outputs, but to do so it has to be calibrated. The process of calibrating consists basically of adjusting the coefficients of the model to fit the available data. After conceiving the 3-D model, the first thing we need to do is, therefore, generate that data.

Design of experiments (DoE) and sampling techniques

One could think that a model would make better predictions when trained on a large dataset. However, more than the quantity of data, what matters is the position of the samples in the design space. If the samples are locally clustered, the model might not be able to make good predictions beyond that specific region.

Once the design space has been defined, several sampling strategies are available to generate samples. Note that the design space is bounded by the upper and lower limits of the input variables being considered. One of the most popular sampling strategies is Monte Carlo (MC) sampling [START_REF] Fishman | Monte Carlo[END_REF], although it can be inefficient because the MC methods rely on repeated random sampling. Besides, they usually require an important number of sample points to converge. The MC methods are widely used to estimate the probability density function (PDF) of a target function, approximate a quantity numerically, and optimize functions under uncertainty. Alternatively, to minimize the computational cost of running expensive simulations, the Design of experiments (DoE) and statistical approximation techniques are often used [START_REF] Simpson | Sampling Strategies for Computer Experiments: Design and Analysis[END_REF].

The DoE theory [START_REF] Montgomery | Design and Analysis of Experiments[END_REF] was initially developed to plan physical experiments. Since there is a random error associated with physical experiments, many of the "classical" experimental designs replicate sample points. In the "classical" designs, a few sample points are left in the middle and the rest are spread around the boundaries of the design space. To this category belong, among others, the well-known factorial and fractional factorial designs, the central composite design, and the Box-Behnken design. The central composite and Box-Behnken designs are also known as response surface designs [START_REF] Montgomery | Response Surface Methods and Designs[END_REF].

In 1989, Sack et al. [START_REF] Sacks | Design and Analysis of Computer Experiments[END_REF] stated that when it comes to deterministic computer experiments, the sample points should be chosen to fill the design space. Moreover, scientist started using "space filling" designs [START_REF] Simpson | Sampling Strategies for Computer Experiments: Design and Analysis[END_REF], such as Latin hypercube designs [START_REF] Park | Optimal Latin-hypercube designs for computer experiments[END_REF]. In Latin hypercube designs, the PDF is divided into disjoint intervals with equal probability, and a sample is randomly drawn from each interval. Thus, the samples are better distributed in the design space. Nevertheless, because of the complexity of industrial facilities, we have decided to focus on "classical" designs.

Factorial designs [START_REF] Montgomery | The 2k Factorial Design[END_REF][START_REF] Montgomery | Introduction to Factorial Designs[END_REF] Experiments can be used to study the behavior of a system. In the experiments, the system is considered as a process that transforms some input into an output. The process involves several factors, and the main objective of the experiments is therefore to determine the influence of those factors on the output response. Hence, given a system, one must choose the factors to be considered and define the ranges over which the factors will be varied. In addition, in factorial designs, the number of levels for each factor must be defined. A level is a value taken by the factor within its range of variation, at which runs will be made. The appropriate number of levels depends on the nature of the system. Thus, non-statistical knowledge of the phenomena is essential to design experiments.

The number of runs in a factorial design is given by all the possible combinations of the levels of the factors considered. This is why they are also known as full factorial designs. Moreover, if we had factors 𝐴, 𝐵, and 𝐶, with two levels, three levels, and four levels, respectively, there would be 2 × 3 × 4 possible combinations → 24 runs.

The variation in the response produced by a change in the level of a factor defines its effect.

One of the most popular factorial designs is that of 𝑘 factors, each at two levels. It is called the 2 𝑘 factorial design, and it is widely used in factor screening experiments because it requires the smallest number of runs with 𝑘 factors. The levels, denoted "-"

and "+", correspond to the lower and the upper limits of the range of variation of each factor, respectively. The high level of any factor is also denoted by the corresponding lowercase letter, and by convention, (1) is used to denote all factors at the low level.

For example, in an experiment with three factors, we would require 2 3 runs. As shown in figure 2-3, the eight combinations of the levels can be displayed geometrically as a cube. The design matrix of the experiment is shown in table 2.2. It is worth noting that 2 𝑘 factorial designs are not always ideal because it is not possible to detect a curvature effect with two levels. If the relationship between

2 𝑎 + - - 3 𝑏 - + - 4 𝑎𝑏 + + - 5 𝑐 - - + 6 𝑎𝑐 + - + 7 𝑏𝑐 - + + 8 𝑎𝑏𝑐 + + +
the response and a quantitative factor is not linear, at least three levels must be considered. Also, if there is a qualitative factor, the levels should match the qualities.

- - - + + + - 2 + - - - - + + 3 - + - - + - + 4 + + - + - - - 5 - - + + - - + 6 + - + - + - - 7 - + + - - + - 8 + + + + + + +
design required to estimate the effects specified by the experimenter.

To illustrate, let us consider a 2 𝑘 fractional factorial design containing 2 𝑘-𝑝 runs,

where 𝑝 is the number of independent generators. This design is called a 2 𝑘-𝑝 fractional factorial design. If 𝑝 = 1, it is also known as a one-half fraction of the 2 𝑘 design. Considering again an experiment with three factors, we would need to reduce the number of runs from 2 3 to 2 3-1 . Hence, we choose the highest-order interaction (𝐴𝐵𝐶) as the generator, and we determine the one-half fractions by grouping together the runs where the 𝐴𝐵𝐶 interaction has the same sign. Supposing that we select the rows with a "+" for the 𝐴𝐵𝐶 interaction as our one-half fraction, we would have the main effects 𝐴, 𝐵, and 𝐶 aliased with 𝐵𝐶, 𝐴𝐶, and 𝐴𝐵, respectively. The rows are highlighted in the extended design matrix of a 2 𝑘 full factorial design shown in table 2.3. Note that the interactions included in the matrix are simply the product of the levels of the factors interacting. Because of the aliasing, we could not estimate the main effects with the 2 3-1 design, we would be estimating 𝐴 + 𝐵𝐶, 𝐵 + 𝐴𝐶, and

𝐶 + 𝐴𝐵.

Central composite design (CCD) [98]

The central composite design (CCD) is a popular design for fitting second-order models. It consists of a 2 𝑘 factorial design with center points, that is augmented with a group of 2𝑘 axial or star points. The axial points allow the quadratic terms to be incorporated into the model. Besides the number of factors (𝑘), two parameters must be specified in the design: the distance from the center of the design space to the axial points (𝛼), and the number of center points (𝑛 𝑓 ). Generally, three to five center points are recommended.

Box-Behnken design [98]

A Box-Behnken design is formed by combining a full factorial design and an incomplete block design. The levels of the factors are positioned in the center and at the midpoints of the edges of the region defined by the upper and the lower limits of each variable. Moreover, the design requires three levels of each factor.

The number of blocks and the number of factors varying in each block depend on the number of variables involved in the design. In an experiment with three factors,

Models and algorithms

Once the data has been generated, we can start building a computational model ℳ to predict the quantity of interest 𝑌 . Let 𝑌 = ℳ(𝑋) be an output random variable, and 𝑋 = [𝑥 1 , . . . , 𝑥 𝑀 ] an input random vector.

Note that the choice of the surrogate models and supervised learning algorithms we can use will be limited by the size of the dataset. The dataset must be large enough to compute the unknown coefficients of the model and validate it. However, there is always a risk of overfitting with large datasets, which occurs when the model fits exactly the training dataset but it is unable to make good predictions for unseen data.

Linear regression [START_REF] Montgomery | The 2k Factorial Design[END_REF] One of the simplest models to implement, if not the simplest, is a linear regression model. The model assumes that there is a linear relationship between the output and input variables. The basic formulation is as follows:

𝑌 = 𝛽 0 + 𝛽 𝑇 𝑋 + 𝜀, (2.6)
where 𝛽 0 is the intercept term, 𝛽 = [𝛽 1 , . . . , 𝛽 𝑀 ] is the vector of regression coefficients, and 𝜀 is the error of the model. This is a first-order model that contains only the main effects on the output response. If the effects are not independent, the model can be improved by including the interaction terms.

The intercept term and all coefficients associated with both the variables and the interaction terms can easily be calculated from a 2 𝑘 factorial design. For example, if we consider the results under the combinations shown in table 2.2, we may estimate the factor effects as follows:

𝐴 = 1 4𝑛 [𝑎 + 𝑎𝑏 + 𝑎𝑐 + 𝑎𝑏𝑐 -(1) -𝑏 -𝑐 -𝑏𝑐], (2.7
)

𝐵 = 1 4𝑛 [𝑏 + 𝑎𝑏 + 𝑏𝑐 + 𝑎𝑏𝑐 -(1) -𝑎 -𝑐 -𝑎𝑐], (2.8
)

𝐶 = 1 4𝑛 [𝑐 + 𝑎𝑐 + 𝑏𝑐 + 𝑎𝑏𝑐 -(1) -𝑎 -𝑏 -𝑎𝑏], (2.9
)

𝐴𝐵 = (1) -𝑎 -𝑏 + 𝑎𝑏 + 𝑐 -𝑎𝑐 -𝑏𝑐 + 𝑎𝑏𝑐 4𝑛 , (2.10) 𝐴𝐶 = (1) -𝑎 + 𝑏 -𝑎𝑏 -𝑐 + 𝑎𝑐 -𝑏𝑐 + 𝑎𝑏𝑐 4𝑛 , (2.11
)

𝐵𝐶 = (1) + 𝑎 -𝑏 -𝑎𝑏 -𝑐 -𝑎𝑐 + 𝑏𝑐 + 𝑎𝑏𝑐 4𝑛 , (2.12
)

𝐴𝐵𝐶 = -(1) + 𝑎 + 𝑏 -𝑎𝑏 + 𝑐 -𝑎𝑐 -𝑏𝑐 + 𝑎𝑏𝑐 4𝑛 , (2.13)
where 𝑛 is the number of times the design is replicated. Experimental designs are generally replicated in physical experiments. Moreover, 𝑛 = 1 because for a given combination, our simulations always result in the same value.

Accordingly, the regression model reads

Ŷ = β0 + β1 𝑥 1 + β2 𝑥 2 + β3 𝑥 3 + β12 𝑥 1 𝑥 2 + β13 𝑥 1 𝑥 3 + β23 𝑥 2 𝑥 3 + β123 𝑥 1 𝑥 2 𝑥 3 = (1) 2 + 𝐴 2 𝑥 1 + 𝐵 2 𝑥 2 + 𝐶 2 𝑥 3 + 𝐴𝐵 2 𝑥 1 𝑥 2 + 𝐴𝐶 2 𝑥 1 𝑥 3 + 𝐵𝐶 2 𝑥 2 𝑥 3 + 𝐴𝐵𝐶 2 𝑥 1 𝑥 2 𝑥 3 , (2.14)
where 𝑥 𝑖 , for 𝑖 = 1, . . . , 𝑀 , is a coded variable that represents a factor. The coding brings all the variables to the same range and produces the familiar ± notation. Thus, 𝑥 𝑖 = -1 when the factor is at the lower limit of its range of variation and, 𝑥 𝑖 = 1 when it is at the upper limit. The effects are divided by two because they are estimated based on a two-unit change (from -1 to 1), whereas the coefficients measure the effect of a unit change in 𝑋 on the mean of 𝑌 .

Polynomial chaos expansion (PCE) [START_REF] Soize | Physical systems with random uncertainties: Chaos representations with arbitrary probability measure[END_REF][START_REF] Marelli | UQLab user manual -Polynomial chaos expansions[END_REF] Another relatively simple model to implement, yet more sophisticated and reliable, is a polynomial chaos expansion (PCE), which is popular in the field of uncertainty quantification. To build a PCE, it is assumed that the output response 𝑌 belongs to the Hilbert space of second-moment random variables, and it can be represented by an infinite series expansion as follows:

𝑌 = ∑ 𝛼∈N 𝑀 𝑦 𝛼 × Ψ 𝛼 (𝑋), (2.15)
where 𝑦 𝛼 are coefficients to be computed, Ψ 𝛼 (𝑋) is a basis of multivariate orthogonal polynomials, and 𝛼 = {𝛼 1 , . . . , 𝛼 𝑀 } is the multi-indices of degree |𝛼| = ∑ 𝑀 𝑖=1 𝛼 𝑖 , which indicates the degree of the multivariate polynomial.

A multivariate polynomial of the basis reads

Ψ 𝛼 (𝑋) = 𝑀 ∏ 𝑖=1 Ψ (𝑖) 𝛼 𝑖 (𝑥 𝑖 ), (2.16)
where

Ψ (𝑖)
𝛼 𝑖 (𝑥 𝑖 ) is a univariate polynomial of degree 𝛼 𝑖 , from the family of polynomials associated with the variable 𝑥 𝑖 .

The family of orthogonal polynomials used to build up the base depends on the input PDF, for example, Legendre polynomials are used for input variables with a uniform distribution. Table 2.4 shows the classical families of orthogonal polynomials and the type of variable they are associated with. Note that classical orthogonal polynomials are defined for reduced variables, e.g. uniform variables 𝒰(-1, 1) and normal variables 𝒩 (0, 1). Thus, an isoprobabilistic transform 𝑋 = 𝒯 (𝜉) is necessary.

Gaussian

Hermite 𝐻𝑒 𝛼 (𝑥)

𝐻𝑒𝛼(𝑥) √ 𝛼! Gamma Laguerre 𝐿 𝑎 𝛼 (𝑥) 𝐿 𝑎 𝛼 (𝑥) √ Γ(𝛼+𝑎+1) 𝛼! Beta Jacobi 𝐽 𝑎,𝑏 𝛼 (𝑥) 𝐽 𝑎,𝑏 𝛼 (𝑥) J 𝑎,𝑏,𝛼
In practice, the infinite series expansion is truncated to limit the number of coefficients to be computed. The number of multivariate polynomials forming the basis depends on the truncation scheme and the number of input variables 𝑀 . If we consider a standard truncation scheme, the basis would be formed by all multivariate polynomials of total degree |𝛼| less than or equal to 𝑝, and there would be (𝑀 +𝑝)! 𝑀 !×𝑝! coefficients to be computed. Nevertheless, when considering a standard truncation scheme, some coefficients turn out to be negligible. Thus, one can also reduce the number of coefficients to be computed by using sparse regression techniques [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF].

These techniques find the relevant terms of the basis from a candidate set of polynomials defined by the truncation scheme. The resultant PCE is called a sparse PCE.

As an example, let us consider a system with two input variables: 𝑥 1 ∼ 𝑈 (𝑎, 𝑏) 𝑥 2 ∼ 𝑁 (𝜇, 𝜎), and a response 𝑌 = ℳ(𝑥 1 , 𝑥 2 ) we want to estimate using a PCE. First, we would have to perform an isoprobabilistic transform to reduce the input variables

𝑥 1 = (𝑎 + 𝑏) 2 + (𝑏 -𝑎) 2 × 𝜉 1 , 𝜉 1 ∼ 𝒰(-1, 1), (2.17
)

𝑥 2 = 𝜇 + 𝜎𝜉 2 , 𝜉 2 ∼ 𝒩 (0, 1). (2.18)
Now, since 𝑥 1 has a uniform distribution and 𝑥 2 has a normal distribution, we need Legendre and Hermite polynomials to build the basis. Hence, a multivariate polynomial of the basis reads

Ψ 𝛼 1 ,𝛼 2 (𝜉 1 , 𝜉 2 ) = 𝑃 𝛼 1 (𝜉 1 ) √ 1 2𝛼 1 +1 × 𝐻𝑒 𝛼 2 (𝜉 2 ) √ 𝛼 2 ! (2.19)
Assuming a second order truncation (𝑝 = 2), the basis would be formed by the polynomials in table 2.5. Moreover, the PCE reads

Ŷ = ℳ 𝑃 𝐶 (𝜉 1 , 𝜉 2 ) = 5 ∑ 𝑗=0 𝑦 𝑗 × Ψ 𝑗 = 𝑦 0 + 𝑦 1 √ 3𝜉 1 + 𝑦 2 𝜉 2 + 𝑦 3 √ 5 2 (3𝜉 2 1 -1) + 𝑦 4 √ 3𝜉 1 × 𝜉 2 + 𝑦 5 1 √ 2 (𝜉 2 2 -1).
(2.20) Table 2.5: Base of the PCE assuming a second order truncation.

𝑗 𝛼 Ψ 𝛼 ≡ Ψ 𝑗 0 [0,0] Ψ 0 = 1 1 [1,0] Ψ 1 = √ 3𝜉 1 2 [0,1] Ψ 2 = 𝜉 2 3 [2,0] Ψ 3 = √ 5 2 (3𝜉 2 1 -1) 4 [1,1] Ψ 4 = √ 3𝜉 1 × 𝜉 2 5 [0,2] Ψ 5 = 1 √ 2 (𝜉 2 2 -1)
It is worth noting that we would require a response sample set of at least 12 runs to compute the coefficients.

Kriging [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Lataniotis | UQLab user manual -Kriging (Gaussian process modeling)[END_REF] Let us consider a linear regression model crossing the 𝑦-axis at zero, i.e., 𝛽 0 = 0. The model would be of the form

𝑌 = 𝛽 𝑇 𝑋 + 𝜀, (2.21)
where, once again, 𝑌 is the output variable, 𝑋 is the input vector, 𝛽 is the vector of regression coefficients, and 𝜀 is the error of the model.

To implement a different kind of regression, the input vector can be projected into an 𝑁 dimensional feature space using a set of mapping functions. The model becomes then

𝑌 = 𝛽 𝑇 𝐹 + 𝜀, (2.22)
where 𝐹 = [𝑓 𝑘 (𝑥 𝑖 )] represents the user-selected basis functions that transform the original vector, 𝛽 is the vector of the functions coefficients, the term 𝛽 𝑇 𝐹 is a deterministic function approximating the mean value of the output, and 𝜀 is the error.

Kriging, also known as Gaussian process regression (GPR) in machine learning (ML), is a popular interpolation method that assumes that the model output is a realization of a stochastic process by considering the error 𝜀 as a zero mean, unitvariance, stationary Gaussian process (GP)

𝑌 = 𝛽 𝑇 𝐹 + 𝜀(𝑋), (2.23)
𝜀(𝑋) ∼ 𝐺𝑃 (0, 𝑘(𝑋, 𝑋 ′ )), 𝑘(𝑋, 𝑋 ′ ) = 𝜎 2 𝑅(𝑋, 𝑋 ′ ; 𝜃), (2.24)
where 𝜎 is the variance of the GP, 𝑅(𝑋, 𝑋 ′ ; 𝜃) is the correlation function, also known as the kernel, which describes the correlation between two sample points in the output space, and 𝜃 is the hyperparameters of the kernel function. The mean value of the output 𝛽 𝑇 𝐹 is called the Kriging trend.

A GP is defined as a collection of random variables such that the examination of a larger set of variables does not change the distribution of a smaller set. Moreover, a

Kriging model predicts the responses assuming that the vector formed by the true responses and the predictions has a joint Gaussian distribution. Hence, the parameters 𝛽, 𝜎 2 , and 𝜃 can be derived by maximizing the likelihood function

ℒ(𝛽, 𝜎 2 , 𝜃; 𝑌 ) = (𝑑𝑒𝑡 𝑅) -1/2 (2𝜋𝜎 2 ) 𝑁/2 𝑒𝑥𝑝 [ - 1 2𝜎 2 (𝑌 -𝐹 𝛽) 𝑇 𝑅 -1 (𝑌 -𝐹 𝛽) ] , (2.25)
where 𝑅 is the correlation matrix with elements

𝑅 𝑖𝑗 = 𝑅(𝑥 𝑖 , 𝑥 𝑗 ; 𝜃), 𝑖, 𝑗 = 1, . . . , 𝑁. (2.26)
To build a Kriging model, it is necessary to select a kernel and a functional basis of the Kriging trend. The most commonly used kernels are summarized in table 2.6.

Depending on the trend type different names have been given to Kriging models:

simple Kriging when the trend is a combination of fully specified functions, ordinary Kriging when the trend has a constant yet unknown value, and universal Kriging when the trend is a linear combination of prescribed arbitrary functions. Typically, the trends are based on a polynomial basis.

Because Kriging uses the spatial correlation between the data to interpolate the values, it is an exact interpolation method at the sample points. Moreover, the prediction error could be minimized by adding sample points with important spatial correlation.

Linear 𝑚𝑎𝑥 ( 0, 1 -|𝑋-𝑋 ′ | 𝜃 ) Exponential 𝑒𝑥𝑝 [ -|𝑋-𝑋 ′ | 𝜃 ] Gaussian 𝑒𝑥𝑝 [ -1 2 ( |𝑋-𝑋 ′ | 𝜃 ) 2 ] Matern 3/2 ( 1 + √ 3|𝑋-𝑋 ′ | 𝜃 ) 𝑒𝑥𝑝 [ - √ 3|𝑋-𝑋 ′ | 𝜃 ] Matern 5/2 ( 1 + √ 5|𝑋-𝑋 ′ | 𝜃 + 5 3 ( |𝑋-𝑋 ′ | 𝜃 ) 2 ) 𝑒𝑥𝑝 [ - √ 5|𝑋-𝑋 ′ | 𝜃 ]
Decision trees [START_REF] Hastie | Tree-Based Methods[END_REF][START_REF] Shalev-Shwartz | Decision Trees[END_REF] A decision tree is a simple yet powerful supervised ML algorithm, which can be used for both regression and classification tasks. It predicts the output associated with an input by traveling from the root node of the tree to a leaf. Each internal node of the tree corresponds to a test on the input attributes, and each leaf specifies a value to be returned. Hence, a "decision" is made by performing a sequence of tests that leads to a leaf. In regression trees, instead of a single value, there is a function of the input attributes at each leaf.

Tree-based methods split the feature space into cells, where each cell corresponds to a leaf of the tree. As an example, consider a continuous output response 𝑌 , two input variables 𝑋 1 , 𝑋 2 , and a partition of the feature into five regions (𝑅 1 , . . . , 𝑅 5 ).

At each region, a value is assigned to the output response. To define the regions, first the space is split into two regions at 𝑋 1 = 1, then the region 𝑋 1 ≤ 𝑡 1 is split at An important advantage of decision trees compared to other ML algorithms is their interpretability. One of their biggest drawbacks is that they have high variance, and therefore, a small change to the data could result in major changes to the model. Decision trees are the fundamental component of bagged trees and a random forest, which are among the most powerful ML algorithms today. To make better predictions, both algorithms combine the predictions from different decision trees constructed using bootstrapped learning sets. The fundamental difference between the algorithms is that all features are considered for splitting a node in bagged trees, whereas in a random forest, the learning algorithm is limited to a random sample of features [START_REF] Géron | Ensemble Learning and Random Forest[END_REF].

𝑋 2 = 𝑡
Neural networks [START_REF] Russell | Artificial intelligence: a modern approach[END_REF][START_REF] Shalev-Shwartz | Neural Networks[END_REF] Artificial neural networks (ANNs) are inspired by the architecture of biological neural networks in the brain. They are simply a collection of basic units (neurons) connected in a communication network by directed links. The properties of the network depend on its topology and the properties of the neurons. If all the connections are in one direction, the network is called a feed-forward network. If there are cycles or loops, it is called a recurrent neural network (RNN).

A numeric weight 𝑤 𝑖,𝑗 is associated to every link, and an activation function, typically a hard threshold or a logistic function, is associated to the neurons. Each neuron receives as input the weighted outputs of other neurons connected to its incoming edges. The neuron computes the sum of its inputs and then, to compute its output, it applies the activation function to the sum. The weights are the unknown coefficients one needs to adjust when training the neural network. Typically, a back-propagation algorithm is used to update the values until the gradient of the loss function with respect to the weights is minimized.

In feed-forward networks, the neurons are generally organized in consecutive layers and the neurons of one layer receive only the weighted outputs of the neurons in the preceding layer (see e.g. figure 2345678). The layer connected to the network's inputs is called the input layer, the layer connected to the output(s) is called the output layer, and all the layers in between are called hidden layers. The input layer always contains a bias neuron, which outputs 1. Often, hidden layers also contain bias neurons.

The rectified linear (ReLU) activation function works well in most cases for the hidden layers. The activation function for the output layer depends on the problem. For binary classification, sigmoid and tanh are generally used, but softmax is usually better for multi-classification. For regression tasks, the output layer does not necessarily have an activation function.

The simplest neural network contains one hidden layer, and it can be sufficient for simple problems. However, to solve problems that are not linearly separable, a substantial number of hidden layers might be required.

Global sensitivity analysis [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]117] Sensitivity analysis is a valuable mathematical tool to study how the variability of the model response is affected by the variability of each input variable. It is useful, among others, to understand the behavior of the model response, reduce the dimensions of the problem, prioritize the uncertainty sources, and conduct robustness analysis.

Numerous approaches to perform sensitivity analysis have been developed throughout the years. They can be classified based on their cost in terms of model evaluations, their assumptions about the model, and their requirements of a priori knowledge about the model. Also, they are typically grouped in local and global methods. In local sensitivity analysis, the impact of small perturbations around the nominal values of the inputs is studied, whereas in global sensitivity analysis, the whole input domain is taken into account. One of the most popular global methods is Sobol' sensitivity indices.

Sobol' indices

Sobol' indices are based on the idea of expanding the computational model into summands of increasing dimensions

ℳ = ℳ 0 + 𝑀 ∑ 𝑖=1 ℳ 𝑖 (𝑥 𝑖 ) + ∑ 1≤𝑖<𝑗≤𝑀 ℳ 𝑖𝑗 (𝑥 𝑖 , 𝑥 𝑗 ) + • • • + ℳ 12...𝑀 (𝑋) = ℳ 0 + ∑ u⊂{1,...,𝑀 } ℳ u (𝑥 u ), (2.27)
with 𝑥 u def = {𝑥 𝑖1 , . . . , 𝑥 𝑖𝑠 }. Thus, the total variance of the model can be seen as the sum of the variances of summands

𝐷 ≡ Var[ℳ(𝑋)] = ∑ u⊂{1,...,𝑀 } Var[ℳ u (𝑥 u )].
(2.28)

𝑆 u def = Var[ℳ u (𝑥 u )] 𝐷 .
(2.29)

The first-order Sobol' indices quantify the additive effect of each input variable (𝑥 𝑖 ) separately:

𝑆 𝑖 = 𝐷 𝑖 𝐷 = Var[ℳ 𝑖 (𝑥 𝑖 )] 𝐷 .
(2.30)

The total Sobol' indices quantify the total effect of the input variable 𝑥 𝑖 , including the interactions with the other variables:

𝑆 𝑇 𝑖 def = ∑ u⊃1 𝑆 u . (2.31)
Multiple-term indices that quantify the effect of the interaction of the variables are referred to as high-order Sobol' indices.

Validation

The accuracy of a computational model can be quantified by estimating its error, which is the deviation of the value it predicts ( Ŷ ) from the "true value" of the quantity of interest (𝑌 ). Once again, in this work, we consider the "true values" as the results obtained from full-wave simulations. Estimating the error of the 3-D model we are using to carry out the simulations is not straightforward. Even if we included in the model all the elements of the electromagnetic environment, there would still be uncertainty associated with the physical phenomena. This is when experiments play an important role. To corroborate certain hypotheses and make sure that the numerical model roughly reflects the behavior of the original system, one could compare simulation results to measurements. Note, however, that the measurement equipment and setup also come with uncertainty.

Estimated error

Several measures of error can be found in the literature.

Measure of error Formula

Empirical error

∑ 𝑁 𝑖=1 (𝑌 -Ŷ ) 2 ∑ 𝑁 𝑖=1 (𝑌 -𝜇 𝑌 ) 2
Mean absolute error (MAE)

1 𝑁 ∑ 𝑁 𝑖=1 |𝑌 -Ŷ | Mean squared error (MSE) 1 𝑁 ∑ 𝑁 𝑖=1 (𝑌 -Ŷ ) 2
Root mean square error (RMSE) √ 𝑀 𝑆𝐸 Normalized root mean square error (NRMSE)

𝑅𝑀 𝑆𝐸 𝑌𝑚𝑎𝑥-𝑌 𝑚𝑖𝑛

We could use a dataset to build a computational model and then estimate an error by applying the formulas in table 2.7. Nevertheless, when comparing the values predicted by the model to the values we used to adjust its coefficients, we are not estimating the prediction error. In that case, we would be estimating the fitting error, also known as the training error.

To estimate the prediction error, we need a validation set formed by an independent set of inputs and outputs. Alternatively, one can build the model and validate it using a 𝑘-fold cross-validation technique [START_REF] Hastie | Cross-Validation[END_REF]. The technique is well-known in ML.

It consists in dividing the dataset into 𝑘 different groups (folds), then estimating the coefficients of the computational model 𝑘 times, each time using the data of 𝑘-1 groups, and calculating the error with the data that was held out. The final values of the coefficients and the error of the model are calculated as the average of the results over the 𝑘 rounds of the 𝑘-fold cross-validation.

The leave-one-out error (𝐸 𝐿𝑂𝑂 ), a technique developed in statistical learning to overcome overfitting, can be seen as the combination of the empirical error of the model over a 𝑁 -fold cross-validation. It consists in using each point of the experimental design as a single validation point for the model built from the remaining 𝑁 -1 points. The 𝐸 𝐿𝑂𝑂 is then calculated as the mean-square average of the resulting 𝑁 empirical errors.

All the regression techniques and algorithms introduced in this chapter will be used in the following chapters. Nevertheless, since the interpretability of some ML models is questionable, we will privilege PCE and Kriging models over other techniques to estimate the peak-values of the magnetic field and the surges.

Chapter 3 Transient magnetic field inside a building struck by lightning

To define protection measures against the LEMP in industrial facilities, it is essential to estimate the magnetic field generated by a direct strike inside the buildings.

However, as we mentioned in the previous chapter, the accuracy of the estimations depends on the representativeness of the models, which relies on the assumptions we make along the way. Ideally, to estimate the transient magnetic field, we would consider all the components of the buildings and carry out a full-wave simulation.

Yet, a full-scale reinforced concrete building can be geometrically very complex; thus, including every detail in the model is virtually impossible. Besides, the simulation would require a large amount of computational resources. Even with the assumptions described in section 2.2, we would find ourselves with thousands and thousands of thin-wires. Carrying out these simulations on a daily basis would not be feasible nor practical.

Simplified 3-D modeling of reinforcing grids

Before attempting to build a surrogate model to estimate the magnetic field inside the buildings, it is worth trying to simplify the representation of the reinforcement in full-wave simulations. If we were able to transform a multi-layered reinforcing into an equivalent single-layered reinforcing grid, we would facilitate the development of the 3-D model and probably reduce the amount of computational resources required.

Distribution of the lightning current in the reinforcement

To simplify the representation of the reinforcement, first, we must analyze the distribution of the lightning current in the reinforcing grid. Let us then consider a 5 m × 5 m × 5 m reinforced concrete building over a perfectly conducting ground plane with a single-layered reinforcing grid. The reinforcing grid is adapted from [START_REF] Maksimowicz | Investigation of models of grid-like shields subjected to lightning electromagnetic field: experiments in the frequency domain[END_REF]. It has squared meshes of 30 cm × 30 cm and is embedded in 30 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.4 mm, and their conductivity is set to 8.33 ×10 6 S/m. The current waveform is defined as a Gaussian function covering a frequency range up to 1 MHz and using TEMSI-FD [START_REF]TEMSI-FD: Time Electromagnetic Simulator -Finite Difference Time Domain[END_REF], the currents in the reinforcement are computed at 1 kHz and 1 MHz. The results, normalized to the current waveform, are shown in figure 3-2. It can be observed that contrary to the implications of the formulas in the IEC 62003 standard [20], and in agreement with [START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF] and [START_REF] Metwally | Reduction of lightning-induced magnetic fields and voltages inside struck double-layer grid-like shields[END_REF], the current is not uniformly distributed in the reinforcement; it is diverted to the edges. Now, let us consider a double-layered reinforcing grid with 30 cm between the layers and hoops every 90 cm. Unexpectedly, figure 34shows that the magnitude and the distribution of the lightning current do not appear to be affected by either the frequency or the number of layers. One could then think that if the currents computed in the outer layer of a multi-layered reinforcing grid are in the same range as the current on a single-layered grid with the same characteristics, it would not be necessary to model the inner layers to compute the surges in underground cables.

Also, if we take a closer look, we can see that there is almost no derivation of the lightning current to the hoops, which suggests that the currents in the inner layer are primarily induced. It may therefore be possible to calculate the characteristics of a single-layered reinforcing grid that would provide comparable shielding effectiveness, at least at a certain distance from the wall.

Surrogate model of single-layered grids

To find the characteristics of an equivalent single-layered reinforcing grid, we propose the methodology in figure 345. The idea is simple: we need a reduced representation of reinforcing grids that will allow us to characterize the reinforcement of the building. We build a surrogate model of the reduced representation of single-layered grids and then, to find the equivalent grid of a multi-layered grid, we use an optimization algorithm to fit the surrogate model to the results obtained from the reduced representation of the multi-layered grid we want to simplify.

Finding a reduced representation of reinforcing grids is not straightforward. Originally, we considered three reduced representations:

• An infinite grid illuminated by a plane wave at normal incidence Infinite grids have been used on numerous occasions to study the transmission and reflection properties of reinforced concrete, see e.g. [START_REF] Casey | Electromagnetic shielding behavior of wire-mesh screens[END_REF][START_REF] Richalot | Electromagnetic propagation into reinforced-concrete walls[END_REF][START_REF] Dalke | Effects of reinforced concrete structures on RF communications[END_REF][START_REF] Pena | Measurement and modeling of propagation losses in brick and concrete walls for the 900-MHz band[END_REF][START_REF] Hyun | Analysis of shielding effectiveness of reinforced concrete against high-altitude electromagnetic pulse[END_REF]]. Yet, they are not appropriate to characterize the reinforcing grid when a building is struck by lightning, because as mentioned in section 2.2, the shielding characteristics of reinforced concrete are different for direct and indirect strikes. Hence, a simplified representation that would be appropriate for studying nearby strikes may be unsuitable for direct strikes and vice versa.

A 2 m × 2 m × 2 m reinforced concrete building is a good representation as long as 1 m is a multiple of the mesh size, which is rarely the case. Otherwise, a rebar would not be passing by the center of the roof, and therefore, we would have to add connection wires, which can induce an error in the results. In addition, the computation time seems to be considerably high to think of carrying out multiple simulations to build the surrogate model. Thus, this representation is also rejected.

Different variations of the remaining option, the 4 m × 4 m grid, are then considered and improved to find the most appropriate representation. The one we finally selected is an 8 m × 8 m grid in a perfectly conducting metallic enclosure (see figure 34567). The grid is connected on one side to the enclosure and on the other side to a perfect voltage source, which is connected to the enclosure. The voltage waveform is defined as a Gaussian function with a maximum amplitude of 100 V, covering a frequency range up to 1 MHz. If there are multiple layers, the first layer of the grid is the one connected to the enclosure. To minimize the effect of the wires interconnecting the grid and the enclosure, their radius is set to 2.5 cm, and they are considered to be perfectly conducting. The conductivity of the rebars is set to 8.33 ×10 3 S/m.

Note that it is the conductivity of steel divided by 10 3 . This change, as well as using a metallic enclosure, is made to reduce the computation time. 3.2 shows that, based on the NRMSE estimated by applying a 10-fold cross-validation [START_REF] Hastie | Cross-Validation[END_REF], the GPR algorithms Fine Tree

• Minimum leaf size: 4

1.835 × 10 -2 make the best predictions. Hence, we focused on obtaining ordinary Kriging models.

To select the kernel, different Kriging models are built and compared using UQLab [START_REF] Marelli | UQLab: A framework for uncertainty quantification in MATLAB[END_REF]. Assuming that the inputs have a uniform distribution and based on the 𝐸 𝐿𝑂𝑂 , the best models are found with the Matern 5/2 kernel. The parameters of these models are estimated using a Hybrid Genetic Algorithm (HGA) with a convergence tolerance of 10 -10 and a maximum number of generations equal to 500. The values estimated for the parameters are given in table 3.3.

Model fitting criteria

Going back to the objective, which is to find the characteristics of an equivalent single-layered grid, what we need is to estimate the components of the input vector to be used in the Kriging models. The input vector should lead to a current and a magnetic field strength similar to the ones resulting from the reduced representation of the multi-layered grid we want to simplify. To estimate the components of the input vector, namely the radius of the rebars and the mesh size of an equivalent single-layered grid, the Kriging models are fitted to the response curves of the reduced representation of the multi-layered grid using least squares. The procedure followed by the fitting algorithm is summarized in figure 3456789.

Note that for the fitting, the tolerance on the model and the values of the coefficients is set to 10 -10 .

We could estimate the components of the input vector by fitting only one of the models to a response curve, as well as by starting the optimization from a random value. However, the solution is not always unique, and in some cases, it might not even be accurate. That is the reason why the fitting algorithm compares the values obtained with two independent models and the error when starting from different values. If there were two solutions, we could find them by starting from the lower and the upper bounds of the models, i.e., the minimum and the maximum values of the factorial design.

We observed that when starting the fitting from the lower bound, the simplification of the multi-layered grid is usually compensated by a smaller mesh size, while starting from the upper bound, it is compensated by a bigger radius. Rarely the second approach is a better fit and the algorithm manages to find an equivalent grid without drastically reducing the mesh size. Even when it does, the first approach always leads to a lower error. The latter means, unfortunately, that an equivalent singlelayered grid might not contribute to reducing the amount of computational resources required in the simulations. For example, if we implement the fitting algorithm to find the characteristics of an equivalent single-layered grid of the double-layered grid considered in section 3.1.1, we would obtain meshes of 10 cm and rebars with a radius of 1.13 mm. Since the mesh size is considerably smaller than the original, the only gain we could achieve would be in terms of memory.

There are ways to make sure that the characteristics of the equivalent grid are convenient and the computation time is reduced. Yet, a choice has to be made between approximating the impedance of the reinforcement and approximating its shielding effectiveness. For example, if the factorial design is expanded to take into account thicker rebars and the conductivity is considered as an additional variable, fitting the models would lead to different results. When fitting the current flowing through the same double-layered grid, we obtain a single-layered grid with a mesh size of 44 cm, rebars with a radius of 1.90 cm, and a conductivity of 8.53 ×10 4 S/m. When fitting the magnetic field strength, we obtain a single-layered grid with a mesh size of 30 cm, rebars with a radius of 2.85 cm, and a conductivity of 1.57 ×10 6 S/m. Also, if instead of expanding the factorial design, we consider the imaginary part of the impedance as an output of the reduced representation of the double-layered grid, the equivalent single-layered grid would have a mesh size of 36 cm and rebars with a radius of 1.34 cm.

To compare these four single-layered grids, the magnetic field strength is computed inside and around the building described in section 3.1.1, considering a lighting strike to the corner. Figure 3-10 shows the results at 1 MHz on a 4 m × 4 m surface, located 1 m away from the roof. Figure 3-11 shows the relative error on the values computed. It can be observed that, in general, the distribution of the magnetic field strength is very similar. As expected, there is a high error in the gap between the two layers of the original grid, which is completely normal. There is also a high error in approximately 5 cm next to the walls because the equivalent grid is placed where the outer layer of the original grid was. Moreover, the inhomogeneity in the field created by the rebars in the inner layer is no longer visible. It would not be reasonable to place the equivalent grid in the middle of the wall or as one of the inner layers, since we already know that most of the lightning current is displaced to the outer grid.

It is worth noting that all four equivalent grids are able to reproduce the field around the building; yet, only the grids found by fitting the magnetic field strength reproduce accurately the field inside. The same conclusions can be drawn from the results on the other surfaces and at lower frequencies. Also, there is always a higher error close to the edges, where the values change fast. An improvement could be achieved by reducing the size of the FDTD cell size.

Certainly, it would be better to use the fitting algorithm to transform the multilayered grid into a single-layered that matches both the field and the impedance.

Nevertheless, it may not always be advantageous in terms of computational resources.

The alternative simplifications could be useful to study separately other indirect effects of lightning, such as the surges induced in the cables. However, to be able to choose the most convenient equivalent grid depending on the case study, we need to understand how the magnetic field strength is influenced by the characteristics of the reinforcement. It would also give us an insight into the approach of the fitting algorithm.

Parametric study on the shielding effectiveness of reinforced concrete

Different studies have shown that the LEMP inside of a reinforced concrete building is attenuated when the mesh size is reduced or the radius of the rebars is increased [START_REF] Maksimowicz | Investigation of models of grid-like shields subjected to lightning electromagnetic field: experiments in the frequency domain[END_REF][START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF]. Also, that the magnetic shielding effectiveness can be enhanced by increasing the number of layers [START_REF] Maksimowicz | Investigation of models of grid-like shields subjected to lightning electromagnetic field: experiments in the frequency domain[END_REF][START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF][START_REF] Metwally | Reduction of lightning-induced magnetic fields and voltages inside struck double-layer grid-like shields[END_REF][START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF][START_REF] Zischank | Laboratory simulation of direct lightning strokes to a modeled building: measurement of magnetic fields and induced voltages[END_REF][START_REF] Kern | Magnetic fields and induced voltages in case of a direct strike-comparison of results obtained from measurements at a scaled building to those of IEC 62305-4[END_REF][START_REF] Mardiguian | Lightning-generated fields in reinforced concrete buildings[END_REF], especially for fast-rising currents [START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF], and by enlarging the distance between the layers [START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF]. In this section, we will conduct a basic parametric study aiming to quantify the variations in the shielding against the magnetic field when the characteristics of reinforced concrete are modified.

Let us then consider once again a 5 m × 5 m × 5 m reinforced concrete building, this time with a foundation buried at a depth of 1 m. The walls, roof, and foundation of the building are made up of a reinforcing grid, the characteristics of which we are going to vary. Based on real configurations, six parameters are chosen to be analyzed independently, starting from the reference values in table 3.4. When one of the parameters is varied, the others are set to their reference values. Note that the internal walls, the columns, and the beams are not taken into account. The meshes are considered squared, the layers are considered to be facing each other, which means that the meshes of the layers are aligned, and the radius of the hoops is equal to the radius of the rebars. The conductivity of the rebars is set to 8.33 ×10 6 S/m.

The concrete is modeled as a lossy dielectric material with variable conductivity and permittivity, in accordance with its moisture content. The soil is considered homogeneous, with a relative permittivity of 10 and a resistivity of 100 Ω•m.

Table 3. The magnetic shielding effectiveness is defined as follows:

𝑆𝐻 = 20 × log ( 𝐻 0 𝐻 ) , (3.1)
where 𝐻 is the magnitude of the magnetic field strength computed inside the reinforced concrete building, and 𝐻 0 is the magnitude of the magnetic field strength computed when the LPS system of the building consists of four down-conductors.

The down-conductors, one in each corner, are connected at the top and the bottom, as shown in figure 3-13. Their radius is set to 6 mm, the reference value of the radius of the rebars. Since the computed values exhibit a minor dependency on the frequency, we decided to plot them against the distance of the points to the foundation and the striking point. It can be observed that the shielding effectiveness increases with the distance to the foundation and as we approach the roof. Lower shielding effectiveness close the foundation could be explained by the electromagnetic field generated by the currents flowing through the foundation as they dissipate into the soil. The interpretation of higher shielding effectiveness close to the roof is less straightforward. It is perhaps higher because the reinforcing grid in the roof provides shielding against the fields radiated by the lightning channel, whereas there is no shielding when the LPS consists of four down-conductors.

Mesh size and radius of the rebars forming the reinforcing grid

As expected, figures 3-15 and 3-16 show that the shielding against the magnetic field improves as the percentage of steel in the structure increases. The effect is therefore consistent with the behavior of reinforced concrete to an incident wave. The shielding effectiveness increases with the radius of the rebars and decreases with the mesh size.

On average, dividing the mesh size by two enhances the shielding against the magnetic field by 9 dB, and multiplying by two the radius of the rebars enhances the shielding by 2.7 dB. This explains why the fitting algorithm presented in section 3.1.3 usually reaches the lower limit of the mesh size before it starts changing the radius.

Number of layers and periodicity of the hoops

Adding an extra layer to the reinforcement would of course enhance the shielding against the magnetic field as well. We saw in section 3.1.1 that the lightning current flows mainly through the outer layer. Thus, an inner layer would work as a shield against an incident electromagnetic field. Accordingly, the number of interconnections between the layers should not influence the shielding effectiveness. In agreement with [START_REF] Metwally | Reduction of lightning-induced magnetic fields and voltages inside struck double-layer grid-like shields[END_REF], the results in figure 3-17 confirm that the effect of increasing the periodicity of the hoops is negligible. Only placing a hoop at every intersection improves slightly the shielding, but it is not by more than 1 dB.

Figure 3-18 shows that on average, the magnetic shielding effectiveness increases by 12.5 with an additional layer. Yet, the increase depends on the mesh size and the distance between the layers. If the distance between the layers is divided or multiplied by two, the shielding effectiveness would increase by 12.5 dB ± 2.7 dB. For example, consider the shielding effectiveness when all the parameters are set to the reference values and the results in figure 3

Moisture content of the concrete

It has already been demonstrated that concrete does not contribute significantly to the magnetic shielding [START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF][START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF]. Figure 3-19 shows that as long as the concrete is modeled, the effect of changing its moisture content is negligible. An effect is noticeable only at low frequency and when the concrete is disregarded. At most, disregarding the concrete leads to an increase of 0.5 dB.

Permeability of steel

Steel is an iron alloy, which has high conductivity and a considerably variable degree of magnetization. The permeability of steel depends on its chemical composition. For example, the relative permeability of carbon steel is on the order of 100, whereas the relative permeability of annealed stainless steel is within 750 and 1800 and can be approximated to 1 for austenitic stainless steel [START_REF] Ellingson | Constitutive Parameters of Some Common Materials[END_REF]. The most common type of steel used in structural applications is carbon steel, which is generally assumed to be nonmagnetic. The IEC 62305-4 standard [START_REF]Electrical and electronic systems within structures[END_REF] suggests considering a relative permeability of 200 for the steel used in reinforced concrete buildings.

A high degree of magnetization is known to contribute to the shielding against an incident field; yet, when a current is injected into the reinforcement, the effects of the permeability are contradictory. It is however worth noting that a good agreement between the effects has been found in [START_REF] Laisné | Eeléments de modélisation de la foudre sur les bâtimentspartie 3 : simulation numérique 3d d'un bâtiment soumis à un impact direct et à une boucle d'injection[END_REF]. When a current is injected into the reinforcement, the skin depth (𝛿 𝑖 ) decreases with both the permeability and the frequency, and when the skin depth decreases, the resistance of the rebars increases and their internal inductance decreases. In addition, the internal inductance of a conductor depends on the degree of magnetization of the material, and therefore, it also increases with the permeability.

By neglecting the skin effect and the internal inductance, we could approximate the internal impedance of the rebars as 𝑍 𝑖𝑛𝑡 = 𝑅 𝑑𝑐 = 0.001 Ω/m, considering a radius of 6 mm. If we consider the skin effect and a worst-case scenario with a relative permeability 𝜇 𝑟 = 1000, the impedance of the rebars would depend on the frequency.

At the frequencies characterizing the different return strokes in [START_REF]Electrical and electronic systems within structures[END_REF], the resistance, approximated as 𝑅 𝑎𝑐 = 𝜌 2𝜋𝑎𝛿 , is equal to 0.0913 Ω/m, 0.289 Ω/m, and 0.578 Ω/m for 25 kHz, 250 kHz, and 1 MHz, respectively. The internal inductance, approximated as 𝐿 = 𝜇𝛿 4𝜋𝑎 , is equal to 0.581 𝜇H/m, 0.184 𝜇H/m, and 0.092 𝜇H/m. Considering these values, we carry out the simulations at each frequency. We vary the resistance and the inductance independently because the internal reactance of cylindrical conductors is often neglected. Figure 3-20 shows that the permeability could considerably reduce the magnetic shielding effectiveness. Nevertheless, the magnitude of the effect depends on the frequency, and it is noticeable only when considering a high degree of magnetization.

If the relative permeability of the type of steel used in reinforced concrete buildings is indeed a few hundred, one may not need to take it into account in the simulations. To summarize, the shielding provided by reinforced concrete does depend on the geometrical configuration of the reinforcement and the type of steel the rebars are made of. The most effective way to improve the magnetic shielding is by adding multiple layers. The further the inner layers are placed from the outer layer, the better the shielding will be. A significant enhancement can also be obtained by reducing the mesh size. The effect of the radius of the rebars and the periodicity of the hoops is insignificant compared to the influence of these two parameters. Also, the variations in the shielding effectiveness observed at all the 29 points seem to be comparable. Hence, one could assume that modifying the geometrical configuration of the reinforcement will affect the magnetic field uniformly inside the building. It would therefore be possible to calculate a correction factor to take into account the modifications on the reinforcement, at least when it comes to estimating the peakvalues of the magnetic field strength. It has been observed that the magnetic field has a longer rise time inside multi-layer grid-like shields [START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF].

On that idea, we could choose a single-layer grid-like shield as the reference shield and divide the problem of estimating the magnetic field strength into two parts. In the first part, we would need to estimate the distribution of the peak-values inside the reference grid-like shield. In the second part, we could calculate a correction factor to adjust the estimated values, depending on the characteristics of the grid-like shield of interest. Note that we are now referring to the reinforced concrete building as a grid-like shield. This is to emphasize that we are not considering all the structural components of the buildings, such as walls, columns, and beams.

Distribution of the peak-values of the magnetic field strength

Calculations according to IEC 62305-4

The IEC 62305-4 standard [START_REF]Electrical and electronic systems within structures[END_REF] suggests the following formula to calculate the maximum value of the magnetic field strength in single-layer grid-like shields:

|𝐻| 𝑚𝑎𝑥 = 𝑘 × 𝐼 0 × 𝑤 𝑚 𝑑 𝑤 × √ 𝑑 𝑟 , (3.2)
where 𝑘 is the configuration factor (typically 𝑘 = 0.01), 𝐼 0 is the maximum value of the lightning current, 𝑤 𝑚 is the mesh width, 𝑑 𝑤 is the shortest distance to the wall, and 𝑑 𝑟 is the shortest distance to the roof. The formula is a first-order approximation of the magnetic field distribution inside the safety volume, the boundaries of which are defined by the safety distance 𝑑 𝑠 = 𝑤 𝑚 from the shield. It was obtained from numerical simulations of three typical grid-like shields over a perfectly conducting To analyze the limitations of the formula, we start by computing the magnetic field under the same considerations inside the 10 m × 10 m × 10 m shield. The radius of the rebars is set to 6 mm and the simulations are carried out using TEMSI-FD. In addition, the IEC formula not only disregards the effect of the radius of the rebars but also assumes a linear relationship between the mesh size and the magnetic field strength. Figure 3-21 shows that is not entirely accurate. To focus on building a model to calculate a correction factor, we could set the mesh size to a fixed value and use the IEC formula to estimate the peak-values. After all, the values calculated with the formula are mostly on the safe side. On the other hand, we could propose a more precise, yet more complex, variant.

Case studies

To propose an improved formula, first, we define the reference shield as a singlelayered reinforcing grid with squared meshes of 10 cm × 10 cm, formed by rebars with a radius of 1 mm. The grid is embedded in 10 cm of plain concrete to avoid direct contact of the foundations with the soil. Then, we compute the peak-values of the magnetic field strength inside the ten different grid-like shields listed in table 3.5.

Note that the dimensions of the shields include the thickness of the concrete.

We could consider the reference shield as a reinforcing grid with more realistic characteristics; yet, since the mesh size and the radius will be variables of interest when computing a correction factor, we decided to set the values to the potential minimum. If the effect of the radius of the rebars is disregarded, we assure by choosing the thinnest rebars that we remain on the safe side. Also, based on the assumptions of the IEC formula, we might be able to evaluate the effect of the mesh size independently. Hence, we choose the smallest mesh size, which defines the greatest safety Note: The dimensions of the buildings include 10 cm of plain concrete.

volume. A safety distance equal to the mesh size must still be considered for the field to be homogenized. At a shorter distance to the reinforcement, the field fluctuates rapidly depending on the position of the rebars.

All the simulations are carried out using TEMSI-FD and executed in parallel. The boundary conditions of the analysis space are defined as PMLs to assume an open space. The volume of the analysis space is variable because the boundaries are defined at least 10 m away from the grid-like shield. Again, the concrete is modeled as a lossy dielectric material with a conductivity of 0.0052 S/m and a relative permittivity of 8.6. The conductivity of the rebars is set to 8.33 × 10 6 S/m. The soil is considered homogeneous, with a relative permittivity of 10 and a resistivity of 100 Ω•m. As usual, the lightning channel is modeled as a monopole antenna, and the current waveform is defined as a Gaussian function covering a frequency range up to 1 MHz. The channel is attached to the corner of the roof because direct strikes to the corner are known to produce the highest magnetic fields [START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF].

Finally, the peak-values of the magnetic field strength are computed at 11 horizontal surfaces inside the shields. The first and the last surfaces are located 10 cm away from the roof and the foundation. The remaining nine are distributed uniformly in height, as shown in figure 3-24. Note that we are considering a Cartesian coordinate system expressed in meters, with the origin at the striking point and the 𝑧-axis oriented towards the foundation. Also, the results are once again normalized to the current waveform and therefore treated as transfer functions.

If we take a look at the results at points (1, 1, 1), (3, 3, 3), and (5, 5, 5) inside B1 and B6 in figure 3-25, we see that the magnetic field strength is almost constant up to a few hundreds of kHz, depending on the size of the building. Then, it decreases with the frequency. The tendency is similar at other points and inside different buildings.

As suggested by the IEC formula in (3.2), we could therefore assume that the peakvalues of the magnetic field strength are proportional to the maximum value of the lightning current (𝐼 0 ), at least for slow-rising currents. Thus, to conceive an improved formula, we choose the results at 25 kHz characterizing the first positive return stroke, which has the highest peak current. According to the LPL I defined in [START_REF]3: Physical damage to structures and life hazard[END_REF], we set 𝐼 0 to 200 kA.

Improved formula

A direct lightning strike generates a transient magnetic field inside a building due to the lightning current flowing along the lightning channel and to the earth-termination system. Since the reinforcement of the building works as a shield against incident waves, we could theoretically neglect the contribution of the lightning channel and calculate the magnetic field as the superposition of the contributions given by the flow of a partial lightning current through each rebar. However, to determine the contributions of the rebars, we would need to estimate the partial currents and evaluate integrals as suggested in [START_REF] Orlandi | Systematic approach for the analysis of the electromagnetic environment inside a building during lightning strike[END_REF]. Alternatively, we can start by supposing that a good approximation of the magnetic field strength can be obtained by superposing the contribution of the rebars in the vertical edges and adding a term (𝑇 𝑤 ) to account for the contribution of the walls, roof, and foundation. We know that the lightning current is always diverted to the edges. And to facilitate the development and application of the formula, we could even consider just the closest vertical edge to the point of observation. Let 𝑇 𝑒 then be the term accounting for the contribution of the closest edge. Note that this simplification causes a discontinuity of the field at the points that are equidistant from two or more vertical edges and fails to reproduce the smallest values, where the fields are canceling each other out. Accordingly, the variant of the IEC formula to calculate the peak-values of magnetic field strength is conceived as

|𝐻| 𝑚𝑎𝑥 = 𝐼 0 × ( 𝑇 𝑒 𝑘 + 𝑇 𝑤 ) , (3.3)
where 𝐼 0 is once again the maximum value of the lightning current, currently set to 200 kA, and 𝑘 is a configuration factor, which depends on the position of the closest edge with respect to the striking point.

To find the composition of the terms 𝑇 𝑒 and 𝑇 𝑤 , we consider from each surface the results along the 𝑥-axis, 10 cm away from the wall, and along the line x = y up to the center. As shown in figure 3-26, the magnetic field strength decreases with the Since the term 𝑇 𝑤 is supposed to be constant when the distance to the wall is fixed, by fitting the curves along the 𝑥-axis up to the center, we find that 𝑇 𝑒 can be approximated as 𝑐 0 ×𝑑 -𝑐 3 𝑒 , where 𝑐 0 and 𝑐 3 are coefficients to be determined. However, whereas 𝑐 3 remains around the same values, 𝑐 0 varies with the distance to the roof (𝑑 𝑟 ). Probably because the distribution of the currents changes in the reinforcement as they flow to the earth-termination system. Thus, we need to adjust 𝑐 0 as a function 𝑑 𝑟 . We approximate the function as 𝑐 1 × 𝑑 -𝑐 2 𝑟 and the composition of the first term of the formula is set.

Next, to find the composition of 𝑇 𝑤 , we fit the curves resulting from the difference between the results along the lines 𝑥-axis and x = y on each surface. We also find that 𝑇 𝑤 can be approximated as

𝑐 0 × 𝑑 -𝑐 6 𝑤 with 𝑐 0 = 𝑐 4 × 𝑑 -𝑐 5
𝑟 . Yet, we observed that the coefficient 𝑐 6 varies slightly with the relative distance to the roof. Close to the walls, it decreases as we move along the 𝑧-axis, then increases again close to the foundation. Hence, we introduce a new coefficient 𝑐 7 that multiplies 𝑐 6 , and we define the dependence on the relative position in height as

⃒ ⃒ 𝑐 7 -2×𝑑𝑟 ℎ ⃒ ⃒
, where ℎ is the height of the shield. The value of 𝑐 7 indicates the point at which the field close to the walls stops decreasing. For example, 𝑐 7 would be equal to 1 if the magnetic fields generated by the flow of the current through the foundations and the roof were comparable.

That would indicate that the field close to the walls will stop decreasing right in the middle. With the composition of the second term of the formula also set, the variant introduced in (3.3) has become

|𝐻| 𝑚𝑎𝑥 = 𝐼 0 × ( 𝑐 1 𝑘 × 𝑑 𝑐 2 𝑟 × 𝑑 𝑐 3 𝑒 + 𝑐 4 𝑑 𝑐 5 𝑟 × 𝑑 𝑐 6 ×|𝑐 7 -2×𝑑𝑟 ℎ | 𝑤 ) , (3.4)
where 𝑐 𝑖 , for 𝑖 = 1, . . . , 7, are coefficients to be determined. Note that the distances are in meters and measured from the center of the rebars forming the grid-like shield.

The formula in (3.4) was conceived assuming that the closest vertical edge is the one in the striking point. If it was not the case, the contribution of the rebar in the vertical edge would increase with the distance to the roof. The rapidity of the increase depends on its distance to the striking point. By fitting the ratios of the peak-values computed at the corners of the horizontal surfaces, we found that the configuration factor k could be approximated as follows:

𝑘 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1,
for the edge in the striking point

1 + 𝑑𝑝 3×𝑑𝑟 , for the adjacent edges 1 + 𝑑𝑝 6×𝑑 1.1

𝑟

, for the opposite edge (3.5) where 𝑑 𝑝 is the horizontal distance between the closest vertical edge and the edge in the striking point.

Finally, to consider the effect of the currents flowing through the roof, which is pronounced in the fields computed at the first surface of each building, we include an exponential decay in the second term. Hence, the formula in (3.4) becomes

|𝐻| 𝑚𝑎𝑥 = 𝐼 0 × ( 𝑐 1 𝑘 × 𝑑 𝑐 2 𝑟 × 𝑑 𝑐 3 𝑒 + 𝑐 4 𝑑 𝑐 5 𝑟 × 𝑑 𝑐 6 ×|𝑐 7 -2×𝑑𝑟 ℎ | 𝑤 × 𝑒 𝑐 8 ×𝑑𝑛 ) , (3.6)
where 𝑒 is the exponential function and 𝑑 𝑛 is the distance to the vertical edge in the striking point. By relative, we mean that it is normalized to the dimensions of the shield as

𝑑 𝑛 = √ ( 𝑥 𝑙𝑒𝑛𝑔𝑡ℎ ) 2 + ( 𝑦 𝑤𝑖𝑑𝑡ℎ ) 2 .
Now that we have the expression of the formula, the next step is to determine the coefficients 𝑐 1 to 𝑐 8 . This time, we can use all the data to guarantee a good fit. In a first attempt, however, we observed that both 𝑐 4 and 𝑐 5 depend on the dimensions of the grid-like shield. The coefficient 𝑐 4 decreases with the size of the shield, probably because in a bigger shield, the current flowing through the rebars is redistributed among more parallel paths. Therefore, we defined 𝑐 4 as a function of the volume of the shield (𝑣) and tried to adjust it to the response. A good fit was obtained with

𝑐 4 = 𝑐 04 3 √
𝑣 . The coefficient 𝑐 5 decreases as the ratio of the cross-sectional area (𝑎) to the height of the shield (ℎ) increases. Perhaps because the currents take longer to be displaced to the corners. Thus, we tried reducing 𝑐 5 as a function

√ 𝑎 ℎ , √ 𝑎 instead of
𝑎 to cancel the units. A good fit was obtained when reducing 𝑐 5 by

√ 𝑎 20×ℎ .
The formula we propose as a more precise variant of the formula given by the IEC standard is then 𝜀 : Estimated error With the formula in (3.7), we were able to estimate the coefficients. We used an optimization algorithm to minimize the sum of the NRMSE with two scales: one with the magnetic fields as they were calculated with the formula and computed with TEMSI-FD, another with all the results in decibels (dB). Adding the two errors helps to prevent high deviations at low values. The algorithm searches for the values of the coefficients that lead to a local minimum of the sum. The search is conducted in a region defined by bound constraints, starting from the middle of the region. The lower bound for all the coefficients is set to zero. The upper bound is set to 0.01 for 𝑐 1 and 𝑐 04 , and to 2 for the other coefficients. The step tolerance is set to 10 -14 .

|𝐻| 𝑚𝑎𝑥 = 𝐼 0 × ( 𝑐 1 𝑘 × 𝑑 𝑐 2 𝑟 × 𝑑 𝑐 3 𝑒 + 𝑐 04 3 √ 𝑣 𝑑 𝑐 05 - √ 𝑎 20×ℎ 𝑟 × 𝑑 𝑐 6 ×|𝑐 7 -2×𝑑𝑟 ℎ | 𝑤 × 𝑒 𝑐 8 ×𝑑𝑛 ) + 𝜀 (3.

Estimated error

To estimate the error of the formula in (3.7), we consider all the points at the 11 horizontal surfaces inside each of the 10 buildings described in section 3. not quite precise. Altogether, the formula proposed is an acceptable alternative to full-wave simulations. Certainly, it is more complex than the IEC formula. Yet, it is still sufficiently straightforward for engineering applications.

Machine learning approach

To reduce the amount of computational resources without compromising the accuracy of the results, we could adopt a more sophisticated approach. The error of the formula in (3.7) seems reasonable until we compare it to the prediction error of ML algorithms.

To illustrate, let us consider a dataset built from the variables defined in (3.7) and the fields computed at the 11 horizontal surfaces in dB. We split the dataset into five groups to train and validate three simple ML algorithms: a feed-forward neural network, a fine regression tree, and a bagged tree. The prediction error is evaluated as the combination of the prediction errors over the five rounds of a five-fold cross- Even with a relatively small dataset, the superiority of the ML algorithms is clear.

Nevertheless, it does not necessarily mean that ML is the best way to estimate the magnetic field strength inside grid-like shields. The choice of the approach would depend on the indicators used to measure the performance of the models and the availability of the resources. A drawback perhaps of ML for industrial applications, is that the final user may require access to the training dataset and a specific system configuration to make predictions.

Correction factor

The distribution of the peak-values of the magnetic field strength inside a single-layer grid-like shield can be estimated with either a full-wave simulation, the IEC formula in (3.2), the improved formula in (3.7), or an ML algorithm. Now, we can move on to calculating a correction factor to adjust the estimated values.

Experimental design

To build the model to calculate a correction factor, we consider modifications to the geometrical configuration of the grid-like shields. The variables that are taken into of the mesh size. The latter means that when 𝑥 4 is equal to 0.5, the actual distance between the layers is equal to half the mesh size (0.5𝑥 1 ), and when 𝑥 5 is equal to 1, the distance between the hoops interconnecting the layers is equal to the mesh size (𝑥 1 ). One could define the distances differently and continue with the same line of reasoning.

A full factorial design for five variables would be suboptimal if we want to consider different grid-like shields. Thus, to limit the number of simulations, we define a Box-Behnken design with a single central point. A Box-Behnken design for five factors results in ten blocks with two factors varying simultaneously, for a total of 41 combinations, including the central point. In the design, the distance between the layers (𝑥 4 ) and between the hoops (𝑥 5 ) is set to 0 m when there is a single layer.

Consequently, four experiments become the same, which reduces the number of runs from 41 to 38.

The full-wave simulation of each of the 38 resulting runs is carried out using TEMSI-FD in shields B1 to B4 (see table 3.5), for a total of 152 simulations. The numerical model of the electromagnetic environment and the points of observation Distance between the layers as a multiple of the mesh size

𝑥 4 0.5 -2
Distance between the hoops as a multiple of the mesh size

𝑥 5 1 - 8
Note: The distances vary within the given ranges when there are at least two layers, otherwise they are set to zero.

remain as described in section 3.3.2. Then, for each run, the peak-values of the magnetic field strength computed with the modified shield are divided by the peakvalues of the magnetic field strength computed with the reference shield. Moreover, for 38 variations in the grid-like shield, we obtain the effect it produces on the peakvalues with respect to the reference shield on 44 surfaces, 11 by building. From the surfaces, we consider only the values at points inside the safety volume, which depends on the mesh size and the number of layers. We convert the values into dB and regroup them by run.

Figure 3-32 shows that even when regrouping the results by shield, the values by run vary considerably. Since the variations are hardly ever lower than 1 dB by run, we fit a normal distribution to the quantified effect. We use the means of the fitted distributions to build a sparse PCE and the standard deviations to estimate the uncertainty of the model.

Sparse PCE and global sensitivity analysis

Let us assume that the mean of the normal distribution fitted to the effect of the variations in the peak-values is the output random variable. If we consider the five factors of the Box-Behnken design as input random variables with a uniform distribution, we can build up the basis of a PCE from the Legendre polynomials 𝑃 𝛼 (𝑥). The family of Legendre polynomials is obtained when the input vector 𝑋 has a uniform

𝑥 𝑖 = 𝑎 𝑖 + 𝑏 𝑖 2 + 𝑏 𝑖 -𝑎 𝑖 2 × 𝜉 𝑖 , 𝜉 𝑖 ∼ 𝒰(-1, 1), (3.8)
where 𝑎 𝑖 and 𝑏 𝑖 are the lower and the upper limits, respectively, defined for each input variable in table 3.8. To transform 𝑥 4 and 𝑥 5 into standard uniform variables, their lower limit is set to 0 because the distances are non-existent when there is only one layer (𝑥 3 = 1). Thus, 𝜉 4 = 𝑥 4 -1 and 𝜉 5 = 2𝑥 5 -8 8 .

Accordingly, a multivariate polynomial of the basis reads

Ψ 𝛼 (𝜉 1 , . . . , 𝜉 𝑀 ) = 𝑀 ∏ 𝑖=1 √ 2𝛼 𝑖 + 1 × 𝑃 𝛼 𝑖 (𝜉 𝑖 ), (3.9)
where 𝑃 𝛼 𝑖 is the univariate Legendre polynomial of degree 𝛼 𝑖 and |𝛼| = ∑ 𝑀 𝑖=1 𝛼 𝑖 is the degree of the multivariate polynomial.

Since we only have 38 runs to compute the coefficients of the multivariate polynomials forming the basis, we could use a second-order standard truncation scheme or find a sparse expansion as proposed in [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF]. We continue with the last option.

Using UQLab [START_REF] Marelli | UQLab: A framework for uncertainty quantification in MATLAB[END_REF], we apply an adaptive algorithm based on least angle regression to identify the optimal basis from a fourth-order truncation scheme and compute the coefficients. The significant terms of the basis and their associated coefficients are listed in table 3.9. Moreover, the correction factor (𝐶𝐹 ) to adjust the peak-values of the magnetic field strength can be calculated as follows:

Ŷ = ℳ 𝑃 𝐶 (𝜉 1 , 𝜉 2 , 𝜉 3 , 𝜉 4 , 𝜉 5 ) = 17 ∑ 𝑗=0 𝑦 𝑗 × Ψ 𝑗 = 0.261 + . . . + 1.971 ( 5𝜉 3 1 × 𝜉 3 -3𝜉 1 × 𝜉 3 ) , (3.10) 𝐶𝐹 = 10 Ŷ 20 , (3.11)
In (3.10), the effect of the modifications on the geometrical configuration of gridlike shields is approximated by an 18-term polynomial series expansion. Note that the Table 3.9: Base and coefficients of the sparse PCE in (3.10)

𝑗 𝛼 𝑦 𝛼 ≡ 𝑦 𝑗 Ψ 𝛼 ≡ Ψ 𝑗 0 [0,0,0,0,0] 0.261 1 1 [0,0,0,1,0] -0.666 √ 3 × 𝜉 4 2 [0,0,1,0,0] -2.212 √ 3 × 𝜉 3 3 [0,1,0,0,0] -2.259 √ 3 × 𝜉 2 4 [1,0,0,0,0] 6.282 √ 3 × 𝜉 1 5 [2,0,0,0,0] -1.691 √ 5 2 (3𝜉 2 1 -1) 6 [1,1,0,0,0] 1.297 3𝜉 1 × 𝜉 2 7 [0,0,3,0,0] -3.379 √ 7 2 (5𝜉 3 3 -3𝜉 3 ) 8 [0,3,0,0,0] -0.466 √ 7 2 (5𝜉 3 2 -3𝜉 2 ) 9
[3,0,0,0,0] 3.077

√ 7 2 (5𝜉 3 1 -3𝜉 1 ) 10 [0,1,0,0,2] 0.369 √ 15 2 (3𝜉 2 × 𝜉 2 5 -𝜉 2 ) 11 [0,1,0,2,0] 0.240 √ 15 2 (3𝜉 2 × 𝜉 2 4 -𝜉 2 ) 12 [1,2,0,0,0] 0.281 √ 15 2 (3𝜉 1 × 𝜉 2 2 -𝜉 1 ) 13 [2,1,0,0,0] -0.836 √ 15 2 (3𝜉 2 1 × 𝜉 2 -𝜉 2 ) 14 [0,0,0,4,0] -0.900 3 8 (35𝜉 4 4 -30𝜉 2 4 + 3) 15 [4,0,0,0,0] -0.977 3 8 (35𝜉 4 1 -30𝜉 2 1 + 3) 16 [1,3,0,0,0] 0.385 √ 21 2 (5𝜉 1 × 𝜉 3 2 -3𝜉 1 × 𝜉 2 ) 17 [3,0,1,0,0] 0.860 √ 21 2 (5𝜉 3 1 × 𝜉 3 -3𝜉 1 × 𝜉 3 )
sparse PCE was built from the effect on the peak-values with respect to the reference shield. Nevertheless, one could use the correction factor to adjust the peak-values estimated with a different grid-like shield if its characteristics are within the ranges defined in table 3.8. The performance of the sparse PCE is shown in figure 3-33.

To perform a global sensitivity analysis, we can analytically calculate the Sobol' indices from the coefficients of the sparse PCE. The first-order and the total indices are shown in figure 3-34. As expected, the variance of the quantified effects is mostly influenced by the mesh size and the number of layers. And, in accordance with [START_REF] Metwally | Reduction of lightning-induced magnetic fields and voltages inside struck double-layer grid-like shields[END_REF] and the parametric study in section 3.2, the effect of the number of hoops interconnecting the layers is negligible. Definitely, the number of layers must be considered in the models to calculate the shielding effectiveness provided by reinforced concrete.

Disregarding an extra layer could lead to overestimating the magnetic field by a factor ------------------up to 6, which may result in costly protection measures.

Since the first-order indices account for most of the contributions to the total variance, we could reduce the number of terms of the sparse PCE by disregarding the effect of the interaction. If the contribution of each variable is considered separately, a fair estimation of the correction factor might still be obtained. Yet, considering the higher-order terms is essential in some cases. For example, figure 3-35 shows that when the mesh size decreases, the influence of the radius of the rebars becomes significant. On the other hand, the distance between the layers does not seem to affect the effect of adding an extra layer.

Uncertainty quantification

The sparse PCE to predict the effect on peak-values ( Ŷ ) was built from the means of the normal distributions fitted to the results of the 38 runs in the experimental design. Thus, even though the RMSE and the leave-one-out error (𝐸 𝐿𝑂𝑂 ) are small (see figure , there could be an important error on the 𝐶𝐹 calculated to adjust the values at an arbitrary point. Certainly, there is a local error in the PCE predictions associated with the deviation of each effect from the mean of the fitted distribution.

We induced that error by disregarding the dependence of the effect on the dimensions of the building and the position of the point of interest. In addition, there is an Then, to estimate the local error of the original PCE, we create 1000 input vectors, different from the ones considered in the experimental design. Each vector corresponds to an alternative modification of the grid-like shield. To create the input vectors, we use again a Monte Carlo method to draw 1000 samples from the input variables within the ranges defined in table 3.8, making sure that 𝑥 3 is an integer and

𝑥 4 = 𝑥 5 = 0 when 𝑥 3 = 1.
Generally, a few hundred samples are sufficient to estimate the basic shape of a PDF [START_REF] Proppe | Equivalent linearization and Monte Carlo simulation in stochastic dynamics[END_REF]. Again, although the sparse PCE was built from the effect on the peak-values with respect to the reference shield, it can be used to predict a correction factor to be applied to the peak-values estimated with a different grid-like shield. In that case, the error within the shield may increase because we would be calculating the effect on peak-values as the difference between two predictions: a prediction considering the characteristics of the shield we have and a prediction considering the characteristics of the shield of interest.

For example, to calculate the peak-values of the magnetic field strength inside the full-scale building described in [START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF], we could assume that the characteristics of its reinforcement correspond to the characteristics of the reference shield and apply the excluding the error of the formula in (3.7). In case B-5, the peak-value computed by the authors in [START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF] at point 𝐴 was 167 A/m. Note that their model included the columns, beams, and floors of the building.

Also, one could calculate a correction factor to adjust the available results, instead of applying the formula in (3.7). If we take the cases B1 and B3 considered in [START_REF] Maksimowicz | Investigation of models of grid-like shields subjected to lightning electromagnetic field: experiments in the frequency domain[END_REF]: a single-layer and a double-layer grid-like shields with meshes of 30 cm and rebars with a radius of 1.4 mm, the sparse PCE would predict 12.335 dB + 𝜀 ∼ 𝒩 (-0.164, 1.894 2 ) and 6.481 dB + 𝜀 ∼ 𝒩 (-0.164, 1.894 2 ), respectively. Thus, the peak-values of the magnetic field inside the single-layer grid-like shield could be adjusted to take into account a second layer by adding -5.854 dB + 𝜀 ∼ 𝒩 (0, 2.679 2 ), which is in good agreement with the authors' findings in [START_REF] Maksimowicz | Investigation of models of grid-like shields subjected to lightning electromagnetic field: experiments in the frequency domain[END_REF]. They observed that adding a layer enhances the magnetic shielding by 4.9 dB experimentally, and by 7.3 dB in full-wave simulations.

Similarly, if we take the cases AI and AIII considered in [START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF]: a single-layer and a double-layer grid-like shields with meshes of 15 cm and rebars with a radius of 3 mm, the sparse PCE will predict -0.053 dB + 𝜀 ∼ 𝒩 (-0.164, 1.894 2 ) and -12.021 dB + 𝜀 ∼ 𝒩 (-0.164, 1.894 2 ), respectively. Thus, one could expect the second layer to improve the shielding by 11.968 dB + 𝜀 ∼ 𝒩 (0, 2.679 2 ). The authors in [START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF] measured an improvement of 3-7 dB and later in [START_REF] Metwally | Reduction of lightning-induced magnetic fields and voltages inside struck double-layer grid-like shields[END_REF] computed an average factor of about 2.6 (∼ 8.3 dB). The difference could be explained by the fact that in both the experimental setup in [START_REF] Metwally | Measurement of magnetic fields inside single-and double-layer reinforced concrete buildings during simulated lightning currents[END_REF] and the computation model in [START_REF] Metwally | Reduction of lightning-induced magnetic fields and voltages inside struck double-layer grid-like shields[END_REF], the ground was a highly conducting plane.

To conclude, there are different alternatives to calculate the peak-values of the magnetic field strength generated by a direct lightning strike inside grid-like shields.

In this chapter, we proposed a variant of the IEC formula to calculate the distribution of the peak-values and a sparse PCE with 18 terms to account for modifications made to the geometrical configuration of the shield. The sparse PCE can be used to adjust the results of either the improved formula or another computational method. The modifications to the configuration of the shield include changes to the mesh size, the radius of the rebars, the number of layers, the distance between the layers, and the distance between the hoops.

Chapter 4

Electrical surges in underground cables due to a direct lightning strike

Now that we have different alternatives to calculate the peak-values of the magnetic field strength, we can move on to calculating the electrical surges in underground cables. We still do not know how the different components of the electromagnetic environment affect the surges; thus, we need to go back to the basis of our 3-D model and conduct a parametric study to evaluate the possibility of reducing its complexity.

Then, as we did for the magnetic field, we will be able to collect data from the simulations and start building surrogate models.

Parametric study

Reference case

First, to study the surges resulting from a direct lightning strike in an industrial facility, a simplified case is defined: two reinforced concrete buildings facing each other, with a height of 20.25 m above the ground and a horizontal cross-sectional area of 20.25 m × 50.25 m. The dimensions include the thickness of the concrete.

The external walls, the roof, and the foundation of the buildings are made up of a single-layered reinforcing grid. No internal walls, columns, or beams are considered.

The reinforcing grid has squared meshes of 50 cm × 50 cm and is 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 at a depth of 5 m.

As shown in figure 4 The equipotential bonding network inside the buildings consists of collective ground As in chapter 3, the lightning channel is modeled as a monopole antenna by a vertical lossy wire, excited at its base by a lumped current source and connected at the top end to a PML [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. The channel is attached to the center of the roof. Aiming to conduct the study in the frequency domain, the current waveform is also defined as a Gaussian function covering a frequency range from 0 Hz to 1 MHz.

Using TEMSI-FD [START_REF]TEMSI-FD: Time Electromagnetic Simulator -Finite Difference Time Domain[END_REF], a non-uniform grid is implemented with cells varying from We saw the advantages of applying the MPM in section 2.2.2. To reduce the computation time, it is also preferable to continue using a Gaussian pulse. Nevertheless, it is important to note that calculating the currents in the time domain from the currents computed using a Gaussian pulse is not always straightforward. There are different parameters, such as the number of poles, that can influence the accuracy of the extrapolation. The process generally leads to comparable results; yet, when verifying, we did sometimes observe differences in the decay time (see e.g. figure 456).

Even in a simplified case such as the reference case, a small variation in the geometry could intensify the surges resulting from a direct lightning strike. Thus, we decided to divide the parametric study into two parts. First, we focus on parameters that are rarely mentioned in lightning protection standards to define a worst-case scenario. Then, the worst-case scenario becomes the new reference case, and we conduct a more generalized study.

First part: Routing, bonding, and grounding

Electrical surges are induced in cables by the transient electromagnetic field and the ground potential rise generated by a lightning strike. Nevertheless, during a direct strike, surges can also be induced by the coupling of the cables to different structures, and by a partial lightning current entering the conductors via their grounding. Hence, the phenomenon is strongly dependent on the impedance of the lines, the conductor routing, and the connections to the bonding network.

Since there are hundreds of cables following different paths in an industrial facility, to estimate the lightning-induced surges, the cables are often considered to be grounded at the entrance to the buildings. By doing so, the uncertainty associated with the cable routing and the bonding network is disregarded. Therefore, the results of the numerical model can be applied to a larger set of arrangements in the facility.

If the conductors in the reference case are connected to the reinforcement at the entrance to the buildings, instead of being connected to the bonding rings, the currents rise significantly. In addition, the resonances are no longer present below 1

MHz. Figure 4-7 shows that the current increases in the insulated conductor, even if it is only the bare conductor that is grounded at the entrance. Moreover, grounding the conductors at the entrance to the buildings appears to be an important overestimation, and one wonders: considering just the influence of the cable routing and the configuration of the earth-termination system, would those currents be representative of a worst-case scenario?

To address the question, we conduct a study on seven parameters: the routing inside the buildings, the configuration of the bonding network, the configuration of the grounding network, the soil resistivity, the characteristics of concrete, the striking point, and the position of the duct. The parameters are varied independently, and their effect is analyzed at low frequency and at the frequency of the first resonance.

Routing inside the buildings

As we saw in the previous chapter, the transient electromagnetic field generated by a direct lightning strike inside the building varies with the distance to the walls, the roof, and the foundation. Thus, the surges are expected to vary if the cable routing is changed. To confirm, first, we tried moving the conductors closer to the walls but keeping them connected to the same bonding ring. Then, we tried connecting the conductors to different bonding rings. Figures 4-8 and 4-9 show the results of the following cases:

Case A: Inside B1, the conductors are moved 4 m closer to the wall on the 𝑥-axis.

Case B: Inside B1, the conductors are moved 5 m further away from the wall on the 𝑥-axis.

Case C: Inside B2, the conductors are moved 5 m further away from the wall on the 𝑥-axis.

Case D: Inside B1, the conductors are connected to the bonding ring on the third floor.

Case E: Inside B1, the conductors are connected to the bonding ring on the first floor.

Case F: Inside B2, the conductors are connected to the bonding ring on the third floor.

Case G: Inside B2, the conductors are connected to the bonding ring on the first floor. The figures show that significant variations are only observed when connecting the conductors to different bonding rings. As long as the conductors are relatively far from the walls, the effect of changing the routing on the horizontal plane is negligible.

On the other hand, when the connection to the bonding network is moved away from the roof, the effect is noticeable at both low frequency and at the frequency of the first resonance. At low frequency, the induced currents increase when the length of the conductors and thus their impedance is reduced. The highest currents are observed when the conductors are connected to the bonding ring on the first floor of B1. At the frequency of the first resonance, the currents decrease when the connection inside B1 is close to the foundation, but increase when the connection inside B2 is close to the foundation. The highest currents are observed when the conductors are connected to the bonding ring on the first floor of B2.

Configuration of the bonding network

Since the bonding rings are interconnected in the corners and the lightning channel is attached to the center of the roof, the effect of changing the routing is probably reduced. Moreover, we tried interconnecting the bonding rings only in one corner and modifying the number of connections from the bonding network to the reinforcement of the building. By forcing the current to take a single path to flow to earth, the position of the connection of the conductors to the bonding network gains influence. It can be observed that when the bonding rings are connected once, rather than at each corner, the frequency of the first resonance is shifted because the length of the path to earth is changed. Additionally, since the currents collected converge to one point, the coupling between the rings and the conductors is stronger in some cases. Reducing the number of connections between the bonding network and the reinforcement of the buildings has almost no effect other than shifting the frequency of the first resonance. An important decrease is observed when there are no connections, but there is always at least one to prevent potential differences. Depending on the part of the facility, the mesh size of the grounding network can be enlarged, or even limited to a single connection between the ground rings.

The smaller the mesh size, the lower the impedance of the network and the potential differences between buildings. Hence, the highest currents are expected when there is not a meshed grounding network, only the ground rings with a single connection between them. Figure 4-14 shows that the differences are significant. An increase in the earth impedance may lead to a higher current flowing through the duct, and therefore, to a stronger coupling to the conductors routing inside.

Soil resistivity

The earth impedance depends on both the configuration of the earth-termination system and the soil resistivity. The higher the resistivity (𝜌), the higher the earth impedance. Thus, as expected, figure 4-15 shows that the currents induced in the conductors also increase when the resistivity increases. material, the phenomenon could perhaps be explained by the dielectric losses that would be disregarded in the model when the reinforcement above ground is surrounded by air.

Striking point

As stated previously, the symmetry of the configuration and the attachment of the channel to the center of the roof can reduce the effect of the routing. It can be observed in figure 4-18 that if the building is struck in corner C11 or C12, the currents rise significantly at low frequency. The rise is expected because C11 is the closest corner to the interconnection of the conductors and the bonding system and C12 is one of the closest corners to the duct. Nevertheless, it appears that the symmetry strengthens the resonances because at the frequency of the first resonance the highest currents are observed in the reference case.

Position of the duct

When a building is struck by lightning, the highest currents flow to earth through the rebars located in the vertical edges. One could assume that those partial currents are equal if the lightning channel is attached to the center of the roof. However, the distribution of the currents in the reinforcement is affected by its connection to other conductive components in the facility, such as the reinforced concrete ducts.

Similarly, the portion of the lightning current that flows through a duct is affected by its distance to the vertical edges. Figure 4-19 shows the results of the following cases:

Case A: The distance between the duct and the vertical edges in corners C12 and C21 is reduced to 1 m on the 𝑦-axis.

Case B: The distance between the duct and the vertical edges in corners C12 and C21 is reduced to 12 m on the 𝑦-axis.

Case C: The distance between the duct and the vertical edges in corners C12 and C21 is reduced to 1 m on the 𝑦-axis and the lightning channel is attached to corner C12.

Case D: The distance between the duct and the vertical edges in corners C12 and C21 is reduced to 12 m on the 𝑦-axis and the lightning channel is attached to corner C12.

It can be observed that the closer the duct is to the edges and to the striking point, the higher the currents. Additionally, changing the position of the duct also changes the path the conductors follow inside the building; thus, at low frequency, the highest currents are obtained in case C because of the proximity of segment S4

Worst-case scenarios

To summarize, the magnitude of the surges generated by a direct lightning strike depends on the configuration of the industrial facility. Aiming to define a worst-case scenario, a study was conducted on parameters associated with the cable routing and the configuration of the earth-termination system. The idea is to compare the currents computed in the worst-case scenario to the currents computed when the conductors are grounded at the entrance to the buildings.

The parametric study showed that it is essential to be as precise as possible when it comes to modeling the concrete and the number of connections between the ground rings. Adding concrete between the reinforcement and the surface leads to an important attenuation of the surges at low frequency. A grounding network with a bigger mesh size or a single connection between the ground rings leads to a higher earth impedance.

Regarding the influence of the other parameters, their effect depends on the position of different elements. Table 4.1 shows the cases in which the highest currents were obtained. Based on these results, two worst-case scenarios are defined, the first to maximize the induced currents for slow excitations and the second to intensify the resonances. The differences between the reference case and the scenarios are described in table 4.2. The currents computed in the worst-case scenarios are shown in figure 4-20. As expected, they are higher than the currents computed in the reference case. Yet, it can be observed that the increase obtained when grounding the conductors at the entrance to the buildings no longer appears to be an important overestimation. At low frequency, the currents computed with or without the conductors grounded at the entrance are comparable. Hence, as long as the grounding network is modeled Table 4.2: Differences between the reference case and the worst-case scenarios.

First scenario Second scenario

Routing

The conductors are connected to the bonding ring on the first floor, 1 m away from the reinforcement of the buildings. Inside B2, the connection is close to C21 to reduce the cable length.

Inside B1, the connection is close to C11 to strengthen the coupling between the conductors and the bonding network.

Inside B1, the connection is close to C12.

Bonding

Inside B1, the bonding rings are interconnected in corner C12.

network The bonding network is connected to the reinforcement of B1 in C11.

The bonding network is connected to the reinforcement of B1 in C12.

Grounding network

There is a single connection between the ground rings instead of a meshed grounding network.

The reinforcement of B1 is connected to the ground ring in C11.

The reinforcement of B1 to the ground ring in C12.

Striking point

Corner C12.

Position of the duct

The distance to the edge is reduced to 1 m.

appropriately, one could get a good idea of the currents expected in the worst-case scenarios by grounding the conductors at the entrance to the buildings. It would reduce the complexity of the model and the number of parameters to be taken into account to estimate the surges. For example, figure 4-21 shows that when the conductors are grounded at the entrance, the number and the position of connections between the ground rings and the reinforcement of buildings have a negligible effect on the results.

However, when the conductors are grounded at the entrance to the buildings, the first resonance is shifted towards a higher frequency and its magnitude is reduced.

Thus, the resonant behavior of the currents is expected to decrease when computing the time response for a fast waveform. To confirm, the current in the worst-case scenarios are calculated for the positive and the subsequent strokes (see figures 4-22 and 4-23). Note that to calculate the currents for the subsequent stroke, the simulations had to be repeated using a Gaussian pulse covering a 10 MHz bandwidth.

The waveforms of the currents calculated for the positive stroke are in good agreement, although the peak-values are higher when the conductors are grounded at the entrance to the buildings. On the contrary, the peak-values calculated for the subsequent stroke are in the same range, but the waveforms are quite different. Since this work is focused on estimating the peak-values of the surges resulting from a direct lightning strike, we decided to continue the study with the conductors grounded at the entrance to the buildings and the configuration of the first worst-case scenario.

This scenario is hereafter referred to as the new reference case.

Second part: Geometry of the structures and characteristics of the components

From the first part of the study, we now have one least complex model to estimate the peak-values of the surges in a worst-case scenario, and three parameters that we may still need to consider: the soil resistivity, the thickness of the concrete, and the position of the duct. An upper view of the model is shown in figure 4-24. Note that segments S2 to S4 and S6 to S8 no longer exist. Segment S5 is still inside the duct, and segments S1 and S9 are the ones connected to the reinforcement of the buildings.

We observed that the position of the duct had an important effect on the surges in the reference case, mainly because changing the position of the duct also changes the path the conductors follow inside the buildings. However, if the conductors are grounded at the entrance to the buildings in the new reference case, does it still have an effect? We repeated the simulations to verify. Figure 4-25 shows the results of the following cases:

Case A: The distance between the duct and the vertical edges in corners C12 and C21 is enlarged. The duct is now at 1 m from the edges in corners C14 and C23.

Case B: The duct is in the middle of B1. The distance between the duct and the vertical edge in corner C21 is still 1 m.

Case C: The duct is in the middle of B2. The distance between the duct and the vertical edge in corner C12 is still 1 m.

It can be observed that the position of the duct with respect to the edges of B2 is not relevant. On the other hand, the position of the duct with respect to the edges of B1 still affects the results. Also, the currents depend on more than the distance to the edges. In the new reference case, where the duct is at 1 m from the edge in corner C12, the currents increase faster with the frequency than in case A. Yet, in case A the currents are higher at low frequency. It could be explained by the fact that in case A the duct is closer to the electrode interconnecting the ground rings.

To confirm, we define three additional cases, where the lightning channel is attached shorter distance between the second duct and the interconnection leads to a stronger coupling, and probably, a higher current passes by the duct to head to the electrode.

But perhaps, a lower current actually flows through the second duct because there is an easier path to earth nearby. Anyhow, since the second duct provides another path for the current to flow, it makes sense to observe lower currents in cases 𝐵 and 𝐶.

We must consider the number of ducts and their position with respect to the edges and striking point. Moreover, observing variations between cases 𝐵 and 𝐶 and between the ducts is a little unexpected. One reason could be that when there are three buildings, a higher current flows through the ducts because the earth impedance is lower, and therefore, there is a stronger coupling to the conductors. In that case, instead of modeling the third building, one could increase the cross-sectional area of B2. To continue with the study, we now go back to the new reference case. In this second part, we focus on parameters associated with the geometry of the structures and the characteristics of the components: the dimensions of the buildings, the distance between the buildings, the number of connections between the duct and the buildings, the geometrical configuration of the reinforcement, the load at the end of the insulated conductor, and the types of conductors. It can be observed that as long as the cross-sectional area remains constant, the variations are small. It is of course expected to observe slightly higher currents if the dimensions of B1 are reduced on the 𝑦-axis because, in that case, the duct is closer to the rebars in the edges. Regarding the height of the buildings and the depth of the foundations, it is interesting to see that they both play an important role.

Dimensions of the buildings

The currents are lower at low frequency if the foundations are buried deeper into the ground. Perhaps because the earthing is better, but it could also be because the foundations are just further away from the duct and therefore, there is less coupling.

The currents are higher if B1 is shorter, probably because the duct is closer to the lightning channel. And the currents are also higher if B2 is taller, most likely because there are more rebars exposed to the electromagnetic field.

Distance between the buildings

A larger distance between the buildings means longer conductors. Thus, as expected, figure 4-30 shows that the frequency of the first resonance is shifted when the distance is changed. Also, the longer the conductors, the higher their impedance and therefore, the lower the current entering via their grounding. That explains why the highest current in the bare conductor is observed when the distance is reduced to 50 m. On the other hand, since the insulated conductor is loaded, the current in the conductor increases with its length. Probably, because the coupling between the conductor and the duct goes on. However, the current steepness in the duct is known to decrease with the distance to the building; hence, after a certain distance, the current induced in the conductor should stop increasing.

Cables

Again, the partial lightning current entering a cable depends on its impedance, which depends on the characteristics and parameters of the cable. Yet, the characteristics and parameters of the cable also influence its electromagnetic coupling to other conductors. For example, increasing the section of the bare conductor decreases its impedance. Thus, as shown in figure 4-37, the current increases in the bare conductor and decreases slightly in the insulated conductor. Nevertheless, the effect is minor.

Similarly, when increasing the number of cores of the insulated conductor, the current increases in the insulated conductor and decreases slightly in the bare conductor. Note that to obtain the results shown in figure 4-38, we modeled a single-core equivalent conductor. The radius of the equivalent conductors was calculated using LAMLIM [START_REF]LAMLIM[END_REF] and considering the arrangements shown in figure 4-39 over a perfectly conducting ground plane.

It can be observed in figure 4-38b that duplicating the number of cores leads to

Parameters to be considered in the numerical model of an industrial facility

Since it is nearly impossible and computationally expensive to consider all the elements of a full-scale industrial facility in full-wave simulations, we started by conducting a parametric study on a simplified case. The study helped us identify some of the parameters that could be disregarded and the parameters that must be considered in the numerical models.

In the first part, we studied the influence of the cable routing inside the buildings, considering the distance to the walls and to the roof. We varied the number of connections between the bonding rings, then between the bonding network and the reinforcement of the buildings, and finally, between the reinforcement of the buildings and the ground rings. We also studied the influence of the earth impedance by changing the type of grounding network, its depth, and the soil resistivity. We observed, once again, that the effect of the moisture content of the concrete is minor; yet, it is essential to be as precise as possible when it comes to its thickness. And, we saw the importance of the position of the duct with respect to the vertical edges and the striking point.

To reduce the complexity of the model, we decided to continue the study with the conductors grounded at the entrance to the buildings. The results obtained in that case give a good idea of the peak-values one can expect in a worst-case scenario. Moreover, it was no longer necessary to consider the cable routing, the bonding network, and the number of connections between the reinforcement of the buildings and ground rings. Also, to reduce the number of parameters and stay on the safe side, we interconnected the ground rings with a single electrode. The mesh size of the grounding network is therefore no longer a parameter to be considered.

In the second part, we repeated some simulations varying the position of the duct, and we included a second duct and a third building to corroborate the results. We saw that the surges are affected by both the dimensions of the buildings and the depth of the foundations. Surprisingly, we observed that the geometrical configuration of the reinforcement of the duct has a significant effect on the surges, whereas only the mesh size of the reinforcement of the buildings seems to be important. We also observed that there is a great dependency on the terminations of the power installation and the type of cable considered.

In short, to obtain a reasonable approximation of the peak-values of the surges resulting from a direct lightning strike, we must consider at least the influence of:

• The same tendency we observed in section 3.1.2 can be observed in figures 4-40 and 4-41. Both the current and the magnetic field strength increase when the radius increases and the mesh size decreases. In addition, as expected, they increase with the conductivity of the rebars. Overall, since we defined a voltage source to characterize the grid in the reduced representation, the current and the magnetic field strength increase when the impedance of the grid decreases. It is almost imperceptible in the figures but the magnetic field strength also has a longer decay time.

Again, we take 65 uniformly spaced samples of the time responses, and we estimate the parameters of the Kriging models using a HGA with a convergence tolerance of 10 -10 and a maximum number of generations equal to 500. The values of the parameters obtained with a Metern 5/2 kernel are given in table 4.4. We did not try to adjust different regression models to the data because we have already obtained good results with a Kriging model and calibrated the algorithms in section 3.1. The Kriging models are implemented in a new fitting algorithm, which also follows the procedure described in figure 3456789. The only difference is that in the new fitting algorithm, the input vector has three components instead of two. In France, most of the rebars used to build the reinforcement of the ducts have a radius of 6 mm or 7 mm. Since we saw in section 4.1.3 that the surges in the cables decrease with the radius of the rebars, we could assume that they all have the same radius and that it is equal to 6 mm. Let us then consider a double-layered grid with meshes of 25 cm × 25 cm, formed by rebars with a radius of 6 mm. The layers are at a distance of 25 cm, and they are interconnected every 50 cm. Using the new fitting algorithm, we would obtain an equivalent single-layered grid with a mesh size of 0.427 m and rebars with a radius of 2.897 cm and a conductivity of 1.06 ×10 6 S/m. However, in FDTD simulations it is preferable to define the wires in the middle of the FDTD cells. Since we had implemented cells of 25 cm, we rerun the algorithm setting the mesh size to 25 cm, and then again, setting the mesh size to 50 cm. With a mesh size of 25 cm, we obtained rebars with a radius of 2.44 cm and a conductivity On the other hand, none of the equivalent grids is able to match the amplitude of the magnetic field strength.

To check if the differences in the magnetic field strength have an important effect, we went back to the new reference case defined in section 4.1.2. We conducted three simulations: First, we considered that the reinforcement of the duct is made of the double-layered grid described above. Then, we replaced the double-layered grid with the single-layered equivalent grids. The currents computed in the conductors are shown in figure 4-43. Clearly, it is important to make sure that the equivalent grid also reproduces the magnetic field strength. However, if one has to choose between setting the mesh size and a small fitting error, it would be advisable to set the mesh size. In that case, to stay on the safe side, the single-layered grid with a mesh size of 50 cm is the best option.

Experimental characterization

We were trying to find a simplified 3-D model of the reinforcement; yet, we do not even know if the reinforcement we are considering as the basis represents appropriately the geometrical configuration of actual constructions.

The impedance of the reinforcement could be affected by the type of connection between the rebars. They could be welded or wire-tied at the intersections or a combination of both. What we know for sure is that good contact cannot be guaranteed at every intersection. Besides, in actual constructions, there is uncertainty regarding the mesh size and the distance between the layers of the reinforcing grid. Nevertheless, in numerical simulations, we generally consider a perfect contact between the rebars and symmetrical reinforcement. Thus, one question remains: how big of an error could we be making in our models? To corroborate some of our assumptions, we computed the impedance of the ducts at the new lightning testing facility and compared the results to the impedance measured on different sections of the ducts.

Lightning testing facility at EDF Power Networks Lab

Impedance measurements

To measure the impedance, we used a vector network analyzer (VNA) and two different methods: One-Port reflection, also known as Shunt measurement, and Shunt-Thru. The One-Port method offers a simple approach to measure the properties of the device under test (DUT), but its sensitivity is limited at high or low impedance values [START_REF] Schuster | Impedance Measurements using the Bode 100[END_REF]. Since the Shunt-Thru method is the one recommended to measure low impedance, we used both methods to compare the results. in the literature (see e.g. [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]). Even if the covering plates are not electrically connected to the reinforcement, they could affect the results.

Regardless of the method and the section of the duct being characterized, the length of the wires used to connect the VNA to the copper plates was constant.

Their impedance was measured independently to remove their contribution from the results. By remove, we mean subtracting the impedance of the wires from the total impedance measured. Also, we repeated all the measurements with the One-Port As expected, the reactance measured is predominantly inductive and does not vary significantly between the ducts. In the inductive region, we can fit the curves with an inductance of approximately 35 𝜇H. However, there are important differences between the ducts at low frequency before the inductance starts to dominate. On average, the resistance of the reinforcement with wire-tied rebars is 302 mΩ, whereas the resistance of the reinforcement with welded rebars is around 16 mΩ, and it might even be lower.

It is clear from the results in figure 4-53b that there is an error at low frequency. The There are also important differences between two sections of the same duct. As shown in table 4.5, the resistance and the impedance calculated from the results vary considerably. Yet, by adding the values of the resistance and the inductance of the sections, we come close to the values calculated with all the sections interconnected.

Note that the values of the resistance and the inductance were adjusted to the average of the impedance measured with the three methods.

It was expected to observe differences between two sections of the same duct because, after all, they are not equally long. Nevertheless, the broad range of variation of the resistance shows that it is difficult to control the quality of the contact between the rebars and therefore, to estimate the impedance of a duct based on its characteristics.

Numerical model

Moving on to the numerical model, several simulations are carried out using TEMSI-FD to calculate the impedance of the ducts at the new lightning testing facility. A uniform grid is implemented with cells of 5 cm, and all the conductors are modeled using, as usual, the thin-wire formalism. The conductivity of the rebars forming the The impedance of the duct was calculated as the voltage of a perfect voltage source over the current injected into the reinforcement, minus the impedance of the return path. However, we observed that even if we try to remove the effect of the return path, the reactance varies slightly with the position of the wires, which might also happen when measuring. The voltage waveform is defined as a Gaussian function covering a frequency range up to 10 MHz. The maximum amplitude of the pulse is set to 100 V. The boundary conditions of the analysis space are defined as PMLs to assume an open space.

The impedance calculated with a return path at 6 m and 8 m from the duct is shown in figure 4-55. Since 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. Surprisingly, an inductance of approximately 35 𝜇H can also be adjusted to the curves in figure 4-55. Also, the fact that the rebars may not always form perfectly symmetrical grids seems to have a minor effect on the results. Hence, we could at least assume that the reactance of a reinforced concrete duct can be well approximated considering the average dimensions of the meshes in full-wave simulations. As expected, the resistance calculated is of course lower because the contact resistance is neglected. We could try to adjust the resistance by, for example, reducing the conductivity of the rebars. The problem is that it would limit the partial lightning current flowing through the duct. If we are not able to measure the impedance of the duct of interest, it would be safer to consider the conductivity of steel. Altogether, disregarding the effect of the non-linearities that a poor contact between the rebars could create, our numerical models are a good representation of actual reinforced concrete ducts.

Estimation of the peak-values of the current in underground cables

We know the parameters that must be considered in the numerical model of an industrial facility (see section 4.2), and we have a simplified model of reinforced concrete ducts (see section 4.3). Now, we can move on to building a surrogate model to estimate the surges. To do so, we need a dataset.

The first step to build a dataset is to define the range of variation of the parameters:

• Since the soil texture in France is mostly loam [128], we consider its resistivity between 100 Ω•m and 1000 Ω•m.

• The distance from the duct to the closest edge and the striking point depends on the dimensions of the buildings. Typically, industrial buildings have a length between 10 m and 100 m, and a width between 5 m and 45 m, which brings us to a cross-sectional area between 50 m 2 and 4500 m 2 .

• Even though a duct can be divided into several ducts that go to different buildings along the way, we consider a single and independent duct connecting the building struck by lightning to each of the other buildings. To limit the volume of the analysis space, we set the number of buildings to two: the building struck by lightning (B1) and a building in front of it (B2). Therefore, we consider only one duct. Also, to facilitate the conception of the numerical model, we consider that the duct is straight.

• For each duct, we would have to define several parameters depending on its position with respect to the vertical edges and the striking point. To simplify, we consider a strike in the corner of the roof and that the duct is closest to the edge that happens to be the one carrying the highest lightning current. Hence, we consider only the relative distance, i.e., the distance divided by the length of the building, to the edge in the striking point.

• Again, to limit the volume of the analysis space, we consider that the distance between the buildings is equal to the length of the ducts and that it is 100 m at the most.

• The mesh size of the reinforcement of the buildings is rarely bigger than 1 m.

However, since a good contact cannot be guaranteed at every intersection, we need to consider meshes up to 5 m.

• Instead of defining the type of cable as a parameter, we consider three cables inside the duct: a bare conductor with a radius of 3.99 mm, an insulated conductor with an overall radius of 11.24 mm, and a coaxial cable with an overall Thevenin or a Norton equivalent circuit [START_REF] Powell | Duals and analogues[END_REF]. All the conductors are grounded at the entrance to the buildings.

To summarize, table 5.5 presents the parameters considered and their range of variation. The other parameters are set to fixed values. The considerations worth noting are the following:

• In an industrial facility, there is a considerable diversity of structures, with heights varying between 3 m and 100 m. However, the height of the buildings is typically between 10 m and 40 m above ground. Since we saw in section 4.2 that the surges increase when reducing the height of the building struck by lightning and the depth of its foundation, we set the height of all the buildings to 10 m and the depth of the foundations to 5 m.

• Regarding the reinforcement of the buildings, we consider a single-layered grid formed by rebars with a radius of 1 cm. We also saw in section 4.2 that the effect of the number of layers and the radius of the rebars is negligible. No internal walls, columns, or beams are considered.

• We model the ducts using the simplified model with a mesh size of 50 cm introduced in section 4.3. The reinforcement of the ducts is connected to the reinforcement of the buildings twice at each end.

• The concrete, as in section 4.1, is modeled as a lossy dielectric material with a conductivity of 0.0052 S/m and a relative permittivity of 8.6. The insulation of the conductors is also modeled as a lossy dielectric material with a conductivity of 1 ×10 -14 S/m and a relative permittivity of 2.9. The conductivity of steel and copper is set to 8.33 ×10 6 S/m and 59.6 ×10 6 S/m, respectively. The relative permittivity of the soil is set to 10.

• We do not consider a meshed grounding system but interconnected ground rings, as in the second part of the parametric study (section 4.2).

Finally, the dataset is built based on a 2 7 full factorial design, not considering the load at the end of the insulated cable as a parameter but as a different case. For each case, 128 simulations are therefore carried out using TEMSI-FD. A total of 256 simulations are carried out to build the surrogate models. The outputs we will try to estimate are the currents in the middle of each cable.

The lightning waveform is defined as a Gaussian function covering a frequency range from 0 Hz to 10 MHz, and a non-uniform grid is implemented with cells varying from 50 cm to 5 m. As in section 4.1.2, the normalized results obtained with the Gaussian function are then used to calculate the currents for the first and the subsequent strokes. The simulations are run up to 100 𝜇s. It is interesting to observe that the highest currents were computed in run 89 and run 121, whereas the lowest currents were computed in run 24 and run 56. The only parameter that changes between runs 89 and 121 and between runs 24 and 56 is the distance between the buildings. The distance is 10 m in runs 24 and 89, and it is 100 m in runs 56 and 121. In all four runs, the duct is in the middle of B1 (a relative distance equal to 0.5). Other than that, the parameters are the complete opposite.

The highest currents were computed when the soil resistivity is 1000 Ω•m, the length of B1 is 10 m, the width of B1 is 5 m, the cross-sectional area of B2 is 45000 m 2 , and the mesh size of the reinforcing grid is 5 m. The lowest currents were computed when the soil resistivity is 100 Ω•m, the length of B1 is 100 m, the width of B1 is 45 m, the cross-sectional area of B2 is 50 m 2 , and the mesh size of the reinforcing grid is 1 m. Also, as expected, the current in the bare conductor and the coaxial cable does not vary significantly between the two cases. Therefore, we will only consider the results of the case in which the insulated conductor is in OC. The average of the peak-values in the bare conductor and the coaxial cable is slightly higher when the insulated conductor is in OC (see table 4.7). In addition, since the currents in both cables are significantly higher for the first stroke, we will not consider the currents for the subsequent stroke. On the contrary, the currents in the insulated conductor are comparable for the first and the subsequent stokes; thus, we have to consider both.

Moreover, we need to build surrogate models to estimate six different outputs:

𝑌 1 :
The current in the bare conductor for the first stroke, when the insulated conductor is in OC.

𝑌 2 : The current in the insulated conductor for the first stroke, when it is in SC.

𝑌 3 : The current in the insulated conductor for the subsequent stroke, when it is in SC.

𝑌 4 : The current in the insulated conductor for the first stroke, when it is in OC.

𝑌 5 : The current in the insulated conductor for the subsequent stroke, when it is in OC.

𝑌 6 : The current in the coaxial for the first stroke, when the insulated conductor is in OC.

Since we have a small dataset, we decided to focus on building a sparse PCE and a Kriging model for each output. The sparse PCEs are built from a fifth-order truncation scheme, but surprisingly, all the significant terms of the basis are of degree equal or less than two. To build the Kriging models, a linear trend and a Matern 5/2 kernel were chosen. Because of the large number of coefficients that had to be computed, we will not list them in detail. The most important information about the models is presented in tables 4.8 and 4.9. As an example, consider the first worst-case scenario introduced in section 4.1.2.

The soil resistivity is 100 Ω•m, the length of B1 is 50 m, the width of B1 is 20 m, the cross-sectional area of B2 is 1000 m 2 , and the length of the duct is 100 m. However, we are out of the range of variation of the other parameters. The closest we could get to the scenario is by considering the mesh size equal to 1 m and the duct at 0.1 × 50 m = 5 m from the striking point. In that case, we would obtain the currents in table 4.10. The models lead to different results; thus, one has the option of selecting the highest current or averaging the values. In any case, it is interesting to observe that the peak-value of the current estimated in the bare conductor corresponds to the peakvalue of the current in figure 4-23b. The bare conductor considered in the models has a radius of about half the radius of the bare conductor in the worst-case scenario introduced in section 4.1.2; yet, the mesh size of the reinforcing grid is bigger, and therefore, the currents are comparable. where 𝐼 𝑂𝐶 is the current calculated when the insulated conductor is in OC, and 𝐼 𝑂𝑆 is the current calculated when the insulated conductor is in OS. Hence, for the first stroke, we would have a Norton resistance of 110.781 kΩ and 104.266 kΩ, for the PCE and the Kriging models, respectively. Consequentially, a current going through a 50

Ω load connected at the end of the conductor between 64.224 A and 61.51 A. The current is about three times higher than the current in figure 4-23a, which is expected because the radius of the insulated conductor is also about three times bigger.

The models already make good predictions. Nevertheless, a way to reduce the error would be by using different techniques to build additional surrogate models and then average the results.

Chapter 5

Conclusions and perspectives

Lightning strikes can have severe consequences in industrial facilities. Reliable methods to estimate their effects are not only necessary to design the LPS but also to demonstrate the safety of the installations. In France, lightning must be considered in the safety demonstration of NPPs. In this work, we focused on proposing alternative models adapted to the particularities of French industrial facilities to estimate the effects of a direct lightning strike to one of the buildings in the facility. Due to time constraints, we limited the scope of our analysis to the transient magnetic field inside the building struck by lightning and the surges induced in the cables routing in reinforced concrete cable ducts.

To deepen the understanding of the phenomena, we started by conceiving a generalized 3-D model of the electromagnetic environment. However, we quickly decided to study the effects separately because variations in the configuration of the facility do not influence the magnetic field and the electrical surges in the same manner. Hence, we conducted several parametric studies that led to two different and less complex 3-D models.

While trying to reduce the complexity of the 3-D model to compute the transient magnetic field, we attempted to find a simplified representation of the reinforcement of the building. We saw that it is possible to replace a multi-layered reinforcing with an equivalent single-layered grid; yet, the characteristics of the equivalent grid can be adjusted to reproduce only one of the features: the peak-values of the magnetic field or the impedance of the grid.

With the less complex 3-D models, we carried out full-wave simulations using TEMSI-FD, a solver based on the finite-difference time-domain (FDTD) method.

Based on the results, we introduced an improved formula as a variant of the formula suggested by the IEC 62305-4 standard. Both formulas can be used to determine the distribution of the peak-values of the magnetic field strength inside a singlelayer grid-like shield. The new variant is more complex; yet, more accurate and still sufficiently straightforward for engineering applications. Nevertheless, a comparison of its error to the prediction error of several ML models showed that there is still room for improvement.

To complement the improved formula, we built a sparse PCE to predict the effect of any modification on the geometrical configuration of the grid-like shield. The PCE allows for extrapolation from the results of either the improved formula or another computational method. It can be used to account for variations in the mesh size, the radius of the rebars, the number of layers, the distance between the layers, and the distance between the hoops.

The Sobols' indices calculated from the coefficients of the PCE showed that when it comes to computing the magnetic field inside a reinforced concrete building struck by lightning, the number of hoops can be neglected but the mesh size and the number of layers must be considered.

Also, based on the results of 256 simulations, we built twelve surrogate models to calculate the peak-values of the currents in three different cables routing in a reinforced concrete duct: a protective conductor represented by a bare conductor, an LV cable represented by an insulated cable, and a control cable represented by a coaxial cable. Half of the models are sparse PCEs, and the other half are Kriging models. The idea is to compare the predictions of the PCE and Kriging models to ensure the accuracy of the results. To have a larger set of predictions to compare, it would be ideal to increase the size of the dataset and build auxiliary surrogate models.

Further research could also consider additional buildings and cable ducts could.

However, since the lightning current distributes among the earth-termination system,

Finally, my biggest 2 - 7

 27 thanks to my family and friends for the unconditional support. Mom and dad, you made me into who I am. I hope I have made you proud.4 Electrical surges in underground cables due to a direct lightning strike 4.1 Parametric study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Reference case . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 First part: Routing, bonding, and grounding . . . . . . . . . . 4.1.3 Second part: Geometry of the structures and characteristics of the components . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Parameters to be considered in the numerical model of an industrial facility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Simplified model of reinforced concrete ducts . . . . . . . . . . . . . . 4.3.1 Equivalent single-layered grid . . . . . . . . . . . . . . . . . . 4.3.2 Experimental characterization . . . . . . . . . . . . . . . . . . 4.4 Estimation of the peak-values of the current in underground cables . Decision tree corresponding to the partition in figure 2-6. Figure adapted from [111]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-8 Feed-forward neural network with three inputs and one hidden layer.

Figure adapted from [ 115 3 - 1 3 - 2 3 - 4 3 - 6 3 - 7 3 - 9 𝑥 3 = 2 , 4 - 1 4 - 2 4 - 3 4 - 5 4 - 6 4 - 7 4 - 9

 1153132343637393241424345464749 Figure adapted from [115]. . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 3

 3 is dedicated to estimating the peak-values of the transient magnetic field. First, a 3-D model is conceived to study the phenomenon, then a parametric study is conducted, and finally, a formula to calculate the peak-values and a polynomial chaos expansion (PCE) to account for modifications made on the reinforcement of the building are presented.

Figure 1 - 1 :

 11 Figure 1-1: Various processes comprising a downward negative ground flash, from the separation of electrical charges in the cloud to the return-stroke process. Figure adapted from [5].

Figure 1 - 2 :

 12 Figure 1-2: External LPS of a house with air terminals on the roof. Figure adapted from [8].

Figure 1 - 3 :

 13 Figure 1-3: Structure protected according to the mesh method. Figure taken from [16].

Figure 1 - 4 :

 14 Figure 1-4: Meshed earth-termination system of an industrial facility. Figure adapted from [17].

Figure 1 - 5 :

 15 Figure 1-5: Computation model of a reinforced concrete building struck by lightning.

Figure 2 - 1 :

 21 Figure 2-1: Waveform considered for the lightning return strokes. Figure adapted from [74].

Figure 2 - 2 :

 22 Figure 2-2: Double-layered reinforcing grid embedded in concrete.

Figure 2 - 3 :

 23 Figure 2-3: Geometric view of a 2 3 factorial design.

Figure 2 - 4 :

 24 Figure 2-4: Geometric view of a CCD with three factors.

Figure 2 - 5 :

 25 Figure 2-5: Geometric view of a Box-Behnken design with three factors.

  2 and the region 𝑋 1 > 𝑡 1 is split at 𝑋 1 = 𝑡 3 , and finally, the region 𝑋 1 > 𝑡 3 is split at 𝑋 2 = 𝑡 4 . This is a predefined set of splitting rules that leads to the partition in figure2-6. The corresponding decision tree is shown in figure2-7.

Figure 2 - 6 :

 26 Figure 2-6: Partition of a twodimensional feature space.Figure adapted from [111].

  Figure 2-6: Partition of a twodimensional feature space.Figure adapted from [111].

Figure 2 - 7 :

 27 Figure 2-7: Decision tree corresponding to the partition in figure 2-6. Figure adapted from [111].

Figure 2 - 8 :

 28 Figure 2-8: Feed-forward neural network with three inputs and one hidden layer. Figure adapted from [115].

Figure 3 - 1 :

 31 Figure 3-1: Striking points considered in the model of a 5 m × 5 m × 5 m reinforced concrete building over a perfectly conducting ground plane.

Figure 3 - 2 :

 32 Figure 3-2: Distribution of the lightning current in the reinforcement of a 5 m × 5 m × 5 m building. (a) Strike in the corner. (b) Strike at 0.9 m from the downconductors. (c) Strike at 1.8 m from the down-conductors. (d) Strike in the center of the roof.

Figure 3 -

 3  shows that even when the building is struck in the center of the roof, the currents in the down-conductors in the edges are different. It can be explained by the asymmetry of the reinforcing grid, given that the length of the building (5 m) is not a multiple of the mesh size (0.3 m).

Figure 3 - 3 :

 33 Figure 3-3: Currents computed at the top of the down-conductors in the edges. (a) Strike in the corner. (b) Strike at 0.9 m from the down-conductors. (c) Strike at 1.8 m from the down-conductors. (d) Strike in the center of the roof.

Figure 3 - 4 :

 34 Figure 3-4: Distribution of the lightning current in the reinforcement of a 5 m × 5 m × 5 m building, considering a double-layered reinforcing grid and a strike in the center of the roof.

Figure 3 - 5 :

 35 Figure 3-5: Methodology proposed to find the characteristics of an equivalent singlelayered grid.

Figure 3 - 6 :

 36 Figure 3-6: Reduced representation of reinforcing grids considered to find the characteristics of an equivalent single-layered grid. (a) Infinite grid illuminated by a plane wave. (b) 2 m × 2 m × 2 m building over a perfectly conducting ground plane. (c) 4 m × 4 m grid over a perfectly conducting ground plane.

  Two variables are considered to characterize reinforcing grids: the current flowing through the grid and the magnetic field it induces. The current is computed in one of the wires connected to the enclosure, and the magnetic field is calculated as the average of the magnetic field strength computed over a 2 m × 2 m surface. As shown in figure3-7, the surface is positioned in the center of the enclosure, one meter away from the first layer of the grid. To keep the dimensions of the grid constant in the reduced representation, the meshes of the last row and the last column are not necessarily squared, they can also be rectangular. This approximation simulates what occurs when the dimensions of the building are not a multiple of the mesh size.

Figure 3 - 7 :

 37 Figure 3-7: Reduced representation of reinforcing grids chosen to find the characteristics of an equivalent single-layered grid: An 8 m × 8 m grid in a perfectly conducting metallic enclosure.

Figure 3 - 8 :

 38 Figure 3-8: Outputs computed in the reduced representation of single-layered reinforcing grids. (a) Current flowing through the grid. (b) Average of the magnetic field strength over a 2 m × 2 m surface.

Figure 3 - 9 :

 39 Figure 3-9: Flowchart of the fitting algorithm used to estimate the components of the input vector.

Figure 3 -

 3 Figure 3-10: Normalized magnetic field strength computed at 1 MHz inside and around a 5 m × 5 m × 5 m building, at 1 m from the roof. (a) Double-layered grid. (b) Single-layered grid with a mesh size of 10 cm. (c) Single-layered grid with a mesh size of 44 cm. (d) Single-layered grid with a mesh size of 30 cm. (e) Single-layered grid with a mesh size of 36 cm.

Figure 3 - 11 :

 311 Figure 3-11: Relative error on the magnetic field strength computed at 1 MHz inside and around a 5 m × 5 m × 5 m building, at 1 m from the roof. (a) Single-layered grid with a mesh size of 10 cm. (b) Single-layered grid with a mesh size of 44 cm. (c) Single-layered grid with a mesh size of 30 cm. (d) Single-layered grid with a mesh size of 36 cm.

4 :Figure 3 - 12 :

 4312 Figure 3-12: Position of the 29 points at which the magnetic field was computed inside the building.

Figure 3 - 13 :

 313 Figure 3-13: Computation model of the building when its LPS system consist of four down-conductors.

Figure 3 -

 3 14b shows that the magnetic shielding effectiveness also increases close to the striking point. It starts decreasing with the distance from the striking point, then fluctuates slightly around the same values, and then decreases. The symmetry of the curve suggests that the shielding effectiveness is affected by the distance between rebars of adjacent walls.

Figure 3 - 14 :

 314 Figure 3-14: Magnetic shielding effectiveness calculated with all the parameters set to their reference value. (a) At the vertically distributed points. (b) At the horizontally distributed points, 1.25 m from the roof.

Figure 3 - 15 :Figure 3 - 16 :

 315316 Figure 3-15: Magnetic shielding effectiveness calculated when changing the mesh size of the reinforcing grid. -------------

- 14 .

 14 The shielding effectiveness with a single-layered grid is around 20 dB in the center of the building. Adding a second layer at a distance of 60 cm from the first would lead to a shielding effectiveness of 35.2 dB. If the second layer is added 30 cm away from the first, the shielding effectiveness would be around 32.5 dB. If a third layer is added at the same distance from the second layer, it would be around 45 dB. Modifying the distance between those three layers from 30 cm to 15 cm would reduce the shielding effectiveness by 5.4 dB, making it 39.6 dB.

Figure 3 - 17 :

 317 Figure 3-17: Magnetic shielding effectiveness calculated when changing the distance between the hoops interconnecting the layers of the reinforcing grid.

Figure 3 -

 3 Figure 3-18: Magnetic shielding effectiveness calculated when additional layers are added to the reinforcing grid. -------------

Figure 3 - 19 :

 319 Figure 3-19: Magnetic shielding effectiveness calculated when changing the moisture content of the concrete.

Figure 3 -

 3 Figure 3-20: Magnetic shielding effectiveness calculated when considering the permeability of steel. (a) Considering the effect on the resistance of the rebars. (b) Considering the effect on the internal inductance of the rebars.

  ground plane: a 10 m × 10 m × 10 m shield, a 50 m × 50 m × 10 m shield, and a 10 m × 10 m × 50 m. The lightning channel was attached to the corner of the roof.

Figure 3 -

 3  shows that at points 𝐴 (𝑑 𝑤 = 𝑑 𝑟 = 5 m) and 𝐵 (𝑑 𝑤 = 𝑑 𝑟 = 3 m) as defined in the standard, the results are in good agreement. In fact, as long as we calculate the field at points that are equidistant from the walls and the roof, and inside the safety volume, the IEC formula in (3.2) appears to be a good compromise between simplicity and accuracy. However, at points that are close to the walls or the striking point, there are major differences between the results. According to the IEC formula, the distribution of the peak-values of the magnetic field strength does not change as we move along the vertical axis. As shown in figure3-22, the distance to the roof affects uniformly the magnitude of the peak-values at points over the same horizontal surface. The fact that the distribution of the peak values is defined exclusively by the distance to the walls contradicts the well-known current displacement phenomenon in grid-like shields. As we saw in section 3.1.1, no matter where the lightning channel is attached, the lightning current is always diverted to the edges. Therefore, as shown in figure3-23 and in agreement with[START_REF] Tatematsu | Analysis of electromagnetic fields inside a reinforced concrete building with layered reinforcing bar due to direct and indirect lightning strikes using the FDTD method[END_REF], the distribution of the peak values changes with the relative position to the roof and the foundation.

Figure 3 - 21 :

 321 Figure 3-21: Comparison between the peak-values of the magnetic field strength computed with TEMSI-FD and the results presented in the IEC 62305-4 standard [17].

Figure 3 - 22 :

 322 Figure 3-22: Peak-values of the magnetic field strength calculated with the IEC formula inside a 10 m × 10 m × 10 m single-layer grid-like shield, considering 𝐼 0 = 200 kA. (a) At 1 m from the roof. (b) At 3 m from the roof. (c) At 5 m from the roof.

Figure 3 - 23 :

 323 Figure 3-23: Peak-values of the magnetic field strength computed with TEMSI-FD inside a 10 m × 10 m × 10 m single-layer grid-like shield, considering 𝐼 0 = 200 kA. (a) At 1 m from the roof. (b) At 3 m from the roof. (c) At 5 m from the roof.

Figure 3 - 24 :

 324 Figure 3-24: Position of the 11 computation surfaces inside B1.

Figure 3 - 25 :

 325 Figure 3-25: Magnetic field strength normalized to the source, computed with TEMSI-FD inside B1 and B6.

Figure 3 - 26 :

 326 Figure 3-26: Maximum magnetic field strength computed with TEMSI-FD inside B1, considering 𝐼 0 = 200 kA. (a) Along the 𝑥-axis with y = 0.1 m. (b) Along the line x = y.

7 )

 7 𝐼 0 : Maximum value of the lightning current 𝑘 : Configuration factor in (3.5) 𝑑 𝑟 : Distance to the roof 𝑑 𝑒 : Distance to the closest vertical edge 𝑑 𝑤 : Shortest distance to the wall 𝑑 𝑛 : Normalized distance to the edge in the striking point 𝑣 : Volume of the grid-like shield 𝑎 : Cross-sectional area of the grid-like shield ℎ : Height of the grid-like shield 𝑐 𝑖 : Coefficients given in table 6

Figure 3 - 27 :

 327 Figure 3-27: Values taken by the variables in (3.6) to calculate the maximum magnetic field strength at point (6,9,3) inside B1.

3 . 2 .

 32 There are around 5.3 million points in total. A histogram of the residuals is shown in figure 3-28. It is worth noting that the error of the formula proposed could be assumed to have a normal distribution 𝜀 ∼ 𝒩 (0.588, 3.941 2 ), whereas a generalized extreme value distribution provides a better fit to the PDF of the residuals obtained with the IEC formula. Besides, the mean of the generalized distribution (7.824 dB) is almost as big as the 95.45% confidence interval of the normal distribution (7.882 dB ∼ factor of 2.478).

Figure 3 - 28 :Figure 3 -

 3283 Figure 3-28: Residuals obtained from the difference between the peak-values of the magnetic field strength computed in the case studies and the values calculated with the formula. (a) Considering the current of the first positive stroke (𝐼 0 = 200 kA). (b) Initially converting the peak-values into decibels.

Figure 3 - 30 :

 330 Figure 3-30: Peak-values of the magnetic field strength calculated with the improved formula inside B1, considering 𝐼 0 = 200 kA. (a) At 1 m from the roof. (b) At 3 m from the roof. (c) At 5 m from the roof.

Figure 3 - 31 :

 331 Figure 3-31: Scatter plot of the peak-values of the magnetic field strength computed in the case studies versus the values predicted with ML models, considering 𝐼 0 = 200 kA. (a) Feed-forward neural network. (b) Regression tree. (c) Bagged tree.

Figure 3 - 32 :

 332 Figure 3-32: Effect on the peak-values of the magnetic field strength calculated in the first run of the experimental design with 𝑥 1 = 10 cm, 𝑥 2 = 1 mm, 𝑥 3 = 2, 𝑥 4 = 1, and 𝑥 5 = 4. (a) Regrouped results and their fitted normal distribution. (b) Results in B1. (c) Results in B2. (d) Results in B3. (e) Results in B4.

Figure 3 - 33 :

 333 Figure 3-33: Scatter plot of the quantified effects on the peak-values of the magnetic field strength versus the effects predicted with the sparse PCE.

Figure 3 - 34 :

 334 Figure 3-34: Sobol' indices calculated from the coefficients of the sparse PCE.-------------------

Figure 3 - 35 :

 335 Figure 3-35: Effect on the peak-values of the magnetic field strength predicted by the sparse PCE for modification made on the geometrical configuration of the shield. (a) 𝑥 3 = 1, 𝑥 4 = 𝑥 5 = 0. (b) 𝑥 3 = 2, 𝑥 4 = 1, 𝑥 5 = 4. (c) 𝑥 2 = 1 cm, 𝑥 4 = 𝑥 5 = 0 if 𝑥 3 = 1, 𝑥 5 = 4 if 𝑥 3 > 1.

Finally, the

  original PCE and the 200 alternative PCEs are evaluated at the 1000 input vectors. The predictions of the original PCE are considered as the reference and subtracted from the predictions of the alternative PCEs. The resulting set of residuals gives us an idea of the deviation we may obtain from the predictions of the original PCE (see figure 3-36). However, since there is already an error on the predictions of the original PCE, we take the 38 residuals calculated as the difference between the prediction of the original PCE and the means of the fitted normal distributions. Then, to each of them, we add the set of 200 residuals. This leads to the histogram of residuals in figure 3-37. The figure shows that PDF of the error can be approximated as a normal distribution with mean 𝜇 = -0.164 dB and standard deviation 𝜎 = 1.894 dB.

Figure 3 - 36 :

 336 Figure 3-36: Residuals obtained from the difference between values predicted by the original PCE and the values predicted by the 200 alternative PCEs.

Figure 3 - 37 :

 337 Figure 3-37: Estimated error of the original PCE when considering the standard deviation of the normal distributions fitted to the calculated effects.

  -1, a 100-meters reinforced concrete duct interconnects the buildings in the middle. The duct has a cross-sectional area of 3.25 m × 2.75 m and is positioned 50 cm above the foundations of the buildings. The single-layered 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.

Figure 4 - 1 :

 41 Figure 4-1: Calculation model of the reference case.

Figure 4 - 2 :

 42 Figure 4-2: Upper view of building B1 in the reference case.

Figure 4 -

 4 3 shows a detailed view of the path the conductors follow inside B1.

Figure 4 - 3 :

 43 Figure 4-3: Detailed view of the bonding network and the path the conductors follow inside building B1.

Figure 4 - 4 :

 44 Figure 4-4: Conductors running in the reinforced concrete duct.

  figure 4-5. Note that they are normalized to the source, i.e., divided by the spectrum of the Gaussian pulse.

Figure 4 - 5 :Figure 4 - 6 :

 4546 Figure 4-5: Currents computed in the reference case. (a) Insulated conductor. (b) Bare conductor

Figure 4 - 7 :

 47 Figure 4-7: Currents computed in segment S1 when the conductors are grounded at the entrance to the buildings. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 8 :

 48 Figure 4-8: Currents computed in segment S1 when the conductors are moved on the 𝑥-axis inside B1. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 9 :

 49 Figure 4-9: Currents computed in segment S1 when the conductors are connected to different bonding rings inside B1. (a) Insulated conductor. (b) Bare conductor.

Figures 4 -

 4 Figures 4-10 and 4-11 show the results of the following cases:

Figure 4 -

 4 Figure 4-10: Currents computed in segment S1 when reducing the number of connections between the bonding rings. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 11 :

 411 Figure 4-11: Currents computed in segment S1 when changing the connections between the bonding network and the reinforcement of the buildings. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 12 :Figure 4 - 13 :

 412413 Figure 4-12: Currents computed in segment S1 when changing the connections between the grounding network and the reinforcement of the buildings. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 14 :

 414 Figure 4-14: Currents computed in segment S1 for different topologies of the grounding network. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 15 :

 415 Figure 4-15: Currents computed in segment S1 when changing the soil resistivity. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 17 :

 417 Figure 4-17: Currents computed in segment S1 when changing the distance between the reinforcement and the surface. (a) Insulated conductor. (b) Bare conductor.

Figure 4 -

 4 Figure 4-18: Currents computed in segment S1 when changing the striking point. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 19 :

 419 Figure 4-19: Currents computed in segment S1 when the duct is positioned closer to the edge. (a) Insulated conductor. (b) Bare conductor.

Figure 4 -

 4 Figure 4-20: Currents computed in segment S1 in the worst-case scenarios. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 21 :

 421 Figure 4-21: Currents computed in segment S1 when changing the connections to the ground rings and when the conductors are grounded at the entrance to the buildings. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 22 :Figure 4 - 23 :

 422423 Figure 4-22: Currents calculated in segment S1 in the worst-case scenarios for the first and the subsequent strokes. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 24 :

 424 Figure 4-24: Upper view of the calculation model of the new reference case.

Figure 4 - 25 :

 425 Figure 4-25: Currents computed in segment S1 when changing the position of the duct in the new reference case. (a) Insulated conductor. (b) Bare conductor.

Figure 4 -

 4 Figure 4-26 shows an upper view of the calculation models and figure 4-27 shows the results. It can be observed that being closer to the electrode interconnecting the ground rings does not necessarily lead to a higher current at low frequency. A

Figure 4 - 26 :

 426 Figure 4-26: Upper view of the calculation model of the cases defined to study the influence of the position of the duct.

Figure 4 - 27 :

 427 Figure 4-27: Currents computed in segment S1 in the cases defined to study the influence of the position of the duct. (a) Insulated conductor. (b) Bare conductor.

  As a remainder, the reinforcement of the buildings has a cross-sectional area of 20 m × 50 m and a height of 25 m. Out of the 25 m, 20 m are above ground. To analyze the effect, we changed the dimensions of the reinforcement of each building independently. Figures 4-28 and 4-29 show the results of the following cases: Case A: The cross-sectional area is the same but the building is rotated 90°. Case B: The cross-sectional area is reduced to 10 m × 50 m. Case C: The cross-sectional area is reduced to 20 m × 25 m. Case D: The height is increased to 45 m but the foundation is still buried at a depth of 5 m. Case E: The height is reduced to 10 m but the foundation is still buried at a depth of 5 m. Case F: The height is increased to 30 m and the foundation is buried at a depth of 10 m.

Figure 4 - 28 :Figure 4 -

 4284 Figure 4-28: Currents computed in segment S1 when changing the cross-sectional area of the buildings in the new reference case. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 30 :

 430 Figure 4-30: Currents computed in segment S1 when changing the distance between the buildings in the new reference case. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 32 :

 432 Figure 4-32: Currents computed in segment S1 when changing the mesh size of the reinforcing grid of the buildings in the new reference case. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 33 :

 433 Figure 4-33: Currents computed in segment S1 when changing the radius of the rebars and the number of layers of the reinforcing grid of the buildings in the new reference case. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 36 :

 436 Figure 4-36: Currents computed in segment S1 when changing the loads at both ends of the insulated conductor in the new reference case. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 37 :Figure 4 - 38 :

 437438 Figure 4-37: Currents computed in segment S1 when changing the section of the bare conductor in the new reference case. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 39 :

 439 Figure 4-39: Arrangements considered to calculate the radius of the equivalent insulated conductors using LAMLIM.

  The soil resistivity • The distance from the duct to the closest vertical edge of the building struck by lightning • The distance from the duct to the striking point • The number of ducts • The number of buildings • The cross-sectional area of the buildings • The distance between the buildings • The mesh size of the reinforcement of the buildings • The geometrical configuration of the reinforcement of the ductTable 4.3: Factors and levels considered in the factorial design.

Figure 4 - 40 :

 440 Figure 4-40: Current computed in the reduced representation of single-layered reinforcing grids. (a) 𝜎 = 0.1 × 𝜎 𝑟𝑒𝑏𝑎𝑟𝑠 . (b) 𝜎 = 𝜎 𝑟𝑒𝑏𝑎𝑟𝑠 . (c) 𝜎 = 10 × 𝜎 𝑟𝑒𝑏𝑎𝑟𝑠 .

Figure 4 - 41 :

 441 Figure 4-41: Average of the magnetic field strength computed in the reduced representation of single-layered reinforcing grids. (a) 𝜎 = 0.1 × 𝜎 𝑟𝑒𝑏𝑎𝑟𝑠 . (b) 𝜎 = 𝜎 𝑟𝑒𝑏𝑎𝑟𝑠 . (c) 𝜎 = 10 × 𝜎 𝑟𝑒𝑏𝑎𝑟𝑠 .

Figure 4 - 42 :

 442 Figure 4-42: Curves fitted to the data of a double-layered reinforcing grid. (a) Current flowing through the grid. (b) Average of the magnetic field strength.

Figure 4 - 43 :

 443 Figure 4-43: Currents computed in segment S1 considering an equivalent singlelayered reinforcing grid in the new reference case. (a) Insulated conductor. (b) Bare conductor.

Figure 4 - 44 :

 444 Figure 4-44: Graphic design of the testing facility.

Figure 4 - 45 :

 445 Figure 4-45: Reinforced concrete ducts covered by steel plates in the testing facility.

Figure 4 -

 4 Figure 4-46: Wire-tied reinforcement.Figure4-47: Welded reinforcement.

Figure 4 - 48 :

 448 Figure 4-48: Front view of the reinforcement of the ducts.

Figure 4 - 49 :

 449 Figure 4-49: Interconnection of two sections of the reinforced concrete duct.

Figure 4 - 50 :

 450 Figure 4-50: Interconnection of the copper plates and the reinforcement of the duct.

Figure 4 - 51 :

 451 Figure 4-51: Connection of the VNA to the copper plates. -----------------

Figure 4 - 52 :

 452 Figure 4-52: Measurement setups.

  method after performing a new open/short/load calibration in which the wires were considered. To differentiate the measurements, we named the One-Port method with the new calibration One-Port*. The results when the three sections of the ducts are interconnected are shown in figure 4-53.

Figure 4 - 53 :

 453 Figure 4-53: Impedance measured at the new lightning testing facility. (a) Duct with wire-tied reinforcement. (b) Duct with welded reinforcement.

Figure 4 - 54 :

 454 Figure 4-54: Computational model of one of the reinforced concrete ducts in the lightning testing facility.

Figure 4 - 55 :

 455 Figure 4-55: Impedance computed using TEMSI-FD.

radius of 13 .

 13 [START_REF] Tatematsu | Three-Dimensional FDTD Calculation of Lightning-Induced Voltages on a Multiphase Distribution Line With the Lightning Arresters and an Overhead Shielding Wire[END_REF] mm. The bare conductor represents a protective conductor of the minimum section one can find on the market (50 mm 2 ). The insulated conductor represents a phase of an LV cable with an insulation thickness of 2.5 mm. The coaxial cable represents a control cable of 48 conductors, each with a section of 1 mm 2 and an insulation thickness of 0.6 mm. The core of the coaxial cable is loaded with 1.04 Ω at both ends as if each of the 48 conductors was loaded with 50 Ω. Its radius, which was calculated using LAMLIM and considering the arrangement shown in figure4-56, is set to 7.862 mm. Note that to be on the safe side, we chose a typical LV cable (3×240+95) and one of the biggest control cables available on the market. One end of the insulated conductor is loaded with a parallel RLC circuit that represents a transformer. The values of the components (R = 100Ω, L = 300 𝜇H, C = 50 nF) were available from the results of different studies that have been conducted at EDF. At the other end of the insulated conductor, we consider a short circuit (SC) and an open circuit (OC) (loads of 0 Ω or 10 MΩ, respectively), which would allow us to calculate a

Figure 4 - 56 :

 456 Figure 4-56: Arrangement considered in LAMLIM to calculate the radius of the core of the coaxial cable.

Table 4 . 6 : 2 -45000 m 2 Relative distance to the closest edge 0. 1 -0. 5 Distance between the buildings 10 m - 100 m 57 and 4 - 58 .

 46221510057458 Parameters considered in the model and ranges within they are varied. Mesh size of the reinforcing grid 1 m -5 m Load at the end of the insulated conductor {0 Ω, 10 MΩ} The peak-values of the currents computed in the cables are shown in figures 4-The scatter plots in figure 4-57 highlight the difference between the resultswhen the insulated conductor is in OC and SC, whereas the histograms in figure4-58 highlight the difference between the results for first and subsequent strokes. Also, table 4.7 presents the details of the maximum and minimum currents computed in each case.

Figure 4 - 57 :

 457 Figure 4-57: Scatter plot of the peak-values of the currents in the cables computed using TEMSI-FD. (a) Bare conductor -First stroke. (b) Bare conductor -Subsequent stroke. (c) Insulated conductor -First stroke. (d) Insulated conductor -Subsequent stroke. (e) Coaxial cable -First stroke. (f) Coaxial cable -Subsequent stroke.

Figure 4 - 58 :

 458 Figure 4-58: Histogram of the peak-values of the currents in the cables computed using TEMSI-FD. (a) Bare conductor (OC). (b) Bare conductor (SC). (c) Insulated conductor (OC). (d) Insulated conductor (SC). (e) Coaxial cable (OC). (f) Coaxial cable (SC).

  2.1 Parameters considered for the different strokes. . . . . . . . . . . . .

	HGA MTLL	Hybrid genetic algorithm Modified transmission line model with linear decay with height
	HIRF NPP	High intensity radiated field Nuclear power plants
	IEC NRMSE	International electrotechnical commission Normalized root mean square error
	Acronyms INB Installation nucléaire de base OC Open circuit
	KTA PCE	Kerntechnischer Ausschuss Polynomial chaos expansion
	PDF	Probability density function
	ABC LAMLIM La méthode des lignes multiconducteurs Absorbing boundary condition PML Perfectly matched layers
	ANN LEMP	Artificial neural networks Lightning electromagnetic pulse
	LIV RMSE	Lightning induced voltage Root mean square error
	CCD	Central composite design
	LOO RNN	Leave-one-out Recurrent neural network
	DoE LPL SC	Design of experiments Lightning protection level Short circuit
	DUT LPS SPD	Device under test Lightning protection system Surge protective device
	EFIE LPZ	Electric field integral equation Lightning protection zone
	TEM	Transverse electromagnetic
	EM LV	Electromagnetic Low voltage
	TEMSI-FD Transient electromagnetic simulator finite difference
	FD MAE TL	Finite difference Mean absolute error Transmission line
	FDTD MC	Finite-difference time-domain Monte Carlo
	FEKO ML VF	Machine learning Feldberechnung für körper mit beliebiger oberfläche Vector-fitting
	FMM MoM VNA	Method of moments Fast multipole method Vector network analyzer
	FT MPM	Fourier transform Matrix pencil method
	GP MSE	Gaussian process Mean squared error
	GPR MTLE	Gaussian process regression Modified transmission line model with exponential decay with height

2.2 Design matrix of a 2 3 factorial design. . . . . . . . . . . . . . . . . . . 2.3 Extended design matrix of a 2 3 factorial design divided into two half fractions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Classical orthogonal polynomials. . . . . . . . . . . . . . . . . . . . . 2.5 Base of the PCE assuming a second order truncation. . . . . . . . . . 2.6 Commonly used kernels. . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Commonly used measures of error. . . . . . . . . . . . . . . . . . . . 3.1 Factors and levels considered in the factorial design. . . . . . . . . . . 3.2 Best adjusted models based on the NRMSE. . . . . . . . . . . . . . . 3.3 Parameters and error of the ordinary Kriging models with the Matern 5/2 kernel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.4 Parameters chosen to be analyzed in the study. . . . . . . . . . . . . 3.5 Dimensions of the grid-like shields considered. . . . . . . . . . . . . . 3.6 Coefficients of the improved formula in (3.7). . . . . . . . . . . . . . .

Table 2 .

 2 2: Design matrix of a 2 3 factorial design.

	Factors

Table 2 .

 2 3: Extended design matrix of a 2 3 factorial design divided into two half fractions.

			Factors			Interactions
	Run	𝐴	𝐵	𝐶	𝐴𝐵	𝐴𝐶	𝐵𝐶	𝐴𝐵𝐶
	1

Table 2 .

 2 4: Classical orthogonal polynomials.

	Type of variable	Orthogonal	Hilbertian basis
		polynomials	Ψ 𝛼 (𝑥)
	Uniform	Legendre 𝑃 𝛼 (x)	𝑃𝛼(𝑥) √ 1
			2𝛼+1

Table 2 .

 2

6: Commonly used kernels. Kernel 𝑅(𝑋, 𝑋 ′ ; 𝜃)

Table 2

 2

.7 summarizes the most commonly used in statistics. Note that 𝑁 is the number of model evaluations and 𝜇 is the sample mean of the quantity of interest.

Table 2.7: Commonly used measures of error.

Table 3 .

 3 2: Best adjusted models based on the NRMSE.

	Model	NRMSE

Table 3 .

 3 3: Parameters and error of the ordinary Kriging models with the Matern 5/2 kernel.

		Model of the current	Model of the magnetic field
			strength
	𝛽 1	36.101	1.337
	𝜎 2	259.245	0.318
	𝜃	[8.534, 7.740, 1.263]	[7.282, 6.924, 1.157]
	𝐸 𝐿𝑂𝑂	4.831 × 10 -6	4.975 × 10 -6

Table 3 .

 3 5: Dimensions of the grid-like shields considered.

	Building	Length (m)	Width (m)	Height (m)
	B1	10.100	10.100	10.100
	B2	10.100	10.100	20.100
	B3	20.100	20.100	10.100
	B4	20.100	10.100	10.100
	B5	5.100	5.100	5.100
	B6	20.100	20.100	20.100
	B7	5.100	5.100	10.100
	B8	10.100	10.100	50.100
	B9	50.100	50.100	10.100
	B10	50.100	20.100	10.100

Table 3 .

 3 6: Coefficients of the improved formula in (3.7).

	Coefficient	Estimated value
	𝑐 1	1.146 × 10 -3
	𝑐 2	0.665
	𝑐 3	1.625
	𝑐 04	5.795 × 10 -3
	𝑐 05	0.585
	𝑐 6	0.520
	𝑐 7	1.491
	𝑐 8	1.413

Table 3 .

 3 7: Prediction error of the ML algorithms compared to the error calculated for the formulas.

	Model

account and the ranges within they are varied are presented in table 3.8. Because of the proportionality we observed in section 3.2, the distances are given as a multiple

Table 3 .

 3 8: Ranges within the characteristics of the grid-like shield are varied.

	Variable	Range

Table 4 .

 4 1: Highest normalized currents observed for each parameter in segment S1.

	Description	Current in the insulated conductor (dB)	Current in the bare conductor (dB)
		10 kHz	First resonance	10 kHz	First resonance
	Reference	-95.50	-45.99	-74.13	-41.08
	Conductors connected to
	the bonding ring on the	+ 4.45	-19.51	+ 5.07	-19.08
	first floor of B1
	Conductors connected to
	the bonding ring on the	-2.24	+ 0.50	-1.10	+ 0.55
	first floor of B2
	Interconnection of the
	bonding rings inside B2	+ 1.17	+ 0.96	+ 1.10	+ 0.34
	in corner C22
	Interconnection of the
	bonding network and	+ 0.66	+ 1.36	+ 0.60	+ 1.09
	the reinforcement of B1
	in C12
	Interconnection of the
	reinforcement of B1 and	+ 5.31	-0.16	+ 4.58	-0.17
	the ground ring in C11
	One connection between	+ 9.40	+ 0.09	+ 8.53	+ 0.08
	the ground rings
	Strike in corner C11	+ 6.69	-12.87	+ 6.08	-12.66
	Duct at 1 m from the
	edge and strike in corner	+ 8.36	-11.74	+ 7.89	-5.78
	C12

Table 4 .

 4 4: Parameters and error of the ordinary Kriging models with the Matern 5/2 kernel.

		Model of the current	Model of the magnetic field
			strength
	𝛽 1	-4.009	-0.008
	𝜎 2	249.559	0.269
	𝜃	[6.711, 9.877, 0.346, 0.835]	[4.632, 9.827, 0.319, 0.779]
	𝐸 𝐿𝑂𝑂	6.218 × 10 -6	5.847 × 10 -6

Table 4 .

 4 5: Resistance and inductance calculated from the impedance measured on different sections of the ducts.

	Section	Welded reinforcement	Wire-tied reinforcement
		R (mΩ)	L (𝜇H)	R (mΩ)	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

Table 4 .

 4 7: Maximum and minimum currents computed in the cables.

	Maximum	Minimum	Average
	Run	𝐼(𝐴)	Run	𝐼(𝐴)

Table 4 .

 4 8: Details and error of the PCEs.

	Output	Maximal	Size of basis	𝐸 𝐿𝑂𝑂
		degree
	𝑌 1	2	16	1.656 ×10 -2
	𝑌 2	2	18	3.055 ×10 -2
	𝑌 3	2	19	5.431 ×10 -2
	𝑌 4	2	17	5.889 ×10 -2
	𝑌 5	2	15	1.366 ×10 -2
	𝑌 6	2	12	1.410 × 10 -1

Table 4 .

 4 9: Parameters and error of the linear Kriging models with the Matern 5/2 kernel. 𝛽 1 [65.406, 3.016, -5.944, 1.278, 3.642, -1.660, -1.851, 0.185] 𝛽 1 [44.896, 2.712, -4.846, -0.568, 1.267, -1.574, -3.389, 3.162] 𝛽 1 [-12.711, 3.343, 12.514, -0.450, 0.704, -1.941, -3.311, 1.692] 𝛽 1 [25.611, 4.431, 9.259, -1.987, 0, -0.934, -6.082, 2.866]

	𝑌 1	𝜎 2 𝜃	2.769 [3.983, 1.920, 2.193, 1.585, 1.144, 1.469, 2.673]
		𝐸 𝐿𝑂𝑂	6.237 × 10 -3
		𝛽 1	[45.394, 2.864, 2.852, 1.116, 3.432, -1.693, -1.839, 0.376]
	𝑌 2	𝜎 2 𝜃	3.530 [3.741, 2.224, 2.766, 1.412, 0.375, 1.938, 4.400]
		𝐸 𝐿𝑂𝑂	7.493 × 10 -3
	𝑌 3	𝜎 2 𝜃	5.805 [1.254, 1.815, 6.340, 2.165, 1.568, 1.483, 2.791]
		𝐸 𝐿𝑂𝑂	6.903 × 10 -3
	𝑌 4	𝜎 2 𝜃	19.088 [1.742, 0.817, 2.395, 0.900, 1.160, 1.194, 0.392]
		𝐸 𝐿𝑂𝑂	4.88 × 10 -2
	𝑌 5	𝜎 2 𝜃	5.486 [1.256, 1.100, 9.993, 5.988, 0.911, 0.959, 3.217]
		𝐸 𝐿𝑂𝑂	4.792 × 10 -4
		𝛽 1	[51.347, 3.110, 2.483, 0.914, 3.250, -1.997, -1.818, 0.253]
	𝑌 6	𝜎 2 𝜃	21.724 [1.106, 2.670, 3.448, 2.302, 1.806, 0.554, 1.664]
		𝐸 𝐿𝑂𝑂	9.778 × 10 -2

Table 4 .

 4 10: Currents estimated using the models in tables 4.8 and 4.9

	Output	Predicted value (A)
		PCE		Kriging
	𝑌 1	594.752		615.124
	𝑌 2	156.722		161.420
	𝑌 3	64.253		61.537
	𝑌 4	0.704		0.635
	𝑌 5	30.977		33.116
	𝑌 6	316.225		289.546
	If we consider the OC resistance equal to 10 MΩ, we can calculate the Norton
	resistance as	𝑅 =	10 MΩ × 𝐼 𝑂𝐶 𝐼 𝑂𝐶 -𝐼 𝑂𝑆	,	(4.1)

Acknowledgments

Fractional factorial designs [START_REF] Montgomery | Two-Level Fractional Factorial Designs[END_REF] In some cases, one cannot afford to run all the combinations of a full factorial design. A full factorial design with 𝑘 factors at 𝑀 levels requires 𝑀 𝑘 runs. Thus, as the number of factors increases, the number of runs required rapidly outgrows the computational resources available. Fractional factorial design can be used to reduce the number of runs by focusing for example on the main effects and the low-order interactions. Note that by doing so, certain high-order interactions are neglected.

Since only a fraction of the design is considered, it is essential to ensure that the effects of interest are not aliased with each other when choosing the generators. There are two fundamentally different approaches to determining the aliasing pattern based on the effects one wants to estimate: the "Clear" approach and the "Distinct" approach [START_REF] Grömping | Clear" and "Distinct": two approaches for regular fractional factorial designs with estimability requirements[END_REF]. In the "Clear" approach, only the main effects and two-factor interactions are considered significant. In the "Distinct" approach, the generators are chosen to guarantee the estimability of the effects that the experimenter thinks are significant. Thus, it is assumed that the experimenter knows the nature of the model. A wellknown algorithm to choose the generators is the Franklin-Bailey algorithm [START_REF] Franklin | Selection of defining contrasts and confounded effects in two-level experiments[END_REF]. The Franklin-Bailey algorithm finds the generators for the smallest fractional factorial

Configuration of the grounding network

As with the configuration of the bonding network, changing the number of connections to the grounding network can also force the currents to take a different path. However, since the foundation of the building also provides a path for the currents to flow to earth, a shift in the frequency of the first resonance is hardly noticeable. Figure 4 Naturally, the currents increase if the building struck by lightning is not connected to the grounding network. But once again, there is always at least one connection.

The analysis is then similar to the one made for the interconnection of the bonding rings: the coupling can be stronger if the number of connections is reduced and placed close to the path of the conductors inside the building.

In general, the earth electrodes forming the ground rings and the meshed grounding network are buried at a depth of approximately 1 m. As shown in figure 4-13, small variations in the depth have a minor effect on the currents. The currents increase slightly at low frequency when the electrodes approach the surface, perhaps because the coupling between the grounding network and the reinforcement of the other structures is stronger.

Characteristics of the concrete

Since the foundations of the buildings are also part of the earth-termination system, but they are not in direct contact with the soil, the moisture content of the concrete and its thickness could also influence the results.

We already know that the concrete does not affect the magnetic field generated inside the building by a direct lightning strike. Figure 4-16 shows that the effect of varying its moisture content on the induced currents can also be neglected. On the contrary, increasing the distance from the reinforcement to the surface leads to significant changes. Figure 4-17 shows the results of the following cases:

Case A: The concrete is not modeled. The smaller the distance between the surface and the reinforcement, the higher the currents. For instance, a significant rise is observed at low frequency when the concrete is not modeled. Yet, the effect of increasing the thickness of the concrete around the buried rebars is also negligible. If we look at the concrete as an insulating

Interconnection of the duct and the buildings

It is clear that the currents induced in the conductors depend on the current flowing through the duct, which depends on the characteristics of the duct itself and its connection to the buildings. We had assumed that the reinforcement of the duct was connected four times to the reinforcement of each building. However, the connection could be made by different means, and in some cases a good continuity is not guaranteed. Figure 4-31 shows the results of the following cases:

Case A: The reinforcement of the duct was connected twice to the reinforcement of each building. The connections are at the top corners.

Case B: The reinforcement of the duct is connected once at each end to the bare conductor, which is also connected to the reinforcement of the buildings.

Case C: There are no connections between the reinforcement of the duct and the reinforcement of the buildings.

Case D: The duct has no reinforcement. As expected, the currents rise significantly if we remove the reinforcement of the duct, because we are reducing the shielding effectiveness it provides for the cables running inside. Also, it can be observed that the increase is comparable to the case in which the reinforcement of the duct and the reinforcement of the buildings are not interconnected. It is important to ensure there is at least one good connection, even if it is made using the protective conductor.

Geometrical configuration of the reinforcement of the buildings

We saw in the previous chapter that it is possible to replace a multi-layered reinforcing grid with an equivalent single-layered grid to calculate the electromagnetic field inside a building struck by lightning. Moreover, we considered a single layer in the reference case. However, one cannot define a unique equivalent grid because the configuration of the reinforcement can vary depending on the building. To analyze the effect of possible changes, we considered the following cases:

Case A: The mesh size of the reinforcing grid is reduced to 25 cm.

Case B: The mesh size of the reinforcing grid is increased to 1 m.

Case C: The mesh size of the reinforcing grid is increased to 5 m.

Case D: The radius of the rebars is reduced to 5 mm.

Case E: The radius of the rebars is increased to 2 cm.

Case F: There is a double-layered reinforcing grid.

As expected from the distribution of the lightning current studied in section 3.1.1, figures 4-32 and 4-33 show that changing the radius of the rebars or adding an extra layer does not affect the results. In contrast, increasing the mesh size, especially of B1, can lead to a significant increase in the currents. The question that follows is whether we can conclude on the effect of the configuration of the reinforcement, considering that the conductors are grounded at the entrance to the buildings. Perhaps the number of layers does have an important effect when the bonding system is modeled.

To verify, we repeated the simulations of the worst-case scenarios defined in the first part of the study but considering a double-layered reinforced grid.

The results are shown in figure 5.34. As expected, there are some differences because adding an extra layer reduces the magnitude of the electromagnetic field typically have a double-layered reinforcing grid with squared meshes of 25 cm and rebars with a radius between 4 mm and 8 mm. Still, we considered the following cases:

Case A: The mesh size of the reinforcing grid is increased to 50 cm.

Case B: There is a double-layered reinforcing grid.

Case C: The radius of the rebars is reduced to 4 mm.

Case D: The radius of the rebars is increased to 16 mm.

Figure 4-35 shows that variations in the reinforcement of the duct have a much more important effect on the results. The effect on the surges is actually similar to the effect of the reinforcement of the building on the magnetic field. The surges increase with the mesh size and decrease with the radius and the number of layers.

We could therefore think about finding an equivalent single-layered grid of a typical duct to use in the simulations.

Loads

Once again, when a building 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 Case C: The insulated conductor is loaded with 3 nF and 10 nF.

Case D: The insulated conductor is loaded with 8 𝜇H and 2.5 nF.

Case E: The insulated conductor is loaded with 800 𝜇H and 100 nF.

As expected, figure 4-36 shows that when the insulated conductor is not loaded, the currents in both the insulated and the bare conductors are in the same range.

In contrast, the current in the insulated conductor decreases significantly when a

• The loads connected at both ends of the cables

• The characteristics of the cables

That is if the conductors are grounded at the entrance to the buildings, and only the ground rings are represented in the model. Otherwise, the configuration of the bonding and the grounding networks must be considered. Being cautious to always remain on the safe side, the other parameters can be set to fixed values. Note however that we did not consider all the components of industrial facilities. Future studies could explore the effect of, among others, the support rails, lines in contact with soil, water pipes, and the grounding arrangement.

Simplified model of reinforced concrete ducts

Equivalent single-layered grid

The geometrical configuration of the reinforcement of the buildings can vary drastically within the facility. In contrast, as already mentioned, the reinforcement of the ducts is basically the same; moreover, the least we can do is try to find a simplified 3-D model of the reinforcement, as we did for the reinforcement of the buildings in section 3.1.

We start once again by considering the reduced representation of a single-layered reinforcing grid as an experiment but with three inputs this time: the radius of the rebars, the mesh size, and the conductivity of the rebars. As the reduced representation, we use the same 8 m × 8 m grid in a metallic enclosure, defined in section 3.1.2.

The outputs of the experiment are still the current flowing through the grid, and the average of the magnetic field it induces over a 2 m × 2 m surface.

Based on a factorial design with mixed levels, 75 simulations are carried out to build a Kriging surrogate model for each output. The levels taken by the factors are shown in table 4.3. Since, theoretically, the conductivity of the rebars only affects the resistance of the grid, we did not think it would be necessary to define a 5 3 factorial design, and therefore, we considered a logarithmic scale to select the levels.

the cable ducts, and soil-contacting conductors leading away from the building during a direct strike, by considering one duct we consider a worst-case scenario.

The measurements of the impedance of the ducts at the new lightning testing facility of EDF Power Networks Lab proved that our representation of reinforced cable ducts in numerical models is appropriate. Yet, those were the first measurements performed at the installation. There is still a great number of tests that could be conducted to validate the representation of other components of the electromagnetic environment in our numerical models.

The error estimated for the surrogate models presented in this work was calculated considering the results obtained from full-wave simulations as the true values. The error of the 3-D models used to carry out the simulations is yet to be quantified.

List of Publications