Emmanuel Roubin

Alexis Vallade

Nathan Benkemoun
email: nathan.benkemoun@univ-nantes.fr

Jean-Baptiste Colliat

Multi-scale failure of heterogeneous materials: a double kinematics enhancement for Embedded Finite Element Method

Keywords: heterogeneous quasi-brittle material, strong discontinuity approach, E-FEM method, EAS method heterogeneous quasi-brittle material, strong discontinuity approach, E-FEM method, EAS method

published or not. The documents may come

Introduction

Fracture in brittle or quasi-brittle materials -such as geological media like clay or cementitious ones -may occur under several loading paths such as tension or compression.

In this paper, following the idea developed in (Benkemoun et al. (2010)) and (Benkemoun et al. (2012)), a new paradigm is presented: the use of the strong discontinuity approach at the mesoscopic scale. In this context, the strong discontinuity approach becomes a tool to model (1) a diffuse process of cracking and (2) the coalescence of micro-cracks leading to one or several macro-cracks and eventually to the failure of the specimen. This massive use of strong discontinuities leads to a large ratio of the number of elements for which strong discontinuities are activated to the total number of elements in the mesh and so to a quite complex and precise shape of the macroscopic crack(s). Moreover the large number of micro-cracks represented also offers a way to model macro-scale crack path continuity and branching as a result of micro-cracks coalescence.

At the mesoscopic scale (e.g. millimeter scale for concrete), heterogeneities play a significant role: for concrete, they are aggregates inclusions embedded within a mortar matrix. In order to explicitly take them into account within a Finite Element context, the authors turn to a method referred as non-adapted meshing process (Moës et al. (2003)). This consists in a unique "homogeneous" mesh with nodes placed independently from the heterogeneous morphology and thus from the interfaces. Hence, for a given micro-structure geometry, a set of elements are split into two parts by a physical interface. These elements are decomposed into two sub-domains that inherit of each material phase properties. In order to introduce this contrast of property in the classical linear elements, an enhancement of the kinematics is performed by means of a jump within the strain field -this being known as a weak discontinuity. The pioneer work on that matter goes back to the end of the 80's (Ortiz et al. (1987)) for the modeling of localization bands. The application of weak discontinuities in the context of material heterogeneities modeling has been introduced later in (Sukumar et al. (2001)).

Dealing with two-phase materials (e.g. matrix-inclusion), the non-adapted meshing method leads to two classes of elements: the ones that are completely within the matrix or within the inclusions, and the elements split into two parts. Only the latter type of elements are enhanced by the mean of weak discontinuities. Hence, the meshing process, referred to as morphological projection, consists in assigning their type to the original elements of the mesh. This methodology can be applied to three-dimensional random shapes, allowing complex geometry (based for example on excursion sets of random fields (Roubin (2013), Roubin et al. (2014)). This point is depicted on figures 1(b) and 1(c) (where matrix, inclusion and interface elements are represented in light, medium and dark grey, respectively). A weak discontinuity (Hautefeuille et al. (2009)) (finite jump in the strain field) is introduced into each element in order to take into account this particular kinematics. This technique provides a meshing process that does not depend on the position and the shape of the inclusions. Cracking and softening behavior observed for quasi-brittle material such as concrete are introduced thanks to a strong discontinuity (Oliver et al. (2002)) (finite jump in the displacement field) in each elements. These discontinuities represent micro-cracks that can occur in any phases (aggregates or mortar matrix) and the debonding at the interface between aggregates and mortar matrix.

The remainder of the paper begins in Section 2 with a description of the kinematics for weak and strong discontinuities with which we work, and is then organized as follows: in Section 3, we introduce those kinematics into a Finite Element model, in Section 4, a general Finite Element discretization of the problem is presented and its resolution is described and applied to 4-node tetrahedron elements in Section 5. Finally, in Section 6, attention is focused on the capability of the model to represent the main features of cementitious materials by means of mechanical loading. First uniaxial tension and compression loading are performed. On the one hand, we show the emergence of the typical asymmetry of the respective macroscopic responses. On the other hand, we make a complete review of the crack patterns. The observations enlighten the close relationship between these two phenomena. In addition, non proportional loadings are performed in order to analyze the induced anisotropic behavior due to the failure process. This study is carried out on the previous uniaxial tension and compression tests by means of damage indicator -linear post-analysis on Young modulus degradation -as well as residual strength -non linear post-analysis on tensile strength degradation.

Kinematics description of weak and strong discontinuities

In this section, we summarize the basic notations employed all along the paper and remind kinematics of both weak and strong discontinuities.

We denote by Ω ⊂ R 3 a domain with smooth boundary ∂Ω and the spatial position of a particle by x ∈ Ω. This domain is discretized by means of standard isoparametric elements such as Ω = ∪ n elm e=1 Ω e . Naturally, position of weak discontinuities are defined by the material heterogeneities surfaces (interfaces) and therefore known prior to any mechanical calculation.

They are represented on figure 2 by the smooth contours of the grey shapes which have to be seen as the heterogeneities. In opposition, strong discontinuities (cracks) are introduced thanks to a stress-based localization criterion. Hence they can be, theoretically speaking, scattered anywhere in the domain (red curves on figure 2). However, in the presented framework, the physically meaningful assumption that they can be present only in the matrix (white zone) or following material heterogeneities surface is made. Whatever the origin of the discontinuity is, it will be called Γ d through this paper for sake of clarity.

Dealing with both weak and strong discontinuities, three cases can be present in an element Ω e : 1) only a weak discontinuity is present, 2) only a strong discontinuity is present, 3) both of them are present. The usual representation of (Jirásek (2000)) on figures 3(a) and 3(b) can be extended in the case of double enhancement as depicted on figure 3(c). In the first case, the shape of the material surface Γ d is defined by the shape of the heterogeneity, the direction of the unit normal n is taken orthogonal to Γ d and a jump of material properties passing through Γ d is present. In the second case, the direction of the unit normal n is defined by a standard principal stress criterion and a displacement jump passing through Γ d has to be considered. In the last case, the shape of the material surface Γ d is defined by the shape of the heterogeneity, the direction of the unit normal n is taken orthogonal to Γ d and a jump of material properties and of the displacement passing through Γ d are present.

It implies that the crack opening will be localized on Γ d or in other term on the physical boundary of each phase. As it is depicted on figure 2, such elements are split into two parts Ω e and Ω e so that Ω e = Ω e ∪ Ω e .

Ω

n Finally having at hand these two discontinuities, the total strain is written in the context of the enhanced strain fields (Simo and Rifai (1990)) such as:

Ω e Γ d Ω e Ω e Ω e Γ d Ω e n Ω e Γ u
ε = ∇ s ū regular + ε weak + ε strong , (1)
where ∇ s ū is the symmetric gradient of the displacement field. As in (Simo and Rifai (1990)), we refer to ε and ε as the enhanced parts of the strain field. The notation • (resp.

•) refers to weak (resp. strong) discontinuity.

Having those considerations in mind, we turn now to kinematics of weak and strong discontinuities.

Weak discontinuity kinematics

In this part, interest is made only on weak discontinuity kinematics thus only ε is considered. In order to determine the form of ε, considerations on its corresponding displacement field labelled ũ are taken. The description of such a displacement field depends naturally on its expression in Ω e and Ω e . Let's call ũ the displacement field in Ω e and ũ the one

in Ω e . A first natural consideration is to impose continuity of ũ passing through Γ d . This yields to

ũ = ũ ∀ x ∈ Γ d . (2)
Whereas ũ itself remains continuous, the displacement gradient shall exhibit a jump across the surface discontinuity. Maxwell's compatibility conditions (Gurtin (1984)) lead to the form,

ũ = Θ n • (x -ξ) ([|ε|] n n + [|ε|] m m + [|ε|] t t) with Θ =    Θ ∀ x ∈ Ω e Θ ∀ x ∈ Ω e , (3)
where ξ represents the position of Γ d , thus letting n • (xξ) act has the signed distance to the surface discontinuity, and Θ a still unidentified piecewise constant shape function with unit jump at discontinuity surface (a specific definition is given in section 3) . The introduction of the Θ function leads to a different definition of ũ weather it is evaluated in Ω e or in Ω e . In (3), [|ε|] n , [|ε|] m and [|ε|] t are three constants that can be identified as the strain jumps. Actually, by considering the symmetric gradient of the displacement field ũ of components (ũ n , ũm , ũt ), ε can be written:

ε = ∇ s (ũ) =      ũn,n 1 2 (ũ n,m + ũm,n ) 1 2 (ũ n,t + ũt,n ) sym ũm,m 1 2 (ũ m,t + ũt,m ) sym sym ũt,t      , (4)
where (• ,n ), (• ,m ) and (• ,t ) represent derivatives in relation with n, m and t.

Then the difference between ε and ε can be expressed as follows

ε -ε =      [|ε|] n [|ε|] m [|ε|] t [|ε|] m 0 0 [|ε|] t 0 0      , (5)
where [|ε|] n , [|ε|] m and [|ε|] t stand for (ũ n,n -ũn,n ), (ũ n,m -ũn,m ) and (ũ n,t -ũn,t ), respectively.

All the specificity of the weak discontinuity kinematics lies on the enhancement displacement field ũ and its related strain field ε. The form of ε is now known: Θ has to be still exactly defined according to additional statical considerations. One can already note that the enhancement brings three parameters (stored in the so-called weak discontinuity vector [|ε|]) that will be future unknowns of the mechanical problem.

Strong discontinuity kinematics

In this part, an element Ω e is only split by a discontinuity in the displacement field.

Thus only ε is present. Since strong discontinuity kinematics has been widely studied in the literature (see Simo et al. (1993), Oliver (1996) and Wells and Sluys (2001b) for instance) herein, focus is placed only on the key points of this method.

The displacement field u of an element Ω e crossed by a strong discontinuity can be mathematically decomposed, following (Simo and Oliver (1994)), by

u = ū + (H Γ d -ϕ e ) [|u|] , (6)
where ū is a regular and continuous part of the displacement field that allows us to impose standard boundary conditions (Oliver (1996), part II), H Γ d is the Heaviside function centred on Γ d such as

H Γ d =    1 if x ∈ Ω e 0 if x ∈ Ω e , (7)
ϕ e is an arbitrary continuous shape function in Ω e of unit value at each nodes in Ω e and null at each nodes in Ω e and finally, [|u|] is a continuous function on Ω e . The components of the displacement jump at the discontinuity are given by [|u|] x∈Γ d .

By taking the symmetric gradient of equation ( 6) and following the development proposed in (Simo and Oliver (1994)), we obtain the underlying strain field such as

∇ s u = ε = ∇ s ū + (H Γ d -ϕ e )∇ s [|u|] -([|u|] ⊗ ∇ϕ e ) s bounded + δ Γ d ([|u|] ⊗ n) s unbounded , (8)
where δ Γ d is the Dirac-delta distribution centred at the surface discontinuity and n its normal vector. Except from the fact that [|u|] represents the components of the displacement jump when evaluated on Γ d , this field will not be explicitly defined. At this stage, the resulting information of this reasoning is the decomposition of ε into a regular and bounded part εb and an unbounded one εu

ε = εb + εu = (H Γ d -ϕ e )∇ s [|u|] -([|u|] ⊗ ∇ϕ e ) s εb + δ Γ d ([|u|] ⊗ n) s εu . (9)
And so the total strain field in presence of a strong discontinuity only has the form

ε = ∇ s ū + ε = ∇ s ū + εb + εu . (10)
We can notice that strong discontinuity approach brings an unbounded part into the strain field. With such properties, physical considerations like continuity of the traction vector at Γ d or bounded stress in Ω e seem rather impossible to respect. Under the name of strong discontinuity analysis, the authors in (Simo et al. (1993)) propose solutions in order to ensure that constitutive models used are still consistent in presence of strong discontinuity. Among them, the so-called Discrete Strong Discontinuity Approach (DSDA) (Oliver (2000)) is used here. Hence, the continuum model equipped with strong discontinuity leads to an underlying discrete model at the discontinuity Γ d . On the one hand, a standard continuum relationship links strain and stress fields in Ω e \Γ d and on the other hand, a traction-separation law linking the traction vector to the crack-opening [|u|] on Γ d is introduced, modeling the failure mechanism.

Remarks

• ε (i.e. Θ) is not defined yet. It will be in Section 3.

• [|u|] is chosen as a constant function thus its symmetric gradient is equal to 0. This leads to an explicit formulation of the bounded part of the strong enhanced strain referred as Kinematically Enhanced Strain (KES). It is defined by

εb = -([|u|] ⊗ ∇ϕ e ) s . (11)
• We remind that the contribution of weak and strong discontinuities can be added if both of them occur, leading to the following strain enhancement

ε = ∇ s ū + ε + ε. (12)
• Additional set of unknowns [|ε|] and [|u|] has to be solved. The resolution will be detailed in Section 4. For the time being, it is sufficient to know that each set is defined at the element level. Hence, a local resolution will be performed for their determination letting the problem size unchanged.

Discontinuity implementation into a three-field variational framework

The Finite Element implementation of the problem is cast into a three-field variational framework. It is shown in this section how such a statement allows us to blend the presented strain enhancement into a Finite Element framework. Both yet undetermined shapes of the enhanced strains are defined following both static and kinematics considerations. Moreover, this formulation leads a handful incorporation of the discrete model using traction-separation law on Γ d .

Hu-Washizu three-field variational formulation

The main idea of the Hu-Washizu (Washizu (1982)) formulation is that, unlike classical displacement formulations, the three fields (ū, ε, σ), respectively the standard displacement field, the standard strain field and the standard stress field are considered independent. We call ( η, γ, τ ) respectively the virtual displacement, the virtual strain and the virtual stress fields. With this notation at hand, the Hu-Washizu formulation reads as follows:

• the classical equilibrium equation (or virtual work) (virtual displacement)

HW ū( ū, ε, σ; η) = Ω ∇ s η : σ dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, (13)
• the kinematics equation (virtual stress)

HW σ (ū, ε, σ; τ ) = Ω τ : (∇ s ū -ε) dΩ = 0, (14)
• the behavior equation (or constitutive model) (virtual strain)

HW ε (ū, ε, σ; γ) = Ω γ : (σ (ε) -σ) dΩ = 0, ( 15
)
where σ is a stress field that verified the behavior law and Ω η • ρb dΩ and Γt η • t d∂Ω the virtual work of the volumic and surfacic loadings. Note that in classical variational formulation (in displacement) we usually have, by definition

   ε = ∇ s ū σ = σ(ε)
which would make ( 14) and ( 15) irrelevant.

Having this three-field variational formulation at hand, we now turn in the next subsection to the so-called Assumed Strain Method developed in (Simo and Rifai (1990)).

Assumed Strain Method

The first idea of the Assumed Strain Method is to enhance both standard and virtual strain fields. Each of them can be decomposed into three parts. For the standard strain field, it gives

ε = ∇ s ū regular + ε weak + ε strong enhanced (16)
and for the virtual strain field

γ = ∇ s η regular + γ weak + γ strong enhanced . ( 17
)
Note that γ respects the form of ε (see equation 9) thus it is decomposed into a bounded part γb and an unbounded part γu .

Substituting equation ( 16) into equation ( 14) and combining equations ( 13), ( 15) and ( 17) we obtain the following modified three-field variational problem:

Ω ∇ s η : σ(∇ s ū + ε + ε) dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, ( 18
) Ω τ : ε dΩ = 0, ( 19
) Ω τ : ε dΩ = 0, (20)
Ω γ : (σ (∇ s ū + ε + ε) -σ) dΩ = 0, (21)
Ω γ : (σ (∇ s ū + ε + ε) -σ) dΩ = 0. ( 22
)
Following the second idea of the Assumed Strain Method, the explicit presence of the stress field is eliminated from the modified three-field variational formulation by choosing the space of the stress field L 2 -orthogonal to the space of the enhanced strains. As a result, equations ( 19) and ( 20) are satisfied and the second term of ( 21) and ( 22) vanish. And so the modified three-field variational formulation becomes:

Ω ∇ s η : σ(∇ s ū + ε + ε) dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, ( 23
) Ωe γ : (σ (∇ s ū + ε + ε)) dΩ = 0, ( 24
) Ωe γ : (σ (∇ s ū + ε + ε)) dΩ = 0. ( 25
)
It can be noted that due to the previous L 2 -orthogonal condition, the stress field σ does not appear in the final statement of the formulation. Moreover, since no inter element (enhanced strain and stress) continuity requirement is imposed, this condition can be applied independently on each n el -typical elements. Henceforth, equations ( 24) and ( 25) are evaluated at a the element level -on Ω e for e = (1...n el ).

The third idea of the Assumed Strain Method is to assure that after having enforced the L 2 -orthogonality condition, the stress field must at least include piece-wise constant functions; the combination of these two conditions (L 2 -orthogonality and piece-wise constant functions) implies satisfaction of the patch test (see Simo and Rifai (1990)

γ dΩ = 0. ( 26
)
Regarding the weak discontinuity, equation (26.1) can be written as

Ω e γ dΩ + Ω e γ dΩ = 0. ( 27
)
By assuming that the surface of discontinuity is flat (n is constant) within an element, solving (27) gives conditions on Θ. If we also assume that both Θ and Θ are constant in Ω e and Ω e , we obtain:

V Θ + V Θ = 0. ( 28
)
The choice retained here is: Θ = V /V and Θ = -V /V where V, V and V are the volumes of Ω e , Ω e and Ω e respectively.

Regarding the strong discontinuity, equation (26.2) can be written as

Ωe γb dΩ + Ωe γu dΩ = 0. (29)
If we choose for γu a form respecting εu namely

δ Γ d ([|η|] ⊗ n) s , we get Ωe γb dΩ + Γ d ([|η|] ⊗ n) s d∂Ω = 0. ( 30
)
The same assumption of constant strain field and flat interface within an element are made.

Thus γb has the form of the so-called Assumed Enhanced Strain (EAS)

γb = - A V ([|η|] ⊗ n) s , (31)
with A the area of the discontinuity.

This method gives naturally an non-symmetric formulation: the virtual strain field space γb based on EAS (equation ( 31)) differs from the standard strain field space εb based on KES (equation ( 11)). It has been shown that this non-symmetric approach leads to superior numerical results in spite of the fact there is a variational inconsistency in the derivation (see Wells and Sluys (2001b)).

In the next section, the Finite Element discretization of equations ( 23), ( 24) and ( 25) is presented.

Finite Element discretization and resolution methodology

In this part, the Finite Element discretization of equations ( 23), ( 24) and ( 25) which correspond to the system to be solved is presented. First the discretization of the standard strain field and the virtual strain field respecting the form obtained in Sections 2 and 3 is shown, in a second time those discretizations are incorporated in equations ( 23), ( 24) and

(25) to get the Finite Element problem. Note that the next Finite Element discretizations are general and valuable for tetrahedron, beam or truss elements.

Note: Unless mentionned, Voigt notation is now used.

Strain fields discretization

Following the form obtained in Sections 2 and 3, the discretization of the standard strain field and the virtual strain field is written as follows:

ε = [ε xx ε yy ε zz 2ε xy 2ε yz 2ε xz ] T = B d ∇ s ū + G w [|ε|] ε + G s [|u|] εb , ( 32a
) γ = [γ xx γ yy γ zz 2γ xy 2γ yz 2γ xz ] T = B δ ∇ s η + G w [|γ|] γ + G * s [|η|] γ , (32b)
where several interpolation matrices introduced correspond to: B (= ∂N ) the standard strain interpolation matrix, G w the standard and virtual fields corresponding to the weak discontinuity, G s the standard field corresponding to the strong discontinuity (bounded part εb ) and G * s the virtual field corresponding to the strong discontinuity (bounded and unbounded part γ). d are nodal displacements and [|ε|] (resp. [|u|]) corresponds to additional unknowns arising from the kinematics enhancement of weak (resp. strong) discontinuity.

Regarding weak discontinuity kinematics, G w can be decomposed into Θ and a constant part H w that only carries information on the interface vector n such as:

G w =    G w = Θ H w = V V H w in Ω e G w = Θ H w = -V V H w in Ω e . ( 33
)
The KES interpolation matrix G s derives from equation ( 11) and can be computed by explicating the arbitrary function ϕ e . Such function can be defined using standard interpolation shape functions as follows:

ϕ e (x) = nen a=1 N a p a with p a =    1 if node number a ∈ Ω e 0 if node number a ∈ Ω e , (34)
where n en is the number of nodes in the element and p a the nodal values of ϕ e . G s is therefore the equivalent symmetric operator (• ⊗ ∇(ϕ e )) s in the Voigt notation of equation (32a).

The EAS interpolation matrix G * s is known from equation is decomposed in a bounded G * s,b and an unbounded G * s,u part that derives from equation ( 31) and (9), respectively. It leads to:

G * s = G * s,b + G * s,u = - A V + δ Γ d H * s , (35)
where H * s is the equivalent symmetric operator (• ⊗ n) s in Voigt notations of the strain field as in (32b).

Having those discretizations and the form of G w , G s and G * s at hand, equations ( 23), ( 24) and ( 25) can now be discretized.

Finite Element discretization of the problem

Injecting equations (32) into ( 23), ( 24), (25) and considering (33) and ( 35), the discretized problem to be solved reads as follows:

n el A e=1 [f e int -f e ext ] = 0, (36a)
h [|ε|] = V V Ω e H T w σ dΩ - V V Ω e H T w σ dΩ = 0 ∀e ∈ [1 .. n el ], (36b)
h [|u|] = - A V Ω e H * ,T s σ dΩ + Ω e H * ,T s σ dΩ + Γ d T d∂Ω = 0 ∀e ∈ [1 .. n el ], (36c)
where all the σ linearly depend on d, [|ε|] and [|u|] and where

f e int = Ω e B T σ dΩ + Ω e B T σ dΩ, (37a)
f e ext = Ωe N T ρb dΩ - Γt N T t d∂Ω and (37b)
T = H * ,T s σ == σ • n Matrix is the traction vector. ( 37c
)
The behavior law written in relation with the spatial position x is As the whole framework fits into the Discrete Strong Discontinuity Approach, two behaviors have to be considered. First, regarding the continuum bodies Ω e and Ω e , the behavior is considered elastic. This rather strong assumption is justified by the general spirit of simple meso-scale modeling. However, more complex behaviors such as plasticity or damage can also be implemented (see Oliver (1996) for details). Then, the second part of the modeling takes place at the discontinuity surface Γ d . A governing law that links the traction vector

σ(d, [|ε|] , [|u|]) =    σ = C ε = C Bd + G w [|ε|] + G s [|u|] if x ∈ Ω e σ = C ε = C Bd + G w [|ε|] + G s [|u|] if x ∈ Ω e , (38)
T = σ| Γ d
• n 1 is defined on the surface Γ d where the crack opening magnitude drives the non-linear failure mechanism. Such laws are often referred as traction-separation laws. In order to fit in the general spirit of simple meso-scale modeling, a very basic single traction criterion and brittle softening is used here. The next subsection describes its main characteristics and shows how it is incorporated in the Finite Element problem (36a, 36b, 36c).

Finally, after linearising these equations, the Finite Element problem is written under a matricial form and the solving strategy is presented.

Discrete constitutive model on Γ d

The discrete model is based upon a relationship between the traction vector T (via an equivalent stress σ eq ) and the crack opening magnitude labelled [u]. The equivalent stress σ eq value has to be tackled with careful attention in order to depict a meaningful representation.

Thus two cases have to be considered. First, if we consider the case where only the strong discontinuity is present (case 1) of Section 1), no geometrical information are given on the crack orientation. The choice is made here to use the larger principal stress component as the equivalent stress σ eq = σ I . This principal component is simply the first eigenvalue of the stress matrix. Furthermore, the corresponding eigenvector n I represents its direction.

The physical meaning of eigenvalue problems naturally leads to the choice of this vector in order to represent the crack orientation n ← n localization

I

. Secondly, if both weak and strong discontinuities are present (case 3) of Section1), it is assumed that crack opens at the interface between each phases. The main difference with the latter case is that n is defined by geometrical characteristics and therefore is independent from any stress state.

Hence, traction vector can be defined prior to localization. In this case, its projection on the interface direction is used such as σ eq = n • T . Attention is drawn to the fact that in this case, the shear components are non-zero but we assume they do not participate to the localization process. This is a major assumption that could be improved by considering more complex criteria. However, in the spirit of multiscale analysis we aimed at first considering the most simple case.

The discrete model is introduced when the equivalent stress σ eq reaches a certain amount σ y , which has to be seen as a meso-scale material characteristic. Mathematically speaking, the strong discontinuity is introduced when the following localization criterion becomes zero:

Φ l = σ eq -σ y . (39)
Then, the failure mechanism is driven by a scalar opening criterion Φ o defined as follows:

Φ o = σ eq -(σ y -q), (40)
where

q = σ y 1 -exp - σ y G f [u] . ( 41
)
It can be seen that a second material parameter G f called the fracture energy that governs the amount of energy necessary for the complete material failure is introduced in (41). Since [u] represents the magnitude of the crack opening [|u|], the criterion is single valued. However, an additional projection of the crack has to be defined in order to fit in the previous framework, introducing a unit jump vector n p :

[|u|] = [u] n p . (42)
Finally, the behavior can be split into two parts: an elastic one in the body Ω e (linear relationship between strain and stress fields out of the discontinuity Γ d ) as drawn on figure 4(a) for a 1D case and a discrete one on the discontinuity Γ d by means of a traction-separation law that links Ω e and Ω e as shown on figure 4(b). Since several choices retained here depend on element kinematics, details on the different ingredients of this model (σ eq , n p . . . ) will be given later.

From equation (36c) and making the assumption of flat interface and constant stress, an easy integration gives an explicit expression of T as a average value of σ and σ weighted by volumes such as:

T = 1 V H * ,T s (V σ + V σ ), (43)
where V, V and V are the volumes of Ω e , Ω e and Ω e , respectively. Considering now (38), equation ( 43) can be written in function of (d, [|ε|] , [u]) such as:

T = 1 V H * ,T s V C + V C B K s * b d (44) + V V V 2 H * ,T s C -C H w K s * w [|ε|] ( 4 5 ) + 1 V H * ,T s V C + V C G s n p K s * s [u]. (46)
Following (Ibrahimbegovic et al. (1998)), the criterion Φ o is now incorporated in the system (36a, 36b, 36c

∆Φ o = ∂σ eq ∂T ∆T + ∂q ∂[u] ∆[u] = ∂σ eq ∂T 1 V H * ,T s V C + V C B K s * b ∆d + ∂σ eq ∂T V V V 2 H * ,T s C -C H w K s * w ∆ [|ε|] + ∂σ eq ∂T 1 V H * ,T s V C + V C G s n p K s * s ∆[u] + σ 2 y G f e -σy[u]/G f Kq ∆[u].
Thus the linearisation of Φ o = 0 gives

K s * b (k) n+1 ∆d (k+1) n+1 + K s * w (k) n+1 ∆ [|ε|] (k+1) n+1 + (K s * s + K q ) (k) n+1 ∆[u] (k+1) n+1 = -Φ o (k) n+1 . (47)
In the next subsection, we present the form of the final problem to be solved under a matricial form and the resolution strategy.

Global system and resolution strategy

Because the linearisation of equations (36a) and (36b) is trivial the mathematical development will be skipped. The global system to be solved in terms of increments of d, [|ε|] and [u] and condensed in a matricial format is

     K bb K bw K bs K wb K ww K ws K s * b K s * w K s * s + K q      (k) n+1      ∆d ∆ [|ε|] ∆[u]      (k+1) n+1 =      - n el A e=1 {f e int -f e ext } -h [|ε|] -Φ o      (k) n+1 , (48)
where

K bb = B T V C + V C B, K bw = V V V B T C -C H w , K bs = B T V C + V C G s n p , K wb = V V V H T w C -C B, K ww = V V V 2 H T w V C + V C H w , K ws = V V V H T w C -C G s n p , K s * b = ∂σeq ∂T 1 V H * T s V C + V C B, K s * w = ∂σeq ∂T V V V 2 H * T s C -C H w , K s * s = ∂σeq ∂T 1 V H * T s V C + V C G s n p , K q = σ 2 y G f e -σy [u] / G f . ( 49
)
Solving system ( 48) is done at two levels. First, following the operator split method spirit,

variables (∆ [|ε|] , ∆[u]
) are determined at the element level (local solving) for a given ∆d by solving: . ∆d is calculated by solving the global equilibrium equation using a static condensation (Wilson (1974)) on the local (known) variables (∆ [|ε|] , ∆[u]). This leads to a matricial system, with a modified stiffness matrix K sc , to be solved such as

   h [|ε|] = 0 Φ o = 0 . ( 50
K sc (k) n+1 ∆d (k+1) n+1 = - n el A e=1 {f e int -f e ext } (k) n+1 , (51)
where

K sc (k) n+1 = K bb -K bw K bs      K ww K ws K s * w K s * s + Kq   (k) n+1    -1   K wb K s * b   (k) n+1 .
Even though the stiffness matrix has been changed due to kinematics enhancement, both its size and sparsity are unchanged. Hence, no matter how many heterogeneities are represented or how many elements have starting to fail, the global size of the problem is preserved. In terms of numerical resources, the memory needed only depends on the mesh size (number of nodes). Naturally, local Newton algorithms slows down the global calculation as the number of strong discontinuity activated increases. By using the static condensation, a standard FE problem is retrieved, where increments of d have to be found in order to respect the global equilibrium equation 51.

The most common method used to solve those problems in case of non linearity are the so-called Newton methods. However, it requires full calculation of full stiffness matrix at each iteration, and since morphological modeling requires rather fine meshes, a quasi-Newton algorithm coupled with an iterative solver is used here. Among the huge diversity of those algorithms, the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm is retained. It makes the number of arithmetical operations to fall down from O(n 3 ) to O(n 2 ). The price to pay is that quasi-linear convergence is obtained (instead of quadratic). Full details are given in the original papers (Broyden (1970a,b); Fletcher (1970); Goldfarb (1970); Shanno (1970)).

Moreover numerical implementation details are in (Matthies and Strang (1979)). Added to the quasi-Newton BFGS, a line-search method is also used, modulating the incremental displacement norm for each iterations by a factor s: k+1) .

û(k+1) = û(k) + s (k+1) ∆d (
See Dahlquist (2003) for details on the computation of s. It is recalled that powerful algorithm adapted to the E-FEM such as (Oliver et al. (2008)) and not implemented here can significantly increase the computation performance.

In the next section, we move to applications of this model in the case of 4-node tetrahedron elements. First the explicit form of the matrices encountered all along this paper is presented for sake of clarity. Then the criterion Φ o is written in the considered case.

Application to 4-node tetrahedron elements

The use of 4-node tetrahedron elements is now presented in the case of this model. This volumic discretization brings several advantages compared to the initial implementation of this framework based on lattice discretization (Benkemoun et al. (2010)). First, an exact representation is made regarding the volumic tessellation of the mesh, leading to exact representation of constant stress problem. Therefore, the problematic of mesh convergence presented in the aforementioned paper with lattices is irrelevant in this context. Furthermore, a complete kinematics of cracks opening can now be represented (mode I, II). Hence, both geometrical construction and mechanical behavior are more accurately depicted.

Interpolation matrices

In this part, attention is drawn to the explicit definition of the matrices encountered through this paper in the case of tetrahedron elements. In addition, for sake of convenience, fields are represented in their Voigt notation. Hence, each matrices are developed in this format.

In the case of an element cut by a discontinuity Γ d , the tetrahedron is split into two sub-domains Ω e and Ω e delimited by the interface Γ d of direction vector n (see figure ( 5)).

It is reminded that this surface is assumed to be flat (n is constant over Γ d ). Numerical implementation of such geometrical construction is not trivial and several cases have be taken into account depending on the surface orientation. Indeed, Ω e and Ω e can be polyhedron of respectively 6 and 4, 5 and 5 or 4 and 6 nodes. And since the formulation involves only volumes V and V (not the discontinuity area, see equation ( 43)), it has to be considered with utmost attention.

n Dealing with weak discontinuity, it is reminded that the interpolation matrix G w can be decomposed into Θ, a piece-wise constant function depending on the considered subdomain and H w a matrix containing information on the discontinuity surface orientation (see equation ( 33)). Following the same Voigt convention, H w can be constructed so that the vector format of the strain enhancement matches its tensor definition in equation ( 33).

Ω e Γ d m t Ω e
The vector n is written in the global coordinate system n = [n x n y n z ] T , leading to the following interpolation matrix format:

H w =               n 2 x n x m x n x t x n 2 y n y m y n y t y n 2 z n z m z n z t z 2n x n y n x m y + n y m x n x t y + n y t x 2n y n z n y m z + n z m y n y t z + n z t y 2n x n z n x m z + n z m x n x t z + n z t x               . ( 52
) G w =    G w = Θ H w = V V H w in Ω e G w = Θ H w = -V V H w in Ω e . ( 53
)
Regarding strong discontinuity, it is reminded that two matrices have to be considered, namely G s and G * s,b , making the formulation non-symmetric. The form of G s is given by

G s = -               ∂ϕe ∂x 0
0 ∂ϕe ∂x               (54)
and the form of

G * s,b by G * s,b = - A V H * s with H * s =               n x 0 0 0 n y 0 0 0 n z n y n x 0 0 n z n y n z 0 n x               . ( 55
)
All components of the enhanced strain field are now explicitly defined. The standard strain field can thus be built respecting the form of equation (32a) and the system (48) solved.

We now turn explanation on the localization and opening criterion.

Localization and opening criterion

Since the discontinuity surface orientation is constructed with geometrical properties for interface element and with stress consideration otherwise, two cases have to be treated separately.

In the case of weak discontinuity within the element, the interface is always defined.

Hence n is known prior to any mechanical calculation. Thus, the traction vector T can be defined before the localization. In order to represent the interface orientation as the weakest direction, the equivalent stress for localization is defined as the projection of the traction vector on it:

Φ l = σ eq -σ y = n • T -σ y = T n -σ y (56)
with

T = 1 V H * ,T s (V σ + V σ ). ( 57
)
On the contrary, if no material discontinuity can define an interface, strong discontinuity appears with stress state consideration. One can note that in this case, a constant stress tensor can be given for the whole element since it is supposed free of material discontinuity.

Its orientation is defined by the principal direction of the stress tensor. If σ I is its eigenvalue then:

Φ l = σ eq -σ y = σ I -σ y . (58)
When localization occurs (Φ l = 0), the corresponding eigenvector n I is recorded and set as the interface orientation: n ← n localization

I

. It is assumed that its value remains constant through time. Afterwards, the traction vector is defined by:

T = H * ,T s σ (59)
which follows its previous definition with σ = σ = σ.

In both cases, after localization the discontinuity surface and its orientation n are defined.

In order to model the same failure mechanism whether an interface element is considered or not, the opening criterion Φ o is assumed to be identical. The equivalent stress is taken to be the projection of the traction vector on n:

Φ o = n • T -(σ y -q). (60)
Furthermore, this definition of the equivalent stress leads to a very simple written expression of the equivalent stress derivative:

∂σ eq ∂T = n. (61)
Finally, by defining the projection vector n p of equation ( 42) as the normal vector n, a mode I opening mechanism is represented. Attention is draw to the fact that, in this case, the evolution of the shear tractions are not driven by the displacement jump. That is a major assumption which represents a flaw in the model. Improvements may come from more complex opening mechanisms such as mode II or mode I+II. Those are both under development and related issues are discussed in the conclusion of the paper.

At this stage, the whole system ( 48) is now explicitly known in the case of tetrahedron elements. Thus we can move to numerical examples in the context of concrete-like material to show the features of the model developed through this paper.

6. Numerical analysis of concrete: from meso to macro-scale

In this section, attention is focused on the capability of the enhanced FE model to represent the main features of concrete by means of mechanical loading. Starting from mesoscale (millimeter scale), uniaxial tension and compression tests are performed in order to

show, on the one hand, the emergence of the typical asymmetry of the respective macroscopic responses and on the other hand, a complete analysis of the crack patterns. Moreover, non proportional loadings are carried out in order to measure the induced anisotropy. This study is based on the previous uniaxial tension and compression tests by means of damage indicator (i.e. linear post-analysis on Young modulus degradation) as well as residual strength (i.e. non linear post-analysis on tensile strength decrease).

Each numerical example given in this chapter is made on a heterogeneous 100 × 100 × 100 mm 3 specimen for which two phases are modeled. Based on unions of excursion sets (Adler (2008), Roubin et al. (2014)), this method explicitly represents aggregates of different sizes melt within a matrix that is, roughly speaking a mortar. The former are modeled with three average diameters of 15, 7 and 5mm representing respectively 25, 50 and 25 % of the total 30 % volume fraction (see figure 6(b)). Random shaped inclusions are yielded by thresholding correlated Random Fields (see Roubin (2013)) and performing the union.

Following the non-adapted mesh spirit, once projected onto the FE discretized space, this morphology is represented by two kinds of elements: those which are completely included within the matrix or an aggregate and those, close to an interface, that are split and thus enhanced by a weak discontinuity. Table 1 summarizes the material properties at the mesoscopic scale. It is reminded that the model contains, for each phase, two elastic and two failure parameters, the Young modulus E, the Poisson ratio ν, the yield stress σ y and the fracture energy G f , respectively.

Furthermore, a yield stress and a fracture energy are also set to define the interface.

It is worth noting that:

• Aggregates are assumed to remain elastic. Thus, no crack may initiate or propagate within those elements.

• Weakly enhanced elements do not require any specific characteristics in the elastic regime. Indeed, as presented in the previous chapter, prior to any debonding (modeled through a strong discontinuity), they model a perfect interface, thus with infinite rigidity.

• For sake of simplicity, the failure mechanisms are taken to be equal for both matrix and interfaces.

However, the difference in elastic properties of the matrix and aggregates leads to stress concentrations at the interface. Since the failure criterion is based on stress considerations, a similar yield stress in both matrix and interface still leads to a weaker behavior of the latter.

Finally, the computations are performed under uniform displacement control along the first spatial axis (X-direction). Moreover two other faces of the domain (normal to the Y and the Z direction) have constant zero value for their normal displacement. Hereafter, a value corresponding to the X-direction is referred as axial whereas any transversal quantities refers to an average value set up on both the Y and Z-directions.

Phase E [GPa] ν [-] σ y [MPa] G f [J.m -2 ] Aggregates 100 0.2 - - Mortar 20 0.2 9 0.1 Interface - - 9 0.1
Table 1: Meso-scale material characteristics of each phase and interface 6.1. Analysis of the macroscopic responses for simple traction and compression 6.1.1. Axial upscaled properties and crack pattern analysis First it is worth noting that, although the failure criterion at the meso-scale is triggered in tension only, macroscopic simple compression leads to the failure of the specimen. Indeed this feature is clearly a consequence of the structural effect that is set up by the explicit representation of the aggregates.

Second, it clearly appears that the macroscopic failures are unsymmetric. This feature is typical of quasi brittle materials such as concrete. Table 2 sums up macroscopic key values extracted from figure 7: we denote by E M the macroscopic modulus, σ M f the macroscopic tensile or compressive strengths, ε M f the corresponding macroscopic failure strain and finally D p the total dissipated energy. The last value is computed through integration of the macroscopic axial force over the axial displacement.. As expected the initial elastic behavior is strictly symmetric. Even though it is quite difficult to observe the end of the purely elastic region on the plots, it can safely be assumed that, for compression, the transition to the non linear behavior occurs for a more important macroscopic stress than for tension. This unsymmetric elastic domain is the first feature

Loading path E M [GPa] ν M [-] σ M f [MPa] ε M f [-] D p [J]
Compression 37.8 0.195 39.2 1.2 10 -3 50 Tension 37.8 0.195 3.9 1.2 10 -4 0.5

Table 2: Macroscopic upscaled material properties for both tension and compression that emerges from the multi-scale analysis.

Regarding the energy D p needed to reach the specimen ruin, it can be noted that both are greater than the mesoscopic fracture energy G f (in J.m -2 ) assigned to the mortar (corresponding values of D p in tension and compression are 50 and 5 000 J.m -2 , respectively).

Furthermore, the fact that this energy is significantly greater in compression reflects the more brittle behavior of concrete when tested in tension. It is naturally linked with the asymmetric strength values and their corresponding failure strains for which the ratio of compression to tension are both 10. A discussion on this meaningful result is given in the main conclusion of this paper. . This corresponds to the elements for which a strong discontinuity has been introduced. It can be seen that those crack patterns are very significant either in tension or compression. First in tension, some micro-cracks are linked in order to set up a single macro-crack that roughly lies in the transversale plane (best seen on Figure 9) . Obviously this macro-crack is tortuous and goes around the aggregates that remain elastic. Second, in compression, it can be observed that several macro-cracks are present (contrary to tension) and that they are roughly parallel to the axial direction.

Aggregates properties

E (-40%) ν (+50%) Compression strength [MPa] 48.2 (+23%) 40.3 (+3%) Tension strength [MPa] 4.0 (constant) 3.72 (-5%)
Compression to tension ratio 12.05 (+20%) 10.8 (+8%) Table 3 shows the influence of the elastic properties of the aggregates on both compression and tension macroscopic strengths. As expected, the ratio between those two quantities is largely influenced by the contrast of Young's modulus. To be more precise, the less this contrast is, the more the compressive strength increases and so the ratio (the tension strength remains almost constant). On the contrary, the Poisson ratio of the aggregates seems to have a small influence on the macroscopic strengths: from ν = 0.2 to ν = 0.3 the compressive strength shows a 3% increase and the tensile a 5% decrease.

On a more general point of view, the question of the geometrical representation of a macroscopic crack using a local method as the Embedded Finite Element is non trivial. For example, in the two-dimensional case, using constant triangular elements, it has been shown (see Jirásek (2012) for details) that the best way to produce a suitable crack trajectory and avoid numerical issues such as stress locking is to combine two methods. First a non-local formulation of the smeared crack approach, giving crack orientation in each element. Then a tracking algorithm to enforce the crack path continuity between each element. The major drawback of this implementation is that the local spirit of the E-FEM (directly inherited from the FEM itself) is lost. Indeed, in addition to non local damage, path continuity enforcement implies, for an element, a crack position that depends of those of its neighbors. Moreover, in the three-dimensional case, continuity of flat plane (crack) is often impossible.

Herein, the E-FEM implementation has to be placed within the multi-scale context. In this case, a single fractured element is not considered to be representative of any specific macroscopic feature. However, it is only when a large number of those activated elements are merging that it may be considered that they model a continuous path at the macroscopic scale. It is for these reasons that, herein, no specific effort has been made in order to enforce any path continuity at the meso-scale. It is reminded that herein, the discontinuity is considered piecewise constant in each element. In the case of weakly discontinuities (interfaces elements), the orientation and position of the strong discontinuity are predefined by the crossing heterogeneity. However, in the case of standard kinematics (matrix element), the orientation of the strong discontinuity is set to be the direction of the larger principal stress at the localization time. In this case it is assumed that the discontinuity path through the centroid of each tetrahedron. This assumption impacts only the construction of the functions ϕ e .

Finally, features like multi-cracking or branching, which usually require a complex local numerical implementation (within an element), are herein omitted at the meso-scale. However, as shown on Figure 11 (which is a zoom made on a subset of Figure 10(b)), it can be retrieved at the macro-scale. This picture shows a crack that splits in two branches; a main branch (on the top) with larger opening values and a second branch that eventually vanishes. Generally, these branchings come from an aggregate "blocking" the way of the crack propagation direction.

Transversal strain analysis

Regarding the transversal behavior of the specimen, several observations can be made either in tension or in compression. It is recalled that the axial direction X corresponds to the imposed displacement direction and transversal values are defined as the average of the values along Y and Z. The results presented here are still based on the same one-dimensional macroscopic tests (tension/compression) mentioned above. The macroscopic Poisson ratio can be determined using the transversal strains by:

ν M = - ε M tr ε M ax (62)
Herein, this property is extended to the diffuse cracking regime in order to illustrate the mechanism that leads to the specimen loss of rigidity -see figure 12 where it is plotted in terms of axial strain for both tension and compression tests. First, the elastic part shows that the macroscopic Poisson ratio is of the same order as for the meso-scale: ν M = 0.195 in both cases. Afterwards those values are diverging. First, in tension, the ignition of diffuse cracking causes local strain to release and thus making the macroscopic strain decrease with the Poisson ratio. On the contrary, in compression, this local strain releasing causes a heightening of the transversal mechanism leading to a significant increase of the apparent Poisson ratio.

Since the meaning of the Poisson ratio is highly contestable with strongly non-linear failure behavior, the post-localization analysis is only based on transversal strains. For that matter, Figure 13(b) and Figure 13(a) show the macroscopic response up to the specimen ruin in terms of the axial strain ε M ax (solid curve) and of the transversal direction ε M tr (dashed curve). Regarding the tension test, during the post-peak phase the transversal strain decreases and tends to vanish (see Figure 13(b)). It represents the unloading that occurs in the specimen -at a macroscopic scale -after the main crack localization. In contrast, during a compression test, the transversal strain still increases after the peak load (see Figure 13 The dilatancy δ of the specimen can also be computed by considering the trace of the macroscopic strain tensor. Thus, the relative variation of the volume, drawn on Figure 14, is:

δ = ∆V V 0 = ε M ax + ε M tr Y + ε M tr Z (63)
Notice that the same magnitude in the dilatancy rate is observed in tension and in compression regarding the elastic region. This result is in conformity with the identical Poisson ratio value. Naturally a tension test produces a volumetric expansion (δ > 0)

while compression first produces a contraction (δ < 0). However, the cracking process increases the dilatation rate for both tests, which is a major feature of many materials such as concrete. Hence, when in tension the slope increases in the diffuse cracking regime, for compression the contraction slows down. The localization is characterized by: in tension a sudden increase in rate and in compression a maximum (considering absolute value) of the dilatancy. Afterwards, the dilatancy sign changes in the post-localization region, which is an interesting feature, also experimentally observed for concrete. In (Torrenti (1987)) the interpretation is made that, in compression, the contraction corresponds to a predominant elastic effect where the following expansion reveals an important cracking stage. However, experiences revealed a volumetric strain which switches sign prior to the localization.

Dissipated energy

By integrating the macroscopic force over the displacement, the total energy can be calculated from the macroscopic response. It is shown on Figure 7 for each calculation step.

Furthermore, considering a fictitious elastic unloading up to a zero force level, a difference can be made between the elastic and the dissipated part. Thus it gives an additional way to understand the failure mechanisms and to compare tension and compression behaviors. in compression can be seen. The fall of the elastic energy to a nearly null value represents the unloading outside the macro cracks region and the fact that all the energy is dissipated close to the macroscopic cracks.

Numerical evaluation of the induced anisotropy

After having shown the macroscopic responses of the model in relation with simple radial and proportional macroscopic loading paths, focus is placed on non proportional cases. Here, our objective is to show the emergence, at macro-scale, of some features related to any anisotropic behaviors. This anisotropy shall be revealed by the non proportionality of the loading path.

Considering a tension or a compression test with monotonic loading, the failure mechanism that leads to the specimen ruin induces a strong anisotropy of cracks pattern. On the basis of this simple observation, an analysis of the macroscopic material propertiese.g. Young modulus, tensile strength, etc. -for each step of the previous monotonic loading tests (in both tension and compression) is now proposed. Hereafter, the two first parts of the macroscopic loading path (tension or compression) are referred as principal calculations. In order to yield residual material properties, additional calculations, that inherit from the principal, are performed. They are referred as secondary computations. Basically, the inheritance from principal to secondary calculation is made through the non linear data, i.e. the whole set of meso-scale cracks (with their orientations and opening values). This second part of the numerical analysis is also made under displacement control. Actually the displacements obtained from the first part are imposed and additional displacements are added. Those displacements correspond to the same kind of boundary conditions along a different direction (Y or Z).

Anisotropic induced damage

Here the residual property of concern is the damage Young modulus, which is computed as the secant modulus on the macroscopic stress -strain curve. In order to catch the anisotropy, from each step of the principal calculation three secondary calculations are performed in the three directions. It leads to three macroscopic secant moduli: an axial ẼM ax value and two transversal ones ẼM tr Y and ẼM tr Z . The results are displayed by defining "damage" variables d ax , d tr Y and d tr Z , respectively. They are built to compare the upscaled residual secant moduli to those corresponding to the initial state E M (which are the same for all direction):

d ax = E M -ẼM ax E M , d tr Y = E M -ẼM tr Y E M and d tr Z = E M -ẼM tr Z E M (64)
Theoretically, in the elastic region of the principal test, those variables are null and then tend to increase along with the specimen failure state. Herein, the secant moduli upscaled values ẼM ax and ẼM tr yield the damage variables d ax and d tr (equation 64). Results for tension are given on Figure 16(b) where the two damage variables are plotted in terms of the macroscopic axial strain of the principal calculation. It can be observed that the axial damage d ax is growing faster than that of the transversal one d tr . Moreover, d ax reaches a value of ≈ 0.85 which corresponds to a highly damaged state when d tr hardly reaches 0.15. Hence, the elastic property is far more spoiled in the axial direction. This result reflects the characteristic morphology of the cracks pattern, splitting the specimen in two by a plane roughly perpendicular to the axial axis. As the macroscopic crack grows, the "link" between each part of the specimen becomes weaker, leading to a decrease of the upscaled secant modulus in this direction. On the contrary, it remains several non-broken paths on the transversal directions that give to the specimen a higher rigidity. Regarding the compression test, Figure 16(a) shows the opposite effect. Indeed, here the transversal damage is more important than the axial one. It may be explained by cracks patterns that form planes parallel to the axial direction, leading to a higher loss of rigidity in the transversal directions. Furthermore, the more diffuse aspect of the crack repartition makes the difference between axial and transversal damage less important.

Finally, for both tension and compression cases, it can be noticed that the two transversal damages are of the same order of magnitude, representing isotropic behavior in these two directions. Somehow, it can be said that the macroscopic elastic behavior is shifting from an isotropic case to a transverse isotropic one. More numerical investigations for the second step of those non proportional loading paths may determine the complete elasticity tensor for this case.

Induced anisotropy for tensile strengths

Still dealing with both tension and compression tests, interest is now taken in the residual tensile strength of the specimen for each direction. As the previous section computed an elastic property, herein, non-linear calculations are performed at each step of the principal test in order to yield these failure properties. The results are plotted in terms of residual strengths defined by the ratio between actual tensile strengths f ax , f tr Y and f tr Z to the initial state one f : However, regarding the transversal directions, it can be noticed that this state seems rather more deteriorated than what experimental results shows.

r ax = f ax f , r tr Y = f tr Y f and r tr Z = f tr Z f (

Discussions and conclusions

This paper first sets a numerical tool for the modeling of the quasi-brittle behavior of heterogeneous materials such as the cementitious ones. This model can be viewed as a multiscale model, aiming at upscaling data from mesoscopic scale (millimeter scale for concrete) to macro-scale.

At meso-scale the numerical implementation is based on the introduction of kinematics enhancements of two kinds within the FE context. The first enhancement, referred to as weak, leads to a non-adapted mesh strategy for heterogeneous morphologies. It can therefore be "simply" projected onto an unstructured mesh, freeing us of any expensive algorithm that aim to match a mesh onto a given morphology. Indeed, according to the nonadapted mesh point of view, this mesh is created regardless of any physical surface (i.e. the interface between each component). The second enhancement, referred to as strong, models discontinuities in the displacement field that are viewed as micro-cracks. Their opening evolution are directly linked to the local fracture energy and tensile strength leading to a model with only two meaningful non-linear parameters. These two enhancements make the model a relevant and efficient tool to represent the failure mechanisms in a continuum region (mortar matrix cracking) as well as on its interfaces (debonding).

Numerical applications of the model have been shown by modeling concrete specimens mainly under uniaxial loadings. The massive use of E-FEM in this case, which consists in the introduction of a high number of strongly enhanced elements (typically more than hundreds of thousands), has lead to complex macroscopic crack patterns. Moreover it has revealed some emergent macroscopic responses that exhibit several features such as asymmetry of the tension-compression stress-strain relationship, which are typical of concrete. However, even though the implementation of full kinematics presented here is a major improvement compared to simpler FE (see Benkemoun et al. (2010)), several weak points are worth noticing regarding the macroscopic responses. They are now discussed.

First, a low ratio of compression to tension strengths and a too brittle behavior in compression are still observed. To the authors point of view, the main reason for that is the choice of rather simple failure mechanisms and cracks opening mode at meso-scale. Hence future work shall include for example, mode II or I+II crack opening. Furthermore, from the authors experience, it has been noticed that high fraction volumes leads to higher ratio of compression to tension strengths. Unfortunately, increasing the fraction volume leads to the need of representing thinner heterogeneity (and therefore finer meshes), thus leading to unreasonable time consuming computations (a computation as presented in this paper takes about 24h of CPU time). Second, as it can be seen in Figure 7, some locking effect can be observed for both tension and compression curves, which exhibit, for large strain, an asymptotic behavior that does not correspond to zero macroscopic force . As a matter of fact, this issue can directly be linked with the two problematics addressed just above. Indeed, it has been said that the orientation of the strong discontinuity at the interface are set to be equal to the weak discontinuities, thus being pre-defined by the heterogeneities (and not based on the principal stress directions). And since the failure is governed by a mode I opening mechanism, stress locking often occurs when the displacement incompatibilities results in shear loading. Furthermore, the ratio of weakly enhanced element increase along with the fraction volume which highlight even more the need of mode II opening. Finally, if not sufficiently discretized by the FE mesh, the morphology can exhibits "elastic paths" within the aggregates, thus explaining partially the stressed state at the end of the loading. Due to the inner complexity of this topological quantification, only a visual validation has been made, thus leading to reasonable doubts as to the disconnected state of all discretized inclusions.

Moreover the model clearly offers a large number of perspectives dealing with multiphysics. This is a major point associated to the durability of concrete structures, for which mass transfers are the cornerstone. Thus improving the physics of the modeling can be made:

(1) through a coupling between the cracks pattern (obtained from a mechanical analysis) to simple flows (such as the Poiseuille flow between two planes). ( 2) according to a better morphological representation of the heterogeneities, e.g. accounting for the largest porosity that can be found within mortar matrix. On a more general point of view, the authors think that the use of the E-FEM in the context of heterogeneous materials shall be applied, in a near future, to smaller scales (typically the micrometer scale dealing with concrete), allowing for a better understanding of the role played by simple failure mechanisms and their propagation to macro-scale. Thus the future development of this model lies in retaining the spirit of the physical significance of fine scales.

Introduction

Fracture in brittle or quasi-brittle materials -such as geological media like clay or cementitious ones -may occur under several loading paths such as tension or compression.

Modeling this phenomenon is still a very active area. The so-called macroscopic approaches are mainly based on the definition of a Representative Volume Element (RVE) which allows for defining macroscopic stress and strain tensors. Usually those quantities lead to deriving rate-independent strain-softening plasticity models (Prévost (1975)) or damage models (Ollivier et al. (2012)). Those models are leading to ill-posedness partial differential equations (Sandler and Wright (1984)) and mesh-dependent results that require techniques introducing intrinsic length scales such as the rate-dependent plasticity (Needleman (1988)) or the non-local approach (Pijaudier-Cabot and Bažant (1987)). In opposition, the approach of interest in this paper, namely the strong discontinuity approach, yields mesh-independent finite element solutions without introducing any intrinsic length scale. The main reason for this is that dissipation is computed over a surface of zero measure, thus being independent from the element size (Simo et al. (1993)). Moreover an important feature of the strong discontinuity approach is that the magnitude of the displacement jump may be viewed as the crack opening value, which thus belongs to the set of unknowns. This is a key point dealing with Civil Engineering structures for which the lifespan is evaluated not only according to the mechanical point of view but mainly through its durability and the corresponding mass transfers (Jourdain et al. (2011)). Considering the literature related to the strong discontinuity approach, it is worth noting that it is mainly used in the context of macroscopic problems (e.g. uniaxial test (Oliver (1996)) or bi-axial test (Wells and Sluys (2001a)). At this scale, only a few macroscopic cracks are to be modeled. Thus the ratio between the number of elements for which strong discontinuities are introduced to the total number of elements in the mesh is quite low. Moreover, because at this scale a relevant crack path modeling is of crucial importance, some authors enforce its continuity by means of numerical techniques such as tracking algorithms (Oliver et al. (2004)) or by solving secondary problems based on a crack path field (Dias (2012)). This implementations often lead to complex non-local computations, thus weakening the attractive local aspect of embedded discontinuity approaches. However, these more accurate modelings of the crack path physical representativity are known to significantly reduce issues related to high strain localization such as stress locking and mesh bias dependency.

In this paper, following the idea developed in (Benkemoun et al. (2010)) and (Benkemoun et al. (2012)), a new paradigm is presented: the use of the strong discontinuity approach at the mesoscopic scale. In this context, the strong discontinuity approach becomes a tool to model (1) a diffuse process of cracking and (2) the coalescence of micro-cracks leading to one or several macro-cracks and eventually to the failure of the specimen. This massive use of strong discontinuities leads to a large ratio of the number of elements for which strong discontinuities are activated to the total number of elements in the mesh and so to a quite complex and precise shape of the macroscopic crack(s). Moreover the large number of micro-cracks represented also offers a way to model macro-scale crack path continuity and branching as a result of micro-cracks coalescence.

At the mesoscopic scale (e.g. millimeter scale for concrete), heterogeneities play a significant role: for concrete, they are aggregates inclusions embedded within a mortar matrix. In order to explicitly take them into account within a Finite Element context, the authors turn to a method referred as non-adapted meshing process (Moës et al. (2003)). This consists in a unique "homogeneous" mesh with nodes placed independently from the heterogeneous morphology and thus from the interfaces. Hence, for a given micro-structure geometry, a set of elements are split into two parts by a physical interface. These elements are decomposed into two sub-domains that inherit of each material phase properties. In order to introduce this contrast of property in the classical linear elements, an enhancement of the kinematics is performed by means of a jump within the strain field -this being known as a weak discontinuity. The pioneer work on that matter goes back to the end of the 80's (Ortiz et al. (1987)) for the modeling of localization bands. The application of weak discontinuities in the context of material heterogeneities modeling has been introduced later in (Sukumar et al. (2001)).

Dealing with two-phase materials (e.g. matrix-inclusion), the non-adapted meshing method leads to two classes of elements: the ones that are completely within the matrix or within the inclusions, and the elements split into two parts. Only the latter type of elements are enhanced by the mean of weak discontinuities. Hence, the meshing process, referred to as morphological projection, consists in assigning their type to the original elements of the mesh. This methodology can be applied to three-dimensional random shapes, allowing complex geometry (based for example on excursion sets of random fields (Roubin (2013), Roubin et al. (2014)). This point is depicted on figures 1(b) and 1(c) (where matrix, inclusion and interface elements are represented in light, medium and dark grey, respectively). 2009)) (finite jump in the strain field) is introduced into each element in order to take into account this particular kinematics. This technique provides a meshing process that does not depend on the position and the shape of the inclusions. Cracking and softening behavior observed for quasi-brittle material such as concrete are introduced thanks to a strong discontinuity (Oliver et al. (2002)) (finite jump in the displacement field) in each elements. These discontinuities represent micro-cracks that can occur in any phases (aggregates or mortar matrix) and the debonding at the interface between aggregates and mortar matrix.

The remainder of the paper begins in Section 2 with a description of the kinematics for weak and strong discontinuities with which we work, and is then organized as follows: in Section 3, we introduce those kinematics into a Finite Element model, in Section 4, a general Finite Element discretization of the problem is presented and its resolution is described and applied to 4-node tetrahedron elements in Section 5. Finally, in Section 6, attention is focused on the capability of the model to represent the main features of cementitious materials by means of mechanical loading. First uniaxial tension and compression loading are performed. On the one hand, we show the emergence of the typical asymmetry of the respective macroscopic responses. On the other hand, we make a complete review of the crack patterns. The observations enlighten the close relationship between these two phenomena. In addition, non proportional loadings are performed in order to analyze the induced anisotropic behavior due to the failure process. This study is carried out on the previous uniaxial tension and compression tests by means of damage indicator -linear post-analysis on Young modulus degradation -as well as residual strength -non linear post-analysis on tensile strength degradation.

Kinematics description of weak and strong discontinuities

In this section, we summarize the basic notations employed all along the paper and remind kinematics of both weak and strong discontinuities.

We denote by Ω ⊂ R 3 a domain with smooth boundary ∂Ω and the spatial position of a particle by x ∈ Ω. This domain is discretized by means of standard isoparametric elements such as Ω = ∪ n elm e=1 Ω e . Naturally, position of weak discontinuities are defined by the material heterogeneities surfaces (interfaces) and therefore known prior to any mechanical calculation.

They are represented on figure 2 by the smooth contours of the grey shapes which have to be seen as the heterogeneities. In opposition, strong discontinuities (cracks) are introduced thanks to a stress-based localization criterion. Hence they can be, theoretically speaking, scattered anywhere in the domain (red curves on figure 2). However, in the presented framework, the physically meaningful assumption that they can be present only in the matrix (white zone) or following material heterogeneities surface is made. Whatever the origin of the discontinuity is, it will be called Γ d through this paper for sake of clarity.

Dealing with both weak and strong discontinuities, three cases can be present in an element Ω e : 1) only a weak discontinuity is present, 2) only a strong discontinuity is present, 3) both of them are present. The usual representation of (Jirásek (2000)) on figures 3(a) and 3(b) can be extended in the case of double enhancement as depicted on figure 3(c). In the first case, the shape of the material surface Γ d is defined by the shape of the heterogeneity, the direction of the unit normal n is taken orthogonal to Γ d and a jump of material properties passing through Γ d is present. In the second case, the direction of the unit normal n is defined by a standard principal stress criterion and a displacement jump passing through Γ d has to be considered. In the last case, the shape of the material surface Γ d is defined by the shape of the heterogeneity, the direction of the unit normal n is taken orthogonal to Γ d and a jump of material properties and of the displacement passing through Γ d are present.

It implies that the crack opening will be localized on Γ d or in other term on the physical boundary of each phase. As it is depicted on figure 2, such elements are split into two parts Ω e and Ω e so that Ω e = Ω e ∪ Ω e .

Ω

n Finally having at hand these two discontinuities, the total strain is written in the context of the enhanced strain fields (Simo and Rifai (1990)) such as:

Ω e Γ d Ω e Ω e Ω e Γ d Ω e n Ω e Γ u
ε = ∇ s ū regular + ε weak + ε strong , (1)
where ∇ s ū is the symmetric gradient of the displacement field. As in (Simo and Rifai (1990)), we refer to ε and ε as the enhanced parts of the strain field. The notation • (resp.

•) refers to weak (resp. strong) discontinuity.

Having those considerations in mind, we turn now to kinematics of weak and strong discontinuities. (2)

Whereas ũ itself remains continuous, the displacement gradient shall exhibit a jump across the surface discontinuity. Maxwell's compatibility conditions (Gurtin (1984)) lead to the form, that can be identified as the strain jumps. Actually, by considering the symmetric gradient of the displacement field ũ of components (ũ n , ũm , ũt ), ε can be written:

ũ = Θ n • (x -ξ) ([|ε|] n n + [|ε|] m m + [|ε|] t t) with Θ =    Θ ∀ x ∈ Ω e Θ ∀ x ∈ Ω e , (3
ε = ∇ s (ũ) =      ũn,n 1 2 (ũ n,m + ũm,n ) 1 2 (ũ n,t + ũt,n ) sym ũm,m 1 2 (ũ m,t + ũt,m ) sym sym ũt,t      , (4)
where (• ,n ), (• ,m ) and (• ,t ) represent derivatives in relation with n, m and t.

Then the difference between ε and ε can be expressed as follows

ε -ε =      [|ε|] n [|ε|] m [|ε|] t [|ε|] m 0 0 [|ε|] t 0 0      , (5)
where [|ε|] n , [|ε|] m and [|ε|] t stand for (ũ n,n -ũn,n ), (ũ n,m -ũn,m ) and (ũ n,t -ũn,t ), respectively.

All the specificity of the weak discontinuity kinematics lies on the enhancement displacement field ũ and its related strain field ε. The form of ε is now known: Θ has to be still exactly defined according to additional statical considerations. One can already note that the enhancement brings three parameters (stored in the so-called weak discontinuity vector

[|ε|]) that will be future unknowns of the mechanical problem.

Strong discontinuity kinematics

In this part, an element Ω e is only split by a discontinuity in the displacement field.

Thus only ε is present. Since strong discontinuity kinematics has been widely studied in the literature (see Simo et al. (1993), Oliver (1996) and Wells and Sluys (2001b) for instance) herein, focus is placed only on the key points of this method.

The displacement field u of an element Ω e crossed by a strong discontinuity can be mathematically decomposed, following (Simo and Oliver (1994)), by

u = ū + (H Γ d -ϕ e ) [|u|] , ( 6
)
where ū is a regular and continuous part of the displacement field that allows us to impose standard boundary conditions (Oliver (1996), part II), H Γ d is the Heaviside function centred on Γ d such as By taking the symmetric gradient of equation ( 6) and following the development proposed in (Simo and Oliver (1994)), we obtain the underlying strain field such as

H Γ d =    1 if x ∈ Ω e 0 if x ∈ Ω e , ( 7
∇ s u = ε = ∇ s ū + (H Γ d -ϕ e )∇ s [|u|] -([|u|] ⊗ ∇ϕ e ) s bounded + δ Γ d ([|u|] ⊗ n) s unbounded , (8)
where δ Γ d is the Dirac-delta distribution centred at the surface discontinuity and n its normal vector. Except from the fact that [|u|] represents the components of the displacement jump when evaluated on Γ d , this field will not be explicitly defined. At this stage, the resulting information of this reasoning is the decomposition of ε into a regular and bounded part εb and an unbounded one εu

ε = εb + εu = (H Γ d -ϕ e )∇ s [|u|] -([|u|] ⊗ ∇ϕ e ) s εb + δ Γ d ([|u|] ⊗ n) s εu . ( 9
)
And so the total strain field in presence of a strong discontinuity only has the form

ε = ∇ s ū + ε = ∇ s ū + εb + εu . ( 10
)
We can notice that strong discontinuity approach brings an unbounded part into the strain field. With such properties, physical considerations like continuity of the traction vector at Γ d or bounded stress in Ω e seem rather impossible to respect. Under the name of strong discontinuity analysis, the authors in (Simo et al. (1993)) propose solutions in order to ensure that constitutive models used are still consistent in presence of strong discontinuity. Among them, the so-called Discrete Strong Discontinuity Approach (DSDA) (Oliver (2000)) is used here. Hence, the continuum model equipped with strong discontinuity leads to an underlying discrete model at the discontinuity Γ d . On the one hand, a standard continuum relationship links strain and stress fields in Ω e \Γ d and on the other hand, a traction-separation law linking the traction vector to the crack-opening [|u|] on Γ d is introduced, modeling the failure mechanism.

Remarks

• ε (i.e. Θ) is not defined yet. It will be in Section 3.

• [|u|] is chosen as a constant function thus its symmetric gradient is equal to 0. This leads to an explicit formulation of the bounded part of the strong enhanced strain referred as Kinematically Enhanced Strain (KES). It is defined by

εb = -([|u|] ⊗ ∇ϕ e ) s . (11)
• We remind that the contribution of weak and strong discontinuities can be added if both of them occur, leading to the following strain enhancement

ε = ∇ s ū + ε + ε. (12)
• Additional set of unknowns [|ε|] and [|u|] has to be solved. The resolution will be detailed in Section 4. For the time being, it is sufficient to know that each set is defined at the element level. Hence, a local resolution will be performed for their determination letting the problem size unchanged.

Discontinuity implementation into a three-field variational framework

The Finite Element implementation of the problem is cast into a three-field variational framework. It is shown in this section how such a statement allows us to blend the presented strain enhancement into a Finite Element framework. Both yet undetermined shapes of the enhanced strains are defined following both static and kinematics considerations. Moreover, this formulation leads a handful incorporation of the discrete model using traction-separation law on Γ d .

Hu-Washizu three-field variational formulation

The main idea of the Hu-Washizu (Washizu (1982)) formulation is that, unlike classical displacement formulations, the three fields (ū, ε, σ), respectively the standard displacement field, the standard strain field and the standard stress field are considered independent. We call ( η, γ, τ ) respectively the virtual displacement, the virtual strain and the virtual stress fields. With this notation at hand, the Hu-Washizu formulation reads as follows:

• the classical equilibrium equation (or virtual work) (virtual displacement)

HW ū( ū, ε, σ; η) = Ω ∇ s η : σ dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, ( 13
)
• the kinematics equation (virtual stress)

HW σ (ū, ε, σ; τ ) = Ω τ : (∇ s ū -ε) dΩ = 0, ( 14
)
• the behavior equation (or constitutive model) (virtual strain)

HW ε (ū, ε, σ; γ) = Ω γ : (σ (ε) -σ) dΩ = 0, ( 15
)
where σ is a stress field that verified the behavior law and Ω η • ρb dΩ and Γt η • t d∂Ω the virtual work of the volumic and surfacic loadings. Note that in classical variational formulation (in displacement) we usually have, by definition

   ε = ∇ s ū σ = σ(ε)
which would make ( 14) and ( 15) irrelevant.

Having this three-field variational formulation at hand, we now turn in the next subsection to the so-called Assumed Strain Method developed in (Simo and Rifai (1990)).

Assumed Strain Method

The first idea of the Assumed Strain Method is to enhance both standard and virtual strain fields. Each of them can be decomposed into three parts. For the standard strain field, it gives

ε = ∇ s ū regular + ε weak + ε strong enhanced (16)
and for the virtual strain field

γ = ∇ s η regular + γ weak + γ strong enhanced . ( 17
)
Note that γ respects the form of ε (see equation 9) thus it is decomposed into a bounded part γb and an unbounded part γu .

Substituting equation ( 16) into equation ( 14) and combining equations ( 13), ( 15) and ( 17) we obtain the following modified three-field variational problem:

Ω ∇ s η : σ(∇ s ū + ε + ε) dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, ( 18
) Ω τ : ε dΩ = 0, ( 19
) Ω τ : ε dΩ = 0, ( 20
) Ω γ : (σ (∇ s ū + ε + ε) -σ) dΩ = 0, ( 21
) Ω γ : (σ (∇ s ū + ε + ε) -σ) dΩ = 0. ( 22
)
Following the second idea of the Assumed Strain Method, the explicit presence of the stress field is eliminated from the modified three-field variational formulation by choosing the space of the stress field L 2 -orthogonal to the space of the enhanced strains. As a result, equations ( 19) and ( 20) are satisfied and the second term of ( 21) and ( 22) vanish. And so the modified three-field variational formulation becomes:

Ω ∇ s η : σ(∇ s ū + ε + ε) dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, ( 23
) Ωe γ : (σ (∇ s ū + ε + ε)) dΩ = 0, ( 24
) Ωe γ : (σ (∇ s ū + ε + ε)) dΩ = 0. ( 25
)
It can be noted that due to the previous L 2 -orthogonal condition, the stress field σ does not appear in the final statement of the formulation. Moreover, since no inter element (enhanced strain and stress) continuity requirement is imposed, this condition can be applied independently on each n el -typical elements. Henceforth, equations ( 24) and ( 25) are evaluated at a the element level -on Ω e for e = (1...n el ).

The third idea of the Assumed Strain Method is to assure that after having enforced the L 2 -orthogonality condition, the stress field must at least include piece-wise constant functions; the combination of these two conditions (L 2 -orthogonality and piece-wise constant functions) implies satisfaction of the patch test (see Simo and Rifai (1990)). Mathematically

speaking, this gives

Ωe

γ dΩ = 0 and Ωe γ dΩ = 0. ( 26
)
Regarding the weak discontinuity, equation ( 26.1) can be written as

Ω e γ dΩ + Ω e γ dΩ = 0. ( 27
)
By assuming that the surface of discontinuity is flat (n is constant) within an element, solving (27) gives conditions on Θ. If we also assume that both Θ and Θ are constant in Ω e and Ω e , we obtain:

V Θ + V Θ = 0. ( 28
)
The choice retained here is: Θ = V /V and Θ = -V /V where V, V and V are the volumes of Ω e , Ω e and Ω e respectively.

Regarding the strong discontinuity, equation (26.2) can be written as

Ωe γb dΩ + Ωe γu dΩ = 0. ( 29
)
If we choose for γu a form respecting εu namely

δ Γ d ([|η|] ⊗ n) s , we get Ωe γb dΩ + Γ d ([|η|] ⊗ n) s d∂Ω = 0. ( 30
)
The same assumption of constant strain field and flat interface within an element are made.

Thus γb has the form of the so-called Assumed Enhanced Strain (EAS)

γb = - A V ([|η|] ⊗ n) s , (31)
with A the area of the discontinuity.

This method gives naturally an non-symmetric formulation: the virtual strain field space γb based on EAS (equation ( 31)) differs from the standard strain field space εb based on KES (equation ( 11)). It has been shown that this non-symmetric approach leads to superior numerical results in spite of the fact there is a variational inconsistency in the derivation (see Wells and Sluys (2001b)).

In the next section, the Finite Element discretization of equations ( 23), ( 24) and ( 25) is presented.

Finite Element discretization and resolution methodology

In this part, the Finite Element discretization of equations ( 23), ( 24) and ( 25) which correspond to the system to be solved is presented. First the discretization of the standard strain field and the virtual strain field respecting the form obtained in Sections 2 and 3 is shown, in a second time those discretizations are incorporated in equations ( 23), ( 24) and

(25) to get the Finite Element problem. Note that the next Finite Element discretizations are general and valuable for tetrahedron, beam or truss elements.

Note: Unless mentionned, Voigt notation is now used.

Strain fields discretization

Following the form obtained in Sections 2 and 3, the discretization of the standard strain field and the virtual strain field is written as follows:

ε = [ε xx ε yy ε zz 2ε xy 2ε yz 2ε xz ] T = B d ∇ s ū + G w [|ε|] ε + G s [|u|] εb , (32a)
γ = [γ xx γ yy γ zz 2γ xy 2γ yz 2γ xz ] T = B δ ∇ s η + G w [|γ|] γ + G * s [|η|] γ , (32b)
where several interpolation matrices introduced correspond to: B (= ∂N ) the standard strain interpolation matrix, G w the standard and virtual fields corresponding to the weak discontinuity, G s the standard field corresponding to the strong discontinuity (bounded part εb ) and G * s the virtual field corresponding to the strong discontinuity (bounded and unbounded part γ). d are nodal displacements and [|ε|] (resp. [|u|]) corresponds to additional unknowns arising from the kinematics enhancement of weak (resp. strong) discontinuity.

Regarding weak discontinuity kinematics, G w can be decomposed into Θ and a constant part H w that only carries information on the interface vector n such as:

G w =    G w = Θ H w = V V H w in Ω e G w = Θ H w = -V V H w in Ω e . ( 33
)
The KES interpolation matrix G s derives from equation ( 11) and can be computed by explicating the arbitrary function ϕ e . Such function can be defined using standard interpolation shape functions as follows:

ϕ e (x) = nen a=1 N a p a with p a =    1 if node number a ∈ Ω e 0 if node number a ∈ Ω e , (34)
where n en is the number of nodes in the element and p a the nodal values of ϕ e . G s is therefore the equivalent symmetric operator (• ⊗ ∇(ϕ e )) s in the Voigt notation of equation (32a).

The EAS interpolation matrix G * s is known from equation is decomposed in a bounded G * s,b and an unbounded G * s,u part that derives from equation ( 31) and ( 9), respectively. It leads to:

G * s = G * s,b + G * s,u = - A V + δ Γ d H * s , (35)
where H * s is the equivalent symmetric operator (• ⊗ n) s in Voigt notations of the strain field as in (32b).

Having those discretizations and the form of G w , G s and G * s at hand, equations ( 23), ( 24) and ( 25) can now be discretized.

Finite Element discretization of the problem

Injecting equations (32) into ( 23), ( 24), (25) and considering (33) and ( 35), the discretized problem to be solved reads as follows:

n el A e=1 [f e int -f e ext ] = 0, (36a)
h [|ε|] = V V Ω e H T w σ dΩ - V V Ω e H T w σ dΩ = 0 ∀e ∈ [1 .. n el ], (36b)
h [|u|] = - A V Ω e H * ,T s σ dΩ + Ω e H * ,T s σ dΩ + Γ d T d∂Ω = 0 ∀e ∈ [1 .. n el ], (36c)
where all the σ linearly depend on d, [|ε|] and [|u|] and where

f e int = Ω e B T σ dΩ + Ω e B T σ dΩ, (37a)
f e ext = Ωe N T ρb dΩ - Γt N T t d∂Ω and (37b)
T = H * ,T s σ == σ • n Matrix is the traction vector. ( 37c
)
The behavior law written in relation with the spatial position x is As the whole framework fits into the Discrete Strong Discontinuity Approach, two behaviors have to be considered. First, regarding the continuum bodies Ω e and Ω e , the behavior is considered elastic. This rather strong assumption is justified by the general spirit of simple meso-scale modeling. However, more complex behaviors such as plasticity or damage can also be implemented (see Oliver (1996) for details). Then, the second part of the modeling takes place at the discontinuity surface Γ d . A governing law that links the traction vector

σ(d, [|ε|] , [|u|]) =    σ = C ε = C Bd + G w [|ε|] + G s [|u|] if x ∈ Ω e σ = C ε = C Bd + G w [|ε|] + G s [|u|] if x ∈ Ω e , (38)
T = σ| Γ d
• n 1 is defined on the surface Γ d where the crack opening magnitude drives the non-linear failure mechanism. Such laws are often referred as traction-separation laws. In 1 Matrix notation order to fit in the general spirit of simple meso-scale modeling, a very basic single traction criterion and brittle softening is used here. The next subsection describes its main characteristics and shows how it is incorporated in the Finite Element problem (36a, 36b, 36c).

Finally, after linearising these equations, the Finite Element problem is written under a matricial form and the solving strategy is presented.

Discrete constitutive model on Γ d

The discrete model is based upon a relationship between the traction vector T (via an equivalent stress σ eq ) and the crack opening magnitude labelled [u]. The equivalent stress σ eq value has to be tackled with careful attention in order to depict a meaningful representation.

Thus two cases have to be considered. First, if we consider the case where only the strong discontinuity is present (case 1) of Section 1), no geometrical information are given on the crack orientation. The choice is made here to use the larger principal stress component as the equivalent stress σ eq = σ I . This principal component is simply the first eigenvalue of the stress matrix. Furthermore, the corresponding eigenvector n I represents its direction.

The physical meaning of eigenvalue problems naturally leads to the choice of this vector in order to represent the crack orientation n ← n localization

I

. Secondly, if both weak and strong discontinuities are present (case 3) of Section1), it is assumed that crack opens at the interface between each phases. The main difference with the latter case is that n is defined by geometrical characteristics and therefore is independent from any stress state.

Hence, traction vector can be defined prior to localization. In this case, its projection on the interface direction is used such as σ eq = n • T . Attention is drawn to the fact that in this case, the shear components are non-zero but we assume they do not participate to the localization process. This is a major assumption that could be improved by considering more complex criteria. However, in the spirit of multiscale analysis we aimed at first considering the most simple case.

The discrete model is introduced when the equivalent stress σ eq reaches a certain amount σ y , which has to be seen as a meso-scale material characteristic. Mathematically speaking, the strong discontinuity is introduced when the following localization criterion becomes zero:

Φ l = σ eq -σ y . (39)
Then, the failure mechanism is driven by a scalar opening criterion Φ o defined as follows:

Φ o = σ eq -(σ y -q), (40)
where

q = σ y 1 -exp - σ y G f [u] . (41)
It can be seen that a second material parameter G f called the fracture energy that governs the amount of energy necessary for the complete material failure is introduced in (41). Since [u] represents the magnitude of the crack opening [|u|], the criterion is single valued. However, an additional projection of the crack has to be defined in order to fit in the previous framework, introducing a unit jump vector n p :

[|u|] = [u] n p . (42)
Finally, the behavior can be split into two parts: an elastic one in the body Ω e (linear relationship between strain and stress fields out of the discontinuity Γ d ) as drawn on figure 4(a) for a 1D case and a discrete one on the discontinuity Γ d by means of a traction-separation law that links Ω e and Ω e as shown on figure 4(b). Since several choices retained here depend on element kinematics, details on the different ingredients of this model (σ eq , n p . . . ) will be given later.

From equation (36c) and making the assumption of flat interface and constant stress, an easy integration gives an explicit expression of T as a average value of σ and σ weighted by volumes such as:

T = 1 V H * ,T s (V σ + V σ ), (43)
where V, V and V are the volumes of Ω e , Ω e and Ω e , respectively.

T = 1 V H * ,T s V C + V C B K s * b d (44)
+ V V V 2 H * ,T s C -C H w K s * w [|ε|] ( 4 5 ) + 1 V H * ,T s V C + V C G s n p K s * s [u]. (46)
Following (Ibrahimbegovic et al. (1998)), the criterion Φ o is now incorporated in the system (36a, 36b, 36c). By injecting equations ( 44), ( 45) and ( 46

∆Φ o = ∂σ eq ∂T ∆T + ∂q ∂[u] ∆[u] = ∂σ eq ∂T 1 V H * ,T s V C + V C B K s * b ∆d + ∂σ eq ∂T V V V 2 H * ,T s C -C H w K s * w ∆ [|ε|] + ∂σ eq ∂T 1 V H * ,T s V C + V C G s n p K s * s ∆[u] + σ 2 y G f e -σy[u]/G f Kq ∆[u].
Thus the linearisation of Φ o = 0 gives

K s * b (k) n+1 ∆d (k+1) n+1 + K s * w (k) n+1 ∆ [|ε|] (k+1) n+1 + (K s * s + K q ) (k) n+1 ∆[u] (k+1) n+1 = -Φ o (k) n+1 . (47)
In the next subsection, we present the form of the final problem to be solved under a matricial form and the resolution strategy.

Global system and resolution strategy

Because the linearisation of equations (36a) and (36b) is trivial the mathematical development will be skipped. The global system to be solved in terms of increments of d, [|ε|] and [u] and condensed in a matricial format is

     K bb K bw K bs K wb K ww K ws K s * b K s * w K s * s + K q      (k) n+1      ∆d ∆ [|ε|] ∆[u]      (k+1) n+1 =      - n el A e=1 {f e int -f e ext } -h [|ε|] -Φ o      (k) n+1 , (48)
where

K bb = B T V C + V C B, K bw = V V V B T C -C H w , K bs = B T V C + V C G s n p , K wb = V V V H T w C -C B, K ww = V V V 2 H T w V C + V C H w , K ws = V V V H T w C -C G s n p , K s * b = ∂σeq ∂T 1 V H * T s V C + V C B, K s * w = ∂σeq ∂T V V V 2 H * T s C -C H w , K s * s = ∂σeq ∂T 1 V H * T s V C + V C G s n p , K q = σ 2 y G f e -σy [u] / G f . (49)
Solving system ( 48) is done at two levels. First, following the operator split method spirit,

variables (∆ [|ε|] , ∆[u]
) are determined at the element level (local solving) for a given ∆d by solving: . ∆d is calculated by solving the global equilibrium equation using a static condensation (Wilson (1974)) on the local (known) variables (∆ [|ε|] , ∆[u]). This leads to a matricial system, with a modified stiffness matrix K sc , to be solved such as

   h [|ε|] = 0 Φ o = 0 . ( 50
K sc (k) n+1 ∆d (k+1) n+1 = - n el A e=1 {f e int -f e ext } (k) n+1 , (51)
where

K sc (k) n+1 = K bb -K bw K bs      K ww K ws K s * w K s * s + Kq   (k) n+1    -1   K wb K s * b   (k) n+1 .
Even though the stiffness matrix has been changed due to kinematics enhancement, both its size and sparsity are unchanged. Hence, no matter how many heterogeneities are represented or how many elements have starting to fail, the global size of the problem is preserved. In terms of numerical resources, the memory needed only depends on the mesh size (number of nodes). Naturally, local Newton algorithms slows down the global calculation as the number of strong discontinuity activated increases. By using the static condensation, a standard FE problem is retrieved, where increments of d have to be found in order to respect the global equilibrium equation 51.

The most common method used to solve those problems in case of non linearity are the so-called Newton methods. However, it requires full calculation of full stiffness matrix at each iteration, and since morphological modeling requires rather fine meshes, a quasi-Newton algorithm coupled with an iterative solver is used here. Among the huge diversity of those algorithms, the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm is retained. It makes the number of arithmetical operations to fall down from O(n 3 ) to O(n 2 ). The price to pay is that quasi-linear convergence is obtained (instead of quadratic). Full details are given in the original papers (Broyden (1970a,b); Fletcher (1970); Goldfarb (1970); Shanno (1970)).

Moreover numerical implementation details are in (Matthies and Strang (1979)). Added to the quasi-Newton BFGS, a line-search method is also used, modulating the incremental displacement norm for each iterations by a factor s: k+1) .

û(k+1) = û(k) + s (k+1) ∆d (
See Dahlquist (2003) for details on the computation of s. It is recalled that powerful algorithm adapted to the E-FEM such as (Oliver et al. (2008)) and not implemented here can significantly increase the computation performance.

In the next section, we move to applications of this model in the case of 4-node tetrahedron elements. First the explicit form of the matrices encountered all along this paper is presented for sake of clarity. Then the criterion Φ o is written in the considered case.

Application to 4-node tetrahedron elements

The use of 4-node tetrahedron elements is now presented in the case of this model. This volumic discretization brings several advantages compared to the initial implementation of this framework based on lattice discretization (Benkemoun et al. (2010)). First, an exact representation is made regarding the volumic tessellation of the mesh, leading to exact representation of constant stress problem. Therefore, the problematic of mesh convergence presented in the aforementioned paper with lattices is irrelevant in this context. Furthermore, a complete kinematics of cracks opening can now be represented (mode I, II). Hence, both geometrical construction and mechanical behavior are more accurately depicted.

Interpolation matrices

In this part, attention is drawn to the explicit definition of the matrices encountered through this paper in the case of tetrahedron elements. In addition, for sake of convenience, fields are represented in their Voigt notation. Hence, each matrices are developed in this format.

In the case of an element cut by a discontinuity Γ d , the tetrahedron is split into two sub-domains Ω e and Ω e delimited by the interface Γ d of direction vector n (see figure ( 5)).

It is reminded that this surface is assumed to be flat (n is constant over Γ d ). Numerical implementation of such geometrical construction is not trivial and several cases have be taken into account depending on the surface orientation. Indeed, Ω e and Ω e can be polyhedron of respectively 6 and 4, 5 and 5 or 4 and 6 nodes. And since the formulation involves only volumes V and V (not the discontinuity area, see equation ( 43)), it has to be considered with utmost attention.

n Dealing with weak discontinuity, it is reminded that the interpolation matrix G w can be decomposed into Θ, a piece-wise constant function depending on the considered subdomain and H w a matrix containing information on the discontinuity surface orientation (see equation ( 33)). Following the same Voigt convention, H w can be constructed so that the vector format of the strain enhancement matches its tensor definition in equation (33).

Ω e Γ d m t Ω e
The vector n is written in the global coordinate system n = [n x n y n z ] T , leading to the following interpolation matrix format:

H w =               n 2 x n x m x n x t x n 2 y n y m y n y t y n 2 z n z m z n z t z 2n x n y n x m y + n y m x n x t y + n y t x 2n y n z n y m z + n z m y n y t z + n z t y 2n x n z n x m z + n z m x n x t z + n z t x               . ( 52
) G w =    G w = Θ H w = V V H w in Ω e G w = Θ H w = -V V H w in Ω e . ( 53
)
Regarding strong discontinuity, it is reminded that two matrices have to be considered, namely G s and G * s,b , making the formulation non-symmetric. The form of G s is given by

G s = -              
0 ∂ϕe ∂x               (54)
and the form of G * s,b by

G * s,b = - A V H * s with H * s =               n x 0 0 0 n y 0 0 0 n z n y n x 0 0 n z n y n z 0 n x               . ( 55
)
All components of the enhanced strain field are now explicitly defined. The standard strain field can thus be built respecting the form of equation (32a) and the system (48) solved.

We now turn explanation on the localization and opening criterion.

Localization and opening criterion

Since the discontinuity surface orientation is constructed with geometrical properties for interface element and with stress consideration otherwise, two cases have to be treated separately.

Finally, by defining the projection vector n p of equation ( 42 Each numerical example given in this chapter is made on a heterogeneous 100 × 100 × 100 mm 3 specimen for which two phases are modeled. Based on unions of excursion sets (Adler (2008), Roubin et al. (2014)), this method explicitly represents aggregates of different sizes melt within a matrix that is, roughly speaking a mortar. The former are modeled with three average diameters of 15, 7 and 5mm representing respectively 25, 50 and 25 % of the total 30 % volume fraction (see figure 6(b)). Random shaped inclusions are yielded by thresholding correlated Random Fields (see Roubin (2013)) and performing the union.

Following the non-adapted mesh spirit, once projected onto the FE discretized space, this morphology is represented by two kinds of elements: those which are completely included within the matrix or an aggregate and those, close to an interface, that are split and thus enhanced by a weak discontinuity. Table 1 summarizes the material properties at the mesoscopic scale. It is reminded that the model contains, for each phase, two elastic and two failure parameters, the Young modulus E, the Poisson ratio ν, the yield stress σ y and the fracture energy G f , respectively.

Furthermore, a yield stress and a fracture energy are also set to define the interface.

It is worth noting that:

• Aggregates are assumed to remain elastic. Thus, no crack may initiate or propagate within those elements.

• Weakly enhanced elements do not require any specific characteristics in the elastic regime. Indeed, as presented in the previous chapter, prior to any debonding (modeled through a strong discontinuity), they model a perfect interface, thus with infinite rigidity.

• For sake of simplicity, the failure mechanisms are taken to be equal for both matrix and interfaces. However, the difference in elastic properties of the matrix and aggregates leads to stress concentrations at the interface. Since the failure criterion is based on stress considerations, a similar yield stress in both matrix and interface still leads to a weaker behavior of the latter.

Finally, the computations are performed under uniform displacement control along the first spatial axis (X-direction). Moreover two other faces of the domain (normal to the Y and the Z direction) have constant zero value for their normal displacement. Hereafter, a value corresponding to the X-direction is referred as axial whereas any transversal quantities refers to an average value set up on both the Y and Z-directions. First it is worth noting that, although the failure criterion at the meso-scale is triggered in tension only, macroscopic simple compression leads to the failure of the specimen. Indeed this feature is clearly a consequence of the structural effect that is set up by the explicit representation of the aggregates.

Phase

E [GPa] ν [-] σ y [MPa] G f [J.m -
Second, it clearly appears that the macroscopic failures are unsymmetric. This feature is typical of quasi brittle materials such as concrete. Table 2 sums up macroscopic key values extracted from figure 7: we denote by E M the macroscopic modulus, σ M f the macroscopic tensile or compressive strengths, ε M f the corresponding macroscopic failure strain and finally D p the total dissipated energy. The last value is computed through integration of the macroscopic axial force over the axial displacement.. As expected the initial elastic behavior is strictly symmetric. Even though it is quite difficult to observe the end of the purely elastic region on the plots, it can safely be assumed that, for compression, the transition to the non linear behavior occurs for a more important macroscopic stress than for tension. This unsymmetric elastic domain is the first feature that emerges from the multi-scale analysis.

Loading path E M [GPa] ν M [-] σ M f [MPa] ε M f [-] D p [J]
Regarding the energy D p needed to reach the specimen ruin, it can be noted that both are greater than the mesoscopic fracture energy G f (in J.m -2 ) assigned to the mortar (corresponding values of D p in tension and compression are 50 and 5 000 J.m -2 , respectively).

Furthermore, the fact that this energy is significantly greater in compression reflects the more brittle behavior of concrete when tested in tension. It is naturally linked with the asymmetric strength values and their corresponding failure strains for which the ratio of compression to tension are both 10. A discussion on this meaningful result is given in the main conclusion of this paper. . This corresponds to the elements for which a strong discontinuity has been introduced. It can be seen that those crack patterns are very significant either in tension or compression. First in tension, some micro-cracks are linked in order to set up a single macro-crack that roughly lies in the transversale plane (best seen on Figure 9) . Obviously this macro-crack is tortuous and goes around the aggregates that remain elastic. Second, in compression, it can be observed that several macro-cracks are present (contrary to tension) and that they are roughly parallel to the axial direction.

Aggregates properties E (-40%) ν (+50%) Table 3 shows the influence of the elastic properties of the aggregates on both compression and tension macroscopic strengths. As expected, the ratio between those two quantities is largely influenced by the contrast of Young's modulus. To be more precise, the less this contrast is, the more the compressive strength increases and so the ratio (the tension strength remains almost constant). On the contrary, the Poisson ratio of the aggregates seems to have a small influence on the macroscopic strengths: from ν = 0.2 to ν = 0.3 the compressive strength shows a 3% increase and the tensile a 5% decrease.

On a more general point of view, the question of the geometrical representation of a macroscopic crack using a local method as the Embedded Finite Element is non trivial. For example, in the two-dimensional case, using constant triangular elements, it has been shown (see Jirásek (2012) for details) that the best way to produce a suitable crack trajectory and avoid numerical issues such as stress locking is to combine two methods. First a non-local formulation of the smeared crack approach, giving crack orientation in each element. Then a tracking algorithm to enforce the crack path continuity between each element. The major drawback of this implementation is that the local spirit of the E-FEM (directly inherited from the FEM itself) is lost. Indeed, in addition to non local damage, path continuity enforcement implies, for an element, a crack position that depends of those of its neighbors. Moreover, in the three-dimensional case, continuity of flat plane (crack) is often impossible.

Herein, the E-FEM implementation has to be placed within the multi-scale context. In this case, a single fractured element is not considered to be representative of any specific macroscopic feature. However, it is only when a large number of those activated elements are merging that it may be considered that they model a continuous path at the macroscopic scale. It is for these reasons that, herein, no specific effort has been made in order to enforce any path continuity at the meso-scale. It is reminded that herein, the discontinuity is considered piecewise constant in each element. In the case of weakly discontinuities (interfaces elements), the orientation and position of the strong discontinuity are predefined by the crossing heterogeneity. However, in the case of standard kinematics (matrix element), the orientation of the strong discontinuity is set to be the direction of the larger principal stress at the localization time. In this case it is assumed that the discontinuity path through the centroid of each tetrahedron. This assumption impacts only the construction of the functions ϕ e .

Finally, features like multi-cracking or branching, which usually require a complex local numerical implementation (within an element), are herein omitted at the meso-scale. However, as shown on Figure 11 (which is a zoom made on a subset of Figure 10(b)), it can be retrieved at the macro-scale. This picture shows a crack that splits in two branches; a main branch (on the top) with larger opening values and a second branch that eventually vanishes. Generally, these branchings come from an aggregate "blocking" the way of the crack propagation direction.

Transversal strain analysis

Regarding the transversal behavior of the specimen, several observations can be made either in tension or in compression. It is recalled that the axial direction X corresponds to the imposed displacement direction and transversal values are defined as the average of the values along Y and Z. The results presented here are still based on the same one-dimensional macroscopic tests (tension/compression) mentioned above. The macroscopic Poisson ratio can be determined using the transversal strains by:

ν M = - ε M tr ε M ax (62)
Herein, this property is extended to the diffuse cracking regime in order to illustrate the mechanism that leads to the specimen loss of rigidity -see figure 12 where it is plotted in terms of axial strain for both tension and compression tests. First, the elastic part shows that the macroscopic Poisson ratio is of the same order as for the meso-scale: ν M = 0.195 in both cases. Afterwards those values are diverging. First, in tension, the ignition of diffuse cracking causes local strain to release and thus making the macroscopic strain decrease with the Poisson ratio. On the contrary, in compression, this local strain releasing causes a heightening of the transversal mechanism leading to a significant increase of the apparent Poisson ratio.

Since the meaning of the Poisson ratio is highly contestable with strongly non-linear failure behavior, the post-localization analysis is only based on transversal strains. For that matter, Figure 13(b) and Figure 13(a) show the macroscopic response up to the specimen ruin in terms of the axial strain ε M ax (solid curve) and of the transversal direction ε M tr (dashed curve). Regarding the tension test, during the post-peak phase the transversal strain decreases and tends to vanish (see Figure 13(b)). It represents the unloading that occurs in the specimen -at a macroscopic scale -after the main crack localization. In contrast, during a compression test, the transversal strain still increases after the peak load (see Figure 13 The dilatancy δ of the specimen can also be computed by considering the trace of the macroscopic strain tensor. Thus, the relative variation of the volume, drawn on Figure 14, is:

δ = ∆V V 0 = ε M ax + ε M tr Y + ε M tr Z (63)
Notice that the same magnitude in the dilatancy rate is observed in tension and in compression regarding the elastic region. This result is in conformity with the identical Poisson ratio value. Naturally a tension test produces a volumetric expansion (δ > 0)

while compression first produces a contraction (δ < 0). However, the cracking process increases the dilatation rate for both tests, which is a major feature of many materials such as concrete. Hence, when in tension the slope increases in the diffuse cracking regime, for compression the contraction slows down. The localization is characterized by: in tension a sudden increase in rate and in compression a maximum (considering absolute value) of the dilatancy. Afterwards, the dilatancy sign changes in the post-localization region, which is an interesting feature, also experimentally observed for concrete. In (Torrenti (1987)) the interpretation is made that, in compression, the contraction corresponds to a predominant elastic effect where the following expansion reveals an important cracking stage. However, experiences revealed a volumetric strain which switches sign prior to the localization.

Dissipated energy

By integrating the macroscopic force over the displacement, the total energy can be calculated from the macroscopic response. It is shown on Figure 7 for each calculation step.

Furthermore, considering a fictitious elastic unloading up to a zero force level, a difference can be made between the elastic and the dissipated part. Thus it gives an additional way to understand the failure mechanisms and to compare tension and compression behaviors. in compression can be seen. The fall of the elastic energy to a nearly null value represents the unloading outside the macro cracks region and the fact that all the energy is dissipated close to the macroscopic cracks.

Numerical evaluation of the induced anisotropy

After having shown the macroscopic responses of the model in relation with simple radial and proportional macroscopic loading paths, focus is placed on non proportional cases. Here, our objective is to show the emergence, at macro-scale, of some features related to any anisotropic behaviors. This anisotropy shall be revealed by the non proportionality of the loading path.

Considering a tension or a compression test with monotonic loading, the failure mechanism that leads to the specimen ruin induces a strong anisotropy of cracks pattern. On the basis of this simple observation, an analysis of the macroscopic material propertiese.g. Young modulus, tensile strength, etc. -for each step of the previous monotonic loading tests (in both tension and compression) is now proposed. Hereafter, the two first parts of the macroscopic loading path (tension or compression) are referred as principal calculations. In order to yield residual material properties, additional calculations, that inherit from the principal, are performed. They are referred as secondary computations. Basically, the inheritance from principal to secondary calculation is made through the non linear data, i.e. the whole set of meso-scale cracks (with their orientations and opening values). This second part of the numerical analysis is also made under displacement control. Actually the displacements obtained from the first part are imposed and additional displacements are added. Those displacements correspond to the same kind of boundary conditions along a different direction (Y or Z).

Anisotropic induced damage

Here the residual property of concern is the damage Young modulus, which is computed as the secant modulus on the macroscopic stress -strain curve. In order to catch the anisotropy, from each step of the principal calculation three secondary calculations are performed in the three directions. It leads to three macroscopic secant moduli: an axial reaches 0.15. Hence, the elastic property is far more spoiled in the axial direction. This result reflects the characteristic morphology of the cracks pattern, splitting the specimen in two by a plane roughly perpendicular to the axial axis. As the macroscopic crack grows, the "link" between each part of the specimen becomes weaker, leading to a decrease of the upscaled secant modulus in this direction. On the contrary, it remains several non-broken paths on the transversal directions that give to the specimen a higher rigidity. Regarding the compression test, Figure 16(a) shows the opposite effect. Indeed, here the transversal damage is more important than the axial one. It may be explained by cracks patterns that form planes parallel to the axial direction, leading to a higher loss of rigidity in the transversal directions. Furthermore, the more diffuse aspect of the crack repartition makes the difference between axial and transversal damage less important.

Finally, for both tension and compression cases, it can be noticed that the two transversal damages are of the same order of magnitude, representing isotropic behavior in these two directions. Somehow, it can be said that the macroscopic elastic behavior is shifting from an isotropic case to a transverse isotropic one. More numerical investigations for the second step of those non proportional loading paths may determine the complete elasticity tensor for this case. The results show approximately the same behavior as those for the elastic moduli. Regarding the tension failure, a more important decrease of the tensile strength is observed in the axial direction than along the transversal directions. The ratios are also of the same order of magnitude. Indeed, when the specimen has lost ≈ 80 % of its strength in the former direction, it has only lost ≈ 20 % in the last two. Regarding the compression test, the specimen seems to follow a rather isotropic behavior. Nevertheless, the transversal resid- However, regarding the transversal directions, it can be noticed that this state seems rather more deteriorated than what experimental results shows.

Discussions and conclusions

This paper first sets a numerical tool for the modeling of the quasi-brittle behavior of heterogeneous materials such as the cementitious ones. This model can be viewed as a multiscale model, aiming at upscaling data from mesoscopic scale (millimeter scale for concrete) to macro-scale.

At meso-scale the numerical implementation is based on the introduction of kinematics enhancements of two kinds within the FE context. The first enhancement, referred to as weak, leads to a non-adapted mesh strategy for heterogeneous morphologies. It can therefore be "simply" projected onto an unstructured mesh, freeing us of any expensive algorithm that aim to match a mesh onto a given morphology. Indeed, according to the nonadapted mesh point of view, this mesh is created regardless of any physical surface (i.e. the interface between each component). The second enhancement, referred to as strong, models discontinuities in the displacement field that are viewed as micro-cracks. Their opening incompatibilities results in shear loading. Furthermore, the ratio of weakly enhanced element increase along with the fraction volume which highlight even more the need of mode II opening. Finally, if not sufficiently discretized by the FE mesh, the morphology can exhibits "elastic paths" within the aggregates, thus explaining partially the stressed state at the end of the loading. Due to the inner complexity of this topological quantification, only a visual validation has been made, thus leading to reasonable doubts as to the disconnected state of all discretized inclusions. Moreover the model clearly offers a large number of perspectives dealing with multiphysics. This is a major point associated to the durability of concrete structures, for which mass transfers are the cornerstone. Thus improving the physics of the modeling can be made:

(1) through a coupling between the cracks pattern (obtained from a mechanical analysis) to simple flows (such as the Poiseuille flow between two planes). ( 2) according to a better morphological representation of the heterogeneities, e.g. accounting for the largest porosity that can be found within mortar matrix. On a more general point of view, the authors think that the use of the E-FEM in the context of heterogeneous materials shall be applied, in a near future, to smaller scales (typically the micrometer scale dealing with concrete), allowing for a better understanding of the role played by simple failure mechanisms and their propagation to macro-scale. Thus the future development of this model lies in retaining the spirit of the physical significance of fine scales.

Figure 1 :

 1 Figure 1: Projection of two-phase morphologies on a 3D mesh

Figure 2 :

 2 Figure 2: Illustration of the two kinds of discontinuous elements used in Ω e . On top right; a so-called matrix element embedded with only strong discontinuity and on bottom right; a interface element with both weak and strong discontinuity

Figure 3 :

 3 Figure 3: From (Jirásek (2000)), element with: a) one weak discontinuity, b) one strong discontinuity and c) both weak and strong discontinuity

  Equation (36a) is the global equilibrium equation of a standard Finite Element problem whereas equations (36b) and (36c) are local equations added by the presence of the enhanced parts of the strain field. It is important to recall that they are solved at the element level, allowing us to dertermine the values of [|ε|] and [|u|] by performing a local resolution.

  Figure 4: Elastic/brittle behavior

Figure 5 :

 5 Figure 5: 4-node tetrahedron element with discontinuity surface

Figure 6

 6 Figure6(a) shows the FE mesh used. Basically, GMSH(Geuzaine and Remacle (2001)) is used to produced the mesh, which is based on the Delaunay triangulation of a set of randomly positioned nodes. In order to catch the geometrical information of the smallest heterogeneities, the mesh used to compute the following examples has about 556 10 3 nodes (1 600 10 3 dof) and 3 500 10 3 elements. Figure 6(b) shows the projection of a typical morphology onto the mesh (only weakly enhanced and aggregates elements are represented in light and dark grey respectively).

Figure 6 :

 6 Figure 6: Projection of a typical excursion set morphology onto the chosen mesh

Figure 7

 7 Figure7shows the macroscopic response obtained through the computation for both a simple tension and a simple compression loading. It plots the macroscopic axial force versus 29

Figure 7 :

 7 Figure 7: Macroscopic response for simple tension and compression tests

  Figure 8: Distribution of the crack openings at the last time step of the simple tension test

Figure 9 :

 9 Figure 9: Macroscopic crack paths for simple tension: 2D slices

Figure 10

 10 Figure 10 plots the crack patterns obtained at the end of the computation for both compression (Figure 10(a)) and tension (Figure 9 and Figure 10(b)). This corresponds to

Figure 10 :

 10 Figure 10: Macroscopic crack paths for simple tension and compression

Figure 11 :

 11 Figure 11: Branching of macroscopic cracks around the aggregates

Figure 12 :Figure 13 :

 1213 Figure 12: Apparent Poisson ratio for tension and compression tests vs. axial strain

Figure 15 :

 15 Figure 14: Macroscopic dilatancy evolutions for tension and compression

Figure 16 :

 16 Figure 16: Axial and transversal damage variables evolutions

Figure 17 :

 17 Figure 17: Evolution of the residual tensile strength during the failure process

Figure 1 :

 1 Figure 1: Projection of two-phase morphologies on a 3D mesh

Figure 2 :

 2 Figure 2: Illustration of the two kinds of discontinuous elements used in Ω e . On top right; a so-called matrix element embedded with only strong discontinuity and on bottom right; a interface element with both weak and strong discontinuity

Figure 3 :

 3 Figure 3: From (Jirásek (2000)), element with: a) one weak discontinuity, b) one strong discontinuity and c) both weak and strong discontinuity

)

  where ξ represents the position of Γ d , thus letting n • (xξ) act has the signed distance to the surface discontinuity, and Θ a still unidentified piecewise constant shape function with unit jump at discontinuity surface (a specific definition is given in section 3) . The introduction of the Θ function leads to a different definition of ũ weather it is evaluated in Ω e or in Ω e . In (3), [|ε|] n , [|ε|] m and [|ε|] t are three constants

  ) ϕ e is an arbitrary continuous shape function in Ω e of unit value at each nodes in Ω e and null at each nodes in Ω e and finally, [|u|] is a continuous function on Ω e . The components of the displacement jump at the discontinuity are given by [|u|] x∈Γ d .

  Equation (36a) is the global equilibrium equation of a standard Finite Element problem whereas equations (36b) and (36c) are local equations added by the presence of the enhanced parts of the strain field. It is important to recall that they are solved at the element level, allowing us to dertermine the values of [|ε|] and [|u|] by performing a local resolution.

  Figure 4: Elastic/brittle behavior

  ) in Φ o , we obtain a criterion depending on (d, [|ε|] , [u]). Since such a criterion leads to a non-linear equation, one has to linearise it. Rewriting it under an incremental form and collecting terms together, Φ o leads to

)

  By developing h [|ε|] = 0 as done previously in (36b), one can note that this equation is linear, thus the non-linear aspect of the local system (50) comes only from the equation Φ o = 0. This one imposes a standard Newton-Raphson procedure implemented within element subroutine. Once system (50) is solved, appropriate values of ∆[|ε|]

Figure 5 :

 5 Figure 5: 4-node tetrahedron element with discontinuity surface

  ) as the normal vector n, a mode I opening mechanism is represented. Attention is draw to the fact that, in this case, the evolution of the shear tractions are not driven by the displacement jump. That is a major assumption which represents a flaw in the model. Improvements may come from more complex opening mechanisms such as mode II or mode I+II. Those are both under development and related issues are discussed in the conclusion of the paper.At this stage, the whole system (48) is now explicitly known in the case of tetrahedron elements. Thus we can move to numerical examples in the context of concrete-like material to show the features of the model developed through this paper.6. Numerical analysis of concrete: from meso to macro-scale In this section, attention is focused on the capability of the enhanced FE model to represent the main features of concrete by means of mechanical loading. Starting from mesoscale (millimeter scale), uniaxial tension and compression tests are performed in order to show, on the one hand, the emergence of the typical asymmetry of the respective macroscopic responses and on the other hand, a complete analysis of the crack patterns. Moreover, non proportional loadings are carried out in order to measure the induced anisotropy. This study is based on the previous uniaxial tension and compression tests by means of damage indicator (i.e. linear post-analysis on Young modulus degradation) as well as residual strength (i.e. non linear post-analysis on tensile strength decrease).

Figure 6

 6 Figure6(a) shows the FE mesh used. Basically, GMSH(Geuzaine and Remacle (2001)) is used to produced the mesh, which is based on the Delaunay triangulation of a set of randomly positioned nodes. In order to catch the geometrical information of the smallest heterogeneities, the mesh used to compute the following examples has about 556 10 3 nodes (1 600 10 3 dof) and 3 500 10 3 elements. Figure 6(b) shows the projection of a typical morphology onto the mesh (only weakly enhanced and aggregates elements are represented in light and dark grey respectively).

Figure 6 :

 6 Figure 6: Projection of a typical excursion set morphology onto the chosen mesh

Figure 7

 7 Figure7shows the macroscopic response obtained through the computation for both a simple tension and a simple compression loading. It plots the macroscopic axial force versus 29

Figure 7 :

 7 Figure 7: Macroscopic response for simple tension and compression tests

Figure 8 :

 8 Figure 8: Distribution of the crack openings at the last time step of the simple tension test

Figure 9 :

 9 Figure 9: Macroscopic crack paths for simple tension: 2D slices

Figure 10

 10 Figure 10 plots the crack patterns obtained at the end of the computation for both compression (Figure 10(a)) and tension (Figure 9 and Figure 10(b)). This corresponds to

Figure 10 :

 10 Figure 10: Macroscopic crack paths for simple tension and compression

Figure 11 :

 11 Figure 11: Branching of macroscopic cracks around the aggregates

Figure 12 :Figure 13 :

 1213 Figure 12: Apparent Poisson ratio for tension and compression tests vs. axial strain

  Three energies -elastic, dissipated and total -are plotted on figure15as a function of the axial imposed displacement for both loadings.These curves give a clear representation of the reversible and irreversible mechanisms that occur during the loadings, which are typical of a softening behavior. It can be observed that, at the beginning, nearly all the energy involved is elastic. Then, at the localization stage, the part of energy dissipated increases significantly. A more brittle failure in tension than -

Figure 14 :Figure 15 :

 1415 Figure 14: Macroscopic dilatancy evolutions for tension and compression

ẼM

  ax value and two transversal ones ẼM tr Y and ẼM tr Z . The results are displayed by defining "damage" variables d ax , d tr Y and d tr Z , respectively. They are built to compare the upscaled residual secant moduli to those corresponding to the initial state E M (which are the same for all direction): the elastic region of the principal test, those variables are null and then tend to increase along with the specimen failure state. Herein, the secant moduli upscaled values ẼM ax and ẼM tr yield the damage variables d ax and d tr (equation 64). Results for tension are given on Figure 16(b) where the two damage variables are plotted in terms of the macroscopic axial strain of the principal calculation. It can be observed that the axial damage d ax is growing faster than that of the transversal one d tr . Moreover, d ax reaches a value of ≈ 0.85 which corresponds to a highly damaged state when d tr hardly

Figure 16 :

 16 Figure 16: Axial and transversal damage variables evolutions

Figure 17 :

 17 Figure 17: Evolution of the residual tensile strength during the failure process

Table 3 :

 3 Influence of elastic properties of aggregates on compression and tension strengths

Table 1 :

 1 Meso-scale material characteristics of each phase and interface

	2 ]

Table 2 :

 2 Macroscopic upscaled material properties for both tension and compression

	Compression	37.8	0.195	39.2	1.2 10 -3	50
	Tension	37.8	0.195	3.9	1.2 10 -4	0.5

Table 3 :

 3 Influence of elastic properties of aggregates on compression and tension strengths

Matrix notation

Emmanuel Roubin a , Alexis Vallade c , Nathan Benkemoun b, * , Jean-Baptiste Colliat c a LMT-Cachan,

In the case of weak discontinuity within the element, the interface is always defined.

Hence n is known prior to any mechanical calculation. Thus, the traction vector T can be defined before the localization. In order to represent the interface orientation as the weakest direction, the equivalent stress for localization is defined as the projection of the traction vector on it:

with

On the contrary, if no material discontinuity can define an interface, strong discontinuity appears with stress state consideration. One can note that in this case, a constant stress tensor can be given for the whole element since it is supposed free of material discontinuity.

Its orientation is defined by the principal direction of the stress tensor. If σ I is its eigenvalue then:

When localization occurs (Φ l = 0), the corresponding eigenvector n I is recorded and set as the interface orientation: n ← n localization

I

. It is assumed that its value remains constant through time. Afterwards, the traction vector is defined by:

which follows its previous definition with σ = σ = σ.

In both cases, after localization the discontinuity surface and its orientation n are defined.

In order to model the same failure mechanism whether an interface element is considered or not, the opening criterion Φ o is assumed to be identical. The equivalent stress is taken to be the projection of the traction vector on n:

Furthermore, this definition of the equivalent stress leads to a very simple written expression of the equivalent stress derivative:

evolution are directly linked to the local fracture energy and tensile strength leading to a model with only two meaningful non-linear parameters. These two enhancements make the model a relevant and efficient tool to represent the failure mechanisms in a continuum region (mortar matrix cracking) as well as on its interfaces (debonding).

Numerical applications of the model have been shown by modeling concrete specimens mainly under uniaxial loadings. The massive use of E-FEM in this case, which consists in the introduction of a high number of strongly enhanced elements (typically more than hundreds of thousands), has lead to complex macroscopic crack patterns. Moreover it has revealed some emergent macroscopic responses that exhibit several features such as asymmetry of the tension-compression stress-strain relationship, which are typical of concrete. However, even though the implementation of full kinematics presented here is a major improvement compared to simpler FE (see Benkemoun et al. (2010)), several weak points are worth noticing regarding the macroscopic responses. They are now discussed.

First, a low ratio of compression to tension strengths and a too brittle behavior in compression are still observed. To the authors point of view, the main reason for that is the choice of rather simple failure mechanisms and cracks opening mode at meso-scale. Hence future work shall include for example, mode II or I+II crack opening. Furthermore, from the authors experience, it has been noticed that high fraction volumes leads to higher ratio of compression to tension strengths. Unfortunately, increasing the fraction volume leads to the need of representing thinner heterogeneity (and therefore finer meshes), thus leading to unreasonable time consuming computations (a computation as presented in this paper takes about 24h of CPU time). Second, as it can be seen in Figure 7, some locking effect can be observed for both tension and compression curves, which exhibit, for large strain, an asymptotic behavior that does not correspond to zero macroscopic force . As a matter of fact, this issue can directly be linked with the two problematics addressed just above. Indeed, it has been said that the orientation of the strong discontinuity at the interface are set to be equal to the weak discontinuities, thus being pre-defined by the heterogeneities (and not based on the principal stress directions). And since the failure is governed by a mode I opening mechanism, stress locking often occurs when the displacement