Clément Cadet

Jacques Besson

Sylvain Flouriot

Samuel Forest

Pierre Kerfriden

Laurent Lacourt

Victor De Rancourt

Strain localization analysis in materials containing randomly distributed voids: Competition between extension and shear failure modes

Keywords: Ductile fracture, Strain localization, Cell computation, Porous materials, Anisotropy

Ductile fracture is studied by simulating the failure of random microstructures.

• Rice's criterion allows to detect localization and failure mode (shear or extension).

• A strong anisotropy of failure behavior is observed when rotating loading.

• This intrinsic anisotropy is stronger in single void cells than in random cells.

• Minimal failure strain curves with low ductility in generalized shear are obtained.

Introduction

Ductile fracture is a major failure mode for industrial metallic parts subjected to monotonic loading. Understanding this type of fracture is crucial for an accurate prediction of the part's resistance and therefore efficient design. Simulating ductile fracture remains a challenging problem [START_REF] Boyce | The Sandia Fracture Challenge: blind round robin predictions of ductile tearing[END_REF][START_REF] Boyce | The second Sandia Fracture Challenge: predictions of ductile failure under quasi-static and moderate-rate dynamic loading[END_REF]. It involves various small scale phenomena such as void nucleation, growth and the final failure by void coalescence or void sheeting [START_REF] Besson | Local Approach to Fracture[END_REF], all strongly dependent on the material behavior and loading type.

In order to accurately simulate ductile fracture, analytical models have been developed, such as [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] seminal one. This model, enriched by [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF], is now widely used for simulating failure. This has allowed the development of a class of models, which continues to be developed by incorporating the effects of other physical mechanisms, such as coalescence [START_REF] Benzerga | Effective Yield Criterion Accounting for Microvoid Coalescence[END_REF]; [START_REF] Morin | A unified criterion for the growth and coalescence of microvoids[END_REF]; [START_REF] Torki | A unified criterion for void growth and coalescence under combined tension and shear[END_REF] among others) or strain localization [START_REF] Keralavarma | Ductile failure as a constitutive instability in porous plastic solids[END_REF]. The effect of void shape (e.g [START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF]) or of random distributions of voids [START_REF] Leblond | A theoretical approach of strain localization within thin planar bands in porous ductile materials[END_REF]Danas andPonte Castañeda, 2009, 2012;[START_REF] Vincent | Porous materials with two populations of voids under internal pressure: I. Instantaneous constitutive relations[END_REF] can also be represented by such models. Following [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF], such models are frequently compared to results from simulations on unit cells, in which an elastoplastic cubic or parallepipedic matrix containing a single void is simulated up to failure. Such simulations can also be used on their own to help identify and distinguish failure mechanisms.

However, single void unit cells correspond to a periodic array of voids, and do not represent accurately interactions between voids, especially at high porosity. This interaction should therefore be investigated by considering more general configuration of pores. Qualitative understanding of the interaction of pores was achieved by simulations involving a cluster of a small number of voids [START_REF] Bandstra | On the influence of void clusters on void growth and coalescence during ductile fracture[END_REF][START_REF] Tvergaard | Effect of void cluster on ductile failure evolution[END_REF][START_REF] Tvergaard | Nucleation from a cluster of inclusions, leading to void coalescense[END_REF][START_REF] Trejo Navas | Void growth and coalescence in a three-dimensional non-periodic void cluster[END_REF]. [START_REF] Shakoor | Ductile fracture of a metal matrix composite studied using 3D numerical modeling of void nucleation and coalescence[END_REF] described precisely the mechanisms of ductile fracture in random populations of voids. However these studies can not provide failure models as they consider too few loading conditions. [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] used simulations on microstructures with a random distribution of voids to identify a Gurson-Tvergaard-Needlemann yield criterion. This approach was extended to other material behaviors [START_REF] Fritzen | Computational homogenization of porous materials of Green type[END_REF] or other void populations [START_REF] Khdir | Computational homogenization of plastic porous media with two populations of voids[END_REF][START_REF] Khdir | A computational homogenization of random porous media: Effect of void shape and void content on the overall yield surface[END_REF]. Yet such studies only deal with yield behavior and not the ultimate failure of the microstructure. Recently, [START_REF] Hure | Yield criterion and finite strain behavior of random porous isotropic materials[END_REF] and [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] simulated random microstructures up to coalescence and investigated the scatter of failure strain due to the random distribution of voids. A significant scatter was found even for large populations. Moreover, by considering the influence of the Lode parameter, [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] identified differences in failure behavior between unit and random cells. The latter can more freely express shear-dominated failure modes. Therefore random microstructures can better represent failure mechanisms and they should be favored for determining failure models.

A significant issue of cell studies is to operationally detect the onset of failure. Many indicators have been proposed, focusing on different physical mechanisms and characteristic features (see review by Zhu et al. (2020a) and discussion by [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]). These different choices are not equivalent and a given microstructure's resistance (strain at failure) can significantly depend on the considered indicator. A traditional approach to failure is If a full response surface for ductile fracture is desired, sufficiently general loading conditions should be used. For an isotropic material model subjected to a proportional loading, stress triaxiality and Lode parameter are mainly used. Yet conflicting results with respect to the influence of the Lode parameter on the failure strain are reported in literature. [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] and [START_REF] Luo | On the prediction of ductile fracture by void coalescence and strain localization[END_REF] found minimal ductility for generalized shear. On the other hand failure strain grows from generalized tension to generalized compression, without a minimum in shear, according to [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF]; Zhu et al. (2020a). In [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF], the behavior was mostly of the former type for random microstructures, and mostly of the latter type for unit cells. The discrepancy, which can be interpreted as a strong indication of the intrinsic anisotropy of the simulation cell, is due to different ways of applying the loading conditions, as pointed out by Zhu et al. (2020a) and [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]. However, a reliable failure model should not present this anisotropy and this dependence on the methodology. In order to avoid this effect, [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], [START_REF] Dunand | Effect of Lode parameter on plastic flow localization after proportional loading at low stress triaxialities[END_REF] and [START_REF] Tekoglu | On localization and void coalescence as a precursor to ductile fracture[END_REF] carried out the simulations by rotating the principal directions of the applied stress. The failure strain of the cell is then found by minimizing over all loading conditions with the same stress invariants. A minimal failure strain for generalized shear is then found. However only one axis of rotation was considered in these studies. This choice is motivated by [START_REF] Rudnicki | Conditions for the localization of deformation in pressuresensitive dilatant materials[END_REF] results on localization, which show that localization bands in homogeneous materials should be orthogonal to the eigendirection associated to the middle principal stress. However it was never verified numerically for unit cells. Moreover the situation for random microstructures could be more complex than for unit cells, which have higher symmetry. Fully general three-dimensional loading orientations, i.e rotations of the principal axes of the applied stress with respect to the cell's axes, should therefore be considered.

The present study therefore aims to investigate ductile fracture and more precisely strain localization behavior in random microstructures, while fully accounting for the effect of the loading orientation. To the authors' best knowledge, this is the first time that localization is studied in cells with a distribution of voids. Moreover, no previous studies had considered the effect of loading orientation with such generality. To this end, we generate cells consisting of a perfectly plastic matrix embedding a random distribution of identical spherical pores. Using the finite element software [START_REF] Zset | Zset[END_REF], simulations are performed in a large strain formulation. Various proportional loading conditions are applied up to the failure of the cell. Following the methodology by Zhu et al. (2020a), this failure is detected by Rice's localization criterion, computed at the global scale of the cell with a macroscopic tangent operator. The failure behavior for random microstructures and unit cells is then compared depending on the loading conditions. It is shown that Rice's criterion successfully describes not only localization onset, as performed by Zhu et al. (2020a), but also failure mode. A region of low ductility is observed in generalized shear. A strong anisotropy of failure behavior is observed for unit cells but is reduced for random microstructures. Minimal failure strains depending on the loading conditions can then be determined.

This paper is organized as follows. We first describe the methodology for generating random microstructures, applying loading conditions, and detecting localization with Rice's criterion. Secondly, the performance of the criterion is assessed on simplified loading conditions, in which the principal axes of the applied stress coincide with those of the cubic cell. Thirdly, more general loading orientations are investigated. Finally, with the help of further computations, we discuss the preceding results, especially the failure indicator used, and the impact of cell anisotropy. Appendices present the numerical validation of the failure indicator and the effect of the simulation parameters (mesh size, and temporal discretization).

An intrinsic notation is used for tensors: vectors are represented as v = v i e i , second order tensors as A ∼ = A i j e i ⊗ e j , where (e i ) is an orthonormal frame. The subscript 0 in the notation A 0 refers to the value of A in the initial configuration at time t = 0. The position of a material point initially at x 0 evolves with time t as x = Φ (x 0 ,t); the deformation gradient is then defined as

F ∼ = ∂ Φ ∂ x 0 .
Quantities decorated with an overlying bar, such as Ā, refer to the macroscopic counterpart (at the level of a cell) of a quantity A defined locally. For instance F ∼ is the average deformation gradient over the cell.

Methodology of micromechanical finite element simulations

This section presents the general methodology for micromechanical simulations, carried out with Zset software [START_REF] Besson | Object-Oriented Programming Applied to the Finite Element Method Part I. General Concepts[END_REF][START_REF] Zset | Zset[END_REF]. The generation and meshing of the random microstructures, the material behavior used for the matrix, the boundary and loading conditions are described. This methodology follows initially [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] one, so it will be here only summarized. The subsection 2.3 presents however an extension of the previous methodology allowing a greater variety of loading conditions when compared to [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] work.

Generation of microstructures and meshing

The microstructures are made of a periodic population of N voids identical nonoverlapping spherical voids of radius r within a cubic matrix of size a cube . The total initial porosity

f 0 = 4π 3 N voids (r/a cube ) 3 (1)
and the number of voids N voids determine together the radius of the spheres. A sphere intersecting the cube is copied on the other side (two copies for an intersected face, four for an edge and eight for a vertex). The periodic microstructure with its spherical voids is thus able to pave space. The positions of the spheres are chosen randomly according to a simple iterative dart-throwing process: A position of the center is chosen uniformly on the cube, and the new sphere is added to the collection if it (or its periodic copies) does not intersect any already accepted sphere.

The microstructure is then meshed by the NETGEN software [START_REF] Schöberl | NETGEN -An advancing front 2D/3D-mesh generator based on abstract rules[END_REF], ensuring finer elements close to the voids. The validity of the meshing parameters for ductile fracture analysis was already verified by [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]. Most of the simulations on random microstructures will be performed on the cells R1 and R2 visible in Fig. 1a and1b. They were constructed with a porosity of f 0 = 6% and N voids = 27 voids and have already been used in [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] study. A unit cell with a single void, is also used for comparison (Fig. 1c). All meshes use quadratic elements with reduced integration.

Material behavior law at finite strain

Simulations use a finite strain framework to account for the large deformation than can take place in the matrix. The strain rate D ∼ and Cauchy stress σ ∼ tensors are convected in a corotational frame [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]:

ė ∼ = R ∼ T D ∼ R ∼ s ∼ = JR ∼ T σ ∼ R ∼ (2)
where R ∼ is a rotation matrix verifying -R ∼

T Ṙ ∼ = Ṙ ∼ T R ∼ = skew( Ḟ ∼ F ∼ -1
) (material spin tensor). The matrix is chosen elastic-plastic with isotropic elasticity and von Mises plasticity without hardening:

ė ∼ = ė ∼ e + ė ∼ p e ∼ e = 1 + ν E s ∼ - ν E (tr s ∼ )1 ∼ s vm = 3 2 s ∼ dev : s ∼ dev f (s ∼ ) = s vm -R 0 ≤ 0 ė ∼ p = ṗ∂ f ∂ s ∼ (3)

Boundary and loading conditions

Periodic boundary conditions are applied on the sides of the cube [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. The displacement field u should therefore have the form:

u = ( F ∼ -1 ∼ ) • x 0 + v (x 0 ) = E ∼ • x 0 + v (x 0 ) (4)
with F ∼ the average deformation gradient, E ∼ = F ∼ -1 ∼ and v a periodic displacement fluctuation field with zero average gradient over the cell. The periodicity of v and the antiperiodicity of traction vectors mean that:

v (x + 0 ) = v (x - 0 ) (5) σ ∼ • n (x + 0 ) = -σ ∼ • n (x - 0 ) (6)
if x + 0 and x - 0 represent two homologous points on opposite sides of the periodic mesh and n (x 0 ) represent the outward-pointing normal to the mesh boundary at x 0 . These periodic boundary conditions are numerically imposed by multi-point constraints on homologous nodes. In this formulation (already used by [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF] and [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]), the degrees of freedom (DOF) are the three components of the displacement fluctuation field for each node of the mesh and the nine components of E ∼ ). The components of the macroscopic deformation gradient therefore appear as explicit degrees of freedom. In order to control rigid body motion, a node of the mesh is fixed, and the tensor F ∼ (or equivalently E ∼ ) is kept upper triangular:

E 21 = E 31 = E 32 = 0 (7)
This condition was investigated by [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] and found to yield equivalent results to the condition imposing F ∼ to remain symmetric [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. In order to facilitate the ulterior computation of Rice's criterion, this condition is applied by penalization. For (i, j) ∈ {(2, 1), (3, 1), (3, 2)}, a high stiffness spring connects the DOF E i j to a fixed one Êi j = 0. The force applied on the DOF E i j is then:

R E i j = k( Êi j -E i j ) (8)
where k is a stiffness (chosen sufficiently high compared to the stiffness of the simulated cube).

The macroscopic Boussinesq (or first Piola-Kirchhoff) S ∼ and Cauchy stress tensors σ ∼ are defined by:

S ∼ = 1 V 0 V 0 S ∼ dV 0 = 1 -f 0 V 0 V matrix 0 S ∼ dV 0 (9) σ ∼ = 1 J S ∼ . F ∼ T (10)
where J = det F ∼ and V 0 is the volume of the cell (matrix and defects) in the initial configuration. The integration on V 0 considers that stress is well-defined and identically zero in the voids. The components of S ∼ also correspond, up to a factor V 0 , to the reaction forces conjugate to the degrees of freedom E i j .

The loading conditions are applied by an extension of [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF] method. A special spring element with ten degrees of freedom is used: Nine are connected to the E i j degrees and the last one, Ê11 , corresponds to the displacement of a ghost node. This element guarantees the following form for the macroscopic Boussinesq stress tensor:

S11 = k( Ê11 -E 11 )/V 0 (11) S ∼ = J σ11 η ∼ F ∼ -T (12)
where k is a stiffness which can be chosen equal to that of Eq. ( 8), and η ∼ is a constant symmetric tensor with η 11 = 1. Eq. ( 12) can be rewritten in terms of E i j , η i j and Ê11 and K by using Eq. ( 11) and explicitly inverting F

∼ . The derivation is straightforward but tedious so the resulting expressions were obtained with the SymPy computer algebra software [START_REF] Meurer | SymPy: symbolic computing in Python[END_REF] and are not reproduced here. With this spring element, the average Cauchy stress tensor remains proportional to the constant tensor η ∼ throughout the simulation:

σ ∼ = σ11 η ∼ (13)
The loading is then driven by E 11 , which essentially coincides with Ē11 due to the high stiffness k. Zhu et al. (2020a) and [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] studied the case where η ∼ is diagonal (i.e. η 12 = η 13 = η 23 = 0). [START_REF] Barsoum | Rupture mechanisms in combined tension and shear-Micromechanics[END_REF], [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] and [START_REF] Liu | Void behaviors from low to high triaxialities: Transition from void collapse to void coalescence[END_REF] used conditions equivalent to η 22 = η 33 and η 12 = η 23 = 0. Following [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF] example, several studies have considered rotating the applied stress σ ∼ around e 2 . The problem is here considered with more generality, with a potentially full η ∼ tensor. As η ∼ is a symmetric tensor, it can be diagonalized in an orthonormal frame and its components can be more usefully expressed as:

[η ∼ ] =   1 η 12 η 13 η 12 η 22 η 23 η 31 η 32 η 33   = Q   σI / σ11 0 0 0 σII / σ11 0 0 0 σIII / σ11   Q T (14)
where σ I ≥ σ II ≥ σ III are the principal stresses and Q is a rotation matrix defining the loading orientation (that is, the orientation of the eigendirections of η

∼

) with respect to the frame canonically associated to the cube.

From the principal stresses, which are linked by the condition η 11 = 1, two stress invariants can be extracted. The stress triaxiality T and the Lode parameter L are here defined as:

T = tr σ ∼ 3 σvm (15) L = 2 σII -σI -σIII σI -σIII ( 16
)
where σeq is the von Mises equivalent stress. In our convention, the values L = -1, L = 0 and L = 1 respectively correspond to states of generalized tension, shear and compression. Therefore a loading condition can be defined by T , L and a 3D orientation matrix Q.

Conversely, given a triplet (T, L, Q) triplet, the correct tensor η ∼ can easily be found by computing the principal stresses (Zhu et al., 2020a), applying the rotation Q, and normalizing so that η 11 = 1.

Formulation and implementation of the macroscopic Rice criterion

The simulations are primarily aimed at determining the random microstructures' resistance to ductile failure. As the material model defined in subsection 2.2 does not account for damage, failure should be investigated at the cell level. An indicator is required to identify failure. In this work, the localization approach is used by applying Zhu et al.'s (2020a) methodology to the random microstructures. Rice's indicator could also be more suited to complex loading conditions, as used in section 4.

Let L be the macroscopic tangent operator linking the macroscopic Boussinesq stress and deformation gradient rate tensors:

Ṡ ∼ = L : Ḟ ∼ (17)
Rice's criterion predicts strain localization when, for some direction n the acoustic tensor along n becomes singular:

det n k Lik jl n l = 0 (18)

In the following, the expression det n L n will refer to the above equation with the correct choices of indices. Whereas the local tangent operator L is directly given by the material behavior, the macroscopic operator L should be computed at the level of the cell in order to link Ṡ ∼ and Ḟ ∼ . It therefore combines the local material behavior and structural effects. Zhu et al. (2020b) compared several numerical techniques to compute this macroscopic tangent operator. A condensation method derived from Zhu et al. (2020a) but tailored to our finite element formulation is here used. The general idea of the condensation technique is to compute a Schur complement on the finite element matrix separating macroscopic and local degrees of freedom.

More precisely, from subsection 2.3, three types of DOF can be distinguished in the simulation : (i) E i j , macroscopic DOF corresponding to the component i j of F ∼ -1 ∼ ; (ii) v n,i , nodal DOF corresponding to the displacement fluctuation in the direction i of the node n; (iii) four ghost nodes Ê11 , Ê21 Ê31 , Ê32 useful to control the applied stress or the rigid body motion. The vector collecting the E i j (resp Êi j , v n,i ) DOF is noted E (resp. Ê, V ).

The stiffness matrix K computed during the Newton iterations of the FEM simulation, can thus be decomposed in blocks so as to verify:

  K Ê Ê K ÊE 0 K E Ê KEE K EV 0 K V E K VV     ∆ Ê ∆E ∆V   =   ∆R Ê ∆R E ∆R V   (19)
The right hand side represents the external forces on each DOF. R V corresponds to the nodal forces, whereas [R E ] i j = V 0 Si j with S ∼ the macroscopic Boussinesq stress tensor applied to the cell, and V 0 the initial volume of the cell. The blocks K ÊV and K V Ê are zero because the ghost nodes are not linked to the displacement fluctuation DOF. Moreover, KEE is obtained by combining a contribution K EE from the tetrahedral elements of the mesh and contributions from the springs:

KEE,i j,kl = K EE,i j,kl + ∑ s spring K (s) E i j ,E kl (20)
where K (s) E i j ,E kl is the elementary stiffness matrix of the spring. Up to a factor equal to the volume of the cell, K EE represents the average local tangent operator on all tetrahedral elements. For each spring controlling the rigid body motion, only one coefficient in the block EE is non-zero:

K (s) E 21 ,E 21 = k for the E 21 spring, for example. Similarly K (s) E 31 ,E 31 = K (s) E 32 ,E 32 = k.
The expressions for the extension of [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF] special spring element are more complicated but always computable, provided Ê11 and F ∼ are known. The tilde on KEE represents the modification of the quantity K EE linked to the physical mesh by contributions from the springs. The effect of Dirichlet boundary conditions and multi-point constraints on K should now be considered. Dirichlet boundary conditions, which are applied on the ghost nodes, lead to the elimination of DOF. For multi-point constraints applying the periodic boundary conditions of Eq. ( 6), rows related to homologous DOF should be summed and combined in a single new row (and the same for columns). Eliminating the Ê DOF due to Dirichlet boundary conditions and combining rows and columns for homologous V DOF is equivalent to modifying the matrix:

K = KEE K EV K EV K VV ∆E ∆V = K -1 ∆R E ∆R V (21)
The block KEE has not been modified and the DOF Ê have disappeared. Had the conditions E 21 = E 31 = E 32 = 0 not been imposed with springs, these DOF would have disappeared from the stiffness matrix. Computing the related coefficients in the macroscopic tangent operator would have been impossible.

Computing the macroscopic tangent operator means determining the linear relationship between ∆ F ∼ and ∆ S ∼ , and therefore between ∆E and ∆R E , when no other force is applied. In this context, the nodes of the meshed cell are free to move so that ∆R V = 0 at equilibrium. In Eq. ( 19), this implies K V E ∆E + K VV ∆V = 0, hence:

∆R E = K∆E with K = KEE -K EV K -1 VV K V E = ([K -1 ] EE ) -1 (22)
In the last expression, the EE block of K -1 is extracted. The macroscopic tangent operator should then be computed from K by removing the contribution from the springs. Otherwise, even for a homogeneous cube, it would not coincide with the local tangent operator. L is then obtained as:

Li jkl = KE i j E kl /V 0 with K = K EE -K EV K -1 VV K V E (23)
K/V 0 is therefore different from K EE /V 0 , which is the average value of the tangent operator over all the elements of the cell. It includes the structural effect -K EV K -1 VV K V E . In practice, at the end of each increment of the Z-set simulation, a Python script with references to Z-set's internal variables is called. This script can access the global stiffness matrix with boundary conditions K . The quantity [K -1 ] EE can then be computed by solving the following equation for the nine E i j DOF:

K α∆E α∆V = α∆R i j E 0 (24)
where all the coefficients of ∆R i j E are zero, except the coefficient related to E i j : R i j E,i j = 1. The quantity α is a very large value, allowing to neglect the modifications of the right hand side made by the Dirichlet boundary conditions. The matrix [K -1 ] EE is then constructed from the nine solutions by taking for each column, only the coefficients from the E block. Solving these nine equations with Z-set's solver is computationally efficient. No inversion of the matrix K is necessary, and its factorization is only performed once, during the increment.

The matrix [K -1 ] EE is then inverted by means of the Scipy solver. This inversion is inexpensive as [K -1 ] EE is only of size 9 × 9. The resulting matrix K can then be stored in a separate file. This process is repeated for each increment of the simulation. Finally, as a post-processing operation, K is computed from K by subtracting the contribution of the springs.

The macroscopic tangent operator L has now been computed, and Eq. ( 18) remains to solve. Rather than solving it directly, the minimum of det n L n is found for n on the unit sphere. By continuity, if this minimum is nonpositive, Eq. ( 18) has a solution, and the cell has exhibited strain localization. The BFGS minimizer from SciPy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF] is used with the multi-start procedure described by Al [START_REF] Al Kotob | A general and efficient multistart algorithm for the detection of loss of ellipticity in elastoplastic structures[END_REF].

The Appendix A validates the methodology for computing L by comparison with theoretical results and other methods. Note than L only computes the consistent macroscopic tangent operator used in numerical computations, and not directly the instantaneous tangent operator corresponding to the material formulation [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. However the validation from Appendix A proves that the difference can be neglected for sufficiently small steps.

Investigation of localization with Rice's criterion for a diagonal loading

Unit cell studies frequently consider the application of diagonal macroscopic stress tensors. To compare with literature, the analysis is first restricted to diagonal η ∼ loading. The loading is then characterized only by T and L. Rice's criterion's is shown in the following to efficiently detect localization onset and failure mode for random microstructures.

Strain at localization depending on T, L loading conditions

The methodology for Rice's criterion is first verified on the simple case of the unit cell loaded with T = 1, L = -1. This loading condition corresponds to generalized uniaxial tension. For such a loading condition, the failure onset directly corresponds to an inflexion point in the stress-strain curve. It is also well understood by applying [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF] criterion: the cell stops thinning and the macroscopic transverse strain E 22 (or equivalently E 33 ) stabilizes.

Fig. 2 thus compares the evolution of three quantities throughout the computation, parametrized by the macroscopic strain E 11 : the transverse strain E 22 , the macroscopic stress component σ11 along the loading axis, and the Rice quantity min det n L n throughout the computation. For E 11 0.28, min det n L n = 0, which corresponds to Rice's criterion activating. This localization is simultaneous with the stabilization of the transverse strain and the inflection point of the strain-stress curve. The three failure indicators therefore yield the same failure onset, which points to the soundness of the methodology using Rice's criterion. This allows defining the strain at localization E loc as:

E loc = arg min E 11 {E 11 | min det n L n ≤ 0} (25)
The dependence of the strain at localization on T and L is depicted in Fig. 3 for the unit cell and the random cells R1 and R2. Each subfigure considers a slice of the (T, L) space: A parameter is fixed whereas the other one varies. At fixed L = -1 (axisymmetric, generalized tension) or L = 0.5, failure strain decreases almost exponentially with stress triaxiality, for all cells. Such a decrease mostly agrees with the evolution E loc ∝ exp(-3T /2) from [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF], although the agreement, especially regarding R2, is slightly reduced compared to [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF], who used a coalescence indicator.

For fixed T = 1 and varying L (Fig. 3c), the situation is more complex and differs between cells. For the unit cell, three zones can be distinguished. For L < -0.5 or L > 0.3, failure strain evolves slowly with L and the plot forms a plateau. A slow increase of E loc with L is frequently found in unit cell studies with diagonal loading [START_REF] Zhu | Combined effects of triaxiality, Lode parameter and shear stress on void growth and coalescence[END_REF](Zhu et al., , 2020a)). On the contrary, between L = 0.5 and L = -0.3, E loc is a convex function of L, minimal for L = 0 where the localization is almost immediate. An early localization in shear is expected for an homogeneous perfectly plastic material. For the random cells, the central zone of reduced ductility with a minimum for L = 0 is wider than for unit cells. The plateaus at high |L| are reduced and even disappear for the cell R2. For the cell R1, the distinction between the plateaus and the reduced ductility zone are still present. A local maximum associated to a strong slope discontinuity is present at L ∼ 0.5 and marks the transition between the central low ductility zone and the upper plateau. For random cells, the reduced ductility zone is therefore dominant. Such an evolution is reminiscent of the results by [START_REF] Barsoum | Rupture mechanisms in combined tension and shear-Micromechanics[END_REF], [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF] or [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF]; [START_REF] Luo | On the prediction of ductile fracture by void coalescence and strain localization[END_REF], who considered non diagonal loading conditions. The difference with Zhu et al. (2020a) who also used a localization criterion and a diagonal loading can be explained by the present use of a cubic cell (instead of a parallelepipedic cell by Zhu et al.) which is more sensitive to shear localization.

The distinction between plateaus and a central reduced ductility zone was already predicted by a coalescence criterion in [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]. However, the localization criterion used in the present study, predicts a wider reduced ductility zone, especially for the unit cell. In the previous article, that zone was limited to the immediate vicinity of L = 0. Moreover, the plateaus at high |L| are not necessarily present with the localization criterion. A more precise comparison between these failure indicators will be performed in subsection 5.1.

Localization bands and failure mechanisms

In the previous section, Rice's criterion was used as an indicator for failure onset. A complex dependence of E loc with respect to L was found with several failure zones. Cadet et al. ( 2021) correlated these zones with differences in failure modes, identified by their plastic strain fields. However the methodology of Rice's criterion can help quantify and distinguish failure mechanisms.

Let n be a unit vector such that det n L n = 0, and call g the unit eigenvector of n L n with eigenvalue 0. Rice's criterion then predicts a localization band of normal n , and the two sides of the band have a relative velocity along g . Two extremal cases can be distinguished. The case where g and n are collinear represents a perfect extension mode: the band opens without tangential movement. If g and n are orthogonal, the two half spaces delimited by the band have a pure tangential displacement: this represents pure shear. Therefore the analysis of the acoustic tensor n L n at localization allows finding the localization mode within the cell. Furthermore the value of g • n distinguishes extension and shear failure modes. This analysis can be performed not only at localization onset, but more generally throughout the computation. To this end, the vector n minimizing det n L n is found, and the eigenvector for the lowest eigenvalue of the associated acoustic tensor is computed. The solution (n , g ) only describes a localization band at the instant of failure.

Nonetheless, this provides insight on likely bands before localization, and the evolution of such bands after localization. Fig. 4, 5, 7 and 8 depict the characteristics of the localization band for the unit cell and the random cell R1 for two loading conditions (T = 1, L = -1) and (T = 1, L = 0). Each figure shows first the plastic strain rate field at the instant of localization. Plastic strain rate is normalized by the strain rate Ė11 , which is required to drive the computations but plays no role for a rate-independent plastic material. The evolution of g • n and of the principal stress σ 11 is then represented for the whole simulation. Finally the evolution of the det n L n landscape is shown at several time steps during the simulation. The minimizing value is shown with a blue triangle (resp. red dot) before (resp. after) localization. White areas correspond to negative values of det n L n . The unit vector n is parametrized as:

n = (sin(φ ) cos(θ ), sin(φ ) sin(θ ), cos(φ )) (26)
For the unit cell loaded with (T = 1, L = -1), localization takes place near E loc ∼ 0.25. The localization band is a plane orthogonal to the first axis, as shown by the plastic deformation rate field. This is consistent with the evolution of the det n L n landscape. Initially, while the cell is still completely elastic, det n L n does not depend on n . As soon as plasticity begins, a minimum is found near (but not exactly) θ = 0[π], φ = π/2, i.e. n = e 1 . The quantity g • n increases until localization and stabilizes shortly after at a value close to 1. The failure mode is thus extension, which is consistent with a localization band orthogonal to the tensile axis. On the contrary, for (T = 1, L = 0), the localization takes place immediately after yield. Two equivalent localization bands at 45 • from the cube's first axis are present. These bands can also be seen on the det n L n landscapes. At the localization onset, g • n ∼ 0.1: this low value is characteristic of the shear failure mode. However, during the rest of the simulation, g • n continues to grow, and the local minima in the det n L n come closer one to another. They even merge for E 11 ∼ 0.3, which also correspond to a stabilization of g • n at a value close to 1. This can be interpreted as a transition from an initial strain localization in shear mode to an extension mode (Fig. 6). For the random cell R1 with T = 1, L = -1, results are globally similar to the situation of the unit cell for the same loading condition. However the localization band has a more complex pattern, as it tries to connect several voids (already discussed in [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]). The overall orientation of this band remains orthogonal to the first loading axis. For T = 1, L = 0, localization takes place early, in shear mode. Contrary to the case of the unit cell, no stabilization of g • n and no transition to extension mode are witnessed. Moreover, due to the lower symmetry of the random microstructures, the det n L n landscape has lost its symmetry with respect to φ = π/2. The oscillations on Fig. 8b correspond to an alternation between two local minima. The link between localization mode and failure strain is shown in Fig. 9. For all cells, the zones of reduced ductility identified in Fig. 3c are systematically associated to shear failure mode. On the other hand, the peaks and plateaus correspond to extension modes. These results are consistent with those obtained by [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] with a coalescence criterion. However, in the previous study, modes were only distinguished by the inspection of the plastic deformation field. On the contrary, the present study's analysis of the acoustic tensor and its eigenvectors allows to distinguish quantitatively the two localization modes. Moreover this analysis is straightforward, as the acoustic tensor has already been computed in order to apply Rice's criterion.

• n ≥ √ 2/2 (resp. g • n < √ 2/2 ) at localization

Influence of the loading orientation with respect to the cube

In the previous section, Rice's criterion was found to effectively predict failure by localization and was able to distinguish between extension and shear modes. Only simple loading conditions with a diagonal η were considered. The evolution of the strain at localization with respect to L in Fig. 3 showed plateaus and cusps due to the hesitation between failure modes. It is questionable whether this evolution still holds for different loading orientations with respect to the cube. General loading conditions with an arbitrary η ∼ are thus now considered. They allow investigating the influence of the loading orientation. The effects of the intrinsic anisotropy of the simulation cells can then be explored.

Parametrisation of the general loading orientation

As shown in Eq. ( 14), the tensor η ∼ , defining the average Cauchy stress tensor up to a proportionality constant, is equivalent to the combination of T , L, and a rotation matrix Q. This rotation represents the orientation of the loading with respect to the cube. In order to define the computations and report on the results, a suitable parametrisation for Q should be adopted. As the rotation group SO 3 is three-dimensional, Q can be represented by three parameters, for instance by the three components of a Rodrigues vector. If Q corresponds to the rotation of angle θ around the direction given by the unit vector m , with ||m || = 1, the associated vector is:

d = tan(θ /2)m (27)
In order to explore all possible loading orientations, the whole rotation group SO 3 is not necessary due to symmetry conditions. The periodic boundary conditions are naturally associated to cubic symmetry. The unit cell display exactly this symmetry so no generality is lost when reducing loading orientations by cubic symmetry. In order to limit the computational expense, we suppose that it is also the case for random microstructures. Equivalently, the image of a random microstructure by a symmetry of the cube can be thought as a distinct microstructure, whose behavior is not studied in this work. On the other hand, the symmetric tensor η ∼ is characterized by three orthogonal axes bearing distinct principal stresses. It can thus be associated to orthotropic symmetry. Accounting for symmetries, the set of all possible loading orientations with respect to the simulation cell is the cubic-orthotropic disorientation space (using [START_REF] Heinz | Representation of orientation and disorientation data for cubic, hexagonal, tetragonal and orthorhombic crystals[END_REF] terminology). Inversion of a principal axis does not change η ∼ , so the disorientation space is only one half of the one found by [START_REF] Heinz | Representation of orientation and disorientation data for cubic, hexagonal, tetragonal and orthorhombic crystals[END_REF]. The appropriate disorientation space can be expressed in Rodrigues formalism:

0 ≤ d i ≤ √ 2 -1 for i = 1, 2, 3 d 1 + d 2 + d 3 ≤ 1 (28)
A sample of loading orientations in this space is represented in Fig. 10. They are repeatedly used in the following section. Red points correspond to the vertices of the disorientation space, whereas green points were randomly sampled. The latter were first sampled on the unit quaternion sphere with a uniform probability distribution, then were converted to Rodrigues vectors. The representative in the disorientation space is then computed with the method by [START_REF] Grimmer | Disorientations and coincidence rotations for cubic lattices[END_REF].

From T , L and a Rodrigues vector d , the macroscopic stress tensor η ∼ can be computed with the methodology described in subsection 2.3. However the 1-axis plays a special role in the spring element applying boundary conditions (Eq. ( 12)). The computations were found to be more stable when η 11 is the maximal component on the diagonal of η ∼ . The principal axes of η are thus cyclically permuted until this condition is reached. This operation is compatible with cubic symmetry and does not change the equivalence class of loading orientations.

The component E 11 does not take into account rotations of the loading axes, and cannot be used to compare simulations in different loading orientations. By extension of the , associated to its largest eigenvalue, E loc can be defined as:

E loc = F ∼ (t loc ) : (e I ⊗ e I ) (29)
This is equivalent to rotating the deformation gradient in the principal frame of the applied stress, and extracting the new component E 11 :

E loc = Q T • (F ∼ loc -1 ∼ ) • Q 11 (30)
where the notation of Eq. ( 14) is used. In the case of diagonal loading, Q is the identity, and the previous definition of E loc is recovered.

Results for general loading orientations

With the preceding formalism, simulations with any loading orientation can be performed. Before exploring the full orientation space, the simpler case of a rotation with respect to one axis of the cube is studied (Fig. 11). For the condition (T = 1, L = -1), simulations are performed on the unit cell and two random cells. The applied stress is rotated with various angles around the second or third axes of the cubes. In the Rodrigues formalism, d is collinear with e 2 or e 3 . As L = -1 is an axisymmetric case, rotation around e 1 has no effect on the applied stress. For the unit cells, results for the second and third axes are identical, due to symmetry. Failure strain decreases as the angle of rotation increases, until a minimum is reached at approximately 35 • . Notwithstanding the angular discretization, this is fully consistent with the results by [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], who found a minimum for an angle of 37 • . On the other hand, the two random microstructures show a more complex dependence of the failure strain with respect to the angle of rotation. The symmetry between the axes 2 and 3 is lost with the random cell. Moreover, results from R1 and R2 are quite different one from another. Minimum failure strain is not necessarily found for an angle of 37 • (10 • for R1, for instance). A clear peak is found for R1 and R2 for a rotation of 40 • around e 2 . The analysis of plastic deformation fields (not presented here) shows a transition of modes between the two sides of the peak. This phenomenon is similar to the cusps in Fig. 3c. In order to explore the effect of orientation for given values of T and L, 150 different rotations were chosen in the disorientation space. Failure strain values were computed for these rotation points on the unit cell and the random cells at (T = 1, L = -1). Results are displayed in the disorientation space (Fig. 12). The case (T = 1, L = 0.5) is also shown for the unit cell. In the rest of this section, a coarser mesh with 7682 nodes (the same as in [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]) was chosen for the unit cell in order to limit the computational expense. For all cases, failure strain strongly depends on the loading orientation. Note that the minimal failure strain is systematically lower than in Fig. 11. In order to find this minimal failure strain, restricting to rotations around the axes is therefore not justified. Values of failure strain appear to depend in a rather smooth way on the loading orientation: Peaks corresponding to hesitation between modes as in Fig. 11a are scarce. Very clear trends can be seen for the unit cell. For instance, at (T = 1, L = 0.5), E loc appears to mainly depend on the value of d 2 . However results from different cells or different values of L do not display the same trends with respect to the loading orientation. After having understood the structure of the dependence with respect to Q for a given value of L, we consider the loading condition as mainly a function of (T, L) and study the scatter which can be attributed to the effect of the loading orientation. Computations were performed on the three cells for various values of L, at fixed T = 1. For each value of L, 30 simulations with different loading orientations (visible on Fig. 10) were performed. Results are shown in Fig. 13 as box plots, representing the minimum, first, second (median), third quartiles and maximum for the 30 different orientations. For each cell, loading orientation is responsible for a large scatter. An order of magnitude can sometimes be found between the minimal and maximum failure strain for a given L. Results from the unit cell are even more scattered, which is the sign of greater anisotropy for the localization behavior. For L = 0, this scatter is however minimal, because the behavior of the material imposes early localization in shear. For random cells, minimal failure strain is mostly a symmetric function of L, with minimal value for L = 0. The curve of the minimal failure strain is however flatter for the unit cells, due to the higher anisotropy. Counter-intuitively, localization is thus found to take place earlier in the unit cell than in the random cells, because the former is more sensitive to the loading orientation.

-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 L 10 2 The localization mode can also be investigated depending on the loading orientation. Fig. 14 shows the plastic deformation rate fields just after localization for the random cells R1 (Fig. 14a to 14d) and R2 (Fig. 14e to 14h). Two loading cases are considered: (T = 1, L = -1) (first and third rows) and (T = 1, L = 0.4) (second and fourth rows). For each loading case, the results for the most (resp. least) resistant orientations associated to the maximal (resp. minimal) failure strain value are shown on the left (resp right). The failure mode for the least resistant orientation is always associated to a shear failure mode. This is even true for the cell R2 at L = -1, for which the localization band is parallel to the sides of the cube. For the most resistant, shear or extension modes are possible. This shows that L is not sufficient to distinguish between the two failure modes. For a given L, the failure mode can change depending on the loading orientation.

The link between early localization and shear failure mode can be quantified with g • n , as in subsection 3.2. Fig. 15 correlates g • n at localization and E loc for all the simulations on random cell R2, seen in Fig. 13b. E loc is almost a linear function of g • n . A small failure strain is thus systematically associated to shear failure mode, whereas late localization is only found with extension mode.

The results of this section show the strong influence of the loading orientation. Yet the response of an isotropic material should only depend on the stress invariants T and L. The dependence with the loading orientation is therefore the sign of the anisotropic behavior of the cubic simulation cells. This anisotropy can not be avoided but could be reduced with larger simulation cells. The random microstructures already display a smaller anisotropy than the unit cell. As proposed by [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF] and in order to obtain conservative results, the application of simulations to isotropic models should focus on the critical loading orientation with the minimal failure strain. This critical orientation may vary with T and L, but is always associated to a shear failure mode.

Discussion

The previous results are now discussed with help from complementary simulations. Rice's criterion is first compared to other criteria from literature. Then the link between cell anisotropy and the scatter of results is further analyzed. The anisotropy in yield strength is studied, and the influence of a reduced porosity on localization is investigated.

Comparison of failure criteria for diagonal loading

In this study, Rice's localization criterion was used to detect the failure of the cells. Its prediction in the case of diagonal loading are here compared to those of two coalescence criteria (Fig. 16). The energy criterion from [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] compares the macroscopic plastic power Ẇp and elastic power Ẇe . In the coalescence state, elastic unloading et al.'s (2021) criterion on the determinant of the deformation gradient rate, and [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] energy criterion. All computations with a diagonal loading, at T = 1. For the random cell, the data points from the δ and energy criteria are directly taken from [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF].

occurs far from the voids, so the ratio Ẇe / Ẇp is minimal. The δ criterion from [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] assumes that for shear or extension coalescence modes, the macroscopic deformation gradient rate Ḟ ∼ becomes singular (rank one for extension, and rank two for shear). The evolution of det Ḟ ∼ can then be compared to that of a homogeneous loading. Coalescence is thus detected when this ratio reaches zero (within a given precision). More precise definition and implementation of these two criteria can be found in [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF].

For both the unit and random cells, localization (from Rice's criterion) takes place earlier than or almost at the same time as coalescence (from the two other criteria). This is consistent with the findings by [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] and Zhu et al. (2020a). For the unit cell, the three criteria give very similar predictions for L < -0.5 and L > 0.3. However, in the region L ∈ [-0.5, 0.3], localization takes place very early, whereas failure strain from the energy criterion remains high and grows with L. The behavior of the δ criterion is intermediate, with results close to that of the energy criterion, except in the immediate vicinity of L = 0, where failure strain is minimal. This division of the E loc -L curve in three zones matches the results from Fig. 9. The three criteria give similar results when the unit cell fails in extension, and quite different results for the shear failure mode. The case L = 0 can be better understood from Fig. 6. There is early strain localization in shear, but the failure mode ends switching to extension mode, which is recognized by the energy criterion. For the random cell, at low and high L, coalescence criteria give the same results. Coalescence and localization are almost simultaneous at high L, whereas localization occurs significantly earlier at low L (for generalized tension). In the region L ∈ [-0.5, 0.4], the energy criterion fails to detect any coalescence (effect already seen in [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]). The failure strain from the Rice and δ criteria both form a U-shaped curve. The curve is wider for localization, as localization takes place earlier than coalescence. In this region, failure takes place in shear mode and no switch to extension mode was seen in Fig. 8. To summarize, localization always precedes coalescence. The three criteria yield similar results for extension modes. For shear modes, the energy criterion fails to detect coalescence whereas the other criteria can recognize this mode. The zone of low ductility is wider for the localization criterion.

Anisotropy of the cells

The results of subsection 4.2 showed a major dependence of the failure strain on the loading orientation with respect to the cell axes. This can be interpreted as an evidence of an anisotropic localization behavior of the simulation cells. In order to precise this anisotropy interpretation, we show that the dependence on loading orientation does not only concern localization strain, but also the initial yield stress. The simulations from Fig. 13 are analyzed again, and the maximum value of the principal stress σI is extracted. Values of σI can be directly compared between simulations at different loading orientations. Fig. 17 shows the scattering of stress values with respect to the orientation for the unit cells and the random cells R1 and R2. This scattering is significant but is not as marked as for failure strain values. It confirms the anisotropic behavior of the simulation cells. Compared to unit cells, anisotropy is significantly reduced for the random microstructures, in agreement with Fig. 13. However, the scatter of maximum stress is maximal near L = 0, whereas the scatter of localization strain was minimal there. The anisotropic behavior, both for localization and for initial yield stress, is due to the geometric anisotropy of the cells. With the periodic boundary conditions, the material is a cubic array of voids and therefore displays cubic anisotropy. The anisotropy is smaller for random microstructures because they are associated to larger volume elements, or equivalently because the period of the cubic array created by boundary conditions is larger (relatively to the size of the voids).

Cubic unit cells thus exhibit major anisotropy. In order to determine an isotropic response from cell simulations, computations should be performed for many loading orientations. For random microstructures, this anisotropy is smaller but remains significant. With larger microstructures containing more voids, the anisotropy could probably be further reduced, at the expense of a higher computational cost. Another method would be to use different boundary conditions. Such a modification was proposed by [START_REF] Coenen | Novel boundary conditions for strain localization analyses in microstructural volume elements[END_REF] to counteract the influence of the periodic boundary conditions, which are here responsible for the anisotropy. The most likely localization band is detected from the strain field. The boundary conditions are then adjusted to allow the development of this band (and only this one). Yet this approach does not seem adapted to the complex localization patterns seen in Fig. 14 or to possible transitions of failure modes (Fig. 6).

Effect of porosity

A single fixed porosity of 6% was considered up to now. In order to assess the influence of the porosity, a unit cell of porosity 1% was generated with the same procedure as previously. The only difference resides in the radius of the embedded spherical void. This cell was subjected to the same loading conditions as in Fig. 13 and the results is presented in Fig. 18. When compared to the 6%-porosity unit cell, a significant influence of the loading orientation is still present. The minimal failure strain still draws a U-shaped curve, with a minimum near L = 0. However the values of failure strain are generally higher. This is logical as lower porosity induces higher resistance to ductile fracture. More interestingly, -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 L 10 2

10 1 the relative scatter with orientation is reduced, and the curve of minimal failure strain is less flat than for the 6%-porosity cell. As the void is less susceptible to interact with the periodic boundary conditions, the anisotropy is reduced for low volume fractions.

E loc

Conclusion

Random microstructures consisting in a perfectly plastic matrix containing a random distribution of voids were simulated up to failure. A wide range of loading conditions, characterized by the stress triaxiality and the Lode parameter, was studied. The orientation of the principal directions of the applied stress with respect to the cubic cell's axes was systematically investigated. Failure was detected with a strain localization indicator: Rice's criterion computed at the macroscopic level. The major findings of this study follow:

• Rice's criterion successfully detects the onset of strain localization, the orientation of the localization band. Moreover it allows distinguishing two types of failure modes: shear and extension modes.

• For diagonal loading, in which the directions of the applied stress coincide with those of the cubic cell, a central zone of reduced ductility is found both for the unit cell and the random cells. Shear failure mode is systematically associated to reduced ductility.

• The loading orientation strongly affects failure strain. For given stress invariants, the values of failure strain in the least and the most resistant orientation can differ by A more complex case is then considered with an elastic porous unit cell. Strain, and therefore the local tangent operator, is no more homogeneous within the cell. The reference results are now computed with a perturbation method (Zhu et al., 2020b). Regularly during the main simulation, the state of the finite element problem is saved. Secondary computations are then launched from the saved states by applying a small perturbation ∆ F ∼ = e i ⊗ e j , with i, j ∈ [ [1,3]]. The evolution of the Boussinesq stress tensor ∆ S ∼ is then observed. The macroscopic tangent operator can then be estimated by Li jkl = ∆ Si j /∆ Fkl . This method is less practical than the condensation method because it requires launching nine secondary computations for each output time of the main computation. More problematically, it is not adapted to multibranch tangent operators, which is automatically the case in elastoplasticity. Perturbations ∆ F ∼ = ±e i ⊗ e j may correspond to either plastic loading or elastic unloading. Distinguishing which one corresponds to the correct evolution in the primary computation is not easily possible a priori. This problem does not affect pure elasticity. For (T = 1, L = -1) on the elastic unit cell, the condensation and perturbation methods give identical results ( Fig. A.20). However only taking the average of the local operators is not enough, for it implies neglecting structural effects.

A last validation is performed by comparison with a theoretical result. The homogeneous cube from the first validation is filled with with a softening material: R soft (p) = R 0 (1exp(500p) (with otherwise the same characteristics as in subsection 2.2). For a von Mises material, according to [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF], the localization should take place when dR/d p = -E/4. This is effectively the case, as shown by This section studies the influence of the temporal and spatial discretization on the localization results (Fig. B.22). Simulations were carried out with diagonal loading conditions, at T = 1 and various values of L. The fine mesh from Fig. 1c and the coarse mesh from [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF] were compared. The effect of fine time steps was also studied. As the simulation is driven by the evolution of E 11 , time steps correspond to the increase ∆E 11 between two iterations of the implicit Euler method. In the main study, increasing time steps from ∆E 11 = 0.0001 to ∆E 11 = 0.005 were used. This allowed a sufficient precision even in the case of early localization. Constant coarse (∆E 11 = 0.005) and fine time steps (∆E 11 = 0.001) are here tried out. For the coarsely and finely meshed unit cells, the overall appearance of the curves are similar with a plateau at low L and a zone of reduced ductility around L = 0. The behavior at high L > 0.3 is however different with a clearer plateau for the fine mesh. This justifies preferring the finely meshed cell in this study. Nevertheless, the unit mesh was used for Fig. 12 and 13 in order to limit the computational cost. For the fine cell, there is almost no effect of the time step (except at the lower end of the central reduced zone). The effect is more pronounced for the coarse cell, especially in the L > 0.3 region. The finer the mesh and the time steps, the narrower the zone of reduced ductility.

∼

  dev the deviatoric part of the rotated Cauchy stress tensor s ∼ , s vm the equivalent von Mises stress and ṗ = 2 3 ė ∼ p : ė ∼ p playing the role of the plastic multiplier. The Young modulus, the Poisson ratio and the yield strength are respectively chosen as E = 200 GPa, ν = 0.3 and R 0 = 500 MPa.

Figure 1 :

 1 Figure 1: Meshes of the cells repeatedly used in the study

Figure 2 :

 2 Figure 2: Evolution of transverse strain, stress component σ11 and Rice criterion for a simulation of a unit cell at T = 1, L = -1. The stabilization of the transverse strain and the localization onset are almost simultaneous. This instant also corresponds to an inflection point of the stress-strain curve.

Figure 3 :

 3 Figure 3: Deformation at localization for unit and random cells (diagonal loading) and slices of the (T, L) space. The triangles indicate a slope of -3/2 in the log-lin plot: E loc ∝ exp(-3T /2)

Figure 6 :

 6 Figure 4: Characteristics of the localization band for the unit cell, at (T = 1, L = -1)

  Figure 7: Characteristics of the localization band for the random cell R1, at (T = 1, L = -1

  Figure 8: Characteristics of the localization band for the random cell R1, at (T = 1, L = 0)

Figure 9 :

 9 Figure9: Link between the deformation at localization and the localization mode (T = 1). Same as Fig.3c, but with localization modes highlighted. The extension mode (resp. shear mode) is shown in red (resp. green) and corresponds to g• n ≥ √ 2/2 (resp. g • n < √ 2/2 ) at localization

Figure 10 :

 10 Figure 10: Disorientation space representing all loading orientations with respect to the cube, along with simulation points frequently used in this study

Figure 11 :

 11 Figure 11: Dependence of the deformation at localization for different rotations around the axes (three cells, T = 1, L = -1)

Figure 12 :

 12 Figure 12: Deformation at localization for 150 different loading orientations, represented in the disorientation space

Figure 13 :

 13 Figure 13: Localization onsets for three cells (T = 1, varying L). Each box plot corresponds to the results of 30 different loading orientations.

Figure 15 :Figure 16 :

 1516 Figure14: ṗ/ Ė11 just after localization, for the orientations maximizing (left) and minimizing (right) E loc , for the two random cells and two loading conditions (rows). To better visualize the failure mode, the deformed mesh is represented, and a factor is applied to the displacements.

Figure 17 :

 17 Figure 17: Anisotropy of the maximum principal stress, for the unit cell and two random ones. Computations at T = 1, varying L, for 30 different loading orientations.

Figure 18 :

 18 Figure 18: Strain at localization for T = 1, varying L and 30 different loading orientations, on a unit cell with porosity 1%

  Figure A.19: Validation of the computation of macroscopic tangent operator on homogeneous volume elements

Fig

  Figure A.20: Validation of the computation of macroscopic tangent operator on an elastic unit cell

Figure A. 21 :

 21 Figure A.21: Verification of the -E/4 theoretical localization result on a homogeneous cube: for a von Mises material with softening R(p), the localization takes place when dR/dp = -E/4

Figure B. 22 :

 22 Figure B.22: Sensitivity of the determination of localization by Rice criterion to the mesh size and to the time step size

(a) ṗ/ Ė11 after localization

more than one order of magnitude. The failure mode may depend on the loading orientation.

• The sensitivity to loading orientation is due to the intrinsic anisotropy of the cell.

Although it affects both random and unit cells, the scatter is reduced for random microstructures, as they are less anisotropic.

• Minimal failure strain on all loading orientations was obtained both for unit cells and random populations of voids. It is a U-shaped function of L, minimal for L = 0. In order to reach this minimum, considering only a rotation along a fixed axis of the cube, as commonly done in literature, is not sufficient and general 3D rotations are required.

In order to develop reliable models of ductile fracture, minimal failure strain should be sought among all loading orientations. Moreover, random microstructures should be favored over unit cells thanks to their lower anisotropy in localization. However the simulation of random microstructures is significantly more computationally expensive than unit cells. An extension of this work should be to obtain a surrogate model for the ductile fracture of random cells. This model should require as small a number of simulations as possible by leveraging the link between random microstructures and unit cells.

Appendix A. Validation of the computation of macroscopic tangent operator

The method for computing the macroscopic tangent operator is here validated by comparison with other methods. Homogeneous cubes of 3 × 3 × 3 elements without voids are first considered. A purely elastic material and an elastoplastic material are considered. The formulation of Eq. ( 3) is used, but R 0 is replaced either by R elas = +∞ or R hardening (p) = R 0 + H p (with H = 2000 MPa). Periodic boundary conditions and a large strain formulation are used. A diagonal loading condition corresponding to (T = 1, L = -1) is applied. As the cubes are homogeneous, the macroscopic tangent operator is equal to the local tangent operator, which is directly computed by the software [START_REF] Zset | Zset[END_REF] within the element. Fig . A.19 compares the evolution of the Rice criterion using macroscopic tangent operators obtained from three methods: (i) the condensation method of subsection 2.4; (ii) the local operator within an element; (iii) taking K EE (in Eq. ( 23)) as a macroscopic tangent operator. This last method is equivalent to averaging the local tangent operator over all elements, and neglecting their interactions. The three methods give identical results for the two types of materials, as should be. However, the averaging method only works here because of the homogeneous cube. For the elastic material, the Rice criterion is not constant because of the large strain formulation.