Nathan Benkemoun

Rachel Gelet

Emmanuel Roubin

Jean-Baptiste Colliat

Poro-elastic two-phase material modeling: theoretical formulation and Embedded Finite Element Method implementation

Keywords: E-FEM, EAS method, Weak Discontinuity, Poroelasticity, Darcy Law, Return Mapping, Operator Split, Static Condensation

This paper presents the formulation of finite element methods for the numerical modeling of a poro-elastic two-phase (aggregates/mixture phase) solid. The displacement and pressure fields are decomposed, following the EAS method, into a regular part and an enhanced part. This leads to discontinuous strain and pressure gradient fields allowing to capture the jump in mechanical and hydrical properties passing through the interface between the aggregates and the mixture phase. All these enhanced fields are treated in the context of the Embedded Finite Element Method through a local enhancement of the finite element interpolations as these jumps appear. The local character of these interpolations lead after a static condensation of the enhanced fields to a problem exhibiting the same structure as common poro-elastic finite element models, but incorporating now the mechanical and hydrical properties of a two-phase solid.

Introduction 2. Introduction

Interest has significantly increased over the last decades in the use of Enriched Finite Element Methods such as the Generalized Finite Element Method (G-FEM, [START_REF] Strouboulis | The design and analysis of the Generalized Finite Element Method[END_REF], [START_REF] Strouboulis | The generalized finite element method: an example of its implementation and illustration of its performance[END_REF]), the eXtended Finite Element Method (X-FEM, [START_REF] Belytschko | Arbitrary discontinuities in finite elements[END_REF], [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], [START_REF] Mariani | Extended finite element method for quasi-brittle fracture[END_REF]) and the Embedded Finite Element Method (E-FEM, [START_REF] Armero | An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids[END_REF], [START_REF] Oliver | Modeling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 2: numerical simulation[END_REF], [START_REF] Borja | A finite element model for strain localization analysis of strongly discontinuous fields based on standard galerkin approximation[END_REF], [START_REF] Benkemoun | Modeling heterogeneous materials failure: 3D meso-scale models with embedded discontinuities[END_REF]). These methods all seek to enrich the basis of the classical Finite Element shape functions by the means of specific enhancements in the displacement/strain fields for mechanical problems and/or in the pressure/pressure gradient, temperature/temperature gradient, concentration/concentration gradient fields for multiphysics problems for instance. It results in higher-performances finite elements capable of capturing physical discontinuity such as cracks, fractures, slip lines and material interfaces.

In the literature, the aforementioned enhancements are in general sorted into two categories:

• strong enhancements -a jump within the displacement (and/or the pressure, the temperature, the concentration,...) fields is introduced by the means of a so-called strong discontinuity ( [START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate independent inelastic solids[END_REF]);

• weak enhancements -a jump within the strain (and/or the pressure gradient, the temperature gradient, the concentration gradient,...) fields is introduced by the means of a so-called weak discontinuity ( [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF]).

On the one hand, due to the flexibility of their numerical implementation and the increase of their numerical efficiency, these tools allow the possibility to treat a broad variety of Civil Engineering problems in reasonable computational time. On the other hand, in many Civil Engineering applications, one has to deal with continuum mechanical problems which cannot be uniquely classified in the context of either "solid mechanics" or "fluid problems".

In general, these applications fall into the categories of porous media (see [12], [START_REF] Wilmanski | A thermodynamic model of compressible porous materials with the balance equation of porosity[END_REF], [START_REF] Schrefler | The finite element method in the static and dynamic deformation and consolidation of porous media[END_REF] for instance). There is a broad variety of problems falling into this category including: geomaterials problems (soil, rock, concrete,...), dealing with the coupled solid deformation and pore-fluid flow behavior of liquid and/or gas into a porous solid skeleton. Thus, modeling these kinds of problems in the context of the aforementioned Enriched Finite Element Methods seems to be an appealing approach. Several applications concerning these coupled problems in porous media referenced in recent literature are worth noting. [START_REF] Alfaiate | On the use of strong discontinuity formulations for the modeling of preferential moisture uptake in fractured porous media[END_REF] develop a model of moisture transport in a fractured porous media, [START_REF] Jourdain | Upscaling permeability for fractured concrete: meso-macro numerical approach coupled to strong discontinuities[END_REF] present a two-scale (mesomacro) model of gas transport in a fractured concrete, [START_REF] Ngo | Model for localized failure with thermo-plastic coupling: Theoretical formulation and ED-FEM implementation[END_REF] a coupled thermo-mechanical failure model in porous construction materials, [START_REF] Callari | Finite element methods for the analysis of strong discontinuities in coupled poro-plastic media[END_REF] and [START_REF] Callari | Strong discontinuities in partially saturated poroplastic solids[END_REF] a fully coupled poro-plastic model in the presence of shear bands in a fully and partially saturated porous solid, [START_REF] Réthoré | A two-scale approach for fluid flow in fractured porous media[END_REF] derive a two-scale (micro-macro) approach for fluid flow in a fractured porous media, [START_REF] De Borst | Discrete vs smeared crack models for concrete fracture: bridging the gap[END_REF] present a coupled model between water and chloride ion transport and stress driven failure and [START_REF] Mohammadnejad | An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model[END_REF] simulate an hydraulically-driven fracture in a porous media. Finally, a classification of this literature is presented in Table 1. It summarizes the models in relation with the performed enhancements in the fields of interest. We denote by (S) a strong enhancement and by (W) a weak one. We call u, p f , q f and T the displacement field, the pressure field, the mass flux and the temperature, respectively .

Models in the literature

Fields of interest Resulting type of enhancement Regarding this literature, all these models, excepted ( [START_REF] Jourdain | Upscaling permeability for fractured concrete: meso-macro numerical approach coupled to strong discontinuities[END_REF]), are used at the macroscopic scale, thus considering the material (cementitious one in general) as an homogeneous medium i.e. a one-phase material with mechanical and permeability or thermal properties. Indeed the explicit representation of heterogeneities (shapes, distribution, sizes and behavior) present at fine scales and having a significant role ( [START_REF] Maruyama | Strain and crack distribution in concrete during drying[END_REF], [START_REF] Yaman | Active and non-active porosity in concrete. part I: experimental evidence[END_REF]) in most of the physical phenomenon (degradation, cracking, transport,...) observed at the macroscopic scale is not taken into account. In this paper, we propose to develop a model at a relevant scale for a reasonable cost in computation time but fine enough to take into account these heterogeneities: the development of a mesoscopic scale (eg. millimeter scale for concrete) model for a fully saturated porous media is proposed. At this scale for concrete, heterogeneities can be considered as aggregates inclusions embedded into a mixture phase and this mixture phase can be seen as a mortar matrix with pores fully saturated by water.

u p f q f T [15] (S) (S)
In order to explicitly take into account these aggregates (shapes and geometry), the authors turn to a method referred as non-adapted meshing process ( [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF]). This consists in a unique homogeneous mesh whose nodes are placed independently from the morphology of the aggregates. A significant amount of computation time is saved at this stage. However this implies the presence of finite elements split into two parts each having different hydromechanical properties. Fig. 1 shows the whole domain Ω and some element split by Γ d into two parts (Ω ⊕ e and Ω ⊖ e ). On Ω ⊕ e , we consider E ⊕ , ν ⊕ , κ ⊕ and on Ω ⊖ e E ⊖ , ν ⊖ ,κ ⊖ which are respectively the Young modulus, the Poisson ratio and the permeability in each parts. In order to take into account this special kinematics in the classical linear elements, a weak enhancement is introduced both in the strain field and the pressure gradient field leading to a W/W enhancement. This paper presents the methodology to introduce this W/W enhancement both in the mechanical problem and the fluid problem in the context of the E-FEM, leading to a mesoscopic poro-elastic fully saturated two-phase model.

It is as follows : In Section 2, the mechanical problem is presented. It is based upon a Hu-Washizu formulation ( [START_REF] Washizu | Variational methods in elasticity and plasticity, 3rd Edition[END_REF]). After that the form of the enhanced fields (enhanced strain field and pressure gradient field) is detailed; their introduction in the strain and the pressure gradient fields is based upon the Enhanced Assumed Strain method ( [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF]). Finally the discretized form of this mechanical problem is obtained as a result of the Finite Element Method. In Section 3, the fluid problem to be solved is introduced. After a little work, the discretized form of this problem is obtained. In Section 4, the solving process to solve this coupled hydro-mechanical problem is presented : it lies on two levels of solving a local one for the enhanced unknowns ( [START_REF] Simo | Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems[END_REF]) and a global one for the regular parts of the displacement and pressure fields. In this section, we first remind the mechanical problem to be solved. Then we show how to adapt it in the context of a porous media.

Ω Ω ⊕ e Γ d Ω ⊖ e n Ω e Γ t t Γ q qf b

Governing equations

As mentioned in the introduction of this paper, some elements can be split into two parts, each part having different mechanical properties. This implies the introduction of a weak discontinuity ( [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF]) in the kinematics of the strain field. It is achieved with the EAS method ( [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF]) considering both the actual and virtual strain fields decomposed into a regular part and an enhanced part. This gives for the actual strain field

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

γ = ∇ s η regular + γ enhanced , (2)
where ∇ s ū (∇ s η) is the symmetric gradient of the actual (virtual) displacement field. As in ( [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF]), we refer to ε and γ as the enhanced parts of the strain fields. The notation (•) refers to weak discontinuity.

Introducing (1) and ( 2) in the Hu-Washizu formulation ( [START_REF] Washizu | Variational methods in elasticity and plasticity, 3rd Edition[END_REF]), one gets a modified Hu-Washizu formulation (see [START_REF] Roubin | Multi-scale failure of heterogeneous materials: a double kinematics enhancement for EFEM method[END_REF] for details) to be solved :

85 Ω ∇ s η : σ dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, ( 3a
) Ωe γ : σ dΩ = 0, (e = 1, 2, ..., n ⊕,⊖ elm ) ( 3b
)
where σ is the total stress tensor, b the imposed volumetric body forces and t the tractions applied at the boundary portion Γ t . We note n ⊕,⊖ elm the number of elements split into two parts Ω ⊕ e and Ω ⊖ e . Remark : Note that the first equation of the problem is global. It has to be solved for all the elements. On the contrary, the second one is local. It has to be solved only for elements 90 split into two parts. This remark is very important for the resolution process described latter.

In the context of saturated porous media ( [START_REF] Tamagnini | Computational Soil Dynamics, aLERT School : Soil -Structure Interaction[END_REF]), the total stress tensor σ is written such

as σ = σ -p f 1, (4)
where σ is the effective stress tensor and p f the fluid pressure (positive in compression).

Combining (3a), (3b) and (4) one obtains

Ω ∇ s η : σ dΩ - Ω p f ∇ • η dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, ( 5a
) Ωe γ : (σ -p f 1) dΩ = 0, (e = 1, 2, ..., n ⊕,⊖ elm ). ( 5b
)
Considering now the fluid problem, some elements can be split into two parts, each part having different properties of permeability. This implies, as for the mechanical problem, the introduction of a weak discontinuity in the kinematics of the pressure field. As for the mechanical problem, the symbol (•) refers to the regular part of the pressure and the symbol (•) to the enhanced part (corresponding to the weak discontinuity) thus the actual pressure p f is such

p f = pf regular + pf enhanced . (6)
Combining (5a), (5b), [START_REF] Armero | An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids[END_REF] and developing terms, one gets

Ω ∇ s η : σ dΩ - Ω pf ∇ • η dΩ - Ω pf ∇ • η dΩ - Ω η • ρb dΩ - Γt η • t d∂Ω = 0, (7a) Ωe γ : σ dΩ - Ωe γ : (p f 1) dΩ - Ωe γ : (p f 1) dΩ = 0, (e = 1, 2, ..., n ⊕,⊖ elm ). ( 7b
)
This represents the mechanical system to be solved. Nevertheless the form of the enhanced fields has still to be explicited. This is the focus of the next subsection.

FE discretization of the fields

In this subsection, we introduce the form of the enhanced fields present in the mechanical problem. Both continuous and discretized forms are introduced.

Continuous form of the enhanced fields

Following the idea developed in [START_REF] Markovič | Modélisation multi-échelles de structures hétérogènes aux comportements anélastiques non-linéaires[END_REF] and [START_REF] Roubin | Multi-scale failure of heterogeneous materials: a double kinematics enhancement for EFEM method[END_REF], the form of the enhanced displacement field is assumed to be

ũ = Θ n • (x -ξ) ([ε] n n + [ε] m m + [ε] t t) with Θ =    Θ ⊕ ∀ x ∈ Ω ⊕ e Θ ⊖ ∀ x ∈ Ω ⊖ e , (8)
where ξ represents the position of the discontinuity Γ d , Θ an unidentified shape function and Applying the symmetric gradient operator to ũ, one gets

([ε] n , [ε] m , [ε] t ) three
ε = Θ [ε] n n ⊗ n + [ε] m 2 (n ⊗ m) s + [ε] t 2 (n ⊗ t) s . (9)
Details are given in Appendix A.

We now use the same reasoning to obtain the form of ∇p f considering the fact that pf is a scalar field whereas ũ was a vectorial field.

We write pf under the following form

pf = Θ n • (x -ξ) (p h f ), (10)
where ξ represents the position of Γ d and Θ an unidentified shape function. The product

n • (x -ξ) is called the distance function.
Note that when x is equal to ξ (in other words 110 when we are on the discontinuity Γ d ), pf is equal to zero thus the pressure is continuous through the discontinuity. ph f is a scalar value computed during the resolution process only for the elements split into two parts. Thus, following ( [START_REF] De Borst | Some observations on embedded discontinuity models[END_REF]), the gradient of pf is given by

∂ pf ∂x 1 = ∂ pf ∂(n T k .(x -ξ) k ) ∂(n T k .(x -ξ) k ) ∂x 1 = Θp h f n 1 , (11)
∂ pf ∂x 2 = ∂ pf ∂(n T k .(x -ξ) k ) ∂(n T k .(x -ξ) k ) ∂x 2 = Θp h f n 2 , (12)
∂ pf ∂x 3 = ∂ pf ∂(n T k .(x -ξ) k ) ∂(n T k .(x -ξ) k ) ∂x 3 = Θp h f n 3 . (13)
Finally, one gets the form of ∇p f such as

∇p f = Θp h f      n 1 n 2 n 3      = Θp h f n. (14)
In the next subsection, we turn to the discretized form of ε and ∇p f suitable for a FE approach.

115

Discretized form of the enhanced fields

Following the idea presented in ( [START_REF] Roubin | Multi-scale failure of heterogeneous materials: a double kinematics enhancement for EFEM method[END_REF]), ε is written under the form

ε = G u w [|ε|], ( 15
)
where ε is written under the Voigt notation and [|ε|] is a vector containing (

[ε] n , [ε] t , [ε] m ).
The form of G u w is explicitly known

G u w =    G u,⊕ w = Θ ⊕ H u w = V ⊖ V H u w in Ω ⊕ e G u,⊖ w = Θ ⊖ H u w = -V ⊕ V H u w in Ω ⊖ e ( 16
)
with

H u w =               n 2 1 n 1 m 1 n 1 t 1 n 2 2 n 2 m 2 n 2 t 2 n 2 3 n 3 m 3 n 3 t 3 2n 1 n 2 n 1 m 2 + n 2 m 1 n 1 t 2 + n 2 t 1 2n 2 n 3 n 2 m 3 + n 3 m 2 n 2 t 3 + n 3 t 2 2n 1 n 3 n 1 m 3 + n 3 m 1 n 1 t 3 + n 3 t 1               . ( 17
)
and

Θ ⊕ = V ⊖ /V , Θ ⊖ = -V ⊕ /V where V, V ⊕ and V ⊖ the volumes of Ω e , Ω ⊕ e and Ω ⊖ e respectively.
By analogy with the enhanced strain field, we write

∇p f = G p w ph f , (18)
where ∇p f is a column vector and ph f a scalar value.

120

The form of G p w is explicitly known

G p w =    G p,⊕ w = Θ ⊕ H p w = V ⊖ V H p w in Ω ⊕ e G p,⊖ w = Θ ⊖ H p w = -V ⊕ V H p w in Ω ⊖ e (19)
with

H p w =      n 1 n 2 n 3      . ( 20
)
Finally this gives for the enhanced pressure field

pf = M p w ph f , (21)
where M p w is a scalar function and ph f a scalar value.

The form of M p w is explicitly known

M p w =    M p,⊕ w = Θ ⊕ n • (x ⊕ -ξ) = V ⊖ V n • (x ⊕ -ξ) in Ω ⊕ e M p,⊖ w = Θ ⊖ n • (x ⊖ -ξ) = -V ⊕ V n • (x ⊖ -ξ) in Ω ⊖ e . ( 22
)
Having the discretized form of the enhanced fields in hand, we now turn to the discretized form of the regular fields.

Discretized form of the regular fields

For the gradient of the regulart part of the virtual strain ∇ s η, we have

∇ s η = B u ηh , (23)
where B u is the classical interpolation matrix (derivatives) for a tetrahedron elements and 125 ηh is a column vector containing the virtual nodal displacement unknowns. Note that ∇ s η is again under the Voigt convention.

For the divergence of the regulart part of the virtual displacement ∇ • η, we get

∇ • η = [1 1 1 0 0 0]∇ s η = m∇ s η = mB u ηh . ( 24
)
For the regular part of the actual pressure field pf , we have

pf = N p pf h , (25)
where N p is a row vector containing the shape functions for the regular part of the actual pressure field and pf h is a column vector containing the actual nodal pressure unknowns.

130

Finally we obtain

∇ s η = B u ηh + G u w [|γ|] , (26)
p f = N p pf h + M p w ph f , (27)
∇p f = B p pf h + G p w ph f , (28) 10
where B p is a matrix containing the derivatives of the shape functions N p . Note that in the expression of ∇ s η, the enhanced part γ has the same form as ε where

[|γ|] = ([γ] n , [γ] m , [γ] t ).

135

Having the expression of ∇ s η, p f and ∇p f in a suitable form for a FE approach, we can now introduce them into the mechanical problem (7a) and (7b) and obtain the discretized mechanical problem to be solved.

FE discretization of (7a) and (7b)

Combining ( 24), ( 26), ( 27) and (7a) and considering the behaviour law for the effective stress such as

σ(ū h , [|ε|]) =    σ ⊕ = C ⊕ ε = C ⊕ B u ūh + G u,⊕ w [|ε|] if x ∈ Ω ⊕ e σ ⊖ = C ⊖ ε = C ⊖ B u ūh + G u,⊖ w [|ε|] if x ∈ Ω ⊖ e , (29)
yield

Ω ⊕ B u,T C ⊕ B u dΩ + Ω ⊖ B u,T C ⊖ B u dΩ ūh + Ω ⊕ B u,T C ⊕ G u,⊕ w dΩ + Ω ⊖ B u,T C ⊖ G u,⊖ w dΩ [|ε|] - Ω B u,T m T N p dΩ pf h - Ω ⊕ B u,T m T M p,⊕ w dΩ + Ω ⊖ B u,T m T M p,⊖ w dΩ ph f - Ω N u,T ρb dΩ - Γt N u,T t d∂Ω = 0. ( 30
)
where C ⊕ and C ⊖ are the elastic tensors for the material in Ω ⊕ and Ω ⊖ .

140

Under matricial format, one gets

K uu ūh + K uw [|ε|] -K u p pf h -K up ph f = f ext , (31)
where

K uu = Ω ⊕ B u,T C ⊕ B u dΩ + Ω ⊖ B u,T C ⊖ B u dΩ, (32)
K uw = Ω ⊕ B u,T C ⊕ G u,⊕ w dΩ + Ω ⊖ B u,T C ⊖ G u,⊖ w dΩ, (33)
K u p = Ω B u,T m T N p dΩ, (34)
K up = Ω ⊕ B u,T m T M p,⊕ w dΩ + Ω ⊖ B u,T m T M p,⊖ w dΩ, (35)
f ext = Ω N u,T ρb dΩ + Γt N u,T t d∂Ω. (36)
Combining ( 24), ( 26), ( 27) and (7b) yields

Ω ⊕ e G u,⊕,T w C ⊕ B u dΩ + Ω ⊖ e G u,⊖,T w C ⊖ B u dΩ ūh + Ω ⊕ e G u,⊕,T w C ⊕ G u,⊕ w dΩ + Ω ⊖ e G u,⊖,T w C ⊖ G u,⊖ w dΩ [|ε|] - Ω ⊕ e G u,⊕,T w m T N p dΩ + Ω ⊖ e G u,⊖,T w m T N p dΩ pf h - Ω ⊕ e G u,⊕,T w m T M p,⊕ w dΩ + Ω ⊖ e G u,⊖,T w m T M p,⊖ w dΩ ph f = 0. ( 37
)
Under a matricial format, one gets

K uw T ūh + K ww [|ε|] -K w p pf h -K wp ph f = 0, (38)
where

K ww = Ω ⊕ e G u,⊕,T w C ⊕ G u,⊕ w dΩ + Ω ⊖ e G u,⊖,T w C ⊖ G u,⊖ w dΩ, (39)
K w p = Ω ⊕ e G u,⊕,T w m T N p dΩ + Ω ⊖ e G u,⊖,T w m T N p dΩ, (40)
K wp = Ω ⊕ e G u,⊕,T w m T M p,⊕ w dΩ + Ω ⊖ e G u,⊖,T w m T M p,⊖ w dΩ. (41)
4. Balance of mass equation : two-phase mixture formulation

Governing equations

In the present work, the model is formulated in a Biot thermodynamic framework ( [START_REF] Biot | General theory of three-dimensional consolidation[END_REF], [12]) employing the variables of solid-skeleton displacements and fluid pressures. Following the development made in ( [START_REF] Tamagnini | Computational Soil Dynamics, aLERT School : Soil -Structure Interaction[END_REF]), the weak form of the mass balance equation for the twophase mixture (mortar matrix fully saturated by water) can be written

Ω ψ ṗf Q dΩ+ Ω ψ∇• us dΩ+ Γ ψv•n d∂Ω+ 1 gρ w Ω ∇ψ•κ∇p f dΩ- 1 g Ω ∇ψ•κb dΩ = 0, ( 42
)
where ψ is the virtual pressure field, vanishing on Γ p where the pressure p f is imposed, Q the ratio between the porosity n and the pore water compressibility K w and us the solid-skeleton velocity. We note the part of the boundary Γ q where the external normal flow v • n q for the 150 outward unit normal n q to d∂Ω, is -q f . Finally we consider Γ d the weak discontinuity surface splitting a given element into two-subdomains Ω ⊕ and Ω ⊖ where in terms of fluxes we

have v•n ⊕ d = -q ⊕ f and v•n ⊖ d = -q ⊖ f . Note that n ⊕ d is pointing toward Ω ⊕ and n ⊖ d toward Ω ⊖ .
With these considerations in hands, (42) can be rewritten

155 Ω ψ ṗf Q dΩ + Ω ψ∇ • us dΩ + Γq ψv • n q d∂Ω + Γ d ψv • n d d∂Ω + 1 gρ w Ω ∇ψ • κ∇p f dΩ - 1 g Ω ∇ψ • κb dΩ = 0. ( 43
)
The part on Γ d can be decomposed into

Γ d ψv • n d d∂Ω = Γ ⊕ d ψv • n ⊕ d d∂Ω + Γ ⊖ d ψv • n ⊖ d d∂Ω = - Γ ⊕ d ψ q⊕ f d∂Ω - Γ ⊖ d ψ q⊖ f d∂Ω. (44)
In the case of a weak discontinuity, we consider that no flux is lost passing through the discontinuity Γ d thus q⊕ f + q⊖ f = 0. This implies that

- Γ ⊕ d ψ q⊕ f d∂Ω - Γ ⊖ d ψ q⊖ f d∂Ω = 0, (45)
and so

Γ d ψv • n d d∂Ω = 0. ( 46
)
Finally combining (43) and ( 46), one has

Ω ψ ṗf Q dΩ+ Ω ψ∇• us dΩ- Γq ψ qf d∂Ω+ 1 gρ w Ω ∇ψ•κ∇p f dΩ- 1 g Ω ∇ψ•κb dΩ = 0. ( 47
)
Following the form of p f , the virtual pressure field ψ is decomposed such as

ψ = ψ regular + ψ enhanced . ( 48
)
The shape functions are the same as for the actual pressure field thus

ψ = N p ψh + M p w ψh , (49)
∇ψ = B p ψh + G p w ψh , (50)
where B p is a matrix containing the derivatives of the shape functions N p .

Having in hands the discretized form of the fields present in (47), we now turn to its discretization.

FE discretization of (47)

To carry out the time integration in (47), a backward Euler method is adopted:

( d(•) dt ) n+1 = (•) n+1 -(•) n ∆t = ∆(•) ∆t ( 51
)
where ∆t is the time increment, while (•) n+1 and (•) n denote the unknowns at time step n + 1 and n, respectively.

Combining (47) with ( 49), ( 50), ( 27) and ( 28), considering the fact that ∇ • us = m ε and using (51) for the time derivatives, one gets ∀ ψh,T

1 Q Ω 1 ∆t N p,T N p dΩ ∆ pf h + 1 gρ w Ω ⊕ B p,T • κ ⊕ B p dΩ + 1 gρ w Ω ⊖ B p,T • κ ⊖ B p dΩ pf h + 1 Q Ω ⊕ 1 ∆t N p,T M p,⊕ w dΩ + 1 Q Ω ⊖ 1 ∆t N p,T M p,⊖ w dΩ ∆p h f + 1 gρ w Ω ⊕ B p,T • κ ⊕ G p,⊕ w dΩ + 1 gρ w Ω ⊖ B p,T • κ ⊖ G p,⊖ w dΩ ph f + Ω 1 ∆t N p,T mB u dΩ ∆ū h + Ω ⊕ 1 ∆t N p,T mG u,⊕ w dΩ + Ω ⊖ 1 ∆t N p,T mG u,⊖ w dΩ ∆ [|ε|] - Γq N p,T qf d∂Ω - 1 g Ω B p,T • κb dΩ = 0, (52)
and ∀ ψh,T

1 Q Ω ⊕ e 1 ∆t M p,⊕ w N p dΩ + 1 Q Ω ⊖ e 1 ∆t M p,⊖ w N p dΩ ∆ pf h + 1 gρ w Ω ⊕ e G p,⊕,T w • κ ⊕ B p dΩ + 1 gρ w Ω ⊖ e G p,⊖,T w • κ ⊖ B p dΩ pf h + 1 Q Ω ⊕ e 1 ∆t M p,⊕,T w M p,⊕ w dΩ + 1 Q Ω ⊖ e 1 ∆t M p,⊖,T w M p,⊖ w dΩ ∆p h f + 1 gρ w Ω ⊕ e G p,⊕,T w • κ ⊕ G p,⊕ w dΩ + 1 gρ w Ω ⊖ e G p,⊖,T w • κ ⊖ G p,⊖ w dΩ ph f + Ω ⊕ e 1 ∆t M p,⊕,T w mB u dΩ + Ω ⊖ e 1 ∆t M p,⊖,T w mB u dΩ ∆ū h + Ω ⊕ e 1 ∆t M p,⊕,T w mG u,⊕ w dΩ + Ω ⊖ e 1 ∆t M p,⊖,T w mG u,⊖ w dΩ ∆ [|ε|] - Γq M p,T w qf d∂Ω - 1 g Ωe G p,T w • κb dΩ = 0, (e = 1, 2, ..., n ⊕,⊖ elm ). ( 53
)
Under a matricial format, (52) and ( 53) are rewritten such as

1 ∆t K T u p ∆ū h + 1 ∆t K T w p ∆ [|ε|] + 1 ∆t M ∆ pf h + H pf h + 1 ∆t L ∆p h f + G ph f = f 1ext , (54)
1 ∆t K T up ∆ū h + 1 ∆t K T wp ∆ [|ε|] + 1 ∆t L T ∆ pf h + G T pf h + 1 ∆t C ∆p h f + B ph f = f 2 . (55)
where

M = 1 Q Ω N p,T N p dΩ, (56)
H = 1 gρ w Ω ⊕ B p,T • κ ⊕ B p dΩ + 1 gρ w Ω ⊖ B p,T • κ ⊖ B p dΩ, (57)
L = 1 Q Ω ⊕ N p,T M p,⊕ w dΩ + 1 Q Ω ⊖ N p,T M p,⊖ w dΩ, (58)
G = 1 gρ w Ω ⊕ B p,T • κ ⊕ G p,⊕ w dΩ + 1 gρ w Ω ⊖ B p,T • κ ⊖ G p,⊖ w dΩ, (59)
f 1ext = Γq N p,T qf d∂Ω + 1 g Ω B p,T • κb dΩ, (60)
C = 1 Q Ω ⊕ M p,⊕,T w M p,⊕ w dΩ + 1 Q Ω ⊖ M p,⊖,T w M p,⊖ w dΩ, (61)
B = 1 gρ w Ω ⊕ G p,⊕,T w • κ ⊕ G p,⊕ w dΩ + 1 gρ w Ω ⊖ G p,⊖,T w • κ ⊖ G p,⊖ w dΩ, (62)
f 2 = Γq M p,T w qf d∂Ω + 1 g Ω G p,T w • κb dΩ. ( 63
)
Remark : As for the mechanical problem, note that the first equation of the problem (54) is global. It has to be solved for all the elements. On the contrary, the second one is local.

It has to be solved only for elements split into two parts. This remark is very important for the resolution process described latter.

Hydro-mechanical coupled problem : Global-local solving process

System to be solved

The hydro-mechanical coupled problem to be solved is obtained by combining ( 31), ( 38), ( 54) and ( 55). This leads to

K uu ūh + K uw [|ε|] -K u p pf h -K up ph f = f ext , (64a)
K uw T ūh + K ww [|ε|] -K w p pf h -K wp ph f = 0, (64b) 1 ∆t K T u p ∆ū h + 1 ∆t K T w p ∆ [|ε|] + 1 ∆t M ∆ pf h + H pf h + 1 ∆t L ∆p h f + G ph f = f 1 ext , (64c)
1 ∆t K T up ∆ū h + 1 ∆t K T wp ∆ [|ε|] + 1 ∆t L T ∆ pf h + G T pf h + 1 ∆t C ∆p h f + B ph f = f 2 , (64d)
where the unknowns are ūh , pf h , [|ε|] and ph f . The solving process is based upon the key idea of [START_REF] Simo | Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems[END_REF]. Since the residual equations (64b) and (64d) hold only for each element e = 1, 2, ..., n ⊕,⊖ elm split into two parts, the most effective solution is to condense out the internal element modes [|ε|] and ph f to arrive at an equivalent displacement model involving only the nodal displacement vector ūh at the global level, where the nodal pressure vector pf h is then post-calculated.

The problem is therefore split into two levels of equations : (1) the global one corresponding to (64a) and (64c), ( 2) the local one corresponding to (64b) and (64d). The term global means that the equations are applied to all the elements of the problem. The term local means that the equations have to be solved only for the elements e = 1, 2, ..., n ⊕,⊖ elm split into two parts. Due to the fact that the E-FEM method consists in adding new unknowns ([|ε|] and ph f in this case) at the element level, the value of these unknowns is discontinuous passing through the boundaries of the elements and so this local property allows a static condensation ( [START_REF] Wilson | The static condensation algorithm[END_REF]) once these new unknowns are solved.

In the rest of this paper and technically speaking, these two levels of equationsglobal and local -correspond to indices called respectively k and l. For a given iterate

∆ū h (k) n+1 , ∆ pf h (k) n+1
, we perform a local sub-iteration on l to solve exactly the non-linear equations (64b) and (64d).

In order to solve the system, one has first to linearize it. This is explained is the next four subsections.

Linearization of (64a)

Let's call the left hand side of (64a), f int thus we have at the kth iteration at the time

step n + 1 f int (ū h , [|ε|] , pf h , ph f ) (k) n+1 -f ext n+1 = 0. ( 65
)
The linearization of (65) at the global level corresponding to k yields

f int (ū h , [|ε|] , pf h , ph f ) (k+1) n+1 -f ext n+1 = f int (ū h , [|ε|] , pf h , ph f ) (k) n+1 -f ext n+1 + ∂f int (X) ∂X (k) n+1 ∆X (k+1) n+1 ≃ 0, (66)
where

∂f int (X) ∂X (k) n+1 ∆X (k+1) n+1 = K uu (k) n+1 ∆ū h (k+1) n+1 + K uw (k) n+1 ∆ [|ε|] (k+1) n+1 -K u p (k) n+1 ∆ pf h (k+1) n+1 -K up (k) n+1 ∆p h f (k+1) n+1 . ( 67
)
Finally one gets

K uu (k) n+1 ∆ū h (k+1) n+1 + K uw (k) n+1 ∆ [|ε|] (k+1) n+1 -K u p (k) n+1 ∆ pf h (k+1) n+1 -K up (k) n+1 ∆p h f (k+1) n+1 = -f int (ū h , [|ε|] , pf h , ph f ) (k) n+1 -f ext n+1 . ( 68
)

Linearization of (64c)

Let's call the left hand side of (64c), f 1 int thus we have at the kth iteration at the time

step n + 1 f 1int (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) (k) n+1 -f 1ext n+1 = 0. ( 69
)
The linearization of ( 69) at the global level corresponding to k yields

f 1 int (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) (k+1) n+1 -f 1 ext n+1 = f 1int (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) (k) n+1 -f 1ext n+1 + ∂f 1 int (X) ∂X (k) n+1 ∆X (k+1) n+1 ≃ 0, (70)
where

∂f 1 int (X) ∂X (k) n+1 ∆X (k+1) n+1 = 1 ∆t K T u p (k) n+1 ∆ū h (k+1) n+1 + 1 ∆t K T w p ∆ [|ε|] (k+1) n+1 + [ 1 ∆t M + H] (k) n+1 ∆ pf h (k+1) n+1 + [ 1 ∆t L + G] (k) n+1 ∆p h f (k+1) n+1 . (71)
Finally one obtains

1 ∆t K T u p (k) n+1 ∆ū h (k+1) n+1 + 1 ∆t K T w p ∆ [|ε|] (k+1) n+1 + [ 1 ∆t M + H] (k) n+1 ∆ pf h (k+1) n+1 + [ 1 ∆t L + G] (k) n+1 ∆p h f (k+1) n+1 = -f 1int (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) (k) n+1 -f 1ext n+1 . ( 72
)

Linearization of (64b)

Let's call (64b), h e (ū h , [|ε|] , pf h , ph f ). The linearization of (64b) is performed localy for a given global iterate k. In other word, for a given k, a local sub-iteration on l is done in order to solve exactly (64b).

For a given k and an iterate l, we have

h e (ū h , [|ε|] , pf h , ph f ) (k,l) n+1 = K uw T ūh (k) n+1 + K ww [|ε|] (k,l) n+1 -K w p pf h (k) n+1 -K wp ph f (k,l) n+1 = 0. (73)
This can be rewritten with the increment of the unknowns under the form

K uw T ūh n + ∆ū h (k) n+1 + K ww [|ε|] n + ∆ [|ε|] (k,l) n+1 -K w p pf h n + ∆ pf h (k) n+1 -K wp ph f n + ∆p h f (k,l) n+1 = 0. (74)
After collecting terms together, one gets

K uw T ūh n + K ww [|ε|] n -K w p pf h n -K wp ph f n + K uw T ∆ū h (k) n+1 + K ww ∆ [|ε|] (k,l) n+1 -K w p ∆ pf h (k) n+1 -K wp ∆p h f (k,l) n+1 = 0, (75)
or under a pithy form

h e (ū h , [|ε|] , pf h , ph f ) n + K uw T ∆ū h (k) n+1 + K ww ∆ [|ε|] (k,l) n+1 -K w p ∆ pf h (k) n+1 -K wp ∆p h f (k,l) n+1 = 0. (76)
Introducing a trial function

h e (ū h , [|ε|] , pf h , ph f ) trial n+1
in the spirit of a return mapping algorithm ( [START_REF] Simo | Computational Inelasticity[END_REF]) such as

h e (ū h , [|ε|] , pf h , ph f ) trial n+1 = h e (ū h , [|ε|] , pf h , ph f ) n + K uw T ∆ū h (k) n+1 -K w p ∆ pf h (k) n+1 , (77)
one gets

h e (ū h , [|ε|] , pf h , ph f ) trial n+1 + K ww ∆ [|ε|] (k,l) n+1 -K wp ∆p h f (k,l) n+1 = 0. (78)
(78) has to be linearized in relation with l, this yields

h e (ū h , [|ε|] , pf h , ph f ) (k,l+1) n+1 = h e (ū h , [|ε|] , pf h , ph f ) (k,l) n+1 + ∂h e (∆X) ∂∆X (k,l) n+1 δ∆X (k,l+1) n+1 ≃ 0, (79)
where

∂h e (∆X) ∂∆X (k,l) n+1 δ∆X (k,l+1) n+1 = K ww (k,l) n+1 δ∆ [|ε|] (k,l+1) n+1 -K wp (k,l) n+1 δ∆p h f (k,l+1) n+1 . (80)
Finally one gets

K ww (k,l) n+1 δ∆ [|ε|] (k,l+1) n+1 -K wp (k,l) n+1 δ∆p h f (k,l+1) n+1 = -h e (ū h , [|ε|] , pf h , ph f ) (k,l) n+1 . (81)
5.5. Linearization of (64d)

200 Let's call (64d), s e (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ).
The linearization of (64d) is performed localy for a given global iterate k. In other word, for a given k, a local sub-iteration on l is done in order to solve exactly (64d).

For a given k and an iterate l, we have (neglecting f 2 )

s e (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) (k,l) n+1 = 1 ∆t K T up ∆ū h (k) n+1 + 1 ∆t K T wp ∆ [|ε|] (k,l) n+1 + 1 ∆t L T ∆ pf h (k) n+1 + G T pf h (k) n+1 + 1 ∆t C ∆p h f (k,l) n+1 + B ph f (k,l) n+1 = 0. (82)
Following the same reasoning as for (64b), one gets

s e (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) trial n+1 + 1 ∆t K T wp ∆ [|ε|] (k,l) n+1 + [ 1 ∆t C + B] ∆p h f (k,l) n+1 = 0, (83)
where

s e (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) trial n+1 = s e (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) n + 1 ∆t K T up ∆ū h (k) n+1 + [ 1 ∆t L T + G T ] ∆ pf h (k) n+1 . (84)
After linearizing (83), one obtains

1 ∆t K T wp (k,l) n+1 δ∆ [|ε|] (k,l+1) n+1 + [ 1 ∆t C + B] (k,l) n+1 δ∆p h f (k,l+1) n+1 = -s e (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) (k,l) n+1 . ( 85
)
Having the linearization of the coupled system in hands, we now turn to its resolution.

As mentioned above, this is achieved at two levels : a local one and a global one. 20

Algorithm for the local system resolution

The local system to be solved is

K ww (k,l) n+1 δ∆ [|ε|] (k,l+1) n+1 -K wp (k,l) n+1 δ∆p h f (k,l+1) n+1 = -h e (ū h , [|ε|] , pf h , ph f ) (k,l) n+1 , (86a)
1 ∆t K T wp (k,l) n+1 δ∆ [|ε|] (k,l+1) n+1 + [ 1 ∆t C + B] (k,l) n+1 δ∆p h f (k,l+1) n+1 = -s e (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆p h f ) (k,l) n+1 , (86b)
in terms of δ∆p h f and δ∆ [|ε|]. We remind that this system has to be solved only for the elements split into two parts.

210

After a little work, the solution of this system is Once convergence is achieved, ∆ [|ε|] and ∆p h f are known. The next subsection presents the global system resolution strategy in the spirit of [START_REF] Simo | Computational Inelasticity[END_REF]. > tol do Evaluate :

δ∆ [|ε|] (k,l+1) n+1 = -K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp -1 h e + s e [ 1 ∆t C + B] K wp (k,l) n+1 , (87) δ∆p h f (k,l+1) n+1 = - s e [ 1 ∆t C + B] + 1 ∆t[ 1 ∆t C + B] K T wp [K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp ] -1 [h e + s e [ 1 ∆t C + B] K wp ] (k,l) n+1 . (88)
ūh (k) n+1 = ūh n + ∆ū h (k) n+1 , (89)
pf h (k) n+1 = pf h n + ∆ pf h (k) n+1 , (90)
• h e (k,l) n+1 = h e trial n+1 + K ww ∆ [|ε|] (k,l) n+1 -K wp ∆p h f (k,l) n+1 ; • s e (k,l) n+1 = s e trial n+1 + 1 ∆t K T wp ∆ [|ε|] (k,l) n+1 + [ 1 ∆t C + B] ∆p h f (k,l) n+1 .
Compute :

δ∆ [|ε|] (k,l+1) n+1 = -K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp -1 h e + s e [ 1 ∆t C + B] K wp (k,l) n+1 , δ∆p h f (k,l+1) n+1 = - s e [ 1 ∆t C + B] + 1 ∆t[ 1 ∆t C + B] K T wp [K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp ] -1 [h e + s e [ 1 ∆t C + B] K wp ] (k,l) n+1 . Update : ∆ [|ε|] (k,l+1) n+1 = ∆ [|ε|] (k,l) n+1 + δ∆ [|ε|] (k,l+1) n+1 , ∆p h f (k,l+1) n+1 = ∆p h f (k,l) n+1 + δ∆p h f (k,l+1)
(k) n+1 = [|ε|] n + ∆ [|ε|] (k) n+1 , (91) ph
f (k) n+1 = ph f n + ∆p h f (k) n+1 , (92)
be the total enhanced parameters values;

iii) evaluate the internal force vectors

f int (ū h , [|ε|] , pf h , ph f ) (k) n+1 and f 1int (∆ū h , ∆ [|ε|] , pf h , ∆ pf h , ph f , ∆ iv) check convergence: if it is satisfied then (•) (k) n+1
is the solution; otherwise, continue; and; v) determine ∆ū h (step (v)). This is detailed here after where, for sake of clarity, superscript and subscript are omitted.

Once the local system is solved thanks to algorithm (1) , one has

K uw T ∆ū h + K ww ∆ [|ε|] -K w p ∆ pf h -K wp ∆p h f = 0, (93a)
1 ∆t K T up ∆ū h + 1 ∆t K T wp ∆ [|ε|] + [ 1 ∆t L T + G T ] ∆ pf h + [ 1 ∆t C + B] ∆p h f = 0, ( 93b
)
where it is now possible to express ∆ [|ε|] and ∆p h f as a function of ∆ū h and ∆ pf h in order to perform a static condensation at the global level.

This gives the following expressions

∆ [|ε|] = -K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp -1 (K uw T + 1 ∆t[ 1 ∆t C + B] K wp K T up )∆ū h + (-K w p + 1 [ 1 ∆t C + B] K wp [ 1 ∆t L T + G T ])∆ pf h , (94)
∆p h f = - 1 ∆t[ 1 ∆t C + B] K T up + 1 ∆t[ 1 ∆t C + B] K T wp (K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp ) -1 (K uw T + 1 ∆t[ 1 ∆t C + B] K wp K T up ) ∆ū h + - 1 [ 1 ∆t C + B] [ 1 ∆t L T + G T ] + 1 ∆t[ 1 ∆t C + B] K T wp (K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp ) -1 (-K w p + 1 [ 1 ∆t C + B] K wp [ 1 ∆t L T + G T ]) ∆ pf h . (95)
For sake of clarity, we set

X = -K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp -1 (K uw T + 1 ∆t[ 1 ∆t C + B] K wp K T up ) , (96)
Y = -K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp -1 (-K w p + 1 [ 1 ∆t C + B] K wp [ 1 ∆t L T + G T ] . (97)
and

W = - 1 ∆t[ 1 ∆t C + B] K T up + 1 ∆t[ 1 ∆t C + B] K T wp (K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp ) -1 (K uw T + 1 ∆t[ 1 ∆t C + B] K wp K T up ) , (98)
Z = - 1 [ 1 ∆t C + B] [ 1 ∆t L T + G T ] + 1 ∆t[ 1 ∆t C + B] K T wp (K ww + 1 ∆t[ 1 ∆t C + B] K wp K T wp ) -1 (-K w p + 1 [ 1 ∆t C + B] K wp [ 1 ∆t L T + G T ]) . (99)
This gives ∆

[|ε|] = X ∆ū h + Y∆ pf h , ( 100
) ∆p h f = W∆ū h + Z∆ pf h . (101)

Static condensation 230

We now come back to the global equation ( 68) and (72) to perform a static condensation ( [START_REF] Wilson | The static condensation algorithm[END_REF]). In (72), we express ∆ pf h as a function of ∆ū h . This gives, using (100) and (101)

1 ∆t K T u p + 1 ∆t K T w p X + [ 1 ∆t L + G]W ∆ū h + [ 1 ∆t M + H] + 1 ∆t K T w p Y + [ 1 ∆t L + G]Z ∆ pf h = -f 1int -f 1ext . (102)
Finally ∆ pf h can be rewritten in terms of ∆ū h only such as

∆ pf h = -[ 1 ∆t M + H] + 1 ∆t K T w p Y + [ 1 ∆t L + G]Z -1 [f 1 int -f 1 ext ] + [ 1 ∆t K T u p + 1 ∆t K T w p X + [ 1 ∆t L + G]W]∆ū h . (103)
For sake of clarity, we set

U = [ 1 ∆t M + H] + 1 ∆t K T w p Y + [ 1 ∆t L + G]Z , (104)
and

V = [ 1 ∆t K T u p + 1 ∆t K T w p X + [ 1 ∆t L + G]W]. (105)
This gives

∆ pf h = -U -1 [f 1 int -f 1 ext ] + V∆ū h . (106)
Inserting equations (100), ( 101) and ( 106) in (68) gives

K uu ∆ū h + K uw X ∆ū h -Y U -1 [(f 1int -f 1ext ) + V∆ū h ] + K u p U -1 (f 1 int -f 1 ext ) + V∆ū h -K up W∆ū h -ZU -1 [(f 1 int -f 1 ext ) + V∆ū h ] = -f int -f ext . (107)
After collecting terms together and reordering the former expression, one gets

K uu + K uw (X -Y U -1 V) + K u p U -1 V -K up (W -ZU -1 V) ∆ū h = -f int -f ext --K uw Y U -1 + K u p U -1 + K up ZU -1 f 1int -f 1ext . ( 108
)
This can be written under the form

K∆ū h = f . (109)
a weak discontinuity in the gradient of the pressure field as performed in this paper and (2) a discontinuous pressure field (see [START_REF] Larsson | Embedded localization band in undrained soil based on regularized strong discontinuity-theory and FE-analysis[END_REF], [START_REF] De Borst | Discrete vs smeared crack models for concrete fracture: bridging the gap[END_REF]) by introducing a strong discontinuity in the pressure field. By keeping the W/W enhancement already done in this work for capturing heterogeneities, two new enhancements can be achieved: (1) a S+W/W if the pressure is continuous and (2) a S+W/S+W if the pressure is discontinuous. Finally, introducing the coupling of fluid transport with the process of cracking is of main interest in Civil Engineering problems such as CO 2 storage or nuclear industry.

Appendix A. details for ε

In order to obtain ε, we apply the symmetric gradient operator to ũ, this yields

ε = ∇ s (ũ) =      ũ1,1 1
where ũ = (ũ 1 , ũ2 , ũ3 ), x = (x 1 , x 2 , x 3 ) and (• ,1 ), (• ,2 ) and (• ,3 ) represent derivatives in relation with x 1 , x 2 and x 3 .

Following the reasoning developed in ( [START_REF] De Borst | Some observations on embedded discontinuity models[END_REF]), we have for ũ1,1

ũ1,1 = ∂ ũ1 ∂x 1 = ∂ ũ1 ∂(n k .(x -ξ) k ) ∂(n k .(x -ξ) k ) ∂x 1 = ũ′ 1 n 1 = Θ([ε] n n 1 +[ε] m m 1 +[ε] t t 1 )n 1 , (113)
where n k .(x -ξ) k = n 1 .(x -ξ) 1 + n 2 .(x -ξ) 2 + n 3 .(x -ξ) 3 .

And so

∂ ũ1 ∂x 2 = Θ([ε] n n 1 + [ε] m m 1 + [ε] t t 1 )n 2 , (114)
∂ ũ2 ∂x 1 = Θ([ε] n n 2 + [ε] m m 2 + [ε] t t 2 )n 1 , (115)
∂ ũ2 ∂x 2 = Θ([ε] n n 2 + [ε] m m 2 + [ε] t t 2 )n 2 , (116)
∂ ũ3 ∂x 3 = Θ([ε] n n 3 + [ε] m m 3 + [ε] t t 3 )n 3 , (117)
∂ ũ1 ∂x 3 = Θ([ε] n n 1 + [ε] m m 1 + [ε] t t 1 )n 3 , (118)
∂ ũ3 ∂x 1 = Θ([ε] n n 3 + [ε] m m 3 + [ε] t t 3 )n 1 , (119)
∂ ũ2 ∂x 3 = Θ([ε] n n 2 + [ε] m m 2 + [ε] t t 2 )n 3 , (120)
∂ ũ3 ∂x 2 = Θ([ε] n n 3 + [ε] m m 3 + [ε] t t 3 )n 2 . ( 121
)
Developing the former expressions gives

ǫ11 = Θ([ε] n n 2 1 + [ε] m m 1 n 1 + [ε] t t 1 n 1 ), (122) ǫ22
t 3 n 3      . ( 128
)
Finally ε can be rewritten under the form

ε = Θ [ε] n n ⊗ n + [ε] m 2 (n ⊗ m) s + [ε] t 2 (n ⊗ t) s . (129)

Figure 1 :

 1 Figure 1: Illustration of the weak discontinuity in an element used in Ω e : Ω ⊖ e represents the aggregate phase and Ω ⊕ e the mixture phase (mortar matrix with pores fully saturated by water). The interface Γ d illustrates the position of the weak discontinuity

  scalar values computed during the resolution process only for the elements split into two parts . The product n • (xξ) is called the distance function. Note that when x is equal to ξ (in other words when we are on the discontinuity Γ d ), ũ is equal to zero thus the displacement is continuous through the discontinuity.

This local system resolution is performed in the more general case when h e

  in this case algorithm (1) is performed. Intermediate cases exist but they are not detailed at this stage.

215 5 . 7 .

 57 Global system resolution strategyWe let (•) (k) n+1 be the value of a variable (•) at the kth iteration during the load step in [t n , t n+1 ]. Accordingly, i) Let ∆ū h nodal displacement and pressure at the kth iteration and let

Algorithm 1 :

 1 Local procedure for a given time step (n + 1) and a given global kth iteration l tol while (s e ) 2 + (||h e ||) 2 (k,l) n+1

n+1 .l

 n+1 ← l + 1 (next local iteration) ; be the total nodal displacement and pressure; ii) Freeze ∆ū h

,

  set k ← k + 1, and go to step (i).The only step in the solution procedure outlined above which remains to be addressed is determining ∆ū h

Table 1 :

 1 Type of enhancement performed in the literature (note that

	[17]	(S)	(W)	S/W
	[18] and [19]	(S)	(S)	S/S
	[20]	(S)	(W)	S/W
	[21]	(S)	(S)	S/S
	[22]	(S)	(W)	S/W

* or (W) † S/S * or S/W † * stands for the G-FEM using and † for the E-FEM). For instance, a strong enhancement in the mechanical problem (jump in the displacement field) and a weak enhancement in the fluid problem (jump in the pressure gradient field) yield an enhancement labeled as S/W

  = Θ([ε] n n 2 2 + [ε] m m 2 n 2 + [ε] t t 2 n 2 ), (123) ǫ33 = Θ([ε] n n 2 3 + [ε] m m 3 n 3 + [ε] t t 3 n 3 ), (n 1 n 2 + n 2 n 1 ) + [ε] m (m 1 n 2 + m 2 n 1 ) + [ε] t (t 1 n 2 + t 2 n 1 )), (n 1 n 3 + n 3 n 1 ) + [ε] m (m 1 n 3 + m 3 n 1 ) + [ε] t (t 1 n 3 + t 3 n 1 )), (n 2 n 3 + n 3 n 2 ) + [ε] m (m 2 n 3 + m 3 n 2 ) + [ε] t (t 2 n 3 + t 3 n 2 )). (m 1 n 2 + m 2 n 1 ) 1 2 (m 1 n 3 + m 3 n 1 ) 1 2 (m 1 n 2 + m 2 n 1 ) m 2 n 2 1 2 (m 2 n 3 + m 3 n 2 ) 1 2 (m 1 n 3 + m 3 n 1 ) 1 2 (m 2 n 3 + m 3 n 2 ) m 3 n 3 (t 1 n 2 + t 2 n 1 ) 1 2 (t 1 n 3 + t 3 n 1 ) 1 2 (t 1 n 2 + t 2 n 1 ) t 2 n 2 1 2 (t 2 n 3 + t 3 n 2 ) 1 2 (t 1 n 3 + t 3 n 1 ) 1 2 (t 2 n 3 + t 3 n 2 )

				 		m 1 n 1	1 2  
	Θ[ε] m	 			 	+
							
					
						t 1 n 1	1 2
			Θ[ε] t	 
					
							(124)
	ǫ12 =	1 2	Θ([ε] (125)
	ǫ13 =	1 2	Θ([ε] (126)
	ǫ23 =	1 2	Θ([ε] (127)
	Collecting terms together yields
								
								n 2 1	n 1 n 2 n 1 n 3	
						ε = Θ[ε] n	  	n 1 n 2 n 2 2	n 2 n 3	  	+
							n 1 n 3 n 2 n 3 n 2 3

n n n

where

and

Once the value of ∆ū h is computed, the value of ∆ pf h is post-calculated using (106).

Conclusion

We have presented in this paper a method to model a poro-elastic fully saturated twophase material at the mesoscopic scale. In particular, a W/W enhancement (discontinuous strain and pressure gradient fields) implementation has been presented in the context of the Embedded Finite Element method with the corresponding discretization of the strain and pressure gradient enhanced fields, giving a contribution to the enhancements presented in Table 1. In parallel, an efficient local-global algorithmic treatment is proposed: the enhanced fields are introduced at the element level as a local enhancement of the strain and pressure gradient fields and this local structure allows the static condensation of the enhanced parameters, leading altogether to an efficient solution of these problems.

Some extensions of the method presented in this paper are also regarded, starting with the extension to partially-saturated two-phase media. This can be relevant in the context of moisture transport problem ( [START_REF] Carmeliet | Determination of the isothermal moisture transport properties of porous building materials[END_REF]) or diffusion problem ( [START_REF] Šavija | Lattice modeling of rapid chloride migration in concrete[END_REF]) in heterogeneous media.

Furthermore, consideration of crack initiation and propagation in the mortar matrix and debonding at the interface between aggregates and mortar matrix is also of our main interest. On the one hand, for the mechanical problem, this can be achieved using a strong enhancement by the means of the Strong Discontinuity Approach (SDA), introducing an Heaviside function in the displacement field. This method is well-known and massively used in the literature (see [START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate independent inelastic solids[END_REF], [START_REF] Wells | Application of embedded discontinuities for softening solids[END_REF] or [START_REF] Roubin | Multi-scale failure of heterogeneous materials: a double kinematics enhancement for EFEM method[END_REF] for instance). On the other hand regarding the literature, for the fluid problem, this enhancement can be achieved in different ways. Thus some authors consider (1) a continuous pressure field (see [START_REF] Alfaiate | On the use of strong discontinuity formulations for the modeling of preferential moisture uptake in fractured porous media[END_REF], [START_REF] Armero | An analysis of strong discontinuities in a saturated poro-plastic solid[END_REF], [START_REF] Réthoré | A two-scale approach for fluid flow in fractured porous media[END_REF]) then introducing


N Benkemoun

J.-B Colliat

A Ibrahimbegovic

Anisotropic constitutive model of plasticity capable of accounting for details of meso-structure of two-phase material

Keywords: heterogeneous materials, meso-scale model, failure surface, embedded discontinuity, biaxial loading

Anisotropic constitutive model of plasticity capable of accounting for details of meso-structure of two-phase composite material

Nathan Benkemoun, A Ibrahimbegovic, J.-B Colliat

Introduction

The non-linear analysis of structural components working under complex loading program is nowadays an indispensable ingredient of a performance based design. Given this Email address: benkemoun@lmt.ens-cachan.fr (N. Benkemoun, J.-B. Colliat and A. Ibrahimbegovic) Preprint submitted to Elsevier October 14, 2016 current trend, it has become crucial to have accurate tools that improve the predictive capabilities of macroscopic constitutive equations and furnish a precise description of nonlinear response in structural analysis. At the structural scale, phenomenological models based on macroscopic quantities such as macroscopic stresses and strains and macroscopic laws placed in the thermodynamical framework are widely used. Considering cement-based materials such as concrete, there is an extensive literature [START_REF] Reynouard | Mechanical behavior of concrete[END_REF] dealing with its phenomenological constitutive behavior modelling ( [START_REF] Lubliner | Plasticity Theory[END_REF], [START_REF] Lemaitre | Mécanique des Matériaux Solides[END_REF]) according to different loading paths, for static or dynamical cases as well as several multi-physics coupling. However, because of their macroscopic nature, these models present difficulties in describing correctly the physical mechanisms (fracture, damage and transport mechanisms) taking place at finer scales and involving macroscopic observations. Moreover in the case of a complex loading program (e.g. non-proportional loading) they require to properly choose the criterion to use in relation with the studied material and the applied loading. This can lead to substantial difficulties when devising a successful identification procedure ( [START_REF] Kucerova | Novel anisotropic continuumdiscrete damage model capable of representing localized failure of massive structures[END_REF]) which would allow to obtain the correct range of all the model parameters. So, their predictive capability is rather limited for very different loading programs with respect to the one which was used in identification. That is why the main goal of this paper is to provide a novel version of a phenomenological constitutive anisotropic model for concrete based upon information coming from finer scales and the corresponding numerical testing. Namely, we first seek to quantify correspondingly the difference in behavior in tension and in compression, as well as the fracture energy accompanying each of different modes of failure. Second, we carry out large number of numerical tests to master the inelastic mechanism's evolution from finer scales and very different loading programs. The resulting constitutive model we propose on numerical testing can be considered as the most appropriate combination of multi-surface models for concrete combining Drucker-Prager for compression stress and Rankine for tensile stress [START_REF] Nechnech | An elasto-plastic damage model for plain concrete subjected to high temperatures[END_REF] with a damage model describing the localized failure of structure ( [START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures[END_REF], [START_REF] Kucerova | Novel anisotropic continuumdiscrete damage model capable of representing localized failure of massive structures[END_REF]). The combination is not done in ad-hoc way, but constructed in accordance with a meso-scale representation of concrete distinguishing aggregates and cement paste [START_REF] Benkemoun | Modeling heterogeneous materials failure: 3D meso-scale models with embedded discontinuities[END_REF] and the result of computations corresponding to the chosen loading program. Therefore, the proposed model can account for several fine scale imprints such as compressive strength increase as a function of aggregate volume fraction or the fracture energy typical of failure mode. In this manner, we obtain a model capable of describing different failure modes, including the phenomena of localized failure which is of great interest for performance based design. This is also the main novelty of our present approach with relation to [START_REF] Benkemoun | Modeling heterogeneous materials failure: 3D meso-scale models with embedded discontinuities[END_REF] and [START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures[END_REF], we now replace a meso-scale model by an equivalent macro-scale model where failure criteria parameters depend directly of meso-scale computations.

The outline of this paper is as follows. In Section 2, we give a brief description of the meso-scale model of concrete first developed in [START_REF] Benkemoun | Modeling heterogeneous materials failure: 3D meso-scale models with embedded discontinuities[END_REF] which is adapted to the present goal.

In Section 3, we turn to providing a detailed description of a failure surface obtained with meso-scale model computations.

Meso-scale model for failure analysis of two-phase quasi-brittle material

In this section, we give a brief description of the meso-model of a typical two-phase composite material, such as concrete. For the readers interested in more details, a complete description of the model, its numerical implementation and a number of illustrative examples of the model predictive capabilities can be found in [START_REF] Benkemoun | Modeling heterogeneous materials failure: 3D meso-scale models with embedded discontinuities[END_REF].

Meso-model features

The numerical tool in [START_REF] Benkemoun | Modeling heterogeneous materials failure: 3D meso-scale models with embedded discontinuities[END_REF] is based upon a two-phase (aggregates melt into a mortar matrix) quasi-brittle finite element model capable of representing the behavior of concretelike materials under complex loading paths. In order to take into account the influence of the shape, the size, the distribution and the mechanical properties of aggregates on the mechanical behavior of concrete, the meso-scale is chosen to be the scale of computation.

This scale has been utilized by others researchers to account for heterogeneities in materials such as concrete [START_REF] Wriggers | Mesoscale models for concrete: Homogenisation and damage behaviour[END_REF] and soils [START_REF] Borja | Critical state plasticity, part VI: Meso-scale finite element simulation of strain localization in discrete granular materials[END_REF]. The meso-scale we work with to model two-phases quasibrittle material is based upon a 3D lattice finite element model [START_REF] Schlangen | Simple lattice model for numerical simulation of fracture of concrete materials and structures[END_REF], [START_REF] Schlangen | Fracture simulations of concrete using lattice models: computational aspects[END_REF], [START_REF] Yip | Automated modeling of three-dimensional structural components using irregular lattices[END_REF] and [START_REF] Lachihab | Aggregate composites: a contact based modeling[END_REF] whose truss elements kinematics is enhanced by two discontinuities.

The first discontinuity is based upon a weak discontinuity (continuous displacement field and discontinuous strain field) [START_REF] Hautefeuille | Failure model of heterogeneous structures using structured meshes and accounting for probability aspects[END_REF]. It is introduced because of the retained meshing process which relies on non-conforming mesh [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] where some truss elements are cut into two parts, each having different elastic properties. Considering a two-phase material example in 2D (Fig. 1(a)) with the mortar matrix in blue, one aggregate in green and the interface in red, the finite element discretization gives three sets of truss elements (Fig. 1(b)): those entirely inside the matrix (in blue) with no weak discontinuity activated, those entirely in the aggregate (in green) with no weak discontinuity activated and those split by a physical interface (in bold red) for which the weak discontinuity is activated. Non-conforming meshes provide the advantage to have a meshing process independent from the microstructure (positions of the aggregates and shapes), as well as the fixed size of the mesh in probability studies [START_REF] Hautefeuille | Failure model of heterogeneous structures using structured meshes and accounting for probability aspects[END_REF] of meso-structure for this kind of material. In the present work, this approach is applied to three dimensional domains. Fig. 2 represents a 0.1 m × 0.1 m × 0.1 m cube with 30 % of aggregates. We find again in blue the mortar matrix, in green the aggregates and in red the interfaces.

The second discontinuity relies upon a strong discontinuity (discontinuous displacement field and unbounded strain field) [START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate independent inelastic solids[END_REF], [START_REF] Wells | Three-dimensional embedded discontinuity model for brittle fracture[END_REF], [START_REF] Oliver | From continuum mechanics to fracture mechanics: the strong discontinuity approach[END_REF]. It is introduced in order to represent micro-cracks that may occur in any of different phases (aggregates or mortar matrix for two-phase materials) and to capture the interface failure (debonding). Moreover, the key point pertains to strong discontinuities capability to model softening behavior without any mesh dependency [START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate independent inelastic solids[END_REF], [START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations[END_REF] which is the major issue dealing with failure of quasi-brittle materials.

The weak discontinuity is present only for the truss elements split into two parts, each having a different Young modulus. The strong discontinuity is introduced by means of a yield function Φ which is triggered only in traction. Thus two constitutive models appear for a truss element : a continuum one (outside the discontinuity : Fig. 3(a)), which is elastic, and a discrete one (at the discontinuity : Fig. 3(b)) which is quasi-brittle. This approach is called the "Discrete Strong Discontinuity Approach" and can be found in [START_REF] Dias-Da-Costa | A discrete strong discontinuity approach[END_REF]. We denote with t Γ the traction vector at the discontinuity and with [|u|] the crack opening. For the discrete model, the softening law is introduced in terms of the internal variable q by considering the exponential form,

q = k([|u|]); k([|u|]) = σ u 1 -exp -[|u|] σ u G f (1)
The latter appears in the yield function which can be written as

Φ = t Γ -(σ u -q) (2)
In summary, there are altogether eight model parameters:

1. the Young modulus E 1 for the mortar matrix and E 2 for aggregates for the continuum model, 2. the ultimate tensile strength before softening, σ u i and the fracture energy, G f i (i = 1, 2, 3 for respectively the mortar matrix, aggregates and interfaces) for the discrete model. We note that G f i as the area under the curve t

Γ -[|u|].
The mathematical framework for the introduction of these discontinuities in a finite element problem is the Hu-Washizu [START_REF] Washizu | Variational Methods in Elasticity and Plasticity[END_REF] three fields variational formulation discretized by using the Incompatible Modes method [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF], [START_REF] Ibrahimbegovic | A modified method of incompatible modes[END_REF]. Among the different possibilities [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], [START_REF] Strouboulis | The design and analysis of the Generalized Finite Element Method[END_REF] to discretized the deformation field, the "Enhanced Finite Element Method" (E-FEM) [START_REF] Ibrahimbegovic | Embedded discontinuity finite element method for modeling of localized failure in heterogeneous materials with structured mesh: an alternative to extended finite element method[END_REF] has been chosen. Practically, this means that the deformation field is enhanced by two functions : G This kind of discretization leads to the following system of equations to be solved : where A denotes the standard assembly operator ( [START_REF] Ibrahimbegovic | Non Linear Solid Mechanics : Theoretical Formulations and Finite Element Solution Methods[END_REF], [START_REF] Bathe | Finite Element Procedures[END_REF]).

               A n elm e=1 f int e (d, [|u|], [|ǫ|]) -f ext e = 0 -σ 1 (d, [|u|], [|ǫ|]) + σ 2 (d, [|u|], [|ǫ|]) = 0 ∀e ∈ [1, n elm ] θσ 1 (d, [|u|], [|ǫ|]) + (1 -θ)σ 2 (d, [|u|], [|ǫ|]) -t Γ = 0 ∀e ∈ [1, n elm ] t Γ -(σ u -q) = 0 ∀e ∈ [1, n elm ] (3) phase 1 phase 2 ℓ θℓ -1 θℓ 1 (1-θ)ℓ G ①/② 1
The first equation is the classical global equilibrium obtained for a finite element problem.

The second and the third ones are respectively the local equilibrium concerning the weak discontinuity enhancement and the strong discontinuity enhancement. Both are written in function of σ 1 and σ 2 which are the stresses in each subdomain of the element and t Γ , the traction vector at the discontinuity. We recall that the scalar value θ parameterizes the discontinuity spatial position. In (3), [|u|] and [|ǫ|] are the interpolation parameters corresponding to both strong and weak discontinuities. Note that [|u|] is nothing but the crack opening. These two parameters are defined independently for each element so that they are local quantities only. It is important to note that only the second equation of ( 3) is a linear equation. Moreover, the strong discontinuity is introduced only upon reaching a fracture criterion, so that the third and fourth equations in system (3) are not always present. Computational procedure used to obtain the solution of ( 3) is explained in details in [START_REF] Benkemoun | Modeling heterogeneous materials failure: 3D meso-scale models with embedded discontinuities[END_REF].

Illustrative validation examples

In this subsection, we present several illustrative examples and predictive features of the meso-scale model.

Exact solution of a 1D problem with weak discontinuity

This first validation concerns the ability of the function

G ①/② 1
to represent a finite jump in the deformation field (weak discontinuity) and the accuracy of the incremental value

∆[|ǫ|] (k+1) computation. Consider a bi-material bar (Ω = Ω 1 Ω 2 ) of length l = 2 with material modulus E 1 in Ω 1 and with E 2 in Ω 2 ,
where Ω 1 = (-1, θl e ) and Ω 2 = (θl e , 1). The interface Γ is located at x = θl e where l e is the length of one finite element. For the pure Dirichlet problem, with u = 0 at x = -1, u = 1 at x = 1, and no body forces, the exact displacement solution is given in [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite element method[END_REF] by :

u(x) =    (1 + x)α , -1 ≤ x ≤ θl e , 1 + E 1 E 2 (x -1)α , θl e ≤ x ≤ 1, (4)
where

α = E 2 E 2 (1 + θl e ) -E 1 (θl e -1)
.

In this example, E 1 = 5 MPa and E 2 = 10 MPa. Numerical results are computed for a bar discretized with a structured mesh containing 20 finite elements with three different locations of the interface : θl e = 0.0025, 0.05 et 0.075. In Fig. 6, the finite element mesh and the interface are shown for θl e = 0.05. The numerical results for the relative error in as :

Ω 1 Ω 2
e = Ω (ǫ -ǫ h )E(ǫ -ǫ h )dΩ (6)
where ǫ and ǫ h are respectively the deformation obtained thanks to the exact solution [START_REF] Kucerova | Novel anisotropic continuumdiscrete damage model capable of representing localized failure of massive structures[END_REF] and the deformation obtained with the finite element computation. E is the Young modulus depending on the position x. By analyzing the numerical results (Tab. 1, Fig. 7 et Fig. 8), one can note that the function

G ①/② 1
is the appropriate enrichment function to capture the jump in the Young modulus values at the interface θl e . This was reasonable to expect because we sought to preserve the equivalence between the structured mesh with E-FEM space and the exact mesh with standard FEM space in which the interface coincides with a node. Namely, our results are capable of reproducing the piecewise linear nature of the exact solution (4) thanks to the chosen enrichment function which is piecewise linear.

Classical homogenization : comparison with Hashin-Shtrikman bounds for 3D case

In this example, we turn now to a weak discontinuity validation in three dimensions by exploring the Hashin-Shtrikman bounds.

At meso-scale for real concrete, we assume that the aggregates embedded in the mortar matrix exhibit an isotropic behavior because of their sizes and the material processing technique. In order to validate the meshing technique with non-conforming mesh used to represent two-phase material with aggregates embedded in a mortar matrix at meso scale, we rely upon the variational principle of Hashin and Shtrikman [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF]. This variational principle provides bounds for the effective properties of a composite that are in general more strict than the classical upper bound of Voigt and lower bound of Reuss [START_REF] Huet | An Integrated Approach of Concrete Micromechanics[END_REF]. These bounds are sensitive to specimen size and are strictly valid only when the body is assumed to be infinite, the microstructure isotropic and the effective responses are isotropic. The RVE used for the numerical analysis is a 0.1 m × 0.1 m × 0.1 m cube made of spherical aggregates melt into a mortar matrix. The aggregates size is ranging from 3 to 19 mm of diameter. Displacement boundary conditions are applied at the boundaries of this cube. On the face X = 0.1 m, a displacement of 10 -5 m is imposed in the X direction and faces X = 0 m, Y = 0 m and Z = 0 m are restrained respectively in the X, Y and Z directions. Tab. 2.2.2 gives the material data of the mortar matrix and aggregates and the effective Young modulus computed from the numerical simulations. Fig. 2.2.2 represents the effective Young modulus for the finite element model and its 3D Hashin-Shtrikman bounds for a volume fraction of aggregates ranging from 0 to 45 %. One can observe that the finite element model response is in-between Hashin-Shtrikman bounds, in spite of the fact that the sample is not of infinitely large size. This confirms its ability to correctly represent isotropic two-phases material by the chosen non-conforming mesh.

Mesh objectivity of computed response in localized failure

In this subsection, a strong discontinuity approach validation is proposed by exploring the mesh objectivity with respect to any choice of finite element discretization.

In order to illustrate the mesh objectivity in relation with the finite element discretiza- The load-imposed displacement response computed for two meshes is shown in Fig. 10. It can readily be seen that the global response is independent from the chosen mesh discretization. Indeed the dissipated energy during the loading program is the same whatever the number (and the size) of chosen elements. One can observe in Fig. 11 the computed displacement field at the end of the computation which is nearly identical for the two meshes.

Moreover, Fig. 12 shows values of displacement discontinuity or crack opening in the broken elements at the end of the computation which is again identical for the two meshes. In order to illustrate the capability of the model to represent failure mechanisms for concrete-like materials, we present two numerical simulations. The first one is a tension test and the second is the standard compression test. The specimen is 0.1 m × 0.1 m × 0.1 m cube with 24, 35 and 45 % of aggregates. Fig. 13 shows the response obtained both in tension and compression. One can observe a significant difference of behavior in tension and compression, which is one of the major characteristics related to quasi-brittle materials like concrete. Both the peak stresses as well as the amount of dissipated energy are much greater in compression than in tension. Moreover, it is important to observe the key role of aggregates in compression : a higher percent of aggregates tend to rise up the ultimate strength whereas in tension the ultimate strength is quasi-independent of this percentage increase (Fig. 13).

Failure surfaces

In this Section, we show how to use this meso-scale model, in order to provide an original approach to obtain a failure surface taking into account heterogeneities and different process of cracking.

Literature review on failure surfaces

The failure surfaces proposed in the literature for a concrete-like material are generally based upon phenomenological models ( [START_REF] Nechnech | An elasto-plastic damage model for plain concrete subjected to high temperatures[END_REF], [START_REF] Pituba | An anisotropic damage model for the concrete[END_REF], [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la ruine du béton de structure[END_REF], [START_REF] Nguyen | A coupled damage-plasticity model for concrete based on thermodynamic principles: Part i: model formulation and parameter identification[END_REF] and [START_REF] Reynouard | Mechanical behavior of concrete[END_REF]). The models of this kind rely upon a rigorous thermodynamics framework [START_REF] Ibrahimbegovic | Non Linear Solid Mechanics : Theoretical Formulations and Finite Element Solution Methods[END_REF], which is very well adapted to structural computation. Nevertheless, they require to properly identify and determine model parameters in relation with the studied material and the applied loading program. Thus the choice of model parameters must be made with a consistent calibration procedure (e.g. see [START_REF] Kucerova | Novel anisotropic continuumdiscrete damage model capable of representing localized failure of massive structures[END_REF]). In the case of a failure surface, this should include the calibration of the initial yield and damage surfaces that are in accordance with experiments, along with other parameters governing the evolutions of these surfaces. Moreover, because of their macroscopic nature, these models present difficulties in taking into account explicitly heterogeneities (e.g. aggregate shape, distribution, sizes and behavior) as well as the corresponding information at finer scales, which plays a fundamental role ( [START_REF] Yaman | Active and non-active porosity in concrete. part I: experimental evidence[END_REF]) in most of the physical phenomena observed at the macroscopic scale. In addition, in order to describe correctly physical mechanisms (as failure mechanisms, for example) taking place at finer scales and responsible of macroscopic observations, phenomenological models require again the identification of several parameters ( [START_REF] Nguyen | A coupled damage-plasticity model for concrete based on ther-modynamic principles: Part ii: non-local regularization and numerical implementation[END_REF]) and remain of fairly limited predictive capabilities outside of the given loading program.

That is why in this section, a new approach is presented: the coupling between the fine scale, here chosen as the meso-scale, where the quasi-brittle behavior is introduced, with the macroscopic scale, where the failure surface is obtained by the numerical computation of ultimate strength, is done naturally. Namely, every points of the surface take into account failure mechanisms taking place at the fine scale and the material heterogeneity. In this sense, a novel version of a phenomenological constitutive anisotropic model for concrete based upon information coming from finer scales is obtained. The proposed model can approximately be considered as the most appropriate combination of multi-surface models for concrete combining Drucker-Prager for compression stress and Rankine for tensile stress [START_REF] Nechnech | An elasto-plastic damage model for plain concrete subjected to high temperatures[END_REF] with a damage model describing the localized failure of structure ( [START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures[END_REF], [START_REF] Kucerova | Novel anisotropic continuumdiscrete damage model capable of representing localized failure of massive structures[END_REF]). However, contrary to this classical combination, the proposed model can also provide the compressive strength related to the volume fraction, as well as the fracture energy corresponding to each particular failure mode.

Numerical loading program

In order to examine the ability of the Finite Element model to provide biaxial failure behavior characteristic of a heterogeneous material specimen, three concrete specimens with a plate form of 100 mm × 100 mm × 10 mm (Fig. 14) and a different volume fraction of aggregates (Table 3) were subjected to biaxial imposed displacements, producing the following stress combinations : compression-compression, tension-tension, compression-tension and tension-compression. Their uniaxial compressive strength values σ c i (i = 1, 2, 3) were taken from a uniaxial compressive test performed on each concrete specimens. In order to numerically cover the four regions (compression-compression, tension-tension, compression-tension and tension-compression), a circular loading program based upon the trigonometrical circle cos(α) 2 + sin(α) 2 = 1 was performed. This leads to couples of imposed displacements (u 1 = cos(α), u 2 = sin(α)) parameterized by the angle α. Displacements u 1 and u 2 are respectively applied to faces X = 100 mm and Y = 100 mm. Faces X = 0 mm, Y = 0 mm and Z = 0 mm are respectively blocked in X, Y and Z directions (Fig. 15).

Specimen i Volume fraction f i (%) Uniaxial compressive Size aggregates range strength σ c i (MP a) (mm)
Tension is chosen as positive. Stress in the principal direction 1 (X direction) is called σ 1 , the one in the principal direction 2 (Y direction) σ 2 . In addition, the ultimate strength in the principal direction 1 is written σ1 and the one in the direction 2, σ2 . Finally, in the region of tension-compression, the value of the ultimate strength decreases in compression when the one in tension increases. In general, the failure pattern is perpendicular to the principal tension direction. Most of these tendencies described in this part are consistent with experimental results of [START_REF] Kupfer | Behavior of concrete under biaxial stress[END_REF], [START_REF] Carpinteri | Fracture Mechanics of Concrete[END_REF] and [START_REF] De | Étude expériementale du comportement mécanique d'un Béton Fibré à Ultra Hautes Performances (BFUHP) en traction biaxiale[END_REF]. For the regions of tension-tension, compression-tension and tension-compression, two phases of computed response can be observed. In the zone of compression-compression, a linear behavior is observed between σ1 /σ c 2 and σ2 /σ c 2 : cracks are prevented from opening (the loading in the Y direction tends to prevent the crack from opening in the X direction and the loading in the X direction tends to prevent the crack from opening in the Y direction), their influences is therefore less important than in the others cases of loading.

(mm)) G f (f 1 ) (J) G f (f 2 ) (J) G f (f 3 ) (J)
To conclude Section 3, one can say that Fig. 16 confirms the interest of concrete plasticity criterion, of multi-surface kind, with the compression described by Drucker-Prager kind of model [START_REF] Dolarevic | A modified three-surface elasto-plastic cap model and its numerical implementation[END_REF] and tension regime described by Rankine model [START_REF] Pearce | On multi-surfaces plasticity and rankine model[END_REF].

Conclusion

In this work, we proposed a novel constitutive model of plasticity suitable for describing the failure phenomena in two-phase composite materials, such as concrete with aggregates melt into a mortar matrix. Contrary to usual phenomenological models, the main features of the proposed model are not obtained from performing real experiments, but rather from numerical experiments carried out with a meso-scale Finite Element model. The mesomodel relies on a non-conforming mesh of spatial truss elements, representing aggregate, cement paste and interface so that some elements must be split into two parts. The latter requires introducing a weak discontinuity into the elements split by the phase interface, which provides an effective method to take into account inclusions without remeshing. Moreover, a strong discontinuity is introduced within each element at a particular instant of localized failure providing the capability to model crack opening. These key points provide an efficient tool to accomplish "virtual testing"; where changes in size, distribution and volume fraction of aggregates can easily be introduced.

Three types of concrete with different volume fractions (15 %, 35 % and 40 %) were tested.

For these three specimens, failure surfaces and failure modes according to different loading combinations were developed. One can note that in the region of compression-compression the ultimate strength rises up significantly with the percentage of inclusions whereas the three others regions are quasi independent of this percentage variation. It is important to stress the fact that the ultimate strength increase in relation to the uniaxial compression strength is dependant of the imposed displacements ratios in the region of compressioncompression, whereas the three others regions remain quasi independent of the same ratio.

These two points are often confirmed as major features of quasi-brittle materials such as concrete.

Finally, the most notable finding is that this approach can provide the sound definitions of the parameters governing the failure process, such as the fracture energy, the ultimate strength taking into account the heterogeneity and different process of cracking and that can be used to identify the parameters of phenomenological models as multi-surface models (e.g.

the combination of Rankine and Drucker-Prager criterion or the three-surface plasticity cap model [START_REF] Dolarevic | A modified three-surface elasto-plastic cap model and its numerical implementation[END_REF]). Once these phenomenological models are clearly identified, a structural analysis could be performed with the potential to provide more predictive results from the classical phenomenological models.

In the future explorations, we plan to study how to account for induced anisotropy in the process of cracking. Thus, the additional computations will be performed for the numerical loading program on pre-cracked specimen where the initial crack can be computed from drying and early age models of concrete [START_REF] Sa | Modélisation à l'échelle mesocopique du comportement hydro-mécanique des matériaux à matrice cimentaire[END_REF].

Figure 1 :

 1 Figure 1: 2D discretization with a non-conforming mesh

Figure 2 :

 2 Figure 2: 3D domain for a typical two phases material: mortar matrix blue, aggregates green, red interfaces

Figure 3 :

 3 Figure 3: Elastic-quasi-brittle behavior

①/② 1 and G 2 .

 12 The first one is a weak discontinuity, introduced as a piecewise linear function over an element capturing the jump of the Young modulus in the deformation field. The second one is split into a constant function and a Dirac function capturing the unbounded nature of the strain field in presence of a strong discontinuity. One can observe these functions on Fig.4and Fig.5thanks to a zoom on a interface truss element (red bold elements on Fig.1(b) and Fig.2).

Figure 4 :Figure 5 :

 45 Figure 4: Enhanced deformation field for a two-phase truss element

Figure 6 :

 6 Figure 6: Bi-material bar problem : mesh and interface position for θl e = 0.05

Interface location ξ

  Relative error in the energy norm 0.0025 1.3 × 10 -10 0.05 3.12 × 10 -10 0.075 3.18 × 10 -10

Figure 7 :Figure 8 :

 78 Figure 7: Exact and computed deformations fields for θl e = 0.0025

Figure 9 :Figure 10 :

 910 Figure 9: Hashin-Shtrikman, Voigt and Reuss bounds

Figure 11 :Figure 12 :

 1112 Figure 11: Displacement magnitude at the end of the computation: a) 260000 elements, b) 690000 elements

Figure 13 :

 13 Figure 13: Macroscopic response in tension/compression

Figure 14 :Figure 15 :

 1415 Figure 14: Specimen and microstructure with 40 % of aggregates

Figure 16 :

 16 Figure 16: Failure surfaces for 15 % (diamond), 35 % (pentagon) and 40 % (triangle) of aggregates

Figure 17 :

 17 Figure 17: Evolution of σ 1 /σ c2 with respect to σ 2 /σ c2 , during the loading until failure (green points), in the case of 40 % of aggregates

Table 1 :

 1 Relative error in the energy norm for the bi-material bar problem

Table 3 :

 3 Specimen data

  The first one corresponding to a linear relation between σ1 /σ c 2 and σ2 /σ c 2 : cracking is not very present (few elements are broken) so a linear relation between the two ratios, driven by a macroscopic Poisson ratio, is still Combined regions (u 1 (mm),u 2 (mm)) σ1 (f 3 )/σ c 2 σ2 (f 3 )/σ c 2

	tension-tension	(0.01,0.0)	0.263	0.087
		(0.005,0.0087)	0.261	0.290
		(0.0087,0.005)	0.291	0.256
		(0.0,0.01)	0.100	0.291
		(0.0096,0.0026)	0.299	0.192
		(0.0026,0.0096)	0.196	0.289
		(0.0099,0.0013)	0.286	0.142
		(0.0038,0.0092)	0.235	0.290
		(0.01,0.01)	0.249	0.244
	compression-compression	(-0.06,0.0)	-1.400	-0.277
		(-0.03,-0.05)	-0.999	-1.372
		(-0.05,-0.03)	-1.327	-0.857
		(0.0,-0.06)	-0.280	-1.383
		(-0.058,-0.015)	-1.376	-0.553
		(-0.06,-0.06)	-1.195	-1.137
		(-0.015,-0.058)	-0.636	-1.449
	compression-tension	(-0.005,0.0087)	-0.266	0.249
		(-0.0087,0.005)	-0.537	0.123
		(-0.0071,0.0071)	-0.391	0.197
		(-0.0061,0.0079)	-0.415	0.224
	tension-compression	(0.005,-0.0087)	0.128	-0.517
		(0.0087,-0.005)	0.247	-0.255
		(0.0071,-0.0071)	0.198	-0.376
		(0.0061,-0.0079)	0.167	-0.446

Table 5 :

 5 Ratios of ultimate strength with respect to σ c2 along X and Y directions under different couples of imposed displacements for 40 % of aggregates

	Combined regions	(u 1 (mm),u 2

Table 6 :

 6 Values of fracture energy under particular couples of imposed displacements for 15 %, 35 % and 40 % of aggregates available. As for the second phase, it corresponds to a non-linear relation between the two ratios : cracking is then important and the behavior is no longer governed by a macroscopic

	Poisson modulus.

Combined regions

(u 1 (mm),u 2 (mm)) σ1 (f 1 )/σ c 2 σ2 (f 1 )/σ c 2 σ1 (f 2 )/σ c 2 σ2 (f 2 )/σ