Yue Sun

Emmanuel Roubin

Jianfu Shao

Jean-Baptiste Colliat

Meso-scale Finite Element modeling of the Fracture Process Zone evolution for concrete

Keywords:

The so-called Fracture Process Zone (FPZ) is well-known as a precursor of macroscopic failure for concrete. Moreover, it is a reasonable assumption that this FPZ extension is partly governed by the underlying random distribution of the material meso-structure. Therefore, meso-scale (millimeter scale) fracture simulations are effective in improving our understanding of this relationship. In this paper, our objective is to simulate the extension of cracks within the FPZ and to investigate the influence of spatial heterogeneity distribution. To this end, we propose to employ the Enhanced Finite Element Method (E-FEM) as our simulation tool. In the proposed simulation method, two types of discontinuities are incorporated within the elements: strong discontinuities serve to simulate actual cracks, while weak discontinuities allow for an explicit representation of heterogeneities. The combination of these two enhancements makes the model well-suited for investigating the FPZ at the meso-scale for being able to simulate complex crack features and interactions between cracks and aggregates. The simulation captures many crack patterns in the FPZ, such as diffused micro-cracks, macro-crack coalescence, multi-branching, and crack closures. The size effect of inclusions in the FPZ, which is widely observed in experiments, is also realized in simulations. The comparison with the experimental observations shows good consistency.

a-AE

. The integration of both discontinuities allows for the simulation of complex crack features, including spontaneous crack initiation and multi-branching, as well as crack-aggregate interactions, such as debonding on the heterogeneity interface.

Introduction

It is a general consensus that the macroscopic failure behaviors of quasi-brittle materials, such as concrete and rocks, are strongly associated with the development of cracks at a fine scale. The formation of a macro-crack is a gradual process that begins with diffused micro-cracking and involves matrix-aggregate debonding, micro-crack coalescence (branching of cracks), and crack closures [START_REF] Bažant | Measurement of characteristic length of nonlocal continuum[END_REF][START_REF] Wittmann | Fracture process zone in cementitious materials[END_REF][START_REF] Mihashi | Correlation between characteristics of fracture process zone and tension-softening properties of concrete[END_REF], see the sketch Fig. 1. This whole collection of mechanisms at the fine scale is dissipative and responsible for stress transfer. The subcritical zone that experiences micro-cracks is referred to as the Fracture Process Zone (FPZ), which can be observed in the vicinity and ahead of the macro-crack. Various experimental techniques have been developed to investigate the FPZ characteristics and track the crack propagation, such as Acoustic Emission (AE) measurements [START_REF] Landis | Micro-macro fracture relationships and acoustic emissions in concrete[END_REF][START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF][START_REF] Muralidhara | Fracture process zone size and true fracture energy of concrete using acoustic emission[END_REF][START_REF] Haidar | Correlation between the internal length, the fracture process zone and size effect in model materials[END_REF], X-ray tomography techniques [START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF][START_REF] Landis | Microstructure and fracture in three dimensions[END_REF], and Digital Image Correlation techniques [START_REF] Choi | Measurement of deformations on concrete subjected to compression using image correlation[END_REF][START_REF] Corr | Digital image correlation analysis of interfacial debonding properties and fracture behavior in concrete[END_REF][START_REF] Tejchman | Calculations of fracture process zones on meso-scale in notched concrete beams subjected to three-point bending[END_REF][START_REF] Wu | An experimental investigation on the fpz properties in concrete using digital image correlation technique[END_REF].

To numerically simulate the failure behaviors of material with cracking, multiple models have been proposed from a macroscopic point of view. In the concept of Representative Volume Element (RVE), the material can be considered as homogeneous. Under this hypothesis, various concepts that consider fracture mechanics have been proposed, such as the Linear elastic fracture mechanics (LEFM) [START_REF] Griffith | Vi. the phenomena of rupture and flow in solids[END_REF][START_REF] Wells | Crack propagation symposium proceedings[END_REF]Rice, 1968;Rice et al., 1968), the cohesive crack model [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF][START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF][START_REF] Veselỳ | An application for the fracture characterisation of quasi-brittle materials taking into account fracture process zone influence[END_REF], the crack band theory [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF], and the cohesive zone model [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF], just mention a few. These macro-scale models are robust and sophisticated for simulating various material behaviors. However, the definition of material parameters, such as fracture energy, is usually strongly related to the FPZ size, which is predominantly determined by the meso-structure of the material. Therefore, meso-scale models are proposed for a detailed analysis of the fracture process in the FPZ. Meso-scale models can be briefly divided into three categories based on their treatment of the displacement field: i) Continuum models: the displacement field is assumed to be continuous; ii) Discrete models: the displacement field is considered to be established by rigid body motion of discrete units; iii) Discontinuous models: discontinuity is introduced into the displacement field.

For continuum models, one can incorporate the meso-structure directly into the finite element framework using randomly placed aggregates [START_REF] Rodrigues | 3d concurrent multiscale model for crack propagation in concrete[END_REF][START_REF] Zhou | 3d mesoscale finite element modelling of concrete[END_REF][START_REF] Tejchman | Calculations of fracture process zones on meso-scale in notched concrete beams subjected to three-point bending[END_REF], random fields [START_REF] Georgioudakis | Stochastic failure analysis of structures with softening materials[END_REF][START_REF] Vořechovskỳ | Computational modeling of size effects in concrete specimens under uniaxial tension[END_REF], or X-ray tomography images. The latter technique makes it possible to incorporate real meso-structure into the finite element framework [START_REF] Trawiński | A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and x-ray 𝜇ct images[END_REF][START_REF] Huang | 3d meso-scale fracture modelling and validation of concrete based on in-situ x-ray computed tomography images using damage plasticity model[END_REF]. In some of these models, one layer of cement-aggregate interface element is used to simulate the Interfacial Transition Zone (ITZ). The mechanical behavior of the material can be simulated using damage-based constitutive models [START_REF] Zhou | Mesoscopic investigation of size effect in notched concrete beams: The role of fracture process zone[END_REF][START_REF] Rodrigues | 3d concurrent multiscale model for crack propagation in concrete[END_REF], plasticity theory [START_REF] Zhang | 3d mesoscale fracture analysis of concrete under complex loading[END_REF], or phase-field models [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF]. However, using the classical finite element method usually requires a conforming mesh strategy, which can lead to high computational costs and a limitation that cracks can only propagate along element boundaries. As a result, the crack trajectory is strongly influenced by the mesh alignment [START_REF] Tijssens | Numerical simulation of quasi-brittle fracture using damaging cohesive surfaces[END_REF][START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: survey of progress[END_REF][START_REF] Zhou | Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency[END_REF].

Multiple discrete models have been proposed to investigate material fracture behavior at the meso-scale. The simplest of these is lattice models [START_REF] Luković | A 3d lattice modelling study of drying shrinkage damage in concrete repair systems[END_REF][START_REF] Pan | Lattice modeling of early-age behavior of structural concrete[END_REF][START_REF] Pan | A review of lattice type model in fracture mechanics: theory, applications, and perspectives[END_REF], where rigid aggregates can ideally be implemented using lattice elements. The heterogeneous meso-structure of the material can be obtained by using the X-ray scan images [START_REF] Nitka | Modelling of concrete fracture at aggregate level using fem and dem based on x-ray 𝜇ct images of internal structure[END_REF], or by proposing spacial variations of material parameters, such as the tensile strength and fracture energy [START_REF] Grassl | Random lattice-particle simulation of statistical size effect in quasi-brittle structures failing at crack initiation[END_REF][START_REF] Eliáš | Stochastic discrete meso-scale simulations of concrete fracture: Comparison to experimental data[END_REF][START_REF] Xenos | Calibration of nonlocal models for tensile fracture in quasi-brittle heterogeneous materials[END_REF]. Recent studies have aimed to simulate the Fracture Process Zone (FPZ) in concrete materials [START_REF] Grassl | Meso-scale modelling of the size effect on the fracture process zone of concrete[END_REF][START_REF] Grassl | Meso-scale approach to modelling the fracture process zone of concrete subjected to uniaxial tension[END_REF]. The methodology employed in these studies is similar to the one used in the current work, where a heterogeneous meso-structure is generated through the placement of randomly distributed spherical aggregates with a minimum inter-aggregate distance. The observed crack patterns in these studies and the present work show similarities, including the presence of diffused micro-cracks in front of and adjacent to the macro-crack, as well as crack closures.

In the context of discontinuous models, two main strategies can be found in the literature that are distinguished by their support of additional degree of freedom: nodal enrichment, such as eXtended Finite Element Method (X-FEM) [START_REF] Kumar | Mixed mode crack growth in elasto-plastic-creeping solids using xfem[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF][START_REF] Benvenuti | Simulation of finite-width process zone in concrete-like materials by means of a regularized extended finite element model[END_REF][START_REF] Li | Fracture process in cortical bone: X-fem analysis of microstructured models[END_REF]; and elemental enrichment, such as Enhanced Finite Element Method (E-FEM). The discontinuous model enables discontinuity in the displacement field, and thus can represent cracks explicitly. In order to numerically investigate the fracture process in the FPZ at the meso-scale, the method used in this study is an E-FEM model.

The E-FEM model has been selected for this study due to its simple yet effective way to model the heterogeneous The structure of this paper is composed of six sections. After the introduction in the first section, Section 2 provides a kinematic description of the weak and strong enhancements. Then Section 3 concentrates on the mechanical behavior of the model at the local scale, where a mode-I crack opening-closing E-FEM model [START_REF] Sun | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF] is selected to simulate the evolution of the FPZ in concrete. This section addresses the four main phases of cohesive behavior on the strong discontinuity. Next, Section 4 outlines the resolution strategy of the model within the finite element framework. Subsequently, in Section 5, the present numerical model is used to simulate a notched cube with different inclusion sizes, with material parameters selected to match the corresponding experimental results by Otsuka [START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF]. The focus is made on the crack feature evolution in the FPZ under different loading levels and the effect of the aggregate size. Finally, a brief conclusion and discussion are given in Section 6 to close the paper.

Nomenclature

𝑩

Kinematics enhancements at the meso-scale

This section starts by summarizing the basic notations that are employed in this paper. Then the kinematic descriptions and weak forms of the enhancements, the weak and strong discontinuities, are presented subsequently.

A sketch of the reference domain Ω is plotted in Fig. 2. It is a heterogeneous material which consists of a heterogeneity (dark gray part) and matrix (light gray part). By applying a non-adapted mesh strategy [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF][START_REF] Kamel | Comparison of advanced discretization techniques for image-based modelling of heterogeneous porous rocks[END_REF], an unstructured mesh is projected to the heterogeneous material. Two kinds of enhancement are generated after the discretization:

• The weak discontinuity is embedded in the elements that are passed by the heterogeneity interface  𝜀 . It brings a jump in the strain field [|𝜺|] to the element because of the difference between the materials' elastic properties, see Fig. 3(a).

• The strong discontinuity is embedded in the elements that are cut by a crack  𝑢 . It brings a jump in the displacement field [|𝒖|] to the element, see Fig. 3(b).

The both discontinuities are oriented. An local orthogonal basis (𝒏, 𝒎, 𝒕) is built based on the normal vector 𝒏 on the discontinuity surface. In the case of a debonding at the heterogeneity interface, the weak and strong discontinuities are independent and can be merged to the same element, see Fig. 3(c). The strain field writes as (Benkemoun, Hautefeuille, Colliat and Ibrahimbegovic):

𝜺(𝒙) ∶= ∇ sym 𝒖(𝒙) + ε(𝒙) + ε(𝒙), (1)
where

∇ sym ∶= 1 2 [∇(•) + ∇ 𝑇 (•)]
is the symmetric gradient operator and ū is the regular part of the displacement field.

Notation • / • presents the weak / strong discontinuity, respectively.

Kinematic description and weak form of the weak discontinuity

As presented previously, the weakly enhanced elements are located at the heterogeneity interfaces of the material.

The normal vector of the discontinuity interface 𝒏 and the corresponding local basis (𝒏, 𝒎, 𝒕) is defined by the time of the mesh projection. Since the displacement field is continuous across the discontinuity surface, the magnitude of the jump in the strain field [|𝜺|] can be entirely defined by three unknown values [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]. The weakly enhanced part of the displacement field ũ writes as [START_REF] Gurtin | The linear theory of elasticity[END_REF][START_REF] Ibrahimbegović | Strong coupling methods in multi-phase and multi-scale modeling of inelastic behavior of heterogeneous structures[END_REF]:

{[𝜀] 𝑛 , [𝜀] 𝑚 , [𝜀] 𝑡 } T (
ũ(𝒙) = Θ𝒏 ⋅ (𝒙 -𝝃)([𝜀] 𝑛 𝒏 + [𝜀] 𝑚 𝒎 + [𝜀] 𝑡 𝒕) with Θ = { Θ + ∀𝑥 ∈ Ω + 𝑒 Θ -∀𝑥 ∈ Ω - 𝑒 , ( 2
)
where 𝝃 represents the position of interface surface, 𝒏 ⋅ (𝒙 -𝝃) can be seen as a signed distance to the discontinuity surface, and Θ is a still undefined function whose explicit expression is given later.

Correspondingly, the enhanced strain field can be calculated by taking the symmetrical gradient of ũ:

ε = ∇ sym ( ũ) = Θ ( [𝜀] 𝑛 𝒏 ⊗ 𝒏 + [𝜀] 𝑚 2 (𝒏 ⊗ 𝒎) sym + [𝜀] 𝑡 2 (𝒏 ⊗ 𝒕) sym
) .

(3)

Kinematic description and weak form of the strong discontinuity

An element is strongly enhanced if it is crossed by a crack. In this case, a jump in the displacement field is embedded in the element, the kinematic description of which writes as (Oliver, 1996b): is considered as a constant function over Ω 𝑒 . The corresponding strain field can be obtained by taking the symmetric gradient of the displacement field Eq. 4 [START_REF] Simo | A new approach to the analysis and simulation of strain softening in solids[END_REF][START_REF] Oliver | On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations[END_REF]:

𝒖 = û + (  𝑢 -𝜑 𝑒 )[|𝒖|], (4)
𝜺 = ∇ sym 𝒖 = ∇ sym û ⏟ ⏟⏟ regular -([|𝒖|] ⊗ ∇𝜑 𝑒 ) sym ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ bounded enhancement ε𝑏 + 𝛿  𝑢 ([|𝒖|] ⊗ 𝑛) ⏟⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞ ⏟ unbounded enhancement ε𝑢 . ( 5
)
The term 𝛿  𝑢 is the Dirac-delta distribution which comes from the gradient of the Heaviside function.

Mechanical behaviors at the meso-scale

In this section, the focus is made on the mechanical behaviors of the model at the local scale, especially the fracture behaviors related to the cracks. As pointed out previously, cracks are presented by strong discontinuities in this model, namely the strong discontinuity analysis [START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids[END_REF]. It aims at ensuring the consistency of the model in the presence of a jump in the displacement field, which will result in an unbounded enhancement as presented in Eq. 5. In this paper, the Discrete Strong Discontinuity Approach (DSDA) [START_REF] Oliver | On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations[END_REF][START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures: Part i: theoretical formulation and numerical implementation[END_REF][START_REF] Kucerova | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures: Part ii: identification from tests under heterogeneous stress field[END_REF] is used, which discretizes the strongly enhanced element into two parts:

• On the discontinuity interface, a traction vector 𝑻 is presented, which is continuous and serves as a bright linking the two sub-volumes through the discontinuity surface.

• Beyond the discontinuity interface, the bulk volumes is pure elastic and can be fully described by the Hooke's law.

In this framework, the nonlinear behaviors are entirely concentrated on the discontinuity. The governing law of these behaviors is depicted by a relationship between the traction vector 𝑻 and the crack opening [|𝒖|], namely the traction-separation law [START_REF] Wells | Discontinuous modelling of strain localisation and failure[END_REF]Dias-da Costa, Alfaiate, Sluys and Júlio, 2009a,b). In existing E-FEM models, several types of cracking modes have been proposed and put into literature, such as the mode-I traction-opening model [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF][START_REF] Sun | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF], the mode-II sliding-opening model [START_REF] Hauseux | Fe modelling with strong discontinuities for 3d tensile and shear fractures: application to underground excavation[END_REF][START_REF] Yue | Fe modeling of concrete with strong discontinuities for 3d shear fractures and comparison with experimental results[END_REF], and the mix-mode of the mode-I and mode-II is also under investigation (Laborin, Malecot, Roubin and Daudeville). In this paper, considering that the mode-I crack mode is fairly a common choice for modeling FPZ, the mode-I model with the crack opening-closing mechanism [START_REF] Sun | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF] is selected. Hence, in the subsequent four parts of the section, the localization criterion, the opening law, and the closing law are presented. These governing laws describe the fracture behaviors of a quasi-brittle material at the local-scale.

Localization criterion

The localization criterion determines the moment of failure initiation and the orientation of the crack. It is a stressbased criterion which writes as:

Φ l = 𝜎 eq -𝜎 y , (6)
where the equivalent stress 𝜎 eq is compared with a critical value 𝜎 y (local material property). A negative value of Φ l means an elastic behavior while a positive value presents that the localization occurs in the element. Depending on the localization takes place in a homogeneous element or a heterogeneous element, two situations are considered:

• In the case of homogeneous element, the equivalent stress is defined as the major principal stress 𝜎 𝐼 . Correspondingly, the crack orientation 𝒏 is determined by the major eigenvector.

• In the case of heterogeneous element, where a weak discontinuity is already present in the element, we consider a hypothesis that the strong discontinuity takes place on the same surface [START_REF] Ibrahimbegovic | Embedded discontinuity finite element method for modeling of localized failure in heterogeneous materials with structured mesh: an alternative to extended finite element method[END_REF][START_REF] Benkemoun | Failure of heterogeneous materials: 3d meso-scale fe models with embedded discontinuities[END_REF]. The equivalent stress is defined as the projection of traction vector on the pre-defined interface: 𝜎 eq = 𝑻 ⋅ 𝒏. This choice is made on the experimental observations that the cracks tend to appear at the interface between mortar and aggregates because the Interfacial Transition Zone (ITZ) [START_REF] Scrivener | The interfacial transition zone (itz) between cement paste and aggregate in concrete[END_REF][START_REF] Jebli | Experimental characterization of mechanical properties of the cement-aggregate interface in concrete[END_REF] usually carries a higher porosity, which make it a weaker zone for concrete.

Traction-opening criterion

After the initiation of the crack, we consider a cohesive traction-opening procedure to the crack, see Fig. 4. As the crack opening value [|𝒖|] increases, the cohesion linking the two bulk sub-volumes decreases. A classical softening criterion is used to depict this mechanical behavior at the discontinuity:

Φ o = 𝜎 eq -𝜎 y exp ( - 𝜎 y  op [𝑢]
)

⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝑞 o . ( 7
)
In this equation, two local material parameters are applied: the critical tensile stress 𝜎 y and the fracture energy  op (the gray zone in Fig. 4). The latter parameter presents the necessary energy to completely open a crack of 1 𝑚 2 . Upon this consideration, a larger fracture energy  op presents a more ductile material behavior.

An analytical solution is available for the non-linear resolution of Eq. 7, readers can refer to section 4.2 for more details.

Closing criterion

In case of a non-proportional loading is applied to an element that carries a opened crack, the tensile stress on the discontinuity interface could reduce progressively or even turns to compression. For a model without closure mechanism, no energy can be dissipated in this phase since the opening value [|𝒖|] remains constant, see Fig. 5(a).

For a model with closure [START_REF] Sun | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF], on the other hand, permits the closing of the opened cracks and the energy dissipation. These are crucial ingredients for the failure behaviors of quasi-brittle materials under non monotonic loadings such as cyclic compression tests. The closing procedure includes several stages, see Fig. 5(b).

As the firs moment of the unloading, the discontinuity interface is still under tension. The first phase therefore corresponds to a decrease of the tensile stress on the discontinuity interface, see stage (1) to stage (2) in the figure.

The opening value [|𝒖|] will remain unchanged before a fully release of the elastic energy that is stored in the bulk volumes.

Subsequently, as the loading continue decreases, a compressive stress is applied to the discontinuity interface, and the closing criterion is triggered. We consider here an exponential form for the closing criterion. It leads to an increasing required compressive load to close the crack. It is impossible to completely close the crack, which means that the opened crack cannot be healed. Upon consideration, the closing criterion writes as:

Φ c = -𝜎 eq +  cl [𝑢] max ln ( [𝑢] [𝑢] max ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝑞 c , ( 8
)
where  cl is the necessary energy for completely closing the crack, and [𝑢] max is the largest opening value. Here, an assumption is adopted that the fracture energy for closing (light gray part in the left in Fig. 5

 cl ∶=  ′ op .
This part of the energy is calculated as:

 ′ op = ∫ [𝑢] max 0 𝜎 eq ([𝑢]) 𝑑[𝑢] =  op ( 1 -exp ( - 𝜎 y [𝑢] max  op )) =  cl . ( 9
)
This assumption makes a larger opened element more difficult to close completely. Furthermore, it saves us from adding additional material parameters.

Re-opening criterion

Let us now consider that the closed element is now reloaded. As presented in the closing phase, the first mechanical behaviors is pure elastic, see stage (3) to stage (4) in Fig. 5(c). The compressive stress on the discontinuity surface progressively decreases and the tensile stress increases until reaches the critical tensile stress. Then the re-opening criterion is triggered (stage (4) to stage ( 5)). In this model, it is considered that the already damaged element cannot completely recovered, even the crack is partially closed. The degradation of the element is permanent and irreversible.

The re-opening criterion is calculated with the lowest residual tensile strength and the residual fracture energy

 re =  op - ′ op .
The re-opening criterion writes as:

Φ o = 𝜎 eq -𝜎 ′ y ( 1 -exp ( - 𝜎 ′ y  re [𝑢]
))

, with

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝜎 ′ y = 𝜎 y exp ( - 𝜎 y ( [𝑢] max -[𝑢] cl )  op )  re =  op exp ( - 𝜎 y ( [𝑢] max -[𝑢] cl )  op
) .

(10)

Resolution scheme of the Finite Element framework with double enhancements

This section aims at presenting the resolution strategy of the E-FEM model. The main difficulty is focused on the treatment of the kinematic enhancements, especially the treatment of the unbounded part, such as the Dirac-delta distribution in the strain field that results from the strong discontinuity. In this paper, both strong and weak disconti- (incompatible) kinematics is solved in a compatible form by applying the Statically and Kinematically Optimal Nonsymmetric Formulation (SKON) [START_REF] Dvorkin | Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions[END_REF][START_REF] Simo | A new approach to the analysis and simulation of strain softening in solids[END_REF]Oliver, 1996a). It leads to a non-invasive form for the E-FEM model, which makes it possible to solve the system in a standard and classical Finite Element framework. Moreover, it allows the model to simulate multi-branching cracks, and the size of the solving system stays constant during the appearance of the new cracks.

Since the kinematic description of the weak and strong discontinuities has been documented many times in the literature [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF][START_REF] Yue | Fe modeling of concrete with strong discontinuities for 3d shear fractures and comparison with experimental results[END_REF][START_REF] Sun | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF], the author choose to keep the writing as light as possible for the sake of reading. This section contains two parts. We show in the first part the discretization of the system. Then the resolution of the non-linear equation at the local scale is presented in the second part.

Discretization of the system

Depending on the element is in the opening procedure or the closing procedure, the linearized element can be written in a matrix form.

• If the element is opening:

⎡ ⎢ ⎢ ⎣ 𝑲 bb 𝑲 bw 𝑲 bs 𝑲 wb 𝑲 ww 𝑲 ws 𝑲 s * b 𝑲 s * w 𝑲 s * s + 𝐾 qo ⎤ ⎥ ⎥ ⎦ (𝑘) 𝑛+1 ⎧ ⎪ ⎨ ⎪ ⎩ Δ𝒅 Δ[|𝜀|] Δ[𝑢] ⎫ ⎪ ⎬ ⎪ ⎭ (𝑘+1) 𝑛+1 = ⎧ ⎪ ⎨ ⎪ ⎩ -(𝒇 𝑒 int -𝒇 𝑒 ext ) - 𝑒 [|𝜺|] -Φ o ⎫ ⎪ ⎬ ⎪ ⎭ (𝑘) 𝑛+1 , (11)
• If the element is closing:

⎡ ⎢ ⎢ ⎣ 𝑲 bb 𝑲 bw 𝑲 bs 𝑲 wb 𝑲 ww 𝑲 ws -𝑲 s * b -𝑲 s * w -𝑲 s * s + 𝐾 qc ⎤ ⎥ ⎥ ⎦ (𝑘) 𝑛+1 ⎧ ⎪ ⎨ ⎪ ⎩ Δ𝒅 Δ[|𝜀|] Δ[𝑢] ⎫ ⎪ ⎬ ⎪ ⎭ (𝑘+1) 𝑛+1 = ⎧ ⎪ ⎨ ⎪ ⎩ -(𝒇 𝑒 int -𝒇 𝑒 ext ) - 𝑒 [|𝜺|] -Φ c ⎫ ⎪ ⎬ ⎪ ⎭ (𝑘) 𝑛+1 , ( 12
)
where the subscript 𝑛 presents the time step, and the superscript (𝑘) denotes the resolution iteration step. The explicit expression of each matrix is given in the appendix.

In Eq. 11 and Eq. 12, the three incremental terms {Δ𝒅, Δ the displacement field can be solved at the global level by making a static condensation on the local variables [START_REF] Wilson | The static condensation algorithm[END_REF]. The solving equation of the assembled global system writes as:

𝕂 sc | | | (𝑘) 𝑛+1 Δ𝒅 | | | (𝑘+1) 𝑛+1 = - 𝑛 e 𝔸 𝑒=1 {𝒇 𝑒 int -𝒇 𝑒 ext } | | | (𝑘) 𝑛+1 , with 𝕂 sc | | | (𝑘) 𝑛+1 = 𝑛 e 𝔸 𝑒=1 𝑲 sc | | | (𝑘) 𝑛+1 . ( 13
)
Providing that the local variables are now computed, the modified stiffness matrix 𝑲 sc can be calculated by writing the internal variables in terms of the displacement:

• If the element is opening:

𝑲 sc | | | (𝑘) 𝑛+1 = 𝑲 bb - [ 𝑲 bw 𝑲 bs ] ( [ 𝑲 ww 𝑲 ws 𝑲 s * w 𝑲 s * s + 𝐾 qo ] (𝑘) 𝑛+1 ) -1 [ 𝑲 wb 𝑲 s * b ] (𝑘) 𝑛+1 , (14)
• If the element is closing:

𝑲 sc | | | (𝑘) 𝑛+1 = 𝑲 bb - [ 𝑲 bw 𝑲 bs ] ( [ 𝑲 ww 𝑲 ws -𝑲 s * w -𝑲 s * s + 𝐾 qc ] (𝑘) 𝑛+1 ) -1 [ 𝑲 wb -𝑲 s * b ] (𝑘) 𝑛+1 . ( 15
)
It can be seen from Eq. ( 14) and Eq. ( 15) that the modified stiffness matrix 𝕂 sc has the same dimension as the elastic stiffness matrix 𝕂 bb . It confirms the non-invasive form of the E-FEM model that the enhancements will not increase the global unknowns of the solving system. However, the resolution of the non-linear equations of the strong discontinuity requires increasing computation time. Some analytical solutions are available for this problem which decrease significantly the computation time, they are presented subsequently.

Resolution of the local non-linear equations

It can be seen from Eq. 11 and Eq. 12 that the non-linearity of the system is only located at the resolution of the strong discontinuity, which is related to the non-linear crack opening/closing criterion Eq. 7 and Eq. 8. Being equations in terms of the crack opening value [𝑢], these two criteria can be reformulated into the following form:

• The crack opening criterion [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

𝑇 𝑒 + 𝑀[𝑢] = 𝜎 y exp ( - 𝜎 y  op [𝑢]
) , with

{ 𝑇 𝑒 = ( 𝑲 s * b -𝑲 s * w 𝑲 -1 ww 𝑲 wb ) 𝒅, 𝑀 = ( 𝑲 s * s -𝑲 s * w 𝑲 -1 ww 𝑲 ws ) . ( 16
)
• The crack closing criterion [START_REF] Sun | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF]:

-𝑇 𝑒 -𝑀[𝑢] =  cl [𝑢] max ln ( [𝑢] [𝑢] max ) . ( 17
)
The analytical solutions can be obtained by using the Lambert function 𝑊 0 [START_REF] Corless | On the lambertw function[END_REF]. For the element in the opening procedure, the crack opening value [𝑢] is solved as:

[𝑢] =  op 𝜎 y ⎛ ⎜ ⎜ ⎜ ⎝ 𝑊 0 ⎛ ⎜ ⎜ ⎜ ⎝ 𝜎 2 y exp( 𝜎 y 𝑇 𝑒  op 𝑀 )  op 𝑀 ⎞ ⎟ ⎟ ⎟ ⎠ - 𝜎 y 𝑇 𝑒  op 𝑀 ⎞ ⎟ ⎟ ⎟ ⎠ . ( 18
)
For the element in the closing phase, the analytical solution writes as:

[𝑢] = - cl 𝑊 0 ⎛ ⎜ ⎜ ⎝ - 𝑀[𝑢] 2 max exp( 𝑇 𝑒 [𝑢] max  cl )  cl ⎞ ⎟ ⎟ ⎠ / ( 𝑀[𝑢] max ) . ( 19
)
The weak discontinuity [|𝜺|] can be determined by the solution of [|𝒖|], the expression gives as:

[|𝜀|] = -𝑲 -1 ww ( 𝑲 wb 𝒅 + 𝑲 ws [𝑢] ) . ( 20
)
Therefore, all local non-linear equations are resolved and the internal unknown variables are determined. Then the assembled modified stiffness matrix 𝕂 sc can be completely defined for the resolution of the global system. This model is implemented in a opened FE code FEAP [START_REF] Taylor | FEAP-ein finite element analysis programm[END_REF]. The non-linear global system is solved by using a quasi-Newton BFGS algorithm [START_REF] Matthies | The solution of nonlinear finite element equations[END_REF].

First Author et al.: Preprint submitted to Elsevier

Meso-scale modeling of the FPZ in concrete -Noteched cube

As previously stated, research has demonstrated that the formation of macro-cracks can result from a complex process starting with the scattered micro-cracks in the Fracture Process Zone (FPZ) [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF]. This section focuses on the investigation of this transition process from micro-cracks at the fine-scale to macro-cracks using an E-FEM model. The simulation is based on an experiment performed by Koji Otsuka in 2000 [START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF]. In this article, experiments are carried out with X-ray and Acoustic emission (AE) techniques to investigate the FPZ behavior in concrete. The effect of the aggregate size is also evaluated.

It is important to note that the numerical simulation has several limitations, which are addressed by making certain simplifications:

• The model resolution is static and does not consider the rate of loading.

• The volume fraction and granulometry of the aggregate in the concrete are not provided in the article, so we chose a volume fraction of 30% and an 'S-shaped' curve for the aggregate granulometry to adapt to a general case of concrete.

• The environmental influences such as temperature and humidity cannot be directly considered in the numerical simulation, but their effect is taken into account in the material parameters.

This section is divided into three parts. In the first part, the numerical simulation setup is presented, including the construction of morphological models, the establishment of boundary conditions, and the determination of material parameters. The second part focuses on the fracture evolution in the FPZ. Finally, the impact of aggregate size is discussed at the end of this section.

Numerical simulation setups

This section is focused on constructing numerical models that can be compared to the experiments in [START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF] in terms of dimensions, boundary conditions, morphological structure, and material properties. It is presented in the following parts how we attain this purpose.

Dimensions and boundary conditions of the specimen

The dimensions of the specimen in experiments are given in Fig. 6(a). It is a notched specimen that is made from a cubic plate. In experiments, the specimen is loaded in tension under displacement control through two steel plates bonded to the specimen using an adhesive. To match the experimental setup, the boundary conditions of the numerical simulation are also set in the same manner, as shown in Fig. 6(b).

Construction of the morphological models

Given that the dimensions of the numerical specimen have now been defined, the next step is to construct the morphological heterogeneity model. In this paper, the specimen heterogeneity is taken into account by considering two material phases: matrix and aggregates. The tool used to construct the heterogeneity structure is the sphere packing method, as described in [START_REF] Stamati | spam: Software for practical analysis of materials[END_REF]. The key point is to generate N spheres of different diameters ([d min , d max ]) and then arrange them randomly in the field while ensuring a minimum rejection length between them to avoid overlaps.

In experiments, three types of concretes are investigated, which are distinguished by their maximum aggregate size d max : 10, 15, and 20 millimeters. However, there is no information available on the volume fraction and granularity curves of the aggregates. Hence, in the numerical simulations of this section, a volume fraction of 30% is chosen for all three types of concrete, and an 'S-shaped' curve is assumed for the granulometry curve to approximate a general situation in concrete. The minimum aggregate size for all three types of concrete is selected as 5 millimeters to eliminate any differences between them.

Within a range of sphere diameters [d min , d max ], the size distribution of the spheres is nearly linear. Therefore, to construct the 'S'-shaped granulometry curves, six classes of spheres are used, as shown in Table . 1. As previously stated, the spheres are randomly arranged in the field, which can result in differences between generations, even with the same settings. The size of the statistical sample is not large enough to reduce this numerical discrepancy. Therefore, this section generates six different fields (noted as r1, r2, r3, r4, r5, and r6) for each type of concrete. The measured granulometry curves are displayed in Fig. 7. It can be seen that there is no obvious difference between the six generations by applying the same setting. It is worth noting that in the E-FEM model of this paper, only one weak discontinuity is embedded in the finite element. Therefore, a finite element can only correctly present morphological information if at most one heterogeneity interface passes through the element. To meet this requirement, the size of the finite element must be smaller than the minimum distance between each sphere. However, the small aggregates and the narrowness between spheres will increase the fin finite element mesh requirement, which will significantly boost the computational time and memory.

Therefore, the choice is made to simulate the specimen in two parts: a coarse part and a refined part (see Fig. 8). The projection of the morphological structure to the mesh is shown in Fig. 9. The coarse mesh part, as shown in Fig. 8, has an average finite element size of 20 mm. It is a homogenized zone, which is considered as a merged homogeneous material of the cement paste and the aggregates. The homogeneous material (note as material 1, see Fig. 9) is considered to be purely elastic based on the fact that the most significant crack pattern is mainly observed in the notched part of the specimen.

The refined part, on the other hand, has an average finished element size of 0.9 mm. It is a heterogeneous area that is composed of cement paste and aggregates. After the projection, three types of elements are generated: material 2 represents the cement paste which contains cement and river sands; material 3 represents the aggregates; and material 4 represents the weakly enhanced elements which carries the geometry information about the heterogeneity interface. It can be seen from Fig. 9 that the weakly enhanced element of one aggregate does not overlap another weakly enhanced element of another aggregate. The fracture behaviors of the specimen occur only in the refined region. The crack pattern is one of the most important points in this section.

For further information, the morphological model comprises approximately 870 thousand nodes and 5.4 million elements, requiring approximately 6 GB of computing memory. As previously stated, the differences between genera-First Author et al.: Preprint submitted to Elsevier tions are a result of the randomly generated morphological structures, such as the positions of the aggregates. Therefore, to avoid falling into particular cases, six different fields are generated for the three types of concretes (d max = 10 mm, 15 mm, and 20 mm).

Identification of material parameters

The material parameters at the meso-scale are determined by conducting an inverse analysis, which involves minimizing the discrepancy between the numerical results, which is the average response of three generations, and the experimental results. This process enables the identification of four parameters for each type of element: Young's module 𝐸, Poisson's ratio 𝜈, critical tensile strength 𝜎 𝑦 , and the fracture energy  op . The former two parameters are elastic parameters and can be defined by comparing numerical and experimental results in the elastic phase. The latter two parameters are critical parameters, which are related to the failure behaviors of the concrete. The load curve in terms of the crack opening displacement (COD) is plotted in Fig. 10, in which, the maximum aggregate size of concrete is 10 mm. The parameter identification process is carried out in three steps.

In the first step, it is assumed that all four types of elements are elastic and have the same parameters as material 1, the homogenized material of the cement paste and aggregates. The purpose of this step is to determine the elastic parameters of material 1, as it is considered to be purely elastic and therefore has no critical parameters to be identified.

In the second step, the elastic parameters of the materials in the refined zone are determined. It is assumed that materials 2, 3, and 4 are purely elastic. Firstly, the Young's modulus and the Poisson's coefficient of the aggregates are selected. Based on the research of Makani [START_REF] Makani | Analytical estimate of the mechanical behavior of rock: Granitic aggregates[END_REF], the Young's modulus of the calcium type aggregates of concrete (material 3) is chosen to be 70 GPa, and the Poisson's coefficient is selected to be 0.2. Then, the elastic parameters of the matrix (material 2) can be identified if the numerical simulation reproduces the same macroscopic elastic modulus as the experimental results.

In the final step, the critical parameters of the materials in the refined zone are identified. These parameters, related to the failure behaviors of the concrete at the macro-scale, can be identified if the average numerical simulation produces the same maximum macroscopic resistance and softening curves as the experimental results, see Fig. 10. In this step, the tensile strength of the aggregate is selected based on Makani's research [START_REF] Makani | Analytical estimate of the mechanical behavior of rock: Granitic aggregates[END_REF], which found that the compressive resistance of calcium aggregates is around 224 MPa. As it is widely known that for geomaterials, the compressive resistance is normally more than 10 times higher than the tensile strength. Therefore, we set the tensile strength of the aggregate at the meso-scale to be 10 MPa.

Considering the two previous steps, the random distribution of aggregates must be taken into account. In this study, we performed these steps on the six independent generations of the d max = 10 mm concrete to obtain an average. Then It is worth noting that our interests are primarily focused on the development and evolution of cracks in the FPZ.

To this end, we have decided not to apply the load until the total failure of the material, but to refine the loading steps at the failure initiation stage. In this way, the evolution of the crack pattern in the FPZ at different loading levels can be well captured. The identified parameters for the four types of elements are summaries in Table . 2.

Crack pattern -formation of the Fracture Process Zone

Given the identification of the simulation configurations in the previous section, this part focuses on the simulation of the crack pattern in the fracture process zone (FPZ). To illustrate the progression and formation of the FPZ, a concrete sample with d max = 10 mm inclusions is selected as an example. Five steps after the maximum resistance are chosen, see Fig. 11. Regarding the opening elements:

• Firstly, we can see from the figure that there is a high density of diffuse micro-cracks surrounding and in front of the macro-crack at each loading stage. Obviously, the orientation and location of these scattered micro-cracks are strongly influenced by the underlying heterogeneous meso-structure.

• Secondly, as the loading increases, some of these micro-cracks eventually transform into macro-cracks, whose path is directly inherited from the micro-cracks at the fine-scale.

• Finally, during macro-crack coalescence, new micro-cracks emerge surrounding and in front of the macro-crack.

This transition procedure continues until the complete post-peak failure.

In this work, the transition from micro-cracks to macro-cracks is gradual, and the distinction between the two is based on the size and distribution of cracks. Micro-cracks refer to scattered, separate cracks with relatively small crack opening values, and stresses are still transferred through the interface of the crack. Macro-cracks, on the other hand, are cracks that form the major crack and actually cause the collapse of the material. They are characterized by a relatively large crack opening value, which results in a stress-free crack and crack closures.

In addition to the open elements, the closed elements can also be observed in the figure. As the loading increases, the crack closure behavior takes place mainly in the vicinity and branches of the macro-crackdue to stress release caused by the fully formed macro-crack. Therefore, no closed element can be observed in the area in front of the macro-crack, since there are only spread micro-cracks exist and no macro-crack has formed yet.

The crack features of generations r2 and r3 in the FPZ can be seen in the annex, depicted in Fig. 17 and Fig. 18 respectively. It can be seen that the crack paths of the three generations are all different, due to variations in their mesostructures. However, the pattern of crack evolution for the three generations is consistent. For further information, the simulation result shows several crack features in common with the study of [START_REF] Grassl | Meso-scale modelling of the size effect on the fracture process zone of concrete[END_REF], such as the tortuous crack path, scattered micro-crack, and the crack closures. [START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF].

Size effect in the Fracture Process Zone

In this section, the focus is made on the size effect of the inclusions. According to the experimental results in [START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF], at the peak load, as the maximum inclusion size increases, the width of the Fracture Core Zone (FCZ) increases, while the length of the FCZ decreases, see Fig. 14. Here, the FCZ is defined as the area that carries 70% of the total energy of all AE events, see Fig. 13(a). This size effect can also be observed at the post-peak load, see Fig. 12.

This section contains two parts. The first part of this section examines the relationship between the inclusion size and the dimensions of the FCZ at peak load, and compares the simulation results with the experimental results. The second part focuses on the crack pattern of different d max concrete at post-peak load. A possible explanation of the size effect is given from the numerical point of view.

Effect of the inclusion size at the peak load

In section 5.1.2, three morphological models are constructed for each type of concrete that carrying different d max .

In this section, the same material parameters in Table . 2 are applied to these models to investigate the size effect of the inclusions.

To make a reasonable comparison between the experimental results and numerical results, it is necessary to present the method used to measure the dimension of the Fracture Core Zone (FCZ) in the simulations. In experiments, the FCZ is defined as the area carrying 70% of the total energy, as seen in Fig. 13(a). However, it is difficult to define such an area in numerical simulations. Therefore, in this study, the elements in the FCZ are chosen as those that have a crack opening value [𝑢] greater than 0.001 mm. The dimensions of the FCZ are then measured by embedding all these elements in a cuboid and measuring the dimensions of the cuboid, see Fig. 13(b). It is worth noting that isolated elements that located far from other cracks will not be included in the FCZ to prevent overestimation of its length and width. This approach is used for a qualitative comparison between the simulation and experimental results, since a quantitative comparison may lack rigorous validity. The comparison is illustrated in Fig. 14.

Regarding the width of FCZ in experimental results, as shown in Fig. 14 inclusion size equal to 10 mm, one test for inclusion size equal to 15 mm, and two tests for inclusion size equal to 20 mm. For the same type of concrete, the maximum discrepancy between two tests of the same concrete can be found at the d max = 10 mm, which is equal to 13.0 mm. The average width increases from 28.9 mm to 45.2 mm as the inclusion size varies from 10 mm to 20 mm. In numerical simulations, there are six simulations for each type of concrete. Similarly, discrepancies can be found between each generation for the same type of concrete, with the maximum difference being 12.42 mm for the d max = 15 mm concrete. The average widths of FCZ are 17.9 mm and 30.25 mm for the smallest and largest inclusion size concrete, respectively.

Regarding the width of FCZ, there are two positive correlations that can be seen between the experimental results and simulation results. Firstly, the average width of the FCZ increases as the maximum inclusion size (d max ) increases.

Secondly, the discrepancies between two samples of the same type of concrete are observed.

For the former, the increase in the average width of the FCZ with the increase in the maximum inclusion size (d max ) is due to the heterogeneous structures in the concrete, as the same material parameters were used for all types of concrete. From the numerical point of view, this effect can be explained by the fact that the cracks are more tortuous as the inclusion size increases, since the cracks tend to deviate through inclusions to avoid obstruction. Readers can refer to the section 5.3.2 for more details. For the latter, clearly, the discrepancy between two samples of the same type of concrete in both experiments and simulations can be explained by the statistical random distribution of aggregates. Now let us move our attention to the length of FCZ, see Fig. 14(b).

The experimental results of the length of FCZ show that the discrepancies between different samples of the same type of concrete are substantially larger than those in the width of FCZ. The largest difference between two tests is 21.3 mm for the d max = 10 mm concrete and 32.4 mm for the d max = 20 mm concrete. For the d max = 20 mm concrete, the longest FCZ measures 47.3 mm which is three times longer than the shortest FCZ of 14.9 mm. The average length decreases from 40.0 mm to 27.2 mm as the inclusion size increases from 10 mm to 20 mm. As for the simulation results, it can be seen that the average length is larger than it is in experimental results. The discrepancy between tests of the same type of concrete is also smaller. There are two possible reasons for these differences. On the one hand, the number of numerical sample (six samples for each kind of concrete) is insufficient to achieve a statistically meaningful conclusion. On the other hand, the length of FCZ in numerical simulations may be affected by the choose of our morphological model. It is presented in Section. 5.1.2 that the numerical sample is divided into two parts, an elastic part and a heterogeneous part. The failure behaviors can only appear at the latter part. Limited by the computation power, the width of this region is chosen to be 50 mm. It is can be noticed in Fig. 11 that the farthest micro-cracks can be observed up to the border of the refined zone, which indicates that the refined zone is obviously insufficient to capture all fracture behaviors in the FPZ. This may lead to a longer FCZ since there is no energy being dissipated outside the 50 mm refined zone, which is not consistent to the experimental observation in Fig. 15. Moreover, This effect is more pronounced for concrete with larger aggregates, as the FCZ of large d max concrete is shorter and wider, and the wider distribution of micro-cracks means reducing the refined zone has a greater impact on the behavior of the FCZ.

Effect of the inclusion size at the post-peak load

In this part, we focus on the crack pattern of the different d max concretes at the post-peak load. The crack pattern of three types of concrete is plotted in Fig. 16. It can be seen that they share many patterns in common, such as the diffused micro-cracks besides and in front of the macro-crack, branching besides the macro-crack, and the closed elements at the branches and neighboring parts of the macro-crack.

Most importantly, it can also be seen from the figure that as the d max increases, the macro-crack becomes more circuitous. It is due to the fact that in case of being obstructed by an inclusion, the crack tend to deviate to continue propagating. Therefore, larger inclusions lead to a more tortuous crack pattern and underlying more dissipated energy in the crack path. Hence, larger inclusion size would lead to a wider and shorter FCZ. It could be a possible explanation for the size effect from the numerical point of view.

Discussions and conclusions

This paper uses a meso-scale approach to model the Fracture Process Zone (FPZ) in concrete. Several notched cube specimens with statistical random distributed aggregates were performed to investigate the formation of the FPZ and its relationship with the inclusion size, namely, the size effect.

At the meso-scale, the specimen is considered as a highly heterogeneous material, and the evolution crack pattern from the fine-scale spread micro-cracks to macro-cracks can be observed. Here, the aggregates are performed explicitly by the sphere packing method [START_REF] Stamati | spam: Software for practical analysis of materials[END_REF] and the embedded weak discontinuities [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF][START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF]. The cracks are depicted by the strong discontinuities [START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids[END_REF][START_REF] Wells | Discontinuous modelling of strain localisation and failure[END_REF]. Providing that the mode-I crack mode is fairly a common choice for the modeling of FPZ, a crack opening-closing E-FEM model [START_REF] Sun | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF] is applied in this paper. By combining the weak/strong enhancements, the presented E-FEM model is capable of simulating complex crack pattern evolution and considering the interaction between cracks and aggregates.

The numerical simulations are compared with relevant experiments that is realized by [START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF].

The comparison results in several conclusions:

• First, the failure process in the FPZ is a complex procedure starting from the diffused local micro-cracks whose location and orientation are closely related to the heterogeneity structure. The coalescence of micro-cracks forms macro-cracks, while new micro-cracks appear in the front of the macro-crack. Meanwhile, crack closures can be observed at the stress release zone, such as at the branches of the major macro-crack.

• Second, we managed to simulate the experimentally observed size effect, that as the inclusion size increases, the FCZ becomes wider and shorter. The most likely explanation for this observation is that a larger inclusion size leads to more tortuous cracks. It is because the crack tends to deviate through inclusions to avoid obstruction.

Thus more circuitous cracks would lead to more energy being dissipated in the crack path.

However, being limited by the computation power, the number of numerical samples is insufficient to achieve a statistically meaningful conclusion. Because of the same limitation, the non-linear zone (zone carrying refined meshes) is chosen to have a width of 50 mm, which is not enough to capture all local-scale cracks in the notched cube specimen.

It could be a possible reason that addresses the difference between the simulation and experiments in the comparison of the FCZ length.

Appendix A -Expressions of the matrix

The explicit expressions of the matrix in Eq. 11 and Eq. 12 are listed in this appendix. The definition of each term in the expressions is given as:

𝑲

  ORCID(s):First Author et al.: Preprint submitted to Elsevier

Figure 1 :

 1 Figure 1: Concept of the crack evolution in the FPZ for heterogeneous media.

Figure 2 :

 2 Figure 2: illustration of the solid that exhibiting heterogeneities and cracks, with their corresponding elements.

Figure 3 :

 3 Figure 3: Three elements carrying different types of enhancements.

First

  Author et al.: Preprint submitted to Elsevier

Figure 4 :

 4 Figure 4: Behavior of the strongly enhanced element at the crack opening phase.

  (b)) equals to the consumed energy in opening (dark gray part in Fig. 5(b)):

  nuities [|𝒖|] and [|𝜺|] are solved inside the element as additional unknown variables at the local scale. The unbounded First Author et al.: Preprint submitted to Elsevier

Figure 5 :

 5 Figure 5: Mechanical behavior of the local constitutive models at the discontinuity interface.

  [|𝜀|], Δ[𝑢]} are solved in two different levels. First, the weak and strong discontinuity (Δ[|𝜺|] 𝑘+1 𝑛+1 , Δ[|𝒖|] 𝑘+1 𝑛+1 ) are solved at the element level as internal variables for a fixed displacement field 𝒅 𝑘 𝑛+1 . It would lead to a null residual of  𝑒 [|𝜺|] and balanced opening/closing criterion Φ o ∕Φ c . Then

Figure 6 :

 6 Figure 6: Geometry information of the notched cube specimen in experiments[START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF] and in simulations.

Figure 7 :

 7 Figure 7: Measured granulometry curves for each type of concrete.

Figure 8 :

 8 Figure 8: Basic geometry information about the mesh with a zoomed view at the refined zone.

Figure 9 :

 9 Figure 9: Projection of the defined field with random morphology to an unstructured mesh.

Figure 10 :

 10 Figure 10: Experiment (Otsuka and Date, 2000) and simulation load curve of d max = 10 mm concrete in terms of the crack opening displacement (COD).

First

  Author et al.: Preprint submitted to Elsevier the average simulation results are used to be compared with the experimental results. The values of material parameters are iteratively adjusted until the simulation results matched the experimental data as closely as possible.

Fig. 11

 11 Fig. 11 presents the front view of our 3D simulation. The aggregates are depicted in a transparent gray color and the cracks are divided into two categories: opening elements and closed elements. The opening elements are displayed in blue and represent elements with continuously increasing crack opening values. The closing elements are displayed in red and correspond to elements with decreasing crack opening values at some point. These two groups of elements are mutually exclusive.

Figure 11 :

 11 Figure 11: Front view of the crack pattern evolution of d max = 10 mm r1 concrete: crack openings and crack closures (aggregates marked in transparent gray).

Figure 12 :

 12 Figure 12: Comparison of FCZ traced from X-ray films of concrete with four different d max[START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF].

  (a), there are four experimental tests for First Author et al.: Preprint submitted to Elsevier

Figure 13 :Figure 14 :

 1314 Figure13: Measurement of the FCZ in experiments[START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF] and in simulations.

Figure 15 :

 15 Figure 15: Relationship between d max and dimensions of FCZ at peak load in numerical simulations.

  (a) Concrete d max = 10 mm (b) Concrete d max = 15 mm (c) Concrete d max = 20 mm

Figure 16 :

 16 Figure16: Relationships between Dmax and dimension of the FCZ at peak load[START_REF] Otsuka | Fracture process zone in concrete tension specimen[END_REF].

•

  𝑩: standard strain interpolation matrix • 𝑉 , 𝑉 + ∕𝑉 -: volume of the element, volume of the upper/lower part of the element First Author et al.: Preprint submitted to Elsevier

Figure 17 :

 17 Figure 17: Crack pattern evolution of d max = 10 mm r2 concrete: crack openings and crack closures.

Figure 18 :

 18 Figure 18: Crack pattern evolution of d max = 10 mm r3 concrete: crack openings and crack closures.

  where the displacement field 𝒖 is decomposed into a continuous part û and a discontinuous part (  𝑢 -𝜑 𝑒 )[|𝒖|]. In this equation,   𝑢 is the Heaviside function that centered at the discontinuity interface, and 𝜑 𝑒 is a continuous arbitrary function which has null value in Ω - 𝑒 , and unit value at Ω + 𝑒 . The magnitude of the jump in the displacement field[|𝒖|]

Table 1

 1 Diameter ranging ([mm]) for three different d max concretes. Six classes are used to reproduce the "S" shaped granulometry curve.

	Class	d max = 10 mm	d max = 15 mm	d max = 20 mm	Relative percentage
	1	[5.0, 5.83]	[5.0,6.67 ]	[5.0,7.5 ]	4%
	2	[5.83, 6.67]	[6.67, 8.33]	[7.5, 10.0]	17%
	3	[6.67, 7.5]	[8.33, 10.0]	[10.0, 12.5]	29%
	4	[7.5, 8.33]	[10.0, 11.67]	[12.5, 15.0]	29%
	5	[8.33, 9.17]	[11.67, 13.33]	[15.0, 17.5]	17%
	6	[9.17, 10.0]	[13.33, 15.0]	[17.5, 20.0]	4%

Table 2

 2 Identified parameters of the four types of elements.

	Material	𝐸 [GPa]	𝜈 [-]	𝜎 𝑦 [MPa]	 op [𝐽 ∕𝑚 2 ]
	material 1	6.0	0.2	-	-
	material 2	2.2	0.2	2.8	1.2
	material 3	70.0	0.2	10.0	1.2
	material 4	-	-	2.8	1.0

  bb = 𝑩 𝑇 ( 𝑉 + 𝑪 + + 𝑉 -𝑪 -) 𝑩 𝑲 bw = 𝑉 + 𝑉 - 𝑉 𝑩 𝑇 (𝑪 + -𝑪 -)𝑯 w 𝑲 bs = 𝑩 𝑇 (𝑉 + 𝑪 + + 𝑉 -𝑪 -)𝐆 s 𝒏 𝑪 + + 𝑉 + 𝑪 -)𝑯 w 𝑲 ws = 𝑉 + 𝑉 - 𝑉 𝑯 w (𝑪 + -𝑪 -)𝐆 s 𝒏. 𝑪 + + 𝑉 -𝑪 -)𝑩 𝑪 + + 𝑉 -𝑪 -)𝐆 s 𝒏

	𝑲 wb =	𝑉 + 𝑉 -𝑉	𝑯 w	𝑇 (𝑪 + -𝑪 -)𝑩
	𝑲 ww = 𝑇 (𝑉 -𝑲 s * b = 𝑉 + 𝑉 -𝑉 𝑯 w 𝜕𝜎 eq 𝜕𝑻 1 𝑉 𝑯 * ,𝑇 s (𝑉 + 𝑲 s * w = 𝜕𝜎 eq 𝜕𝑻 𝑉 + 𝑉 -𝑉 𝑯 * ,𝑇 s (𝑪 + -𝑪 -)𝑯 w
	𝑲 s * s = s (𝑉 + 𝐾 qo = 𝜕𝜎 eq 𝜕𝑻 1 𝑉 𝑯 * ,𝑇 𝜎 2 𝑦  op 𝑒 -𝜎 y [𝑢]∕ op
	𝐾 qc =	[𝑢]	 cl

max

[𝑢]

Notwithstanding these limitations, the study certainly adds to our understanding of the crack evolution inside the FPZ and its relationship with the material composite structure. Further research in this field would be of great help in increasing calculating efficiency and reducing required computation power.


Torridon M Millan

C Pozarycki

A Mcadam

S An- Drejkovicova

P Mahaffy

D Glavin

A Buch

C Szopa

C Freissinet

A Srivastava

S Teinturier

C Malespin

R Williams

A Williams

J Eigenbrode

R Navarro-Gonzalez

S S Johnson

OPTIMIZATION OF THE SAMPLE ANALYSIS AT MARS WET CHEMISTRY EXPERIMENT FOR THE DETECTION OF ORGANICS IN GLEN

Introduction:

As Curiosity ascends Mount Sharp, it is analyzing stratigraphic rock layers to find clues to Mars' environmental history and habitability [1]. The Sample Analysis at Mars (SAM) instrument suite onboard Curiosity includes a pyrolyzer coupled to a gas chromatograph-mass spectrometer (pyro-GCMS) mostly dedicated to the search for organic molecules on Mars. SAM is able to perform in situ molecular analysis of gases evolved from heating solid samples collected by Curiosity up to ~900°C. SAM can then detect, separate, and identify volatiles inorganic and organic compounds released from the gas phase of the solid samples.

SAM also carries nine sealed wet chemistry cups to allow for a new capability to be employed at Mars' surface, opening the possibility for a larger set of organics to be detected. Seven of the cups contain a 0.5 mL mixture of N-methyl-N-(tert-butyldimethylsilyl) trifluoroacetamide and dimethylformamide (MTBSTFA:DMF 4:1) for derivatization and the two others contain tetramethylammonium hydroxide (TMAH), 25% in methanol mixture for thermochemolysis [2]. The cups also contain standards of organic molecules (36.2 nmol of 3-fluorovaline and 24.2 nmol of pyrene in the MTBSTFA cups) (Figure 1) for calibration and validation purposes.

The derivatization with MTBSTFA allows the detection and identification of complex, polar and/or refractory molecules. In this silylation method, the labile hydrogen of the targeted molecule is replaced by the MTBSTFA silyl group, transforming molecules into volatile derivatives easily amenable to GCMS analysis (Figure 1). Thus MTBSTFA reacts with organics of astrobiological interest-such as nucleobases, amino acids, carboxylic acids, and sugars-while preventing degradation of their chemical structure. Derivatize benzoic acid and ammonia derivatized were also detected by both their mass spectra and retention times, above instrument background. However, no amino acids or fatty acids were detected.

Optimization of the GCMS parameters: Following the wet chemistry experiment on the OG sample, laboratory experiments were conducted to improve the detection of derivatized organic molecules with SAM in future runs. Tests included the analysis of multiple organic mixtures (e.g. amino acids, carboxylic acids, etc.) under SAM-like operating conditions (e.g. temperature, helium flow) and varying the operating temperature of the GC columns. These runs showed that 1897.pdf 51st Lunar and Planetary Science Conference (2020) by measuring the signal of compounds that elute while the column is cooling down, we can extend the effective range of heavy molecules that can be detected in comparison to previous SAM measurements (Figure 2). Laboratory experiments have also now shown that successive column cleanups after the sample GCMS experiment can release heavier compounds that may not have been released in the original experiment.

Wet chemistry on glycine-clays laced analogs:

To prepare the next chemistry run in Glen Torridon and improve the detection of organics, laboratory experiments were performed on glycine-laced montmorrilonite analogs (Figure 3). The goals were to evaluate the extraction and derivatization efficiency of amino acids, here glycine, present in the phyllosilicates interlayers and evaluate the limitations of the SAM-flight experiment in order to optimize it. Pyrolysis parameters were varied (oven, trap and transfer line temperatures) and the glycine abundance recovery was calculated. Results showed that glycine detection is limited in SAMlike pyrolysis but can be greatly improved by slightly increasing the transfer line and pyrolysis rate temperatures within the range allowed by SAM.

Figure 1 .

 1 Figure 1. The MTBSTFA derivatization cup (top) and an example of the MTBSTFA reaction with an amino acid that displaces the labile hydrogens with a tertbutyldimethylsilyl group to produce a volatile amino acid derivative (bottom).

1

  st wet chemistry on Ogunquit Beach: EGA and GCMS of the Ogunquit Beach scooped sample have been analyzed and showed the detection of compounds of medium to high molecular masses, including derivatized molecules as well as varying amounts of various MTBSTFA by-products, such as N-methyl-2,2,2trifluoroacetamide (TFMA) and 1,3-bis(1,1dimethylethyl)-1,1,3,3 tetramethyldisiloxanetetramethyldisiloxane (bi-silylated water), along with the whole MTBSTFA molecule. These results indicate that the derivatization reagent was in excess, as reactions clearly occurred and produced derivatized compounds.

Figure 2 .

 2 Figure 2. Gas chromatogram obtained with a sp are GC4 column of a 22 amino acids standards mixture analyzed in the laboratory. Only alanine and glycine would be detectable within the current SAM-flight operating conditions used for the GC4 column on the SAM flight model instrument. However, this laboratory experiments shows that recording the signal while the temperature is cooling down on Mars could allow the detection of heavier molecules such as the amino acids targeted for this study.

Figure 3 .

 3 Ground Clay, packed between quartz wool in a quartz pyrolysis tube and spiked with MTBSTFA derivatization reagent before pyrolysis. Pyrolysis & 2 nd wet chemistry on Glen Etive (Glen Torridon): Both a pyrolysis only and wet chemistry experiments were performed on two separate Glen Etive drilled samples from the Clay-Bearing Unit (Glen Torridon). EGA and GCMS have been analyzed and showed the efficiency of the lessons learned from OG through a much less saturated chromatogram in chemistry reagents. High molecular weight compounds are still detected with some abundance higher than those present in the OG sample although less reagent was introduced. The highest abundance of Sulfurbearing organics ever measured by the SAM instrument was detected in both the pyrolysis and wet chemistry runs. Pathways of formation of those sulfur compounds are under investigations through laboratory experiments to understand their potential source and origins.


Longxiang Huang
email: huangl@coria.fr

Diego Ferrando

Chetankumar S Vegad

François-Xavier Demoulin

Benjamin Duret

Julien Reveillon

Curvature-surface analysis applied to a commercial and an academic injection system

Keywords: Curvature, Spray, Primary break-up, Simplex-swirl-atomizer, Drop-Size-Distribution

 proposed an analysis based on the curvature and surface area of the early atomization process to calculate the final DSD. This analysis was further extended by Ferrando [2] adding the capability to calculate the velocity of the droplets. The curvature-surface analysis has been applied to a commercial simplex swirl atomizer [3] using Large-Eddy Simulation (LES) where the DSD was well captured in comparison with experiments [4]. However, the numerical simulation overestimated the joint velocity distribution of the droplets. This particular injector has a complex geometry and several physical processes are involved. The model and discretization schemes chosen play a role in solving the flow, therefore, the overestimation of the droplet velocity may be due to the flow downstream, and not the curvature analysis itself. To further study this analysis, this work uses a simpler academic case with a simpler geometry. A Direct Numerical Simulation is performed using the ARCHER code, in order to avoid using any turbulence model that may affect the final results. Unlike previous works which used the interFoam solver (volume of fluid-based solver), this study uses a more precise ICM called the level-set method. The curvature-surface analysis method is applied to evaluate its effectiveness.

Introduction

In traditional practice, the measurement of spray characteristics, such as droplet diameter or droplet velocity, is performed at a considerable distance from the nozzle orifice to ensure the completion of the break-up process and the formation of fully spherical liquid structures. Experimental methodologies, such as phase Doppler anemometry (PDA), exhibit an enhanced validation ratio as they are implemented further away from the nozzle orifice [START_REF] Vegad | 31th Conference on Liquid Atomization and Spray Systems[END_REF]. Similar limitations are encountered in computational fluid dynamics (CFD) simulations when counting the number of droplets or converting them from an Eulerian formalism to a Lagrangian formalism [START_REF] Warncke | [END_REF]7]. In such cases, droplets are required to possess a spherical morphology, which requires enlarging the numerical domain, and thus, increasing the simulation costs. In recent years, a novel technique for investigating the near flow field to the nozzle has been proposed by [1] through a multiphase simulation using a Volume of Fluid (VOF) approach. The principal objective of this study was to scrutinize the curvature of liquid structures created near the nozzle to estimate the diameter distribution of the droplets. The investigation was carried out in an academic airblast atomizer and validated by comparing it with experimental results [START_REF] Gepperth | 23rd Annual Conference on Liquid Atomization and Spray Systems[END_REF]. Subsequently, some studies have adopted this new approach to obtain spray characterization. For instance, [START_REF] Ferrando | 15th Triennial International Conference on Liquid Atomization and Spray Systems[END_REF] retrieved the velocity of the droplets, together with the Droplet Size Distribution (DSD), in another airblast atomizer. The latest advancements in this field are concentrated on adapting the analysis to the experimental domain [START_REF] Vegad | 31th Conference on Liquid Atomization and Spray Systems[END_REF]. A significant disparity between simulations and experiments arises due to the fact that, in the former, the threedimensional flow field is well-known in space and time. However, as stated in subsequent sections of this article, the data obtained from experiments only provides two-dimensional projections of liquid structures. Consequently, the analysis must be adapted to operate with 2D data rather than 3D data. The objective of this study is to explain the mechanism of the curvature-surface analysis, its utilization in Computational Fluid Dynamics (CFD) and experiments, and to conduct a review of the validations that have been performed on a specific injector.

Curvature-surface analysis

As proposed by Palanti [1], it is plausible to derive the Droplet Size Distribution (DSD) of a spray that has completed the first break-up stage by calculating the curvature at the surface interface between gas and liquid. The distribution of curvatures at the surfaces of liquid structures, weighted by the amount of surface interface, is then used to compute the drop size distribution [START_REF] Canu | [END_REF]. Initially, studies were carried out employing a Volume of Fluid-based solver (interFoam [10]). Consequently, the methodology was developed within the VOF framework. However, to enable the application of this analysis beyond this specific framework, several parameters needed to be defined. The Phase Indicator (PI) quantifies the amount of liquid present in a discretized space. It is analogous to the Liquid Volume Fraction (LVF) or α in VOF solvers. Thus, within the VOF framework, it directly corresponds to the LVF. However, in a 2D projection (Projected Phase Indicator or PPI), it represents the ratio of liquid to gas within a discretized area. Other variables can be obtained from this parameter. For the sake of simplicity, PI will be called α in equations from now on. The curvature κ is defined as the inverse of the droplet diameter, κ = 4/D. In this study, the curvature is determined as the divergence of the normal vector to the surface interface, which is given by the following equation:

κ = ∇ • ∇α |∇α| (1)
However, there exist alternative definitions of curvature that can be employed using a framework different from PI. These alternative definitions must be explored and studied further.

The amount of surface interface between gas and liquid can be obtained as:

Σ = |∇α| (2)
In 3D, it corresponds to the ratio of the amount of surface to the control volume. In 2D, it represents the ratio of the amount of perimeter of the droplet to the control surface.

Equivalent Sauter Mean Diameter

To investigate the evolution of the break-up process, an equivalent length to the Sauter Mean Diameter (SMD or D32) is defined as l32. l32 is equivalent to the SMD when applied to a cloud of spherical droplets and provides an estimation when applied to a cloud of non-spherical liquid structures. l32 is defined as follows:

l32 = 6 V l A l = 6 T S (αu • nds)dt T S (Σu • nds)dt (3)
l32 is the ratio of the liquid volume flux to the surface interface flux, both integrated in time and space in a given plane. The factor of 6 establishes the equivalence between the classic definition D32 [START_REF] Lefebvre | Atomization and Sprays[END_REF] and this particular and more general definition.

D32 = i NiD 3 i i NiD 2 i (4)
When the analysis is applied to a 2D projected field, D32 cannot be estimated. Therefore, the characteristic diameter l21 is estimated using the same approach. The formula for l21 is as follows:

l21 = 4 A l P l = 4 T S (αu • nds)dt T S (Σu • nds)dt (5)
In this case, the factor that makes both equations equivalent is 4.

It has been validated that l32 is a good estimation of the final l32 obtained from experiments [1,[START_REF] Ferrando | 15th Triennial International Conference on Liquid Atomization and Spray Systems[END_REF][START_REF] Ferrando | Atomization and Sprays[END_REF] when the primary break-up is completed.

Surface Curvature Distribution (SCD)

The derivation of the drop size distribution (DSD) requires the initial recovery of the curvature distribution, which is subsequently weighted by the surface interface. As the approach proposed by [START_REF] Canu | [END_REF], the DSD can be converted into the SCD. In this analysis the opposite conversion is made: from the curvature distribution weighted by the surface interface, the DSD is obtained.

It is worth noting that the methodology used to recover the curvature and surface interface density data may differ between studies. For numerical simulations, measurements have been taken by evaluating κ and Σ in each numerical cell within a predetermined box (as in the case of Palanti [1]), or in the face cells of a designated plane (as in the case of Ferrando [START_REF] Ferrando | Atomization and Sprays[END_REF]), which were later integrated over time. On the other hand, in experiments, due to the 2D projection of the data and the uncorrelated nature, on time, of snapshots, the methodology employed involved in measuring the data is like the one used in the first case, a predefined box.

Drop Size Distribution (DSD)

To recover DSD from SCD, certain assumptions and estimations are necessary. This analysis assumes that the surface interface recovered corresponds to a spherical droplet, and its diameter is a function of the curvature. Consequently, the number of droplets can be calculated as a function of diameter using the following equation:

D(κ) = 4/κ (6)
The diameter of the spherical droplet is directly proportional to the inverse of the curvature. Specifically, the diameter can be calculated as four times the inverse of the curvature.

n(κ) = ST (κ) πD(κ) (7)
being ST (k) the total surface recovered for a particular value of κ.

The number of droplets can be determined by dividing the total surface interface, as a function of curvature, by the surface area of a spherical droplet with a diameter of D(k). This assumption of spherical droplets corresponding to the surface interface in near-field regions is only an approximation. Figure 1 shows the SCD, DSD, and an iso-contour of α in 3D coloured by the curvature (blue κ < 0, red κ > 0 and white κ = 0). There is some surface interface corresponding to negative curvatures, which corresponds to bubbles or concave liquid structures, rather than spherical droplets. Therefore, the initial correction to this procedure is to eliminate the surface interface of negative curvatures and thus, the "droplets" with negative diameters. Figure 2 shows the SCD, DSD, and the iso-contour of Figure 1 with the initial correction already applied. In this instance, the iso-contour is coloured grey in the regions where the surface corresponds to negative curvatures. The measurement plane is represented by the red line, nearly all surfaces are recovered and not removed. However, as evident from the DSD, there exist a small number of excessively large droplets corresponding to near-zero curvatures. These droplets are unphysical since they correspond to almost flat elements that are not close to spherical droplets. To resolve this issue, the D32 calculated with l32 is imposed in the DSD since it is a reliable value, as previously demonstrated. As a result, the largest droplets are removed from the distribution until the D32 calculated using Equation 4 in the DSD and the l32 calculated with Equation 3 are identical. Figure 3 shows the final outcome of the approximations. The corrections made have been verified against experimental results and empirical equation [1,[START_REF] Ferrando | 15th Triennial International Conference on Liquid Atomization and Spray Systems[END_REF][START_REF] Ferrando | Atomization and Sprays[END_REF].

To apply it in 2D projected cases some adaptations are done:

1. The diameter is 2 times the inverse of the curvature: D(κ) = 2/κ 2. The number of droplets is calculated with the perimeter of one droplet:

n(κ) = PT (κ) πD(κ)
3. The corrections on the DSD are done using the l21 instead of l32.

Simplex Swirl Atomizer Configuration

The present study focuses on investigating the spray formed by the Coria Rouen Spray Burner (CRSB) injector, which is a simplex swirl atomizer. This injector has been the subject of several previous works [START_REF] Verdier | [END_REF][START_REF] Marrero | [END_REF][START_REF] Vegad | 31th Conference on Liquid Atomization and Spray Systems[END_REF], where PDA measurements were conducted to characterize its behaviour.

To obtain the necessary data for the curvature-surface analysis, both the experimental and numerical setups are described in this section.

Experimental set-up

The CORIA Rouen spray burner (CRSB) was adopted to atomize n-heptane fuel. The burner comprises the fuel injection system and plenum. The fuel injector is located at the top of the injection system and surrounded by a co-flow air stream guided through a convergent annular passage. The downstream end of the annular passage is 20 mm outer diameter and 10 mm inner diameter. While its upstream end is connected to the plenum housing 18 radial vanes of zero-degree angle. The purpose of providing radial vanes is to break large flow structures. Further details on burner geometry are available in [START_REF] Verdier | [END_REF]. In Fig. 4, Z and R represent the axial and radial directions, respectively. The spray formation is in the downward (positive Z) direction. An industrial pressurized swirl injector (Danfoss make) of 80 • hollow cone, 1.35 kg h -1 flow rate and 200 µm orifice diameter was used for n-heptane atomization. An internal geometry of the injector is such that an extent of fuel pressure is utilized to form a rotating liquid film leading towards the orifice, while the remaining fuel pressure is applied to force the liquid through an orifice as a conical sheet into the air. N-heptane mass flow rate was 0.28 g s -1 and regulated using Coriolis mass flow controller (Bronkhorst, CORI-FLOW, 0-2 g s -1 ). The air mass flow rate was set to 6 g s -1 and controlled by a thermal mass flow controller (Bronkhorst, EL-FLOW, 0-15 g s -1 ). The temperature of n-heptane and air were measured using K-type thermocouples (T1 and T2), and their pressures were measured using pressure transducers (P1 and P2).

To perform 2p-LIF experiments, a regenerative amplified Ti:Sapphire (Coherent Inc.) laser system of 800 nm wavelength, 3.5 mJ pulse -1 energy, 1 kHz repetition rate, and 100 fs pulse duration was used. The laser beam exhibits a 5 mm full width at half maximum (FWHM) Gaussian profile. A cylindrical lens was employed to form a ∼50 µm thin laser sheet illuminating fuel-tracer solution spray. After absorption, the excited tracer molecules jump to higher energy states and decay to lower states by emitting fluorescent photons. An ICCD camera (PiMAX 4, 1024 × 1024 pixels 2 ) was orthogonally positioned to collect fluoresced light (emission) from the spray. Microscopic imaging was preferred in the present study [START_REF] Durst | [END_REF] as the primary breakup of the n-heptane conical sheet is at about 1 mm axial distance from the nozzle exit. Thus, a long-distance microscope (LDM, QM-100, Questar) was equipped with a camera to capture a 2.5 × 2.5 mm 2 field of view at the nozzle exit.

A direct measurement of drop size was carried out using a phase Doppler anemometry (PDA) system (Dantec make). The purpose of PDA measurement is to compare drop size distribution (DSD) predicted from the curvature analysis of numerical and experimental data. In PDA arrangement, green (514. It was challenging to find a tracer dye for n-heptane fuel in terms of 2p-LIF technique. Common tracers for water or amines are fluorescein, eosin, Rhodamine 6G, and many [15,16]. None of these dyes are soluble in non-polar alkanes, i.e., n-heptane. Further, very few tracer dyes for n-heptane fuel are noted in the existing literature. Adding a small amount of organic compound such as alcohol may help to dissolve polar (i.e., eosin) or slightly polar (i.e., Nile red) tracer into non-polar fuel (i.e., n-heptane). Durst et al. [START_REF] Durst | [END_REF] added 20% propanol to dissolve eosin Y salt (60 mg l -1 ) in n-heptane. Adding a small percentage of ethanol into n-heptane and Nile red solution results in a maxima shift in the emission spectrum [START_REF] Durst | [END_REF]17]. Thus, adding organic compounds can increase the tracer solubility or emission, but it may also change fuel properties. However, this is not the case here. No organic compound was added during the solubility test of various dyes in n-heptane. The solubility test for solvent n-heptane was performed for the fluorescein-548 (F-548), pyrromethene-597 (PM-597), pyrromethene-567 (PM-567), and coumarin-504T (C-504T). All the dyes tested here are in the solid phase and nontoxic. In comparison, coumarin-504T did not dissolve entirely as some solid particles were observed at the bottom of the solution. However, the other two dyes pyrromethene-567 and pyrromethene-597 were seen to be completely soluble in n-heptane.

For each tracer dye, a fuel-tracer solution was prepared with 0.5 % by weight dye concentration in the n-heptane.

The solution was kept in a quartz cuvette and excited with 400 nm and 800 nm laser lights to measure fluoresce characteristics of 1p-LIF and 2p-LIF respectively. The 400 nm laser light was produced by placing a frequencydoubling crystal in an 800 nm laser light path. PM-567 is used as a tracer dye for n-heptane fuel in the present study. The 1p-LIF and 2p-LIF emission spectrum for PM-567 and n-heptane solution is shown in Fig. 5. The measured peak emission (λ max,f l ) for 2p-LIF is ∼576 nm at peak absorption (λ max,abs ) of 800 nm. Similarly, in the case of 1p-LIF, the peak emission ∼571 nm is noted for 400 nm excitation. Both the peaks are marked in Fig. 5. It is to be noted that a part of the excitation wavelengths was directed toward the probe to show the absorption wavelengths in Fig. 5. The theory of two-photon absorption is available in [18,19].

Simulation set-up

Previous works on the Coria Rouen Spray Burner (CRSB) injector have included numerical modelling, with studies conducted by [START_REF] Ferrando | 15th Triennial International Conference on Liquid Atomization and Spray Systems[END_REF][START_REF] Ferrando | 31th Conference on Liquid Atomization and Spray Systems[END_REF][START_REF] Ferrando | Atomization and Sprays[END_REF] presenting detailed analyses of the injector's multiphase modelling. Interested readers are directed to consult these references for further information.

In this article, a brief summary of the numerical simulation performed is provided. Experimental measurement techniques were employed to measure the internal 3D geometry of the injector's metallic pieces. Then, a hybrid tetrahedral-hexahedral mesh is generated. This mesh allows for high-resolution modelling in the zone of interest while using a tetrahedral approach in areas where mesh parameters are insufficient, thereby minimizing numerical instabilities.

The internal multiphase flow of the injector is solved using a Volume of Fluid (VOF) based method within the OpenFOAM toolbox, specifically the interFoam solver. The solved equations are:

∂ ūj ∂xj = 0 (8)
∂ ρūi ∂t + ∂ ∂xi (ρūj ūi) = - ∂ P ∂xi + ∂ ∂xj (τij + τtij) + ρgi + fσi (9)
where ū represents the velocity vector, P the pressure, ρ the density of the cell mixture (defined as ρ = ᾱρ l + (1ᾱ)ρg, τij and τtij the viscous and turbulent stresses including all sub-grid scale effects. Additionally, gi the gravitational acceleration, whereas fσ represents the surface tension forces, calculated as:

fσ = σκ ∂ ᾱ ∂xj ( 10
)
where σ is the surface tension between the liquid and the gas, α is the liquid volume fraction and κ is the curvature of the surface interface.

In order to determine the interface position, an equation for the liquid volume fraction α advection is solved:

∂ ᾱ ∂t + ∂( ᾱūj) ∂xj + ∂(U r j ᾱ(1 -ᾱ)) ∂xj = 0 (11)
In addition, a Large Eddy Simulation (LES) framework was adopted, with the WALE sub-grid scale model retained from the internal simulation. However, the atomization process was not adequately captured due to the low mesh resolution above the nozzle.

To address this, a second simulation was conducted on a smaller but more refined numerical domain beginning at the nozzle orifice. In this simulation, the velocity vector u and the liquid volume fraction α were recovered from the first simulation and mapped onto the refined domain to achieve better qualitative agreement with experimental results.

In the second configuration, which involves external flow simulation, the curvature κ, surface interface density Σ, and liquid volume fraction α are obtained to apply the curvature-surface analysis.

Figure 6 shows a central slice of the external flow simulation domain. On it, the black dashed lines represent the planes where the flux of α and Σ are stored to later calculate l32 with Equation 3. The red line represents the plane where κ and σ are stored to calculate first, the SCD and derive from it the DSD.

Results

In Figure 7, l21 (Equation 5) is computed along the axial distance to the injector using 2p-lif images. Each l21 value is computed on one window of the image for a size of 128x1024 pixels, counting a total number of 15 windows with a step size of 64 pixels for each time we step towards the bottom of the image. In Figure 7 b, the drop size distribution is computed on the last window of 2p-lif images with the same parameters used in while calculating l21. Both, the drop size distribution and l21 computed on the last window match the PDA measurement results. However, this correspondence depends on the treatment of the raw images before applaying the curvature-surface analysis. Figure 8 shows Equation 3 applied to the numerical computation in the planes reported in Figure 6. At z = 0mm, l32 gives a high value not corresponding to a physical D32. This value is obtained at the nozzle where the break-up process has not begun yet and, hence, the calculated value corresponds to a corona-like surface. When the liquid sheet begins breaking-up, more ligaments and smaller liquid structures appear generating more surface interface. l32 is then going down up to a point where it reaches an almost constant value (z ≈ 0.2mm). At this point, most of the break-up process is completed, and the remaining process that is yet to finish will not generate a considerable amount of surface interface that will vary the l32 value. Therefore, from this point to the last measured plane (z = 1.3mm), the approximation to D32 is considerably accurate. The figure shows a small increase tendency of l32 from z = 0.2mm to z = 1.3mm. This effect was noticed in experiments by [START_REF] Marrero | [END_REF]. The value at z = 1.2mm, where the DSD is calculated, is l32 = 28µm which is reasonably similar to PDA measurements D32 = 26µm.

Conlcusions

Curvature-surface analysis has been applied in several configurations and both, in experiments and numerical simulations, giving promising results. The analysis estimates accurately characteristics diameters such as D32 and D21 and it is capable of reconstructing a similar drop size distribution in the early stages of the atomization process. So far, the analysis has been applied over VOF simulations and 2pLIF pictures following a similar approach to VOF. The methodology could be adapted to other formalisms and high-accuracy numerical methods to enhance its capabilities.

Thus, more accurate interface capturing methods, such as Level-Set method (LS), may be used to calculate the curvature and surface interface density. Higher-order discretization schemes and Direct Numerical Simulations (DNS) should be tested in order to reduce any numerical errors to corroborate the great benefits of using this technique to calculate the DSD in a stage where the droplet cloud is not fully formed.

In experiments, the treatment of the raw 2p-lif images has demonstrated an influence on the final results. Guidelines to apply the correct filters and parameters should be proposed in the future to be consistent with the methodology.

Figure 1 .

 1 Figure 1. Left to right: Surface Curvature Distribution, Drop Size Distribution, Iso-Contour (α = 0.5) with the curvature κ coloured on it. Raw results after applying the analysis.

Figure 2 .

 2 Figure 2. Left to right: Surface Curvature Distribution, Drop Size Distribution, Iso-Contour (α = 0.5) with the curvature κ coloured on it. Results after removing negative curvature values.

Figure 3 .

 3 Figure 3. Left to right: Surface Curvature Distribution, Drop Size Distribution, Iso-Contour (α = 0.5) with the curvature κ coloured on it. Results after removing large droplets.

Figure 4 .

 4 Figure 4. The schematic of the experimental setup demonstrating the instruments and optics arrangement for 2p-LIF technique.The numbers 1 and 2 correspond to n-heptane and air, MFC-mass flow controller, LDM-long distance microscope, T-thermocouple, P-pressure transducer, Z-axial direction, and R-radial direction.

  5 nm) and blue (488 nm) argon laser beams (2.2 mm diameter and 38 mm spacing) were utilized for measurements. The focal lengths of the transmission and receiver lenses were 350 mm and 310 mm, respectively. Beam crossing resulted in a probe volume of 0.148 × 0.148 × 3.10 mm 3 . The receiver probe was set at θ = 70 • off-axis angle in front scattering mode. It is close to Brewster's angle for n-heptane (refractive index, n = 1.392). It helps to reduce the effect of trajectory ambiguity while collecting reflected and refracted light by the receiver probe. The acquisition was started from Z = 1 mm downstream of the injector orifice and limited by 30 seconds or 1,20,000 droplets to achieve convergent statistics.

Figure 5 .

 5 Figure 5. 1p-LIF and 2p-LIF emission spectrum of n-heptane fuel and PM-567 tracer dye solution. The gray and black lines correspond to 1p-LIF and 2p-LIF techniques, respectively. The peak emissions of 1p-and 2p-LIF when excited with 400 nm and 800 nm laser lights are ∼570 nm (gray line) and ∼576 nm (black line).

Figure 6 .

 6 Figure 6. The liquid volume fraction contour is presented with dashed black lines indicating the planes where the surface interface and liquid volume fraction data are stored for computing the equivalent Sauter Mean Diameter (SMD). The red line denotes the region where the probes are setup to obtain measurements of the Droplet Size Distribution (DSD)

Figure 7 .

 7 Figure 7. On the left: l 21 in function of the axial coordinate calculated using 2p-lif images. On the right: drop size distribution of the 2p-lif images

Figure 8 .

 8 Figure 8. On the left: l 32 evolution on axial coordinate obtained by numerical simulation. On the right: Drop size distribution comparison between PDA and curvature analysis applied on numerical simulation and 2p-lif imagesIn Figure8the drop size distribution obtained with three different methods is shown. The PDA measurements in grey are the reference value to validate the curvature-surface analysis applied to the studied of this work. The CFD results are in agreement with both, PDA measurements and DSD obtained through 2p-lif images.

Acknowledgements

Authors from Université de Rouen Normandie thankfully acknowledge the computer resources at IDRIS, TGCC and CINES under the allocation A0112B06153 and A0092B06153 made by GENCI (Grand Equipement National de Calcul Intensit) and also for the computing time at CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie) under the scientific project No. 2006011.