François Bignonnet
email: francois.bignonnet@univ-nantes.fr

Ghazi Hassen

Luc Dormieux

Fourier-based strength homogenization of porous media

Keywords: Strength Homogenization, Fast Fourier Transform, Von Mises, Third invariant, Porous media

come

Introduction

Strength homogenization is a valuable tool to study the impact of the microstructure morphology of heterogeneous materials on their macroscopic strength properties. Heterogeneous materials may indeed exhibit a great variety of morphologies, among others: matrixinclusion morphologies such as voided metals or clay with silica and carbonate inclusions; porous granular materials such as compacted metallic powders, sandstone or chalk. This study is a contribution to the understanding of the effect of the microstructure morphology on the macroscopic strength properties.

A large number of studies have focused on the use of micromechanics to derive constitutive models of the ductile failure of porous materials. The celebrated Gurson criterion [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth. Part I: Yield criteria and flow rules for porous ductile media[END_REF], based on the exterior kinematic approach of limit analysis on a hollow sphere, was designed to study matrix-inclusion like voided metallic materials. Nonlinear homogenization approaches based on the variational principle introduced by [START_REF] Ponte-Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF] or the estimation technique by [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli approach and its link with Ponte-Castaneda's nonlinear variational procedure[END_REF] have also been proposed to study ductile porous media (see also reviews in [START_REF] Bilger | Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis[END_REF][START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF]). These types of micro-mechanical estimates have been applied to various combinations of morphologies and local strength properties for instance by [START_REF] Dormieux | Résistance d'un polycristal avec interfaces intergranulaires imparfaites[END_REF][START_REF] Maalej | Micromechanical approach to the failure criterion of granular media[END_REF][START_REF] He | Cohesive Mohr-Coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF][START_REF] Dormieux | Micromechanical models of the strength of a sandstone[END_REF][START_REF] Barthélémy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF][START_REF] Maghous | Micromechanical approach to the strength properties of frictional geomaterials[END_REF][START_REF] Barthélémy | A micromechanical approach to the strength criterion of Drucker-Prager materials reinforced by rigid inclusions[END_REF][START_REF] Sanahuja | Résistance d'un milieu poreux à phase solide hétérogène[END_REF][START_REF] Dormieux | Microporomechanics[END_REF][START_REF] Traxl | Multi-level homogenization of strength properties of hierarchical-organized matrixinclusion materials[END_REF]. Most of these estimates do not however model third invariant effects at the macroscopic scale, except works by [START_REF] Revil-Baudard | New three-dimensional strain-rate potentials for isotropic porous metals: Role of the plastic flow of the matrix[END_REF][START_REF] Benallal | An assessment of the role of the third stress invariant in the Gurson approach for ductile fracture[END_REF], who revisited the limit analysis problem on the Gurson hollow sphere, and by [START_REF] Danas | A homogenization-based constitutive model for isotropic viscoplastic porous media[END_REF][START_REF] Mbiakop | An analytical model for porous single crystals with ellipsoidal voids[END_REF][START_REF] Mbiakop | A homogenization based yield criterion for a porous Tresca material with ellipsoidal voids[END_REF] using a "second-order" non-linear homogenization method.

These closed-form or analytic, homogenizationbased models have early been compared to numerical homogenization results on simplified 2D representations of voided porous media [START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF][START_REF] Tvergaard | On localization in ductile materials containing spherical voids[END_REF], which lead to the Gurson-Tvergaard-Needleman model. Simple 3D representations of the porous medium have then been considered, such as pores in simple cubic, body centered cubic or face centered cubic periodic configurations [START_REF] Richelsen | Dilatant plasticity or upper bound estimates for porous ductile solids[END_REF][START_REF] Mcelwain | Yield criterion of porous materials subjected to complex stress states[END_REF] as well as the hollow sphere [START_REF] Pastor | 3D-FEM formulations of limit analysis methods for porous pressure-sensitive materials[END_REF]. More representative microstructures involving multiple pores have only recently been considered [START_REF] Wicklein | Numerical investigations of the elastic and plastic behaviour of an open-cell aluminium foam[END_REF][START_REF] Bilger | Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis[END_REF][START_REF] Bilger | Bounds and estimates for the effective yield surface of porous media with a uniform or a nonuniform distribution of voids[END_REF][START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF][START_REF] Fritzen | Computational homogenization of porous materials of Green type[END_REF][START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF][START_REF] Traxl | Multi-level homogenization of strength properties of hierarchical-organized matrixinclusion materials[END_REF][START_REF] Mbiakop | An analytical model for porous single crystals with ellipsoidal voids[END_REF][START_REF] Mbiakop | A homogenization based yield criterion for a porous Tresca material with ellipsoidal voids[END_REF], but are limited to approximately 50 pores. Most of the simulations where discretized using the Finite Elements Method (FEM) and solved either by solving an incremental elasto-plastic problem or directly by resorting to second-order programming solvers. Exceptions to these are the FFT-based simulations carried out in [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF][START_REF] Bilger | Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis[END_REF][START_REF] Bilger | Bounds and estimates for the effective yield surface of porous media with a uniform or a nonuniform distribution of voids[END_REF][START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF] based on the efficient numerical method initially introduced in linear elasticity by [START_REF] Moulinec | A fast numerical method for computing the linear and non linear properties of composites[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] Most of these numerical studies are restricted to the exploration of loading modes comprising a hydrostatic part combined either to a pure shear mode or to an axi-symmetric shear mode, thus involving only two or three values of the Lode angle. Studies exploring complex stress states to reconstruct a 3D yield surface in the stress space are scarce [START_REF] Richelsen | Dilatant plasticity or upper bound estimates for porous ductile solids[END_REF][START_REF] Wicklein | Numerical investigations of the elastic and plastic behaviour of an open-cell aluminium foam[END_REF][START_REF] Mcelwain | Yield criterion of porous materials subjected to complex stress states[END_REF][START_REF] Mbiakop | An analytical model for porous single crystals with ellipsoidal voids[END_REF][START_REF] Mbiakop | A homogenization based yield criterion for a porous Tresca material with ellipsoidal voids[END_REF] and often restricted to basic periodic cells or microstructures with a limited number of pores. The work of [START_REF] Wicklein | Numerical investigations of the elastic and plastic behaviour of an open-cell aluminium foam[END_REF], involving FEM computations on actual aluminum foam microstructures obtained by micro-tomography, featuring both many pores and complex loading directions is an exception.

The lack of studies combining both complex microstructures and complex loading modes is mainly due to the high computational cost of the involved non linear problems. To overcome this limitation, a slight adaptation of the mesh-free, FFT-based method by [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] using the Augmented Lagrangian and Uzawa scheme is proposed. The adaptation relies on the use of refined discretizations of the Green operators, following [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Brisard | Combining Galerkin approximation techniques and the principle of Hashin and Shtrikman to improve two FFT-based numerical methods for the homogenization of composites[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF]. The use of these recently proposed operators allows to directly deal with the singularity in the non-linear problem, without the need to resort to a regularizing term as initially proposed by [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF]. The non linear problem can thus be solved directly for a specified loading direction (i.e. it does not require load increments), which results in a faster resolution. The proposed method does thus combine two major assets: as a FFT-based method, it is mesh free and can deal easily with complex microstructures; and it is a direct method as it does not require load increments. However, the method is currently restricted to specific local strength properties.

Thanks to this efficient method, various complex morphologies are studied; including not only voided, matrix-inclusion like materials (as in most of the above mentioned numerical studies), but also granular materials. Hundreds of realizations of a random Boolean model of spheres are considered, covering the whole porosity range. Some reconstructions of the macroscopic yield surface in a 3D stress space are proposed to assess and compare the role of both the Lode angle and the morphology.

The article is organized as follows. To start with, the background for strength homogenization is recalled in Sec. 2. The FFT-based strength upscaling method is then presented in Sec. 3 and validated in Sec. 4.2 by comparison to analytical results and finite element elasto-plastic simulations for simple microstructures with either a Green or a Von Mises strength criterion with voids or rigid inclusions. The method is finally applied in Sec. [START_REF] Bilger | Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis[END_REF].3 to study the influence of both the microstructure morphology and the Lode angle on the strength of porous materials with a Von Mises solid phase. The study namely includes a random Boolean model of spheres which aims at representing disordered granular media. Numerical results are then compared to some existing closed-form homogenization based models. Finally, an interpolation of the obtained yield surfaces with respect to the Lode angle is proposed.

Background: strength homogenization

In this section, classical definitions and results on strength homogenization are briefly recalled. The reader is referred to [START_REF] Suquet | Analyse limite et homogénéisation[END_REF][START_REF] De Buhan | Approche fondamentale du calcul à la rupture des ouvrages en sols renforcés[END_REF][START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF] for a comprehensive presentation.

Description of the strength properties

Static description

The strength properties of a material may be described by a domain G of allowable stress states σ. The domain G is assumed convex and must include the free stress state σ = 0. It can be characterized by a strength criterion f (σ):

σ ∈ G ⇔ f (σ) 0, (1)
where the function f is convex.

Kinematic description

In an equivalent and dual way [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF][START_REF] Salençon | Seismic bearing capacity of circular footings: a yield deisgn approach[END_REF], the domain G may be described by its socalled support function π defined for any direction d of the space of symmetric second order tensors by:

π(d) = sup σ∈G σ : d. (2)
The variable d can be physically interpreted as a strain rate tensor, and the support function π as a dissipated power. An allowable stress state is equivalently characterized by [START_REF] Barthélémy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF] or:

σ ∈ G ⇔ ∀d, σ : d ≤ π(d). (3)
Geometrically, this means that the convex domain G is located in the half-space σ : d ≤ π(d) for any direction d such that π(d) < +∞. By construction, the support function π has values in R + , is convex and positively homogeneous of degree 1, i.e.:

∀t ∈ R + , π(td) = tπ(d). (4)
Boundary of the strength domain If a stress state σ ⋆ belongs to the domain G and to a hyper-plane σ : d = π(d), then it belongs to the boundary ∂G of the domain G in a point where d is outward normal to ∂G :

σ ⋆ ∈ G σ ⋆ : d = π(d) ⇒    σ ⋆ ∈ ∂G ∃λ > 0, d = λ ∂f ∂σ (σ ⋆ ) (5)
In turn, a consequence of ( 4) and ( 5) is that for any direction d such that π(d) is finite, the second order tensor σ ⋆ = ∂π/∂d(d) is a stress state which belongs to ∂G :

σ ⋆ = ∂π ∂d (d) ∈ ∂G, (6)
provided that π(d) is regular enough for ∂π/∂d(d) to exist. This property will turn out to be particularly interesting to determine the macroscopic strength properties of a heterogeneous media.

Homogenization of the strength properties

Representative volume element and boundary conditions

A representative volume element (RVE) Ω of a heterogeneous material is now considered. At each point z of Ω the local strength properties are described by a local strength domain G(z), which is characterized equivalently by the strength criterion f (σ, z) or its support function π(d, z).

The macroscopic strength properties will be defined by means of a limit analysis problem on the RVE Ω. In this paper, we focus on periodic boundary conditions with imposed macroscopic strain rate D:

∀z ∈ ∂Ω, v(z) = D • z + v per (z) with v per (z) periodic σ(z) • n(z) anti-periodic ( 7
)
where v is the velocity field, σ(z) is the stress field, ∂Ω is the boundary of Ω and n(z) the outward unit normal to ∂Ω at point z. The set of velocity fields which are kinematically admissible with D is denoted K(D) and defined by:

K(D) = u | piece wise C 1 , u(z) -D • z periodic ∀z ∈ ∂Ω} . ( 8
)
where C 1 denotes the set of continuous fields with continuous derivatives. The strain rate field d = grad s u associated to a velocity field u kinematically admissible with D readily verifies the strain averaging rule:

d = D where • = 1 |Ω| Ω •(z) d V. (9)
Note that the set K(D) comprises velocity fields with discontinuities. If a field u ∈ K(D) is discontinuous across a surface S with unit normal n, then the strain rate d(u) has a singularity on S. Denoting u the velocity jump across S, the strain rate must be decomposed in the sense of distributions as:

d(u) = {d(u)} + u s ⊗ nδ S ( 10
)
where {d(u)} is the regular part of d(u) and δ S the Dirac distribution of the discontinuity surface S. The average in [START_REF] Budiansky | On the elastic moduli of some heterogeneous materials[END_REF] then involves a volume contribution from the regular part {d(u)} and a surface contribution from the discontinuity u :

d(u) = 1 |Ω| Ω {d(u)} d V + S u s ⊗ n d S . (11)
Similarly, the set of statically admissible stress fields σ with a macroscopic stress Σ is denoted S(Σ) and defined by:

S(Σ) = {σ | divσ = 0, σ = Σ, σ(z) • n anti-periodic ∀z ∈ ∂Ω} . ( 12
)
More generally, the set of statically admissible fields is denoted by S = ∪ Σ S(Σ).

Static definition of the homogenized strength

The homogenized strength domain G hom is the set of macroscopic stress states Σ such that it is possible to exhibit a microscopic stress field σ(z) which is statically admissible with Σ and allowable with the local strength domain G(z) at any point z of the RVE Ω [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF][START_REF] De Buhan | Yield design of reinforced earth walls by a homogenization method[END_REF] :

Σ ∈ G hom ⇔ ∃σ, σ ∈ S(Σ) ∀z ∈ Ω, σ(z) ∈ G(z) (13)

Kinematic definition of the homogenized strength

Similarly to (2) and ( 3), the support function Π hom (D) of the homogenized strength domain G hom is defined such that:

Σ ∈ G hom ⇔ ∀D, Σ : D ≤ Π hom (D) = sup Σ∈G hom Σ : D. (14)
Interestingly, as a consequence of the equivalence between the kinematic and static approach of yield design [START_REF] Frémond | Analyse limite. Comparaison des méthodes statique et cinématique[END_REF], it turns out that the homogenized support function further verifies:

Π hom (D) = inf u∈K(D) π(d(u), z). (15)
In the case of velocity discontinuities such as in [START_REF] Danas | A homogenization-based constitutive model for isotropic viscoplastic porous media[END_REF], the average of the support function in [START_REF] Dormieux | Résistance d'un polycristal avec interfaces intergranulaires imparfaites[END_REF] must be understood as:

π(d(u), z) = 1 |Ω| Ω π({d(u)}, z) d V + S π( u s ⊗ n, z) d S (16)
In order to realize the minimum of (15) in K(D), velocity fields such that π(d) is infinite must be avoided. A candidate velocity field is called pertinent if the associated strain rate ensures that π(d) < ∞ at any point. In practice, the optimization of ( 15) is thus carried out on the set of kinematically admissible and pertinent velocity fields.

The optimal velocity field which realizes the minimum of (15) corresponds to the free plastic flow mechanism within the framework of perfect plasticity. In what follows, it is also somehow abusively referred to as a failure mechanism.

Classical methods to bound the homogenized strength domain

The so-called interior static and exterior kinematic approaches of the yield design theory [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF] enable to build bounds on the homogenized strength domain G hom , as follows:

Interior static approach Based on the direct definition of G hom (13), the interior static approach consists in finding a statically admissible stress field σ ∈ S(Σ) and allowable with the local strength properties, which implies that the macroscopic stress Σ ∈ G hom . If this is repeated for several macroscopic stress states, their convex hull is included in G hom and thus provides an interior approximation to its boundary ∂G hom .

Exterior kinematic approach Based on the dual definition [START_REF] Dormieux | Microporomechanics[END_REF], the exterior kinematic approach consists in finding a velocity field which is pertinent and kinematically admissible with a macroscopic strain rate D. This allows to determine via [START_REF] Dormieux | Résistance d'un polycristal avec interfaces intergranulaires imparfaites[END_REF] an upper bound on the homogenized support function Π hom (D), and thus a halfspace of macroscopic stress states which comprises the entire domain G hom . If this is repeated for several directions D, the intersection of all the corresponding halfspaces is a convex domain which comprises the entire homogenized strength domain G hom and thus provides an exterior approximation to G hom . Efficient numerical implementations of these two approaches have been proposed in [START_REF] Krabbenhoft | Formulation and solution of some plasticity problems as conic programs[END_REF][START_REF] Pastor | 3D-FEM formulations of limit analysis methods for porous pressure-sensitive materials[END_REF][START_REF] Gueguin | An optimization method for approximating the macroscopic strength criterion of stone column reinforced soils[END_REF][START_REF] Bleyer | Yield surface approximation for lower and upper bound yield design of 3D composite frame structures[END_REF] to mention a few. In the next section, an alternative approach is presented.

Determination of the boundary of the homogenized strength domain

The purpose of this section is to retrieve any macroscopic stress state on the boundary ∂G hom of the homogenized strength domain as the average of a stress field solution to a fictitious non-linear viscous problem.

Fictitious non-linear viscous problem

The second definition of the homogenized support function [START_REF] Dormieux | Résistance d'un polycristal avec interfaces intergranulaires imparfaites[END_REF] indicates that Π hom (D) can be obtained by the resolution of a minimization problem of π(d(u), z) on kinematically admissible velocity fields u ∈ K(D). This minimization problem has formally the same structure as a non-linear viscous problem in which the support function is interpreted as the strain potential. Under the assumption that the local support function π is regular enough, it is appealing to consider the following fictitious viscous problem in which the unknowns are the velocity field v and the stress field σ:

         σ ∈ S v ∈ K(D) and d = grad s v ∀z ∈ Ω, σ(z) = ∂π(d, z) ∂d (17)
Following [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF], the solution to the problem [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] realizes the minimum in [START_REF] Dormieux | Résistance d'un polycristal avec interfaces intergranulaires imparfaites[END_REF]. In other words, the injection of the solution velocity field in the average over the RVE of the local support function yields the homogenized support function for the macroscopic direction D. Consequently, the average of the stress field solution to the fictitious viscous problem (17) fulfills the following properties:

   Σ = ∂π ∂d (d) ∈ ∂G hom D is exterior normal to ∂G hom at point Σ (18)
Hence, the resolution of the fictitious viscous problem [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] provides a strategy to build the boundary of the homogenized strength domain ∂G hom .

Note that the support function π is not differentiable in d = 0. The fictitious constitutive law σ = ∂π/∂d(d) in problem [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] degenerates to the condition σ ∈ G at any point where d = 0 [START_REF] Barthélémy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF].

FFT-based method

The proposed method relies on the approximate resolution of the fictitious non-linear problem (17) by a numerical method on a regular discretization grid. The strain rate field computed numerically will be used to approximate the homogenized support function and hence the homogenized strength criterion.

The proposed scheme is directly based on the Augmented Lagrangian iterative scheme proposed by [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] for the efficient resolution of non-linear elastic problems using FFTs, with minor adaptations. This FFT-based Augmented Lagrangian scheme has also recently been used by [START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF] to homogenize strength properties. However, the present method is based on the direct resolution of the fictitious non-linear viscous problem [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF], whereas the method presented in [START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF] involves the resolution of an elasto-plastic problem with load increments until the limit load, which is more time consuming.

Augmented Lagrangian

As shown previously, the strength homogenization problem may be reformulated in the minimization problem [START_REF] Dormieux | Résistance d'un polycristal avec interfaces intergranulaires imparfaites[END_REF]:

Π hom (D) = inf u∈K(D) π (d(u)). ( 19
)
This minimization problem can be relaxed by the introduction of the linear relation e(z)d (u(z)) = 0 at any point z of Ω. The Lagrange multiplier associated to this relation is denoted λ(z). The classical Lagrangian would then read:

L(d(u), e, λ) = π (e) + λ : (d(u) -e) (20)
However, it is numerically more stable to consider an alternate function, the so-called Augmented Lagrangian, defined as follows:

L 0 (d(u), e, λ) =π (e) + λ : (d(u)e)

+ 1 2 (d(u) -e) : C 0 : (d(u) -e), (21)
where the homogeneous stiffness C 0 which has been introduced is associated to a positive definite quadratic form. The initial problem is retrieved at the stationarity of the Augmented Lagrangian L 0 : L 0 (d(u), e, λ), [START_REF] He | Cohesive Mohr-Coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF] where the exchange between inf and sup is justified when π(d) is convex and has sufficient growth at infinity [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF].

Let us focus on the right hand side of [START_REF] He | Cohesive Mohr-Coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF]. The inner minimization problem over u ∈ K(D) refers to a linear problem defined on the domain Ω with a uniform stiffness C 0 and a prestress τ = λ -C 0 : e (see [START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF] and ( 55)). The solution to this problem is formally known and involves the Green operator Γ 0 of the uniform medium with stiffness C 0 and periodic boundary conditions (see [START_REF] Van Der Marck | Network approach to void percolation in a pack of unequal spheres[END_REF]). The solution strain rate field to the inner minimization problem over u ∈ K(D) is hence

d = D -Γ 0 * τ .
Subsequently, the stationarity of the Augmented Lagrangian [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth. Part I: Yield criteria and flow rules for porous ductile media[END_REF] with respect to λ and e leads to the following non-linear system:

∀z ∈ Ω,          d(u(z)) = D -[Γ 0 * (λ -C 0 : e)] (z) d(u(z)) -e(z) = 0 ∂π(e, z) ∂e -λ(z) + C 0 : (e(z) -d(u(z))) = 0 (23)
At the solution to the system [START_REF] Hill | A self consistent mechanics of composite materials[END_REF], e (resp. λ) is equal to the strain rate (resp. stress) field solution to the initial fictitious non-linear problem [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF].

Discretization

Discretization grid

We now turn to the numerical resolution of the problem (19) using a FFT-based method which were initially introduced in linear elasticity [START_REF] Moulinec | A fast numerical method for computing the linear and non linear properties of composites[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. The Lagrange multiplier fields e and λ involved in [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth. Part I: Yield criteria and flow rules for porous ductile media[END_REF] are discretized on the space of voxel-wise constant fields :

e(z) = β∈I e β χ β (z) ; λ(z) = β∈I λ β χ β (z) ( 24
)
where χ β is the function equal to 1 at any point of the voxel Ω β (and its periodic reproductions) and equal to 0 elsewhere. Hence, the discretized fields e and λ have both 6N degrees of freedom in three dimensions or 3N degrees of freedom in plane strains which correspond to the data series (λ β ) β∈I and (e β ) β∈I .

Discrete non-linear system

As in the continuous case presented in the previous section, the stationarity of the Augmented Lagrangian ( 21) is first formally written on the space of kinematically admissible velocity fields, involving the Green operator. The stationarity is then written on the two discrete data series (λ β ) β∈I and (e β ) β∈I . The arising discretized non-linear system is then:

∀β ∈ I,          d β = D -Γ 0 * (λ -C 0 : e) β d β = e β ∂π β (e β ) ∂e + C 0 : e β = λ β + C 0 : d β (25)
with the following notations

   d β = d(u) β π β (e) = π(e, z) β where • β = 1 |Ω β | Ω β •(z) d V. (26)
It is already clear that the proposed discretization is not tailored to efficiently take into account velocity fields with discontinuous mechanisms, which are sometimes optimal in strength homogenization. Indeed, the strain rate field d = D -Γ 0 * (λ -C 0 : e) is by construction the solution to an elasticity problem with a uniform stiffness C 0 and a voxel-wise constant prestress λ -C 0 : e and is thus continuous.

The first line of the discretized system (25) requires the evaluation of d β which involves the volume average over a voxel of the Green operator applied to a voxel-wise constant polarization field τ = λ -C 0 : e. This term can be evaluated exactly by resorting to the consistent discretized Green operator introduced in [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] and recalled in (61), instead of the truncated operator initially proposed in [START_REF] Moulinec | A fast numerical method for computing the linear and non linear properties of composites[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. However, the evaluation of the consistent operator involves infinite sums and is of limited practical use. Instead, simulations are carried out in the present work using the filtered, non consistent discrete Green operator latter introduced in [START_REF] Brisard | Combining Galerkin approximation techniques and the principle of Hashin and Shtrikman to improve two FFT-based numerical methods for the homogenization of composites[END_REF] and recalled in (63). The latter discrete operator has indeed been shown to provide an excellent approximation to the consistent one, while remaining computationally efficient. Note that alternative propositions for the discretization of the Green operator have recently been introduced in [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF], which could likely be substituted to the proposed choice.

Numerical resolution

Provided some requirements are met on the local support functions π β which will be discussed in Sec. 3.5, the discretized non linear system (25) can be solved by the iterative algorithm of Uzawa as proposed by [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF]. Each iteration of the Uzawa algorithm involves the three following steps:

iteration i : λ i-1 β β∈I and e i-1 β β∈I are known -Step 1. compute d i for the prestress λ i-1 -C 0 : e i-1 : d(u i ) = D-Γ 0 * λ i-1 -C 0 : e i-1 and set d i β = d(u i ) β (27)
-Step 2. ∀β ∈ I compute e i β solution to the local non linear equation:

∂π β (e i β ) ∂e + C 0 : e i β = λ i-1 β + C 0 : d i β (28)
-Step 3. ∀β ∈ I update λ i β :

λ i β = λ i-1 β + C 0 : d i β -e i β = ∂π β (e i β ) ∂e (29)
As a stopping criterion, we chose that the relative difference between two successive iterations of

β∈I π β (d i β
) is less than a tolerance (e.g. 10 -5 ). The practical implementation of Step 1 classically involves the evaluation of the convolution product in the Fourier space as follows:

1. ∀β ∈ I, compute the polarization τ i-1 β = λ i-1 β - C 0 : e i-1 β 2. evaluate the Discrete Fourier Transform (DFT) τ i-1 b b∈I of τ i-1 β β∈I 3. ∀b ∈ I apply the Fourier transform of the discrete Green operator ηi-1 b = Γ 0,b : τ i-1 b 4. evaluate the inverse DFT η i-1 β β∈I of ηi-1 b b∈I 5. ∀β ∈ I, update the strain rates d i β = D -η i-1 β
where b is the multi-index of the frequencies in the Fourier space. We recall that in the present work, simulations are carried using the filtered, non consistent discrete Green operator Γ fnc 0,b introduced by [START_REF] Brisard | Combining Galerkin approximation techniques and the principle of Hashin and Shtrikman to improve two FFT-based numerical methods for the homogenization of composites[END_REF] and defined in (63).

Since all the non-linearity is dealt with on a voxelper-voxel basis at Step 2 by ( 28), the non-linear Uzawa algorithm can be efficiently implemented. Note that an alternative FFT-based method has also been proposed in [START_REF] Gélébart | Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials[END_REF] to deal with non-linearities using a Newton-Raphson algorithm.

Further, as far as memory efficiency is concerned, the numerical implementation only requires to store in memory a single copy of the unknown vectors (λ β ) β∈I and (e β ) β∈I as well as some information on the discrete support function π β of each voxel β.

Resolution of the non-linear equation (28) in specific cases

At Step 2 of the Uzawa algorithm, the non-linear tensorial equation ( 28) can be solved by a numerical method in the general case. In this section, useful specific cases for which the non-linear equation ( 28) has analytical solutions are considered.

Case of local strength criteria depending on the two first invariants

Let us focus on the case where the local strength criterion f (σ) depends only on the two first invariants of the stress tensor:

f (σ) = f (σ m , σ d ) with σ m = tr σ/d ; σ d = √ σ d : σ d , (30)
where σ d = σ -σ m 1 (with 1 denoting the second order identity tensor) is the deviatoric part of the stress tensor and σ m is the mean (or hydrostatic) stress. The support function is then a function of the two first invariants of the strain rate tensor:

π(d) = sup{σ m d v + σ d d d | f (σ m , σ d ) ≤ 0} = π(d v , d d ) with d v = tr d ; d d = d d : d d , (31)
where

d d = d -d v 1/d
is the deviatoric part of the strain rate tensor and d v is the volume strain rate. The fictitious non-linear constitutive law may then be decomposed as:

σ = ∂π ∂d (d v , d d ) = ∂π ∂d v (d v , d d )1 + ∂π ∂d d (d v , d d ) d d d d = C(d v , d d ) : d (32)
where C(d v , d d ) is the isotropic secant rigidity tensor defined by:

C(d v , d d ) = dk(d v , d d )J + 2µ(d v , d d )K with        k(d v , d d ) = 1 d v ∂π(d v , d d ) ∂d v 2µ(d v , d d ) = 1 d d ∂π(d v , d d ) ∂d d (33)
and where k and µ are fictitious secant bulk and shear modulii, d=2 or 3 is the dimension of the problem and the volume and deviatoric projection tensors J and K are defined as J = 1 d 1 ⊗ 1 and K = I -J with I denoting the fourth order identity tensors.

The decomposition ( 32) is not unique. For example, the constitutive law could be written in a prestressed manner σ = C ′ (d) : d + σ p (d). To keep C ′ (d) isotropic, the prestress must be hydrostatic σ p (d) = σ p (d)1. Such alternative choices may prove useful to ensure that the quadratic form associated to C ′ (d) is definite positive.

Using the secant formulation (32) of the fictitious constitutive law, the non-linear equation [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF] reads:

(C(e v , e d ) + C 0 ) : e = s with s = C 0 : d u i + λ i-1 . (34)
In the case where the reference stiffness C 0 = dk 0 J + 2µ 0 K is isotropic, this tensorial equation can be projected on its volume and deviatoric parts:

(k(e v , e d ) + k 0 )e v = s m ; (2µ(e v , e d ) + 2µ 0 )e d = s ⋆ d (35)
where s ⋆ d is the positive or negative scalar such that

e d = e d s d /s ⋆ d and s d = |s ⋆ d |.
Let us now consider some cases of practical interest.

Green criterion

The Green strength criterion corresponds to:

f (σ) = σ 2 m a + σ 2 d b -1 ≤ 0. ( 36
)
In the (σ m , σ d ) plane, the boundary of the corresponding strength domain G is an ellipsis centered at the origin, with half-axes √ a along σ m and √ b along σ d . The support function of G is:

π(d) = ad 2 v + bd 2 d . (37)
The deriving fictitious isotropic secant modulii are then given by:

k(d v , d d ) a = 2µ(d v , d d ) b = 1 ad 2 v + bd 2 d . (38)
The non-linear system of equations ( 35) leads to a single scalar polynomial equation of the fourth order. Provided that the reference stiffness is chosen such that k 0 /a = 2µ 0 /b, the fourth order polynomial equation reduces to a second order one whose solution is:

e d = s d 2µ 0 1 - s 2 m a + s 2 d b -1/2 ; e v = e d bs m as d if k 0 a = 2µ 0 b (39)
The sign of e d is usually positive since s is by construction the sum of λ, which is on the boundary of the local strength domain G, and of C 0 : d, which must be directed towards the exterior of G at convergence.

Von Mises criterion

The Von Mises strength criterion corresponds to the asymptotic case a/b → ∞ of the Green criterion [START_REF] Barthélémy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF][START_REF] Dormieux | Microporomechanics[END_REF]:

f (σ) = σ 2 d b -1 ≤ 0. ( 40
)
The support function of the Von Mises criterion is:

π(d) = +∞ if d v = 0 √ bd d if d v = 0. ( 41
)
Pertinent strain rates must then be incompressible (d v = 0). This corresponds to an infinite fictitious secant bulk modulus. Without specific restrictions on the reference stiffness, the solution to the tensorial equation [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF] is:

e i (z) =      1 2µ 0 1 - √ b s d s d if s d ≥ √ b 0 if s d ≤ √ b (42)
which is pertinent since it is purely deviatoric.

Pores and rigid inclusions

The case of a voxel Ω β made up of an infinitely rigid material is simply handled by setting e i β = 0 as a solution to [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF]. In the case of a voxel Ω β made up of a pore, e i β solution to ( 28) is

e i β = d i β + C -1 0 : λ i-1 β .
Due to the use of a properly discretized Green operator, the convergence of the Augmented Lagrangian iterative scheme proposed by [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] is ensured for these infinite contrasts.

Convergence requirements

In the case of a local support function depending only on the two first invariants, an isotropic secant stiffness tensor C with bulk modulus k and shear modulus µ can be defined.

The Uzawa algorithm can then be seen as a fixed point algorithm [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF]. Indeed, at iteration i, λ i and e i are known and related by λ i = ∂π(e i )/∂e = C(e i ) : e i from [START_REF] Maalej | Micromechanical approach to the failure criterion of granular media[END_REF]. At the next iteration, equation [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF] implies:

C(e i+1 ) + C 0 : e i+1 = C : e i + C 0 : D -Γ 0 * (C(e i ) -C 0 ) : e i (43)
In the linear case, based on the eigenvalues of the linear iterative operator in [START_REF] Priour | Percolation through voids around overlapping spheres: A dynamically based finite-size scaling analysis[END_REF], the Augmented Lagrangian scheme has been shown in [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF] to converge if k 0 > 0, µ 0 > 0 and k ∈]0; +∞[, µ ∈]0; +∞[. However, if k and k 0 or µ and µ 0 have opposite signs, the convergence is not ensured. It has been further observed that the convergence can be achieved for infinite contrasts (presence of pores or rigid inclusions) when specific discretizations of the Green operator are used. These include both the consistent and filtered non consistent discretized operators introduced in [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Brisard | Combining Galerkin approximation techniques and the principle of Hashin and Shtrikman to improve two FFT-based numerical methods for the homogenization of composites[END_REF] and the discretized Green operator introduced in [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF].

This feature is of great interest for the proposed application since in a region of strongly localized deformation the fictitious secant modulii derived from the support function tend to be very soft while in regions without deformation they tend towards infinity. This is mathematically due to the fact that the support function is positively homogeneous of degree 1 (see ( 4)).

Heuristically, the Uzawa algorithm is conjectured to converge in the non linear case provided that the secant bulk and shear modulii are positive and that properly discretized Green operators are used. The proposed scheme is therefore suited to the study of heterogeneous materials governed at the micro-scale by the above mentioned Green and Von Mises criteria, as well as the Gurson criterion for example.

On the contrary, the proposed method also faces a major limitation. For example, the fictitious secant bulk modulus associated to the Drucker-Prager criterion may have negative values since the pertinent strains must have a positive volume strain while negative mean stress are allowable with the criterion. Hence, the generalization of the present method to frictional strength criteria is still an open question.

Approximation of the homogenized strength domain

At each iteration i, the continuous (i.e. non-discretized) field d(u i ) of Step 1 ( 27) is kinematically admissible with D by construction. Consequently, an upper bound Π u.b. (D) on the homogenized support function Π hom (D) could be obtained:

Π hom (D) ≤ Π u.b. (D) = π(d(u i )). ( 44
)
However this continuous strain field cannot be evaluated easily. The consistent discretized Green operator (61) would only give access to the volume average

d i β = d(u i (z))
β of the strain rate over each pixel β. Further, the use of the filtered non-consistent discretized Green operator (63) only gives access in practice to an approximation to this average. Provided the discretization is sufficiently fine, at convergence, the following approximation is used:

π(d(u i )) ≈ Π app (D) = 1 N β∈I π β (d i β ). ( 45
)
It must be emphasized that, strictly speaking, this approximation unfortunately does not preserve the upper bound status of a rigorous external kinematic approach of the yield design theory. However, this approximation behaves similarly to an external kinematic approach in so far as increasing the refinement of the discretization grid is numerically observed to provide a decreasing approximation of the homogenized support function.

In practice, estimates of exterior envelopes to the boundary of the homogenized strength domain G hom are build as follows. Several simulations are carried out for various directions of macroscopic strain rates D. The output of each simulation is the approximate macroscopic support function Π app (D). This allows to define a hyper-plane Σ : D = Π app (D) in the macroscopic stress states space whose normal is D and whose distance to the origin (or free stress state) is Π app (D). The interior convex hull of the hyper-planes obtained for all explored directions D provides then an exterior approach (up to the discretization error [START_REF] Richelsen | Dilatant plasticity or upper bound estimates for porous ductile solids[END_REF]) to the boundary ∂G hom of the homogenized strength domain. This procedure is illustrated Fig. 4 where all the hyperplanes are displayed. For all the other results presented in this document, the same reconstruction procedure has been carried out, but only the convex hull of the hyperplanes is displayed so as to improve the readability of the figures.

Unfortunately, the proposed method cannot easily be used to derive an interior static approach to the boundary of the homogenized strength domain. The latter would be desirable to quantify the distance of the numerical approximation to the exact macroscopic strength domain. Interestingly, the elasto-plastic FFTbased method proposed in [START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF] is claimed to provide an interior static approach (up to discretization errors). In the present paper, this complementary method has not been used as it is computationally more intensive.

Numerical results

This section aims both at validating the FFT-based method and at studying the effect of the microstructure morphology on the homogenized strength properties.

Considered microstructures

Morphologies

Square pore in plane strain (2D) The first considered microstructure is a square pore in a matrix obeying a Von Mises criterion [START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF] (see Fig. 1a). Plane strains are assumed for this two dimensional problem.

Face Centered Cubic and Primitive Cubic morphologies (3D) We next turn to simple periodic three dimensional microstructures: Primitive Cubic (PC) and Face Centered Cubic (FCC) lattices (see Fig. 2a). The solid phase obeys a Von Mises strength criterion [START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF]. Two cases are considered: either the interior of the spheres is the solid phase, similar to a granular material; either the interior of the spheres in the pore phase, similar to a matrix-inclusion composite. The percolation threshold, which corresponds to the maximum porosity for which the solid phase is connected, is different for each configuration: φ ≈ 26% for the FCC granular morphology, φ ≈ 48% for the PC granular, φ ≈ 96% for the FCC matrix-inclusion and φ ≈ 96.5% for the PC matrix-inclusion. Note that in the case of high porosities, the terminology "matrixinclusion" may be misleading since the pore space is interconnected: the microstructure is then closer to a sponge.

Additionally, for the matrix-inclusion morphologies, we considered the case where the matrix is made of a Green criterion [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] and the inclusions are either pores or rigid particles.

Boolean sphere model (3D) Finally, the Boolean model of spheres is considered (see Fig. 3a and3c). This model is based on a random Poisson process. Realizations of the microstructure are obtained by randomly choosing a determined number of points within the periodicity cell. These points correspond to the centers of potentially overlapping spheres. To be more specific, the domain comprising all the points which lie within a distance smaller than a given radius R of any of these sphere centers (including their periodic reproductions) is attributed to a material phase 1. The complementary domain is attributed to a material phase 2.

We restrict the study to the case where one of the material phase (solid phase) follows a Von Mises criterion [START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF] whereas the other phase is the pore space. Two situations are considered: when the material phase 1 is the solid phase, the morphology is considered as a granular material (Fig. 3c). Conversely, when the material phase 1 is the pore space, the morphology is considered as a matrix-inclusion composite (Fig. 3a).

For the granular case, the maximum porosity for which the solid phase is connected is φ ≈ 71% [START_REF] Rintoul | Precise determination of the critical threshold and exponents in a three-dimensional continuum percolation model[END_REF]. For the matrix-inclusion case, it has been estimated to φ ≈ 97% [START_REF] Van Der Marck | Network approach to void percolation in a pack of unequal spheres[END_REF][START_REF] Priour | Percolation through voids around overlapping spheres: A dynamically based finite-size scaling analysis[END_REF] to or φ ≈ 95% [START_REF] Jeulin | Percolation d'agrégats multiéchelles de sphêres et de fibres: Application aux nanocomposites[END_REF].

Note that the Face Centered Cubic (FCC) and Primitive Cubic (PC) lattices may be seen as very specific realizations of the Boolean model of spheres. However the main difference lies in the fact that the FCC and PC lattices are deterministic and organized whereas the Boolean model of spheres is random and disorganized.

Local strength properties and reference stiffness

The optimal choice of the reference stiffness is still an open question in the non linear case. It is observed to affect the convergence rate of the Uzawa algorithm [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF]. The use of the filtered non consistent discrete Green operator [START_REF] Brisard | Combining Galerkin approximation techniques and the principle of Hashin and Shtrikman to improve two FFT-based numerical methods for the homogenization of composites[END_REF] is also known to introduce a dependence on the reference stiffness of the fields at convergence. As an indication, we provide below the values chosen for the presented simulations, which empirically appeared as a good compromise.

In the case of a Von Mises solid phase with pores, the reference stiffness C 0 is chosen such that the shear modulus µ 0 = 2 √ b and the Poisson's ratio ν 0 = 0.49. In the case of a Green matrix phase [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] with pore (resp. rigid) inclusions, the reference stiffness C 0 is chosen such that the shear modulus µ 0 = √ b/5 (resp. µ 0 = 10 √ b) and the Poisson's ratio ν 0 such that 2µ 0 /k 0 = b/a.

Loading modes

Plane strain problem For the square pore problem, the study is restricted to the following imposed macroscopic strain rates:

D = cos(ψ d ) e 1 ⊗ e 1 + e 2 ⊗ e 2 2 +sin(ψ d ) e 1 ⊗ e 1 -e 2 ⊗ e 2 2 ( 46
)
These directions allow to study the boundary ∂G hom of the macroscopic strength domain in the (Σ 11 , Σ 22 ) plane.

Three dimensional problems The study is restricted to the following imposed macroscopic strain rates:

D = cos(ψ d ) 3 1 + 2 3 sin(ψ d ) -cos θ d + π 3 e 1 ⊗ e 1 -cos θ d - π 3 e 2 ⊗ e 2 + cos (θ d ) e 3 ⊗ e 3 = cos(ψ d ) 3 1 + sin(ψ d ) [cos(θ d )s 1 + sin(θ d )s 2 ] (47)
where the following orthogonal unit deviatoric tensors have been introduced:

s 1 = 2 3 e 3 ⊗ e 3 - e 1 ⊗ e 1 + e 2 ⊗ e 2 2
s 2 = e 1 ⊗ e 1 -e 2 ⊗ e 2 √ 2 (48)
The tensors (1/ √ 3, s 1 , s 2 ) constitute a basis of the space (e 1 ⊗e 1 , e 2 ⊗e 2 , e 3 ⊗e 3 ) which is orthonormal for the double contracted product " : ". The macroscopic strain rate D in [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF] is such that the volume strain rate is D v = cos(ψ d ) and the deviatoric strain rate is

D d = | sin(ψ d )|.
The angle ψ d is thus related to the triaxiality ratio of the macroscopic strain rate. In turn, the angle θ d is similar to a Lode angle and is related to the third invariant of the macroscopic strain rate J 3 (D) by:

J 3 (D) = I 3 (D d ) = det(D d ) = 3 2 D 3 d 9 cos(3θ d ) (49)
For a morphology which is isotropic or which has three planes of symmetry with normals e 1 , e 2 and e 3 , the symmetry properties of the homogenization problems for the loading modes given by [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF] imply that the macroscopic stress tensor Σ, at the point of the boundary of the homogenized strength domain with outward normal D, is diagonal in the (e 1 ⊗ e 1 , e 2 ⊗ e 2 , e 3 ⊗ e 3 ) basis [START_REF] Maghous | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique[END_REF], hence also in the (1/ √ 3, s 1 , s 2 ) basis. Consequently, the macroscopic stress state can be decomposed as follows:

Σ = Σ m 1 + Σ d1 s 1 + Σ d2 s 2 = Σ m 1 + Σ d (θ) (cos(θ)s 1 + sin(θ)s 2 ) with Σ d1 = Σ : s 1 = 2 3 Σ 33 - Σ 11 + Σ 22 2 Σ d2 = Σ : s 2 = Σ 11 -Σ 22 √ 2 ( 50
)
where θ is the Lode angle of the macroscopic stress state, which corresponds to the polar angle in octahedral planes (as in Fig. 12 and13).

Since the presented FFT-based method is formulated so as to impose a macroscopic strain rate, the parametrization of the loading direction involves the Lode angle of the macroscopic strain rates. However, after reconstruction of the failure envelope by taking the convex hull of all the hyperplanes resulting from simulations in several loading directions, the results are expressed in the macroscopic stress space. Thus all the present results are discussed in terms of effect of the Lode angle of the macroscopic stress tensor. This is an important point since in the general case the macroscopic stress and strain Lode angles differ. Note that due to the symmetry properties of the Face Centered Cubic and Primitive Cubic microstructures, the macroscopic stress and strain Lode angle are then equal for the values θ = nπ/6, where n is an integer.

Symmetry properties of the considered three dimensional problems

The symmetry properties of the local strength criterion and of the microstructure imply symmetry properties of the macroscopic strength domain. In the present case, these properties may be summarized as follows [START_REF] De Buhan | Approche fondamentale du calcul à la rupture des ouvrages en sols renforcés[END_REF][START_REF] Mcelwain | Yield criterion of porous materials subjected to complex stress states[END_REF].

First, irrespective of the microstructure geometry, the homogenized strength domain G hom inherits the following symmetry property of the local strength domains G(z):

(∀z ∈ Ω, σ ∈ G(z) ⇒ -σ ∈ G(z)) ⇒ Σ ∈ G hom ⇒ -Σ ∈ G hom , (51)
or, for the loading directions (47):

Π hom (D(ψ d , θ d )) = Π hom (D(ψ d + π, θ d )).
Second, when the local strength domain is isotropic and the directions e 1 , e 2 and e 3 are equivalent for the microstructures presented in Fig. 2a, the boundary of G hom has a 2π/3 periodicity with respect to the Lode angle θ. Equivalently, for the loading directions (47):

Π hom (D(ψ d , θ d )) = Π hom (D(ψ d , θ d + 2π/3)).
Hence, the boundary ∂G hom of the homogenized strength domain in the (Σ 11 , Σ 22 , Σ 33 ) space can be build while exploring only one sixth of the (ψ d , θ d ) directions on a unit sphere, e.g. ψ d ∈ [0; π] and θ d ∈ [0; π/3]. Further, the intersection of ∂G hom with the half-plane θ = π/6 is symmetric with respect to the plane Σ m = 0 and

Σ(Σ m , Σ d , θ = π/3) ∈ ∂G hom ⇒ Σ(-Σ m , Σ d , θ = 0) ∈ ∂G hom .

Validation of the FFT-based method

Comparison to literature results

Von Mises matrix with square pore in plane strain (2D) To start with, the numerical method is compared to exact analytical results. For the plane strain square pore problem with a Von Mises solid, in the case where the pore volume fraction (or porosity) φ is greater than 25%, the macroscopic strength domain G hom has been derived exactly in [START_REF] Maghous | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique[END_REF]:

G hom = {Σ, sup (Σ 11 , Σ 22 , |Σ 11 -Σ 22 |) ≤ (1 -φ) √ 2b (52)
The optimal kinematic mechanisms of [START_REF] Maghous | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique[END_REF] exhibit velocity discontinuities. The studied problem is thus challenging for the FFT-based method since it features three numerical difficulties that may hinder its performances:

1. due to the velocity discontinuity, the optimal strain rate is singular; 2. in the pore, the fictitious secant stiffness is zero (first infinite stiffness contrast); 3. in the matrix, the fictitious bulk modulus is infinite (second infinite contrast).

The directions of the macroscopic strain rate D (46) are discretized with a π/64 angular stride on ψ d . Note that the symmetry analysis of this square pore problem 52) (black) and exterior envelopes by the FFT-based method (red: 32 2 pixel grid, blue: 512 2 pixel grid) for the 2D square pore model (see Fig. 1a). The light red lines correspond to the hyper-planes used to build the envelope for the 32 2 pixel grid. by [START_REF] Maghous | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique[END_REF] can be expressed as follows with the notations of ( 46):

Π hom (ψ d ) = Π hom (π -ψ d ) = Π hom (-ψ d ),
so that it is sufficient to explore the angular range ψ d ∈ [0; π/2]. The exact [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli approach and its link with Ponte-Castaneda's nonlinear variational procedure[END_REF] and FFT-based approximate boundary of the macroscopic strength domain are presented in Fig. 4 for a porosity of 25%. Two discretization grids are compared: a "coarse" grid with 32 2 pixels and a "fine" grid with 512 2 pixels.

For both discretizations, the numerical envelopes reproduce effectively the angular shape of the boundary of the macroscopic domain. The relative overestimation of the strength is less than 11% (resp. 5%) for the 32 2 pixel (resp. 512 2 ) discretization grid. The 5% overestimation for the 512 2 pixel grid is somewhat disappointing as compared to the refinement of the discretization. This is attributed to the impossibility to take efficiently into account the velocity discontinuity of the optimal failure mechanism. The localization of the strain rate field is illustrated Fig. 1b for the macroscopic strain rate D = 1/2(e 1 ⊗ e 1 + e 2 ⊗ e 2 ) which corresponds to an isotropic traction in the (e 1 , e 2 ) plane. The convergence rate with grid refinement will be discussed in detail in Sec. 4.2.2.

Primitive Cubic with Von Mises matrix and pore inclusion (3D) The present FFT-based method has then been compared to the numerical results presented in [START_REF] Richelsen | Dilatant plasticity or upper bound estimates for porous ductile solids[END_REF][START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF] for the Primitive Cubic matrix-inclusion morphology with a porosity φ = 4% and a Von Mises matrix. Projections of ∂G hom on (Σ m , Σ d ) planes for stress Lode angles θ = 0 and θ = π/6 are presented in Fig. 5. [START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF] in the solid phase and 4% of porosity : comparison of the present method (lines) to the numerical results of [START_REF] Richelsen | Dilatant plasticity or upper bound estimates for porous ductile solids[END_REF] and [START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF].

The results are presented only for Σ m > 0 owing to the symmetry properties discussed in Sec. 4.1.4. The simulations have been carried out on a 64 3 -voxel discretization grid and show a good agreement with literature results, with a slight overestimation. The effect of the Lode angle is clear and will be discussed latter.

Convergence with grid refinement

The convergence of the proposed method with respect to a refinement of the discretization grid is quantitatively assessed for all of the considered morphologies.

Two dimensional results

For the two dimensional square pore problem of Fig. 1, the convergence is studied in Fig. 6a for two loading modes : a plane isotropic tensile load D = e 1 ⊗ e 1 + e 2 ⊗ e 2 and a shear load D = e 1 ⊗ e 1e 2 ⊗ e 2 for discretization grids ranging from 8 2 to 2048 2 . The effect of the choice of the reference medium is also assessed: while the reference Poisson ratio is kept to 0.49, two values of the reference shear modulus are compared: µ 0 = 2 and µ 0 = 4 (the parameter of the Von Mises criterion is here set to b = 1). Comparison elasto-plastic Finite Element Method (FEM) simulations have been carried out using the Cast3m code. The same regular discretization grids than for the FFT-based method have been used, together with quadrangular elements with either Q4 or Q8 interpolation.

Although the morphology of the square pore problem is very simple, the numerical resolution is very challenging for numerical methods which do not explicitly take into account discontinuities of the velocity field (see e.g. Fig. 1b). Consequently, both FFT and (continuous) FEM methods lead to a significant overestimation of the yield load at coarse discretization levels (see Fig. 6a). For low to medium sized discretization grids, the relative error on the yield load for the FFT simulations is comprised between the Q4 and Q8 FEM ones. For very fine discretization grids however, the FEM results are more accurate than the FFT ones for this square pore problem. The choice of the reference medium slightly affects the convergence of the FFT based method: for µ 0 = 4 the convergence rate is faster than for µ 0 = 2, but the error for low discretizations is higher.

Three dimensional results

We next study the effect of the discretization on the macroscopic isotropic tensile strength Σ max m (obtained for ψ d = 0 in ( 47)) for:

the four cubic lattice morphologies described in Fig. 2a with sphere radii such that a porosity φ = 20% is achieved, -two specific realizations of the Boolean model of spheres: a granular one with 400 spherical grains (see Fig. 3c) and a matrix-inclusion one with 50 spherical pores (see Fig. 3a). For these two realizations of the Boolean model of spheres, the sphere radii are set to 10% of the size of the periodicity cell and the achieved porosities are close to 20%.

For these six fixed geometries, the proposed FFTbased method is applied for various refinements of the discretization grid. The converged values of Σ max m are estimated so as to obtain good power-law fits of the estimated relative errors for finite discretizations (see Fig. 6b). The convergence w.r.t. discretization is the fastest for the matrix-inclusion like microstructures, with exponents ranging from 1.5 to 2.1 in the power-law fits. On the contrary, the convergence is slower for the granular microstructures, with exponents ranging from 0.8 to 1.1 in the power-law fits. The convergence rates of the FFT-based method for the 3D granular morphologies are yet better than for the previous square pore case, and compare to the convergence rate of the FEM method in the square pore case. Inspection of the local dissipation fields (see Fig. 2b, 3b and3d) reveals that in the granular cases the failure mechanisms involve clear velocity discontinuities. On the contrary, the failure mechanisms are rather bulk in the matrix-inclusion cases. Clearly, the poor convergence w.r.t. discretization of the 3D granular cases and the 2D square pore is directly related to the already discussed limitation of the FFT-based method concerning discontinuous mechanisms. Fig. 7: Adopted mesh for one eight of the unit cell of the Primitive Cubic matrix-inclusion morphology for the comparison elasto-plastic FEM simulations.

Closer inspection of the dissipation fields reveals that for the granular Boolean model of spheres, which is disorganized, the failure mechanism features both intergranular and intra-granular localized slip bands. For the FCC and PC granular morphologies, which are organized, the presence of localized slip bands depends on the direction of the load (see Tab. 1). Additionally, in most cases, there are regions of the solid phase where the strain rate of the failure mechanism is almost null, which means that the local yield limit is not reached in every point of the solid phase.

Comparison to elasto-plastic FEM : effect of Lode angle

As a last validation, the FFT-based method is compared to elasto-plastic simulations performed with a dedicated FEM code specially designed to handle elliptic criteria. The Primitive Cubic matrix-inclusion morphology is considered. The matrix phase follows a Green criterion [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] with b = 1 and a = 1 or 10. The volume fraction of the inclusions is set to φ = 10%. Two extreme kind of inclusions are considered: pores or rigid particles. The FFT-based simulations are carried out on a 64 3 discretization grid.

The FEM simulations are carried out on one eighth of the microstructure, accounting for symmetry properties. The macroscopic strain ( 47) is classically divided into small loading steps and applied incrementally to the unit cell. A kinematically admissible velocity field and a statically and plastically stress field associated to the velocity field are derived at each loading step. The ultimate load is characterised by the occurrence of a failure mechanism of the unit cell, for which the macroscopic stress tensor remains constant whereas the loading can be arbitrarily increased. The adopted mesh is composed of 375 20-noded hexahedral elements and 1991 nodes (see Fig. 7). The angular stride on the angles ψ d and θ d involved in the macroscopic strain rate ( 47) is set to π/60. The projection of the reconstructed boundaries of the homogenized strength properties in octahedral planes are presented Fig. 8. The agreement is excellent in all cases. Note however that in the case of pores, the strength properties derived by FEM are slightly greater than those derived by FFT. More precisely, the FEM results mostly over-estimate the macroscopic mean stress. The FEM simulations are indeed carried out on a rather coarse mesh, whereas by analogy with the Gurson analysis, the optimum radial velocity field in isotropic compression roughly behaves as 1/r 2 where r is the distance to the center of the sphere. This behavior is difficult to interpolate precisely with only 5 quadratic elements in the radial direction, as used in the FEM simulations (see Fig. 7). This slight difference between the FFT and FEM method in the case of pore does not appear on Fig. 8 since we compared the octahedral projections of the simulations by each method (FFT and FEM) at Σ m =0, 60, 80 and 95 % of the maximum macroscopic isotropic tensile strength Σ max m of the corresponding method, thus not exactly at the same intrinsic value of Σ m .

Study of microstructure effects for porous materials with Von Mises solid phase

In this section, the considered microstructures comprise a solid phase governed by the Von Mises criterion and a porous phase with volume fraction φ. All the kinds of three dimensional microstructure morphologies presented in Sec. 4.1.1 are compared.

For the simulations which have been carried out on random realizations of the Boolean model of spheres in this section, the spheres have a radius R = 24 pixels and the discretization grid has 240 3 voxels (this concerns Fig. 9, 10, 11b and 13). From the convergence results w.r.t. discretization presented in Sec. 4.2.2 it may be inferred that for the choice R = 24 pixels, the obtained results should lead to less than 1% of relative error for the matrix-inclusion Boolean model of spheres, but around 10% of overestimation or more for the granular Boolean model of spheres. Note that the error level likely depends on the porosity as well. Approximately 170 realizations of this random model have been generated to cover almost the whole range of porosity for both the granular and matrix-inclusion cases. Note that for this model one cannot impose a priori the porosity, but rather the intensity of the Poisson process, i.e. the number of sphere centers in the periodicity cell. For this reason, the porosity of the realized microstructures is not evenly distributed.

Porosity dependence of the maximum isotropic and deviatoric strength

Numerical results

The effect of microstructure morphology on the homogenized isotropic tensile strength and deviatoric strength is presented in Fig. 9 and 10 with respect to porosity. For the sake of clarity, the results for matrix-inclusion and granular morphologies are separated.

For the 3D cubic lattices, the discretization grid is 256 3 voxels, except for the granular Face Centered Cubic where a 512 3 -voxel grid had to be used to face convergence issues (see Fig. 6b).

For the boolean sphere model, the results present a moderate dispersion which is due to the finite size of the simulated microstructures (see e.g. [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]). As opposed to uniform stress (resp. strain) boundary conditions which introduced a softening (resp. stiffening) bias, the periodic boundary conditions are not biased for strength homogenization. Consequently, up to convergence issues due to the discretization, the results are deemed representative of the macroscopic strength properties.

The main result is that the sensitivity to morphology can be significant. As already observed by [START_REF] Kushch | Effect of microstructure on yield strength of porous solid: A comparative study of two simple cell models[END_REF] for the Primitive Cubic model, the "matrix-inclusion" and "granular" morphology classes have significantly different macroscopic strengths at the same porosity. The present results show that the differences may be marked even within a same "matrix-inclusion" or "granular" class. For instance for a matrix-inclusion composite with 30% of porosity, the isotropic tensile strength is 35% higher for the Face Centered Cubic morphology than for the boolean sphere model; and the deviatoric strength is 20% higher for the Primitive Cubic morphology than for the boolean sphere model. For the granular morphologies, the difference is even more pronounced. Although this was expected due to the marked difference in percolation threshold of the three granular morphologies, some results are counter-intuitive at first sight: the granular boolean sphere model has a percolation threshold at a larger porosity than the granular Primitive Cubic model, but its macroscopic strength is lower than the Primitive Cubic model one for porosities below 45%. A possible explanation is that, for the explored loading modes, the failure mechanism exhibits velocity discontinuities for the granular boolean sphere model whereas it is rather a bulk mechanism for the granular Primitive Cubic model.

Comparison to micromechanical models

The above numerical results are next compared in Fig. 9 and 10 to the closed-form estimates based on micromechanics presented in B. These estimates include the Gurson model (64) and Eshelby-based non-linear homogenization schemes (66). The Gurson model (64) and its modified Gurson-Tvergaard-Needleman (GTN) version [START_REF] Tvergaard | On localization in ductile materials containing spherical voids[END_REF] as well as the estimates based on the Mori-Tanaka (67) and differential (68) schemes are known to be best suited to matrix-inclusion composites, whereas the estimates based on the self-consistent scheme (69) are commonly used to deal with poly-crystalline like (here granular) composites. In Eqs. ( 67), ( 68) and ( 69), spherical inclusions are considered in the underlying Eshelby problems. For the self-consistent scheme, resorting to spherical inclusions to describe both the solid and pore phases leads to a percolation threshold φ = 1/2 from (69). Alternatively, spheroidal inclusions with an isotropic distribution of orientation have been considered in the underlying Eshelby problems for the selfconsistent scheme as an attempt to account for the difference in morphology between the solid and pore phases. The results presented in Fig. 9 and 10 illustrate that resorting to prolate inclusions with aspect ratio ω s ≥ 1 for the solid phase and oblate inclusions with aspect ratio ω p ≤ 1 for the pore phase modifies the percolation threshold of the self-consistent scheme1 .

-Considering the deviatoric strength of the matrixinclusion morphologies, the estimate based on the Mori-Tanaka scheme provides a close estimation to the deviatoric strength of the Primitive Cubic model for a wide range of porosities. The estimate based on the differential scheme is closer to the Face Centered Cubic and boolean sphere models, although it overestimates them by 8% for φ = 0.3. Both models predict the asymptotic behavior at low porosity Σ d / √ b = 1-4/3 φ+o(φ) which is in good agreement with the numerical results for the three matrixinclusion morphologies. The Gurson model systematically overestimates the deviatoric strength, but the GTN model with q 1 =1.5 is accurate for porosities up to 30% for the matrix-inclusion boolean model of spheres.

-Considering the isotropic tensile strength of the matrix-inclusion morphologies, the numerical results for the three morphologies asymptotically tend to the Gurson model at very low porosities. For the present results, the numerical results are well fitted by the GTN model in the porosity range 0 ≤ φ ≤ 0.3 with q 1 = 1.5 for the boolean sphere model, q 1 = 1.10 for the Face Centered Cubic model and q 1 = 1.18 for the Primitive Cubic model. The hydrostatic strength estimates based on the Mori-Tanaka and differential schemes are known to be ill-behaved at low porosity, but the estimate based on the differential scheme reproduces well the behavior of the Primitive Cubic and Face Centered Cubic models at large porosities. -Considering the deviatoric strength of the granular morphologies, the self-consistent scheme with spheres (ω p = ω s = 1) is in good agreement with the Primitive Cubic model. However, it fails to predict the strength of the granular boolean sphere model. In order for the self-consistent scheme to fit the granular boolean sphere model results, the aspect ratio of the prolate spheroidal solid inclusions has to be adjusted to ω s = 9.5 and the one of the oblate spheroidal inclusions to ω p = 1/7.5. Similarly for the granular Face Centered Cubic model, these aspect ratios have to be adjusted to ω s = 1 (sphere) and ω p = 1/11. Although this indicates that an information on the difference in morphology of the two phases has to be input in the self-consistent scheme, these values chosen a posteriori are surprisingly far from the spherical geometry which seemed a priori appropriate in this case. This is particularly striking for the boolean sphere model because the self-consistent scheme precisely aims at dealing with such totally disordered microstructures. Hence, the self-consistent scheme here faces a difficulty in capturing the effect of the spatial organization of the solid grains. -Considering the isotropic tensile strength of the granular morphologies, the self-consistent scheme with the above mentioned adjustments tends to overestimate the numerical results, except for the Face Centered Cubic model (which is likely a simple coincidence). As previously mentioned for the Mori-Tanaka and the differential schemes, the overestimation at low porosities is due to the use of an Eshelby-based scheme together with the modified secant method [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli approach and its link with Ponte-Castaneda's nonlinear variational procedure[END_REF]. This method indeed approximately deals with the non-linearity with phase-wise constant effective secant modulii, whereas the fluctuations of the secant modulii are important within the solid phase in the case of an imposed macroscopic isotropic traction.

Lode angle dependence

The entire boundary ∂G hom of the homogenized strength domain is now reconstructed in the (Σ 11 , Σ 22 , Σ 33 ) space for each of the four cubic lattices presented in Fig. 2a with a porosity of 20% and for a single realization of the granular boolean sphere model with a porosity close to 30%. For the cubic lattices, computations are carried out on a 128 3 -voxel discretization grid. For the granular boolean sphere model, the discretization grid has 240 3 voxels. The angular stride on ψ d and θ d for the loading direction D given by ( 47) is set to π/30. Additionally, we present for comparison the macroscopic strength domain resulting from of an exterior kinematic approach on the Gurson hollow sphere. As proposed in [START_REF] Benallal | An assessment of the role of the third stress invariant in the Gurson approach for ductile fracture[END_REF], the trial velocity fields used in this kinematic approach are solutions to incompressible linear elasticity problems defined on the hollow sphere. The approximation of ∂G hom is obtained by construction of the convex-hull of the hyper-planes Σ : D = Π app (D) (see Sec. 3.6). As the direct comparison of the obtained surfaces in the three dimensional space (Σ 11 , Σ 22 , Σ 33 ) is uneasy, we present only their intersection with half-planes θ=constant ((Σ m , Σ d ) plots in Fig. 11) or with planes Σ m =constant (octahedral planes, (Σ d1 , Σ d2 ) plots in Fig. 12 and13) with the notations defined in [START_REF] Sanahuja | Modelling elasticity of a hydrating cement paste[END_REF].

The simulated microstructure of the granular boolean sphere model is not large enough for the mechanical response to be isotropic. This is for example illustrated in Fig. 13a by the fact that the intersection of ∂G hom with octahedral planes does not feature the 2π/3 symmetry on the Lode angle. To ease the comparison with the other microstructures, an "isotropized" version of the results is also presented in Fig. 13a and11b, where the isotropy is enforced on the computed support functions by setting:

Π hom iso (D(ψ d , θ d )) = 1 3 Π hom (D(ψ d , θ d )) +Π hom (D(ψ d , θ d + 2π/3)) + Π hom (D(ψ d , θ d + 4π/3)) (53)
Projections of ∂G hom in (Σ m , Σ d ) meridional planes for Lode angles θ = 0, π/6 and π/3 are presented Fig. 11. The effect of the Lode angle on the macroscopic strength depends on the microstructure morphology: it is almost negligible for the matrix-inclusion Face Centered Cubic morphology while significant for the matrix-inclusion Primitive Cubic and granular boolean sphere model morphologies. For example, at Σ m = 80% of the isotropic tensile strength Σ max m , the macroscopic stress deviator on ∂G hom is 40% (resp. 25%) greater for the Lode angle θ = π/3 than for θ = 0 for the matrix-inclusion Primitive Cubic morphology (resp. for the granular boolean sphere model, isotropized). These apparently high differences have however to be moderated by a more appropriate measure of the effect of the Lode angle. Indeed, while considering radial loading path (imposed stress directions), the maximum relative difference of the norm between stress states on ∂G hom for θ = π/3 or θ = 0 is only around 10%.

A general trend for all morphologies is that the macroscopic strength is lower for θ = 0 than for θ = π/3 when the mean stress is positive (i.e. isotropic traction). In other words, the combination of an isotropic traction with an axisymmetric deviator due to an axial traction and a lateral compression is more harmful than the combination of an isotropic traction with an axisymmetric deviator due to an axial compression and a lateral traction. Further, the macroscopic strength for a stress deviator in the pure shear case θ = π/6 is comprised within the extremal cases θ = 0 and π/3 for high stress triaxiality ratio (i.e. mostly hydrostatic), whereas at low stress triaxiality ratio (i.e. mostly deviatoric) it is lower than for θ = 0 and π/3. For deviatoric loadings, the lower strength at θ = π/6 than at θ = 0 and π/3 is likely linked to the observation by [START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF] that the failure mechanisms are more localized in the former case than in the latter case.

These general trends are confirmed by the inspection of the projections of ∂G hom in octahedral planes for Σ m =0, 60, 80 and 95% of the maximum macroscopic isotropic tensile strength Σ max m presented Fig. 12 and13. As mentioned by [START_REF] Mcelwain | Yield criterion of porous materials subjected to complex stress states[END_REF], octahedral cuts of ∂G hom must be 3n-lobed shapes (where n is integer). The convexity of G hom then imposes that for a given mean stress Σ m , the values of the stress deviator Σ d (θ) can vary with Lode angle effects by a maximum of two. The computed projections of ∂G hom in octahedral planes evolve from 3-lobed shapes for high triaxiality ratio to 6-lobed shapes for low triaxiality ratio in order to preserve the symmetry conditions discussed in Sec. 4.1.4. This general shape of the numerical criteria is consistent with recent yield design analysis on the Gurson hollow sphere [START_REF] Revil-Baudard | New three-dimensional strain-rate potentials for isotropic porous metals: Role of the plastic flow of the matrix[END_REF][START_REF] Benallal | An assessment of the role of the third stress invariant in the Gurson approach for ductile fracture[END_REF] aiming at capturing third-invariant effects, as well as with the results of the micromechanics based variational procedure presented in [START_REF] Danas | A homogenization-based constitutive model for isotropic viscoplastic porous media[END_REF]. Interestingly, Fig. 8, 12 and 13 illustrate that projections of ∂G hom in octahedral planes can be fitted very precisely by combinations of a conic, as detailed in C. This could prove useful for the proposition of strength criteria depending on the three invariants of the stress tensor, similarly to the constitutive modeling presented in [START_REF] Mcelwain | Yield criterion of porous materials subjected to complex stress states[END_REF] based on the fitting of numerical results. This also implies that in the present case it is sufficient to compare cuts of the boundary of the macroscopic strength domain by meridional planes for Lode angles θ = 0, π/6 and π/3 to completely assess the effect of the Lode angle.

Conclusion

An efficient numerical method has been proposed for the homogenization of the strength of heterogeneous materials with periodic boundary conditions. This FFT-based method is fast and mesh free but relies on a regular discretization grid. As such, it can be directly interfaced with three dimensional imaging techniques. However, the counterpart of the regular discretization is that velocity discontinuities are not efficiently accounted for. Further, relying on the Uzawa iterative algorithm restricts the use of the method to materials with local strength domains whose associated fictive non linear viscous behaviors have positive secant modulii.

The numerical method has been validated on simple geometries by comparison to closed-form results or FEM simulations. The comparison showed an overall good behavior of the proposed method, but its convergence with respect to spatial discretization may be impaired in the case of optimal failure mechanisms which involve velocity discontinuities.

Finally, a thorough study of the effect of microstructure morphology and Lode angle of the macroscopic stress evidences that even within a same class of mor-phology such as "matrix-inclusion" or "granular" with similar constituent shapes (spheres), the precise morphology (random boolean model of spheres, Face Centered Cubic or Primitive Cubic lattices) has a more significant effect on the homogenized strength (up to ≈ 30% of variability) than the Lode angle does (up to ≈ 10% of variability).

Future works could focus on an extension of the present method to deal with frictional geo-materials (e.g. with Drucker-Prager or Mohr-Coulomb local strength criterion).

A Green operator

A.1 Continuous operator

Consider the following problem of prestressed linear elasticity: find ξ ∈ K(E), σ ∈ S such that ∀z ∈ Ω, σ(z) = C 0 : ǫ(z) + τ (z) with ǫ = grad s ξ [START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF] In this problem, ξ is the displacement field, ǫ the associated strain field and σ the stress field. The set K(E) (resp. S) is defined by (8) (resp. ( 12)). The stiffness C 0 is homogeneous, E is the imposed macroscopic deformation and τ is the socalled polarization field which is an imposed arbitrary loading parameter corresponding to a prestress field.

The displacement field ξ solution to the prestressed problem (54) minimizes the potential energy among kinematically admissible displacement fields ξ ′ ∈ K(E):

ξ = arg min ξ ′ ∈K(E) Ω 1 2 ǫ ′ : C 0 : ǫ ′ + ǫ ′ : τ d V (55)
By definition, the Green operator Γ 0 is the operator such that the strain field solution to (54) is:

ǫ(z) = E -(Γ 0 * τ ) (z)
where (Γ 0 * τ ) (z) = Ω Γ 0 (x, y) : τ (y) d V y [START_REF] Van Der Marck | Network approach to void percolation in a pack of unequal spheres[END_REF] In the Fourier space, let k denote the wave vector, k = |k| its norm and n = k/k its normalized direction. The Fourier transform of a Ω-periodic function f is denoted f and is computed using the following convention:

f (k) = 1 |Ω| Ω f (z) exp (-ık • z) d V z . (57)
The FFT-based method benefits from the property on the Fourier transform of a convolution product:

Γ 0 * τ = Γ 0 : τ (58)
When the stiffness tensor C 0 is isotropic with shear modulus µ 0 and Poisson's ratio ν 0 , the Fourier transform of the Green operator is:

Γ 0 (k) = k s ⊗ Ĝ0 (k) s ⊗ k. ( 59
)
where Ĝ0 is the Fourier transform of the Green function:

Ĝ0 (k) =    1 µ 0 k 2 1 - 1 2(1 -ν 0 ) n ⊗ n if k = 0, 0 if k = 0. (60)
The Green operator is null for k = 0 consistently with the property Γ 0 * τ ∈ K(0).

A.2 Discretized operators

The discretization of the fictitious non-linear problem carried out in Sec. 3.2 leads to compute in [START_REF] Kushch | Effect of microstructure on yield strength of porous solid: A comparative study of two simple cell models[END_REF] the average over any voxel of the strain rate d solution to (54) for a voxel-wise constant polarization field. As shown in [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF], the term involving the Green operator can be rigorously computed using the Discrete Fourier Transform as detailed in Sec. As the consistent discrete Green operator Γ c 0 involves infinite sums, it is uneasy to implement. As an alternative, a socalled filtered, non consistent Green operator Γ fnc 0 has been introduced in [START_REF] Brisard | Combining Galerkin approximation techniques and the principle of Hashin and Shtrikman to improve two FFT-based numerical methods for the homogenization of composites[END_REF] as a good approximation to the consistent discrete Green operator:

This approximate operator can be readily computed at each iteration and does not require to be stored in memory.

B Estimates of the homogenized strength of porous material with Von Mises solid phase B.1 Criteria based on an exterior kinematic approach of yield design on a hollow sphere

Gurson criterion

The Gurson criterion corresponds to an exterior kinematic approach of the yield design theory on a hollow sphere with a porosity φ and a Von Mises solid matrix [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth. Part I: Yield criteria and flow rules for porous ductile media[END_REF]. Using the notations [START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF] for the Von Mises criterion, the Gurson criterion reads:

Σ 2 d b(1 + φ 2 ) + 2φ 1 + φ 2 cosh 3 2b Σ m -1 ≤ 0 (64)
As this criterion does not take into account interaction between pore, [START_REF] Tvergaard | On localization in ductile materials containing spherical voids[END_REF] has proposed to replace φ by the correcting term q 1 φ with q 1 around 1.5.

Criterion based on trial fields from linear elasticity

While using the kinematic fields solution to the same hollow sphere with an incompressible linear elastic matrix which are richer than the trial fields used by Gurson, the derived criterion improves the Gurson criterion (see e.g. [START_REF] Benallal | An assessment of the role of the third stress invariant in the Gurson approach for ductile fracture[END_REF]). Although the expression are no longer known in closed-form, they only require a numerical volume integration and feature effects of the third invariant of the macroscopic stress tensor.

B.2 Eshelby-based criteria

Non-linear techniques relying on Eshelby-based homogenization schemes have been used to estimate the solution to the fictitious problem [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF].

In case of heterogeneous material comprising a Von Mises solid phase and pores, these techniques rely on the estimate of the homogenized stiffness C hom of a linear elastic composite with pores and incompressible solid phase with shear modulus µ s . When C hom is isotropic, C hom = 3k hom J + 2µ hom K with k hom = 2µ s x 2µ hom = 2µ s y (65)

where x and y are functions depending on morphological parameters (porosity, aspect ratio of the constituents, ...).

  plastic power for D = 1/2(e 1 ⊗ e 1 + e 2 ⊗ e 2 ) with Von Mises solid

Fig. 1 :

 1 Fig. 1: (a) Square pore in matrix in 2D with 25% of porosity. (b) Local values of the support function π(d) for the failure mechanism obtained by the FFT method. The exact optimal failure mechanism only involves velocity discontinuities.

  (a) matrix-inclusion Boolean sphere model (b) plastic power for D = 1/3 with Von Mises matrix (c) granular Boolean sphere model (d) plastic power for D = 1/3 with Von Mises grains

Fig. 3 :

 3 Fig. 3: Boolean sphere models with porosity φ ≈ 20%.

Fig. 4 :

 4 Exact boundary ∂G hom of the macroscopic strength criterion (

6 Fig. 5 :

 65 Fig. 5: (Σ m , Σ d ) cuts of the boundary of the homogenized strength criterion for the Primitive Cubic matrixinclusion model with a Von Mises criterion[START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF] in the solid phase and 4% of porosity : comparison of the present method (lines) to the numerical results of[START_REF] Richelsen | Dilatant plasticity or upper bound estimates for porous ductile solids[END_REF] and[START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. full-field simulations[END_REF].

56 FFT

 56 hydrostatic load D = e 1 ⊗ e 1 + e 2 ⊗ e 2 deviatoric load D = e 1 ⊗ e 1e 2 ⊗ e 2 fit: + 0.61N -0.46 × 1.04N -0.on macroscopic hydrostatic strength number n of pixels per sphere radius fit:• 2.23n -0.85 2.12n -0.90 10.8n -1.06• 0.70n-1.47

Fig. 6 :

 6 Fig. 6: Convergence with grid discretization for porous microstructures with Von Mises solid phase (a) 2D square pore model (see Fig. 1a) (b) 3D models for the loading direction D = 1/3.

33 Fig. 8 :

 338 Fig. 8: Homogenized strength properties of the Primitive Cubic matrix-inclusion microstructure Fig. 2a with a Green criterion (36) with b = 1 and a = 1 or 10 in the matrix phase and 10% of porosity (resp. rigid inclusions): octahedral plane cuts for Σ m = 0, 60, 80 and 95 % of the computed maximum macroscopic isotropic tensile strength Σ max m (resp. of √ a). Blue × : FEM results. Red + : FFT results. Black line: interpolation of the FFT results by combinations of a conic (see Appendix C).

Fig. 9 :

 9 Fig. 9: Homogenized deviatoric strength (for ψ d = π/2 and θ d = 0 in (47)) for various morphologies a porous materials with a Von Mises solid phase. The lines correspond to closed-form micro-mechanical models recalled in Appendix B; the dots correspond to the FFT-based simulations.

Fig. 10 :

 10 Fig. 10: Homogenized isotropic tensile strength (for ψ d = 0 in (47)) for various morphologies a porous materials with a Von Mises solid phase. The lines correspond to closed-form micro-mechanical models recalled in Appendix B; the dots correspond to the FFT-based simulations.

3 Fig. 11 :

 311 Fig. 11: Homogenized strength properties for different morphologies of porous materials with a Von Mises criterion (40) in the solid phase and comparison to closed-form estimates. (a) periodic cubic lattices of Fig. 2a. (b) granular boolean sphere model. The green, cyan and yellow lines correspond to different values of the Lode angle θ and have been obtained by the present FFT-based method.

Fig. 12 :

 12 Fig. 12: Homogenized strength properties of the periodic cubic lattices Fig. 2a with a Von Mises criterion (40) in the solid phase and 20% of porosity: octahedral plane cuts for Σ m = 0, 60, 80 and 95 % of the maximum macroscopic isotropic tensile strength Σ max m . Red + : FFT results. Black line: interpolation by combinations of a conic (see Sec. C).

22 ( 3 - 22 ( 3 Fig. 13 :

 22322313 Fig. 13: Homogenized strength properties of the granular boolean sphere model Fig. 3c with a Von Mises criterion (40) in the solid phase and ≈ 30% of porosity: octahedral plane cuts for Σ m = 0, 60, 80 and 95 % of the maximum macroscopic isotropic tensile strength Σ max m . Red + : FFT results. Black line: interpolation by combinations of a conic (see Sec. C). (a) raw results. (b) isotropized results using the procedure (53).

2 ×

 2 Γ 0 (k b 1 +n 1 N 1 ,...,b d +n d N d ), (61) where for any multi-index b = b 1 , ..., b d of Z d the wave vector k b is defined by:k b = 2πb 1 L 1 e 1 + ... + 2πb d L d e d . (62)with the notations introduced in Sec. 3.2.

2 ×

 2 Γ 0 (k b 1 +n 1 N 1 ,...,b d +n d N d ),

  , ..., L d in the d directions of an orthonormal basis (e 1 , ..., e d ). A regular discretization grid of size N = N 1 × ... × N d dividing the domain Ω is introduced. The grid is made up of N identical voxels (or pixels in 2D) which each occupy a domain Ω β where β = (β 1 , ..., β d ) is a multi-index with β i ∈ {0, 1, ..., N i -1}. The set I = {0, 1, ..., N 1 -1} × ... × {0, 1, ..., N d -1} is a subset of N d comprising all the multi-index β of all voxels.

	The domain Ω is
	supposed to be a rectangular parallelepiped of the d-
	dimensional space (d = 2 : plane strains, d = 3 : three
	dimensional problem). Its dimensions are denoted by
	L 1

Table 1 :

 1 ⊗ e 3 + e 3 ⊗ e 2 intra & inter-granular localization inter-granular localization e 3 ⊗ e 3e 2 ⊗ e 2 inter-granular localization bulk 2e 3 ⊗ e 3e 2 ⊗ e 2e 1 ⊗ e 1 inter-granular localization Effect of the loading direction on the type of failure mechanism of the cubic lattices with a Von Mises solid phase and pores.

	macroscopic strain rate D	FCC granular morphology	PC granular morphology
	e 2 bulk
	1	bulk + corner localization	bulk + corner localization

The specific case ω s = ω p = 1 corresponds to spherical inclusions for both phases and is given by (69).

Based on the modified secant method for non-linear homogenization techniques [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli approach and its link with Ponte-Castaneda's nonlinear variational procedure[END_REF], the homogenized strength properties may then be estimated by a Green strength criterion [START_REF] Dormieux | Microporomechanics[END_REF]:

Values of x and y are recalled below for classical homogenization schemes (see [START_REF] Dormieux | Microporomechanics[END_REF] for details):

-Mori-Tanaka scheme [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] with spherical pores

-Differential scheme [START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF] with spherical pores

-Self-consistent scheme [START_REF] Hill | A self consistent mechanics of composite materials[END_REF][START_REF] Budiansky | On the elastic moduli of some heterogeneous materials[END_REF] with spherical pores and solid grains

The expressions of x and y for the self-consistent scheme with spheroidal solid grains with aspect ratio ω s and spheroidal pores with aspect ratio ω p having both an isotropic distribution of orientation (see e.g. [START_REF] Sanahuja | Modelling elasticity of a hydrating cement paste[END_REF]) are implicit functions of φ, ω s and ω p , too long to be reproduced here.

C Interpolation of octahedral projections of ∂G hom by combinations of a conic

Projections of the boundary of the macroscopic strength criteria by octahedral planes Σ m =constant can be excellently approximated by a combination of conics (ellipse or hyperbola). In such (x = Σ d1 , y = Σ d2 ) plane (refering to (50)), the equation of an arbitrary ellipse or hyperbola of center (x 0 , y 0 ) is:

The conic ( 70) is an ellipse if AB -C 2 < 0 and AB > 0.

The parameters A, B, C, x 0 and y 0 are deduced from the values of the stress deviator Σ d at three Lode angles θ: R 14). Imposing that the tangent to the conic in θ = 0 and π/3 is orthogonal to the radial vector in polar coordinates, geometrical considerations lead to the unique expression of the conic parameters given in equation (71).

For Lode angles θ ∈ [0; π/3], the interpolation of the macroscopic stress deviator is then explicitly given by equation (72).

Due to the symmetry properties exposed in Sec. 4.1.4, the whole cut of the boundary of the macroscopic strength criteria by the octahedral plane Σ m =constant is constructed first by symmetry with respect to the Σ d2 axis and then by periodic reproduction with a 2π/3 period. The convexity of the strength domain imposes the conditions:

In the case where R 2 is chosen as a function of R 1 and R 3 such that y 0 vanishes and C → ∞, the expression (72) corresponds to the Lode angle dependence of the Willam-Warnke criterion [START_REF] Willam | Constitutive models for the triaxial behavior of concrete[END_REF]. In this case, a circle is obtained if R 1 = R 3 or an equilateral triangle (as for the Rankine criterion) is obtained if

, the expression (72) corresponds to a regular hexagon and hence to the Lode angle dependence of the Tresca and Mohr-Coulomb criteria. The case R 1 = R 2 = R 3 corresponds to a circle.

Excepted in the limit cases R 1 = R 3 = 2R 2 / √ 3 (hexagonal shape) and 2R 1 = √ 3R 2 = R 3 (triangular shape), the obtained criterion is smooth: the tangent is continuously varying.