Laurent Brochard
email: laurent.brochard@enpc.fr

Sabri Souguir

Karam Sab

Scaling of brittle failure: strength versus toughness

Keywords: strength, toughness, scaling law, graphene PACS 62.20.-x

We study the scaling of strength and toughness in function of temperature, loading rate and system size, to investigate the difference between tensile failure and fracture failure. Molecular simulation is used to estimate the failure of intact and cracked bodies while varying temperature, strain rate and system size over many orders of magnitude, making it possible to identify scaling laws. Two materials are considered: an idealized toy model, for which a scaling law can be derived analytically, and a realistic molecular model of graphene. The results show that strength and toughness follow very similar scalings with temperature and loading rate, but differ markedly regarding the scaling with system size. Strength scales with the number of atoms whereas toughness scales with the number of cracks. It means that intermediate situations of moderate stress concentrations (e.g., notch) can exhibit not obvious size scaling, in-between those of strength and toughness. Following a theoretical analysis of failure as a thermally activated process, we could rational-

Introduction

Understanding mechanical failure is essential for the reliability of structures. Yet, the modeling of failure initiation remains quite debated. Many initiation criteria have been proposed in the literature but no consensus exists [START_REF] Brochard | Scaling laws and renormalization groups for strength and toughness of disordered materials[END_REF]. The various theories aim at reconciling two well understood limit cases, namely tensile strength and fracture toughness. Tensile strength, which describes the failure of flawless materials, is captured by a stress criterion (failure surface). In contrast, fracture toughness, which describes the failure of cracked bodies, is captured by an energy released criterion (linear elastic fracture mechanics). But a general theory of failure must be able to handle intermediary situations, that is to combine stress and energy in a more general criterion. So far, initiation theories have been postulated; e.g., cohezive zone models [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF], non-local approaches [START_REF] Novozhilov | On a necessary and sufficient criterion for brittle strength[END_REF] or finite fracture mechanics [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF]. But a rigorous approach, necessary to reach a consensus, requires to understand the elementary physics at the heart of mechanical failure.

The most elementary scale of mechanical failure is the scale of atoms and molecules where atomic bonds break and crack surface is created. While the physical evolution of macroscopic systems is described by continuum thermodynamics, at the atomic scale, fluctuations become significant and need to be accounted for, which is the purpose of statistical physics. In particular a general theory of failure should be able to relate the failure stress to the temperature-induced fluctuations through general scaling laws. Such a theory, called kinetic theory of strength, was proposed initially by Zhurkov (see [START_REF] Regel | The kinetic nature of the strength of solids[END_REF] or [START_REF] Zhurkov | Kinetic concept of the strength of solids[END_REF]). This theory captures how failure stress σ f ailure depends on both temperature T and loading time τ through a general scaling law that has been verified for a wide variety of materials : τ ∝ exp ∆E(σ f ailure ) kT

, with ∆E (σ f ailure ) > 0 the energy barrier to reach the transition state to failure. One can recognize in this law the Boltzmann factor characteristic of thermally activated processes. Confrontation with experiments shows that, for most materials, ∆E ≈ E 0 -γσ, hence :

σ f ailure σ f ailure (0K) = 1 - T T cr ln τ τ cr (1)
with τ cr and T cr two constants with dimensions of time and temperature respectively that depend on the material. The derivation of Equation 1 assumes that failure is the consequence of the occurrence of a transition state (e.g., breaking of a bond) associated with an energy barrier ∆E and with no reverse process (e.g., reformation of the bond), which is valid in the limit of small temperatures (∆E reverse -∆E kT ). Interestingly, thermally activated scaling laws are not confined to brittle failure but is encountered in many other failure mechanisms, although the expression of energy barrier ∆E may vary from one mechanism to another leading to different exponents on T ln (τ ) in Eq. ( 1) : plasticity of metallic glasses (Argon, 1979;[START_REF] Johnson | A Universal Criterion for Plastic Yielding of Metallic Glasses with a (T/Tg)ˆ(2/3) Temperature Dependence[END_REF][START_REF] Schuh | Mechanical behavior of amorphous alloys[END_REF], dislocation nucleation in ductile crystals [START_REF] Zhu | Temperature and Strain-Rate Dependence of Surface Dislocation Nucleation[END_REF], yield by bond rotation in carbon nanotube [START_REF] Wei | Tensile strength of carbon nanotubes under realistic temperature and strain rate[END_REF][START_REF] Dumitrica | Symmetry-, time-, and temperature-dependent strength of carbon nanotubes[END_REF]. A particular case of interest here is fracture failure, which has been investigated theoretically in the framework of thermal activation by Petrov andOrlov (1975, 1976). A key consideration pointed out by Petrov and Orlov is the size effect. Indeed, many 'competing' transitions states can lead to failure. For instance, in a body with many identical micro-cracks, each micro-crack can be the source of a macro-crack, but only one does actually become a macro-crack. This points to a difference in the scaling between fracture toughness and tensile strength, since the presence of flaws concentrates the transition states in the vicinity of the flaws.

In the literature, scaling and size effects of failure have attracted attention mostly in the context of distribution of flaws (Carpinteri, 1994;[START_REF] Carpinteri | Are scaling laws on strength of solids related to mechanics or to geometry[END_REF] and of large scale yielding [START_REF] Bazant | Scaling of Structural Failure[END_REF][START_REF] Brochard | Capturing material toughness by molecular simulation: accounting for large yielding effects and limits[END_REF]. In contrast, the combined temperature-loading rate-size scaling of the kinetic theory has attracted only little attention. Part of the difficulty stems to the logarithmic scaling which requires very large variations in loading rate and system size for proper investigation. Only the temperature scaling is easily achieved experimentally. In this respect, molecular simulations techniques offers an interesting alternative to experiments since one can easily vary the loading rate and the system size over orders of magnitudes. [START_REF] Slutsker | Atomic-Level Fluctuation Mechanism of the Fracture of Solids (Computer Simulation Studies)[END_REF] propose a detailed study of Zhurkov's theory by molecular simulation but, to the best of our knowledge, the confrontation between strength and toughness scaling has not been studied so far. The purpose of this paper is to investigate with molecular simulation techniques the scalings of strength (tensile failure) and toughness (fracture failure) and confront the two scalings. We show, for a toy model (triangular lattice) and for graphene, that their failure is a thermally activated process exhibiting the expected Zhurkov-like scaling. Confronting strength and toughness, we highlight the essential difference in size scaling between fracture toughness and tensile strength. Thus, even though tensile failure and fracture failure arise both from the occurrence of a transition state, the distribution of these states differs radically. This study can serve as a starting point for the elaboration of a physically-based theory of initiation.

In section 2, we present the systems considered and the methodology used. In section 3, we present the results and we propose a theoretical analysis of the observed scaling in section 4.

Methods

We use molecular dynamics (MD) simulations to investigate failure at the atomistic scale. Molecular dynamics techniques give access to a space and time resolution hardly accessible experimentally. In particular, in this work, we focus on the scaling of the stress at failure with respect to temperature, system size and loading rate. Investigating such scaling requires to vary all those parameters over several orders of magnitude. Even though conventional MD is limited to millions of atoms simulated over nanoseconds, one can reasonably investigate 4 to 5 decades of scaling, which is enough for the present purpose.

We study two different materials: a toy model 2D solid and a realistic model of graphene. The respective atomic structures of these two models are displayed in Figure 1. The toy model is a crystalline material made of a regular triangular lattice in which the inter-atomic potential is limited to harmonic interactions between nearest neighbors:

U = i<j U harmo (r ij ) with U harmo (r) = K 2 (r -r 0 ) 2 -U cr if r < r cr 0 otherwise (2)
where r ij is the distance between atoms i and j, K, r 0 and r cr > r 0 are the force constant, equilibrium distance and critical distance of the potential. The shift term U cr = K 2 (r cr -r 0 ) 2 ensures continuity of the energy at r = r cr . The potential is made irreversible, that is, once a bond is broken (r > r cr ) then it is removed from the system and can never reform even if the distance between the two atoms returns to less than r cr . Doing so, we avoid any reversibility of the failure process (i.e., no 'healing' of bond).

The interest for such a toy model is that the results can be easily interpreted and confronted to theory. Toy models (including the triangular lattice) have been widely studied in the literature since the 1970's to understand the physics of failure. [START_REF] Thomson | Physics of Fracture[END_REF] treated this type of toy model via Green's functions and, notably, identified energy barriers to failure similar to the expressions obtained in the present paper. [START_REF] Slepyan | Models and Phenomena in Fracture Mechanics[END_REF] found analytical solutions for crack dynamics and considered in particular the triangular lattice. Many pursued these investigations in the following decades. Let us mention [START_REF] Marder | Origin of Crack Tip Instabilities[END_REF]; [START_REF] Thomson | Lattice imperfections studied by use of lattice Greens functions[END_REF]; [START_REF] Pechenik | Steady-state mode I cracks in a viscoelastic triangular lattice[END_REF] among others, who worked on the triangular lattice. To the best of our knowledge, the scaling study we report in this paper has never been done before.

For the triangular lattice we consider, the elastic mechanical behavior and tensile strength can be anticipated as follows, under small strain hypothesis at 0K: elastic behavior:

       σ xx = √ 3K 3εxx+εyy 4 σ yy = √ 3K εxx+3εyy 4 σ xy = √ 3 2 Kε xy (3)
tensile strength:

       σyy √ 3 +σxy K ≤ rcr r0 -1 σyy √ 3 -σxy K ≤ rcr r0 -1 3σxx-σyy 2 √ 3K ≤ rcr r0 -1 (4)
According to equations (3), one expects a simple linear elastic isotropic behavior with bulk modulus √ 3K/2 and shear modulus √ 3K/4. The tensile strength (Eq. 4), however, exhibits a non trivial surface. In particular, the strength is clearly anisotropic. The relations (2), ( 3) and ( 4) call for considering the following dimensionless quantities:

r * = r/r 0 , U * = U/ Kr 2 0 , σ * -= σ -/K, T * = k B T / Kr 2 0 (k B is Boltzmann constant).
The results presented throughout this paper refer to those quantities.

The theoretical behavior given by Equations ( 3) and (4) assumes static equilibrium of the lattice which is the limit of low temperature (0K) and quasi-static loading. At finite temperatures and finite loading rate, thermal agitation and lattice vibrations (phonons) give rise to deviations from this theoretical prediction. In the following sections, we investigate those deviations. We perform MD simulations of periodic systems with r * cr = 1.1 and compute the tensile strength for various loading orientations, system sizes, temperatures and loading rates. More precisely, we consider 84 different loading orientations ranging from purely volumetric to purely deviatoric. We consider sizes ranging from N = 30 to N = 491520 atoms, temperatures from T * = 10 -4 to T * = 1, and loading rates from ε * = 6.87 • 10 -11 to ε * = 2.06•10 -5 . Here, we define the dimensionless loading rate ε * = dε dt * with respect to the dimensionless time t * = t/ K/m where m is the mass of the atom. The Cauchy stress tensor is computed with the usual virial estimate [START_REF] Allen | Computer simulation of liquids[END_REF]. All these simulations are performed under isothermal conditions imposed by a Langevin thermostat and a velocity-Verlet integration scheme [START_REF] Frenkel | Undestanding Molecular Simulation: From Algorithms to Applications[END_REF]. In this work, we preferred to use a Langevin thermostat rather than the more conventional Nose-Hoover thermostat because the Nose-Hoover thermostat is known to be non-ergodic in the particular case of the harmonic oscillator [START_REF] Legoll | Non-Ergodicity of the NoséHoover Thermostatted Harmonic Oscillator[END_REF], which is an issue with respect to the toy model we study. Considering isothermal conditions is valid in this work because the system considered exhibit no irreversibilities untill the onset of failure. In the general case, the temperature field is not uniform, for instance, plasticity at a crack tip induces a strong heating [START_REF] Rice | Local Heating by Plastic Deformation at a Crack Tip[END_REF], which can dramatically affect the failure behavior [START_REF] Zehnder | On the temperature distribution at the vicinity of dynamically propagating cracks in 4340 steel[END_REF][START_REF] Flores | Enhanced toughness due to stable crack tip damage zones in bulk metallic glass[END_REF][START_REF] Ponson | Cleaved surface of ¡math display="inline"¿ ¡mrow¿ ¡mi¿i¡/mi¿ ¡mtext¿¡/mtext¿ ¡mi mathvariant="normal"¿Al¡/mi¿ ¡mi math-variant="normal"¿Pd¡/mi¿ ¡mi mathvari-ant="normal"¿Mn¡/mi¿ ¡/mrow¿ ¡/math¿ quasicrystals: Influence of the local temperature elevation at[END_REF][START_REF] Wang | Instability of crack propagation in brittle bulk metallic glass[END_REF].

Besides tensile failure, we also study the fracture failure of the toy model. One can estimate an approximate theoretical value of the toughness of the toy model by considering that fracture failure occurs when the first bond at a crack tip breaks. According to fracture mechanics, for a crack in the x direction and a loading in the y direction, the stress field is singular at the tip with an asymptotic stress ahead of the tip σ yy (r) ≈ K Ic √ 2πr , where r is the distance to the crack tip and K Ic is the toughness. The first bond at tip breaks when the force it supports exceeds its carrying capacity. As a first approximation, one can estimate the stress it support by the integration of the stress singularity over the the length supported by the bond: r0/2 0 σ yy (r) dr = K Ic r 0 /π (for a crack in the x direction, the first bond is oblique and supports the load in y over a length r 0 /2). The carrying capacity is given by the first two criteria in Equation ( 4) for the oblique bonds: r0 2 √ 3K rcr r0 -1 . We thus obtain the following estimate of the toughness of the toy model:

K Ic ∼ K 3π 4 r 0 r cr r 0 -1 (5)
Note that, unlike Equations (4) for strength, this is only an estimate that provides an order of magnitude. To study the fracture failure by molecular simulations, we simulate flawed systems with initial cracks initiated inside the material as displayed in Figure 2. A crack is initiated by removing bonds along a line. As for tensile failure, we consider periodic boundary conditions for the simulation of fracture failure. This choice is made to avoid spurious interface effects at non periodic boundaries. As a consequence, the cracked system is surrounded by periodic replica of the crack. The presence of a periodic crack has to be accounted for carefully in order to estimate the fracture toughness from the stress at failure. According to linear elastic fracture mechanics, when the system is submitted to a remote loading in the direction orthogonal to periodic cracks, the stress intensity K I in mode I at the cracks' tips is :

K I = Σ √ πaC 2a L , H L (6)
where Σ is the average stress in the direction orthogonal to the crack, a is the half crack length and L and H are the x and y dimension of the periodic cell. One can identify the usual stress intensity for a finite crack in an infinite body (Σ √ πa), and the coefficient C 2a L , H L is a correction due to the presence of periodic cracks.

In particular, as 2a/L → 0, C → 1 and one recovers the non-periodic case; and, as 2a/L → 1, C → +∞ since periodic cracks merge in this limit. There exists no analytic expression of C 2a L , H L but several numerical estimate have been proposed. In this work, we use the numerical estimate of [START_REF] Karihaloo | Doubly periodic arrays of bridged cracks and short fibre-reinforced cementitious composites[END_REF]. According to fracture mechanics, failure occurs when the stress intensity reaches a critical value K Ic , called toughness: Accordingly, the toughness and the average stress Σ cr at failure are related according to:

Σ cr = K Ic √ πaC 2a L , H L (7)
Equation ( 7) can be used to estimate the toughness of the material from the result of a molecular simulation of the periodic system with initial crack. Note however, that Equation ( 7) is based on linear elastic fracture mechanics which assumes small scale yielding [START_REF] Anderson | Fracture Mechanics: Fundamentals and Applications[END_REF]). Therefore, it is valid only if the crack length and the distance between periodic replica is sufficiently large compared to the characteristic length r pl = (K Ic /σ strength )

2 of the process zone at the crack tip.

In practice, r pl is often much larger than the size of an atom, and using Equation ( 7) at the atomic scale is questionable. We investigated this issue in previous work [START_REF] Brochard | Capturing material toughness by molecular simulation: accounting for large yielding effects and limits[END_REF][START_REF] Brochard | Scaling laws and renormalization groups for strength and toughness of disordered materials[END_REF]. Regarding the present study, the toy model has a very low process zone: using the toughness estimation (Eq. 5), we obtain r pl = πr 0 /4. With the simplistic pair potential considered, fracture failure arises from a single bond breaking, that is r pl is about the size of an atom. So we expect Equation (7) to be reproduced precisely at the atomic scale.

In contrast, the characteristic length r pl is about 1.5 nm for graphene and capturing the toughness with Equation ( 7) is possible for system at least 10 nm large which starts to be computationally expensive (see [START_REF] Brochard | Scaling laws and renormalization groups for strength and toughness of disordered materials[END_REF]). So, we limit ourselves to the study of tensile failure for graphene. We simulate fracture failure of the toy model for various system sizes, crack lengths, temperatures and loading rates. In all the simulations, cracks are initiated in the direction x and the systems are loaded in the y direction. Therefore, one obtains only the mode I toughness for crack propagation in the x direction. We consider crack lengths ranging from 2a/L = 0 (i.e., tensile failure) to 2a/L = 1 (i.e., fully separated). The range of system sizes and temperatures considered is the same as for tensile failure. Regarding loading rates, we consider rates one order of magnitude smaller than for tensile failure. Indeed, as will be discussed in section 4, excessive rates tend to introduce bias in our results and this bias appears more significant for the toughness estimate than for the strength estimate.

In addition to temperature, size and loading rate, we also investigate the scaling with respect to the number of initial crack in the periodic cell. To do so multiple cracks are initiated periodically in the system as displayed in Figure 3. The system being periodic, one recovers a situation similar to that of Figure 2. But the two situations are not statistically equivalent: in the original system (Fig. 2) all the periodic cracks experience an identical time evolution, whereas in the system of Figure 3 the time evolutions of the cracks in the primary cell are different. In particular, in the original system, failure occurs simultaneously at all the cracks, whereas in the system with multiple cracks, only one of the cracks in the primary cell fails. As a consequence, at same reduced crack length 2a/L, temperature, system size and loading rate, the estimated toughness differs between the two systems. This points to a specific scaling with respect to the number of initial cracks in the periodic cell that we address together with the other factors. In the simulations, we consider a number of crack ranging from 1 to 256 (i.e., 16 by 16).

In addition to the toy model, we study the more realistic situation of graphene. Graphene is a crystalline material with honeycomb structure (Fig. 1) which has attracted a lot of attention because of its exceptional physical properties. Graphene failure has been the focus of many studies in the literature (see, for instance, the review of [START_REF] Zhang | Fracture of graphene: a review[END_REF] and references therein). It appears that graphene is brittle, i.e., exhibits high strength (∼ 130 GPa) but relatively low toughness (∼ 4 MPa. √ m). Many molecular simulation studies [START_REF] Omeltchenko | Crack Front Propagation and Fracture in a Graphite Sheet: A Molecular-Dynamics Study on Parallel Computers[END_REF] wapriya et al ( 2014) to mention a few) investigated the failure properties of graphene from its tolerance to flaw to the effects of polycrystalinity. In the present paper, we consider graphene as a realistic case study in contrast with the case of the toy model.

Several inter-atomic potentials exist to simulate graphene. We consider the REBO potential of [START_REF] Brenner | Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films[END_REF]. This potential is known to lead to an excessive strength due a spurious bump in the aromatic carbon-carbon interaction force as the material is stretched [START_REF] Oa | Atomistic modeling of the fracture of polycrystalline diamond[END_REF][START_REF] Belytschko | Atomistic simulations of nanotube fracture[END_REF]. A way to circumvent this problem consists in modifying the switching function performing a smooth cut-off of the C-C interaction between r min = 1.7 Å and r max = 2.0 Å. The distance r min is set to 2.0 Å so that the function performs a sharp cut-off of the energy at 2.0 Å [START_REF] Yazdani | Failure criterion for graphene in biaxial loadinga molecular dynamics study[END_REF]. This modification effectively suppresses the bump and leads to an elasticity and strength that compares well with ab-initio calculations [START_REF] Liu | Ab initio calculation of ideal strength and phonon instability of graphene under tension[END_REF] and experimental measurements [START_REF] Lee | Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene[END_REF]. Note, however, that this modification introduces a discontinuity in the original formulation of the energy which must be handled carefully. In the present work, all the implementation of the molecular dynamics rely on the atomic forces only and the associated energy potential is the continuous integration of the force field. The general methodology is the same as for the simulation of the toy model except that we investigate tensile strength only. The graphene systems studied contain 2508 to 2618 atoms (about 8 × 8 nm 2 ) and are loaded with a strain rate of 0.001 ps -1 at temperatures ranging from 0K to 2700K. Even though graphene has a two dimensional structure, the failure behavior is affected by the third dimension. To investigate this question we perform both simulations in 3D and simulations in 2D by constraining the material in the plane.

T/T 0 =0.01 T/T 0 =0.03 T/T 0 =0.1 T/T 0 =0.3 T/T 0 =1 T/T 0 =3 T/T 0 =10 T/T 0 =30
Fig. 4 Tensile strength obtained for the toy model for various temperatures and loading directions. No shear is applied to the material so that σ xy = 0 for all the points represented.

All the calculations are performed with LAMMPS software (http://lammps.sandia.gov, [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF]).

Results

In this section we present the results of the molecular simulations. A detailed analysis of the scalings is presented in section 4.

Effect of temperature on strength

First of all, we focus on the effect of temperature on tensile strength. We show in Figure 4 the failure surface of the toy model at various temperatures. For sake of generality, the displayed results are reduced with respect to the characteristic tensile strength σ * 0 = σ 0 /K = r * cr -1. Moreover, the temperature scale is given relative to the characteristic failure temperature

T * 0 = k B T0 Kr 2 0 = (r * cr -1)
2 . At 0K, the theoretical failure surface is given by Equations ( 4). This theoretical surface is drawn with a black line in Figure 4. Note that we display the Piola I stresses (first Piola-Kirchhoff stress) and not the Cauchy stresses as is usually done. The reason is that the theoretical failure surface of Equations ( 4) is calculated for the undeformed configuration, whereas the Cauchy stress relates to the deformed configuration. In order to confront the results of molecular simulation to the 0K theory, we convert the Cauchy stress tensor (σ) into the Piola I stress tensor (π) as follows:

π = det F σ t • F -t (8)
where F is the deformation gradient tensor. In the limit of small strains, the Piola I stresses are almost equal to the Cauchy stresses. But, for the case of the toy model, the deformations are of the order of r * cr -1 = 10% which is large enough to induce significant differences between the two tensors. In terms of Piola I stresses, the simulations results at low temperatures are in very good agreement with the 0K theory.

One can readily interpret the shape of the tensile failure surface. It exhibits two modes of failure. The horizontal branch (π yy constant) corresponds to the failure of the bonds at ±30 • of the y axis in the crystal structure (Fig. 1). The oblique branch (3π xx -π yy constant) corresponds to the failure of the bonds along the x axis in the crystal structure. As temperature is increased, one observes a decrease of the tensile strength. Interestingly, this decrease appears almost independent of the loading direction.

In Figure 5, we display the tensile failure surface of graphene for various temperatures. As for the toy model, we represent the Piola stress which differs significantly from the Cauchy stress since the deformations reach ∼ 20% at failure. Again, the tensile strength decreases with temperature. However, this decrease is no more independent of the loading direction: the strength decreases faster in the zigzag ('zz') direction than in the armchair ('ac') direction. Therefore, the orientational invariance observed for the toy model is a peculiarity of this system. One can clearly identify two regimes in the tensile failure surface of graphene: an horizontal branch (π ac almost constant) corresponding to the failure of the bonds oriented in the armchair direction in the crystal structure (Fig. 1); and a vertical branch (π zz almost constant) corresponding to the failure of the bonds at ±30 • of the zigzag direction in the crystal structure. Of course, this bond failure interpretation is an idealization since the actual atomic interactions in

T/T 0 =1 T/T 0 =3 T/T 0 =10 T/T 0 =30 T/T 0 =0 T/T 0 =0.01 T/T 0 =0.03 T/T 0 =0.1 T/T 0 =0.3
Fig. 6 Failure stress of the pre-cracked systems in function of the initial crack length 2a/L at various temperatures. The marks are the simulation results, whereas the lines are the fits with fracture mechanics (Eq. 7) providing estimations of the toughness. For clarity, the curves are arbitrarily shifted vertically. In an exact representation, all the curves converge at 0 for 2a/L = 1. The 0K case is obtained by energy minimization and not by molecular dynamics.

graphene are many-body and much more complex than bond interactions.

Effect of temperature on toughness

In this section, we focus on the toughness estimation. We perform molecular simulations of pre-cracked systems as illustrated in Figure 2. When changing the initial crack length, one expects Equation (7) of fracture mechanics to apply, thus providing a way to estimate the toughness. We perform simulations at different temperatures and we display the stress at failure in function of the initial crack length in Figure 6. Then, we fit each curve with Equation ( 7) to derive the evolution of the toughness with the temperature. For sake of clarity, the results of Figure 6 are shifted vertically. In an actual representation, all curves would end at 0 for 2a/L = 1. Fits with the theory of fracture mechanics are precise, except at the largest temperatures were deviations are observed for small crack lengths. At 0K, we obtain a toughness K * Ic = K Ic K √ r0 = 0.198. This value compares well with the theoretical estimate (Eq. 5): K * Ic = 3π/4 (r * cr -1) = 0.153. The tensile strength in the same loading direction (y) at 0K is given by Eq. 4:

σ * strength = σ strength K = √ 3 (r * cr -1) = 0.1 √ 3.
Therefore the characteristic length of the process zone at 0K is r * pl = r pl r0 ≈ 1.31. As expected, the process zone is about the size of an atom, and fracture mechanics is expected to capture well the failure behavior even at the atomic scale. This is indeed observed in Figure 6. Nevertheless, at high temperature (T /T 0 = 30), the same calculation gives a much larger process zone: r * pl = 4.69. Hence the deviation from fracture mechanics at high temperatures.

With this method, we can estimate the evolution of the toughness with temperature. The next section confront the temperature scalings of strength and toughness.

Trends with temperature

We compare in Figure 7 the scalings of strength and toughness with temperature. The error bars have different meanings for the two quantities. For strength, it represent the standard deviation between the various loading directions (Fig. 4), whereas, for the toughness, it represents the standard deviation associated to the estimations from equation ( 7) at the various values of 2a/L. In particular, the error bars for toughness increases with the temperature because of the discrepancy discussed in the previous section. In the following of the paper, all errors bars refer to these definitions.

As temperature is increased, both strength and toughness decrease. Interestingly, the toughness decreases less rapidly than the strength. The ratio between toughness and strength is also the size of the process zone (r pl = (K Ic /σ strength )

2 ) which, therefore, increases significantly. Apart from this difference, strength and toughness exhibit very similar scalings. We draw empirical scaling laws in Figure 7 (dashed lines) which depends on the square root of temperature:

σ strength σ strength (0K) or K Ic K Ic (0K) = 1 - √ αT (9)
α is the single fitting parameter and takes a different value for strength and toughness. It is notable that the same expression captures both strength and toughness trends. It appears that this scaling is not linear as predicted by Zhurkov's theory (Eq. 1), but we show in the following section that this scaling law is consistent with Zhurkov's approach. Likewise, we display in Figure 8 the scaling of graphene strength with temperature. In contrast with the toy model, the decrease of strength depends on the orientation. So we distinguish the evolution for the armchair direction to that for the zigzag direction. We also distinguish the results for graphene constrained in 2D to that in 3D. We observe that strength decreases almost linearly with temperature, in both directions and irrespective of 2D constraint, as expected from Zhurkov's theory. The difference of scaling with the toy model shows that the scaling with temperature is closely related to the nature of the atomic interactions. Since the linear scaling is common in real materials, one can suspect that some aspects of the toy model are unrealistic. We discuss this point later on. The strength of graphene decreases more rapidly in 3D than in 2D. This is expected since constraining the material in 2D prevents transition paths to failure that would involve out-ofplane deformations of the atomic structure. The energy barrier to failure is expected to be larger in 2D than in 3D, hence the difference in temperature scaling.

Effect of loading rate and system size

According to the theory of thermally activated processes, the scaling of failure should combine temperature, loading rate and system size. Following Equation ( 1) with N T S degenerate transition states, the failure stress is expected to be a function of T ln N T S ε/ ε cr , with ε cr a constant of the dimension of a strain rate. One expects that this rate-temperature-size equivalence is characterized by a single master curve. So we completed the toy model scaling by performing a large series of simulations to investigate the effect of system size and loading rate, for both fracture failure and tensile failure. The results are displayed in Figure 9 for tensile failure and in Figure 10 for fracture failure.

We observe that both strength and toughness increase with the loading rate. However, while strength decreases with the number of atoms, the toughness appears almost independent of it. Our interpretation is that the number N T S of transition states scales with the number of atoms in a flawless material, and with the number of cracks in the case of a cracked material. Accordingly, we investigated the scaling with the number of cracks as explained in the methods (section 2). Doing so, we observe a decrease of the toughness with the number of initial cracks.

In Figure 9, we observe that all the strengths align very well on a single master curve validating the ratetemperature-size equivalence with the number of atoms as a proxy for the number of transition states. On the top chart, all the simulations are performed at the same temperature (T * = T * 0 = (r * cr -1)

2 ). A few results corresponding to large systems submitted to high loading rates deviates from the master curve. We attribute this discrepancy to insufficient equilibration that prevents sampling correctly the phase space of atomic configurations. On the bottom chart, all the simulations are performed with the same number of atoms (N = 270) and the strengths relative to 0K are plotted in function of the quantity ln N ε/ ε cr T T0 . By fitting appropriately the constant ε cr , the data points align very well along a master curve.

In Figure 10, we observe that the toughness results also align along a single master curve considering the number of cracks as a proxy for the number of transition states. The chart on top represents the toughness relative to 0K for systems simulated at the temperature T * = T * 0 . As for strength, we observe a deviation of the molecular simulation results from the scaling law at high loading rates, and the deviation increases with the system size. The deviation is more pronounced than for strength: for a system of given size, reaching the scaling law requires to consider loading rates about two orders of magnitude smaller than for the strength simulations. As before, we attribute this deviation to an insufficient equilibration which introduces a bias in the sampling of the phase space. Interestingly, cracked bodies seem to require much more time to equilibrate than flawless systems. The bottom chart in Figure 10 displays the results for systems of identical sizes with 1 crack Relative strength Relative strength

σ yield /σ yield ( T =0K ) T/T 0 =1 N =30 N =270 N =1080 N =4320 N =17280 N =69120
σ yield /σ yield ( T =0K ) N =270 T/T 0 =0.01 T/T 0 =0.1 T/T 0 =1 T/T 0 =10

Scaling law

Fig. 9 Scaling of the strength of the toy model with respect to the combination of system size, temperature and loading rate. The top chart displays the results of systems all at temperature T * = T * 0 = (r * cr -1) 2 but of various sizes and loading rates. The strength relative to 0K is plotted in function of the quantity ln (N/ ε * ). The bottom chart displays the results of systems all of the same size (N = 270) but at various temperatures and loading rates. The strength relative to 0K is plotted in function of the quantity ln

N ε/ ε cr T T 0
. In both chart, the results align very well on a general master curve which is accurately fitted by the proposed scaling law (Eq. 12).

at various temperatures and loading rates. We plot the toughness relative to 0K in function of the quantity ln N crack ε/ ε cr T T0 . Here also, with an appropriate choice of the constant ε cr , the results align well along a master curve. And, again, the results of systems submitted to the highest loading rate deviate from the scaling law.

Analysis and discussion

One can interpret the observed scaling by adopting a statistical physics description akin to Zhurkov's theory. We develop such an analysis in this section. Relative toughness Relative toughness

K Ic /K Ic ( T =0K ) T/T 0 =1 N crack =1 N crack =4 N crack =16 N crack =64 N crack =256
K Ic /K Ic ( T =0K ) N crack =1 T/T 0 =0.01 T/T 0 =0.1 T/T 0 =1 T/T 0 =10
Scaling law . The bottom chart for systems with N crack = 1 at various temperatures and loading rates. We plot the toughness relative to 0K in function of the quantity ln N crack ε/ ε cr T T 0 . In both chart, the results exhibit the scaling behavior which is well captured by the proposed scaling law (Eq. 12).

Toy model scaling

We consider failure as a thermally activated process with no reverse process. Following Boltzmann statistics, the probability of occurrence of a transition state leading to failure is proportional to exp

-∆E T S k B T
, where ∆E T S > 0 is the energy barrier that has to be overcome. In the toy model, failure occurs when one of the bonds has a length exceeding the critical distance r cr . At 0K, all atoms are in static equilibrium and failure occurs precisely at the theoretical surface (Eq. 4). At finite temperature though, failure can occur below this limit because of thermal agitation. Let us consider the material loaded below the theoretical surface (Eq. 4). The ground state (state of minimum energy) at this loading is the system at 0K. A transition states leading to fail-ure is such that one of the bonds reaches r = r cr . Since all bonds are harmonic, the energy barrier to reach the transition state is quadratic1 :

∆E T S = K 2 κ (r cr -r) 2 (10)
with κ > 0 a constant. This energy barrier can be rewritten in term of stress by taking advantage of the linearity of the mechanical behavior: the length of a bond depends linearly on the strain and stress. We have:

∆E T S = k B T cr 1 - σ f ailure σ f ailure (0K) 2 (11)
where k B T cr = K 2 κ (r cr -r 0 ) 2 = κ 2 k B T 0 is a caracteristic temperature of failure. Expression 11 is valid for both tensile failure and fracture failure. A difference appears regarding the number of transitions states. For tensile failure, by periodicity of the atomic structure, the number of transition states is proportional to the number of atoms. For fracture failure, it is proportional to the number of cracks. Therefore, the probability of failure which is the cumulative probability of all the degenerate transition states takes the form

∝ N T S exp -∆E T S k B T
, with N T S = N for tensile failure and N T S = N crack for fracture failure. This Boltzmann statistics can also be interpreted as the inverse of the characteristic time of recurrence of this event during a time evolution. Equivalently, the probability of failure over a given time scale is inversely proportional to the loading rate. We finally obtain the scaling relation:

σ f ailure σ f ailure (0K) = 1 -ln N T S ε/ ε cr T T cr (12)
where ε cr is a constant which has the inverse dimension of a time. The ratio σ f ailure σ f ailure (0K) is the ratio of strength for tensile failure and the ratio of toughness for fracture failure. This scaling law is precisely of the form of the empirical law we suggested in the previous section (Eq. 9). In Figures 9 and10

Graphene scaling

In the previous section, we derived a scaling law for the toy model in which the failure depends on the square root of temperature. This scaling differs markedly from the linear scaling we obtained for graphene and which is commonly observed for real materials (Eq. 1). This difference is a direct consequence of the peculiar expression of the energy barrier for the toy model (Eq. 11). This expression holds for perfectly harmonic atomic interactions (Eq. 2). This idealized potential is convenient for theoretical derivations such as that of the scaling law, but it exhibits a discontinuity in its derivative at the critical length r cr which, of course, is not representative of real atomic interactions. Instead, significant non linearity is expected in the atomic interactions when approaching failure. This is indeed the case for graphene simulated with the REBO potential. How this difference impacts the energy barrier ∆E T S is hard to anticipate since a formal demonstration is no more possible. Nevertheless, one can estimate numerically ∆E T S by searching for the minimum transition path between an initial configuration (intact material) and a final configuration (broken material).

We can estimate numerically the energy barrier for graphene failure. The minimum transition path is found by increasing the length of a bond while relaxing all other degrees of freedom. We display in Figure 11 the energy barrier in function of the stress for different loading directions (armchair, zigzag and spherical). The xintercept corresponds to the failure at 0K (σ f ailure (0K)).The trend is clearly different from the energy barrier of the toy model (Eq. 11). Instead the energy barrier depends almost linearly on the stress, that is:

∆E T S = k B T cr 1 - σ f ailure σ f ailure (0K) (13)
with T cr the temperature corresponding to the energy barrier at zero stress. Combining this expression with the statistics for thermally activated process, one recovers a linear scaling law for failure, consistent with the results of Figure 8:

σ f ailure σ f ailure (0K) = 1 -ln N T S ε/ ε cr T T cr (14)
This linear scaling is expected to be more relevant for real materials than that obtained for the toy model (Eq. 12) since a linear scaling with temperature is commonly observed experimentally [START_REF] Zhurkov | Kinetic concept of the strength of solids[END_REF].

Conclusion

The analysis of scaling provides a general rate-temperaturesize equivalence for both tensile failure and fracture failure. The main difference between strength and toughness stems to the number of transitions states N T S which scales with the number of atoms for tensile failure and with the number of cracks for fracture failure. Failure always originates from the occurrence of a transition state, but the distribution of transition states and the associated energy barriers is strongly affected by stress concentration and the presence of flaws. Tensile failure can occur at any location in the material, fracture failure is limited to a small number of transition states at the crack tip.

There are several important implications to this study:

-The rate-temperature-size equivalence relates the failure properties under very different conditions. One can easily relate failure at the atomistic level, typically ∼ 10 3 -10 6 atoms loaded at ∼ 10 10 -10 7 s -1 , and at the macroscopic level, typically ∼ 10 23 atoms loaded at rate ∼ 1s -1 . Proper confrontation of failure properties at different scales can be achieved by taking advantage of the rate-temperature-size equivalence.

-Knowing the effect of temperature on failure is enough to predict the scaling with loading rate and system size. This is of practical interest, in particular for molecular simulation studies, in which it is easy to vary temperature with little or no additional computational cost; whereas increasing the system size or decreasing the loading rate over orders of magnitudes comes with very significant computational costs, if not prohibitive.

-Stress concentration at a crack tip strongly localizes the transition states. The transition from fracture failure to tensile failure through flaws of more moderate stress concentrations (e.g., notch, holes etc.) involves the transition from a regime of very few possible transition states, to a regime of lots of transition states. How the number of transition states and their energy barrier depends on the stress concentration is thus an important question in the formalization of a theory of failure initiation.

The scaling law we propose here is not directly applicable to the prediction of failure of real materials since one would have to account for the presence of existing flaws of various length and density. Extending the type of scaling laws we propose to distribution of flaws is a perspective.

Fig. 1

 1 Fig. 1 Crystal structures of the toy model 2D solid (left) and of the graphene (right). The basis vectors give the name of the orientations used throughout this paper.

Fig. 2

 2 Fig.2System considered for the study of fracture failure: an initial crack is initiated by removing bonds along a line. Due to periodicity, the same crack exists in the periodic replica of the system.

  ;[START_REF] Belytschko | Atomistic simulations of nanotube fracture[END_REF];[START_REF] Khare | Coupled quantum mechanical/molecular mechanical modeling of the fracture of defective carbon nanotubes and graphene sheets[END_REF];[START_REF] Zhang | Flaw insensitive fracture in nanocrystalline graphene[END_REF];[START_REF] Moura | Tearing of free-standing graphene[END_REF]; De-

Fig. 3

 3 Fig.3Type of system considered to investigate the scaling of toughness with the number of cracks: several periodic cracks are initiated within the elementary cell actually simulated. Periodic boundary conditions makes this geometry similar to that of Figure2, but the two systems are not statistically equivalent; in particular, the corresponding toughnesses are different.

Fig. 5

 5 Fig. 5 Tensile strength results for graphene simulated in 2D (results in 3D are similar). All points are in the plane σ xy = 0. The results at 0K are obtained by minimization (static equilibrium) and not by molecular dynamics.

Fig. 7

 7 Fig. 7 Scaling of strength and toughness of the toy model with temperature. The error bars for strength correspond to the standard deviation associated to the loading direction. The error bars for the toughness correspond to the standard deviation associated to the initial crack length.

Fig. 8

 8 Fig.8Scaling of graphene strength with temperature. The strength in the armchair (ac) and zigzag (zz) directions are plotted separately. We also present results for the graphene simulated in 3D and for graphene constrained in 2D. The dashed lines are trend lines.

Fig. 10

 10 Fig. 10 Scaling of the toughness of the toy model with respect to the combination of number of cracks, temperature and loading rate. The top chart displays the results for systems simulated at T * = T * 0 with various number of cracks N crack and loading rates ˙ . We plot the toughness relative to 0K in function of the quantity ln N crack ε *

Fig. 11

 11 Fig. 11 Energy barrier to reach the transition state for graphene failing in the armchair direction. The linear trend suggests a relation ∆E T S ∝ -(σ -σ strength (0K)), hence a strength scaling σ strength σ strength (0K) = 1 -ln N T S ε/ ε cr T T cr consistent with the results of Figure 8.

& 8 avenue Blaise Pascal,

Marne-la-Vallée Champssur-Marne, France. K. Sab Laboratoire Navier, UMR 8205, École des Ponts, IFSTTAR, CNRS, UPE 6 & 8 avenue Blaise Pascal, 77455 Marne-la-Vallée Champssur-Marne, France.

∆E T S is the integration of the repelling force to exert on the atoms of the failing bond to bring it to failure, while maintaining all other atoms at static equilibrium. By linearity of the bond network, this repelling force depends linearly on the bond length, hence the quadratic energy barrier.

Acknowledgements We gratefully acknowledge funding from the Labex MMCD provided by the national program Investments for the Future of the French National Research Agency (ANR-11-LABX-022-01)


Jean Dobrowolski

Petr Korba
email: korb@zhaw.ch

Walter Sattinger
email: walter.sattinger@swissgrid.ch

Felix Rafael

Segundo Sevilla

Centralized Wide Area Damping Controller for Power System Oscillation Problems

Keywords: Inter-area Oscillation, Particle Swarm Optimization, Linear Quadratic Gaussian, Dynamic Stability, Modal Analysis

In this paper, three different centralized control designs that vary on complexity are presented to damp interarea oscillations in large power systems. All the controls are based on phasor measurements. The first two proposed architectures use simple proportional gains that consider availability of measurements from different areas of the system and fulfill different optimization functions. The third controller is based on a more sophisticated Linear Quadrartic Gausian approach which requires access to the state space model of the system under investigation. The novelty of the proposed scheme resides in designing a single control to command the most influence group of machines in the system. To illustrate the effectiveness of the proposed algorithms, simulations results in the IEEE New England model are presented.

I. INTRODUCTION

With the general trend of massive penetration of renewable sources in power systems and the expected decommissioning of nuclear based energy as in Switzerland [START_REF]Energy Strategy 2050[END_REF], the future evolution of power grids is a current topic of interest worldwide. Particularly, there is the question if future power systems should evolve in the form of micro-grids, which means to have multi-autonomous small geographical areas where load and generation are integrated locally or continuing developing to even larger super-grids. Each bifurcation has benefits and limitations, for instance micro-grids lack of long distance interconnections having as result absence of inter-area oscillations, however such systems are prone to frequency and other issues if the balance between generation and demand is not meet. On the other hand, super-grids are more robust networks and also significantly more complex to operate. In Europe, the interest for super-grids is evident as demonstrated with the interconnection of Turkey in 2010 to the already complex network of continental Europe [START_REF] Nassar | System analysis of the turkish power system for interconnection with continental europe[END_REF]. In addition, there is the The authors acknowledge the Swiss National Science Foundation (SNSF) under the program Ambizione Energy Grant (AEG). This research is also part of the activities of SCCER FURIES, which is financially supported by the Swiss Innovation Agency (Innosuisse-SCCER program). potential to grow the interconnection to include areas such as the Nordic region or the UK grid, profiting from the existing HVDC connections and building even more [START_REF]Trans-European Infrastructure -Projects of Common Interest[END_REF]. However, despite the robustness and reliability that extra large power systems provide, seldom events like those reported by ENTSO-E in December 2016 on continental Europe, trigger large interarea oscillation causing severe negative effects on the system [START_REF]Analysis of ce inter-area oscillations of 1st december 2016, entsoe[END_REF]. For this reason, damping of inter-area oscillations is a topic of interest within the transmission system operators in Europe, particularly in Switzerland as power system centrally located in the hart of Continental Europe.

Facing inter-area oscillations is not a new topic and most of the solutions proposed can be divided in two directions: improvement of Power System Stabilizer (PSS) by tuning [START_REF] Panda | Coordinated design of tcsc controller and pss employing particle swarm optimization technique[END_REF], [START_REF] Stativȃ | Optimal tuning and placement of power system stabilizer using particle swarm optimization algorithm[END_REF] or redesigning [START_REF] Dobrowolski | Inter-area Oscillation Control Based on Eigensystem Realization Approach[END_REF] its parameters and developing new central schemes such as Wide Area Damping Controllers (WADC) [START_REF] Zamani | Design of an h pid controller using particle swarm optimization[END_REF], [START_REF] Zhang | Design of wide-area damping controllers for interarea oscillations[END_REF] to command generators based on their participation in such events. Designing WADC can vary from simple PID to sophisticate H 2 /H ∞ controllers. In this paper, based on the assumption that readjusting parameters on every existing PSS in a large power grid is ideal but unrealistic, the effects of a new centralized WADC to damp inter-area oscillations is proposed. Three designs that go from a simple proportional gain to a more sophisticated Linear Quadratic Gaussian controller are proposed and explained in detail in section II. The methods for tuning the proposed algorithms and system identification are presented in section III and finally, section IV presents a validation of the methodologies.

II. PROPOSED CONTROL APPROACH

To deal with the problem of inter-area oscillations, as introduced in Section I, an algorithm based on the idea of adding a centralized WADC to the control loop, as initially suggested in [START_REF] Zhang | Design of wide-area damping controllers for interarea oscillations[END_REF], is presented. The proposed WADC receives measurements from Phasor Measurement Units (PMUs) and defines a new input for the excitation system, which is added to the existing controls as PSS. The proposed control sequence of one area, which is used in a multi-area grid is presented by Fig. 1. To illustrate how the proposed approach works, consider first a simple system with only two areas such as the New England IEEE 39 bus system [START_REF] Dutta | Optimal rescheduling of generators for congestion management based on particle swarm optimization[END_REF], or the widely known Kundur system models [START_REF] Klein | A fundamental study of inter-area oscillations in power systems[END_REF], respectively. Then, three control variations for the design of a WADC are derived based on different assumptions related to the availability of PMU measurements to feed the WADC system as follows and full details are given in the subsequent subsections:

A. Description of Control A (C A )

In the first control variation, the PMUs are measuring variables only in the area of interest. In this case, the WADC controller will be an additional control that minimizes the gap between the measured frequency and its reference during an event. To achieve this goal and at the same time deal with the original objective of simplification in power grid control, a simple proportional gain as WADC is used. The control architecture of variation A is presented on Fig. 2. A conven- tional approach for tuning the parameters k ai , i ∈ [1, n] in the WADC is doing manually, however, this is not straightforward nor simple task for the reason that each generator in the power system under investigation has different participation factor during inter-area oscillations. Thus, more sophisticated algorithms such as the Particle Swarm Optimization (PSO) are available to achieve this task. PSO is a widely used method for power systems regulators tuning [START_REF] Panda | Coordinated design of tcsc controller and pss employing particle swarm optimization technique[END_REF], [START_REF] Stativȃ | Optimal tuning and placement of power system stabilizer using particle swarm optimization algorithm[END_REF], [START_REF] Zamani | Design of an h pid controller using particle swarm optimization[END_REF], [START_REF] Gaing | A particle swarm optimization approach for optimum design of pid controller in avr system[END_REF] and therefore, in this work is used. An introduction and full overview of PSO approach is presented in Section III.

B. Description of Control B (C B )

In the second design, it is assumed that there are PMU measurements available from both areas. Thus, the control has the capability to minimize the differences between generators by minimizing the frequency variation between the two areas. To deal with the restrictions of (C A ) aforementioned in II-A, a proportional controller depicted on Fig. 3 is proposed. To find the most appropriate value of the parameters k bi , i ∈ [1, n], PSO approach is also used.

C. Description of Control C (C C )

In the third control variation, in addition to the assumptions described for (C B ), the following state-space representation of the system is also considered to be available:

x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) (1)
From this representation, a Linear Quadratic Regulator (LQR) controller is proposed based on its flexibility to adapt to the state space system (1). Since this control scheme requires the knowledge of the state vector x(k), which in the case of power systems correspond to generator speeds and angles as well as magnitudes and angles of bus voltages, the association with a Linear Quadratic Estimator (LQE) is required to form the so called Linear Quadratic Gaussian (LQG) control. As demonstrated in [START_REF] Dobrowolski | Inter-area Oscillation Control Based on Eigensystem Realization Approach[END_REF], the solution of the LQG controller has two components:

• One control sequence u = -Kx with K as the solution of the quadratic minimization problem;

• One Kalman gain G such that (2), for x(k) = x(k) then x(k + 1) = x(k + 1) else, lim n→+∞ (x(k + n) -x(k + n)) → 0 x(k + 1) = Ax(k) + Bu(k) + G(y(k) -ŷ(k)) ŷ(k) = C x(k) (2)
The block diagram of the third design of WADC is depicted in Fig. 4. Subsequent section III introduce the process for identification of matrices A, B and C, respectively.

III. METHODOLOGY

After the three different variations to designs a WADC have been defined, the PSO methodology and its application in power systems to tune controllers is presented. Then, subsection III-C presents the identification of a state-space representation of the interactions between two areas in the IEEE New England system. The model identification process is required to calculate the gains K and G introduced on Section II

A. Particle Swarm Optimization theory

The Particle Swarm Optimization is one of the Swarm method for solving optimization problems. It is an optimization technique, first presented by Eberhat and Kennedy in 1995 [START_REF] Kennedy | Particle swarm optimization[END_REF], inspired by the natural movement of bird and fish groups. The idea of this algorithm is to define a number of particles n where every particle represents a candidate solution for the defined optimization problem J, which is represented by an objective function to minimize in this case. Every particle i is characterized by three different factors:

1) Position P i and Velocity V i ; 2) A measure of quality J(P i );

3) A best solution since the launch of the algorithm P besti If all particles are considered, an extra parameter defined as P best is introduced, which correspond to the global or best possible position that satisfy the inequality J(P best ) ≤ J(P besti )∀i from the launch of the algorithm.

During the first iteration, particles are randomly placed in the searching space restricted by limitations in the position and a its velocity is randomly defined. Then, at each iteration k the position P i and velocity V i of the i th particle is updated using equation (3):

       V i (k + 1) = H × V i (k) +c 1 × r × (P i (k) -P besti ) +c 2 × r × (P i (k) -P best ) P i (k + 1) = P i (k) + V i (k + 1) (3)
where H is the inertia, c 1 and c 2 are weighting coefficients and r is a random number that changes value at each iteration. From ( 3) is possible to understand that PSO approach operates according to three principles as follows:

1) Momentum, who corresponds to the influence of the current motion related to the inertia H; 2) A cognitive part, which corresponds to the attraction to the particle with best solution P besti ; 3) A social section, related to the amount of attraction to the best position of the tribune (or swarm) P best . In [START_REF] Panda | Coordinated design of tcsc controller and pss employing particle swarm optimization technique[END_REF], a good representation of the PSO behavior is proposed as depicted on Fig. 5. Since PSO formulation has some random variables, its performance depends mainly on the number of particles defined. Indeed, the more particles are used, the more space will be cover after the first iteration and the more likely it is to find the global optimum. Nevertheless, is also worth mentioning that the more particles used, the longer will take to the algorithm to converge. One simulation per particle is required at each iteration and thus, more computational time demanded. To ensure the stability of the algorithm, parameters c 1 and c 2 must be appropriately tuned to control the velocity of the particles.

B. Particle Swarm Optimization applied in power systems

PSO has been widely used for power systems controls tuning [START_REF] Panda | Coordinated design of tcsc controller and pss employing particle swarm optimization technique[END_REF], [START_REF] Stativȃ | Optimal tuning and placement of power system stabilizer using particle swarm optimization algorithm[END_REF], [START_REF] Zamani | Design of an h pid controller using particle swarm optimization[END_REF], [START_REF] Gaing | A particle swarm optimization approach for optimum design of pid controller in avr system[END_REF]. Its convergence and stability in case of multidimensional space has been proven in [START_REF] Clerc | The particle swarm-explosion, stability, and convergence in a multidimensional complex space[END_REF]. For this reason, to tune the proposed controllers defined in Sections II-A and II-B, respectively, PSO is applied defining every particle as a n dimensional particle (with n positions) corresponding to the n parameters k (a/b)i , i ∈ [1, n]. In this form, every particle corresponds to a WADC configuration. The particularity of PSO and the quality of its results arises from the definition of the measure of quality J(P i ) introduced on Section III-A.

Considering as a main priority improving the effects of low frequency (0.1-0.9 Hz) inter-area oscillations in the system [START_REF] Kundur | Power system stability and control[END_REF], two different methods for the definition of optimization function J are considered:

1) First, as proposed in [START_REF] Stativȃ | Optimal tuning and placement of power system stabilizer using particle swarm optimization algorithm[END_REF], the eigenvalue (λ i = σ i ± jω) based multi-objective function defined in equation ( 4) is considered. Given that improving the damping (ξ) for a particular mode may cause a negative effect on other sensitives modes, the group Θ denotes the critical modes (within frequency range of interest), which are selected following a simulation without WADC. The final objective is then to drive the selected eigenvalues towards a desired location in the complex plain, represented by shifting the real part σ 0 of the eigenvalue λ 0 and consequently improving its damping ratio ξ 0 .

J 1 = i∈Θ [σ 0 -σ i ] 2 + β i∈Θ [ξ 0 -ξ i ] 2 (4)
where ξ i the damping ratio of the i th eigenvalue and β a scaling factor. 2) The second method is based on measuring the quality of the actual measurements. Since the proposed controls C A and C B are based on generator speed measurements, the objective function in equation ( 5) is used to minimize the differences between generators rotor speed (ω r ) in the two areas while maintaining the regulation to the nominal value ω ref and minimizing the number of oscillations N osc .

J 2 = i∈a2 [ω a1 -ω i ] 2 + i∈a2 [ω i -ω ref ] 2 + βN osc (5)
where ω a1 is the average speed in area one and β a scaling factor.

C. Identification of inter-area interactions

In this section, identification of the state-space representation of the interaction of the inter-area oscillations is proposed. In our previous work [START_REF] Dobrowolski | Inter-area Oscillation Control Based on Eigensystem Realization Approach[END_REF] it has been proved that accurate dynamic identification of the system allows to design an adaptive control and then to improve the damping of interarea oscillations. Due to the fact that in this work power system simulations have been performed using the commercial software DigSilent PowerFactory, the proposed adaptive control requires co-simulation with Matlab, which has been proved to be computational expensive, specially for large power systems. Considering that the dynamic identification approach presented in [START_REF] Dobrowolski | Inter-area Oscillation Control Based on Eigensystem Realization Approach[END_REF], which is referred as Eigensystem Realization Algorithm (ERA), has the drawback of relying on an impulse as control input. For this reason, in this work a more conventional identification method to design the third control architecture introduced on Section II-C was selected. In order to find an good representation of the inter-area oscillation interactions, one possibility is to consider the state-space (1) with y ∈ n corresponding to variations between generator rotor speeds of the n generators in area 2 with the average speed of area 1 and u ∈ n as the PSS output, corresponding to u pss in Fig. 1. To achieve this goal, it is possible to use statespace identification subroutines such as the MATLAB/N4SID.

The IEEE-39 bus system depicted on Fig. 6 has 5 inter-area modes in the range of [0.59 -0.97] Hz. In this work, special attention is placed to the mode with smallest frequency since it is the weakest and represents the interaction between G1 (Area 1) gainst rest of the generators (area 2). In this system, the following event was simulated: A 29 sec simulation was performed, where at 1 sec the line 16 was suddenly tripped and the WADC controller was disable. After performing the simulation, the frequency of G1 was compared against the frequency of the nine remaining machines and the vector y(k) in ( 1) was formed. In the commercial software PowerFactory, is also possible to have access to the PSSs output and build the control vector u(k). Since there is no information a priory about the actual order of the system, a model reduction process is performed and the results are depicted on Fig. 7-(a). From this subplot, it can be observed that a new system of order 20 provides the lowest error between the reduced model and the original system. Figs. 7-(b-j) are a comparison between the simulation of the estimated state-space and the real measurement in the rest of the generators.

IV. POWER SYSTEMS APPLICATION & MAIN RESULTS

To demonstrate the effectiveness of the proposed control schemes presented on Section II, the IEEE New England system depicted on Fig. 6 has been used. The power system has been subject to a randomly selected 3-phase short-circuit at line 16 after 1 sec of simulation and the fault was cleared after 180 msec. Note that the proposed controls C A , C B , and C C have different methodology and thus, the software configuration depends of the control to be applied.

A. Performance Evaluation of Controls

C A & C B
To find the optimal parameters for controls C A and C B , respectively, PSO method was applied. To integrate the control objectives into the DigSilent PowerFacrtory environment, the Application Programming Interface (API) that link python with the power system software was operated. The principle is to launch a python code that controls parameters within PowerFactory whenever is requited, and at the same time to have access into the stored results. Hence, a subroutine for a application of PSO with p number of particles for i iterations was developed and the flow chart of the subroutine is depicted on Fig. 8. In this figure, blue instructions represent special actions for the evaluation of the optimization function J 1 , while red instructions represent special actions for the evaluation of the optimization function J 2 and black instructions are common actions for both optimization functions . From 250 tests, the optimal configuration for tuning of variables in equation ( 3) was found and the results are summarized on table I. Using this configuration, the tuning of the proposed controls C A and C B for implementation of the objective functions J 1 and J 2 , was carried on and the final results are depicted on Table II. For comparisons purposes, the value of the objective functions without WADC are J 1 = 550.651 and J 2 = 62154.2, respectively. The control C C has fewer constraints than controllers C A and C B and therefore no co-simulations are required, which speed-up simulation time and hence only PowerFactory is required. The main challenge related to this controller is a limitation from the power system software related to how matrix calculations are handled. To overcome this issue, the reduced order model of the IEEE New England system of 20 states was built from scratch as a new DSL model within DigSilent Power Factory [START_REF]Power Factory Users Manual[END_REF].

C. Main Results

In this section, the performance of the proposed controllers are presented for different optimization functions (OF -J i ), as presented on Table II. In total, five possible combinations using the three different controllers were investigated:

• Case study 1: C A for OF -J 1 ;

• Case study 2: C A for OF -J 2 ;

• Case study 3: C B for OF -J 1 ;

• Case study 4: C B for OF -J 2 ;

• Case study 5: C C . To compare the performance of the proposed controllers under different case studies, generators speed measurements from both areas in the IEEE New England system during RMS simulations are presented. For sake of simplicity, Fig. 9 (a)-(d) compares the speed variation of one generator from each area: G1 vs G6 (randomly selected) following a 3-phase fault of 180 msec applied to line 16 as in the previous case. Fig. 9 (a) depicts the simulation result without WADC. Fig. 9 (b) presents the results of the case studies 1 and 3. Fig. 9 (c) show the simulation results for case studies 2 and 4 and finally, Fig. 9 (d) depicts the simulation results for case study 5, which is the C C LQG controller. From the results presented on Fig. 9, the following remarks are derived:

• When no control is used, the inter-area mode with a frequency of 0.59 Hz is clearly visible in the response of the generators, which generates active power fluctuation in the order of 500 MW that can severely affect the objectives defined in the optimizing problem. • From subplot 9 (b) it can be observed that C A minimize the amplitude of the oscillations, but do not improve the stability of the inter-area mode in G1. The measure of J 1 for C A in Table II confirm this, with no significant improvement compared to the case without WADC. On the other hand, C B presents significant improvement, damping the oscillation fast (approximately 2.5 seconds after the fault). Since J 1 is designed only to improve the stability of selected eigenvalues, C B has no internal characteristics to automatically reduce the amplitude of the oscillations and hence the final amplitudes are not significantly reduced for this example. • From subplot 9 (c), it can be seen that the results are similar to subplot 9 (b). The main difference is that, the objective function J 2 described on equation ( 5), includes a component to minimize the amplitude. For this reason C A and C B provide better performance than in the previous case. Another remark is that the first element of J 2 , is used to minimize the difference between generator speeds in both areas and based in the simulation results, it can be concluded that this effect is as efficient as the eigenvalue based multi-objective function. • Subplot 9 (d) demonstrates the effectivenes of C C for case study 5. The controller reduce the amplitude of the oscillations and provide a significant damping in about 3 seconds.

The proposed controllers present different behavior from which C A provides the weakest performance to damp interarea oscillations. C A can positively influence the amplitudes of the oscillations, but cannot improve the stability of the system, which is the main objective in this work. On the other hand, C B and C C are the most efficient control schemes in terms of minimize oscillation amplitude while significantly improving the damping of the system, the challenge to implement the proposed architectures is that they require global measurement and a global communication with PMUs, which can be prone to communication problems.

V. CONCLUSIONS

Three different designs for WADC to damp inter-area oscillations have been presented. The first two architectures are based on proportional gains that consider availability of PMU measurements from different areas of the system and fulfill different optimization functions. The third controller is based on a more sophisticated Linear Quadrartic Gausian approach, which requires access to the state space model of the system under investigation. Performance of control C A indicates that measurements from one area are not sufficient to damp interarea oscillations. Similarly, simulation results based on C B and C C demonstrate that with measurements from more than one area is possible to significantly improve damping of interare modes. Moreover, from the results obtained with control C B based on OF -J 1 and OF -J 2 it can be noticed that is possible to enhance stability of the system as good as with approaches based on modal analysis.

Fig. 1 .

 1 Fig. 1. Control sequence of an area with n Generators (Gen) in a multiple area power systems. The Average Voltage Regulator (AVR) includes the excitation system. The Governor (GOV) is in charge of speed regulation.

•

  Control A (C A ): Measurements from only one area; • Control B (C B ): Measurements from both areas; • Control C (C C ): Measurements from both areas and a model of the interactions between the two areas is available.

Fig. 2 .

 2 Fig. 2. Design A of WADC. The PMUs are measuring variables of the area of interest and the WADC is an additional control to stabilize frequency (speed).

Fig. 3 .

 3 Fig. 3. Design B of WADC. The PMUs are measuring variables in both areas and the WADC is an additional control to minimize the gap between frequencies (speed) of both areas.

Fig. 4 .

 4 Fig. 4. Design C of WADC. The PMUs are measuring variables in both areas and a state-space representation of interactions between the two areas is available. The WADC is based on the Linear Quadratic Gaussian theory, K is the regulator gain and G is the Kalman gain to estimate states x.

Fig. 5 .

 5 Fig. 5. Representation of Particle Swarm Optimization behavior from [5]. Every particle move following three part: momentum, cognitive and social

Fig. 6 . 3 Fig. 7 .

 637 Fig. 6. Representation of Particle Swarm Optimization behavior from [5]. Every particle move following three part: momentum, cognitive and social

Fig. 8 .

 8 Fig.8. Subroutine to apply PSO with python and PowerFactory. The blue actions are specials for J 1 evaluation, the red ones for J 2 and the black one are common actions.

Fig. 9 .

 9 Fig. 9. (a) No control, (b) C A and C B to minimize OF -J 1 , (c) C A and C B to minimize OF -J 2 and (d) LQG control. M 1 and M 2 are respectively maximum amplitude of ω G 6 and ω G 1 .

TABLE II TUNING

 II OF CONTROLS C A AND C B FOR DIFFERENT OBJECTIVE FUNCTIONS (OF)

	k (a/b) i	C A	C B	C A	C B
	OF-J i	OF-J 1	OF-J 1	OF-J 2	OF-J 2
	k (a/b) 1	0.01	94.600	68.804	91.231
	k (a/b) 2	20.768	25.736	37.954	29.245
	k (a/b) 3	39.475	97.177	20.969	96.021
	k (a/b) 4	93.687	58.6967	23.569	49.568
	k (a/b) 5	49.713	6.841	57.285	10.254
	k (a/b) 6	79.819	85.303	1.116	82.187
	k (a/b) 7	26.635	92.851	0.090	96.329
	k (a/b) 8	78.989	11.678	0.01	9.548
	k (a/b) 9	87.928	72.948	11.857	68.245
	OF	516.885 274.804 50857.4 40533.6
	B. Performance Evaluation of Controls C C


Nicolas Pinchon

Olivier Cassignol

Frédéric Bernardin

Adrien Nicolas

Patrick Leduc

Jean Philippe Tarel

Roland Bremond

Emmanuel Bercier

Johann Brunet

J-P Tarel

R Brémond

All-weather vision for automotive safety: which spectral band?

Introduction

In the automotive industry, New Car Assessment Programs (NCAP) are increasingly pushing car manufacturers to improve performances of Advanced Driver Assistance Systems (ADAS), and especially autonomous emergency braking on vulnerable road users (VRU). For instance the 2018 Euro NCAP roadmap is moving towards pedestrians and pedal cyclists detection in day and night conditions. This trend matches accidentology figures, like those provided by French Road Safety Observatory (Table 1). 1. 2014 French accidentology data in adverse conditions [START_REF] Onisr | Les accidents corporels de la circulation 2014 -Recueil de données brutes[END_REF] In the longer term, after automated parking and highway driving, all weather and city driving will be the main technical challenge in the automated driving roadmap.

Injury casualties Fatalities

Current ADAS sensors as visible cameras or Lidars are fitting functional requirements of VRU and obstacle detection in normal conditions (day or night). However, these technologies show limited performances in adverse weather conditions such as fog or rain.

Automotive industry is thus facing this new challenge of detecting vehicle environment in all conditions, and especially in poor visibility conditions, such as night, fog, rain and snow. This topic has been addressed in the framework of the AWARE French public funded project, aiming at the development of a sensor enabling a vision in all poor visibility conditions. This paper presents an experimental comparative study of four different spectral bands: Visible RGB, Near-Infrared (NIR), Short-Wave Infrared (SWIR) and Long-Wave Infrared (LWIR). Sensors and field tests are described in section 2 and 3. Experimental results are detailed in section 4, focusing on pedestrians, vehicles, traffic signs and lanes detection.

Sensors

In this project, we only focus on cameras technologies, and not on distance measurement systems like LIDARs or RADARs. But it is well-known that both technologies are complementary and necessary to bring redundancy for improving the detection system's reliability and accuracy [START_REF] Premebida | Lidar and vision-based pedestrian detection system[END_REF]. Four cameras were tested during the project. The LWIR sensor is an array of microbolometers. It detects the thermal radiation in the spectral band extending from 8 µm to 14 µm. Any object emits radiations which depend on its temperature. For a human or an animal at ambient temperature, the maximum of emission corresponds to a wavelength close to 10µm. LWIR is used for the detection of a temperature contrast and do not require an illuminator. Mid-Wave Infrared (MWIR) has not been added in the study for reasons of cost and compacity, due to the cooling system required for the detectors.

Field tests

Outdoor test campaign

The four cameras were installed during one month in 2015 along the French motorway A75/E11 at the site named La Fageole. A weather station located near the cameras was equipped with a diffusiometer (meteorological visibility), a rain gauge, a luxmeter (ambient light) and a temperature and humidity sensor (ambient air).

Tunnel test campaign

The tests were carried in the CEREMA fog tunnel in Clermont-Ferrand (30m length, 5.5m width and 2.5m height). Artificial fog and rain are reproduced and controlled: fog and rain drop size, meteorological visibility of fog and rain intensity. Two fog classes are available: unimodal drop size distribution (DSP) centred around 1 micron and bimodal DSP centred around 1.5 and 10 microns. 74 scenarios were performed by varying the type of adverse weather conditions (fog and rain), scene illumination (night, day, automotive lighting, aircraft lighting) and glare (front vehicle). The scene in front of the cameras contained road pavement, road marking, pedestrians and lights (see Fig. 2).

Fig. 2.

CEREMA fog tunnel

Field tests

In this section, we describe the detection and recognition range performances that were measured in this study. It is important to keep in mind that these results reflect not only the intrinsic characteristics of the spectral bands but also the capability of the chosen cameras. The cameras were selected to be representative of the typical current state-of-the-art.

In order to prevent any detection algorithm artefact, a visual analysis has been performed by two different human observers. As expected, exact detection range values differed from one observer to the other, but the relative values were consistent. In all cases, brightness and contrast were carefully tuned in order to optimize ranges. With respect of each camera's channels created by the four spectral bands, a video database has been created to remotely record videos of relevant scenes for each listed scenario. The Fig. 3 provides a sample of the video database (outdoor campaign): Objects moving in the fog generate important transmission inhomogeneities. In order to avoid errors due to this effect, we recorded films of human test subjects standing still at the end of the tunnel while the fog cleared over time. Detection was declared successful when the outline of the chest of the test subject became visible against the background. This way, detection ranges were measured at the height at which the transmission meter was set. The test subjects were ~ 25 m away from the cameras. One of the test subjects wore a high visibility jacket and another one wore dark clothes. As expected, Visible, NIR and SWIR detection performances were better for the subject wearing the high visibility jacket (even though this improvement is less pronounced for thicker fog). For this study, the case of the subject in dark clothes was deemed more relevant.

The following table gives the fog density, expressed as standard visibility ranges, at which the pedestrian becomes visible. A reduced visibility range indicates a successful pedestrian detection in a thicker fog, and hence a better capability to see through fog. Cases with glare are not included. Error ranges mostly reflect the dispersion between the different scenarios used in the study. Conclusions are the following: ► The LWIR camera has a better capability to see through fog than the NIR and SWIR ones. The visible camera has the lowest fog piercing capability. ► The LWIR camera is the only one that allows pedestrian detection in full darkness. ► The LWIR camera also proved more resilient to glare caused by facing headlamps in the fog. Other cameras sometimes missed a pedestrian because she or he was hidden by the glare.

Camera Fog density for pedestrian detection

Examples of detection and recognition

Distances were calibrated within the cameras' field of view using the road markings (which are clearly visible for all cameras in images recorded on a sunny day) and an aerial map of the area. That method gives a distance measurement precision on the order of a few meters simply by noting the vehicle position within an image. The maximal distance that could be reliably assessed by this method was on the order of 150 m. In VIS, NIR or SWIR, detection was performed using the vehicle headlamps. The case of a vehicle driving with headlamps off while in adverse conditions was not encountered in this study. Should it happen, however, detection ranges in VIS, NIR and SWIR would be on the order of the recognition ranges. In LWIR, detection relied on the observation of hot vehicle parts: the wheels, the motor or the exhaust system. In some foggy instances, recognition proved difficult or even impossible in VIS, NIR or SWIR because the vehicle remained entirely hidden by the glare of its own headlamps. When this happened, the vehicles were not taken into account in the average ranges given in Fig. 9. Images illustrating this phenomenon are given in the Fig. 10. SWIR here presents more glare than the other bands but it is only due to the camera settings, it is not an intrinsic characteristic of the spectral band. The visual observation of the videos also confirmed the well-known fact that the exploitation of movement by the human visual system greatly increases detection capabilities: the success of the detection task is much higher when performed while watching a film than when performed by observing still images taken from the very same film. This is due to the human visual cortex implementing advanced spatiotemporal denoising and should inspire detection software developers.

In some recordings, wild animals are visible on the side of the road in the LWIR band. These animals are visible in none of the other spectral bands.

The conclusions are the following: ► For vehicles with headlights on, detection in adverse conditions is better in SWIR. NIR comes next and then VIS and LWIR. ► For the recognition task in the same conditions, the conclusions are inverted. ► The VIS, NIR and SWIR bands are sensitive to glare in foggy conditions, making the recognition task impossible in many cases. ► For both tasks, the LWIR camera gives much more reproducible results than all the others. In particular, its performance is independent of the vehicle's headlights being on or off. ► The LWIR camera was also the only one allowing the detection of hot blooded animals on the side of the road in adverse conditions.

Road marking detection

Detection ranges of road markings were evaluated in the outdoor test campaign by using the road lines of the highway. The same calibration as for vehicle detection was used to measure the distances. The Table 5 gives the average maximal detection distance for the four scenarios of interest: Conclusions are the following: ► Observation in the LWIR spectral band is only depending on road marking thermal emissivity. In 15 of the 21 scenarios recorded in LWIR, road markings are not visible. The visibility depends on the weather: rain is cleaning the lines while sun exposure enhances them (see Fig. 11). Thus road marking observation in the LWIR spectral band is not relevant. ► Visible, NIR and SWIR all required additional lighting to detect road marking at night. Detection in NIR and SWIR are equivalent, and slightly better than VIS. This is due in part to: a broader overall spectral band, monochromaticity (no RGB filter), larger pixels and a higher bit depth.

Scenario

Traffic signs recognition

The traffic sign located at 38 m from the cameras on La Fageole test site has been used to realize the comparative study. The Figure 12 shows the sign in two different weather conditions. Conclusions of the analysis of all the scenarios are the following: ► The SWIR sensor does not allow the identification of the traffic sign, even during daylight. This is due to the fact that the difference of reflectivity between the letters and the background is very low in the SWIR band. Traffic signs are indeed designed for visible spectrum and not for the other bands. ► The NIR camera provides a good vision and a better SNR than the RGB visible camera, especially in adverse weather conditions. This is due to the same reasons as for road markings detection. ► As expected, the traffic sign is never identifiable in the LWIR band. Letters and background have indeed the same temperature and emissivity.

Conclusion

As they are complementary to other distance measurement systems like LIDARs or RADARs, the AWARE project only focused its experiments on cameras technologies that are necessary to bring redundancy and complementary characteristics for improving ADAS detection system's reliability and accuracy.

In order to detect pedestrian, vehicle, road marking or recognize traffic signs, the relevance of four different spectral bands has been evaluated under adverse weather conditions: Visible RGB, Near-Infrared (NIR), Short-Wave Infrared (SWIR) and Long-Wave Infrared (LWIR).

The Table 6 Results of experiments clearly state that: ► In addition to the visible spectral band, only the LWIR spectral band provides outstanding interests. Targets can be detected with or without additional light. LWIR detection is not sensitive to any dazzle.

► NIR and SWIR provide equivalent performances, mainly due to the fact that experiments used extended visible to NIR and visible to SWIR cameras. However, as scattering in fog is reduced when wavelength increases, a NIR-only or a SWIR-only camera would certainly have provided much better performances than extended spectral bands. ► Visible RGB extended to NIR (or Red-Clear sensors) combined with LWIR provide the best spectral bands combination to improve ADAS performances of detection such as vehicle, pedestrian, bicycle, animals or road marking, and recognition such as traffic signs. ► Caution shall be considered while using LED headlight technology to provide additional light. LED pulsed technology could reduce detection reliability of system based on Visible, NIR and SWIR cameras.

Fig. 1 .

 1 Fig. 1.Outdoor test campaign More than 33 different scenarios have been set depending on weather conditions (rain type and intensity, ambient temperature), ambient light, and human eye visibility distance. Table3shows the most interesting recorded scenarios, with a class 3 fog (visibility distance < 100 m):

Fig. 3 .

 3 Fig. 3. Example of snapshots recorded by LWIR (top left), Visible (top right), SWIR (bottom left) and NIR (bottom right) cameras

Fig. 4 .

 4 Fig. 4. Pedestrian detection test setup into Cerema fog tunnel, LWIR (top left), Visible (top right), SWIR (bottom left) and NIR (bottom right) cameras

  Visible RGB Moderate (visibility range = 47 ± 10 m) Extended NIR High (visibility range = 28 ± 7 m) Extended SWIR High (visibility range = 25 ± 3 m) LWIR Extreme (visibility range = 15 ± 4 m)

Fig. 5 .

 5 Fig. 5.Example of images recorded in the fog tunnel with the four different cameras

Fig. 7 .

 7 Fig. 7.Reference of range based on site map analysis and T1 lanes type of road marking

Fig. 8 .Fig. 9 .

 89 Fig. 8.Vehicle detection ranges

Fig. 10 .

 10 Fig. 10. Images showing the glare effect in the VIS, NIR and SWIR spectral bands

Fig. 11 .

 11 Fig. 11. Road marking observation in the LWIR spectral band during day sunny condition (left) and day rainy condition (right)

Fig. 12 .

 12 Fig. 12. Images of traffic sign acquired in different spectral bands. In the first line, weather conditions is a sunny day and in the second line is from scenario 7 (fog class 3 with snow)

Table

	Night	32%	41%
	Wet road	20%	20%
	Adverse weather	21%	23%

Table 2

 2

	Camera	Sensor	Spectral band	Optics
		CMOS -SXGA		HFOV = 54°
	Visible RGB	(1280 x 966)	0.4-0.65 µm	VFOV = 40°
		3x8 bits, pitch 4.2 μm		F-number = 2
		CMOS -SXGA		HFOV = 39°
	Extended NIR	(1280 x 1024)	0.4-1 µm	VFOV = 31°
		10 bits, pitch 5.3 μm		F-number = 2.9
	Extended SWIR	InGaAs -VGA (640 x 512) 14 bits, pitch 15 μm	0.6-1.7 µm	HFOV = 39° VFOV = 31° F-number = 1.8
		Microbolometer -VGA		HFOV = 44°
	LWIR	(640 x 482)	8-12 µm	VFOV = 33°
		14 bits, pitch 17 μm		F-number = 1.2

shows their characteristics. The visible RGB CMOS camera is used here as a reference for the test. Extended NIR camera uses monochrome CMOS photodiodes with a cut-off wavelength close to 1 µm. It detects the reflective visible and NIR light from the scene. It thus requires an illumination by sun, moon or night glow or an illuminator positioned on the vehicle.

Table 2 .

 2 Cameras characteristics

	Extended SWIR camera is based on InGaAs III-V material and extends from a
	wavelength of 0.6 µm, red to human eye, to 1.7 µm in the SWIR infrared band.
	SWIR spectral band is typically used for active (reflective) vision in very dark
	condition with a good contrast as SWIR light is generally more reflective than vis-
	ible light.

Table 3

 3

	Scenario	Weather	Visibility distance	Ambient light	Tempe-rature
	7	Day, heavy fog, light snow	75 m	5032 lux	-1°C
	11	Day, heavy fog and snow	69 m	1096 lux	-1°C
	22	Day low light, heavy fog and snow	75 m	104 lux	-0.8°C
	26	Night, heavy fog	99 m	0 lux	+1.5°C

shows the most interesting recorded scenarios, with a class 3 fog (visibility distance < 100 m):

Table 3 .

 3 Most interesting recorded outdoor scenarios

Table 4 .

 4 Fog thickness for pedestrian detection at 25 m with the different cameras

Table 5 .

 5 Road lines detection ranges

Table 6 .

 6 and Table7summarize the relative performances of ADAS function for the four different cameras. Comparison of performances of ADAS functions using Visible, NIR, SWIR or LWIR camera technologies, Night fog, headlights

			Camera Spectral band
	ADAS FUNCTION		Night Fog, using headlights
		VIS	NIR	SWIR	LWIR
	Pedestrian, Bicycles, Animals detection	+	++	++	++++
	Vehicle shape recognition	-	-	--	++
	Vehicle lights detection	+	++	+++	--
	Traffic signs recognition	+	++	-	--
	Road marking detection	+	++	++	-
			Camera Spectral band
	ADAS FUNCTION		Day Fog, using headlights
		VIS	NIR	SWIR	LWIR
	Pedestrian, Bicycles, Animals detection	+	++	++	++++
	Vehicle shape recognition	+	++	++	++
	Traffic signs recognition	+	++	-	--
	Road marking detection	+	++	++	-

Table 7 .

 7 Comparison of performances of ADAS functions using Visible, NIR, SWIR or LWIR camera technologies, Day fog, headlights

Acknowledgement

The authors acknowledge the contribution of their colleagues to this work: P. Morange, J-L. Bicard and all the pedestrians from CEREMA, A. Picard from Sagem and B. Yahiaoui from Nexyad.