Yue Sun

Emmanuel Roubin

Jianfu Shao

Jean-Baptiste Colliat

Strong discontinuity FE analysis for heterogeneous materials: the role of crack closure mechanism

Keywords:

We present in this paper a Finite Element based model, which is devoted to describing the failure mechanics of quasi-brittle materials, such as the unilateral effect and the cyclic behaviors with a low number of cycles. The material will be studied at the meso-scale, and considered as a heterogeneous medium. The model is formulated in a framework of the Enhanced Finite Element Method (E-FEM). As internal enhancements, two kinds of discontinuities are performed. On the one hand, strong discontinuities aim to illustrate cracks and fractures. On the other hand, weak discontinuities are used to describe the heterogeneity. In addition to the initiations and propagations of cracks, the closure of cracks is also taken into account. The proposed model is employed to reproduce the mechanical responses of heterogeneous material by adding the closure mechanism to the mode-I strong discontinuity. We show the ability of the model to simulate some of the main characteristics of such materials, such as the emerged asymmetric traction/compression behavior, the unilateral effect, and the hysteresis phenomenon.

Introduction

The process of modeling the inelastic behavior of heterogeneous materials such as concrete is strongly related to the choice of the observation scale. At the macroscale, providing any Representative Volume Element can be defined, the usual approach relies on FE models based on phenomenological theories such as plasticity theories [START_REF] Menetrey | Triaxial failure criterion for concrete and its generalization[END_REF][START_REF] Hofstetter | Computational mechanics of reinforced concrete structures[END_REF][START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] Grassl | Concrete in compression: a plasticity theory with a novel hardening law[END_REF] or damage models [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure[END_REF][START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF]. Those models are based on macroscopic quantities that are space averages of the local corresponding heterogeneous fields and mainly describe stress-strain relationships. Here, the use of classical smeared crack method usually requires an intrinsic length scale to overcome the difficulties related to constitutive laws with softening (mesh-dependent solution). Because of this feature, many regularisation theories have been developed and put into literature such as the non-local approach (Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: survey of progress[END_REF] and the gradient-enhanced damage model [START_REF] De Borst | Gradient-dependent plasticity: formulation and algorithmic aspects[END_REF][START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF]. For cement-based materials, this point of view has been very fruitful during the last decades, leading to a large range of applications related to the behavior of Civil Engineering structures along several loading paths [START_REF] Reynouard | Comportement mécanique du béton[END_REF]. Today, under the consideration of long-term loadings, the macroscopic observation scale is still the only valid choice to model complex and large structures, such as dams, bridges, or nuclear power plants.

Yet, most of the observations made on those kinds of structures or in laboratory tests show that the origin of the material non-linearities (mainly cracking) can be related to some fine-scale mechanisms. Regarding cracks, it is now well-established that macro-cracks result from a complex process starting with fine-scale spread micro-cracks within the so-called Fracture Process Zone [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF]. Although these micro-cracks are on average isotropic and diffuse, the micro-structure and meso-structure (aggregates, voids, etc.) have a strong influence on the crack pattern at this stage. Hence, with the increase in the macro-stress level, this set of micro-cracks eventually turns to macro-cracks in which the orientation and path directly inherit from the underlying meso-structure. This point has led to some strong efforts during the past few years to propose efficient and reliable tools in order to model concrete at the so-called meso-scale, which here stands for the millimeter scale. At this scale, a large part of the aggregates, as well as the largest voids or porosity, can be represented explicitly. Among the most recent models, the reader can refer to Sukumar, Chopp, [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF]; [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF]. Here, following [START_REF] Hautefeuille | Failure model for heterogeneous structures using structured meshes and accounting for probability aspects[END_REF] and [START_REF] Oliver | From continuum mechanics to fracture mechanics: the strong discontinuity approach[END_REF], we develop a 3D strategy that relies mainly on two ingredients: i) incompatible FE meshes in order to set up a simple and efficient way to model heterogeneous materials with the FEM, ii) the strong discontinuity analysis in order to model strain localization and cracking at the meso-scale [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], and get rid of the mesh-dependence difficulties presented above. The key point is that both of those two ingredients can fit into a single method for implementation, which consists of embedding discontinuities within the elements: i) "weak" or strain discontinuity in order to build incompatible meshes, and ii) "strong" or displacement discontinuity in order to model cracks at the meso-scale. Using the emergence property from the meso-scale to macro-scale, some researchers have reproduced some typical mechanical behaviors of concrete-like materials, such as the asymmetric response in tension/compression, and the progressive loss of stiffness, by using mode-I opening discontinuity [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF], or mode-II sliding discontinuity [START_REF] Hauseux | Fe modelling with strong discontinuities for 3d tensile and shear fractures: application to underground excavation[END_REF]. In recent decades, the use of the X-ray micro-tomography technique makes it possible to build morphological structures based on real material heterogeneity in 3D [START_REF] Hashemi | A tomographic imagery segmentation methodology for three-phase geomaterials based on simultaneous region growing[END_REF][START_REF] Lu | X-ray microtomographic studies of pore structure and permeability in portland cement concrete[END_REF]Stamati, Roubin, Andò and Malecot, 2018a). Capable of predicting the crack pattern, the study of [START_REF] Stamati | Tensile failure of micro-concrete: from mechanical tests to fe meso-model with the help of x-ray tomography[END_REF] shows the important role of the morphological structure in order to simulate failure.

To the author's acknowledge, in existing E-FEM studies, the applied mechanics on the discontinuity at the finescale consider only the initiation and propagation of cracks, and the performance of the model is limited at simple and proportional mechanical loadings. Within this context, this paper is mainly focused on the mechanism of crack closure at fine scale. Our main concerns are obviously related to non-proportional loading paths and we show how adding this simple mechanism at the meso-scale leads to a large improvement in the macro-scale solution, allowing it to exhibit a behavior that is consistent with the so-called unilateral effect for concrete. Moreover, we show that the closure mechanism also leads to some improvements at the macro-scale considering simple loadings as tension and compression. The main effect here is clearly the ability of a large set of meso-cracks to close after the localization of a single macro-crack.

The remainder of this paper is organized as follows. Section 2 recalls the main idea related to the choice of kinematics enhancements. Our interest in this paper is focused on mode-I cracks. After the general formulation of the weak and strong discontinuities, the governing equations describing the local behaviors at the strong discontinuity are presented in section 3. The chosen mechanism for the meso-scale closure of cracks is detailed in this section. Then section 4 sets up the Finite Element approximations of the developed framework as well as the numerical solution. Subsequently, section 5 presents several 3D cases in order to illustrate the main abilities of this model. The numerical application is firstly made on a specimen with a very simple morphological structure, which helps us to investigate clearly the crack pattern evolution and illustrate preliminary the performance of the model. Then the model is applied to a concrete-like specimen with a more complex morphological structure. When different mechanical loadings are applied to the specimen, many typical mechanical behaviors of quasi-brittle materials can be observed, such as the asymmetric responses in traction/compression, the unilateral effect, and the hysteresis phenomenon. Among these typical behaviors, many of them are found closely related to the closure of cracks. Finally, a relevant conclusion is summarised and some discussions are made in section 6 to close the paper.

Kinematics enhancements at the meso-scale

It is presented in this section the basic notations that are employed in this paper. Let us consider the reference domain Ω ⊂ ℝ of a solid exhibiting heterogeneities and cracks, see Fig. 1. It is surrounded by a closed boundary Γ ⊂ ℝ -1 , which can be divided into two disjointed boundaries: the displacement boundary Γ and the traction boundary Γ . Considering the framework of Finite Element Method, the domain can be discretized by a set of elements Ω with Ω = ∪ =1 Ω . Within the solid, two kinds of discontinuities are defined:

• the weak discontinuities  (jump in the strain field), such as the material interfaces

• the strong discontinuities  (jump in the displacement field), such as fracture and cracks By coupling both weak and strong discontinuities on the same geometrical interface, decohesion (or debonding) is represented. As depicted in Fig. 2, three cases can occur to an enhanced element: the element carries a weak discontinuity; the element carries a strong discontinuity; and the element carries both a weak and a strong discontinuity. Finally, it is considered that the kinematic behavior of the strong/weak discontinuity are independent, and the strain field admits an additive form

( ) ∶= ∇ sym ( ) + ̃ ( ) + ̂ ( ), (1)
where

∇ sym ∶= 1 2 [∇(•)+∇ (•)
] denotes the symmetric gradient operator. Notation • / • represents the weak / strong discontinuity (Benkemoun, Hautefeuille, Colliat and Ibrahimbegovic). In Fig. 2,[| |] and [| |] represent the ingredients of the weak and strong discontinuity, respectively. They are presented in the following sections.

Kinematics and governing equations of weak discontinuity

The kinematics of the weak discontinuities [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF][START_REF] Markovic | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior[END_REF] is first presented. As stated above, in this FE context, a weak discontinuity is used within an element split by a material interface (aggregate/mortar matrix). As it presents in Fig. 2(a), an element exhibiting a weak discontinuity carries two kinds of materials, a jump in the strain field is emerged from the different material parameters, which can be defined as

[| |] ∶= + ( ) --( ).
Based on the normal vector of the discontinuity surface , an orthonormal basis ( , , ) can be built.

Considering the condition that the displacement field is continuous across the discontinuity surface, and the derivative of the displacement field has a discontinuity only along the direction perpendicular to the interface but no discontinuity in the direction along the interface, it can be deducted that the jump in the strain field [| |] can be entirely defined by three unknown values {[ ] [ ] [ ] } T [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]. The enhanced displacement field is obtained as [START_REF] Gurtin | The linear theory of elasticity[END_REF][START_REF] Ibrahimbegovic | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior[END_REF]:

̃ ( ) = Θ ⋅ ( -)([ ] + [ ] + [ ] ) 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 of Ω , the explicit expression is given latter. The enhanced strain field can then be obtained by taking the symmetrical gradient of ̃ :

̃ = ∇ sym ( ̃ ) = Θ [ ] ⊗ + [ ] 2 ( ⊗ ) sym + [ ] 2 ( ⊗ ) sym .
(3)

Kinematics and governing equations of strong discontinuity

In the present FE context, the strong discontinuities are used within an element to represent a micro crack (jump in the displacement field). Their kinematics has been well documented in the last decades [START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids[END_REF]Oliver, 1996a;[START_REF] Wells | Discontinuous modelling of strain localisation and failure[END_REF][START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]. Based on previous studies, the displacement field of an element exhibits a strong discontinuity expresses as:

= ̂ + (  -)[| |], (4)
where ̂ is a smooth and continuous displacement field that we can impose standard boundary conditions (Oliver, 1996b), [| |] is a continuous function representing the displacement jump,   is the Heaviside function centered at the discontinuity interface, and is a continuous arbitrary function (Oliver, 1996b) which has null value in Ω -, and unit value at Ω + . 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]:

= ∇ sym = ∇ sym ̂ ⏟ ⏟⏟ regular + (  -)∇ sym ([| |]) -([| |] ⊗ ∇ ) sym ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ bounded enhancement ̂ + Γ ([| |] ⊗ ) ⏟⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞ ⏟ unbounded enhancement ̂ , ( 5
)
where is the Dirac-delta distribution which emerges from the gradient of the Heaviside function. The bounded enhancement ̂ can be simplified by assuming [| |] as a constant function, referred as Kinematically Enhanced Strain (KES) [START_REF] Oliver | On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations[END_REF]. It leads to a null value of ( Γ -)∇ sym ([| |]), and ̂ can be obtained as

-([| |] ⊗ ∇ ) sym .

Admissible discontinuity model with cracks closure at fine-scale

As presented in the previous section, it can be noticed that the strong discontinuity brings into the strain field an unbounded enhancement, which makes it difficult to respect the traction vector continuity at the interface. In order to ensure the consistency of the constitutive model with the existence of the strong discontinuity, the DSDA (Discrete Strong Discontinuity Approach) [START_REF] Oliver | On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations[END_REF][START_REF] Brancherie | Modéles continus et" discrets[END_REF] is used in this study. In this framework, the continuum model exhibiting a strong discontinuity induces an underlying discrete implementation at the discontinuity interface. Therefore, beyond the discontinuity, the bulk volumes Ω ∖ can be described by standard continuum stressstrain relationship; and on the discontinuity surface , the traction vector continuity is ensured, the traction vector is served as a bridge linking the two bulk sub-volumes. In the case of exhibiting weak discontinuity, the constitutive model can be employed in the same way [START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids[END_REF][START_REF] Oliver | On the use of strain-softening models for the simulation of strong discontinuities in solids[END_REF][START_REF] Benkemoun | Failure of heterogeneous materials: 3d meso-scale fe models with embedded discontinuities[END_REF], and the weak and strong enhancement admit an additive form, see Eq. 1. A so-called 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) is used to associate the traction vector with the crack-opening, which depicts the failure and softening behaviors of the element. The focus of this section is mainly made on this relationship.

In existing E-FEM, many enriched discontinuity modes have been proposed to describe various behaviors. In our case, the failure mechanism at the level of the element is driven by a cohesive law with which accounts for the local quasi-brittle behavior of concrete. The localization criterion along the principal direction (tensile) and the mode-I opening criterion are describe in the next two sub sections. Then a new closure criterion is proposed.

Failure criterion and discontinuity orientation

To initiate the failure, a localization criterion is used which compares an equivalent stress eq (defined below) with an ultimate stress y (local material property). It reads

Φ l = eq -y , ( 6
)
where a negative value of Φ l means an elastic behavior while a positive value leads to localization. The equivalent stress is defined as the projection of the traction vector onto : eq = ⋅ , thus it depends on the orientation of the surface. The hypothesis is made here that in case of weak discontinuities (i.e. the element contains a material interface), the surface of the morphology defines the interface orientation . It can be noted that in this case, the orientation is known a priori [START_REF] Benkemoun | Contribution aux approches multi-échelles séquencées pour la modélisation numérique des matériaux à matrice cimentaire[END_REF][START_REF] Benkemoun | Failure of heterogeneous materials: 3d meso-scale fe models with embedded discontinuities[END_REF]. In the other case without weak discontinuities (i.e. the element is either within the matrix or an inclusion), the orientation is the principal direction (eigenvector corresponding to the highest principal stress). It can be noted that, in this case, the equivalent stress corresponds to the principal stress.

In opposition to rotating cracks models (Jirásek and Zimmermann, 1998a,b), once the localization is reached, the orientation of the interface is fixed throughout the computation. It is justified by the local aspect of the strong discontinuities in this model that have to been seen are micro cracks.

Traction separation law -Opening criterion

The opening criterion is brought into use after the failure occurs. As it is mentioned in the previous section, the discontinuity surface is considered flat, thus we have [ ] equals to the normal value of [| |], refers to as the opening value. The equivalent stress and the opening value of the crack are linked through a classical softening criterion at the discontinuity interface:

Φ o = eq -y exp - y  op [ ] ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ o . ( 7
)
A positive value of the opening criterion Φ o means an non-equilibrium state of the element and the crack needs to go further until it reaches an equilibrium state, (i.e. Φ o equals to zero). A negative value leads to an elastic loading, unloading, or closure of cracks, which are discussed later. The softening function in Eq. ( 7) is an exponentially decreasing function in terms of the opening value [ ], tending to zero as the crack opens. The fracture energy  op is considered as a local material parameter, which represents the necessary energy to create a fully opened crack. Physically, a higher value of the fracture energy represents a more ductile material. The softening function is non-linear and continuous, an analytical solution is available for solving the function, readers can refer to the section 4.4 for more details.

Fig. 3 shows an example of the local constitutive model on a one-dimensional element. Naturally, as the imposed displacement Δ increases, the crack continues to open at the opening phase. The mechanical behavior at the discontinuity respects the equation Eq. 7 while the bulk volumes Ω ∖ always remain elastic.

Traction separation law -Closing criterion

As depicted in Fig. 4(a), without considerations of crack closures, energy cannot be dissipated in case of local unloading since the crack opening does not decrease. In order to be able to represent non monotonic behaviors (such as cyclic compression tests), a closure law that follows the same context as the opening one is now proposed.

• Hypothesis 1. The closing mechanism is non-linear and exponential, the crack opening value tends to zero for an infinite compressive stress but the crack will never totally close, see Fig. 4(b).

The exponential form of the law leads to an increasing required compressive stress to close the crack. It depicts a closing getting increasingly harder to close. Moreover, the crack will never wholly closed because the infinite required stress. Finally it leads to a similar form as the opening law, thus benefiting from its simplicity such as the analytical local resolution (described in section 4.4).

• Hypothesis 2. The required energy to completely close the crack equals to the energy dissipated during the opening phase.

This assumption only impacts the value of the closing fracture energy (which can also be seen as a material property). It is depicted in Fig. 4(b) and means

 cl ∶=  ′ op (8)
• Hypothesis 3. In the reopening phase, the opening criterion is computed based on a residual tensile strength and a residual fracture energy.

This hypothesis correspond to a loading phase in an already damaged element. Herein, the degradation of the tensile strength of a damaged material (even if the crack is partially closed) is here accounted for.

As shown in Fig. 4(c), this residual strength is set to be the lowest tensile stress of the previous loading stage.

As results, the damage that have occurred is irreversible and permanent. The element will become more fragile as the opening-closing process continues.

Before going into the details of the equations of the proposed closure law, the next three subsections go through the main phases of the mechanical process.

Unloading procedure Let us first consider the unloading of an element that has already a crack opening. As always, the element will obey an elastic behavior outside of the discontinuity, while the opening value [ ] remains constant. This is because the tensile energy is still stored in the bulk volumes, and the discontinuity interface is under traction before a total release of the elastic energy. It is worth noting that the elastic unloading procedure will end before the imposed displacement Δ reaches 0. More specifically, the elastic unloading procedure will end when Δ = [ ] (the value during the opening procedure), see Fig. 5.

At the state of null stress (and null strain), see point in Fig. 5, the bulk volume has completely released its elastic energy.

Crack closure Following the unloading procedure, the crack will close if a compressive stress is applied to it.

As depicted in Fig. 6, when the imposed displacement is null (Δ = 0) the crack still carries a positive opening value. At this point we have [ ] < [ ] max and the residual fracture energy Reloading procedure Let us now apply another loading displacement to the element.

 re =  op - ′ op . [| |] Δ (a)
In a first stage, the normal stress moves from a negative to a positive value with a constant crack opening [ ] cl . It corresponds to an elastic unloading from its compressive state of the bulk of the element.

Then, when the equivalent stress reaches the residual tensile strength ′ y , the previous opening law drives again the crack opening with a residual fracture energy  re defined just above, see Fig. 7.

Governing equations of the closure mechanism As introduced in previous parts, the closure law is non-linear and exponential, which depicts a "symmetrical" form as the opening law. The choice made here is :

Φ c = -eq +  cl [ ] max ln [ ] [ ] max ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ c , (9)
A positive value of Φ c means the crack should continue closing to a further position while a null value represents an equilibrium stress state. A negative value of Φ c means an elastic loading or opening. It worth reminding that the opening process and the closing process will never happen at the same time, i.e.

Φ o ⋅ Φ c < 0, because [ ] max is always superior or equal to [ ].
The hypothesis is made to have the same required energy for completely closing the crack as the dissipated energy in the opening phase. Hence, regarding the closure law, the dissipated energy in the closing phase should respect the following equation:

 cl = ∫ 0 [ ] max eq ([ ]) [ ] =  ′ op . ( 10
)
And the dissipated energy during the opening phase is calculated as:

 ′ op = ∫ [ ] max 0 eq ([ ]) [ ] =  op 1 -exp - y [ ] max  op . ( 11
)
It can be seen that the dissipate energy can be completely defined by  op , y , and the largest opening [ ] max . Filling Eq. 11 into Eq. 9, it can be noticed that the closure law can also been fully defined by these parameters.

In the reopening phase, the parameters of the opening criterion are taken as the residual tensile strength (marked as ′ y in Fig. 7(c)) and the residual fracture energy (marked as  re in Fig. 6(c) and Fig. 7(c)). The updated opening criterion thus reads:

Φ o = eq -′ y 1 -exp - ′ y  re [ ] , with ⎧ ⎪ ⎨ ⎪ ⎩ ′ y = y exp -y ([ ] max -[ ] cl )  op  re =  op exp -y ([ ] max -[ ] cl )  op . ( 12
)

Finite Element framework and resolution scheme

In the E-FEM model of this study, both weak and strong discontinuities are introduced in elements, which bring to the system additional unknown variables. Therefore, additional equations are required compared to the resolution of a standard FE system. The main purpose of this section is to solved the enhanced system within a classical Finite Element Method framework. The formulation of the problem is based on a three-field variational principle [START_REF] Washizu | Variational methods in elasticity and plasticity[END_REF], which proposes two sets of fields, a standard one and a virtual one. By using the Statically and Kinematically Optimal Non-symmetric 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), these two fields are not enhanced in the same way, which allows the zero mean condition and suitable representation of the discontinuity kinematics. Then the three-field variational formulation is further developed using the Enhanced Assumed Strain method [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF], where the strain field is decomposed into compatible and incompatible parts. By taking the 2 orthogonal condition between the space of the stress field and enhanced strain field, and the property of the Dirac-delta function (∫ ( )  = ∫ ( = ) ), the unbounded enhancement can be included into the variational statement, and the incompatible system can be solved in a special compatible form [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]. We present in the following sections the resolution of the E-FEM model.

Finite Element interpolation

We present in this part the discretization of the system. Since the Finite Element framework of the E-FEM model with double enhancements has been well documented in pioneer works [START_REF] Benkemoun | Failure of heterogeneous materials: 3d meso-scale fe models with embedded discontinuities[END_REF][START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF], it is chosen here to present only the essential statements.

The discretized solving system gives as [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

 ∶= e =1 ̄ int -̄ ext = e =1 ̄ ext -∫ Ω ̌ ( , [| |], [| |]) Ω = 0, (13a)  [| |] = ∫ Ω w ̌ ( , [| |], [| |]) Ω = 0 , (13b)  [| |] = ∫ Ω * , s ̌ ( , [| |], [| |]) Ω = ∫ Ω ∖ * , s,b ̌ ( , [| |], [| |]) Ω + ∫  Ω = 0, ( 13c
)
where

=1
represents the standard assembly operator which groups all e elements into the computational domain Ω, the vector ̄ ext and ̄ int denote the external and internal force, and * s the strong enhancement interpolation matrix of the virtual field which contains the bounded and unbounded part * s ∶= * s,b + * s,u [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

* s,b = - * s , * s,u =  * s . ( 14
)
In this equation, the term * s is equivalent to the operator (• ⊗ ) sym , and can be seen as the Voigh form of the normal vector .

It is recalled that the bulk volumes always respect linear elastic behavior, i.e., Hooke's law, thus the behavior law in Eq. 13 reads [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

̌ ( , [| |], [| |]) = + + = + ( + + w [| |] + [| |]), ∀ ∈ Ω + , --= -( + - w [| |] + [| |]), ∀ ∈ Ω -, (15)
where the several matrices are: +|-the stiffness tensor in Ω +|-, (= ) the standard strain interpolation matrix, w the enhanced strain interpolation matrix corresponding to the weak discontinuity, and s the strong enhancement interpolation matrix for the actual field. Note that the strain fields are all expressed in Voigt notation. The set of unknowns are the classical nodal displacements , and discontinuity variables

[| |] [| |].
Given the zero mean condition, the expression of w is obtained [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

w = ⎧ ⎪ ⎨ ⎪ ⎩ + w = Θ + w , with Θ + = -in Ω + - w = Θ - w , with Θ -= - + in Ω -, ( 16
)
where w is a matrix that carries only information on the normal vector , +|-is the volume for Ω +||-. And s is the equivalent symmetric operator (• ⊗ ∇ ) in Voigt notation, with is an explicitly defined arbitrary function to separate nodes at Ω + from Ω -by setting unit value to nodes at Ω + and null at other part:

( ) = e ∑ =1 with = 1 if node number ∈ Ω + 0 if node number ∈ Ω -. ( 17
)
With Eq. 13, Eq. 15, and the explicit expressions of the interpolation matrix at hand, the resolution of the discrete system is presented in the following parts.

Linearization of the discretized system

The discretized system Eq. 13 is formulated in two levels: Eq. 13a presents the equilibrium of the global system, while Eq. 13b and Eq. 13c are local equations for each element, corresponding to the weak and strong discontinuity. Regarding Eq. 13c, it can be seen that the traction separation law is incorporated in the system through the traction vector [START_REF] Ibrahimbegović | Classical plasticity and viscoplasticity models reformulated: theoretical basis and numerical implementation[END_REF]. Assuming that the discontinuity surface is flat, the traction vector in Eq. 13c can be obtained by calculating the average value of ̌ and weighted by volumes:

= 1 * , s ( + ̌ + + -̌ -). (18)
Depending on the element exhibiting a strong discontinuity is in the opening phase or in the closing phase, the linearization of Eq. 13c should be treated separately.

If the element is in the opening phase, the linearization equation can be obtained by injecting Eq. 18 into the opening criterion Φ o (Eq. 7). Write it in an incremental form, we have [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

ΔΦ o = eq Δ + [ ] Δ[ ] = s * b Δ + s * w Δ[| |] + s * s Δ[ ] + qo Δ[ ]. (19)
And if the element is in the closing phase, the linearization equation is calculated by injecting Eq. 18 into the closing criterion Φ c (Eq. 9). The equation gives as:

ΔΦ c = - eq Δ + c [ ] Δ[ ] = -s * b Δ -s * w Δ[| |] -s * s Δ[ ] + qc Δ[ ]. (20)
In summary, depending if the element is opening or closing, the linearization of the discrete system Eq. 13 can be written in a matrix form:

⎡ ⎢ ⎢ ⎣ bb bw bs wb ww ws s * b s * w s * s + qo ⎤ ⎥ ⎥ ⎦ ( ) +1 ⎧ ⎪ ⎨ ⎪ ⎩ Δ Δ[| |] Δ[ ] ⎫ ⎪ ⎬ ⎪ ⎭ ( +1) +1 = ⎧ ⎪ ⎨ ⎪ ⎩ -( int -ext ) - [| |] -Φ o ⎫ ⎪ ⎬ ⎪ ⎭ ( ) +1 , (21)
or

⎡ ⎢ ⎢ ⎣ bb bw bs wb ww ws -s * b -s * w -s * s + qc ⎤ ⎥ ⎥ ⎦ ( ) +1 ⎧ ⎪ ⎨ ⎪ ⎩ Δ Δ[| |] Δ[ ] ⎫ ⎪ ⎬ ⎪ ⎭ ( +1) +1 = ⎧ ⎪ ⎨ ⎪ ⎩ -( int -ext ) - [| |] -Φ c ⎫ ⎪ ⎬ ⎪ ⎭ ( ) +1 , ( 22
)
where the subscript ( ) states for the time step, and superscript ( ) is the number of iteration. For the sake of clarification, the explicit expression of the assembled stiffness matrix are presented in the appendix.

Solving the system

The linearized system Eq. 21 and Eq. 22 are solved at two levels. The internal variables

(Δ[| |] +1 +1 , Δ[| |] +1 +1
) are firstly solved at the element level for the given displacement at the iteration ( ).

Once (Δ[| |] +1 +1 , Δ[| |] +1 +1
) are solved, with a null value residual  [| |] and opening/closing criterion (Φ o = 0 and Φ c = 0), the global equation can then be solved using a static condensation [START_REF] Wilson | The static condensation algorithm[END_REF]

: sc | | | ( ) +1 Δ | | | ( +1) +1 = - e =1 { int -ext } | | | ( ) +1 , with sc | | | ( ) +1 = e =1 sc | | | ( ) +1 . ( 23
)
Depending on the element is in the opening phase or closing phase, the assembled matrix sc is calculated separately: From Eq. ( 24) and Eq. ( 25), we can see that the modified stiffness matrix sc has the same size of matrix bb . It means that the size of the problem equals to the size of a classical FE formulation without enhancements, no matter the number of heterogeneities or cracks, thus easing the needed computation power to solve meso-scale problems. However, it is recalled that at the local level, an element carrying a strong discontinuity requires the resolution of a non linear equation, which in turns can increase the computation time. The next section describe an analytical solution to these local equations.

Resolution of the traction separation law

The cracks start opening after the localization occurs. It contains a softening function o to simulate the failure mechanism for a quasi-brittle material. Being an equation in terms of the opening value [ ], the opening criterion can be reformed into [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

+ [ ] = y exp - y  op [ ] , (26) with
= s * b -s * w -1 ww wb , = s * s -s * w -1 ww ws . ( 27
)
In this equation, s * b , s * w , ww , wb are constant matrix. is constant within the local resolution which depends proportionally on . As for , it is determined by the time traction vector is defined. A solution can only be solved if > y and < 0, see Fig. 8. With the help of the Lamber W function [START_REF] Corless | On the lambertw function[END_REF] 0 , an analytical solution gives as [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

[ ] sol =  op y ⎛ ⎜ ⎜ ⎜ ⎝ 0 ⎛ ⎜ ⎜ ⎜ ⎝ 2 y exp( y  op )  op ⎞ ⎟ ⎟ ⎟ ⎠ - y  op ⎞ ⎟ ⎟ ⎟ ⎠ . ( 28
)
A closing process in the strong discontinuity of mode-I represents a decreasing value of the crack. The governing law for the closure mechanism is non-linear:

Φ c = -eq + c , with c =  cl [ ] max ln [ ] [ ] max . ( 29
)
In order to solve the equation at the local scale of Eq. 22, let us assume that the local system has already been solved, with and

 [| |] = 0 and  [| |] = 0: wb + ww [| |] + ws [ ] = 0, (30)
-s * b -s * w [| |] -s * s [ ] +  cl [ ] max ln [ ] [ ] max = 0. ( 31
)
By replacing all terms of [| |] by and [ ], we can obtain the following equation:

-s * b + s * w -1 ww wb + -s * s + s * w -1 ww ws [ ] +  cl [ ] max ln [ ] [ ] max = 0. (32)
Hence, we can reform this function into this shape:

-

-[ ] =  cl [ ] max ln [ ] [ ] max , ( 33
)
where and have the same definition as in Eq. 27, and the analytical solution for [ ] can be solved as:

[ ] sol = - cl 0 ⎛ ⎜ ⎜ ⎝ - [ ] 2 max exp( [ ] max  cl )  cl ⎞ ⎟ ⎟ ⎠ [ ] max . ( 34
)
The resolution of the crack opening [ ] sol is plotted in Fig. 9. Once the opening value [ ] sol is determined, the value of weak discontinuity [| |] sol can be analytically calculated as [START_REF] Roubin | Modélisation EF et morphologique de milieu hétérogènes à l'échelle mésoscopique: applications aux matériaux à matrices cimentaire[END_REF])

[| |] sol = --1 ww wb + ws [ ] sol ( 35
)
This framework is implemented in the FE code FEAP [START_REF] Taylor | FEAP-ein finite element analysis programm[END_REF]. Solving the nonlinear global equations is performed through a quasi-Newton BFGS resolution scheme [START_REF] Markovic | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior[END_REF].

Test on a single element

In order to verify the validation of the proposed model, we consider here an example with one single tetrahedral element, see Fig. 10. It is a homogeneous element. Thus it can have only the strong discontinuity. Applying a uniaxial tension at the node 2 towards to direction , a discontinuity interface is obtained at the centroid of the element once the localization is triggered, see Fig. 10(b). The normal vector of the discontinuity surface is determined by the major principal stress, in this case, following the direction of (0, 0, 1).

The elastic parameters and failure parameters of the element are given in Table . 1. In order to investigate the performance of the closure law, we apply a cyclic tensile loading with increasing amplitude to the model with and without closure law.

It can be seen from Fig. 11 that first, for the same loading path, the two models have the same opening value [ ] during the first opening procedure. This is reasonable because they have the same opening criterion. Secondly, focus on their most evident difference: the closure procedure. Owing to the closure law, the model with closure mechanism Table 1 Considered material parameters for the single element. shows an ability to close the crack and dissipate an amount of energy at the closing phase. The decreasing value of [ ] does not exist in the model without closure mechanism, which leads to a consequence that the element carries a much bigger value of the compressive stress on the discontinuity interface. This value of the compressive stress comes from the response of the bulk part (Ω ∖) which is considered as always elastic itself. In this context, the model with closure law will dissipate more energy not only in the closure procedure but also in the opening procedure.

E [GPa] [-] [MPa]  op [ ∕ 2 ]

Illustration of the performances of the model to heterogeneous materials

The heterogeneous morphological structure of the concrete at the mesoscopic scale is complex which exhibits hard aggregates and macro-pores. This section begins by applying the model to an example which carries a simple morphological structure. The used morphological model consists of one single spherical heterogeneity in the matrix. The interests are made on the emergent phenomena in the macroscopic scale and the crack patterns observed in the specimen. Then the model is applied to a material that carries a more complex heterogeneous morphological structure.

Note that for the following numerical examples, the uniaxial loadings are displacement-controlled and applied in the z-direction. The lateral faces (normal to x-direction and y-direction) are set as free-stress. The heterogeneities are considered as purely elastic for the sake of simplicity. For the interface elements between the heterogeneities and the matrix, the critical tensile strength is defined as the same as the matrix. Since the elastic parameters are defined differently for matrix and heterogeneities, there will be stress concentrations at the interface. Using the same failure parameter in the matrix and the interface elements will lead to a weaker behavior at the latter.

Cube with an internal sphere

In this section, a cube with one single spherical heterogeneity is proposed, see Fig. 12. The specimen has a length of 100 millimeters and carries an internal sphere of 30 millimeters' radius. Two cases are considered: the internal sphere is harder compared to the matrix, representing an aggregate; or the internal sphere elements are deleted from the original system to simulate a macro-pore. The parameters of the material are given in Table . 2.

Table 2

Considered material parameters for the cube and its internal sphere.

Phase

E [GPa] [-] [MPa]  op [ ∕ 2 ] Matrix 20.0 0.2 4.0 0.1 Interface - - 4.0 0.1 Aggregate 100.0 0.2 - - Macro-pore - - - -
A uniaxial loading is applied to the sample. The reaction of the specimen is plotted by the macroscopic stress versus the strain in Fig. 13. We can see that first, the tensile and compressive loading have the same macroscopic Young's module, which equals to 23.9 GPa for the cube with hard aggregate, and 15.6 GPa for the cube with a macro-pore. As expected, this value is between the proposed Young's module of the matrix (20.0 GPa) and aggregate (100.0 GPa) or pore (0.0 GPa), which fits the Hashin-Shtrikman range [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF].

Second, it is worth noting that for a pure homogeneous material, the simple compressive loading won't cause any failure, because the maximum principal stress is negative and will never exceed the critical stress y . However, in this case, with a great number of elements and explicit representation of heterogeneity, the macroscopic failures can be observed in both traction and compression. And clearly, an asymmetrical macroscopic response for tensile and compressive loading is obtained, which is a typical propriety for quasi-brittle materials. Moreover, a reasonable ratio between the compressive and tensile strength can be estimated between 12.2 (cube with hard sphere) and 9.9(cube with a macro-pore). Then focus is made on the effect of the closure mechanism. The differences between the model with or without closure law are illustrated in Fig. 13. In the aspect of the macroscopic responses, there are little differences between the two models in traction. Whereas in compression, it can be observed that the model with closure law behaves a little more resistant feature at the post-pic stage. This phenomenon is related to the appearance of crack closures, which allows the model to dissipate more energy. The performances of the cube with a hard aggregate and a macro-pore are discussed separately in the next part.

Case of a cube with hard sphere

In this part, our focus is made firstly on the performances of the specimen in uniaxial traction loading. The evolution of the crack pattern is shown in Fig. 14.

We can see from the figure that first, outside the two major cracks, there are no micro-cracks diffuse in the cube. This is because the stress concentration is highly localized above and below the internal sphere. And also, the fracture energy of the material is very small ( op = 0.1 ∕ 2 ), meaning that the material is quite brittle.

Next, it can be seen that at the beginning of the failure procedure, the two major cracks are rather symmetric. Then, since the elements will become more brittle as the crack opening value increases, the balance between the two symmetric major cracks is broken at some point during the softening procedure. And one of the macro-crack meets a catastrophic failure, which leads to stress release to another macro-crack. Therefore, there is only one macro-crack across the specimen at the end of the loading.

In order to illustrate the degree of closing, we show in Fig. 15 the percentage of closure and the dissipated closing energy. The percentage of closure is calculated by

# percentage of crack closure = [ ] max -[ ] [ ] max × 100%, (36)
with [ ] max is the maximum crack opening values, and [ ] is the value at the current loading step. It can be observed that the crack above the sphere is almost totally closed, with the closing percentage roughly reaches to 100%. As for the dissipated energy in closing phase, as expected, it is always smaller than the fracture energy which equals to 0.1 ∕ 2 . In the next part, our attention is moved to the unilateral effect. It is a common sense that the unilateral effect is one of the typical characteristics of the concrete observed in experimental tests. It represents a phenomenon that when the loading switches from tension to compression, the material has its stiffness recovered after been damaged. The proposed specimen is now applied to a loading-unloading movement. The macroscopic response of the model is plotted in Fig. 16.

From Fig. 16, it can be observed that the model with crack closures has its stiffness recovered after been damaged at the traction loading, i.e. the unilateral effect can be observed. At the unloading phases, it can be seen that the stiffness of the material for the model with crack closures progressively recovers, and leads to a small amount of plastic deformation at the point (2). This plastic deformation corresponds to the residual crack openings remain in the material at the stage of null stress. As for the model without crack closures, the elastic stiffness of the material seems like has never been damaged by the development of cracks. This is because, for the E-FEM model with strong discontinuities, damages and cracks take place at discontinuity interfaces, while the sub-domains always carry a pure elastic behavior. Hence, the unloading phase is purely elastic, and higher residual openings lead to more significant plastic deformation. We present the crack patterns of the material of three different loading steps in Fig. 28. In order to illustrate the effect of the crack closure mechanism, we put the model with crack closures in the left and the model without closures in the right for comparing the differences between them. The chosen three loading steps are marked in Fig. 16, with the macroscopic strain equals to 0.95 × 10 -4 , 0.085 × 10 -4 , and -0.36 × 10 -4 .

It can be seen from Fig. 28(a) that at the end of the traction loading, the material is heavily damaged with two major -0.5 0 0.5 1

⋅10 -4 -1 0 1 2 (1) (2)
(3) macroscopic cracks traverse the cube. The crack patterns of the model with or without crack closures are quite similar. And in Fig. 28(b), for the model with closure mechanism, there are a group of localized elements that have closed while some elements still opening when the unloading phase arrives at the null stress stage. The plastic deformation at the state of unloading is the consequence of these still opening elements. And in Fig. 28(c), the fractured elements are almost completely closed under the compressive loading, the stiffness of the material is fully recovered at this point, the axial stress is also turned into compression. Unlike the model with crack closure mechanism, the model without crack closures produces a crack pattern that barely changed during the unloading and compressive loading. The already created cracks are not decreased. Speculated in the macroscopic response in Fig. 16, it leads to a high value of the residual deformation at the state of null stress, and also higher compressive stress than the model with closure mechanism.

Strain [-] Stress [MPa]

Case of a cube with macro-pore

In this part, the elements of the internal sphere are deleted from the original system to represent a macro-pore in the cube. We can see from Fig. 17(a) and Fig. 17(b) that only one macro-crack apparent in the middle of the cube, the crack pattern is highly similar for the model with and without closure mechanism. And it is can be seen from Fig. 17(c) that there are almost no elements obey a crack closure for the model with closure mechanism. This is due to the fact that there is only one major crack occurs in the specimen, hence no crack is located at the stress release zone.

The above observations lead to the conclusion that the closure of cracks exists not only in cyclic loading but also in simple uniaxial loading. The crack closure is related to the interaction between cracks, that the stress concentration between them will be released as the cracks develop and propagate.

Test for a concrete-like material

Let us now consider in this section a morphological more complex sample. The heterogeneous morphological structure is represented by numerous spherical aggregates dispersed randomly in the cube and may touch the boundaries, see Fig. 18. Two sets of spheres are used to represent a total volume fraction of 20%, which are respectively 60% between 3 and 5 millimeters, and 40% between 8 and 12 millimeters. The main idea of the used method [START_REF] Vallade | Modélisation multi-échelles des shales: influence de la microstructure sur les propriétés macroscopiques et le processus de fracturation[END_REF] is to generate N spheres of different radius and place them randomly in the field while guaranteeing a minimum distance between each sphere. The material's parameters are listed in Table . 3, which are chosen to represent concrete at the mesoscopic scale with the matrix plays the part of the cement paste, and the heterogeneities are hard aggregates.

Applying uniaxial loadings to the cube, the responses of the specimen are displayed by the macroscopic stress in Table 3 Considered material parameters for the concrete-like specimen.

Phase

E [GPa] [-] [MPa]  op [ ∕ 2 ] Matrix 22.0 0.2 4.0 1.0 Interface - - 4.0 1.0 Aggregate 78.0 0.2 - -
terms of the strain, see Fig. 19. The model with and without closure mechanism are compared in this figure to illustrate the effect of crack closures. The observed macroscopic Young's module of this sample is equal to 27.4 GPa.

In the previous example of the cube with an internal sphere, several typical behaviors of concrete can be observed. These behaviors can also be seen in this case of the cube exhibiting a more complex heterogeneous morphology. The asymmetric failure response in traction/compression is observed with a ratio equals to 13.49. The differences between the model with or without crack closures are very slight for traction, and relatively more obvious for compression. This is due to the fact that the model with the closure mechanism will dissipate more energy in crack closures. Since the specimen presents a lot more cracks in compression than in traction, the more ductile behavior is therefore more obvious in compression than in traction.

To analysis the evolution of the crack pattern, we choose here four loading stages along with the traction loading

⋅10 -4 0 1 2 3 (1) (2)
(3) (4)

Strain [-]
Stress [MPa] (b) at the tensile responses which corresponds to the state before failure, the state after the failure, relatively damaged state, and very damaged state, see Fig. 19(b) and Fig. 29. Fig. 29(a) plots the crack patterns at the loading stage number (1), the macroscopic response at this stage is just before the maximum resistance of the material. A number of micro-cracks can be observed dispersed over the cube, especially beside the inclusions because of the stress concentration. There are no crack closure exhibits in the cube at this phase.

Then the cracks at the post-failure stage are shown in Fig. 29(b), which corresponds to the stage (2). The failure behaviors of the specimen are closely related to the development of the macro-cracks. It can be seen that the number of cracks heavily increases, and the micro-cracks link with each other and form some macroscopic cracks. A very small number of crack closures apparent inside the material.

From stage (3) to stage (4), the macroscopic response of the damaged material tends to steady, and the crack closures broadly develop inside the material. For the sake of clarity, we label three major cracks, see Fig. 29(d). It has been introduced in the previous part that the crack closures are triggered by stress releases, and the stress releases are closely related to the interaction between major cracks. Depending on the evolution of major cracks, we can roughly classify the closing pattern into two main categories.

The first crack closing patterns can be represented by crack 1 and 2, as shown in Fig. 29(c). These two cracks are closed because there is another crack that develops faster than them and results in stress releases. From Fig. 29(b) to Fig. 29(d), it can be seen that crack 1 and 2 are almost entirely closed. As for the model without closure law, the stress releases will not produce any crack closures. However, the consequences of the variations of the stress field can also be observed in Fig. 20. We draw in this figure the increasing value of the crack opening from the stage (3) to the stage (4), it can be seen that crack 1 and 2 do not continue propagating as the loading increases.

The major crack number 3 belongs to another category of the closing pattern. Since elements at the center of cracks tend to develop faster than the others, the stress release will induce the crack closures for the bordered elements of a macro-crack. As for the model without closure law, it can be seen from Fig. 20 that the increasing fractures are also centered in the middle of the crack 3. At the border of these major cracks, the crack opening values are barely enlarged. Compare the crack pattern at the end of traction between the two models, it can be noticed that the model without crack closures (Fig. 20) produces a crack pattern that is more dispersed over the material, while the cracks for the model with crack closures (Fig. 29(d)) are deeper and thinner.

We plot the variations of crack opening values in traction for the model with and without closure law in two histograms, see Fig. 21. For the model with crack closures, it can be seen that after the crack closures apparent in the material, i.e. after the stage (2), the number of the element which carries a bigger value of crack opening increases, and so does the element which carries a smaller value of crack opening. The former corresponds to the continue-opening elements which are located at the center of major cracks, they carry a relatively higher value of crack opening (for example, the major crack 3 in Fig. 29(c)). And the latter corresponds to the closed cracks induced by stress releases, for example, the major crack 1 and 2 in Fig. 29(d). As for the model without closure law, we can see from Fig. 21(b) that it only carries an increasing number of the element which exhibits bigger crack opening value, and the number of the element which carries a relatively smaller crack remain constant.

Next, we apply a loading-unloading displacement to the specimen to illustrate the ability of the model in the aspect of the unilateral effect. The macroscopic responses of the model is plotted in Fig. 22(a) and Fig. 22(b). It can be seen in Fig. 22(a) that, though the material is very damaged in traction, the elastic stiffness under compression is fully recovered. These observed behaviors are linked with the closure of cracks, see the diagram in Fig. 22(c). It is worth noting that the curve in Fig. 22(c) illustrates the response of a concrete beam under cyclic bending, thus the curve itself is not comparable with our simulation results. The resources of the stiffness recovery and the plastic deformation, though, should be the same. It can be observed in Fig. 22(a) that the stiffness of the material progressively recovered as the imposed tensile displacement reduces and changes to compressive displacement. Also, at the state of null stress, the material carries an amount of permanent opening. Compared with the model exhibits only one sphere that we presented in section 5.1, the plastic strain is much more significant. This is because the fracture energy of the material ( op = 1.0 ∕ 2 ) is bigger than the fracture energy of the cube with an internal sphere ( op = 0.1 ∕ 2 ) in the previous section. Therefore, the material is more ductile and requires more energy to close the cracks. Besides the stiffness recovery and the plastic deformation, a hysteresis loop can also be seen in Fig. 22(a), which is also consistent with experimental observations. This is related to the additional dissipated energy and the stiffness recovery which are both produced by the crack closure mechanism.

Interested in the crack patterns of the model with the closure mechanism along with the unilateral loadings, we present the cracks of the material at different states in Fig. 23. It can be noticed that the evolution of the crack patterns in Fig. 23 shows a good consistency to the behaviors of concrete in Fig. 22(c), that: • At the point (1), a large amount of crack openings apparent in the specimen, accompanied with the loss of stiffness of the material.

• At the point (2), the state of null stress, the cracks are partially closed while leading to an amount of plastic strain, the resource of which is the residual openings.

• At the point (3), the cracks are almost entirely closed, leading to the recovery of the stiffness of the specimen. Now let us move our attention to the performance of the model under compressive cyclic loadings. The trajectory of displacement gives in Fig. 24. It is a cyclic loading with increasing amplitude. The macroscopic responses are shown in Fig. 25.

We can see from Fig. 25 that the first unloading occurs in the elastic phase. Then in the subsequent cyclic loadings, the curves are all enveloped in the curves of monotonic loadings. The volumetric dilatation is observed in both models while the model with closure mechanism has less dilatation effect. This is because the expansion of the material is related to the "opening" behavior, and the closure of the cracks will restrict this effect.

The effect of the closure mechanism can be illustrated by comparing the responses of the two models in Fig. 25 and Fig. 26. Among several differences between them, the most remarkable observation is the hysteresis phenomenon. The emergence of the hysteresis loops at the macroscopic scale takes its origin from several constituents on the mesoscopic scale:

• The elements dissipate additional energy as the closing phase, which corresponds to the area of the hysteresis loops at the macro-scale.

• The model with the closure mechanism loses less of the stiffness of material compared with the model without closure mechanism owing to the closing of cracks. Drawing the macroscopic stiffness of the material in Fig. 26 in dashed lines, we can see that the model with closure law has 42.3% of the initial stiffness at the end of the loading, while the model without crack closure only has 18.5% of the initial stiffness.

However, it can also be noticed that the hysteresis loops are relatively small regarding the experimental observations. We believe that this is due to several reasons:

• The tensile strength and the fracture energy are eroded during the cyclic loading, and this procedure is irreversible. This leads to a result that the additional dissipated energy which forms the hysteresis loops will become increasingly smaller and tend to zero after several cycles.

• The mode-I discontinuity is only capable of depicting the "opening" cracks. It is common sense that friction plays an important role in failure behaviors of concrete in compression. Hence, the insufficient hysteresis phenomenon is one of the consequences of lacking necessary mechanisms.

In summary, the occurrence of the hysteresis phenomenon is closely related to the added closing mechanism. Yet, it can be noticed that the model proposed in this study cannot produce enough hysteresis loops, some other mechanisms could be indispensable such as the consideration of the friction slip between the lips of the cracks.

As supplementary, we draw the variation of the number of localized elements and closing elements in Fig. 27. First, as we can see from the figure, the model with and without crack closures both have an amount of fractured element at the first cycle, even though the material is still at the elastic phase. This is because, at this stage, there are micro-cracks scattered over the specimen, but not coalescent to form macroscopic cracks. Second, the figure shows a similar increasing tendency in the number of localized elements during the cyclic loading for both models. Such curves are observed because the localization criterion, which decides the appearance of the cracks, is formulated as the same for both models. Finally, it can also be observed that the number of closing elements for the model with crack closures increases and decreases during the cyclic loading, while the number of localized elements stay constant at the several last cycles. This declares that not all of the localized elements will turn to a closing state when unloading displacement is applied, the closing and reopening elements are the same groups of elements.

Discussions and conclusions

In this paper, the aim is to assess a FE model to simulate the macroscopic behavior of concrete by upscaling from the mesoscale. The used method is the Enhanced Finite Element method (E-FEM). Based on the existing E-FEMs, the model incorporates two kinds of enhancement at the meso-scale: the weak enhancement and the strong enhancement. The former enables the model to explicitly simulate the heterogeneities. The latter reproduces jumps of the displacement field in the material, representing cracks at fine scale. These two kinds of enhancements are independent and additive, making the model capable of simulating accurate and complex cracks, such as the debonding of the interfaces and the branching of cracks.

Existing E-FEM models have shown their ability to simulate many of the main characters of concrete, such as the asymmetric behaviors in traction/compression, the volume dilatation, and the progressive loss of stiffness. In this paper, a closure mechanism is integrated within the model. The closure criterion is proposed under physical and meaningful considerations and requires no additional parameter. A number of numerical applications of the model have been shown through the modeling of concrete-like specimens. Compared with the model without the closure mechanism, some emergent macroscopic responses can be obtained:

• It can be observed that the crack closures exist in both cyclic and monotonic loadings. The existence of the closure mechanism in the monotonic loading will result in deeper and thinner crack patterns.

• The unilateral effect is observed, and the evolution of the crack pattern shows a good agreement with the experimental observations.

• The hysteresis phenomenon emerges from crack closures in the cyclic loading, but the hysteresis loops are not very significant.

From the simulations in section 5, a limitation of this study can be noticed that the material is incapable of reaching sufficient softness, especially under the compressive loading, namely the stress locking effect. After the macroscopic failure of the material, its macroscopic stress is not equal to zero, even under a substantial-high loading. Among several reasons, the most important one should be related to a set of interface elements that cannot meet the failure criterion. This is because the orientation of the discontinuity is predefined by the geometry of heterogeneity, and the separation is pure mode-I. Therefore, the stress locking effect will manifest itself when the loading results in shear stress, especially for a high volume fraction material, which raises the ratio of the interface elements.

Another limitation of the model is the insufficient hysteresis loops in the cyclic loading. This is related to the lack of consideration of the friction. The friction between the lips of micro-cracks plays an essential role in producing hysteresis phenomena of quasi-brittle materials.

In order to overcome the limitations mentioned above, a thinner mesh can be used to reduce the stress locking effect, especially for the target solid with high volume fraction. However, a thinner mesh means an increase in the computation cost, requiring more necessary memories and longer calculation time. Therefore, adding the mode-II opening can be a good solution to the two problems presented above.

Figure 1 :

 1 Figure 1: Reference body exhibiting heterogeneities and cracks and their corresponding enhanced elements.

Figure 2 :

 2 Figure2: Three cases of discontinuities on a tetrahedral element[START_REF] Roubin | FE mesoscopic modelling of a micro-concrete based on x-ray scanmorphologies[END_REF].

Figure 3 :

 3 Figure 3: Local constitutive model at the continuous part and at the discontinuity interface at the opening stage.

Figure 4 :

 4 Figure 4: Hypothesis for local constitutive models with non-linear closing laws and damaged tensile strength at the discontinuity interface.

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Local constitutive model at the continuous part and at the discontinuity interface as the imposed displacement decrease.

Figure 8 :

 8 Figure8: Resolution of the strong discontinuity equation in the opening procedure[START_REF] Roubin | Modélisation EF et morphologique de milieu hétérogènes à l'échelle mésoscopique: applications aux matériaux à matrices cimentaire[END_REF].

Figure 9 :

 9 Figure 9: Resolution of the strong discontinuity equation in the closing procedure.

Figure 10 :

 10 Figure 10: Basic geometric information and boundary conditions for the studied single tetrahedral element.

Figure 11 :

 11 Figure 11: Comparison between the model with/without closure law for a single tetrahedral element.

  (a) Morphology of the cube with one single internal sphere (b) Projection of the morphological structure to an unstructured mesh

Figure 12 :

 12 Figure12: Morphology of the cube with an internal sphere and its projection to an unstructured mesh.

  Figure 13: Comparison of the macroscopic responses between the model with/without crack closures for uniaxial loadings.

Figure 14 :

 14 Figure 14: Variation of the macroscopic crack patterns for simple traction applying to the cube with hard internal sphere (the cube is split into half in the lateral direction).

Figure 15 :

 15 Figure 15: Illustration of the degree of closing and the dissipated energy in closing phase for all localized elements at the loading step = 10 × 10 -5 .

Figure 16 :

 16 Figure 16: Macroscopic responses of the cube with a hard internal sphere under unilateral loadings for the model with/without crack closures.

Figure 17 :

 17 Figure 17: Observed crack pattern of the cube with internal macro-pore at the tensile loading of = 10 × 10 -5 .

Figure 18 :

 18 Figure 18: Illustration of the heterogeneous morphological structure.

Figure 19 :

 19 Figure19: The comparison between the model with/without closure law for the cube under monotonic loading.

Figure 20 :

 20 Figure 20: Cracks for the model without closure law at the stage (4), and the increasing values of crack opening from the stage (3) to the stage (4).

Figure 21 :

 21 Figure 21: Histogram of crack opening values for the model with and without closure mechanism in traction.

Figure 22 :

 22 Figure 22: Illustration of the unilateral behaviors of concrete.

Figure 23 :Figure 24 :

 2324 Figure 23: Variation of the cracks of the material during the unloading process in a unilateral test, the points are marked in Fig. 22(a).

Figure 25 :Figure 26 :

 2526 Figure 25: Macroscopic response of monotonic/cyclic tests for model with/without adding closure law, in terms of axial (note Ax) transversal (note Tr) and volumetric strain (note Vol).

  Localization elements for model with closureLocalization elements for model without closure Closing elements for model with closure

Figure 27 :

 27 Figure 27: Accumulate number of localized elements and closing elements along the 4 cycles for model with/without crack closure mechanism.

Figure 28 :

 28 Figure 28: Illustration of the crack pattern and the axial stresses (shown in the demi-cube) of the cube with a hard internal sphere at three different stages of unilateral loading, with the model with crack closures in the left and the model without closures in the right.

Figure 29 :

 29 Figure 29: Cracks (in the middle) and closing percentages (in the right) at four different loading stages for simple traction, they are overlaid in the left with the rigid inclusions are marked in gray.

First Author et al.: Preprint submitted to Elsevier

Appendix -Expressions of the assembled matrix

This appendix aims at presenting the explicit expressions of the assembled matrix shown in Eq. 21 and Eq. 22. Readers can refer to [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF] for details of developments. s * b = eq 1 * , s ( + + + --)