Modesar Shakoor
email: modesar.shakoor@mines-paristech.fr

Marc Bernacki

Pierre-Olivier Bouchard

Ductile fracture of a metal matrix composite studied using 3D numerical modeling of void nucleation and coalescence

Keywords: ductile fracture, void nucleation, void growth, void coalescence, metal matrix composites

come

Ductile fracture of a metal matrix composite studied using 3D numerical modeling of void nucleation and coalescence

Introduction

Among recent efforts in materials science, the development of more predictive and general constitutive models based on micromechanical criteria has gained an increasing interest. Thanks to progresses in experimental means, the in situ observation of ductile fracture's void nucleation, growth and coalescence phenomena in full 3D is now possible [START_REF] Hannard | Characterization and micromechanical modelling of microstructural heterogeneity effects on ductile fracture of 6xxx aluminium alloys[END_REF][START_REF] Hosokawa | Void growth and coalescence in model materials investigated by high-resolution X-ray microtomography[END_REF][START_REF] Kaye | A new parameter for modelling three-dimensional damage evolution validated by synchrotron tomography[END_REF][START_REF] Morgeneyer | In situ 3-D observation of early strain localization during failure of thin Al alloy (2198) sheet[END_REF][START_REF] Seo | Three-Dimensional Investigation of Void Coalescence in Free-Cutting Steel using X-ray Tomography[END_REF]. Incorporating these observed micromechanisms into constitutive models requires a thorough analytical and/or numerical analysis and comparisons with experiments [START_REF] Benzerga | Ductile Fracture by Void Growth to Coalescence[END_REF][START_REF] Pineau | Failure of metals I: Brittle and ductile fracture[END_REF]. However, both analytical and computational models do not yet address conditions that are fully comparable to what is observed in experiments.

Very restrictive assumptions of a single spherical void embedded in a perfectly plastic von Mises matrix under specific loading conditions were used by Gurson in his well-known fundamental work [START_REF] Gurson | Plastic Flow and Fracture Behavior of Ductile Materials Incorporating Void Nucleation, Growth, and Interaction[END_REF]. Many analytical developments have generalized Gurson's yield criterion to an ellipsoidal void shape and more general loading conditions [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-Case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Castañeda | Constitutive models for porous materials with evolving microstructure[END_REF][START_REF] Nahshon | Modification of the Gurson Model for shear failure[END_REF][START_REF] Torki | Theoretical and numerical analysis of void coalescence in porous ductile solids under arbitrary loadings[END_REF], and to complex matrix behavior [START_REF] Benzerga | Plastic potentials for anisotropic porous solids[END_REF][START_REF] Besson | An extension of the Green and Gurson models to kinematic hardening[END_REF].

Similarly, some FE analysis studies have considered microstructures consisting of a single void, or a periodic array of voids [START_REF] Borg | Size effects on void growth in single crystals with distributed voids[END_REF][START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF][START_REF] Dahl | Effect of Contact Conditions on Void Coalescence at Low Stress Triaxiality Shearing[END_REF][START_REF] Fabrègue | A constitutive model for elastoplastic solids containing primary and secondary voids[END_REF]. However, metal alloys feature particles and considering them as voids means neglecting their effect on the load carrying capacity, and the influence of void nucleation mechanisms.

The influence of the particles and their debonding from the matrix has been considered in Refs. [START_REF] Charles | Modelling the competition between interface debonding and particle fracture using a plastic strain dependent cohesive zone[END_REF][START_REF] Ganguly | Rearrangement of local stress and strain fields due to damage initiation in a model composite system[END_REF][START_REF] Giang | Influence of carbide particles on crack initiation and propagation with competing ductile-brittle transition in ferritic steel[END_REF], with a restriction to periodic arrays of particles.

The competition between particle debonding and fragmentation for these periodic arrangements is neglected in most studies, with the exception of Ref. [START_REF] Giang | Influence of carbide particles on crack initiation and propagation with competing ductile-brittle transition in ferritic steel[END_REF].

Idealized microstructures such as periodic arrays of voids or particles are not representative of real ductile materials, which feature random arrangements of voids and particles. Random arrangements of particles with a debonding model have been addressed in Ref. [START_REF] Mcveigh | An interactive micro-void shear localization mechanism in high strength steels[END_REF], while a fragmentation model has been used in Ref. [START_REF] Eckschlager | A unit cell model for brittle fracture of particles embedded in a ductile matrix[END_REF]. Random arrangements of voids have also been considered in the literature [START_REF] Bandstra | On the influence of void clusters on void growth and coalescence during ductile fracture[END_REF][START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF].

In spite of this substantial literature, a number of assumptions still have to be investigated. The number of particles modeled in most 3D FE analysis studies is generally not large enough to statistically represent the random distribution of particles found in real materials [START_REF] Steven Greene | Computational uncertainty analysis in multiresolution materials via stochastic constitutive theory[END_REF][START_REF] Matouš | A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials[END_REF]. The competition between the two void nucleation mechanisms of particle debonding and fragmentation, and their interaction with void coalescence at large plastic strains have not been investigated yet for large and random 3D arrangements [START_REF] Pineau | Failure of metals I: Brittle and ductile fracture[END_REF]. There is definitely an interest for a parallel numerical framework enabling for simulations of large 3D microstructures with a representative number of randomly distributed particles [START_REF] Matouš | A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials[END_REF], as well as for numerical methods that can be used to model particle debonding and fragmentation, and void coalescence [START_REF] Giang | Influence of carbide particles on crack initiation and propagation with competing ductile-brittle transition in ferritic steel[END_REF]. An attempt to develop such framework has been proposed in a previous work by the authors, and successfully applied to the 2D full field modeling of a representative number of randomly distributed voids [START_REF] Roux | A level-set and anisotropic adaptive remeshing strategy for the modeling of void growth under large plastic strain[END_REF][START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF] and particles [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF]. Particle debonding and fragmentation criteria were introduced [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF], as well as a void coalescence criterion [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF]. These developments nevertheless suffered from several limitations that did not enable for large scale 3D simulations.

In this paper, these limitations are addressed thanks to recent numerical developments. Then, a FE micromechanical analysis is proposed that accounts for both particle debonding and fragmentation, and the growth and coalescence of the nucleated voids, at large plastic strains and for large 3D microstructures.

The main objective of this analysis is to show the relevance of a simultaneous modeling of both void nucleation and coalescence, and to reveal interactions between these two phenomena, with a major influence on ductility.

The existing numerical framework is described in Section 2, while new methods and micromechanical models introduced for the purpose of the present analysis are presented in Section 3. This improved FE analysis tool is then applied to show the influence of void nucleation on localization and void coalescence in Section 4.

Numerical framework

Constitutive modeling

The domain is a 2D or 3D Representative Volume Element (RVE) generally composed of a matrix, particles, and voids. In the present framework all phases are meshed, with the use of a conform FE meshing strategy at boundaries between different phases. A particularity of the present framework is that voids are also meshed, meaning that there are mesh elements inside cracks. These elements are only used for remeshing purposes, and measurement of void volume fraction, but have no stiffness or constitutive behavior. Their effect is hence no different than actual voids (that would not be meshed).

Hooke's law is used to model elasticity in both matrix and particles. In addition, a von Mises yield criterion is used for the matrix, with power hardening law [START_REF] Ludwik | Elemente der technologischen Mechanik[END_REF]

σ 0 (ε) = σ y + K (ε) n (1)
where ε is the equivalent (von Mises) plastic strain, σ y the yield stress, K the plastic consistency and n the hardening exponent. Particles are considered elastic brittle and hence assumed to debond or fragment before yielding can occur [START_REF] Lee | Stress concentration induced by an elastic spheroidal particle in a plastically deforming solid[END_REF][START_REF] Eckschlager | A unit cell model for brittle fracture of particles embedded in a ductile matrix[END_REF][START_REF] Ganguly | Rearrangement of local stress and strain fields due to damage initiation in a model composite system[END_REF][START_REF] Mcveigh | An interactive micro-void shear localization mechanism in high strength steels[END_REF]. Fracture criteria are detailed in Section 3.

To reach large strains, an updated Lagrangian formulation is used where a static load is applied progressively using small increments. Each increment consists in solving continuum mechanics equations and then applying the computed displacements to move mesh nodes. Because of incompressible plastic behavior in the matrix, and the use of triangular (2D) and tetrahedral (3D) elements, a mixed velocity-pressure formulation solved using a P 1 + /P 1 element is adopted to avoid locking [START_REF] Boffi | Mixed Finite Elements, Compatibility Conditions, and Applications[END_REF]. The nonlinearity of this formulation, due to the plastic behavior of the matrix, is solved implicitly using a Newton-Raphson scheme [START_REF] Wagoner | Metal Forming Analysis[END_REF].

Mesh motion and adaption

The key aspects of the present framework are linked to remeshing and the mesh motion step that is applied at the end of each load increment. These algorithms enable for simulations featuring large void growth ratios and complex topological events such as void coalescence. The mesh motion algorithm first tries to move mesh nodes to their final positions by applying the computed displacement directly. If element flipping is detected (by computing element volumes), nodes are moved back to intermediary positions by dividing the displacement by two. This subdivision is applied recursively until the displacement is small enough so that no flipping is detected [START_REF] Shakoor | An adaptive level-set method with enhanced volume conservation for simulations in multiphase domains[END_REF]. Remeshing is then operated with the expectation that the new mesh will render a motion towards final positions possible. If not, the same procedure is applied again until final positions are reached.

Both this mesh motion and the mesh adaption algorithm used at each motion iteration were introduced in previous works [START_REF] Shakoor | An efficient and parallel level set reinitialization method -Application to micromechanics and microstructural evolutions[END_REF][START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF][START_REF] Shakoor | An adaptive level-set method with enhanced volume conservation for simulations in multiphase domains[END_REF], and are hence only summarized herein. The main properties of the mesh adaption algorithm are that:

• It will try to maximize the volume of each element to reduce the risk of element flipping, while also accounting for a mesh size criterion [START_REF] Gruau | Coupez, 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF], and preserving phase boundaries [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF]. A volume conservation constraint is also added to preserve at best the volume of each phase, though this constraint can be relaxed to eliminate degenerated elements of very low volume [START_REF] Shakoor | An adaptive level-set method with enhanced volume conservation for simulations in multiphase domains[END_REF].

• The mesh size criterion is isotropic but varies locally based on the local maximum principal curvature at phase boundaries. The aim is to capture fine features (e.g., small particles, crack tips) as well as coarse ones [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF]. The maximum principal curvature of a boundary can be shown to be linked to the second derivatives of the distance functions to that boundary [START_REF] Quan | Anisotropic mesh adaptation with optimal convergence for finite elements using embedded geometries[END_REF]. These derivatives are computed using an efficient and parallel distance computation algorithm [START_REF] Shakoor | An efficient and parallel level set reinitialization method -Application to micromechanics and microstructural evolutions[END_REF], then followed by a parallel second order superconvergent patch recovery technique [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF].

Lastly, fracture indicators presented in Section 3 require the computation of element-wise tags to distinguish between different particles and particle fragments. This is not a trivial operation since the FE mesh is distributed among multiple computing units, and both mesh and particles distributions evolve due to remeshing and particle fragmentation. This tracking of particle tags is performed using an automatic and parallel connected components identification algorithm developed in a previous work [START_REF] Scholtes | New finite element developments for the full field modeling of microstructural evolutions using the level-set method[END_REF].

3. Improvements for full field 3D simulations of void nucleation 3.1. Void nucleation modeling 3.1.1. Crack meshing method Thanks to the numerical framework detailed in Section 2, RVE simulations can be performed. In a previous work [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF], a region deletion technique was developed where quite large cracks with rounded crack tips were introduced to model crack initiation and propagation. In this section, a novel crack modeling method is proposed to capture fine cracks with sharp crack tips. It has the added advantage that mass loss issues linked to the region deletion technique are drastically reduced.

The proposed method relies on Level-Set (LS) functions and automatic interface meshing. The LS function associated to a closed surface Γ s = ∂Ω s is defined as the distance field to ∂Ω s , multiplied by -1 for points that are out of Ω s [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. The use of LS functions for crack modeling is inspired from the X-FEM framework [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF]. The crack surface, which is generally an open surface, is linearly extended to the whole computational domain and an LS function φ c is associated to this extended closed surface Γ c . This LS function φ c is positive on one side of Γ c , and negative on the other side. Another surface Γ f delimits the region of the domain Ω f occupied by the actual crack (generally an open surface), so that its associated LS function φ f is positive inside Ω f , and negative outside. The automatic interface meshing algorithm consists in adding to the FE mesh all intersections between Γ c and mesh edges that have a positive φ f value on at least one of their ends, as well as all intersections between Γ c and Γ f . This algorithm is an extension of the interface fitting algorithm introduced in a previous work [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF].

To avoid numerical issues in the FE code due to superimposed nodes, cracks are initially modeled with a small opening c , which is at least one order of magnitude lower than mesh size1 . Thus, the crack surface Γ c is in the implementation defined by two crack faces Γ + c and Γ - c that are not initially perfectly superimposed. By convention, the LS function associated to Γ + c is φ + c = φ c , and the LS function associated to Γ - c is φ - c = -(φ c + c χ f ), where χ f (x) = 1 if φ f (x) > 0, 0 otherwise. This definition ensures that the two crack faces join on Γ f in order to form the crack tip.

The main advantage of the proposed approach over the X-FEM is that cracks are explicitly meshed; hence unstructured mesh adaption can be used to track the deformation of crack faces and to model large void growth ratios and void coalescence. This computational fracture mechanics approach is illustrated in the following in the case of particle debonding and fragmentation.

(a) (b) ϵ c ϵ c Γ f Γ c Γ c Γ f l D ϵ f void/crack matrix particle + - Γ c + Γ c - Figure 1:
Definitions of the crack faces LS functions φ + c and φ - c , the filter LS function φ f , and numerical parameters for void nucleation modeling by: (a) particle debonding, (b) particle fragmentation. Both c and f are exaggerated on this figure for sake of illustration.

Particle debonding

A simple critical interface stress criterion is used to model particle debonding. While the interface stress was defined as a linear combination of normal and tangential stresses in a previous work, the effect of tangential stresses is not considered herein due to an absence of experimental data. Particle debonding is hence predicted when the local normal stress σ n at any mesh face of the particles/matrix interface reaches a critical value σ D c . To avoid any mesh size dependency, crack initiation is not restricted solely to that mesh face but to a neighborhood in a radius D . The crack LS function φ c = φ + c to the crack surface is hence defined at mesh nodes as negative inside the concerned particle, and positive inside the matrix. The filter LS function φ f is defined as positive inside the sphere of radius D centered on the concerned mesh face, and negative outside this sphere. This definition of the crack geometry is illustrated in Fig. 1(a). Elements with positive φ + c , φ - c = -(φ c + c χ f ), and φ f are added to the void phase, and thus removed from the matrix phase. It is reminded that due to the very low value of c (exagerated in Fig. 1(a)), this loss of mass is negligible, as opposed to previous work [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF].

Particle fragmentation

The particle fragmentation criterion requires a more thorough analysis. First, particle fragmentation is assumed as instantaneous, as crack propagation inside particles of ≈ 10 µm size is a phenomenon that would require an analysis at a spatial scale finer than the microscale, and an extremely fine time scale. Second, an instantaneous fragmentation criterion requires a fracture indicator representative of the whole particle. In a previous work [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF], a local critical first principal stress σ 1 criterion was used. The crack surface was defined as orthogonal to the first principal direction. Because σ 1 often reached large values close to the interface, this criterion led to non physical fragmentation cracks normal to the interface [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF]. To avoid this issue, a non local criterion is proposed [START_REF] Moorthy | A Voronoi Cell finite element model for particle cracking in elastic-plastic composite materials[END_REF]. The first principal stress is averaged on the whole particle, and the first principal direction is taken at the σ 1 -weighted center G p of the particle Ω p

G p = 1 Ωp σ + 1 (x)dx Ωp σ + 1 (x)xdx. (2)
The positive part σ + 1 (x) is defined as σ + 1 (x) = σ 1 (x)χ(σ 1 (x)). In the present framework with linear FEs, σ 1 is computed by diagonalization of the Cauchy stress tensor σ, which is defined element-wise. The coordinates x in Eq. ( 2) are hence the coordinates of the barycenter of each element. The characteristic function χ(σ 1 (x)) filters out elements with zero or negative first principal stress, which represent parts of the particle that are under compression and not carry any tensile load. Similarly, the particle fragmentation indicator σ 1 is defined as

σ 1 = 1 Ωp χ(σ 1 (x))dx Ωp σ + 1 (x)dx. (3)
Fragmentation is predicted if σ 1 reaches a critical value σ F c and G p belongs to an element of Ω p with positive σ 1 . The first principal direction at that element is chosen as the normal vector to the crack surface, which is defined as a plane. The crack LS function φ c = φ + c is defined as positive on one side of this plane, and negative on the other side. The filter LS function φ f is defined as positive inside the particle, and negative outside. To ensure that all mesh elements and nodes of the two particle fragments are effectively separated, this filter LS function is extended with an artificial thickness f . This definition of the crack geometry is illustrated in Fig. 1(b). Elements with positive φ + c , φ - c = -(φ c + c χ f ), and φ f are added to the void phase, and thus removed from the particle. Because of the artificial fragmentation crack extension length f , some elements are removed from the matrix phase. A similar criterion with artificial crack extension can be found in Ref. [START_REF] Moorthy | A Voronoi Cell finite element model for particle cracking in elastic-plastic composite materials[END_REF].

Summary of fracture parameters

As a conclusion, two stress-based fracture criteria have been formulated.

Debonding is predicted based on the normal stress at the interface, and leads to a crack at a part of the matrix/particles interface. This model adds two parameters σ D c and D . Fragmentation is predicted based on the average first principal stress per particle, and leads to a crack across a whole particle. This model adds a parameter σ F c . The numerical method used to insert these cracks through LS functions and mesh adaption requires two numerical parameters c and f . All parameters are summarized in Tab. 1. Thanks to these models, the two modes of void nucleation can be predicted and modeled during RVE simulations. This evaluation of the criteria and the subsequent modifications on the microstructure and the mesh are performed at the beginning of each load increment, resulting in a weak coupling between the nucleation models and constitutive equations. Time step sensitivity is investigated in the Appendix. criterion material parameter(s) numerical parameter(s)

debonding σ D c (MPa) D (µm), c (µm) fragmentation σ F c (MPa) f (µm), c (µm)
Table 1: Material and numerical parameters for the void nucleation criteria.

Boundary conditions

As widely used in RVE simulations [START_REF] Benzerga | Ductile Fracture by Void Growth to Coalescence[END_REF], periodic boundary conditions are imposed 2 . A cubic RVE is used where:

• A constant normal tensile velocity v α is applied on the upper face.

• Symmetry boundary conditions are applied on the lower, left, and back faces.

• A constant normal compressive velocity v β is applied on the remaining faces.

This RVE is illustrated in Fig. 2, along with a similar setting for 2D plane strain simulations. The word RVE is used to denote the meshed domain, although in strict sense the unit cell that is actually being periodically repeated is the whole domain in Fig. 2. This is evidenced in Fig. 3.

2 Strain localization patterns may be influenced by the choice of periodic boundary conditions at late stages of loading, close to final failure. Thus, the present analysis will focus on early localization and void coalescence events. More appropriate choices of boundary conditions are being investigated [START_REF] Buljac | Numerical validation framework for micromechanical simulations based on synchrotron 3D imaging[END_REF][START_REF] Shakoor | On the choice of boundary conditions for micromechanical simulations based on 3D imaging[END_REF]. While the velocity v α is fixed, the value of v β is controlled in order to maintain a constant macroscopic stress state within the RVE. The macroscopic Cauchy stress tensor Σ is defined as the average of the Cauchy stress tensor σ over the RVE [START_REF] Benzerga | Ductile Fracture by Void Growth to Coalescence[END_REF]. To identify the stress state, three invariants of Σ and its deviatoric part Σ are defined

(a) (b) v α v β v α v β v β v β v β v β v α v α
Σ m = Σ ii 3 , Σ = 3 2 Σ ij Σ ij , J 3 = det(Σ ), (4)
where Σ m is the macroscopic mean stress, Σ the macroscopic von Mises equivalent stress, and J 3 the third invariant of the macroscopic deviatoric stress tensor. The first two invariants Σ m and Σ are used to define the macroscopic

stress triaxiality ratio η η = Σ m Σ . (5)
A detailed presentation of these invariants and their values for some widely used specimen geometries can be found in Ref. [START_REF] Bao | On fracture locus in the equivalent strain and stress triaxiality space[END_REF]. In particular, the axisymmetric tension condition used herein is identified by J 3 = 0. A triaxiality ratio η 1 3 means increasing tensile stresses in the directions perpendicular to the loading direction, while a triaxiality ratio η -1 3 means increasing compressive stresses in the directions perpendicular to the loading direction. The remaining range -

1 3 ≤ η ≤ 1 3 corresponds to shear-dominant conditions.
All these settings are interesting as they can be met locally in various material forming processes and also during extreme events, and no model existing in the literature succeeds in predicting damage change for any J 3 and any η [START_REF] Bao | A Comparative Study on Various Ductile Crack Formation Criteria[END_REF][START_REF] Bao | On fracture locus in the equivalent strain and stress triaxiality space[END_REF].

RVE calculations with more and more accurate microscopic models are hoped to give some clues on how to build such a model.

As stated above, J 3 = 0 due to loading axisymmetry. A linear search algorithm has been implemented to find the value of v β that would yield a desired macroscopic stress triaxiality ratio η.

While the problem is well-defined initially, crack initiation and propagation events may split the domain so that boundary conditions are not completely defined for some regions. This is avoided herein by automatically removing these isolated regions, as illustrated in Fig. 4. This treatment may have an influence for nonproportional loading paths (e.g., cracks opened under tension may close under shear or compression) but is acceptable for the present study where only monotonic tensile loading is considered.

(a) (b)

Application

The main objective of the following simulations is to show that void nucleation induces an early localization mechanism before the onset of void coalescence, reducing drastically ductility. Sensitivity of the results with respect to numerical parameters, which include the mesh refinement parameters, the load step, c and f , is investigated in the Appendix. Using the reference numerical parameters chosen in the Appendix, reliable results are obtained in terms of averaged quantities, but discussions regarding local phenomena such as the crack propagation path can be only qualitative. A systematic experimental-numerical validation framework will be considered in future work to enable quantitative comparisons [START_REF] Buljac | Numerical validation framework for micromechanical simulations based on synchrotron 3D imaging[END_REF][START_REF] Shakoor | On the choice of boundary conditions for micromechanical simulations based on 3D imaging[END_REF].

Material

An interesting experimental study on void nucleation can be found in Ref. [START_REF] Babout | Characterization by X-ray computed tomography of decohesion, porosity growth and coalescence in model metal matrix composites[END_REF], where ductile fracture of a particle reinforced aluminum alloy was investi-gated. This material was composed of an aluminum matrix reinforced with 20% of Zirconia-Silica (ZS) particles, and particle fragmentation was observed to be the dominating void nucleation mechanism.

In Ref. [START_REF] Babout | Characterization by X-ray computed tomography of decohesion, porosity growth and coalescence in model metal matrix composites[END_REF], material parameters were identified for the matrix by carrying tensile tests on purely matrix material specimens, and for the particles by using nanoindentation. Some uncertainties may arise due to the limitations of these identification techniques. First, the forming process is different for the purely matrix material and the reinforced one, so the behavior of the matrix may also be different. Second, nanoindentation can only be applied to the particles at the surface of the specimen, hence some differences may apply due to polishing for instance.

Nevertheless, such advanced identification is not commonly found in the literature. In particular, the fracture properties of the particles were also determined and are given in Tab. 2. These properties were determined by comparing X-ray tomography data acquired in situ with RVE calculations. The latter considered a single particle with periodic boundary conditions. The macroscopic equivalent plastic strain was identified by measuring the section reduction of the specimen, while the macroscopic stress triaxiality was determined using Bridgman's formula. • a nucleation simulation, taking into account both void nucleation mechanisms of particle debonding and fragmentation;

Material behavior

• a fragmentation simulation, taking into account only the particle fragmentation mechanism;

• a pre-fragmentation simulation, taking into account initially fragmented particles instead of modeling void nucleation during the simulation;

• a void simulation, taking into account voids instead of particles 3 .

The four RVEs are shown in Fig. 5 before and after 15 % of RVE elongation. On the one hand, the load is nearly equivalently shared by all voids in the void and pre-fragmentation simulations. As a consequence, void growth is nearly homogeneous in the whole RVE for these two simulations, apart from some void coalescence and linkage events. Thus, the stress and plasticity localization mechanisms that could take place in some regions of the RVE before the onset of void coalescence are totally neglected.

On the other hand, plasticity is localized in the top region of the RVE in the fragmentation and nucleation simulations. The microstructure fails due to the three aligned particles at the top of the RVE. The stress first decreases significantly due to the failure of the two big particles (≈ 8 % of RVE elongation), and then later a minor decrease is noticed due to the third particle (≈ 10 % of RVE elongation). In fact, nucleation due to debonding or fragmentation of a particle leads to a stress localization phenomenon that is likely to trigger further nucleation events in neighboring particles [START_REF] Ganguly | Rearrangement of local stress and strain fields due to damage initiation in a model composite system[END_REF].

Results with a modeling of nucleation mechanisms (only fragmentation or both) during the simulation are hence qualitatively very different from results with an initial population of voids or pre-fragmented particles. A more quantitative idea is given by the porosity change and macroscopic von Mises equivalent stress change curves in Fig. 6. These differences are explained when looking at the stress curves in Fig. 6(b).

Due to the presence of undamaged particles, the nucleation and fragmentation

RVEs have an ultimate strength twice higher than in simulations that do not consider void nucleation criteria.

The difference between the nucleation and fragmentation simulations is mainly concentrated in the first nucleation event, which is debonding in the former while it is fragmentation in the latter. Debonding has a less significant effect on the load carrying capacity (at ≈ 8 % of RVE elongation in Fig. 6(b)). This difference does not seem to have a significant effect after 20 % of RVE elongation, in the void growth and coalescence stages.

An important aspect is that while the void and pre-fragmentation simulations show very different void growth curves, with different slopes, the stress/strain curves correspond quite well. This reveals that the difference between pennyshaped voids and spherical voids is not significant under tensile loading. It may nevertheless have a major influence under different loading conditions, such as shear. The presence of particle fragments in the pre-fragmentation configuration may also have a significant influence under such loading conditions.

Additionally, the artificial thickness f may also reduce significantly the influence of particle fragments in these simulations. This may also certainly explain why the macroscopic von Mises equivalent stress in the pre-fragmentation result is lower than that in the void result in Fig. 6

(b).

As a conclusion, this 2D FE analysis shows that while the particle debonding mechanism may be neglected for the studied material and loading conditions, the fragmentation mechanism has to be accounted for. In particular, ultimate strength is significantly underestimated if particles are assumed to be already fragmented in the initial RVE, or if they are replaced by voids.

Three-dimensional analysis

The importance of void nucleation modeling is investigated further in 3D

with an RVE of size 165 µm (33 particles), and varying macroscopic stress triaxiality ratio η. For these 3D simulations, due to a large computational cost, only the two extreme cases are considered, namely the nucleation and void configurations of Paragraph 4.2.1. For instance, the 3D nucleation simulation at η = 0.33 takes 46 hours to run on a 2.5 GHz Intel Xeon Linux cluster. Due to mesh adaption and fracture events, the number of elements varies throughout the simulation (from 0.5 to 1.5 million), and so does the number of parallel CPUs (from 30 to 128). Results are presented in Fig. 7 for the nucleation RVEs and Fig. 8 for the void RVEs. Independently of whether particles and void nucleation mechanisms are taken into account, larger voids and more numerous void coalescence and linkage events are visible in the results corresponding to η = 0.66. In particular, most voids in the results with η = 0 have regular and nearly spherical shapes, which shows that coalescence is not the dominating mechanism yet, even at 50 % of RVE elongation.

Although this increase of porosity is well captured by the void simulations, there are significant differences between the void and nucleation simulations. These differences seem to be increased at lower triaxiality ratios. With η = 0.66, void coalescence and linkage have led to the localization of damage in a region of the RVE, while void growth out of this region is moderate (Fig. 8). This is not true for the simulation with η = 0, where all voids seem to have flattened homogeneously in Fig. 8, while some undamaged particles are still visible in Fig. 7 A remarkable result is that for η = 0, no void growth is observed in the void simulation (Fig. 9(a)). Voids only flatten for η = 0 in Fig. 8, but the overall void volume does not vary. Anisotropic versions of the Gurson model have been proposed in the literature to account for this shape change effect [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-Case of axisymmetric prolate ellipsoidal cavities[END_REF].

The absence of void growth explains that the macroscopic von Mises equivalent stress for η = 0 does not decrease during the simulation in Fig. 9(b). However, void flattening in a direction could have an influence if the loading direction varies.

Moreover, the void growth curves for η = 0 reveal the major influence of the void nucleation mechanisms. Void growth is observed in the nucleation simulation (Fig. 9(a)), with a small loss of load carrying capacity (Fig. 9(b)). The ultimate strength is also higher due to the presence of particles.

At higher triaxiality ratios (Fig. 9(a)), the differences between the nucleation and void simulations regarding porosity curves are reduced, as the slope of these curves are similar at large plastic strain. The macroscopic von Mises equivalent stress curves also show similar changes for large plastic strain (Fig.

9(b)).

This is however not true in the first part of loading, before 10 % of RVE elongation. In these first 10 %, there is still a difference regarding the ultimate strength, although this error seems to decrease when the triaxiality ratio increases. This is explained by the fact that void coalescence is the dominant damage mechanism at high triaxiality ratio, and not void nucleation. Small drops in the macroscopic von Mises equivalent stress curve in Fig. 9(b) are only visible for an RVE elongation below 15 % when η = 0.66. After 15 % of RVE elongation, the smoothness of the curve indicates that void nucleation is negligible.

It is important to underline that in all Gurson-type micromechanical models, void nucleation is considered based on the porosity, which considers only void volume (and sometimes also particle volume). Thus, cracks that nucleate but do not open and grow significantly cannot be described with a Gurson-type micromechanical model, while their influence on the load carrying capacity may not be negligible. In fact, even if this influence was negligible for a given monotonic loading, this would most certainly not be true if the loading direction was then modified, as in most material forming processes.

The plastic localization zone appears more clearly in Fig. 10, where plastic strain for the uniaxial tension configuration is shown. The range of RVE elongation between the first failure step (due to nucleation) and the second one (due to coalescence) is larger for low triaxiality ratios, and smaller for high triaxiality ratios. While in classical Gurson-type models, an additive decomposition of void nucleation and coalescence terms is assumed, these results indicate that for high triaxiality ratios this assumption is not valid anymore, and both phenomena should be accounted for simultaneously in a unified criterion. This has been done recently in the case of void growth and coalescence [START_REF] Morin | A unified criterion for the growth and coalescence of microvoids[END_REF], but still has to be studied for void nucleation.

Conclusions

The numerical methods and micromechanical models proposed in the present work enable for simulations of ductile fracture after the onset of void nucleation, up to void coalescence (Fig. 10). Applications to large and complex 3D RVEs with particles of random size and position are possible, but restricted to a ductile material with high particle volume fraction in this work. Multiple void linkage events occur in these RVEs and are captured by the robust mesh adaption procedure (Figs. 7 and8). This purely plasticity-driven void coalescence and linkage model will have to be coupled to matrix cracking criterion in order to model ductile materials with lower particle and void volume fraction. The initiation of multiple cracks during simulations due to particle debonding and fragmentation is also modeled in this work. To the author's best knowledge, this is the first time that simulations taking into account all these micromechanisms at large plastic strains are conducted. Future work will include further improvements of the numerical methods to eliminate the remaining artificial parameters and model the dissipated energy rate. Experimental-numerical validation is being investigated to identify fracture parameters directly at the microscale and enable quantitative comparisons with 3D data.

When compared to less advanced approaches where particles and/or nucleation mechanisms are neglected, significant differences are observed, hence

showing the importance of modeling particles and the subsequent void nucleation mechanisms. It is obvious that these simulation tools are too costly to be applied directly to model plasticity and mechanical behavior at industrial scale.

For these applications, classical Gurson-type models are now well accepted to be the most promising approaches. However, these models usually consider the three aspects of void nucleation, growth, and coalescence mechanisms assuming additive decomposition, based on micromechanical or phenomenological considerations. In order to improve predictions, these assumptions should be reconsidered.

It is very difficult to engineer materials with given microstructures, for instance with a given initial arrangement of particles or voids. In numerical methods, however, this is greatly simplified. Computational approaches also allow to compute accurately relevant mechanical variables and study relations between them. In the future, the present numerical framework should be used to assess the capabilities of classical Gurson-type models, especially regarding non proportional loading paths. Such comparisons between full field and mean field models should give insights into the limitations of mean field models, and how they could be improved.

Appendix

In this appendix, numerical parameters such as meshing parameters and the load step are chosen based on a sensitivity analysis. The phrasing load step is used instead of time step as all calculations are performed in static loading conditions, and material constitutive laws are rate-independent (Subsection 2.1).

The same material and boundary conditions as defined in Subsection 4.1 are used. A global picture of the methodology developed in Section 3 and also in previous works to model ductile fracture at the microscale is presented in Algorithms 1 and 2. The prescribed macroscopic stress triaxiality ratio is denoted η user and, unless otherwise mentioned, fixed to 0.33.

Algorithm 1 Micromechanical simulation methodology with references to previous work and the new developments introduced in the present paper. procedure MainAlgorithm

Compute distance functions

Ref. [START_REF] Shakoor | An efficient and parallel level set reinitialization method -Application to micromechanics and microstructural evolutions[END_REF] Compute particle tags Ref. [START_REF] Scholtes | New finite element developments for the full field modeling of microstructural evolutions using the level-set method[END_REF] v β ← 0 Paragraph 3.2

for Each load increment do Update node coordinates (mesh motion) Ref. [START_REF] Shakoor | An adaptive level-set method with enhanced volume conservation for simulations in multiphase domains[END_REF] Write output files end for end procedure

The 3D computational domain has a size of 55 µm, and contains a single particle of radius 20 µm. The chosen reference numerical parameters are:

• an isotropic mesh size between h max = 6 µm and h min = 2 µm, depending on the local maximum principal curvature at phase boundaries (Subsection 2.2),

• a crack thickness c = 0.2 µm and an artificial fragmentation crack extension length f = 10 µm (Subsection 3.1),

• a load step of 0.1 % RVE elongation, which is decreased by a factor of two every 20 iterations of the Newton-Raphson algorithm if convergence was not reached (Subsection 2.1).

A sensitivity analysis is conducted by varying numerical parameters and investigating the influence of this variation on relevant mechanical variables.

Because the present framework relies extensively on remeshing operations and local fracture criteria (Algorithm 2), all parameters that depend on mesh size are expected to have a significant influence. Results may also be sensitive to the load step, because the fracture criteria are explicit (Algorithm 1).

For the analysis conducted in this paper, only averaged quantities such as the porosity f and the macroscopic von Mises equivalent stress Σ are of interest.

However, for future comparisons with experiments, an accurate prediction of local variables is also important. This is verified based on the maximum normal stress at the matrix/particle interface σ n and the particle-wise average maximum stress σ 1 . The present sensitivity analysis includes both averaged and local quantities. Sensitivity to mesh size is assessed using a coarse mesh setting with h max , h min , and c multiplied by a factor of 2, and a fine mesh setting with a factor of 0.5. The three meshes are shown in Figure 11. Sensitivity to the load step is assessed using load steps of 1 % and 0.01 %, in addition to the reference one of 0.1 %.

coarse mesh reference mesh fine mesh The results obtained using the reference numerical parameters are shown in Fig. 12. Only particle fragmentation is predicted, which is in agreement with experimental observations [START_REF] Babout | Characterization by X-ray computed tomography of decohesion, porosity growth and coalescence in model metal matrix composites[END_REF]. Both σ n and σ 1 increase up to values that are close to the critical thresholds, but the fragmentation threshold is reached first, at nearly 3 % of RVE elongation.

Regarding the fragmentation crack, it is a horizontal fracture plane, which is also in agreement with experimental observations [START_REF] Babout | Characterization by X-ray computed tomography of decohesion, porosity growth and coalescence in model metal matrix composites[END_REF]. Because of the small RVE size, due to the fact that the particle volume fraction is initially fixed to 20 %, the plane is already in contact with the RVE boundary at its nucleation. This is not physical and is linked to the numerical parameter f , whose influence is investigated in Subsection 6.2.

Additionally, fragmentation cracks could be accompanied with a debonding of the matrix/particle interface close to the fragmentation crack tip [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF]. Although this is not the case for the present material [START_REF] Babout | Characterization by X-ray computed tomography of decohesion, porosity growth and coalescence in model metal matrix composites[END_REF], the presence of a quite large part of the fragmentation crack inside the matrix would prevent this phenomenon.

Regarding the crack thickness parameter c , its influence cannot be seen directly in Fig. 12(b). Since the fracture criteria are considered at the beginning of each load increment, the porosity shown in this figure includes both the artificially inserted porosity due to c and f , and the porosity due to crack opening.

In order to distinguish the two, the artificial porosity is computed separately during the simulation. It is equal to 0.30 % at the onset of particle fragmentation. It then decreases because no additional artificial porosity is inserted, while the global void volume increases due to void growth. This is a major improvement over previous work [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF][START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF], where a mesh size at least ten times finer was required to reduce the artificial porosity to similar values [START_REF] Roux | A new finite element approach for modelling ductile damage void nucleation and growth-analysis of loading path effect on damage mechanisms[END_REF]. Moreover, this artificial porosity is negligible when compared to the final porosity, which is close to 20 %.

The porosity curve in Fig. 12 Based on this reference result, a sensitivity analysis is conducted regarding numerical parameters. In Fig. 13, it is verified that reference parameters enable an accurate prediction of the onset of particle fragmentation.

The fact that σ 1 is null for the first loading step is simply due to the fact that the fragmentation criterion is computed at the beginning of each load increment.

The threshold for particle fragmentation is never reached using the coarse load step (i.e., particle debonding is predicted instead). This means that this load step is too large to enable for the modeling of the competition between the two void nucleation mechanisms. This result is hence discarded in the remaining of Opposed to simulations with varying load steps, particle fragmentation is predicted in all these simulations. However, the onset of particle fragmentation is underestimated using coarse meshing parameters. The absolute uncertainty on RVE elongation at the onset of particle fragmentation for the reference mesh remains small, namely 1.0 %.

The influence of this uncertainty on void growth is shown in Fig. 16.

The results are similar to those obtained in the load step sensitivity analysis. Particle fragmentation is delayed but this does not affect void growth, as the slopes of the curves correspond well between the reference result and the result with fine mesh. The absolute uncertainty on the porosity at 30 % of RVE elongation is 0.04 %. The reference results are hence well-suited for the present material. Moreover, the fact that this uncertainty is very small indicates that the crack thickness parameter c has, in the end, a very low influence on the results.

As a conclusion, absolute errors regarding the prediction of the onset of par- ticle fragmentation close to the percent are obtained using reference numerical parameters. These uncertainties mainly delay void growth, but do not affect significantly the slope of the porosity change curves. The absolute uncertainty on the final porosity is hence also close to the percent.

Particle fragmentation crack extension length f

As mentioned in Subsection 6.1, the numerical implementation of the particle fragmentation criterion includes a length f that leads to non physical results. In this paragraph, simulations with the reference numerical parameters but varying f are conducted to study the influence of this parameter. Results are reported in Fig. 17. Due to a small distance between the particle and the RVE boundary, the sensitivity to f is very high, both regarding the load carrying capacity and the porosity, due to an overestimation of void coalescence. Even for the simulation with f = 2.5 µm, it can be seen that the remesher is already initiating void linkage with the RVE boundary, due to a low number of elements between the fragmentation crack tip and the RVE boundary. The present sensitivity analysis indicates that the numerical parameter f should be exactly zero in order to obtain an accurate and physical modeling of particle fragmentation.

In the future a fragmentation model replacing crack tips inside the matrix by crack tips along the particle/matrix interface will be considered.

Particle debonding

For particle debonding to be predicted instead of fragmentation, the particle fragmentation criterion is deactivated. Debonding is then triggered, but for larger strains, as shown in Fig. 18.

While particle fragmentation in Subsection 6.1 is predicted at an RVE elongation of nearly 3 %, particle debonding does not occur until 5 %. An interesting study is to compare the effect of particle debonding on the load carrying capacity to the effect of particle fragmentation. For instance, the macroscopic von Mises equivalent stress curve in Fig. 18(b) can be compared to that in Fig. Sensitivity analysis results regarding the load step with deactivated particle fragmentation are shown in Fig. 19. The absolute error on the onset of particle debonding using the reference load step is 1.8 %. It is important to check the influence of this error on porosity change (Fig. 20).

The same remarks as for particle fragmentation apply. In particular, the delay regarding the onset of void nucleation does not seem to affect void growth, as the slopes of the porosity change curves correspond for all used load steps.

The absolute error on the porosity at 50 % of RVE elongation is 0.7 %. These results validate the use of the reference load step for further FE analysis.

The dependence on mesh size is assessed in Fig. 21 The absolute uncertainty regarding the prediction of the onset of particle debonding using the reference mesh is very low (0.01 %). However, porosity change curves in Fig. 22 show a significant dependence on mesh size. A possible explanation is the sensitivity of the particle debonding site with respect to numerical parameters. The calculation of the particle debonding indicator is accurate, but it reaches values very close to the critical threshold at multiple points of the interface. The debonding of one of these points instead of the other is purely linked to numerical errors. This effect is illustrated in Fig.

0%

23.

The result with a coarse mesh overestimates void coalescence due to a small number of elements between the particle and lateral RVE boundaries. This explains the quite large difference between the results obtained using a coarse mesh and the other results in Fig. 22.

coarse mesh reference mesh fine mesh

Conclusion

Except for the artificial crack extension length, the reference numerical parameters chosen in the introduction of this appendix, and used in the remainder of this paper, give accurate predictions regarding fracture criteria and averaged quantities. In particular, the absolute error on the onset of both particle fragmentation and debonding is estimated to be below 2 %, as well as the absolute error on the porosity. The effect of the artificial crack thickness c and the amount of artificial porosity introduced because of this parameter is also shown to be negligible, which is a major improvement over previous work.

The artificial crack extension length f , however, should be removed. It has a major influence on both local and averaged results. Additionally, the energy dissipation rate at the onset of fragmentation or debonding should be modeled using a progressive traction-separation law. It can nevertheless be assumed that the conclusions raised in this paper regarding the importance of void nucleation modeling with competing particle debonding and fragmentation mechanisms would not be changed with a more physical fragmentation model.

Figure 2 :

 2 Figure 2: RVE under axisymmetric tension in: (a) 2D, (b) 3D. Thanks to symmetry boundary conditions, the part in light color is not meshed.

Figure 3 :

 3 Figure 3: Example of 2D periodic microstructure with symmetries, the part in light color is not meshed.

Figure 4 :

 4 Figure 4: Elimination of matrix regions and particle fragments that could lead to rigid body modes during mechanical solution: (a) initial 2D microstructure with rigid body modes, (b) resulting microstructure after the numerical treatment.

Figure 5 :

 5 Figure 5: Microstructures before and after 15 % of RVE elongation for the assessment of the importance of void nucleation modeling in 2D plane strain. Particles are shown in red, and voids in white.

Figure 6 :

 6 Figure 6: Importance of void nucleation modeling in 2D plane strain.

Figure 7 :

 7 Figure 7: Particles (red) and voids (light gray) at different RVE elongations for the 3D nucleation computations at different macroscopic stress triaxiality ratio η.

Figure 8 :

 8 Figure 8: Voids at different RVE elongations for the 3D void computations at different macroscopic stress triaxiality ratio η.

  Macroscopic von Mises equivalent stress change for the three simulations.

Figure 9 :

 9 Figure 9: Importance of void nucleation modeling in 3D.

Figure 10 :

 10 Figure 10: Equivalent plastic strain for the 3D computation taking into account nucleation mechanisms at macroscopic stress triaxiality ratio η = 0.33 and RVE elongation of: (a) 25 %, (b) 50 %.

FractureCriteria Algorithm 2 repeat

 2 Mechanical solution using v α , v βRef.[START_REF] Roux | A level-set and anisotropic adaptive remeshing strategy for the modeling of void growth under large plastic strain[END_REF] Compute η Paragraph 3.2 Update v β (linear search) Paragraph 3.2 until |η -η user | < 0.01 Paragraph 3.2

Figure 11 :

 11 Figure 11: Inside view of the three meshes used for analysis of mesh size sensitivity.

Figure 12 :

 12 Figure 12: Reference result in 3D featuring only particle fragmentation.

Figure 13 :

 13 Figure 13: Sensitivity of the particle fragmentation indicator σ 1 during loading with respect to the load step.

Figure 15 :

 15 Figure 15: Sensitivity of the particle fragmentation indicator σ 1 during loading with respect to mesh size.

  Microstructures at the onset of particle fragmentation, with particle in red and voids in white. Evolution of the macroscopic von Mises equivalent stress Σ and the porosity f .

Figure 17 :

 17 Figure 17: Sensitivity of the results with respect to the fragmentation crack extension length f in 3D.

17

 17

  (b) using a particle fragmentation crack extension length f = 2.5 µm. The two curves are very similar, regarding the beginning and the end.However, particle fragmentation leads to a more severe and instantaneous drop of the load carrying capacity (Σ ≈ 250 MPa at 3 % of RVE elongation) while this is less significant and more progressive in the case of debonding (Σ ≈ 300 MPa at 5 % of RVE elongation). Void growth is also slowed down in the case of debonding, as the porosity reaches 10 % at 50 % of RVE elongation (Fig.18(b)), while this occurs at only 30 % of RVE elongation (Fig.17(b)) in the particle fragmentation case with f = 2.5 µm. Thus, while particle fragmentation leads to an instantaneous void nucleation mechanism, particle debonding is more progressive. This effect is not necessarily due to the fact that debonding cracks may still be propagating at the onset of void coalescence. Only one void nucleation mechanism appears clearly in Fig.18(b). Debonding cracks are propagating during the whole simulation (Fig.18(a)), but only the first debonding event at 5 % of RVE elongation affects clearly the load carrying capacity. Later debonding Microstructure at different RVE elongations, with particle in red and voids in white. equivalent stress porosity (b) Evolution of the particle debonding indicator σn (threshold at 1060 MPa), the macroscopic von Mises equivalent stress Σ, and the porosity f .

Figure 18 :

 18 Figure 18: Reference result in 3D featuring only particle debonding, the fragmentation criterion being deactivated.

Figure 19 :

 19 Figure 19: Sensitivity of the particle debonding indicator σn during loading with respect to the load step.

Figure 20 :

 20 Figure 20: Sensitivity of the porosity during loading with respect to the load step.

Figure 21 :

 21 Figure 21: Sensitivity of the particle debonding indicator σn during loading with respect to mesh size.

Figure 22 :

 22 Figure 22: Sensitivity of the porosity during loading with respect to mesh size.

Figure 23 :

 23 Figure 23: Microstructures using various meshing parameters, with particles in red and voids in white.

Table 2 :

 2 Material parameters for the Al2124 matrix and the ZS particles, including fracture criteria parameters[START_REF] Babout | Characterization by X-ray computed tomography of decohesion, porosity growth and coalescence in model metal matrix composites[END_REF], to which is added the numerical parameter D . Young's modulus is denoted E and Poisson's ratio ν, while plastic properties refer to Eq. (1).

		Phase	E (MPa)	ν	σ y (MPa) K (MPa)	n
		Matrix		72000	0.33	270	580	0.54
		Particles		123000	0.23	∞
	Fracture criteria	Phase		σ D c (MPa)	D (µm) σ F c (MPa)
		Particles	1060		10	700
	In the present study, particles are considered as spherical, with radii dis-
	tributed according to a normal law of mean value 20 µm and standard deviation

Unless otherwise mentioned, boundary conditions as defined in Section 3.2 are imposed and the macroscopic stress triaxiality ratio η is maintained to 0.33.

4.2. Importance of void nucleation modeling

4.2.1. Two-dimensional analysis It has been discussed in the introduction that most studies on ductile fracture often disregard the void nucleation mechanism by considering an initial microstructure where particles are replaced by voids or are initially already debonded and/or fragmented. A 2D RVE of size 240 µm (9 particles) is meshed, and comparisons are proposed between:

  at 50 % of RVE elongation. As a conclusion, these results confirm that early stress and plasticity localization phenomena due to void nucleation are neglected in the void simulations. Additionally, the effect of this assumption is increased at low triaxiality ratios.

	These effects are investigated further in terms of global porosity and load car-
	rying capacity in Fig. 9.
		45%	nucleation η=0,00		void η=0,00
		40%	nucleation η=0,33		void η=0,33
		35%	nucleation η=0,66		void η=0,66
		30%
	Porosity	20% 25%
		15%
		10%
		5%
		0%
		0%	10%	20%	30%	40%	50%
				RVE elongation

  .

	14	%
			load step 0.01%
	12	%	load step 0.1%
	10 %	load step 1%
	8 %
	6 %
	4 %
	2 %
	0 %
		0 %	10 %	20 %	30 %	40 %	50 %
				RVE elongation

Mesh adaption would merge nodes that are superimposed[START_REF] Shakoor | An adaptive level-set method with enhanced volume conservation for simulations in multiphase domains[END_REF]. A proper computational contact mechanics implementation will be investigated in the future to solve this issue.

Due to the absence of a contact model, this configuration is here equivalent to initially debonded particles.

Compute mesh size criterion from distance functions Ref. [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF] Optimize mesh quality through local mesh modifications Ref. [START_REF] Shakoor | An adaptive level-set method with enhanced volume conservation for simulations in multiphase domains[END_REF] Update distance functions on optimal mesh Ref. [START_REF] Shakoor | An efficient and parallel level set reinitialization method -Application to micromechanics and microstructural evolutions[END_REF] Update particle tags on optimal mesh Ref. [START_REF] Scholtes | New finite element developments for the full field modeling of microstructural evolutions using the level-set method[END_REF] Transfer history variables from old mesh to optimal mesh Ref. [START_REF] Shakoor | An adaptive level-set method with enhanced volume conservation for simulations in multiphase domains[END_REF] Replace old mesh with optimal mesh end procedure this analysis.

The differences between the curves using the reference and fine load steps are small. The measured absolute uncertainty of 0.08 % of RVE elongation regarding the prediction of particle fragmentation is acceptable using the reference load step.

In Fig. 14, the influence of this uncertainty on porosity change is assessed. No significant influence of the load step on void growth can be observed.

The slope of the curve obtained using the reference load step corresponds well to the slope obtained using the fine load step. The main difference is due to the delay regarding the prediction of the onset of particle fragmentation. The corresponding absolute error on the porosity at 30 % of RVE elongation is 0.1 %. This is acceptable for the present material with high porosity, but should be investigated further for industrial materials where the porosity is lower.

The influence of mesh size on the prediction of the onset of particle fragmentation is assessed in Fig. 15.

Acknowledgments

The supports of the Carnot M.I.N.E.S Institute and the French Agence Nationale de la Recherche (ANR) through the COMINSIDE project are gratefully acknowledged.