Xavier Brajer
email: xavier.brajer@saint-gobain.com

François Hild
email: francois.hild@lmt.ens-cachan.fr

Stéphane Roux
email: stephane.roux@lmt.ens-cachan.fr

ON THE DYNAMIC FRAGMENTATION OF GLASS: A MESO-DAMAGE MODEL

Keywords: FE simulations, fragmentation regime, non-local damage, Poisson-Weibull model

This paper presents an anisotropic damage model to deal with the fragmentation induced by impact loadings on glass samples. As small-scale (i.e., sub-element) damage is described as well as cracks extending above the element scale, an approach referred to as "meso-damage" is developed. The latter, which is based on the knowledge of random distributions of initiation sites, predicts different regimes such as single or multiple fragmentation. The experimental opening crack pattern obtained in edge-onimpact test is reproduced numerically.

Introduction

Impact resistant materials are being developed for armor applications [START_REF] Cagnoux | Déformation et ruine d'un verre pyrex soumis à un choc intense : étude expérimentale et modélisation du comportement[END_REF][START_REF] Den Reijer | Impact on Ceramic Faced Armor[END_REF]. For example, glass is used against ballistic threats (e.g., in windshields). For the latter application, multi-layered armor materials are used (with polycarbonate, PC, as back layer and polyurethane, PU, between glass plies) to optimize the ballistic performance. The use of PU and PC is essential for maintaining the fragments when glass is impacted, especially for the back layer that ensures the structural integrity and absorbs part of the impact energy. The work presented herein aims at modeling the fragmentation of sodalime silicate glass impacted by a bullet of .44 Magnum type. Edge-on-impact (EOI) configurations [START_REF] Hornemann | Experimental Investigation of Wave and Fracture Propagation in Glass -Slabs Loaded by Steel Cylinders at High Impact Velocities[END_REF][START_REF] Brajer | Behavior of soda-lime glasses impacted by a soft bullet: experimental and numerical investigations[END_REF][START_REF] Brajer | The role of surface and volume defects in the fracture of glass under quasi-static and dynamic loadings[END_REF][START_REF] Strassburger | High-speed transmission shadowgraphic and dynamic photoelasticity study of stress wave and impact damage propagation in transparent materials and laminates using the edge-on impact method[END_REF] allow for the visualization of cracking during the experiment by using a high-speed camera, or thereafter when a sarcophagus configuration is utilized [START_REF] Denoual | Analysis of the Degradation Mechanisms in an Impacted Ceramic[END_REF]. Different cracking patterns are observed. The present work intends to describe numerically the cracking features observed experimentally.

To predict crack patterns in brittle materials dynamic loading conditions, several approaches are followed [START_REF] Ravi-Chandar | Dynamic fracture of nominally brittle materials[END_REF]. Some computational modeling directly uses discrete descriptions of fracture [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Camacho | Computational Modelling of Impact Damage in Brittle Materials[END_REF][START_REF] Mastilovic | High-Velocity Expansion of a Cavity within a Brittle Material[END_REF][START_REF] Zhou | Stochastic fracture of ceramics under dynamic tensile loading[END_REF], or models based upon Smoothed Particle Hydrodynamics [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF][START_REF] Libersky | High Strain Lagrangian Hydrodynamics[END_REF] to simulate the behavior of brittle materials under shock loading [START_REF] Hiermaier | Numerical Simulation of Failure in Brittle Materials using Smooth Particle Hydrodynamics[END_REF]. On the other side of the modeling spectrum, continuum models have been used for glass (e.g., Ref. [START_REF] Johnson | An Improved Computational Constitutive Model for Brittle Materials[END_REF]) but the relationship with the microstructure is remote, and the description of damage is implicit. To account for multiple microcracking under dynamic loading conditions, the early models are based upon the analysis of cracked systems [START_REF] Budiansky | Elastic Moduli of a Cracked System[END_REF][START_REF] Margolin | Elasticity Moduli of a Cracked Body[END_REF] on a continuum level at which constitutive equations are written. The models developed by Taylor et al. [START_REF] Taylor | Microcrack-Induced Damage Accumulation in brittle Rocks under Dynamic Loading[END_REF] and Rajendran [START_REF] Rajendran | Modeling the Impact Behavior of AD85 Ceramic under Multiaxial Loading[END_REF] use the previous framework. This approach usually leads to an isotropic description of damage. Cracking is essentially anisotropic and the model needs to account for it to be realistic in terms of actual damage predictions [START_REF] Denoual | On the Identification of Damage during Impact of a Ceramic by a Hard Projectile[END_REF][START_REF] Denoual | Dynamic Fragmentation of Brittle Solids: A Multi-Scale Model[END_REF]. Combining discrete features and continuum descriptions are also possible [START_REF] Espinosa | A Finite Deformation Continuum/Discrete Model for the Description of Fragmentation and Damage in Brittle Materials[END_REF][START_REF] Denoual | Dynamic Fragmentation of Brittle Solids: A Multi-Scale Model[END_REF].

The latter approach will be followed herein.

Based upon a dynamic fragmentation analysis, an anisotropic damage model is introduced. The latter uses characteristic parameters that depend on material properties (i.e., distributions of surface and bulk initiation sites in glass [START_REF] Brajer | The role of surface and volume defects in the fracture of glass under quasi-static and dynamic loadings[END_REF]) and loading conditions. 3D-FE simulations are performed and analyzed with respect to experimental observations. Within a finite element modeling framework, a procedure is proposed to describe on the same footing micro-cracks (responsible for local damage) and mesocrack propagation. The mesh size that allows for this distinction has to be made transparent to this multiscale description. We will use the term "meso-damage" to refer to this modeling. These developments are implemented in an explicit finite element code to model edge-on-impacts. The transition between multiple fragmentation close to the impact zone and single fragmentation far from impact is studied with the meso-damage version of the model.

Experimental results

Bullets are fired by a gun and impact a glass target of surface area 100 × 100 mm 2 and of thicknesses varying between 8 and 15 mm. The projectile speed is measured by two optical cells one meter apart. An "open" configuration is used enabling for in-situ observations by utilizing a high-speed camera (Ultranac FS 501, time of exposure: 1 µs).

When the bullet reaches the second cell, flashlights are triggered. The bullet impact on the target triggers the camera to take pictures. According to the European standard EN1063 [START_REF] En | Glass in building Security glazing Testing and classification of resistance against bullet attack[END_REF], one of the threats for armors is a "soft" bullet (e.g., .44 magnum) further referred to as "BR4," traveling at a (moderate) speed of 430 m/s. BR4 bullets have a soft core made of lead and a brass envelope (Figure 1). Typical results obtained by using the high-speed camera are shown in Figure 1.

When the bullet impacts the glass sample and induces cracking a few microseconds thereafter, dark zones appear on the pictures. They correspond to damaged zones where glass is no longer transparent (caused by the light scattering by crack surfaces).

One observes the damage zone propagating in the sample and the crack front is faithfully identified by this type of inspection and consequently propagation velocities are determined with a good accuracy. After impact, fragments of the impacted sample are scattered and the BR4 core (made of lead) is fully sublimated.

A "sarcophagus" configuration that maintains fragments together on a glass target of size 100×100×10 mm 3 is also used (Figure 2). After impact, a hyper-fluid (EPOFIX) resin infiltrates (in vacuum for 12 h) the cracks and maintains the fragments together to observe the impacted sample. The projectile does not penetrate the target. There are different zones that characterize the interaction between the projectile and the target. A first zone in the vicinity of the impact is totally comminuted. A second area displays a high density of radial and hoop cracks. A third zone appears in which only few long cracks propagate. These zones correspond to different cracking features in brittle materials:

multiple fragmentation in the first zone with mainly closed cracks (glass is totally comminuted; the size of thinner fragments is estimated to be a few micrometers), this fragmentation is caused by high shear stresses just under the contact zone; multiple fragmentation under mode I loading in the second zone, due to high tensile stresses induced by the radial motion of the initial compressive wave; single fragmentation in the third zone (consisting mainly of radial macrocracks).

Fragmentation analysis

When a bullet impacts a glass target, a compressive stress wave propagates ahead of the impacted area, with the speed of longitudinal acoustic waves [START_REF] Graff | Wave Motion in Elastic Solids[END_REF]. The radial motion also induces tensile stresses that nucleate microcracks. It is assumed that the population of initiation sites leading to damage and failure in glass is identical when the material is subjected to quasi-static and dynamic loading conditions.

Poisson point process

The material is assumed to contain initiation points with random critical stresses σcr.

For each applied stress level σ, points whose critical stress is less than σ are potential initiation sites. Their density is denoted by λ in (or equivalently, the statistical distribution of initiation sites follows a Poisson point process of intensity λ in [START_REF] Gulino | Weibull Strength Statistics for Graphite Fibres Measured from the Break Progression in a Model Graphite/Glass/Epoxy Microcomposite[END_REF][START_REF] Jeulin | Modèles morphologiques de structures aléatoires et changement d'échelle[END_REF][START_REF] Ravi-Chandar | On the role of microcracks in the dynamic fracture of brittle materials[END_REF][START_REF] Denoual | A Damage Model for the Dynamic Fragmentation of Brittle Solids[END_REF]).

Consequently, the probability of finding δ points within a uniformly loaded domain Ω of (surface or volume) measure Z reads

P (δ, Ω) = (Λ in ) δ δ! exp (-Λ in ) (1)
so that Λ in corresponds to the average number of points within Ω. By definition, Λ in is related to the density λ in by Λ in = λ in Z for a uniformly loaded domain. The larger the applied stress σ, the greater the number of points for which σcr < σ. One possible choice to account for this trend is given by a power law function of the applied stress σ

λ in (σ) = λ 0 σ σ 0 m ( 2
)
where m and σ m 0 /λ 0 will be interpreted as the Weibull parameters when single fragmentation occurs, and • denotes the Macauley brackets indicating that initiation is governed by tensile stresses.

Single fragmentation

In a 3D setting, the previous analysis holds for a stress σ corresponding to maximum principal stress σ 1 . Under quasi-static loading condition, a weakest link hypothesis usually holds [START_REF]Statistical Approach to Brittle Fracture[END_REF]. Therefore, the failure probability P F is the probability of finding at least one initiation site within Ω when σ 1 ≡ σ F > 0

P F = P (δ ≥ 1, Ω) = 1 -P (δ = 0, Ω) = 1 -exp -λ 0 Z σ F σ 0 m (3)
when a uniform stress is applied. If the stress field is heterogeneous, Λ in is related to

λ in by Λ in (σ F ) = Ω λ in [σ 1 (x)]dx (4)
and the failure probability is written as

P F = 1 -exp -λ 0 Z eff σ F σ 0 m (5)
where σ F denotes the failure stress

σ F = max Ω σ 1 (x) (6)
and Z eff the effective volume or surface [START_REF] Davies | The Statistical Approach to Engineering Design in Ceramics[END_REF]

Z eff = Ω σ 1 (x) σ F m dx (7)
associated with the maximum principal stress σ 1 when positive. The mean failure stress σ F is then expressed as

σ F = σ 0 (λ 0 Z eff ) -1 m Γ 1 + 1 m (8)
and describes the fact that the larger the volume (or surface), the smaller the mean failure stress. Γ (a) denotes the (Euler) Gamma function

Γ (a) = +∞ 0 t a-1 exp(-t)dt (9)
The microstructure model accounts for random distributions of initiation sites. The Poisson-Weibull model enables one to relate the Weibull parameters to microstructural properties describing the population of initiation sites. The latter is the key for understanding probabilistic features related to the fragmentation of brittle materials.

In this first case, the initiation conditions are the only ones that control the failure properties. The Weibull parameters used for glass are either surface [START_REF] Oakley | An Empirical Study of the Effect of Stressed Area on the Strength of Float Glass Surfaces[END_REF] or bulk [START_REF] Gy | Characterization of a mode of rupture of glass at 610 • C[END_REF] parameters. Both will be considered in the sequel.

Multiple fragmentation

Under dynamic loading conditions, a propagating crack may shield other initiation sites because of the dynamic growth of a domain of measure Zo where stresses are relaxed. The cracking pattern then results from a competition between obscuration (or shielding) phenomena and the loading rate (Figure 3) initiating new cracks. The density of initiation sites λ in is divided into two parts, namely, λ b the crack density, and λo the density of obscured sites. The obscuration probability, which will be used to relate these three densities, is expressed as an extension of the classical Weibull law [34]

Po = 1 -exp -Ẑoλin(σ) with Ẑoλin(σ) ≡ T 0 Zo(T -t) dλ in dt dt = Zo λin (10)
where Ẑo is the measure of the mean obscuration zone, the convolution product, ċorresponds to the first order time derivative, and Zo the size of the obscuration zone.

The latter is dependent on the crack propagation velocity kC 0 , a shape factor S and the space dimension n

Zo(t) = S(kC 0 t) n H(t) ( 11
)
where H is the Heaviside step function. If the loading rate is small enough, Ẑo becomes equal to the entire representative volume or surface, and yields the Weibull failure probability [START_REF] Brajer | The role of surface and volume defects in the fracture of glass under quasi-static and dynamic loadings[END_REF]. For example, when n = 3, the growth law of the obscuration probability is derived from Equation ( 10)

- d 3 dt 3 ln(1 -Po) = ... Zo λ in = 6S(kC 0 ) 3 λ in (12)
Let us assume that the microscopic stress rate σ is constant. By noting that, in pure tension, the mesoscopic stress Σ is related to the microscopic (or effective) stress σ by σ = Σ/(1 -D) [START_REF] Rabotnov | Creep Problems in Structural Members[END_REF], the ultimate tensile strength (dΣ/dσ = 0), denoted by Σmax, is then expressed as

Σmax σc = 1 e Γ (m + n) n!Γ (m + 1)
1 m+n [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF] where σc is the characteristic stress associated with dynamic fragmentation

σc = σ m 0 σn λ 0 S (kC) n 1 m+n (14)
The normalized ultimate strength only depends upon the Weibull modulus m, and the space dimension n. The ultimate strength itself is then proportional to σn/(m+n) . This result is in agreement with experimental data of oil shale [START_REF] Grady | Continuum Modeling of Explosive Fracture in Oil Shale[END_REF], microconcrete [START_REF] Brara | An Experimental Method for Dynamic Tensile Testing of Concrete by Spalling[END_REF], ceramics and glass [START_REF] Hild | On the Probabilistic-Deterministic Transition Involved in a Fragmentation Process of Brittle Materials[END_REF].

It is worth noting that in the present case, not only are the Weibull parameters needed, but also a crack propagation velocity, in addition to the considered stress rate.

The multiple fragmentation regime is the result of a competition between new initiations and crack propagations, both of them described by the obscuration probability.

Transition

Equations ( 8) and ( 13) define two different regimes. The first one is obtained when a weakest link hypothesis is made. It corresponds to single fragmentation. The second one assumes multiple fragmentation. The transition between "quasi-static" and "dynamic" strength is estimated by the intersection between the weakest link and the multiple fragmentation solutions (Figure 4)

σ F = Σmax( σ). (15)
A threshold stress rate σsm allows for the description of the transition between single and multiple fragmentation regimes

σ       
< σsm single fragmentation ≥ σsm multiple fragmentation [START_REF] Johnson | An Improved Computational Constitutive Model for Brittle Materials[END_REF] with σsm = σ 0 kC 0 (λ 0 S)

1 n (Z eff λ 0 ) m+n mn e Γ (m + 1)n! Γ (m + n) Γ m+n m + 1 m 1 n (17)
This transition does not only depend on material parameters, but also involves the measure Z eff of the considered examination domain Ωe. The response of a large structure is "dynamic" for low stress rates even if the material follows a weakest link hypothesis for the same loading applied on a smaller domain.

These two regimes are depicted in Figure 4 when the surface is endowed with

Weibull parameters ms = 7, σ 0s = 100 MPa, 1/λ 0s = 100 cm 2 [START_REF] Oakley | An Empirical Study of the Effect of Stressed Area on the Strength of Float Glass Surfaces[END_REF], whereas the bulk is characterized by mv = 30, σ 0v = 3 GPa, 1/λ 0v = 10 -6 mm 3 [START_REF] Gy | Characterization of a mode of rupture of glass at 610 • C[END_REF]. For low strain rates, damage and fracture are dominated by surface initiations (i.e., the ultimate strength corresponding to surface initiations is less than that for initiations in the bulk). Conversely, for very high strain rates (i.e., greater than 10 9 s -1 ), initiations in the bulk may dominate the fragmentation regime. This result is consistent with experimental observations of glass samples [START_REF] Pugh | Glass Cracking Caused by High Explosives[END_REF][START_REF] Brajer | The role of surface and volume defects in the fracture of glass under quasi-static and dynamic loadings[END_REF].

The framework proposed herein is thus valid for describing in a unified way both fragmentation regimes. One of its advantages then lies in the fact that quasi-static experiments may be used to determine the material (Weibull) parameters. Once they are known, they can be used in a situation in which single fragmentation no longer occurs (e.g., dynamic fragmentation), provided the crack propagation velocity is known.

Microcracking analysis

The previous fragmentation results are now incorporated in an anisotropic damage model [START_REF] Denoual | A Damage Model for the Dynamic Fragmentation of Brittle Solids[END_REF]. The level of stress and stress rate determine the type of fragmentation, namely single or multiple fragmentation. A so-called multi-scale model [START_REF] Denoual | Dynamic Fragmentation of Brittle Solids: A Multi-Scale Model[END_REF] is used to determine the transition between single and multiple fragmentation regimes. This fragmentation model accounts for quasi-static and dynamic loading conditions in each examination volume or surface. The formulation chosen here would strictly require that the principal stress directions remain fixed. In practice it is believed that this requirement is not much limiting.

Anisotropic damage model

Up to now, only the maximum principal stress direction was considered. In the following, anisotropy has to be implemented to account for the directionality of cracking as observed experimentally. By averaging over an examination domain, the mesoscopic

principal stress Σ i (i = 1, n) is equal to σ i (1 -D i )
, where σ i is the local (i.e., microscopic and unobscured) principal stress. The damage variable D i is equal to Po associated with a given direction of the eigen stress i (this definition follows the usual assumption of D i = 0 for an undamaged domain and D i = 1 for a fully damaged, or obscured, domain [START_REF] Lemaitre | A Course on Damage Mechanics[END_REF]), and the growth law of each damage variable D i for bulk initiations is derived from Equation ( 12)

d 2 dt 2 1 1 -D i dD i dt = 6S(kC 0 ) 3 λin (σ i ) (18)
when σ i > 0 and dσ i /dt > 0, and with the modified density associated with the Weibull model expressed as [START_REF] Denoual | Dynamic Fragmentation of Brittle Solids: A Multi-Scale Model[END_REF]

V FE λin (σ i ) =        0 if σ i < σ F max V FE λ 0v σi σ0v mv , 1 otherwise (19)
where V FE is the volume of the considered finite element (i.e., coincident with the examination domain Ωe of the fragmentation analysis of Section 3.3). The value σ F is a stress obtained by random selection of a real number R lying between 0 and 1, and inverting P F (σ F ) = R when V eff ≡ V FE [see Equation [START_REF] Brajer | The role of surface and volume defects in the fracture of glass under quasi-static and dynamic loadings[END_REF]]. For surface initiations, only two "in-plane" damage parameters are needed (i.e., one principal stress is always equal to 0) and their growth law read

d dt 1 1 -D i dD i dt = 2S(kC 0 ) 2 λin (σ i ) (20)
when σ i > 0 and dσ i /dt > 0, and

S FE λin (σ i ) =        0 if σ i < σ F max S FE λ 0s σi σ0s ms , 1 otherwise (21)
where S FE is the external surface of the considered finite element. Equations ( 19) and ( 21) express the average number of initiation sites within a volume or on the surface that is equal to 0, 1 or greater than 1. The integration of Equations ( 18) and ( 20)

usually assumes zero initial conditions.

The compliance tensor S is written in the principal frame as [22]

S = 1 E                      1 1-D1 -ν -ν 0 0 0 -ν 1 1-D2 -ν 0 0 0 -ν -ν 1 1-D3 0 0 0 0 0 0 1+ν (1-D2) α (1-D3) α 0 0 0 0 0 0 1+ν (1-D3) α (1-D1) α 0 0 0 0 0 0 1+ν (1-D1) α (1-D2) α                      (22)
where α = S ta /Sno ≈ 1/2, with S ta and Sno are shape factors that describe the obscuration zone of tangential and normal stresses, E Young's modulus of the undamaged material and ν the corresponding Poisson's ratio.

Simulation of EOI experiments

The edge-on-impact experiment is now simulated by using an explicit finite element code. Materials of the bullets are assumed to behave elasto-plastically. Glass is assumed to be elasto-plastic (E = 70 GPa, ν = 0.22) with a Hugoniot elastic limit (HEL) evaluated at 6 GPa. Fragmentation will be described by the previous damage model. The FE mesh accounts for surface properties (when needed) endowed with Weibull parameters ms = 7, σ 0s = 100 MPa, 1/λ 0s = 100 cm 2 , whereas the bulk is characterized by mv = 30, σ 0v = 3 GPa, 1/λ 0v = 10 -6 mm 3 .

Damage is mainly nucleated on external surfaces of the samples since the maximum strain rates are significantly less than the value 10 9 s -1 , which corresponds to the transition from multiple fragmentation induced by surface and bulk initiations (Figure 4). Figure 5 shows that multiple fragmentation appears in the first half of the target and along the side. This 3D computation is expected to describe accurately the initial stage of bullet impact. It does not account for the propagation of long cracks but only for crack initiation as illustrated in Figure 5. This limitation will be overcome by the introduction of the following "meso-damage" approach where crack extension will be allowed to extend beyond single finite elements.

A "meso-damage" approach

The present section aims at extending the previous damage model to account for the propagation of long cracks.

Extension of the damage model

New damage parameters ∆ i are introduced in order to describe cracks that are longer than the elements so that the mesoscopic principal stress Σ i is now equal to σ i (1 -

D i )(1 -∆ i ).
Fracture propagation is described by the obscuration zone that grows at a constant speed kC 0 = 1500 m/s [START_REF] Brajer | The role of surface and volume defects in the fracture of glass under quasi-static and dynamic loadings[END_REF] from element to element. Again, the damage variables ∆ i are equal to the ratio between the obscuration volume and the total volume in one element (Figure 6) associated with the three cracking directions. The growth of ∆ i is obtained by direct derivation of S(kC 0 t) 3 /V FE (when n = 3)

d 2 dt 2 d∆ i dt = 6S(kC 0 ) 3 λin (σ i ) ( 23
)
where λin is equal to 1/V FE , which means that a single mesocrack propagates in the considered element. One may note that this growth law is identical to that used in the case of multiple fragmentation ( 18) except for the factor 1/(1 -D i ) that accounts for overlappings of obscuration zones. All elements where ∆ i = 0 are listed and they numerically correspond to a crack whose length is that of all elements traversed by the latter. In practice, the local damage variables start to grow, i.e., D i > 0 (Figure 6a).

When D i reaches a critical value Dc corresponding to coalescence in the considered element, ∆ i starts to grow towards the neighboring element (Figure 6b). The size of the crack is followed by computing the size of the obscuration zone. Two cases may occur:

the stress intensity factor is high enough to let it propagate (i.e., σ i (1

-D i ) √ a >
Kc, where Kc denotes the toughness of the material) and the element begins to crack and the kinetics of ∆ i is given by Equation ( 23);

the stresses are not high enough; no propagation occurs in the neighborhood of the initially damaged element. For the sake of simplicity, the propagation and arrest conditions are associated with the same value of Kc.

This process is repeated as the stress wave invades the material. The macrocrack propagation may therefore be the result of a series of "microscopic" coalescences when the damage variables D i are high enough (D i > Dc). These multiple events induce a velocity of the crack front that can be greater than kC 0 and a multiple fragmentation regime is observed. Conversely, the propagation conditions may be satisfied (i.e.,

σ i (1 -D i ) √ a > Kc)
even though the local values of D i are small. A single fragmentation regime is then observed.

Effect of element size

An in-house FE code is used to simulate the propagation of a cylindrical divergent wave [START_REF] Denoual | On the Identification of Damage during Impact of a Ceramic by a Hard Projectile[END_REF]. Two element sizes are used when a cavity in a slab is subjected to a pressure step of -6 GPa (Figure 7). The space-time diagram of the damage variable ∆ 1 is

shown when using the meso-damage model. The white zone corresponds to values of ∆ 1 greater than Dc = 0.85. Two different mesh densities are used and the spacetime diagrams are virtually identical. Except for the scatter induced by the multi-scale damage model, the results are quasi mesh-independent. The damage variable grows up to a given point where it stops when the propagation condition is no longer satisfied.

Divergent cylindrical wave

The previous 3D simulations (Section 4.2) help us in determining a contact force between the .44 Magnum bullet and the surface of the target (Figure 8). During the first microseconds after impact, the initial shock wave propagates followed by the collapse of the projectile so that the contact force increases. For the sake of simplicity, the propagation of a cylindrical divergent wave is assumed.

Figure 9 shows the results obtained by using the damage variable ∆ 1 . In the first part of the interaction, the stress level induced by the compressive wave are high enough to initiate several cracks in each finite element and the damage front velocity is related to the speed at which the induced tensile (hoop) stress develops. The zigzagged profile corresponds to multiple initiations and propagations. The apparent propagation velocity of the damage front is greater than kC 0 since this phenomenon is dominated by initiation events. Later, as observed (Figure 1), long cracks begin to propagate at a constant speed kC 0 = 1500 m/s with no further initiation. In that regime, it is described by propagation conditions. This transition could not have been assessed by using a microscopic damage model alone. Last, the dashed line in Figure 9 shows the result for the experimental location of the crack front evaluated from Figure 1. The two regimes are found and the transition is well captured.

Edge-on impact

Figure 10 compares the simulations of the edge-on-impact test obtained with the microand mesoscopic versions of the damage model. Multiple fragmentation in the first half of the target and along the side is followed by single fragmentation as already shown in Figure 9. This 3D computation describes more accurately the initial stage of bullet impact and the propagation of long cracks (Figure 1) than with the local damage model.

The last zone to be better accounted for is that close to the impact point in which damage develops under highly confined conditions. This damage regime is only indirectly included in the elastoplastic description associated with the chosen HEL value, and is therefore not depicted in Figure 10 since only mode I cracks are explicitly described in the present model. However, it was shown that as a first approximation this type model is able to account for confined damage in rocks [START_REF] Grange | On the Dynamic Fragmentation of two Limestones using Edge-on Impact Tests[END_REF].

Conclusion

Edge-on-impact tests have been performed on sodalime silicate glass by using a soft bullet. These experiments allow for a visualization of the crack pattern by using a highspeed camera or a sarcophagus configuration. Fragmentation of impacted sodalime silicate glass is mainly caused by the growth of cracks nucleated on the surface of glass.

Far from the impact surface, propagation of long cracks is always observed.

Numerically, the zone of multiple fragmentation is reproduced with the help of an anisotropic damage model based on a probabilistic analysis of multiple crack activations. Since the population of initiation sites is the key parameter that governs the cracking pattern, Weibull parameters representative of the surface and the bulk properties have been used. However, this model is not sufficient to account for the propagation of long cracks. An extension was therefore proposed.

As the propagation of long cracks is a non-local phenomenon for classical FE simulations, the zone obscured by one macrocrack is used to follow its propagation. This is achieved through a "meso-damage" model coupled with the previous one. 3D simulations prove that the different zones observed in edge-on-impact tests where mode I cracks initiate and propagate are reproduced by the present approach. A local approach would only give the initiation sites, and would not be able to capture the mesocrack formation. The different regimes of damage front velocities observed experimentally are then well captured in the simulations with the "meso-damage" model.

Even though the model proposed herein was used to simulate the dynamic fragmentation of glass, its domain of application may not be restricted to glass. The microscopic damage model was successfully used to analyze the fragmentation of SiC ceramics [START_REF] Denoual | A Damage Model for the Dynamic Fragmentation of Brittle Solids[END_REF][START_REF] Denoual | Dynamic Fragmentation of Brittle Solids: A Multi-Scale Model[END_REF][START_REF] Hild | On the Probabilistic-Deterministic Transition Involved in a Fragmentation Process of Brittle Materials[END_REF][START_REF] Forquin | Effect of Aluminum Reinforcement on the Dynamic Fragmentation of SiC Ceramics[END_REF], high performance concrete [START_REF] Forquin | Dynamic Fragmentation of an Ultra-High Strength Concrete during Edge-On Impact Tests[END_REF], and limestones [START_REF] Grange | On the Dynamic Fragmentation of two Limestones using Edge-on Impact Tests[END_REF]. In the two last cases, there were situations for which the crack length was greater than the element size, thereby calling for the same type of extensions as those proposed herein.

Acknowledgements

This work has been supported by a Saint-Gobain grant.

List of Figures

Fig. 1 Fig. 2

 12 Fig. 1 Sequence of pictures of the impact of a BR4 bullet (velocity = 430 m/s) on a 100 × 100 × 8 mm 3 glass target. The dark areas correspond to damaged zones. The dashed line shows the main crack used in Figure 9 for comparison between predicted and experimental results.

Fig. 3

 3 Fig.3Examination domains in tension. The dots depict points that are potential initiation sites. Left) The tensile stress rate is low; one single crack has initiated, leading to the entire obscuration of the element. Right) The tensile stress rate is high; there is a competition between the stress increase and the obscuration of the element. Several cracks are initiated.

Fig. 4 Fig. 5 FE

 45 Fig.4Average ultimate tensile strength vs. strain rate for a 10 mm 3 element made of sodalime silicate glass (the surface is endowed with Weibull parameters ms = 7, σ 0s = 100 MPa, S eff = 100 cm 2[START_REF] Oakley | An Empirical Study of the Effect of Stressed Area on the Strength of Float Glass Surfaces[END_REF], whereas the bulk is characterized by mv = 30, σ 0v = 3 GPa, V eff = 10 -6 mm 3[START_REF] Gy | Characterization of a mode of rupture of glass at 610 • C[END_REF]).

Fig. 6 Fig. 7 Fig. 8 Fig. 9 Fig. 10 FE

 678910 Fig. 6 Schematic view of a propagation of a long crack via the meso-damage model. The first step (left) corresponds to distributed damage modeled by damage variables D i . The second step (right) corresponds to coalescence in one element described by meso-damage variables ∆ i .

  Comparison of the space-time diagrams of the damage variable ∆ 1 for two different meshes: 4 elements per mm (a) and 16 elements per mm (b). The white zone corresponds to ∆ 1 greater than 0.85. . . . . . . . . 31 8 Contact force vs. time for an impact of a BR4 bullet on a sodalime silicate target determined by a 3D FE simulation. . . . . . . . . . . . . . 32 9 Space-time diagram of the damage variable ∆ 1 obtained with the mesodamage model. The white zone corresponds to ∆ 1 greater than 0.85. It is compared with an experimental crack front velocity evaluated from Figure 1 (dotted line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 10 FE simulations of one quarter of an EOI configuration. Top: the damage variable ∆ 1 associated with the first (maximum) principal stress shows the location of cracks. Bottom: the damage variable D 1 associated with the first principal stress shows the location of crack initiations. . . . . . 34

	7
	1	Sequence of pictures of the impact of a BR4 bullet (velocity = 430 m/s)
		on a 100 × 100 × 8 mm 3 glass target. The dark areas correspond to
		damaged zones. The dashed line shows the main crack used in Figure 9
		for comparison between predicted and experimental results. . . . . . . . 25
	2	Cracking pattern of impacted glass (44 Magnum bullet, V = 330 m/s)
		obtained with a sarcophagus configuration. Near the impact zone, the
		material is fully comminuted. Far from impact, long cracks are observed. 26
	3	Examination domains in tension. The dots depict points that are poten-
		tial initiation sites. Left) The tensile stress rate is low; one single crack
		has initiated, leading to the entire obscuration of the element. Right)
		. 29
	6	Schematic view of a propagation of a long crack via the meso-damage
		model. The first step (left) corresponds to distributed damage modeled
		by damage variables D i . The second step (right) corresponds to coales-
		cence in one element described by meso-damage variables ∆ i . . . . . . . 30

The tensile stress rate is high; there is a competition between the stress increase and the obscuration of the element. Several cracks are initiated.

27 4

Average ultimate tensile strength vs. strain rate for a 10 mm 3 element made of sodalime silicate glass (the surface is endowed with Weibull parameters ms = 7, σ 0s = 100 MPa, S eff = 100 cm 2 [32], whereas the bulk is characterized by mv = 30, σ 0v = 3 GPa, V eff = 10 -6 mm 3 [33]). 28 5 FE simulations of one quarter of an EOI configuration. The damage variable D 1 associated with the first (maximum) principal stress shows the location of crack initiations. . . . . . . . . . . . . . . . . . . . . . .

"! #!


Julien Réthoré
email: julien.rethore@insa-lyon.fr

Stéphane Roux

Franc ¸ois Hild

Hybrid analytical and extended finite element method (HAX-FEM): A new enrichment procedure for cracked solids

Keywords: X-FEM, Williams' series, energy coupling, stress intensity factors

A new numerical method is proposed to estimate directly stress intensity factors at the tip of a crack in an elastic body. In the same spirit as for the extended finite element method, the approximation of the displacement field is enriched in the vicinity of the crack tip. Yet the method proposed herein differs by the way the enrichment is introduced. Instead of using partition of unity concepts, a two-description formulation is implemented. The first one uses standard finite elements while the second one, in the vicinity of the crack tip, resorts to a few purely analytical expressions. These two descriptions are then coupled by partitioning the energy with an overlapping zone. The performance of this new enrichment technique is illustrated and compared with existing techniques by means of two examples.

INTRODUCTION

Since the pioneering papers of Black and Belytschko [START_REF] Black | Elastic crack growth in finite elements with minimal remeshing[END_REF] and Moës et al. [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], the extended finite element method (X-FEM), initially dedicated to elastic solids, was extensively developed to deal with many applications. Concerning the simulation of failure in solids, 3D crack propagation [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets. Part I: mechanical model[END_REF][START_REF] Gravouil | Non-planar 3D crack growth by the extended finite element and level sets. Part II: level set update[END_REF], dynamic crack propagation [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF][START_REF] Duarte | A generalized finite element method for the simulation of three-dimensional dynamic crack propagation[END_REF][START_REF] Réthoré | An energy conserving scheme for dynamic crack growth with the extended finite element method[END_REF][START_REF] De Borst | Mesh-independent numerical representations of cohesive-zone models[END_REF][START_REF] Grégoire | Dynamic crack propagation under mixed-mode loadingcomparison between experiments and X-FEM simulations[END_REF][START_REF] Menouillard | Efficient explicit time stepping for the extended finite element method[END_REF] were simulated, non-linear models [START_REF] Elguedj | Appropriate extended functions for X-FEM simulation of plastic fracture mechanics[END_REF][START_REF] De Borst | Discrete vs smeared crack models for concrete fracture: bridging the gap[END_REF], and multi-grid solvers [START_REF] Rannou | A local multigrid X-FEM strategy for 3-d crack propagation[END_REF] were implemented, and coupling with molecular models [START_REF] Gracie | On a new extended finite element method for dislocations: core enrichment and nonlinear formulation[END_REF][START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF] was shown to be possible. Even full field measurement techniques such as digital image correlation [START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF][START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Réthoré | Experimental investigation of localized phenomena using digital image correlation[END_REF][START_REF] Réthoré | Extended three-dimensional digital image correlation (X3D-DIC)[END_REF] now use enriched kinematics.

X-FEM rests on different enrichments, namely, displacement discontinuities across the crack faces, and (standard crack) singular functions at the crack tip. Ironically, the latter enrichments cannot be trusted to provide reliable estimates of stress intensity factors (SIFs) [START_REF] Nicaise | Optimal convergence analysis for the extended finite element method[END_REF]. Hence, most of these analyses require a post-processing stage to extract SIFs from the simulated displacement field. This stage is usually performed by using interaction integrals [START_REF] Black | Elastic crack growth in finite elements with minimal remeshing[END_REF] or least-squares techniques [START_REF] Duarte | A generalized finite element method for the simulation of three-dimensional dynamic crack propagation[END_REF].

In References [START_REF] Liu | XFEM for direct evaluation of mixed mode sifs in homogeneous and bi-materials[END_REF][START_REF] Xiao | Direct evaluation of accurate coefficients of the linear elastic crack tip asymptotic field[END_REF], an alternative enrichment basis is proposed where SIFs K I and K II are directly evaluated as parameters of the discretization basis. It is referred to as KX-FEM. It is shown in these references that nodal enriched degrees freedom corresponding to K I and K II must be constrained to have the same value. This strategy is generalized in Reference [START_REF] Chahine | Crack-tip enrichment in the XFEM method using a cut-off function[END_REF] by the use of a cut-off function. Other requirements to get accurate SIFs are to include in the basis not only the usual asymptotic terms (with K I and K II as pre-factors) but also higher-order terms [START_REF] Liu | XFEM for direct evaluation of mixed mode sifs in homogeneous and bi-materials[END_REF].

Having implemented enriched kinematics in a digital image correlation framework [START_REF] Réthoré | Experimental investigation of localized phenomena using digital image correlation[END_REF] spurred the present authors to work on the improvement of the robustness of SIF extraction techniques. For instance, noise robust interaction integrals [START_REF] Réthoré | Noise-robust stress intensity factor determination from kinematic field measurements[END_REF] but also a general formalism of the extraction problem (Réthoré et al., submitted 2008) were recently proposed. In the latter reference, one of the conclusions was that the best determination of SIFs is obtained when they are directly the unknowns of the correlation problem as in integrated approaches of digital image correlation [START_REF] Roux | Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches[END_REF]. One of the driving idea of the present paper is to propose a numerical methodology that allows for direct extractions of SIFs with no additional data processing.

Although the method proposed by Liu et al. [START_REF] Liu | XFEM for direct evaluation of mixed mode sifs in homogeneous and bi-materials[END_REF] and Xiao and Karihaloo [START_REF] Xiao | Direct evaluation of accurate coefficients of the linear elastic crack tip asymptotic field[END_REF] is very appealing, it will be shown herein that the underlying partition of unity that is kept active over the enriched region reduces drastically the potential of the method. It can be traced back to the conditioning of the system [START_REF] Nicaise | Optimal convergence analysis for the extended finite element method[END_REF][START_REF] Chahine | Crack-tip enrichment in the XFEM method using a cut-off function[END_REF][START_REF] Chahine | Spider-XFEM, an extended finite element variant for partially unknown crack-tip displacement[END_REF]. By reducing the number of degrees of freedom of the enrichment leads to a better conditioning, and hence a more accurate evaluation of SIFs. An alternative to the partition of unity enrichment is proposed herein. An analytical model written in Williams' series formalism [START_REF] Williams | On the stress distribution at the base of a stationary crack[END_REF] is used in the vicinity of the crack tip where the finite element model is not activated. The analytical model includes rigid body translations, rotation, classical asymptotic terms (and T -stress and higher-order term can also be included). Then both descriptions are matched using an energy coupling procedure or Arlequin method [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF].

Section 2 briefly recalls the family of functions describing the displacement of a cracked elastic solid. In Section 3, the new formulation is proposed and its discretization is discussed in Section 4. In Section 5, the performance of the proposed enrichment technique is illustrated with two examples.

CRACK DISPLACEMENT FIELD

Let us consider a homogeneous body with isotropic elastic behavior, and a 2D setting, the displacement field u is conventionally represented by its complex writing, u = u x +i u y . It was expanded by Williams [27] for a straight crack as a double series

u(r, ) = i=I,II n c n i / n i (r, ) (1)
where c n i are the influence parameters, and / n i reference fields that form the complete set of linear elastic fields satisfying a zero traction condition along the crack path

• even n / n I (r, ) = r n/2 e in /2 - n 2 e i(4-n) /2 + n 2 +1 e -in /2
(2)

/ n II (r, ) = ir n/2 e in /2 + n 2 e i(4-n) /2 - n 2 -1 e -in /2
(3)

• odd n

/ n I (r, ) = (-1) (n+1)/2 r n/2 e in /2 - n 2 e i(4-n) /2 + n 2 -1 e -in /2 (4) / n II (r, ) = i(-1) (n-1)/2 r n/2 e in /2 + n 2 e i(4-n) /2 - n 2 +1 e -in /2 (5)
where is the Kolossov constant, namely, = (3-)/(1+ ) for plane stress or = (3-4 ) for plane strain conditions, being Poisson's ratio. The partition between modes I and II is performed on the basis of the following symmetries

/ n I (z) = / n I (z) / n II (z) = -/ n II (z) (6)
All these fields are homogeneous functions of the distance to the crack tip r . With this convention, all elastic displacement fields with finite strain energy densities are generated with 0 n<∞. All odd indices label fields with a discontinuity across the crack mouth, whereas even indices correspond to continuous fields. To mention just a few examples, n = 0 corresponds to translations either parallel (hence termed I) or perpendicular (labeled II) to the crack path, n = 1 are the classical mode I and mode II displacement fields whose amplitudes are proportional to the corresponding SIFs.

The case n = 2 gives access to either the T -stress component for mode I, or rigid body rotation about the crack tip for mode II (no discontinuity is involved here since 2 is even). Larger indices correspond to 'subsingular' or higher-order fields that may capture the remote heterogeneity of the loading, but do not affect the fracture behavior at the crack tip. This family of fields is thus the appropriate basis function to describe the displacement field for a traction-free crack in an elastic solid.

PROBLEM FORMULATION

Owing to the difficulties described in the introduction, one wants to enrich a finite element discretization using the family of functions introduced above, but not in the standard extended finite element formalism. The underlying partition of unity must be deactivated over the region where the enriched part of the displacement basis is adopted. For this purpose, the energy coupling (or Arlequin) procedure proposed by Ben Dhia and Rateau [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF] is chosen.

For each description of the displacement field, the derived elastic stress must satisfy the balance of momentum (body forces are neglected herein)

divr = 0 ( 7
)
where r is Cauchy's stress tensor, and div the divergence operator. The domain is now subdivided into two overlapping sub-domains 1 and 2 . In each of these sub-domains, a description u 1 , respectively, u 2 , of the displacement field is adopted. The boundary conditions are held by the sub-domain 1 , namely, tractions t d and displacements u d are prescribed on two complementary parts * 1t and * 1u of its boundary.

From the strong form of the equilibrium equation ( 7), two weak forms (one for each description) are derived, and are weighted by a partition of unity over the overlapping region 12

= 1 ∩ 2 1 (x)+ 2 (x) = 1 ( 8
)
where 1 and 2 are the respective weights given to each energy contribution. The definition of these weighting functions is extended over the entire domain so that the description 1 (resp. 2) holds in \ 2 only (resp. \ 1 ). In this setting, the linear l i and bilinear a i weighted energy forms are obtained for i = 1 or 2

a i (u i , v * i ) = i i e(u i ) : C: e(v * i ) d (9) l 1 (v * 1 ) = * 1t 1 t d •v * 1 dS ( 10
)
where C is Hooke's tensor, e the infinitesimal strain tensor, and v * i admissible displacement fields. Note that l 2 is 0 since the crack faces are traction free.

In addition to this partition of energy, a kinematic coupling must be added to ensure the compatibility of the two displacement fields in the sense of a scalar product between the displacement gap and Lagrange multipliers k. If is the L 2 scalar product over the coupling region 12 , we have

(u 1 -u 2 , l * ) = 12 l * •(u 1 -u 2 ) d = 0 ( 11
)
where l * are admissible Lagrange multipliers. The total formulation thus reads: find

(u 1 , u 2 , k) such that for all admissible (v * 1 , v * 2 , l * ) fields, a 1 (u 1 , v * 1 )+a 2 (u 2 , v * 2 )+ (u 1 -u 2 , l * )+ (v * 1 -v * 2 , k) =l 1 (v * 1 )+l 2 (v * 2 ) (12)
Using this weak formulation, a discretization scheme is being adopted for the two descriptions.

DISCRETIZATION AND NUMERICAL ELABORATION

After the weak formulation is chosen, the discretization schemes for the two descriptions and some details on the numerical elaboration are detailed in the sequel.

Outer domain (Model 1)

For the outer domain, model 1, a finite element discretization is adopted. is divided into a given number of elements, a subset being 'active' as they define the discretized domain h 1 and the complementary subset being 'inactive'. Further, a discontinuous enrichment is added to the interpolation to account for the displacement jump across the part of the crack faces that lies in 1 . The interpolation of the displacement for the outer model thus reads

u 1 (x) = i∈N 1 N i (x)d i + i∈N cut N i (x)H(x)b i = D T U 1 ( 13
)
where N i are standard finite element shape functions supported by the set of nodes N 1 included in h 1 , d i the corresponding degrees of freedom. In this equation, H(x) is the discontinuous symmetrized Heaviside step function that enables with its associated degrees of freedom b i for the description of displacement discontinuities for the nodes in the subset N cut of N 1 that have their support cut by the crack. This interpolation relies on a usual extended finite element discretization based on the partition of unity property of N i [START_REF] Babuska | The partition of unity method[END_REF]. Matrix D collects the value of the discretization functions and vector U 1 the degrees of freedom of the outer domain.

In the following, only regular meshes of quadrangular elements will be used. For the numerical integration of the discontinuous function in the weak form, the quadrangles are subdivided into n sub ×n sub quadrangular integration cells containing one Gauss point. Standard Gauss quadrature rule with four integration points per quadrangle is used for the standard terms.

Inner domain (Model 2)

For the inner domain, model 2, the functions are selected in the basis presented in Section 2. The displacement field u 2 (x) reads

u 2 (x) = 1 2 √ 2 n∈[0;n max ] (/ n I (x) p n +/ n II (x)q n ) = U T U 2 (14)
where is Lamé's coefficient, p n and q n the degrees of freedom associated with mode I and II functions, and n max an integer that defines the maximum order considered in the interpolation. Matrix U collects the values of the discretization functions and vector U 2 the degrees of freedom of the inner domain. Note that for n = 1, p 1 and q 1 enable for a direct evaluation of K I and K II .

The discretized domain h 2 is defined as the area of the 'inactive' finite elements plus those in the coupling zone h 12 . The numerical integration of the weak form a 2 is performed in the same way for the enriched terms of a 1 , namely, by using the underlying finite elements that are subdivided into n sub ×n sub quadrangular integration cells containing one Gauss point. The discretization functions adopted for the inner domain satisfy Equation [START_REF] Réthoré | An energy conserving scheme for dynamic crack growth with the extended finite element method[END_REF] and thus the integration of a 2 could be split into two domain integrals (the first one defined on the domain for which 2 = 1 vanishes, and the second one in h 12 ), and a contour integral of the normal tractions on the inner boundary of h 12 . Up to the influence of the weighting function, this contour integral is exactly the one that is required for the hybrid crack element (see e.g. [START_REF] Xiao | Implementation of hybrid crack element on a general finite element mesh and in combination with XFEM[END_REF]). In practice, it is more convenient to integrate the initial form of a 2 .

Coupling

The discretized coupling domain is the subset of elements in h 12 . Lagrange multipliers are discretized using standard finite element shape functions

k(x) = i∈N 12 N i (x)l i = L T K ( 15
)
where N 12 is the set of nodes in h 12 , N i standard finite element shape functions, and l i their associated degrees of freedom. Matrix L collects the values of the discretization functions and vector K the degrees of freedom of Lagrange multipliers. Since the coupling terms involve discontinuous functions and basis functions of description 2, the above-mentioned quadrature using cells with one Gauss point is used.

Numerical elaboration

After the discretization schemes have been introduced in the weak form, a linear system is obtained

⎡ ⎢ ⎢ ⎣ K 1 0 C T 1 0 K 2 C T 2 C 1 C 2 0 ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ U 1 U 2 K ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ F 1 F 2 0 ⎤ ⎥ ⎥ ⎦ (16)
with

K 1i j = h 1 1 e(D i ) : C: e(D j ) d ( 17
)
K 2i j = h 2 2 e(U i ) : C: e(U j ) d (18)
C 1i j = h 12 L i •D j d (19)
C 2i j = - h 12 L i •U j d ( 20
)
are elementary blocks of the global matrix and

F 1i = * t 1 t d •D i dS (21)
F 2i = * t 2 t d •U i dS (22)
the external force vectors. Note that, in practice, F 2 is equal to zero. For practical reasons, the weighting functions i are computed numerically by solving a Laplacian problem ( i = 0) using a finite element method in the coupling zone h 12 . For example, 1 is computed using the appropriate boundary conditions on the boundary of h 12 and 2 is then directly obtained using Equation (8).

EXAMPLES

In the following, two examples are discussed. The first one corresponds to a closed-form solution of a crack in an infinite domain, i.e. the extension of the inner field to infinity, and it will test the compatibility of the entire global scheme. The second one is more informative as it corresponds to a problem that has no exact (closed-form) solution.

Infinite plate with a semi-infinite crack

The first example is dedicated to a convergence study in the idealized case of an infinite plate with a semi-infinite crack. The numerical model consists of a 39×39 quadrangular elements with a horizontal centered crack whose tip is located in the middle of the plate. To fall within the hypothesis of an infinite plate, asymptotic fields for mode I, with K Io as a prescribed SIF, and mode II, with K IIo as a prescribed SIF, are imposed on the outer boundary of the mesh. The plate is considered to be of area 1 m 2 , the material properties are set to 210 MPa for Young's modulus and 0.3 for Poisson's ratio. The prescribed SIF values are 1 MPa √ m for K Io and K IIo . 16×16 quadrature cells in each quadrangles are used where needed. Let us focus on the influence of the size of the analytical model (that contains indices ranging between 0 and 2, i.e. n max = 2) and also on the size of the coupling zone.

The size of h 2 is defined as R outer = hr outer (h being the element size), and that of 'inactive' elements R inner =r inner h. The size of the coupling region is h overlap = R outer -R inner . As an example, Figure 1 shows the weighting function 1 of the finite element model for an analytical model of 10-node layer around the crack tip, r inner = 10 and an overlapping zone of 5-node layer, overlap = 5. Figure 2 shows the two components of the displacement field for the two models when r inner = 3 and overlap = 3.

In order to study the influence of the size of the inner domain and its coupling with the FE zone (outer domain), r inner is varied in the interval 1 r inner 10 and 1 overlap 5. Figures 3 and4 summarize the results. In Figure 3, the relative errors on K I and K II are plotted as functions of r inner for different coupling zone sizes overlap . From the plots using a linear scale, a good convergence is observed. It is to be noted that estimated K I s are first less than the prescribed value, then higher and finally converge to K Io itself. For mode II, only higher values are observed. The log-scale plots show the same trends for mode I and mode II, namely, a linear decrease of the absolute relative error is obtained, the size of the coupling zone acting as a multiplicative pre-factor. For large domain sizes r inner = 10, (i.e. 1 4 of the size of the mesh), absolute errors on K I as low as 10 -4 and 2×10 -5 on K II are obtained for overlap = 5 coupling zones. The effect of the size of the coupling zone is shown in Figure 4. The relative error on SIFs follows the same type of dependence with overlap . Even though increasing the size of the analytical model and using larger coupling zone reduces the error on estimated SIFs, the values obtained with as small parameter values r inner = overlap = 1, yields 2.5×10 -3 errors for K I and 5.0×10 -3 for K II , which are very low values.

Cracked plate in remote tension

In the previous example, because of the idealized boundary conditions, only classical asymptotic terms were activated. This is not the case when realistic boundary conditions are used. The present example aims at evaluating this influence. A cracked plate whose geometry is parameterized by its width w, its length L and the length of the crack a is considered. Under remote tension, the

In the numerical model w is set to 7×10 -3 m, a is equal to 0.5 w, and L is equal to 17×10 -3 m. Young's modulus is equal to 200 MPa, and Poisson's ratio to 0.3. The applied stress is 10 MPa, which gives 2.98 MPa √ m for K Io the reference value for K I . The finite element mesh consists of 49×119 quadrangular elements that are subdivided into 8×8 integration cells when needed. Figure 5 shows the displacement field for r inner = overlap = 1.

Owing the use of realistic boundary conditions, the normalized values of K I obtained for n max = 2 shown in Figure 6 lead to an error increase when the domain of influence of the analytical model is enlarged (by increasing its size or the size of the coupling zone). The traction-free edge in front of the crack tip activates higher-order term in Williams' series. However, using the smallest domain of influence (i.e. r inner = overlap = 1), the normalized value of K I is 1.0146 (i.e. less than 1.5% error, which is very low).

In Figure 7, r inner = 10, overlap = 3 and the effect of the number of terms n max is investigated. A stabilization of the normalized SIF is observed for a value to 1.002, when seven terms are considered. For different values of n max , the normalized value of K I is plotted as a function of r inner in Figure 8(a). One notes that the strong dependence observed for n max = 2 is decreasing when more terms are considered. For n max = 7, the normalized SIF K I remains almost constant, varying between 0.999 and 1.0012. Let us note that for the smallest domain of influence, r inner = overlap = 1, the normalized value of K I is 0.999 (i.e. less than 0.1% error, which is again very low). Figure 8(b) shows the results obtained using KX-FEM under the same conditions. The cut-off function is exactly the weighting function 2 and the same number n max of higher-order terms are considered (only odd terms are introduced in KX-FEM because the underlying partition of unity is assumed to capture the even terms that are not discontinuous across the crack faces). The numerical integration is also performed using the same quadrature. For KX-FEM, an optimal value exists for r inner ≈ 4-5 for a normalized value of K 1 of 0.93 for n max = 1 and 0.983 for n max = 5. Increasing n max , the results are improved but a critical dependence on r inner is still observed, namely, for too small r inner , the domain of influence of the enrichment function is not large enough for a robust SIF estimation but for too large r inner the conditioning of the system is increasing [START_REF] Nicaise | Optimal convergence analysis for the extended finite element method[END_REF], which also reduces the robustness of the SIF estimation.

For different values of n max , Figure 9(a) shows the effect of the mesh size h while the ratio of the number of elements along the length and width of the sample is kept constant and r inner = overlap = 1. For all values of n max , a standard convergence is obtained with a rate between 1.5 and 1.25 that decreases when n max increases. Note, that the convergence rate obtained with the present enrichment procedure is faster than that of KX-FEM (Figure 9(b)). Even in its optimal version [START_REF] Nicaise | Optimal convergence analysis for the extended finite element method[END_REF], the latter gives a convergence rate of about 0.5. Here, for comparison purposes, the cut-off function is equal to 2 and thus the outer radius of the domain of influence of the enrichment decreases when the element size decreases (r inner = overlap = 1). Further 2 is numerically computed using first-order finite elements. It is thus continuous but with discontinuous derivatives, this is not the optimal choice for KX-FEM (Nicaise et al. [START_REF] Nicaise | Optimal convergence analysis for the extended finite element method[END_REF] use an analytically defined fifth-order polynomial cut-off function). However, this crude form of 2 seems to be appropriate for the technique proposed herein. We also performed the same convergence study for KX-FEM using r inner = 4 and overlap = 1. The convergence is improved but no rate faster than 1 was obtained. This example shows that for realistic boundary conditions, higher-order terms must be considered to obtain low errors on estimated SIF values. As in the previous example, using the proposed technique with the smallest analytical domain and the smallest coupling zone yields already very low errors. This is all the more important if one wants to deal with noisy data when experimental analyses are preformed (e.g. by resorting to DIC (Réthoré et al., submitted 2009)).

CONCLUSION

The present paper proposes an alternative to partition of unity methods to incorporate enriched interpolation functions. Cracked solids are modeled by using two descriptions that are coupled over a non-zero area region using a partition of energy, namely, using a bridging domain or Arlequin method. The first description consists of standard finite elements in the outer domain, and the second one relies on Williams' series solution in the inner domain. The numerical implementation is described and two examples are proposed to assess the validity of the proposed approach, and to evaluate its performance when seeking stress intensity factors in comparison with usual X-FEM approaches. High accuracy and robustness with respect to the enriched domain size and the overlap width are obtained. Very small inner domains can be used and thus crack propagation may also be investigated without any modification of the enriched basis due to the crack curvature. As a possible extension of the proposed method, analytical crack tip fields with a cohesive law such as those derived by Karihaloo and Xiao [START_REF] Karihaloo | Asymptotic fields at the tip of a cohesive crack[END_REF] may be introduced using the proposed formalism to simulate the cohesive crack growth. Last, the developments discussed herein from a numerical standpoint are currently implemented in an experimental technique to measure displacement fields of crack bodies (by digital image correlation). In that context, the robustness of SIFs extraction is critical since pictures are noisy and the displacement field measurement itself is an ill-posed problem. The fact that the two kinematic bases (numerical and experimental) are identical allows one to avoid interpolation errors and further measurement uncertainties.

Figure 1 .Figure 2 .

 12 Figure 1. Weighting function of the finite element model 1 for the square plate with an edge crack.

Figure 3 .

 3 Figure 3. Normalized error on K I (left) and K II (right) as functions of the size of the analytical domain r inner . Top figures are plotted using a linear scale for the y-axis whereas the bottom figures use a log scale.

Figure 4 .

 4 Figure 4. Normalized on K I (a) and K II (b) as functions of width of the coupling zone overlap .

Figure 5 .

 5 Figure 5. Horizontal (a) and vertical (b) displacement field (in m) for the cracked plate in remote tension for r inner = overlap = 1.

Figure 6 .

 6 Figure 6. Normalized K I as function of the size of the analytical domain r inner (a), the size of the overlapping domain overlap (b) for n max = 2.

Figure 7 .

 7 Figure 7. Normalized K I as function of n max for r inner = 10 and overlap = 3.

Figure 8 .

 8 Figure 8. Normalized K I as function of the size of the analytical domain for different values of n max obtained using the proposed HAX-FEM (a) and KX-FEM (b).

Figure 9 .

 9 Figure 9. Normalized relative error on K I as function of the mesh size for different values of n max obtained using the proposed HAX-FEM (a) and KX-FEM (b) with r inner = overlap = 1.

Copyright q 2009 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2009) DOI: 10.1002/nme

ACKNOWLEDGEMENTS

This work was funded by the CETIM Foundation grant entitled PROPAVANFIS: 'Advanced methods for the experimental and numerical analyses of crack propagations under complex loadings'.

Contract/grant sponsor: CETIM Foundation


domain: spi spi

text:
Mohammadreza Aghaei
email: mohammadreza.aghaei@ntnu.no

Amin Moazam

Silvia Erba
email: silvia.erba@polimi.it

Mohammad Hosseini
email: mohammad.hosseini@ntnu.no

Italo Aldo

Campodonico Avendano

Muhammad-Salman Shahid
email: muhammad-salman.shahid@g2elab.grenoble-inp.fr

Ignacio Torrens-Galdiz

Giuseppe Mastandrea
email: giuseppe.mastandrea@energyatwork.it

Runar Solli
email: runar.solli@emsystemer.no

Peter Riederer
email: peter.riederer@cstb.fr

Amin Moazami
email: amin.moazami@ntnu.no

Gloria Bevilacqua
email: gloria.bevilacqua@geonardo.com

Kavan Javanroodi
email: kavan.javanroodi@byggtek.lth.se

Panayiotis Papadopoulos
email: p.papadopoulos@cyi.ac.cy

Salvatore Carlucci
email: s.carlucci@cyi.ac.cy

Vahid M Nik
email: vahid.nik@byggtek.lth.se

Collective Intelligence for Energy Flexibility -Collectief: An EU H2020 Project for Enhancing Energy Efficiency and Flexibility in Existing Buildings

Keywords: COLLECTiEF, Collective intelligence, Energy flexibility, Climate resilience, Urban energy system, Demand side management I

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The European Union (EU) has been at the forehand of international efforts to tackle the global challenge of climate change and emissions of carbon dioxide (CO2) impacts, and to deploy affordable, reliable, and modern energy services as well as to increase the share of renewable energy, according to the Sustainable Development Goals no. 7 and 13 [START_REF]7th Sustainable Development Goal (SDG7)[END_REF], [START_REF]13th Sustainable Development Goal (SDG13)[END_REF]. Rapid urbanization is exacerbating energy consumption in urban areas [START_REF] Mauree | A review of assessment methods for the urban environment and its energy sustainability to guarantee climate adaptation of future cities[END_REF]. Building and construction sector is responsible for approximately 65% of global primary energy consumption, making them the largest single energy consumer [START_REF]Clean energy for all Europeans package[END_REF], [START_REF] Nik | Using collective intelligence to enhance demand flexibility and climate resilience in urban areas[END_REF]. Improving energy performance in buildings can contribute significantly to achieving the energy and climate goals since the urban areas play a crucial role in the energy transition and the path towards sustainability.

To this end, it is crucial to search for strategies and technological solutions that reduce energy use and to cover the remaining energy demand through sustainable approaches such as increasing the penetration of variable renewable energy (VRE) [START_REF] Nik | Using collective intelligence to enhance demand flexibility and climate resilience in urban areas[END_REF]- [START_REF] Perera | Quantifying the impacts of climate change and extreme climate events on energy systems[END_REF]. Both VRE and urban energy demand are highly affected by climate conditions, causing a mismatch in demand and generation profiles [START_REF] Perera | Challenges resulting from urban density and climate change for the EU energy transition[END_REF], [START_REF] Hosseini | Improving Climate Resilience and Thermal Comfort in a Complex Building through Enhanced Flexibility of the Energy System[END_REF]. Enhancing the energy flexibility on both supply and demand sides can boost the movement towards sustainable and resilient urban energy solutions [START_REF] Hosseini | Improving Climate Resilience and Thermal Comfort in a Complex Building through Enhanced Flexibility of the Energy System[END_REF], [START_REF] Perera | Redefining energy system flexibility for distributed energy system design[END_REF].

The demand flexibility of buildings can become a major source of energy flexibility since buildings account for a large proportion of energy consumption [START_REF] Jensen | IEA EBC annex 67 energy flexible buildings[END_REF]. However, buildings with flexible loads work best if they are connected into a cluster that responds to the grid signal. In addition, larger buildings normally have a higher level of building automation coming from legacy equipment than small homes. Thus, creating a system that allows up-smarting the legacy equipment for larger buildings will have a greater impact than conventional smart home approaches targeting private homeowners. The Collective Intelligence for Energy Flexibility (COLLECTiEF) project aims to develop an energy management (EM) system based on Collective Intelligence (CI) that allows easy and seamless integration of legacy equipment into a collaborative network within and between existing buildings and urban energy systems. Fig. 1 depicts the conceptual design of the COLLECTiEF project. The main objectives and contributions of the COLLECTiEF project are summarized as follows: The main core of the COLLECTiEF project is shaped around enhancing the existing CI algorithm [START_REF] Nik | Using collective intelligence to enhance demand flexibility and climate resilience in urban areas[END_REF] and developing it further to address the complex interactions between buildings at the urban scale and to generate optimum control strategies at the building and cluster levels. The CI algorithm (Obj.1) will shape the design of the sensor network (Obj.4), the user inputs/interactions and the digital dashboards (Obj.5), and the system components at the edge and cluster nodes (Obj.2). The solutions will be tested in the pilots and ameliorated during the demonstration phase (Obj.3). The final goal is to have the system qualified (TRL8) and ready for commercialization with a new business model (Obj.6). Objective's relations are depicted in Fig. 2.

II. METHOD AND IMPLEMENTATION

COLLECTiEF implements and evaluates an interoperable and saleable EM system based on Collective Intelligence (CI). The concept of CI connects the following three pillars in the COLLECTiEF project: (1) Resilience, demand-response, and energy flexibility on the level of multiple buildings; (2) Smart buildings and smart components; (3) Occupant-centric operation.

A. Methodology

COLLECTiEF develops software and hardware packages to smart up buildings of different typologies common in Europe and their legacy equipment, sort out consumer-centric and industry-centric challenges, and maintain simple and robust communication with the energy grid. The simple communication acts as a stimulus to help buildings sense the grid's need for flexibility and accordingly within the limits of their local constraints and users' preferences. The overall methodology of the project includes the following steps and approaches: (1) Development of control strategies, algorithms, and frameworks for demand side management; (2) Adaptation of the technical components to the use case; (3) Evaluation of the COLLECTiEF algorithms and control strategies as well as the complete COLLECTiEF solutions on the District MOdeller and SIMulator (DIMOSIM) platform [START_REF] Garreau | District MOdeller and SIMulator (DIMOSIM)-A dynamic simulation platform based on a bottom-up approach for district and territory energetic assessment[END_REF]; (4) Testing in a small-scale demonstration site; (5) largescale demonstration of the COLLECTiEF solutions in 13 pilot buildings located in three countries namely, Norway, Italy, and Cyprus; (6) Development of business models and exploitation of the results; [START_REF] Perera | Quantifying the impacts of climate change and extreme climate events on energy systems[END_REF] Capacity building and stakeholder engagement. Fig. 3 shows the value chain of COLLECTiEF and the involved partners in each stage.

B. Implementation and work packages

The COLLECTiEF consortium consists of 15 beneficiaries from universities, institutions, manufacturing and construction companies, and municipal sectors from six countries across Europe. COLLECTiEF comprises seven work packages (Fig. 4) aligned with each other which aim at developing, implementing, testing, and evaluating the proposed EM system based on CI. WP1 coordinates the COLLECTiEF project management and administration as well as supervises the activities carried out by all partners with defined tasks and according to the Grant Agreement and planned action plan. WP2 focuses on the development of algorithms and communication logic of the COLLECTiEF EM system. The partners are developing, testing and enhancing existing and new solutions using a virtual test bed, DIMOSIM [START_REF] Garreau | District MOdeller and SIMulator (DIMOSIM)-A dynamic simulation platform based on a bottom-up approach for district and territory energetic assessment[END_REF], and in a real environment at the G2Elab. WP3 aims at developing and implementing the prototype of the COLLECTiEF network, including hardware and software technologies for demonstration in a real smallscale environment in the G2Elab. WP4 focuses on the Edge Node, Cluster Node, occupant-centric fusion sensor network, IoT Operating System, and Human-Building interface. WP4 has worked on the template for identifying user parameters and will adapt it to COLLECTiEF requirements. WP5 deals with the definition of the measurement and verification protocol and the evaluation of the impacts given by the implementation of the COLLECTiEF solutions on building smartness, energy saving and flexibility, indoor environmental quality, and climate resilience. WP5 interfaces with all the work packages for the eventual completion and qualification of the COLLECTiEF system at TRL 8. WP6 aims to carry out an effective exploitation process of the project results, emphasizing results with commercial and industrial exploitation potential. For this, all COLLECTiEF partners are involved. Finally, WP7 focuses on the project's proper visibility and spreads pertinent information on its goals, activities, and results to the relevant stakeholders and scientific communities, thereby fostering the engagement of the target audience including end users and stakeholders in the COLLECTiEF activities.

C. Project structure

COLLECTiEF project originates from the identification of available knowledge and technologies within the consortium that have been validated and demonstrated in previous projects and which fall under the technology readiness level (TRL) 5 or 6. These technologies allow the creation of a complete CI-based energy network and within the 4-year timeline of the project, the main aim is to qualify the COLLECTiEF solutions in four real applications, hereafter called "DEMO", which will bring the overall TRL to 8. Therefore, project management is divided into two main phases, namely (1) before the start of the DEMO, (2) after the start of the DEMO.

This classification aims at organizing and monitoring the defined tasks as well as at providing a better overview of the project's progress. The main activities are summarized in Table I.

A. Activities performed and objectives progress

The COLLECTiEF consortium performed numerous dedicated tasks and planned activities as well as obtained several achievements according to milestones and action plan by month 24. Additionally, the WPs and Tasks leaders submitted up to 22 deliverables to EU Commission by month 24. The main developments and achievements are described as follows:

The first version of the EM algorithms at the building (Edge node) and cluster (Cluster Node) levels has been developed based on combining COLLECTiEF Intelligence and Reinforcement Learning with promising results. The algorithms are being fine-tuned to be implemented at the edge and cluster nodes. The major indicators in developing the algorithms have been energy demand and indoor thermal comfort so far. Depending on the need, the indicators and objectives will be updated in future versions of the algorithms. Representative weather data sets have been generated considering several future climate scenarios. The weather data sets have been used in running the algorithms and setting the control strategies to assess their impact in adapting to climate variations and decreasing energy demand while guaranteeing adequate indoor thermal comfort conditions for occupants. Different approaches for modelling and assessing indoor thermal comfort have been developed, which will be integrated further into the control algorithms. Moreover, the algorithms are being developed to include other factors/indicators, such as energy price. A functioning platform for data sharing and integration, algorithm development, and co-simulation has been created in GitHub, to which all the relevant users have access. The virtual test bed using DIMOSIM for algorithms and controls has been developed and tested at G2ELab. The Border Router has been adapted on the Raspberry Pi 4 open board to have a platform capable of supporting the computational load foreseen by the project. In addition, the legacy and already-in-place devices in the different pilot buildings were investigated, and the smart devices for interfacing with the field, such as smart plugs and smart thermostats, have been identified. Significant steps forward have been made to achieve this goal, mainly from the hardware point of view. Particularly, the reference open board (Raspberry Pi 4), capable of integrating the edge computing functions envisaged by the project, has been selected based on which the entire hardware platform of the Edge node-BRiG (Border Router + iGateway) was built. The specifications for the use of Raspberry Pi 4 were defined, and the porting of the Border Router firmware (previously based on the lessperforming BeagleBone open board) on Raspberry Pi 4 was implemented. The electronic schematic and PCB of BRiG were developed in this first period. Afterward, based on them, the different components were mounted on the BRiG prototypes and subsequently functionally tested.

A Measurement & Verification (M&V) plan has been defined to evaluate the implementation of the COLLECTiEF solutions targeting not only energy savings and indoor environmental quality (IEQ) but also including the methods to measure and verify the load reduction and management improvement, resulting from the utilization of energy flexibility strategies, and the enhanced climate resilience. The Deliverable D5.1 -Performance Measurement & Verification Protocolconcepts and methods for performance evaluation of COLLECTiEF solutions describes the process of M&V in COLLECTiEF, presents the methodology, and the key performance indicators (KPI) selected to assess the impacts on the smart readiness of buildings, energy performance, indoor environmental quality, energy flexibility, and climate resilience. Furthermore, Deliverable D5.2 -Ongoing performance evaluation of the COLLECTiEF system implemented in the pilot cases presents in detail the methodologies to assess the indoor environmental quality in the buildings and shows the progress on the monitoring and preliminary analysis of the data collected in the pilot buildings (see Fig. 6), pointing out lessons learnt and mitigation strategies. The smart-ready technologies (e.g., smart meters, smart thermostats, sensors, and IoT devices) to be installed in the different pilot buildings have been selected, and it has been initiated the activities related to the development of communication, development, and testing in simulated environments of the algorithms able to connect to the use of digital technologies, and the algorithms for energy, comfort optimization, and flexibility management are guaranteeing the protection of privacy and security of smart buildings.

We finalized the implementation of the occupant-centric fusion sensor network in the pilot buildings. For most parts, the implementation has been successful, and the communication between the sensors and the BRiG device is operational and running. The sensor placement has been executed in a manner that minimizes interference with the user environment, while effectively characterizing the indoor environments see Fig. 7. Besides, the development of the Human-Building Interface is under further development, see Fig. 7. The comprehensive market analysis has been finalized, providing a foundation for the development of business models that can meet the market requirements for the effective exploitation of the COLLECTiEF solutions. Particularly, R2M provided the analysis of the contexts where the business models operate through a PESTLE analysis (see Fig. 9) and the assessment of the market readiness across various EU countries. Additionally, the analysis also included evaluations of successful examples and potential competitors in the COLLECTiEF market space, highlighting their main functionalities, provided services, commercialization strategies, and adopted business models. This information is valuable for benchmarking and serves as a reference point for developing business models for exploiting COLLECTiEF solutions. Further discussions have been carried out to guide the development of the business.

From the beginning of the project, we have been actively engaged in efforts to enhance the visibility of the project. Regarding dissemination, first, the visual identity was designed. It was followed by the launch of the website. A set of materials such as leaflets, rollups, banners, and posters were created. A stakeholder database was set up for the partners to mobilize their network. When it comes to the engagement with the pilot sites, the first workshops were held physically in Milan, Italy, and in Ålesund, Norway. Social media accounts on popular platforms such as Twitter, LinkedIn, and YouTube were launched. As of M12, WP7 started to publish a newsletter within the networking activities. In addition, the COLLECTiEF project first collaborated with the SMART2B project and then with three additional Horizon2020 projects, creating a cluster called EU4BET. EU4BET participates in the Horizon Results Booster service. As part of the capacity building and sustainability actions, the project joined forces with the SmartBuilt4EU project and participated in the workshop about Smart Buildings in Nice, France.

B. Expected impacts and progress

The COLLECTIEF project foresees different impacts in relation to the various domains which are affected by the system implementation (building smartness, energy, indoor environmental quality, climate flexibility and resilience). In particular, the project aims to reach the following six main results:

• Impact 1: Upgrade existing buildings to higher smartness levels

The COLLECTiEF system aims at upgrading the existing pilot buildings with at least one level of smartness according to Smart Readiness Indicator (SRI) scoring methods, independently from the starting smart readiness level. The assessment of the smartness level requires carrying out visual inspections of the pilots and collecting technical data and information. The SRI is calculated in each pilot building before and after the implementation of the systems.

• Impact 2: Reduction in energy consumption and costs

Depending on the building's thermal quality and climate zone, 9-27% relative energy reduction is estimated for buildings after installation of the COLLECTiEF system and upgrading by one level of smartness. The assessment will be carried out in a hybrid way: (1) by the use of measurement of energy consumption and environmental variables (i.e. independent variables such as weather) at the pilot buildings and, thus limited to the boundary of the selected apartments/rooms. (2) By the use of calibrated models of the pilot sites on which the COLLECTiEF algorithms will apply their optimized operation. This will allow extending the impact assessment from the zone to the building/cluster of buildings level. Further, energy cost savings will be evaluated.

• Impact 3: Primary Energy Savings triggered by the project (in GWh/year) An average of 16% reduction in primary energy use of the pilot buildings, equivalent to 0.93 (GWh/year), is estimated after the implementation of the COLLECTiEF systems. For the evaluation of primary energy savings, updated primary energy factors (PEF) will be considered.

• Impact 4: Investments in sustainable energy triggered by the project (in million EUR)

The implementation of the COLLECTiEF system is expected to generate up to a 52% increase in the selfconsumption rate. Among the 14 pilot buildings, only two buildings include renewable energy production from photovoltaics systems. Thus, to get a broader view of the potential impact of increasing self-consumption, some calibrated pilot building models will be virtually equipped with REN energy systems.

• Additional Impact 5: User satisfaction

The project aims at increasing user satisfaction with respect to some of the components of the indoor physical environment by at least 15%, thanks to the implementation of smart occupant-centric sensoring and self-learning control algorithms. User satisfaction will be assessed through two questionnaires: (i) the Post-Occupancy Evaluation (POE) questionnaire and (ii) the Satisfaction questionnaire, that aim to acquire information about occupant's perception of living conditions and to collect data about their satisfaction. A comparison will be conducted between the baseline period, where no control action has been taken on the pilot buildings, and the implementation period, where the COLLECTiEF system will be applied.

• Additional Impact 6: Increased climate flexibility and resilience in urban areas, inducing larger integration of renewable generation.

Fluctuations in the energy demand are greatly influenced on climate variations, particularly during extreme climate events. Climate change can exacerbate this situation by causing more frequent and severe extreme weather events. This results, for example, in having more frequent heatwaves and cold snaps. Supplying enough energy for such times is challenging, since the demand profiles change significantly with considerable loads on the supply system. With the energy management approach that is suggested in COLLECTiEF, it is possible to pass these extreme conditions safely since the COLLECTiEF system enables enhancing the energy flexibility in the network through a collaborative response of buildings and energy systems to climate shocks, resulting in improved climate resilience in urban areas.

IV. CONCLUSION

COLLECTiEF is an EU-funded H2020 project, running for 4 yearsfrom 2021 to 2025. COLLECTiEF aims to enhance, implement, test, and evaluate an interoperable and saleable energy management system based on CI that allows easy and seamless integration of legacy equipment into a collaborative network within and between existing buildings and urban energy systems with reduced installation cost, data transfer, and computational power while increasing data security, energy flexibility and climate resilience.

COLLECTiEF project commits addressing reliable and practical solutions for the challenges of climate change impacts and renewable energy penetration by enhancing the energy flexibility and climate resilience through a collective intelligence (CI) approach.

The COLLECTiEF consortium performed numerous dedicated tasks and planned activities as well as obtained several achievements according to COLLECTiEF's milestones and action plan by month 24. Furthermore, COLLECTiEF's partners submitted up to 22 deliverables to EU Commission by month 24.

During the next two years, the COLLECTiEF consortium will focus on further developing the software/hardware and initiating the large-scale demonstration as well as testing and qualifying COLLECTiEF's solutions. We will also work on the business model and the exploitation plan in order to roll up the commercialization after completing the project.

Fig. 1 .

 1 Fig. 1. The conceptual design of the COLLECTiEF project

Fig. 2 .

 2 Fig. 2. The objectives relations in the COLLECTiEF project

Fig. 3 .

 3 Fig. 3. The COLLECTiEF value chain

Fig. 4 .

 4 Fig. 4. The overview of work packages in the COLLECTiEF project

Fig. 5

 5 demonstrates the schematic diagram of the implementation of algorithm. The flowchart of the code is presented on top right of the figure. The G2ELab's dashboard is used to visualize the data (i.e. sensor and actuator values) in real-time.

Fig. 5 .

 5 Fig. 5. Workflow of the small-scale test in G2ELab. bottom-left: BMS panel in the rooms, top-right: flowchart of the developed algorithm, bottomright: real-time visualization platform to demonstrate the performance of the code.

Fig. 6 .

 6 Fig. 6. Example of data representation regarding TVOC concentration in one of the monitored pilot buildings

Fig. 7 .Fig. 8 .

 78 Fig. 7. Picture of the installed IEQ monitoring system and Post Occupancy Evaluation QR code in one of the pilot buildings in Norway Fig. 8. Picture of the installed IEQ monitoring system and Post Occupancy Evaluation QR code in one of the pilot buildings in Norway. The concept and architecture of the Human-Building Interface were defined and described in Deliverable D3.8: Fully Integrated Dashboard for Cluster Node and Human-Building Interface for the Edge Node (first version). Moreover, the prototype has been developed by including different functions, data insights, and authentication mechanisms.

Fig. 8 .

 8 Fig. 8. The COLLECTiEF's Human-Building Interface homepage

Fig. 9 .

 9 Fig. 9.PESTLE analysis (main topics covered)

TABLE I .

 I THE MAIN ACTIVITIES AND MILESTONES

Before the Start of the DEMO (M1-M24) After the Start of the DEMO (M24-M48)

	Definition of the
	Measurement and	Testing in small-scale test bed
	Verification Protocol	completed
	completed
	Involvement of the users and occupants of the pilot buildings conducted	Demonstration phase in large-scale completed
	Project communication kit	Involvement of the users and
	developed	occupants of the pilot buildings
	Data management plan	Stakeholder workshops conducted
	Control algorithms ready for testing	Business model defined, showing commercial feasibility
	Edge Node, Cluster Node,
	human-building interface, and occupant-centric sensor network ready for	Capacity training material ready for use
	testing
	Installation in pilots completed	Evaluation of the impacts. System complete and qualified (TRL8)
	Market Analysis
	(Stakeholder and
	competitors analysis, Regulatory framework,	Project completed
	market environment
	research) completed
	III. PROJECT PROGRESS AND ACHIEVEMENT

ACKNOWLEDGMENT

COLLECTiEF project has received research funding from European Union's H2020 research and innovation programme under Grant Agreement No 101033683.