Laurent Charpin
email: laurent.charpin@enpc.fr

A Ehrlacher

Keywords: Micromechanics, Cracked elastic solid, Interaction Direct Derivative, finite element simulation

We present a comparative study of the ability of some micromechanics estimates to predict the overall properties of heterogeneous materials. We focus mainly on cracked materials, for which this task is difficult and many estimates fail. We study particularly the Interaction Direct Derivative estimate, proposed by Zheng and Du, which is an approximation of the generalized self-consistent scheme but has the very convenient property to be always explicit. A modified version of this estimate, called Full Range IDD by Zheng and Du, yields good results when comparing all poromechanical coefficients predicted by the estimate to finite element simulations of a 2d cracked material in plane strain, up to crack density factors of 1 for aligned cracks and 0.60 for randomly oriented cracks. The accuracy of finite element computations of the overall moduli is also commented by plotting the convergence of the average of the properties as well as the confidence intervals on these averages.

Introduction

The field of micromechanics is dedicated to the prediction of the overall properties of composites, starting from information about the volume fraction, morphology, and mechanical properties of the various phases. Numerous estimates for the mechanical properties have been proposed, and some of them have become very classical, such as the dilute scheme (directly derived from Eshelby's solution [START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF]), the self-consistent scheme [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF], the differential scheme [START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF], and the Mori-Tanaka scheme [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]. These estimate give more or less good results depending on the volume fractions of the phases and the morphology types considered. For example, Mory-Tanaka's estimate is considered to be particularly adapted for matrix/inclusion morphologies.

Part from these very classical estimates, more refined estimates exist. The general self-consistent estimate seems the most complete [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF]. It is, as the mentioned estimates, based on a simple auxiliary problem on which the localization tensors on the inclusionary phases of the real problem are estimated. It has the advantage of taking the distribution of the inclusions into account, in addition to their shape. This is done through an ellipsoid containing the matrix material, embedding the inclusion, and itself embedded into the homogenized material which is sought for. However, it has a closed-form solution in only few cases: isotropically distributed spherical inclusions (this solution was given by Christensen and Lo [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF]) and some cases of cracked medium [START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF]. For other cases, it is necessary to use approximations (such as the one proposed by Hori and Nemat-Nasser in the same article [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF]), or to simplify the problem (as is done in the two-phase model, where the matrix material also occupies the infinite medium embedding the matrix ellipsoid [START_REF] Shige | Effective elastic moduli of inhomogeneous solids by embedded cell model[END_REF]). It is also possible to evaluate the general self-consistent solution numerically, as is well known and will be illustrated in this article.

More recently, two improvements have been achieved. First, Ponte-Castañeda and Willis showed precisely the role of the distribution of the inclusions in the framework of Hashin-Strickman principle [START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF][START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF]. This work, which can be used to provide bounds or an estimate for the elasticity tensor, clarifies the underlying microstructures more or less explicitly considered when establishing the classical estimates. A major point is that the distribution functions of pairs of inclusion families need to be symmetric. In most estimates which consider the distribution of the inclusions, only one distribution is chosen for each inclusion, neglecting this point. However, when the various inclusion families have different distributions, the estimate is not explicit and difficult to use. Second, Zheng and Du proposed an estimate based on the generalized self-consistent model called the Interaction Direct Derivative estimate, which has the advantage of being fully explicit in all cases, even with multiple inclusion families, each having a different space distribution [START_REF] Zheng | An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution[END_REF].

Concerning the accuracy of these estimates, the classical estimates give good results with simple morphologies such as discs in 2d or spheres in 3d [START_REF] Miled | Effective elastic properties of porous materials: Homogenization schemes vs experimental data[END_REF][START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF][START_REF] Klusemann | Homogenization methods for multi-phase elastic composites: Comparisons and benchmarks[END_REF]. However, when moving to cracked microstructures, it is much more difficult to predict the overall properties of a Representative Elementary Volume of the material. There has been a number of studies in this direction in 2d (see [START_REF] Huang | The numerical calculation of two-dimensional effective moduli for microcracked solids[END_REF][START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF]), but it seems some improvements remain possible. We won't tackle the problem of 3d cracked structures which requires large computational tools for a rigorous treatment, but we will attempt to assess the efficiency of Interaction Direct Derivative estimate of Zheng and Du for cracked structures in 2d. We will study its efficiency by comparison to finite element computations, using monodisperse cracks, to the highest crack density factor we were able to reach using a Random Sequential Addition algorithm [START_REF] Torquato | Morphology and effective properties of disordered heterogeneous media[END_REF], which was roughly 1 in the case of aligned cracks and 0.6 in the case of randomly oriented cracks. Of course this limit could be improved by using multiple crack sizes, but we won't discuss this point in this article.

We believe that our article brings interesting information in addition to Shen and Li [START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF] where the case of parallel and aligned cracks in 2d is also dealt with. First, since we deal with cracks, we think that it is important to foresee the application to poromechanics of the homogenization procedure, by computing and displaying the overall poromechanical coefficients such as the Biot coefficient and Biot modulus [START_REF] Coussy | Poromechanics[END_REF], which we do systematically when comparing the simulations to the estimates. Second, we use the IDD estimate which is very versatile, unlike the solutions developed in some particular cases of aligned or randomly oriented cracks (see Shen and Li [START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF]), and unlike the differential scheme and the self-consistent scheme in complicated situations. Third, we take some time to comment the error made in the finite element simulations, so that the reader knows exactly to what extent our results can be trusted.

The authors used the IDD estimate for the prediction of degradation of concrete submitted to Alkali-Silica Reaction. This estimate was used to improved a linear elastic fracture mechanical model proposed in [START_REF] Charpin | A computational linear elastic fracture mechanics-based model for alkali-silica reaction[END_REF]. Therefore, more details about the comparison of the IDD estimate and finite element simulations can be found in the author's Ph.D. thesis, as well as applications to Alkali-SIlica Reaction [START_REF] Charpin | Modèle micromécanique pour l'étude de la réaction alcali-silice[END_REF].

The article is organized as follows. First, we give a brief explanation of the variational formulations as well as the boundary conditions we use in the finite element code. Second, we describe the procedure of determination of the overall properties. We then present the various microstructures that we study in this article, and explain how we describe the error made in the finite element simulations. We then move to the micromechanics by recalling the expressions and main features of the estimates we use, and finally compare them to the results obtained in the simulations, showing that the IDD estimate is well suited to compute the overall properties of cracked materials.

Plane strain variational formulation

All computations were done in FreeFem++ [START_REF]FreeFem++[END_REF]. In this software the user implements its own variational formulation, which will be shortly presented here. We decided two work in 2d, thinking that in the view of assessing the efficiency of micromechanics estimates, we could obtain first results. Of course to be able to test all possibilities in terms of orientation of inclusions, for example, 3d simulations would be preferable. We will also present shortly the boundary conditions that were used for the computations.

Our variational formulation is written starting from the virtual work principles, as done classically, in a static case [START_REF] Bonnet | Analyse des solides déformables par la méthode des éléments finis[END_REF]. We call Ω the continuum solid of local tensor of elasticity C(x), and boundary ∂Ω. The external boundary is split between the external boundary ∂Ω ext and an internal boundary ∂Ω int which separates the solid from pores embedded in the solid. The internal boundary is submitted to pressure loadings, while the external boundary can be submitted to imposed tractions T D and imposed displacements u D .

Since we are in plane strain, the displacement is assumed to be zero in direction 3 and hence, the related components of the strain are also zero:

ε 13 = ε 23 = ε 33 = 0 (1)
Let us quickly describe the boundary conditions we use.

Imposed average strain

The boundary conditions we call imposed strain are such that the imposed displacement on the external boundary is u D (x) = E.x , ∀x ∈ ∂Ω ext , where E is a second order, symmetric, plane strain tensor.

Calling C(u D ) the space of displacements admissible with u D , and assuming a pressure p exerts on the internal boundary, the variational formulation is (calling ε the strain tensor):

Find u ∈ C(u D ) such that ∀w ∈ C(0), Ω ε[u] : C : ε[w]dV + ∂Ωint pn.wdS = 0 (2)

Imposed average stress

The boundary conditions we call imposed stress are such that the tractions on the external boundary are equal to T D (x) = Σ.n , ∀x ∈ ∂Ω ext where Σ is a second order symmetric tensor. Let us remark that the average stress is only partly imposed because we are still in the plane strain framework. Therefore, only the components (Σ 11 , Σ 12 , Σ 22 ) are imposed. Calling C the space of admissible plane strain displacements, still considering the pressure p is applied in the pores, the variational formulation is:

Find u ∈ C such that ∀w ∈ C, Ω ε[u] : C : ε[w]dV + ∂Ωint pn.wdS - ∂Ω T ext (Σ.n).wdS = 0 (3)

Periodic boundary conditions

The boundary conditions we call periodic are such that the displacement on the external boundary is equal to u = E.x + u ′ , ∀x ∈ ∂Ω ext where the field u ′ is periodic and E again is a plane strain second order symmetric tensor. Let us decompose the strain as ε[u] = E + ε[u ′ ] and call C per the space of periodic place strain displacement fields. The variational formulation is:

Find u ′ ∈ C per such that ∀w ∈ C per , Ω ε[u ′ ] : C : ε[w]dV + Ω E : C : ε[w]dV + ∂Ωint pn.wdS = 0 (4)
These various types of boundary conditions in place strain allow us testing our microstructure to determine their average properties. Let us see how we determine these properties from the numerical results.

3 Determination of the average properties First, we need to be able to compute strain and stress averages. For a quantity f and a domain V , the volume average of g over V is < g > V = V gdV . Let us remind these quantities can be computed as surface integral which is sometimes convenient. For example for the average strain:

< ε > Ω = 1 Ω ∂Ωext u ⊗ s ndS (5)
Where ⊗ s is the symmetrized tensorial product. Let us also recall how the pressure appears in the average stress, assuming the volume fraction of pores is f = Ω p /Ω:

< σ > Ω = -f p½ + 1 Ω Ω\Ω p σdV (6)

Overall constitutive law

Knowing the local elastic behaviour of the solid, and for various imposed macroscopic strains and pressures in the pores, we want to be able to write the macroscopic constitutive law of our porous medium. Following the presentation of Dormieux et al. [START_REF] Dormieux | Microporomechanics[END_REF], assuming in this work that there is only one pore family (which does not mean that they are all identical, but that we study their volume changes and the effect of an applied pressure in them as a whole), we can write the macroscopic constitutive law linking the macroscopic stress Σ and the porosity variation of the porous phase (φ -f ) to the average strain E and the pressure in the pores p:

Σ = C hom : E -pB φ -f = B : E + pM (7)
Where C hom is the homogenised stiffness tensor, B is the Biot coefficient, while we call M the Biot compliance, inverse of the usual Biot modulus N as presented in Dormieux et al. [START_REF] Dormieux | Microporomechanics[END_REF]. These homogenized properties can all be computed using the solution under imposed strain alone, without solving the problem with a pressure in the pores (this is specific to the case of equation 7 where all cavities are considered together). The strain localisation tensor is a fourth order tensor with the minor symmetry A (such that when a macroscopic strain E is imposed, the local strain writes, by virtue of the linearity of the problem, A (x) : E). Assuming that the solid matrix is homogeneous of tensor elasticity C 0 and compliance S 0 , the microporomechanics theory then allows writing the homogenized poroelastic properties as:

   C hom = (1 -f )C 0 :< A > Ω\Ω p B = f ½ :< A > Ω p M = ½ : S 0 : B -f ½ (8)
3.2 Precise determination of the overall properties from strain and stress averages

In this section we provide the expressions we used for the determination of the properties, in the general case where the overall medium is orthotropic. Some of these expressions simplify if the medium is transversely isotropic but we won't detail this case. Let us mention that we use the Voigt notation, with the convention that all symmetric second order tensors (particularity the stress, strain, and Biot coefficient) are written:

A v = A 11 A 22 A 33 √ 2A 32 √ 2A 31 √ 2A 21 t (9)
The fourth order are written accordingly:

D v =         D 1111 D 1122 D 1133 √ 2D 1123 √ 2D 1113 √ 2D 1112 D 2211 D 2222 D 2233 √ 2D 2223 √ 2D 2213 √ 2D 2212 D 3311 D 3322 D 3333 √ 2D 3323 √ 2D 3313 √ 2D 3312 √ 2D 2311 √ 2D 2322 √ 2D 2333 2D 2323 2D 2313 2D 2312 √ 2D 1311 √ 2D 1322 √ 2D 1333 2D 1323 2D 1313 2D 1312 √ 2D 1211 √ 2D 1222 √ 2D 1233 2D 1223 2D 1213 2D 1212         (10)
Second, let us remind that we work in plane strain to reproduce a 3d material which would be invariant and infinitely extended in direction 3. The matrix is of Young's modulus and Poisson's ratio (E, ν). This analogy to the infinite medium case is useful when coefficients relative to direction 3 are needed.

It can be shown that regardless of the shape of the voids, at volume fraction f , the Young's modulus and Poisson's ratio in direction 3 write ( [START_REF] Nemat-Nasser | Micromechanics: overall properties of heterogeneous materials[END_REF], p.85):

E 3 = (1 -f )E (11)
ν 32 = ν 31 = ν (12)
While the Biot coefficient in direction 3 writes:

B 33 = f [1 + ν(A v 11 + A v 12 + A v 21 + A v 22 )] (13)
Let us then show, as an example, how we determine the poromechanical coefficients in the case where the macroscopic strain is prescribed (imposed strain or periodic boundary conditions). Three simulations are performed (we don't use equation 8, but determine the coefficients directly from equation 7). During the first one, the following macroscopic strain is imposed:

(E I ) v = 1 0 0 0 0 0 t (14)
And average stresses Σ I 11 , Σ I 22 et Σ I 12 are measured. Second, the macroscopic strain

(E II ) v = 0 1 0 0 0 1 t ( 15
)
is imposed, while average stresses Σ II 11 , Σ II 22 et Σ II 12 are measured. We can then write the following expressions for the elastic properties:

E 1 = Σ I 11 Σ I 22 2 - Σ II 11 Σ II 22 2 1 (Σ I 22 ) 2 Σ I 11 + ν 2 E3 Σ I 11 2 -Σ I 22 2 + 1 (Σ II 22 ) 2 Σ II 22 + ν 2 E3 Σ II 22 2 -Σ II 11 2 (16)
E 2 = Σ II 22 Σ II 11 2 - Σ I 22 Σ I 11 2 1 (Σ I 11 ) 2 Σ I 11 + ν 2 E3 Σ I 11 2 -Σ I 22 2 + 1 (Σ II 11 ) 2 Σ II 22 + ν 2 E3 Σ II 22 2 -Σ II 11 2 (17)
ν 12 = E 1 Σ II 22 Σ II 11 E 1 - ν 2 E 3 Σ II 11 + Σ II 22 ; µ 12 = Σ II 12 2 (18)
Concerning the Biot compliance and coefficient, they are determined imposing zero average strain and a unit pressure in the pores. Using this solution, we then have by definition (see equation 7):

B 11 = -< σ 11 > Ω B 22 = -< σ 22 > Ω (19)
M = φ -f = f < tr ε > Ω p (20)
The case of imposed stresses is dealt with similarly, but the expressions are not detailed here.

Studied microstructures

Let us now describe the type of microsctructures whose overall properties we will study. Let us first talk about the generation of the random microstructures, and then describe the periodic cells we used.

Random microstuctures built by Random Sequential Addition

This process allows creating matrix/inclusion morphologies easily [START_REF] Torquato | Morphology and effective properties of disordered heterogeneous media[END_REF][START_REF] Huang | The numerical calculation of two-dimensional effective moduli for microcracked solids[END_REF][START_REF] Franciosi | Anisotropy features of phase and particle spatial pair distributions in various matrix/inclusions structures[END_REF][START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF]]. We show an example and monodisperse spheres on Fig. 1. Let us briefly sum up the procedure. First, a domain size, inclusion

Fig. 1 One realisation of a microstructure with isotropically distributed circular pores size and volume fraction f are chosen. Then the coordinates of the center of the particle are randomly generated (and, for example, its orientation if necessary). If the particle has no intersection with the external boundary of the domain nor any previously added particle, it is kept. This operation is then repeated until the desired volume fraction is reached, which is more accurate if there is a large number of inclusions.

The first difficulty of the method is that for particle shapes other than aligned ellipsoids, testing the intersection of two particles can be difficult and can require the discretisation of the boundary of the particle. Second, even if for discs in 2d the maximal density is roughly 0.9, this method, for monodisperse inclusion sizes, creates empty spaces between particles that are too small to be filled by a new particle. The density is therefore limited, for spheres, to roughly 0.5. Finally, for meshing reasons, the particles cannot touch each other. The small exclusion distance reduced even more the density we can reach.

Periodic microstructures

Periodic simulations are convenient and easy to implement. They are based on the choice of a periodic pattern. In our case, for spheres, we followed the choice of Zheng and Du [START_REF] Zheng | An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution[END_REF], placing the centers of the spheres on a hexagonal pattern. This allows attaining, in 2d, optimal densities for monodisperse spheres. In the case of aligned ellipsoids, there are two interesting possible choices. The first possibility is to keep the ellipsoids centered on the same patterns as the spheres, representing an isotropic space distribution of the ellipsoids, but then the volume fraction is geometrically very limited. Another idea is to flatten the hexagonal pattern on which the ellipsoids are placed by the same factor as the aspect ratio of the ellipsoids. An example of obtained pattern in this case is shown on Fig. 2. In this case we can reach the same high density as in the case of spheres.

Remarks concerning the accuracy of the results

Concerning the size of the elements of the meshes used for the computation, an automatic mesh refinement procedure based on the geometrical error in the representation of the microstructure was used, as well as a limitation of the elongation of the elements. The results are converged with respect to the mesh size.

Another issue is the fact that since our inclusions don't intersect the external boundary of the domain, the inclusion density is lower close to the boundary. However, since the convergence of the poromechanical properties is verified with respect to the ratio between the inclusion size and domain size, the final values which are used for the comparisons with the estimates are not affected by this phenomenon.

Determination of the REV size in a simple case, with Kanit's approach

To obtain representative values of the homogenized properties thanks to the numerical simulations, we must use domains which are large enough compared to the inclusion sizes. The convergence of these properties with respect to the domain size can be checked for the different types of boundary conditions.

Let us first recall that the obtained tensors of elasticity using our simple boundary conditions always verify the following relation for any second order symmetric tensor E:

E : C hom Σ : E ≤ E : C hom general : E ≤ E : C hom E : E (21)
Where C hom Σ is determined using the imposed stress boundary conditions, C hom E is determined using the imposed strain boundary conditions, while C hom general refers to any boundary conditions. This property of the different stiffness tensor that can be determined on the same microstructure is well known [START_REF] Nemat-Nasser | Micromechanics: overall properties of heterogeneous materials[END_REF][START_REF] Hori | On two micromechanics theories for determining micro-macro relations in heterogeneous solids[END_REF].

As a consequence, the difference between the moduli obtained with imposed stresses or strains also give an insight about the error which is made.

Description of the error on the estimation of the mechanical properties

Considering a domain of calculus of size L and a characteristic size of the inclusions R, we study the convergence of the mechanical properties, for the three types of boundary conditions, when R L → 0. There has been many studies of this kind, see for example Gusev [START_REF] Gusev | Representative volume element size for elastic composites: A numerical study[END_REF]. Two sources of error appear: first, a bias due to the the lack of representativity of small domains, which cannot be suppressed even with an infinite number of realizations, and which is more or less pronounced depending on the boundary conditions used. Second, the insufficient number of realizations, which lead to imperfect estimation of the average. Following [START_REF] Anantharaman | Introduction to numerical stochastic homogenization and the related computational challenges: some recent developments[END_REF] and [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF], we compute the average and the variance of the properties we are interested (say A i ) for N realizations, by using the classical formula for the average µ

N (A) = 1 N N i=1 A i and for the variance σ N (A) = 1 N -1 N i=1 (A i -µ N (A))
2 . Hence we can plot, for each poroelastic quantity, the empirical average and the confidence interval at 95 % (µ N (A) ± 1,96 σ N (A)/N ) as a function of the ratio R L .

Example: circular inclusions

We study the three types of boundary conditions, for a fixed volume fraction of circular cavities f , hence inducing transversely isotropic overall poromechanical properties. We plot the in-plane shear modulus µ 12 , a plane strain two-dimension bulk modulus we define as:

k 2ddp = 1 2 Σ 11 + Σ 22 E 11 + E 22 (22)
from the macroscopic strains and stresses, as well as the Biot compliance M and Biot coefficient B. Looking at Fig. 3, we observe that the error bars, representing the uncertainty to which we know the empirical average, get smaller with increasing R/L, at fixed number of realizations N . The imposed strain and imposed stress solutions are biased, which means that even at large number of realizations, we don't get the converged value of the poromechanical properties. Finally, the periodic boundary conditions show almost no bias, which is encouraging because it means it is possible to use much smaller ratios of the domain size to heterogeneity size to obtain good estimates of the poromechanical properties than with the other boundary conditions. Therefore, from now on, only the periodic boundary conditions will be used.

Let us recall the expressions of the micromechanics estimates we will compare to these finite element results.

Micromechanics estimates

For each estimate, we give the expression of the strain localization tensor on the inclusionary phase, as well as the homogenized stiffness tensor, except for the general self-consistent estimate that is not analytical. The classical micromechanics results from Dormieux et al. [START_REF] Dormieux | Microporomechanics[END_REF] reminded in equation 8 allow writing the Biot modulus and coefficient easily once the localization tensor is known.

Dilute scheme (dil) The dilute scheme is not shown in the comparisons but is used to build other estimates. Considering N families of ellipsoidal heterogeneities, embedded in a matrix of stiffness tensor C 0 , each caracterized by their stiffness tensor C i , Hill tensor P 0 i and volume fraction f i , the dilute localization tensor of family i writes:

A dil i = I + P 0 i : C i -C 0 -1 (23)
It is the direct application of Eshelby's result [START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF] and it is used in the construction of all analytical micromechanics estimates. The stiffness tensor is then simply obtained by the following general expression for the stiffness tensor:

C dil = C 0 + N i=1 f i (C i -C 0 ) : A dil i ( 24
)
Two-phase estimate (2ph) The two phase estimate, when there is only one heterogeneity family, is built by embedding the heterogeneity ω of stiffness tensor C i in an infinite elastic medium of same stiffness tensor as the matrix of the composite C 0 [START_REF] Shige | Effective elastic moduli of inhomogeneous solids by embedded cell model[END_REF]. Then an ellipsoid D embedding the inclusion ω is defined to represent the spatial distribution of the inclusion rather than its shape, and so that the volume fraction of the inclusions in the real composite f is equal to ω/D. Both ellipsoids are characterized by their Hill tensor. This defines a new localization tensor, linking the average strain on ω to the average strain on D when a strain is applied at infinity:

A 2ph ω = A dil ω : I -f P 0 D : C i -C 0 : A dil ω -1 (25)
Mori-Tanaka estimate (MT) We then build the famous Mori-Tanaka estimate from the two-phase estimate, by assuming that the atmosphere D is similar to the inclusion ω, which means we use the same Hill tensor, and use equation 25 [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF][START_REF] Weng | The theoretical connections between Mori-Tanaka's theory and the Hashin-Shtrikman-Walpole bounds[END_REF]. The moduli of cracked materials were originally computed by Zhao et al. [START_REF] Zhao | Elastic moduli for a class of porous materials[END_REF].

Self-consistent estimate (SC)

The self-consistent estimate assumes that the inclusion is embedded in the homogenized material rather than the matrix material [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF][START_REF] Budiansky | Elastic moduli of a cracked solid[END_REF]. Therefore, the self-consistent localization tensor is similar to the dilute localization tensor except for this difference:

A SC i = I + P SC i : C i -C SC -1 (26)
Where C SC is unknown, but can be expressed as a function of the localization tensor using equation 24:

C SC = C 0 + f (C i -C 0 ) : A SC i ( 27
)
Which yields an implicit expression for C SC , which can be difficult to solve when the homogenized medium is not isotropic due to much more complicated expressions of the Hill tensor appearing in the equation. We will only use this estimate in a simple case where we have a closed form solution to this equation. If there is one family of spherical holes in an homogeneous material characterized by a shear modulus and Poisson's ratio (µ 0 , ν 0 ), then the evolution of the in-plane shear modulus and plane strain two-dimensional bulk modulus as defined in equation 22 write (see Nemat-Nasser and Hori [START_REF] Nemat-Nasser | Micromechanics: overall properties of heterogeneous materials[END_REF], p. 103):

µ AC 12 µ 0 = 1 -3f 1 + f (1 -4ν 0 ) ; k AC 2ddp k 2ddp 0 = 2 (3 -4ν 0 )(1 -ν 0 )(1 -3f ) 3 -4ν 0 1 -f (1 -4ν 0 ) (28)
Interaction Direct Derivative estimate (IDD) This estimate is built on the ideas of the generalized self-consistent problem [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF], which is recognized to be very good for matrix/inclusion morphologies, but has no closed-form solution except in some simple cases. The auxiliary mechanical problem which is considered to determine the localization tensor for inclusion family i is the the following: inclusion i is embedded in its atmosphere of ellipsoidal shape D i so as to respect the volume fraction f i = ω i /D i and of stiffness tensor that of the matrix. The atmosphere is itself embedded into an infinite medium which has the mechanical properties we are looking for, as in the self-consistent model. While in Hori and Nemat-Nasser [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF], the meaning of the shape of the atmosphere was not clear, later developments such as in Ponte-Castañeda and Willis [START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF] made it clear that it is related to the spatial distribution of the inclusions, independently from their shape. From this generalized self-consistent set-up, several assumptions are made, which lead to the IDD estimate, as characterized by the localization tensor of inclusion i:

A IDD i = A dil i : I - k f k P 0 Dk : A dil k -1 (29)
This estimate is identical to the two-phase estimate 2ph when there is only one heterogeneity family, and both these estimates also become identical to the Mori-Tanaka MT when the atmosphere is taken similar to the inclusion itself.

It is also identical to the Ponte-Castañeda and Willis estimate in the case where the distribution of the various phases is identical [START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF]. Therefore, the properties of the PCW estimate explained in [START_REF] Weng | A dynamical theory for the Mori-Tanaka and Ponte Castañeda-Willis estimates[END_REF] can be extended to the IDD estimate when all atmospheres are identical, such as the dynamical basis of the estimate in some simple situations of aligned and randomly oriented cracks. Also, the drawbacks of the PCW estimate underlined in [START_REF] Hu | Some reflections on the Mori-Tanaka and Ponte Castañeda-Willis methods with randomly oriented ellipsoidal inclusions[END_REF][START_REF] Hu | The connections between the double-inclusion model and the Ponte Castañeda-Willis, Mori-Tanaka, and Kuster-Toksoz models[END_REF][START_REF] Weng | A dynamical theory for the Mori-Tanaka and Ponte Castañeda-Willis estimates[END_REF], such as the fact that flat ellipsoids cannot keep an isotropic spatial distribution at high volume fraction, because it leads to a violation of Hashin-Strickman bounds [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF], apply to the IDD estimate.

However, the PCW estimate is not explicit when these distributions are different for each inclusionary phase. The PCW estimate is rigorous in the sense that imposing the symmetry of pair-wise distribution functions is required by the Hashin-Strickman framework they use. In the IDD estimate, independent distribution functions can be chosen, which might violate this requirement, but makes the use of the estimate much simpler (as also stated in [START_REF] Klusemann | Homogenization methods for multi-phase elastic composites: Comparisons and benchmarks[END_REF], even if in the cases presented in this paper, these estimates coincide). This estimate won't be used directly here, but we acknowledge that it is of great importance to mention it because it allows understanding the drawbacks of other estimates, for example concerning the symmetry of distribution functions of pairs of families of different inclusions, which can be violated using the MT or IDD estimate [START_REF] Kanaun | Elastic properties of hybrid composites by the effective field approach[END_REF].

Generalized self-consistent estimate (GSC) This estimate is not in closed form. The GSC set-up explained when presenting the IDD estimate can be dealt with numerically, instead of using approximations such as those of Zheng and Du which yields the IDD estimate or Hori and Nemat-Nasser's [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF]. It is solved iteratively: at each step of calculus we have the inclusion material in ω, the matrix in D\ω, and an approximation of the homogenized stiffness tensor outside D. A new approximation of the stiffness tensor is determined by imposing a strain at infinity and finding the fourth order tensor that relates the average strains and stress on D. At next step, this new stiffness tensor is used for the infinite medium. It has been shown that this algorithm converges, and the limit stiffness tensor is the GSC estimate [START_REF] Bornert | A generalized pattern-based self-consistent scheme[END_REF].

Comparisons between the estimates and the finite element simulations

We now have all elements to test the validity of the estimates in different cases. We will begin with discs, then aligned and randomly oriented ellipses.

Discs

It is the simplest situation, in which the homogenized medium is isotropic, and the distribution of the voids is similar to their shape (see Fig. 1). The study of that type of morphology is classical, whether experimentally [START_REF] Miled | Effective elastic properties of porous materials: Homogenization schemes vs experimental data[END_REF] (in this study the differential scheme performs better than the MT scheme) or numerically [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF], where Torquato's third order estimate [START_REF] Torquato | Effective stiffness tensor of composite media: II. applications to isotropic dispersions[END_REF] is preferred to the GSC estimate, while the MT estimate is considered the poorest. Note that we use monodisperse discs, which limits the volume fraction of voids, but was shown for example by Ma (in the case of cracks) not to have a big influence on the mechanical properties [START_REF] Ma | Numerical analysis of interaction and coalescence of numerous microcracks[END_REF] and by Kurukuri for spheres [START_REF] Kurukuri | Homogenization of Damaged Concrete Meso-structures using Representative Volume Elements -Implementation and Application to SLang[END_REF].

We plot four types of results (see Fig. 4). First, simulation results obtained on the microstructures generated by Random Sequential Addition that we label alea. The periodic boundary conditions were used, and convergence of the properties was checked (see § 4.1). Second, simulations on the periodic cell (see § 4.2), labelled per. Third, two classical estimates: self-consistent SC and Mori-Tanaka MT, which is in this case identical to the 2ph and IDD estimates since there is only one void family and its distribution is similar to its shape. Finally, the GSC which was obtained numerically.

The SC scheme, when compared to the alea, as well as per results, performs poorly. It overestimates the decrease of the mechanical properties. This confirms the general opinion that this estimate is not suitable for matrix/inclusion morphologies, but for percolating morphologies, of the kind we could have obtained by using microstructures generated by an algorithm which allows for the interpenetration of the generated voids. The MT=IDD=2ph estimate gives good results, as well as the GSC numerical estimate. Interestingly, except for the shear modulus, the estimates are closer to the periodic simulation than to the random simulation. Comparison with classical estimates (SC and MT=IDD=2ph), computations on a periodic cell per, and the GSC numerical estimate.

Flat aligned pores

In this case the material is no more transversely isotropic as in the case of discs. Its overall properties are described by four coefficients in the plane. We choose two Young's moduli E 1 and E 2 , the in-plane shear modulus µ 12 and the Poisson's ratio ν 12 . In addition to that, we have two Biot coefficients and the Biot compliance (b 11 , b 22 , M ). As was mentioned, the cracks do not intersect with the outer boundary of the domain, as in simulations found in Huang et al. [START_REF] Huang | The numerical calculation of two-dimensional effective moduli for microcracked solids[END_REF]. Huang has concluded in his work that for crack densities up to 0.6, the GSC estimate is the best, which they only prove it in the case of randomly oriented cracks. However, it seems the computations domains he used were too small compared to the inclusion sizes [START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF]. It is however encouraging for us, because the IDD estimate is built to approach the GSC estimate. In Shen and Li [START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF], the differential scheme is considered satisfactory up to a crack density of 0.6, while in Yin and Ehrlacher [START_REF] Yin | Size and density influence on overall moduli of finite media with cracks[END_REF] it is shown efficient up to 0.3, but it is not tested above this values.

Fig. 5 Flat pores of aspect ratio X=1/10. The domain has the same aspect ratio as the pores, so that the error on the estimation of the poromechanical properties is similar in both in-plane directions.

We plot three sorts of curves. First, the results on the microstructures generated by Random Sequential Addition alea. Second, two sorts of simulations on periodic cells, as explained in § 4.2. perEE denotes the simulations in which the distribution of the ellipses is similar to their shape (see Fig. 2), while perEC denotes the simulation in which the ellipses remain distributed isotropically. Finally, two classical estimates: the 2ph estimate, in this case giving the same result as the IDD estimate, in which the distribution of the pores has been chosen as isotropic, and the MT estimate.

The poromechanical properties of this orthotropic medium are plotted on Figs. 6 and7. Note that the results are displayed as a function of the volume fraction of pores f , but when the ellipsoids are of aspect ratio X, a link can be made with the classically used crack density factor ǫ = f /( 4π 3 X), as in Budiansky and O'Connell [START_REF] Budiansky | Elastic moduli of a cracked solid[END_REF]. First, the 2ph=IDD estimate does not perform well at all. The Young's modulus in direction 2 is estimated as zero at a low volume fraction of 0.13, which does not match the simulation results. Similarly, the Biot coefficient b 2 becomes greater than 1, which is also thermodynamically impossible. The Poisson's ratio is not well estimated either. This failure is related to the fact that when putting the inclusion ω in the atmosphere D and respecting volume fractions, the first needs to fit in the second, which when the inclusion is flat and the atmosphere is not becomes impossible when increasing the volume fraction, exactly when we reach f = X. To solve this problem, we need to work on the choice of the shape of the atmosphere. Compared to this failure of the 2ph=IDD estimates, MT performs better, but the estimate is still poor. Concerning the simulations on the periodic cells, the two simulations yield very different results. We see that alea is close to perEC at low volume fractions, while getting closer to perEE at large ones. This is related to the process of construction of the microstructure. While generating the inclusions randomly in the domain, no anisotropy is created. However, the rejection of a new void when it intersects with a previously placed void creates anisotropy by making the distribution of the voids in space closer and closer to their shape. This can be measured [START_REF] Franciosi | Anisotropy features of phase and particle spatial pair distributions in various matrix/inclusions structures[END_REF], and tells us what to improve in the 2ph=IDD estimate to avoid the failure occurring when the distribution of flat voids becomes anisotropic at high volume fractions.

As the volume fraction increases, the assumption of isotropic distribution of the inclusions which was made using the 2ph estimate becomes less realistic than the one which is embedded in the MT estimate, in which the distribution of inclusions is similar to their shape. We need to imagine a transition from one regime to the other, inspiring ourselves of the case where the inclusions are placed on the periodic grid (see [START_REF] Zheng | An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution[END_REF][START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF][START_REF] Bary | Estimation of poromechanical and thermal conductivity properties of unsaturated isotropically microcracked cement pastes[END_REF]).

Let us come back to the auxiliary problem of the inclusion ω in the atmosphere D, and try to bring a solution to this problem in 2d. The inclusion is of aspect ratio X, the atmosphere is of aspect ratio X d , and the volume fraction of inclusion is respected f = ω/D. At low f , we start with X d = 1. When we reach f = X, the inclusion does not fit in the atmosphere anymore. Notice that the total volume fraction is used, not partial volume fractions according to the orientations. One way of modifying the atmosphere is to choose X d = X/f , so that the volume is respected thanks to a flattening of the atmosphere. The atmosphere is then, for each crack, flatten according to this rule and coaxial with the crack. This concept was already proposed by Ponte-Castañeda and Willis [START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF] in the case of aligned cracks, and Zheng and Du [START_REF] Zheng | An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution[END_REF] in the case of randomly oriented cracks, under the name Full-Range IDD (FRIDD ), and yields an estimate which lies in between the IDD and MT estimates. It is used for example in [START_REF] Stora | Influence of inclusion shapes on the effective linear elastic properties of hardened cement pastes[END_REF].

Let us see how this estimates performs with flatter ellipses, we will call cracks.

Aligned cracks

A similar analysis is made with flat ellipses of aspect ratio 1/100, which we use to represent cracks in the medium. We represent on Figs. 8 and9 the results obtained with the MT, IDD, and FRIDD estimates. The results obtained with the two first estimates have the same defaults as when we compared them to results with 1/10 aspect ratio ellipses. The FRIDD gives an intermediary solution, which is not perfect but preferable to the other estimates since it is thermodynamically correct (for example, the Young's moduli remain positive) and performs better than the MT estimate which generally overestimates the stiffnesses. However, even using this FRIDD estimate, the estimation of the Poisson's ratio is still not satisfactory.

Flat pores of random orientation, or in two orthogonal families

In this section we test two possibilities. First, the orientation of the cracks is random, as one can see on Fig. 10. Second, there is two families of orthogonal cracks as on Fig. 11. When the microstructure contains isotropically oriented cracks, the overall material is transversely isotropic. The accessible volume fractions with the RSA algorithm are smaller than in the case of aligned cracks because some spaces are too small to contain a crack of any orientation but still represent large volume fractions (see Fig. 10). The overall properties of such a microstructure are compared to the three estimates MT , IDD, and FRIDD. When the cracks have two orthogonal orientations (Fig. 11), the results are very similar, as can be seen on Fig. In this situation, the results obtained with the FRIDD estimate are very good for all the poromechanical coefficients, as we can see on Figs. 13 and14, except for the Poisson's ratio. The difficulty of estimating the Poisson's ratio is also illustrated by Segurado and Llorca [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF]. However, one can argue that the precise value of a Poisson's ratio is of little importance when the stiffness is close to zero, as it is the case at large volume fractions of cracks. Let us comment the error made on the empirical averages of these quantities. We have plot, for 3 different sizes of voids A (which is the great semi axis of an ellipse of aspect ratio 0.1) compared to domain of size L = 1, and a constant volume fraction of voids f = 0.25, the evolution of the empirical average and the standard deviation, as explained in § 5.1. We only plot the case of periodic boundary conditions. As can be seen on Fig. 15, the poromechanical properties are not fully converged when the domain size is 80 times larger than the larger dimension of the void. The rate of convergence seems similar for all quantities. However, the error made on the estimation of the empirical average, given by the error bars, is not the same for all quantities. It seems that very large domains are necessary to obtain converged values of the overall properties. It might be a drawback of the fact that there is always a zone close to the external boundary where the void density is lower, due to a wall effect. Even if the size of this zone decreases with decreasing particle size, its presence might affect the rate of convergence of the homogenized properties.

Discussion on the comparison of the estimates to the simulations

Our main conclusion is that the modified version of the IDD estimate, which was originally proposed by the authors under the name Full Range IDD, provide a good estimate of the overall poromechanical properties, when compared to simulations of 2d microstructures generated by an algorithm of Random Sequential Addition, that is in the case of matrix/inclusion morphologies. The estimate is very good in the case of randomly oriented cracks (Figs. 13 and14) for all poromechanical coefficients except for the Poisson's ratio. It also performs well in the case of aligned cracks (Figs. 8 and9), but it is less accurate, particularly concerning the Young's modulus and Biot coefficient in the direction normal to the cracks.

A side conclusion of our work is the very close results obtained when estimating the mechanical properties of cracked media containing isotropically oriented cracks on the one hand, and cracks distributed in two orthogonal families on the other hand.

Conclusion

We have built a finite element model in 2d which allows computing the overall properties of porous solids with various sorts of voids, shaped as discs, or flat ellipses. Our objective was to compare the results of these computations to classical micromechanics estimates, mainly Mori-Tanaka's and the Interaction Direct Derivative estimate of Zheng and Du, but also in some cases the Self-Consistent estimate and Generalized Self Consistent estimate. We have shown that the modified version of the IDD estimate, called Full Range IDD estimate, gives very good results in all the tested cases. This conclusion is of course restricted to the case of microstructures generated with a Random Sequential Addition Algorithm which, by forbidding the interpenetration of inclusion, creates structures which clearly have matrix-inclusion morphologies. The estimates where tested up to volume fractions as large as was available using the RSA algorithm with monodisperse voids (1 in the case of aligned cracks and 0.60 in the case of isotropically oriented cracks). Of course, to extend the validity of the estimates, one should test larger volume fractions by mixing various sizes of inclusions, as well as various sorts. Nevertheless, we think that our simple computations clearly demonstrate the efficiency of the Full Range IDD estimate for cracked media in 2d.

Therefore, we see our results as complementary to Shen and Li [START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF], since even if on some situations the estimates they uses might be more precise, the IDD estimate from Zheng and Du we chose to use here is very versatile. In our opinion it is an important property that makes it easier to use.

Fig. 2

 2 Fig. 2 Periodic cell, inclusion and ditribution of aspect ratio 0.1

Fig. 3

 3 Fig.3Evolution of the poromechanical properties with the ratio of the size of cavities to calculus domain R/L, for a volume fraction f =0,4, Young modulus E = 1 Pa and Poisson's ratio ν = 0,25, for a number of N = 50 realizations for each case, as well as the 95% confidence interval

Fig. 4

 4 Fig. 4 Evolution of the poromechanical properties as a function of the volume fraction of inclusions of voids shaped as discs, with E = 1 Pa and ν = 0.25, for N = 50 realizations for each volume fraction for alea. Comparison with classical estimates (SC and MT=IDD=2ph), computations on a periodic cell per, and the GSC numerical estimate.

Fig. 6

 6 Fig. 6 Evolution of the overall elastic properties as a function of the volume fraction of inclusions of flat pores of aspect ratio X=1/10, with E = 1 Pa and ν = 0.25, for N = 50 realizations for each volume fraction for alea. Comparison with classical estimates MT and IDD=2ph), computations on a periodic cell with isotropic distribution perEC and flat distribution perEE. Note that f=0.4 corresponds to a crack density of ǫ = 0.95.

Fig. 7

 7 Fig. 7 Evolution of the overall poroelastic properties as a function of the volume fraction of inclusions of flat pores of aspect ratio X=1/10, with E = 1 Pa and ν = 0.25, for N = 50 realizations for each volume fraction for alea. Comparison with classical estimates MT and IDD=2ph), computations on a periodic cell with isotropic distribution perEC and flat distribution perEE. Note that f=0.4 corresponds to a crack density of ǫ = 0.95.

Fig. 8 Fig. 9

 89 Fig.8Cracks of aspect ratio 1/100. Simulations, IDD estimate with circular atmosphere, MT estimate, FRIDD estimate with evolving atmosphere. Elastic properties. Note that f=0.04 corresponds to a crack density of ǫ = 0.95.

Fig. 10

 10 Fig. 10 Microstructure with isotropically oriented cracks Fig. 11 Microstructure with two orthogonal crack families

12

 12

  .

Fig. 12

 12 Fig.12Cracks of aspect ratio 1/10, split in two families of orientations 0 and π/2, or isotropically oriented.

Fig. 13 Fig. 14

 1314 Fig.13Cracks of aspect ratio 1/10, of isotropically distributed orientations, or in two families at 0 and π/2. Simulations, IDD estimate with circular atmosphere, MT estimate, FRIDD estimate with evolving atmosphere. Elastic properties. Note that f=0.25 corresponds to a crack density of ǫ = 0.60.

Fig. 15

 15 Fig. 15 Cracks of aspect ratio 1/10, of isotropically distributed orientations, volume fraction f = 0.25. Study of the convergence of the poromechanical properties determined by finite element simulations.

Acknowledgements

This work was funded by the "Chaire Durabilité des Matériaux et des Structures pour l'Énergie" at ENPC.


O Desmaison
email: olivier.desmaison@mines-paristech.fr

G Guillemot
email: gildas.guillemot@mines-paristech.fr

M Bellet

A MULTI-PHYSIC LEVEL SET APPROACH FOR THE SIMULATION OF THE HYBRID LASER / GMAW PROCESS

A new hybrid welding technique has recently been developed in order to answer to the industry needs to gather high thickness steel sheets. The combination of a laser beam and a gas metal arc enables to develop welds with a high added metal rate and low porosity. Moreover lower deformation and residual stresses are observed after cooling. Nevertheless the mastery of this technique is still in development due to the various physical phenomena which occur and interact during the welding process.

A three dimensional finite element model has been developed to simulate this welding process. Metal and air gas domains are both meshed. A Eulerian-Lagrangian approach is used in which the interface between the metal and the surrounding air or plasma is defined by a level set function. Fluid flow phenomena and temperature evolution in the weld pool are simulated.

Two moving heat sources are considered at the surface of the metal. As the arc plasma and laser beam are not modeled in the level set framework, the 'Continuum Surface Force' approach is used: a volumetric heat source distribution is applied to the immediate neighbourhood of the interface. The added metal represents an additional heat source. The Navier-Stokes equations are solved in the weld pool regarding the surface Marangoni force, the volumetric buoyancy. After solving the momentum conservation equation, metal / air interface is tracked through the resolution of a convection equation with the calculated velocity field as input.

As the hybrid welding technique is usually a multi-pass process applied on high thickness piece with high chamfer, it is of primary importance to develop a correct modeling of solid / solid contact interfaces. It is shown that the present level set approach coupled with an adaptative remeshing tool enables to follow these interfaces to simulate the entire process steps. Weld geometry on a specific steel alloy is shown and compared to the expected result.

INTRODUCTION

Welding process is largely used in the heavy industry to gather steels sheets. Recently the need to weld high thickness industrial pieces has developed. New welding conditions are needed for such pieces. For instance several passes are required in order to fill the chamfer and to develop the entire weld bead. High mechanical properties, low deformation and residual stress are expected in such process. A new hybrid welding technique has recently been proposed to answer to these needs. It combines a laser beam and a gas metal arc (Fig. 1). This hybrid laser / GMAW process enables to develop weld with a high added metal rate due to the added heat source. Moreover the final microstructure shows low defaults like porosities. This process is still in development due to the difficulties to estimate the better process parameters and to the physical phenomena interacting during the weld bead development.

A three dimensional finite element modelling is proposed to simulate this welding process. The aim of this new model is to improve the mastery of this hybrid process. A level-set approach has been developed. In a first part, this model will be presented. Its advantage for the modelling of multipasses welding will be shown. The conservation equations governing the physical phenomena are detailed. Secondly some results are presented. The first one corresponds to the full modelling. This initial result will lead to do some assumptions in the second calculation in multipasses hybrid laser / GMAW welding. This result will be detailed for various passes. The weld geometry will be discussed and compared with the expected geometry.

MODELLING

A Eulerian-Lagrangien approach is proposed to model the hybrid Laser / GMAW process. A single meshed volume is built for both metal and gas domains. A level set (LS) function is then defined corresponding to the distance, ϕ, from the metal / gas interface (Fig. 2). This approach enables to overcome the numerical difficulties induced by the metal / metal contacts in a multi pass welding. In this Level-Set approach, the entire materials properties are averaged using the distance function and a single mixed law corresponding to the Heaviside function (Eq. 1).

            + + = ε ϕ π π ε ϕ ϕ sin 1 1 2 1 ) ( H (1)
Thus, any materials property α, is defined in the transition domain as an average of local gas and metal properties, α gas and α metal as:

[ ] gas metal H H α ϕ α ϕ α ) ( ) ( 1 + - = (2)

Heat conservation equation

The temperature field evolution is calculated using the fully implicit resolution of the heat conservation equation. In non-steady state, this equation is written as:

Q T T v t T C p & r = → ∇ → ∇ - → ∇ + ∂ ∂                         λ ρ . . (3)
ρ, C p and λ are the volume mass, the equivalent specific heat and the thermal conductivity. v is the local fluid velocity and Q & is the heat source. The latent heat of transformation is taken into account following Morgan's model. Indeed the equivalent heat capacity is calculated from the enthalpy evolution as:

dt t t dt t t t p T T H H C - - - - = (4)
where t is the current time and dt is the fixed time step. The enthalpy evolution as a function of the temperature, H(T), is previously calculated for a fixed solidification path model. A program has been developed using thermodynamic equilibrium calculation based on the software Thermo-Calc (Thermocalc SA). The Gulliver-Scheil solidification path [START_REF] Gulliver | The quantitative effect of rapid cooling upon the constitution of binary alloys[END_REF][START_REF] Scheil | Bemerkungen zur Schichtkristallbildung[END_REF] has been chosen for solid fraction evolution and the corresponding enthalpy evolution has been calculated.

The volumetric heat source Q & is considered as the sum of several heat sources. The arc plasma and the defocused laser beam supply two heat sources at the metal / gas interface,

p Q & and L Q & .
Theses heat sources are modelled by a Gaussian distribution with specific parameters deduced from the process. Heat convection and radiation effects correspond to the third boundary condition at the metal / gas interface, T Q & . It is noticeable that the Level Set approach does not allow imposing specific boundary conditions at the metal / gas interface. Indeed this interface is not discretized and corresponds to the continuous curve ϕ = 0 (ie. LS = 0) (Fig. 2) which evolves through the elements. Consequently it is not possible to impose a specific heat source as a Neumann condition on this interface. The Continuum Surface Force method developed by Brackbill [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] allows to transform a boundary condition into a volumetric source term [START_REF] Ki | Modeling of laser keyhole welding: Part i. mathematical modeling, numerical methodology, role of recoil pressure, multiple reflections, and free surface evolution[END_REF][START_REF] Ki | Modeling of laser keyhole welding: Part ii. simulation of keyhole evolution, velocity, temperature profile, and experimental verification[END_REF]. Regarding a surface heat source S Q & , we can integrate this source on the area of the interface S in order to obtain the total heat:

( ) ∫ ∫ Ω = Γ V S S S d Q d Q ϕ δ & & ( 5
)
where V is the volume of the model and δ(ϕ) is the derivative of the Heaviside function, H(ϕ) (Eq. 1). This function evolves in the domain [-ε : + ε] and is equal to 0 outside this domain. This integration transforms any surface condition in volume condition. Consequently this transformation enables to impose a metal / gas interface condition in a Level-Set approach on a finite element mesh.

It should be mentioned that the average equation ( 2) is modified in the transition domain

[-ε : + ε]
. Indeed a shift of -ε/4 is introduced for the ρ C p and λ parameter evolution in order to prevent a sharp increase of the temperature. This phenomenon is due to high difference of these parameter between the two domains.

The impingement of molten droplets of filler metal in the weld pool corresponds to the last heat source. These droplets weld at the tip of the filler metal and fall into the liquid bath with a high velocity. Lancaster [START_REF] Lancaster | The physics of welding[END_REF] has initially proposed to impose this added source, D Q & , to the nodes included in a cylinder below the torch beam. This model has later been improved by Kumar and Bhaduri [START_REF] Kumar | Three-dimensional finite element modeling of gas metal-arc welding[END_REF].

D Q & and P Q & correspond to the two parts of the heat source of the GMAW process, GMAW Q & : D P GMAW Q Q I U Q & & & + = = η (6)
where U and I are the tension and the intensity of the GMAW torch and η is the efficiency of the process. The D Q & heat source is first estimated as a function of the temperature of droplets. The geometry of the corresponding cylinder depends of GMAW process parameters and liquid properties as detailed by Bellet and Hamide [START_REF] Bellet | Direct modelling of material deposit and identification of energy transfer in gas metal arc welding[END_REF]. p Q & is then estimated with the heat source GMAW Q & (Eq.6). Finally, the total volumetric heat source is equal to the sum of the various sources:

[ ] D T L P Q Q Q Q Q + - + = ) (ϕ δ (7)

Mechanical behaviour

We propose to model the fluid flow into the weld pool during the welding process. The liquid phase velocity is computed using the Navier-Stokes equation coupled with the mass conservation equation. Two forces are considered into this equation. The buoyancy force, B F r , is computed regarding the density evolution as a function of the temperature. The Marangoni force, M F r , is a surface force defined as a function of the surface tension gradient:

( ) [ ]      ∂ ∂ × ∂ ∂ = - - = s T T F T T g F M B r r r r γ β ρ 0 0 1 (8)
where g r is the gravity vector, ρ 0 is the density for the temperature T 0 and β is the thermal expansion coefficient. γ is the surface tension and s r is the vector tangent to the surface. In the level-set approach, the surface Marangoni force, M F r , is transformed in a volumetric force, vol M F r , using the Continuum Surface Force method applied in the [-ε : + ε] domain:

( ) ϕ δ M vol M F F r r = (9)
A Multi-Physic Level Set Approach for the Simulation of the Hybrid Laser / GMAW Process Finally, the set of equations is written as:

( )      = ∇ + = ∇ + ∇ + ∇ ∇ -         ∇ + ∂ ∂ 0 . . . v F F p v v v v t v vol M B T r r r r r r r r η ρ (10)
where is the liquid velocity, p is the pressure and η is the viscosity.

Surface formation

The surface evolution at the metal / gas interface is modelled with the mass conservation equation and the Navier-Stokes equation resolution regarding the gravity, g F r , and the surface tension, γ F r , forces. The forces γ F r and g F r are expressed as:

     = = g n F g r r r r γ κ γ ( 11
)
where κ is the local curvature of the interface and n r is the vector normal to the surface.

Parameters κ and n r are calculated from the distance function field, ϕ :

     ∇ ∇ = -∇ = ϕ ϕ κ n n r r . ( 12
)
Regarding these forces, the gravity is applied on the entire liquid volume domain. The surface tension should be applied on the liquid / gas interface. In the level-set approach, this force is also transformed in a volumetric force, vol F γ r , using the Continuum Surface Force method (Eq. 9). We finally solve the set of momentum and mass conservation equations as:

( )

     = ∇ + = ∇ + ∇ + ∇ ∇ -         ∇ + ∂ ∂ θ η ρ θ γ θ θ θ θ θ & r r r r r r r r v F F p v v v v t v vol g T . . . ( 13
)
whereθ & is the volume expansion source term. This source term corresponds to the added metal rate and is a function of the welding process parameter and metal properties (filler metal velocity, droplets diameter, detachment frequency of the droplets …). The nodes influenced by this source term are the same as previously defined for the heat source

D Q & .
Thus a cylinder is estimated [START_REF] Bellet | Direct modelling of material deposit and identification of energy transfer in gas metal arc welding[END_REF] in which an homogeneous material supply source term is introduced. . θ v r is the liquid velocity corresponding to the added metal rate. This velocity is specifically used to solve the convection equation of the level set [START_REF] Ville | Convected level set method for the numerical simulation of fluid buckling[END_REF] :

0 . = ∇ + ∂ ∂ ϕ ϕ θ v t r (14)
The new level set field enables to define the new liquid / metal interface (ϕ = 0) at current time regarding the added metal rate. A reinitialization method is then developed in order to respect the Eikonal property on the field ϕ. Indeed the norm of the distance function gradient has to be kept equal to one (|| ∇ ∇ ∇ ∇ ϕ || = 1) after computation of the new distance function field. This condition is usually obtained by solving a modified convection equation coupled with the transport equation [START_REF] Ville | Convected level set method for the numerical simulation of fluid buckling[END_REF]. It is proposed to develop a new topological algorithm in the present model. Indeed the coupling of the modified convection equation with equations ( 13) has not provided valid results. The present algorithm is based on the work of Elias et al. [START_REF] Elias | Simple finite element-based computation of distance functions in unstructured grids[END_REF] and is purely geometrical. The new mesh is propagated from the computed zero isovalue level set on the whole mesh step by step. An enhancement has been proposed to adapt this algorithm to heterogeneous mesh.

The resolution of (Eq. 8-13) enable to calculate both fluid flow velocity and level set position as a function of the forces that promote these fields. As presented it is proposed to separate the two resolutions. The momentum conservation equation without any added metal is first solved (Eq. 8-10) in order to estimate the fluid flow velocity. This latter velocity is then used in the heat conservation equation (Eq. 3) in the convection term.

Then the weld geometry is updated (Eq. 11-13) according to the added metal rate and surface tension force. The velocity θ v r corresponds to the weld geometry movement and the progress of the metal / gas interface into the chamfer. Figure 3 represents the successive steps of the generalize algorithm.

RESULTS AND DISCUSSION

The present approach has been developed to model first the welding process on a steel sample corresponding to grade 18MnNiMo5 in a simple case to show fluid flow development. As detailed previously the material properties have been obtained from Thermo-Calc. The thermal conductivity in solid metal, liquid metal and gas are respectively equal to 25 W.m -1 .K -1 , 35 W.m -1 .K -1 and 0.025 W.m -1 .K -1 .

Figure 4 shows the meshed domain used for the modelling of the GTAW process. The sample is a steel sheet of 15 mm thickness, 100 mm long and 50 mm wide. The total meshed domain is 20 mm height. The electrode moves in the y direction with a velocity of around 1 m.min -1 . Heat source is modelled as a moving Gaussian condition with a total power of 9500 W. Figure 5 shows the fluid flow induced in the weld pool superimposed with liquid fraction field (a-b) and temperature field (c-d). The Marangoni and buoyancy forces (Eq. 10) induce a liquid movement in the weld pool of the order of a few centimetres per second. This high velocity compared with the liquid thermal diffusivity and the small size of the weld pool corresponds to a high Peclet number into the liquid bath. Consequently the heat flow through the liquid domain is essentially due to the convection for these process parameters. It should be emphasized that a high value of the air dynamic viscosity, of 100 Pa.s is used. This value has been introduced in order to prevent shear stress at the gas/metal boundary in the weld pool. For the chosen value of the surface tension derivative (ie. 10 -4 N.m -1 .K -1 ), an inward flow is observed as shown on Fig 5 a). This flow digs the weld pool and tends to develop deeper weld beads than the one corresponding to outward flows. This preliminary result leads us to model the heat conduction into the liquid domain with an enhanced value. This approach will model the liquid mixing into the weld pool and induce a decrease of the temperature gradient in this domain, comparable to the one shown in the calculations done with the momentum conservation equation resolution. In this condition, the Navier-Stokes equation resolution with mass conservation (Eq. 10) is skipped. This approach is also proposed to save time with a low influence on the temperature field in the solid domain and the weld development. It is proposed to use an enhancement factor, f λ , of 20 as proposed by Bellet and Hamide [START_REF] Bellet | Direct modelling of material deposit and identification of energy transfer in gas metal arc welding[END_REF] compared with the liquid conductivity.

An industrial application is then proposed corresponding to the welding of a steel sample for the same composition. The piece has a rectangular chamfer of 50 mm high and around 10 mm wide. The whole piece is 90 mm high, 120 mm long and 90 mm wide. Half of the piece is shown on Figure 6. A mesh is superimposed on the sample and a Level-Set approach is developed in order to model a multipasses process in the chamfer direction. The materials properties are constants in this second computation. Nevertheless the estimated properties [START_REF] Coret | Experimental study of the phase transformation plasticity of 16mnd5 low carbon steel induced by proportional and nonproportional biaxial loading paths[END_REF][START_REF] Coret | Experimental study of the phase transformation plasticity of 16mnd5 low carbon steel under multiaxial loading[END_REF][START_REF] Vincent | Thermo-mechanical consequences of phase transformations in the heat-affected zone using a cyclic uniaxial test[END_REF] are closed to the ones computed with Thermo-Calc. A hybrid heat source moves along the chamfer direction to develop the weld bead. The power of the laser beam is around 5000 W. This first heat source is located ahead the GMAW electrode. This latter heat source has a power of around 10 000 W and is divided in plasma and droplets heat source (Eq. 6). In this computation 30 % of its power is used for the melting of the electrode metal into droplets and 70 % for the development of the arc plasma source. Both heat sources move at the previous heat source velocity of around 1 m.min -1 . The weld development is shown on figure 7 for the first pass. The liquid pool is formed ahead the weld and is progressively solidified for temperature lower than T L . The temperature field moves with the heat sources. The temperature range is close to the expected one. A refinement of the mesh ahead the weld bead has been developed in order to catch the liquid / gas interface corresponding to the isovalue zero of the LS. Mesh is unrefined on the walls of the chamfer which explains the rough appearance of the other surface on figure 7. Three passes are then shown on figure 8 for the same time of 1.5 s. These passes correspond to the first ones that are formed at the bottom of the chamfer. Each layer corresponds to two passes. Thus a pass is developed on each side of the chamfer (left then right) before the next layer is formed. As it has been previously explained, the LS approach enables to manage the metal / metal contact in the multipasses welding. The overlapping of these three passes as well as the weld bead / chamfer contact is mastered. Figure 9 shows the overlapping of the same three passes in an orthogonal section. This section is closed to the experimental ones as it has been observed. The total thickness of these passes is 4 mm which is close to the measured height on experiments. A full campaign of instrumented welding is planned in this research project in order to allow an experimental validation of temperature field and weld shape for several passes.

CONCLUSION

A three-dimensional Level-Set approach is proposed to model the hybrid Laser / Gas Metal Arc Welding process. A large domain is considered including both initial steel sheet and the gas above the sample. This approach enables to consider the metal / metal interactions and to develop a modelling of the multipasses welding. The plasma arc heat source and droplets heat source are both considered in the GMAW heat source. The momentum conservation equation is first solved with buoyancy and Marangoni forces to calculate the liquid movement in the weld pool. This equation is then solved with the gravity and surface tension forces in order to estimate the level-Set velocity and the progress of the weld bead geometry. A continuum surface force method is used to transform surface force into volumetric force for the level-set approach. It is shown that the calculation of the liquid velocity in the weld pool can be neglected.

The efficiency of this Level-Set approach is then shown on the simulation of the welding of a steel sheet with a deep chamfer. This model enables to simulate the process with excepted results for weld bead development and weld geometry. The temperature field is closed to the expected one. Future experiments will be developed in the research project to validate these results.

Fig. 1 :

 1 Fig. 1 : Description of the hybrid Laser -GMAW process with the various heat sources.

Fig. 2 :

 2 Fig. 2 :The meshed gas / metal domain and the corresponding level set function. The level set function is equal to 0 on the boundary between gas and metal (ie. ϕ = 0).

Fig. 3 :

 3 Fig.3: Schematic of the present algorithm proposed to model the hybrid Laser -GMAW welding process in a Level-Set approach. Superscript t and t+dt refer to the value of the fields at computation time t and t+dt.

Fig. 4 :

 4 Fig. 4: Meshed domain developed to model fluid flow during the planar welding of a 18MnNiMo5 steel sheet: The temperature field is superimposed on the Level Set 0 (metal/gas interface).

Fig. 5 :

 5 Fig. 5: a-b) Fluid flow movement during the welding of a 18MnNiMo5 steel sheet superimposed with liquid fraction field for transversal and longitudinal view. c-d) For the same time, temperature field and corresponding weld pool border.

Fig. 6 :

 6 Fig. 6: Half of the piece with grade 18MnNiMo5 weld in hybrid Laser-GMAW process.

Fig. 7 :

 7 Fig. 7: Weld bead development and corresponding temperature field distribution for various time: a) 1.5s, b) 3 s, c) 4.5s. The white line corresponds to the liquidus temperature.

Fig. 8 :

 8 Fig.8: The first a), second b) and third c) passes formed for the same time t = 1.5 s after each pass has started. The temperature field is superimposed at the weld bead surface. The white line corresponds to the liquidus temperature.

Fig. 9 :

 9 Fig. 9: Weld bead shape in an orthogonal section after three passes. The interface of each pass is shown and superimposed on the other.

& & & & &

ACKNOWLEDGEMENTS

This work is supported by the ANR project SISHYFE (ANR-09-MAPR-0019).