**3. Model**

#### *3.1. Gurson–Tvergaard–Needleman Damage Model*

The plastic behavior of the open-cell alloy foam is directly influenced by the properties of the solid phase. The fracture behavior of the solid aluminum phase is governed by ductile damage, namely the nucleation, growth, and coalescence of cavities. Nucleation is mainly due to microcrack initiation (e.g., from pores or defects in the material [30]), second phase particle decohesion, or fracture in the alloy [31]. The growth of cavities occurs by plastic yielding of the matrix surrounding the cavity. Coalescence occurs when neighboring cavities or cracks merge together [32].

In this study, a standard Gurson–Tvergaard–Needleman (GTN) model [33–37] is used to represent the effect of void nucleation and growth in the simulation. The GTN model renders the effect of porosity on the yield locus and its sensitivity to the hydrostatic component of loading. It reduces to a standard isotropic von Mises yielding criterion in the absence of porosity. The governing equation is defined as:

$$\Phi\left(\sigma\_{t\eta}, \sigma\_{\mathcal{Y}}, \sigma\_{\mathcal{H}}, f\right) = \left(\frac{\sigma\_{t\eta}}{\sigma\_{\mathcal{Y}}}\right)^2 + 2fq\_1 \cosh\left(\frac{3q\_2\sigma\_{\mathcal{H}}}{2\sigma\_{\mathcal{Y}}}\right) - \left(1 + q\_3f^2\right) = 0,\tag{1}$$

where Φ is the yield function, *q*1, *q*2, and *q*3 are calibrating parameters, *<sup>σ</sup>eq* is the von Mises equivalent stress, *<sup>σ</sup>y* is the yield stress, *σH* is the hydrostatic stress, and *f* is the void volume fraction (VVF) in the matrix. In addition, starting from an initial void volume fraction *f*0, the total change in *f* , noted ˙ *f* , is defined as [34]:

$$f = f\_{\mathcal{S}^{\mathcal{T}}} + f\_{\text{mucl}} = (1 - f) \operatorname{tr}(\dot{\varepsilon}^{pl}) + \frac{f\_N}{s\_N \sqrt{2\pi}} \exp\left[ -\frac{1}{2} \left( \frac{\varepsilon\_{eq}^{pl} - \varepsilon\_N}{s\_N} \right)^2 \right] \dot{\varepsilon}\_{eq\prime}^{pl} \tag{2}$$

where ˙ *fnucl* is the contribution of nucleating voids; ˙ *fgr* is the void growth rate, which is based on mass conservation and directly proportional to the hydrostatic component of plastic strain rate tensor *tr*(*ε*˙ *<sup>p</sup><sup>l</sup>*); and *εpl eq* is the equivalent plastic strain. *εN* and *sN* are the mean value and standard deviation of the normal nucleation distribution. *fN* is the volume fraction of the nucleated voids. The power law work hardening used in this study is defined as [38]:

$$\frac{\sigma\_{\mathcal{Y}}}{\sigma\_{0}} = (\frac{\sigma\_{\mathcal{Y}}}{\sigma\_{0}} + \frac{3G}{\sigma\_{0}} \varepsilon^{pl})^{N},\tag{3}$$

where *σ*0 is the initial yield stress, *N* is the hardening exponent, *G* the elastic shear modulus and *εp<sup>l</sup>* is the plastic strain. In addition, the initiation of necking in tension takes place according to a Considère criterion when the work hardening rate converges to the current yield stress:

$$\frac{\partial \sigma\_y}{\partial \varepsilon^{pl}} = \sigma\_y. \tag{4}$$

#### *3.2. Mesh Generation*

In this study, the initial 3D tomographic volumes of the foam samples were used to generate 3D image-based FE meshes. First, a surface mesh was generated from the solid phase boundaries in the volume. Next, the surface mesh was simplified and remeshed to reduce the number of triangles while preserving a proper description of the surface. Finally, the solid volume was filled by first-order tetrahedra [13] (C3D4 elements in the Abaqus software), allowing for explicit non-linear simulations including damage and fracture. The whole meshing procedure was performed using the commercial Avizo-R software [39]. Each foam sample was meshed with four different mesh sizes in order to investigate the mesh size sensitivity of the results, as detailed in Section 3.4. The number of tetrahedra and nodes of the reference volume meshes are detailed in Table 4.

**Table 4.** Characteristics of the reference volume meshes for the two foam samples.


Two classes of FE models were considered in this work: (i) *homogeneous* (or microstructure "blind") models, where the local constitutive behavior corresponds to the nominal average behavior of the aluminum alloy everywhere; and (ii) *heterogeneous* models, where the local constitutive behavior depends on the local microstructure (intermetallic fraction in the present case) as informed from the detailed initial 3D tomographic scans of the samples.

The procedure used to generate the *heterogeneous* FE model is described in [40] and briefly recalled here. For each element of the 3D image-based FE mesh, a python script retrieves the voxels of the segmented tomographic volume located in the element interior. The number of white (intermetallic) voxels is then evaluated in order to compute the local intermetallic volume fraction *fIM* in the element. As a result, the tetrahedra of the volume mesh are tagged depending on the volume fraction of intermetallic particles regrouped into 100 equidistant classes with *fIM* varying from 0.0 to 1.0. The distributions of the intermetallic volume fraction in the tetrahedra are illustrated in Figure 3 for the 20 PPI and 30 PPI foam samples.

**Figure 3.** Distributions of intermetallic particles volume fraction in the finite element (FE) tetrahedra for the 20 PPI and 30 PPI foam samples.

#### *3.3. Identification of the Constitutive Model Parameters*

The identification of constitutive model parameters is a challenging task, especially when it comes to the local behavior of the solid (metallic) phase in macroporous structures like the present foam samples. Several studies demonstrated that the tensile mechanical behavior of individual struts extracted from the foam might be significantly different from the behavior of the corresponding bulk metal [41–43], even subjected to the same hardening heat treatment. The problem becomes even more complicated if one would like to take the local microstructure (e.g., local intermetallic fraction here) of foam struts into account in the constitutive modeling. The present work proposes to rely on axisymmetric FE unit cell calculation to compute the local constitutive behavior of the different aluminum matrix composites with various intermetallic volume fractions that are available in each tetrahedron of foam sample mesh.

The 2D axisymmetric cell model is illustrated in Figure 4a. Varying the intermetallic (green) radius *r* allows to cover intermetallic volume fractions *fIM* from 0.0 to 0.6 in order to compute the resulting composite behavior to be considered for the corresponding elements in the distribution of Figure 3. A linear elastic behavior with a Young's modulus of *E* = 160 GPa and a Poisson's ratio *ν* = 0.33 was considered for the intermetallic properties [20]. The power-law hardening of Equation (3) was used for the aluminum phase, with *E* = 70 GPa, *ν* = 0.33 , an initial yield stress *σ*0 = 97 MPa, and a hardening exponent *N* = 0.052 . This initial yield stress corresponds to the lower range of individual strut yield stresses measured by [41] on a similar foam. The resulting tensile stress–strain curves obtained for increasing intermetallic volume fraction *fIM* are illustrated in Figure 4b. They were extracted along with the corresponding Young's moduli (ranging from 70 GPa for pure aluminum to 115 GPa for *fIM* = 0.6 ) to define the elasto-plastic properties of the foam tetrahedra in the different classes of Figure 3.

**Figure 4.** (**a**) 2D axisymmetric model of the aluminum matrix material (red) with a given intermetallic (green) volume fraction *fIM*. The right exterior line remains straight and vertical during loading, but is free to move in the radial direction. (**b**) Resulting hardening behavior for increasing intermetallic volume fraction.

Damage model parameters were also identified for each element class (each value of *fIM*). Following [36], the calibrating parameters of the GTN yield function in Equation (1) were chosen as *q*1 = 1.5, *q*2 = 1, and *q*3 = 2.25. In the nucleation part of Equation (2), the volume fraction of nucleating voids was taken as *fN* = 0.04 [20,44]. The initial void volume fraction *f*0 was taken as zero, which is consistent with the negligible initial porosity of the solid phase measured in the tomographic volumes (Section 2.1). *εN* and *sN* were gradually decreased with the increasing value of *fIM* to mimic the transition towards a brittle behavior with the increase in intermetallic fraction in the material (Table 5). Figure 5 illustrates the resulting tensile stress–strain behaviors, including GTN damage of the corresponding aluminum matrix materials. The increase in volume fraction of intermetallic particles clearly results in a stiffer, stronger, but more brittle behavior. This will impact the constitutive response of the tetrahedra with a high intermetallic fraction in the foam FE mesh.


**Table 5.** Parameters of the Gurson–Tvergaard–Needleman (GTN) nucleation model for increasing volume fraction of intermetallic particles.

#### *3.4. Simulation Conditions*

The FE approach described in the present study consists of a uniaxial monotonous tensile loading applied to the 3D image-based FE meshes generated in the previous section. Two nodes, one on the top surface and another on the bottom surface of the foam samples, were assigned as master nodes. *z* displacement degree of freedom of each node on the top surface was constrained to be the same

as the one of the top master node, which was assigned a positive displacement *uz*. Similarly, the *z* displacement degree of freedom of each node on the bottom surface was constrained to be the same as the one of the bottom master node, which was blocked. Furthermore, nodes located on two central lines on the top and bottom faces were constrained by imposing *ux* = 0 to prevent the potential rotation of the sample about the tensile axis, which is inhibited experimentally by the machine grips.

**Figure 5.** Tensile stress–strain curves of aluminum matrix materials including GTN damage, for increasing fraction of intermetallic particles.

The sensitivity of the simulation results with respect to the element size is illustrated in Figure 6a,b with the predicted apparent Young's modulus obtained with the homogeneous non-porous models. It converged to about 125 MPa with decreasing element size in the case of the 20 PPI foam sample and to about 210 MPa in the case of the 30 PPI foam sample. The characteristic element size of 150 μm was therefore chosen as the reference element size for both samples in Table 4. However, one can note that these predicted values of Young's moduli were overestimated with respect to the experimental measures based on the (3D) image-based macroscopic tensile strain. Two modeling choices in the present study might contribute to this situation: (i) The use of first-order tetrahedral elements imposed by the complexity of the structure and the compatibility with the Explicit solver of the Abaqus FE software. First-order tetrahedra are known to be too stiff in bending-dominated problems [38]. (ii) The uniaxial tension boundary conditions applied in the simulations slightly differ from the effective boundary conditions in the experiments.

**Figure 6.** Apparent Young's modulus convergence tests for (**a**) 20 PPI and (**b**) 30 PPI foam samples in the case of the homogeneous non-porous models.

#### **4. Results and Discussion**

Figures 7 and 8 respectively illustrate the progressive deformation of the 20 PPI and 30 PPI foam samples during the in situ tensile tests. Spatially uniform straining was observed until the first fracture happened in a strut. Both plastic stretching and plastic bending of the struts might be observed during the tensile deformation of foam blocks. Struts approximately parallel to the loading direction were mostly stretched, while the others were deformed first by bending and then potentially by stretching once they were aligned with the tensile direction. It was observed here that the majority of the struts were not parallel to the loading axis. Therefore, the deformation of the foam blocks was dominated by bending for small levels of elongation. The deformation then shifted to a stretching-dominated mode after significant elongation and re-alignment of the struts. It should be noted that bending-dominated deformation tended to lower the apparent stiffness and strength of the foam samples, while stretching dominated deformation rather tended to increase them.

Fracture in the foam blocks initiated in struts that were parallel to the loading direction and stretched. Subsequent propagation of fracture was achieved in the neighboring struts, mostly reclined with respect to the loading direction and exhibiting first bending and then stretching deformation prior to collapse. Once some struts broke, the fractured area expanded to the adjacent cells.

**Figure 7.** Deformation of the 20 PPI foam sample during the in situ tensile test at nominal strain: (**a**) 0.000; (**b**) 0.019; (**c**) 0.040; (**d**) 0.064. The blue arrow indicates the loading direction.

**Figure 8.** Deformation of the 30 PPI foam sample during the in situ tensile test at nominal strain: (**a**) 0.000; (**b**) 0.031; (**c**) 0.085; (**d**) 0.142. The blue arrow indicates the loading direction.

Figure 9a illustrates the five biggest void-cells of the 20 PPI foam sample, highlighted in red. Figure 9b shows the final step of the fractured foam sample during the in situ tensile test. The fracture

area was located in the vicinity of the three big neighboring cells, and passed through two of them (numbers 1 and 3, Figure 9b). In the case of the 30 PPI foam sample, the five biggest cells (Figure 10a) were more scattered than those of the 20 PPI sample. The fracture area was located between the five biggest cells and ended into two of them (numbers 1 and 3, Figure 10b). Furthermore, the 30 PPI foam sample, which was elongated in the transversal direction, was contracted significantly in the *x* and *y* directions.

**Figure 9.** (**a**) The five biggest void-cells and (**b**) the final fractured state of the 20 PPI foam sample. The blue arrow indicates the loading direction. Cells of interest are labeled 1 to 3.

**Figure 10.** (**a**) The five biggest void-cells and (**b**) the final fractured state of the 30 PPI foam sample. The blue arrow indicates the loading direction. Cells of interest are labeled 1 to 3.

In order to study the effect of intermetallic particles in the solid phase, four FE simulations were performed for each foam sample: (i) *homogeneous* non-porous J2 model, (ii) *heterogeneous* non-porous J2 model, (iii) *homogeneous* porous GTN model, and (iv) *heterogeneous* porous GTN model. J2 models here correspond to standard isotropic (von Mises) plasticity.

Figure 11a,b illustrate the macroscopic tensile curves obtained with the 20 PPI and 30 PPI foam samples in the experiments and simulations. The black points correspond to experimental measurements (note that only the first point laid in the linear elastic regime for each foam sample). The dashed lines correspond to the simulations with isotropic non-porous plasticity (J2 plasticity). The bold lines correspond to the simulations with the GTN damage model. The blue curves correspond to the *homogeneous* simulations with an aluminum matrix without intermetallic particles included. The red curves correspond to the *heterogeneous* simulations where the local presence of intermetallic particles with a given volume fraction *fIM* is taken into account in the constitutive behavior of the

elements. It is clear that taking the presence of intermetallic particles into account in the FE simulations *did not* affect the calculated macroscopic stress–strain curves significantly. However, the non-porous models did not capture the last stages of foam deformation properly, overestimating the reaction stress in the absence of damage.

The tensile curves of Figure 11 exhibited some macroscopic hardening behavior and a peak stress attained after a few percents of macroscopic deformation. The hardening arose from both the geometrical rearrangemen<sup>t</sup> of the struts during loading and the constitutive strain hardening of the aluminum solid phase. The 20 PPI foam sample exhibited a 55 MPa Young's modulus and 0.70 MPa tensile strength. The 30 PPI foam sample exhibited a 138 Pa Young's modulus and 0.87 MPa tensile strength, which was both stiffer and stronger than the 20 PPI sample due to the smaller size of the cells. One can note also the higher apparent ductility of the 30 PPI foam sample, accepting a larger level of macroscopic deformation before collapse.

**Figure 11.** Macroscopic experimental and FE tensile curves of the (**a**) 20 PPI and (**b**) 30 PPI foam samples.

Figure 12a,b illustrate the presence of intermetallic particles in a subregion of the 20 PPI foam sample, before the test and after fracture of the struts. It shows that strut numbers 1 and 5 broke in the vicinity of clusters of intermetallic particles. This was also the case of strut number 4, ye<sup>t</sup> with a smaller amount of particles. On the contrary, other regions with important clusters of intermetallic particles did not fail, illustrating the complexity of the situation leading to fracture. Several factors contributed and competed to trigger fracture, including the loading mode (or misalignment of the struts with the tensile load), local geometrical stress concentration (e.g., due to local reduction in cross section), and the presence of hard and brittle intermetallic particles.

The contour plots of void volume fraction computed in the FE simulations with the homogeneous and the heterogeneous GTN models are illustrated in Figure 12c,d. Figure 12d with the heterogeneous model shows, for example, a similar VVF to the homogeneous model (Figure 12c) in strut number 5 at the place of fracture, despite the presence of a cluster of particles that was not taken into account in the homogeneous model. Strut number 2 illustrates that similar levels of VVF could also be reached in the absence of particles due to geometrical plastic strain localization, which was predicted in both homogeneous and heterogeneous models. Strut number 1 illustrates that such localization of plastic flow due to local cross section reduction can be the critical point even in the presence of particles, in which case both models again provided comparable local VVF.

**Figure 12.** (**a**) Distribution of intermetallic particles in the solid phase of the 20 PPI sample. (**b**) Fractured struts of the foam. FE contour plots (at 6.4% strain) of VVF with the (**c**) homogeneous and (**d**) heterogeneous GTN models. Corresponding contour plots of equivalent stress with the (**e**) homogeneous and (**f**) heterogeneous GTN models. Struts of interest are labeled 1 to 5.

Similar comparisons can be made in terms of equivalent stress distribution in the struts in Figure 12e,f). Stress concentration arose either from local reduction of cross section, in which case both models provided comparable results, or with an additional contribution of the local presence of hard particles. One can note, however, that the sole prediction of stress concentration would fail in systematically discriminating the potential fracture zones at the microscopic level, justifying the use of the GTN porous plasticity model in the present work.

This study illustrates that stress analysis in metal foams is not straightforward due to the non-uniform stress distribution arising from the non-uniform geometry [45] and microstructure. The distribution of von Mises stress in the foam under loading depends on size, orientation, and spatial arrangemen<sup>t</sup> of the cells in a complex manner [17], in addition to the local presence of hard particles. It was however observed here that tensile loading of the foams promoted stress concentration in the struts of the biggest void-cells. The subsequent fracture of these struts formed the fracture plane.

Other studies have also suggested that the fracture mode of ERG foam struts depends on the type of precipitates in the Al matrix. The failure of the struts was observed to begin by realignment and ductile transgranular fracture of struts in the fracture plane [25]. The fracture mode later shifted to major brittle intergranular and minor ductile transgranular failure of the remaining struts due to the presence of α-AlFeSi (Al8Fe2Si) and β-AlFeSi (Al5FeSi) precipitates in the grain boundaries. The experimental protocol of the present study using laboratory tomography does not allow us to easily distinguish these different types of precipitates. The use of in situ synchrotron tomography might be an interesting perspective in this context.

As concerns modeling, current limits of the present approach might be overcome in the future by: (i) using a modified FE simulation environment allowing to alleviate the requirement for first-order tetrahedral elements in highly non-linear simulations with extensive damage development, (ii) using more realistic boundary conditions (e.g., directly mapped from the 3D in situ images), (iii) better identifying the local strut constitutive behavior. Indeed, in this study and in the recent one of Petit et al. [20], direct measures of yield stresses and hardening behavior on single-strut micro-tensile tests [41,43] had to be adjusted to properly reproduce the macroscopic foam response. A perspective might consist of taking the crystallographic orientation of the grains in the struts into account for the calculation of their plastic behavior. This could have significant consequences in this material, since some struts are made of a limited number of grains [21], and could exhibit rather anisotropic plastic behavior. Image-based crystal plasticity FE simulations [46] might be particularly interesting in this context.

As a final comment, the present study shows that *quantitatively* taking the local intermetallic particles into account in the prediction of the foam mechanical behavior does not really affect the macroscopic response, but allows a rather good discrimination of the critical zones for fracture in the structure. This confirms the *qualitative* conclusions and prospects of Petit et al. [20]. Nevertheless, it also tends to indicate in this particular case that accounting for intermetallics in such a homogenized GTN damage simulation framework is only of secondary importance for the prediction of fracture, after the proper accounting for local geometry and loading mode of the nodes/struts. Different conclusions were recently found in the study of Amani et al., [40] where the presence of local heterogeneous microporosity in additively manufactured lattice structures had more significant consequences on the prediction of fracture behavior, ye<sup>t</sup> with a higher global volume fraction and a different type of defect (process-induced cavities) as compared to the present study. Note also that direct full-field modeling of the small intermetallics in the struts would allow a more accurate prediction of local stress heterogeneity and fracture onset. However, such full-field simulations would require a severe increase in computational resources and would hardly be applicable on real-size samples. Moreover, the identification of the local constitutive behavior in the struts remains a challenge. The micromechanical approach developed in this study is based on cell-calculation and the GTN model to derive the constitutive behavior of local particle-rich regions of the foam. It can constitute the basis for a wider range of applications to complement 3D image-based FE studies of macroporous materials, which are now well established. This might include, for example, any multiscale investigation of architectured materials where both (i) a global description of the macroscopic structure, but also (ii) local image-based microstructural details are expected to be required to understand and predict the macroscopic behavior of the material.
