A micro-mechanics based extension of the GTN continuum model accounting for random void distributions

Randomness in the void distribution within a ductile metal complicates quantitative modeling of damage following the void growth to coalescence failure process. Though the sequence of micro-mechanisms leading to ductile failure is known from unit cell models, often based on assumptions of a regular distribution of voids, the effect of randomness remains a challenge. In the present work, mesoscale unit cell models, each containing an ensemble of four voids of equal size that are randomly distributed, are used to find statistical effects on the yield surface of the homogenized material. A yield locus is found based on a mean yield surface and a standard deviation of yield points obtained from 15 realizations of the four-void unit cells. It is found that the classical GTN model very closely agrees with the mean of the yield points extracted from the unit cell calculations with random void distributions, while the standard deviation S varies with the imposed stress state. It is shown that the standard deviation is nearly zero for stress triaxialities T ≤ 1 / 3 , while it rapidly increases for triaxialities above T ≈ 1 , reaching maximum values of about S /σ 0 ≈ 0 . 1 at T ≈ 4 . At even higher triaxialities it decreases slightly. The results indicate that the dependence of the standard deviation on the stress state follows from variations in the deformation mechanism since a well-correlated variation is found for the volume fraction of the unit cell that deforms plastically at yield. Thus, the random void distribution activates different complex localization mechanisms at high stress triaxialities that differ from the ligament thinning mechanism forming the basis for the classical GTN model. A


Introduction
Ductile failure processes lead to a loss of load-carrying capacity and occur through either localization of plastic flow (Tekoğlu et al., 2015;Guo and Wong, 2018;Liu et al., 2019) or by macroscopic yielding (Hure, 2021).The complex sequence of micromechanisms controlling failure has been studied for decades (Huang and Hutchinson, 1989;Tvergaard, 1990;Thomson et al., 1998;Benzerga and Leblond, 2010) and it is widely accepted that ductile damage predictions in component-sized structures require continuum modeling for computational efficiency.Thus, numerous homogenized yield criteria have been developed for ductile failure modeling, with the most famous being the model by Gurson (1977), which explicitly accounts for the porosity, f , of the material.However, the Gurson model was early on adjusted by Tvergaard (1981) and Tvergaard and Needleman (1984) to improve the model accuracy (through the Tvergaard-constants q 1 and q 2 ) and account for accelerated void growth at microscopic localization (through the coalescence model imposed through a critical porosity f * ).This widely used combined model is nowadays known as the Gurson-Tvergaard-Needleman (GTN) model.However, many other extensions have been proposed over the years (see Benzerga and Leblond, 2010, and references therein), and while most focuses on obtaining homogenized material properties, the present study will quantify the effect of random void distributions and illustrate how the variation on meso-scale properties can be included in a modified GTN model.The choice of the GTN model is grounded on its appropriateness for spherical pores.However, the framework can be readily applied to other models, and this is particularly relevant when investigating localization, where the GTN model is limited relative to other models (Thomason, 1985; J o u r n a l P r e -p r o o f Journal Pre-proof Pardoen and Hutchinson, 2000;Benzerga, 2002;Leblond and Mottet, 2008;Tekoglu et al., 2012), or size effects (Niordson, 2008;Nguyen et al., 2020).
At a scale comparable to the void spacing, the flow strength is highly sensitive to void distribution (Needleman and Kushner, 1990).Realistic void distributions are far from periodically arranged, and both Scanning Electron-microscopy and X-ray tomography experiments of voided materials have demonstrated either random or clustered configurations (Buffiere et al., 1999;Lecarme et al., 2014;Hannard et al., 2017) -a property not built into the usual continuum-based models.Despite this, the void distribution effects have proven important to the microscopic localization process (Dubensky and Koss, 1987;Magnusen et al., 1988) and the macroscopic fracture in ductile plate tearing (Tekoglu and Nielsen, 2019;Andersen et al., 2020;Çelik et al., 2021).Generally, a lower ductility is obtained for random distributions due to the triggering of local plastic yielding and loss of load-carrying capacity.For example, Becker (1987) performed finite element simulations of a 2D model for an elastic-viscoplastic material obeying the Gurson-Tvergaard yield criterion with an inhomogeneous distribution of voids.The results showed that plastic flow concentrates into bands in areas with large volume fractions of voids.Perrin and Leblond (1990) analyzed a composite sphere of two porous plastic materials, each obeying the Gurson criterion with individual porosities.The macroscopic yield stress was proven to be different from the one derived from the average homogeneous porosity under hydrostatic stress, indicating that the porosity distribution will affect the effective properties of voided materials.
Unit cell calculations with a single void representing a regular distribution of voids give important fundamental insight into yield properties and damage on a macroscopic scale.However, these material properties are also highly influenced by the spatial void distribution.Many studies have been carried out to investigate details of plastic flow and damage in materials with realistic void distributions.Fritzen et al. (2012) studied statistical effects in large ensembles of voids with volume fractions of up to 30%, leading to a proposed extension of the Gurson-Tvergaard model in terms of a volume fraction dependence of the correction parameters introduced by Tvergaard (1981).Khdir et al. (2015) demonstrated that the same model could be used for differently shaped voids assuming that the RVE is of sufficient size.Significant dispersion of failure strains, even for large RVEs, was reported through a numerical investigation in Cadet et al. (2021Cadet et al. ( , 2022)).While small RVEs with a single void may not provide a sufficiently realistic representation of macroscopic properties in terms of plasticity and damage, very large RVEs may yield results that do not adequately represent the microscopic variation in material properties on a small scale.On the other hand, intermediate RVEs with a few voids may be used to represent the spatial statistical variations on the scale of individual integration points in a numerical model, thus leading to an appropriate representation of spatially varying mesoscale properties.This work aims to understand and quantify the effect of the spatial void distributions in terms of the macroscopic yield stress and its dependence on stress triaxiality.
To achieve this, three-dimensional representative volume element calculations with periodic boundary conditions containing four spherical voids of equal size distributed randomly are carried out.The voids are embedded in an elastic-perfectly plastic material.Several randomizations are considered for constant initial void volume fraction and stress triaxiality to bring out the statistical variation of the yield locus.In this way, the aim is not to achieve a homogenized response by pursuing a sufficiently large unit cell but rather to understand how void distribution affects the dispersion around a mean yield locus.A subsequent statistical analysis gives input for a proposed yield surface accounting for the dispersion of the material yield point for different void configurations.A modified GTN model is proposed including statistical variations, and the model compares well to the unit cell simulation results.Finally, a procedure for implementing the new yield surface into a large-scale calculation is presented.
The paper is organized as follows: Section 2 presents the problem formulation in terms of a unit cell with random void distributions, the modeling approach, and the fundamental quantities for the discussion of results.The results from the unit cell study, alongside a comparison to the classical GTN model and a new extension, are presented and discussed in Section 3. The work is concluded in Section 4.

J o u r n a l P r e -p r o o f
Journal Pre-proof

Problem formulation and modeling approach
This work considers a limit load-type analysis of a porous metal with random distributions of discretely modeled microvoids to determine the statistical variation in the yield surface characteristics.Attention is on axisymmetric stress states in the full range of positive stress triaxiality, and the imposed condition is kept constant for each point evaluated on the yield surface (ensuring proportional loading).The matrix is modeled as a rate-independent, perfectly-plastic von Mises material (J 2 -flow theory), and a small strain finite element formulation is adopted to mimic the limit load of the material.The matrix material is characterized by the parameters: σ 0 /E = 0.001 and ν = 0.3, where σ 0 is the yield stress, E is Young's modulus, and ν is the Poisson ratio.The unit cell configuration is described in Section 2.1, the modeling approach is outlined in Section 2.2, while the fundamental quantities in the statistical analysis and for the discussion of results are outlined in Sections 2.3 and 2.4.

Unit cell configuration
Figure 1 shows the cubic unit cell setup, where a 0 denotes the side length along the three coordinate axes x i (i = 1, 2, 3).Each unit cell contains four (N = 4) spherical voids defined by their center coordinates and initial radius r 0 .Here, the non-dimensional total porosity f 0 = 4π 3 N r 3 0 V (or unit cell void volume fraction) of the unit cell with volume V determines the void radius, and the results are presented for three different f 0 -values in Section 3. The spatial location of the voids is determined using an ad-hoc algorithm implemented by means of the Abaqus2Matlab software (see Papazafeiropoulos et al., 2017).The algorithm generates a given number of 3D spheres inside a 3D domain, with the radii and the positions of the spheres being uniformly random.The algorithm, simplified to equally sized voids, involves the following steps: (i) A 3D point grid of potential void centers is created within the unit cell.
(ii) The order of the grid points is randomly permuted.
(iii) The mutual distances between all center points are calculated and a new center point is defined if the distance is smaller than a minimum distance 2r 0 + L, with J o u r n a l P r e -p r o o f Journal Pre-proof r 0 and L being the void radius and minimum ligament size, respectively.
In this way, the algorithm introduces a set of non-overlapping spherical voids, and a representation of three random configurations is shown in Fig. 2 in terms of finite element meshes.
For comparison purposes, the present work also considers a regular Face-Centered Configuration (FCC) void distribution loaded along the cubic axes.Figure 3a illustrates how the configuration can be modeled considering only four voids when translating the unit cell in the positive x 1 -, x 2 -and x 3 -direction.The corresponding finite element mesh for the FCC unit cell is shown in Fig. 3b.It is worth noticing that the FCC unit cell response is independent of the translation along the coordinate axes.

Numerical modeling approach
The present work adopts a small strain finite element formulation and exploits the commercial software package Abaqus (2020).Thus, the model setup cannot account for the softening owing to void growth.Instead, the load-carrying capacity of the unit cell represents the limit load and, thereby, a point on the yield surface when plotted in stress space.Throughout, axisymmetric stress states with σ 2 = σ 3 are considered such that the stress state is defined by the stress ratio where σ 1 is the stress along the main loading axis, and ρ is kept constant and prescribed for each individual analysis of a point on the yield surface.Thus, the von Mises equivalent stress σ e , and the overall mean stress σ m are while the stress triaxiality T is related to the stress ratio through In the finite element calculations, the ratio ρ between the transverse and axial stress components is kept constant using multiple-point constraints (MPCs in Abaqus).This

J o u r n a l P r e -p r o o f
Journal Pre-proof is achieved by introducing an extra set of degrees of freedom to the finite element mesh in terms of three dummy nodes N i (i = 1, 2, 3) placed outside the finite element mesh as depicted in Fig. 4. The dummy nodes are connected with spring elements (SPRING2 in the Abaqus) to three master nodes M i (i = 1, 2, 3), which are part of the unit cell mesh.In this way, the displacement of the dummy nodes is related to the forces F i acting on the faces of the unit cell (along its normal) through where k i are the spring element constants.Moreover, the forces F i are related to the macroscopic stresses through where A i are the surface areas of the unit cell.Thus, combining Eqs. ( 1), (4), and ( 5), and solving for the displacement of the dummy nodes (for constant ρ) gives u Here, u N 1 1 is the prescribed displacement in the main loading direction x 1 , while the remaining displacements of the dummy nodes are calculated in the MPC subroutine.
Furthermore, to ensure the periodicity of the unit cell, a set of linear constraint equations is imposed according to Faces: Edges: Vertices: Here, u i is the displacement in the ith direction (i = 1, 2, 3) of the nodes at the faces, edges, and vertices of the unit cell, as shown in Fig. 4. The nodes at the corners B ′ , C, and D ′ are the master nodes.

J o u r n a l P r e -p r o o f
Journal Pre-proof

Statistical variation
The random void distributions naturally introduce a statistical variation in the model output.The yield point extracted at the limit load of the unit cell depends on the localization process between voids (or lack thereof) and, thus, on the intervoid distance and the location of the voids.Thus, a mean µ and standard deviation S are introduced to characterize the span of yield points obtained for a specific stress state and initial porosity.Note that proportional loading is imposed using a number of constant stress ratios ρ, such that all the determined yield points for a given f 0 and ρ will be located on a straight line through the origin in the σ e − σ m stress space.Let s i = σ 2 e + σ m 2 be the distance from the origin to the yield point of the ith unit cell calculation for a specific value of ρ.The mean of n unit cell calculations with different random void distributions is then given by and the corresponding standard deviation of the mean distance is

Characterisation of deformation mechanism
To discuss the mechanism of plastic deformation leading to the loss of load-carrying capacity of the unit cell, a plastic index PI is introduced as where V m is the volume of the matrix material, and V p is the volume of the unit cell undergoing plastic yielding.The plastic index provides a way to distinguish between macroscopic yielding and localization of plastic flow (or homogeneous versus inhomogeneous yielding in the terminology of Hure, 2021).At the limit load, the major part of the unit cell undergoes macroscopic yielding when PI → 1, while small values of PI indicate localization in part of the unit cell volume.

J o u r n a l P r e -p r o o f
Journal Pre-proof

Numerical results and discussion
In the following, yield surfaces are constructed in the von Mises versus mean stress space for three values of the total porosity f 0 = 0.00085, 0.017, and 0.034 to investigate the influence of the void distribution.The results are obtained by imposing stress states corresponding to nine different values of the stress ratio ρ = −0.5, 0, 0.4, 0.625, 0.73, 0.8, 0.85, 0.9, and 0.99, spanning the range of positive stress triaxialities (T = 0, 0.33, 1, 2, 3, 4.33, 6, 9.33, 99.33).The purely hydrostatic state of stress ρ = 1 is here omitted due to convergence issues.The calculations for each stress state and porosity are repeated using 15 different randomizations of the void distribution to form a statistical basis for the discussion of results.In addition, the FCC distribution is investigated.Thus, a total of 432 combinations of stress state, void distributions, and void volume fraction are considered.
3.1.Yield surfaces extracted as the unit cell limit-load show results for unit cells with random void distributions, while the square markers are the mean value given by Eq. ( 8) for a specific combination of porosity and triaxiality.
Here, solid lines indicate the yield surface represented by the mean values.In addition, triangular markers indicate the corresponding results for the unit cell with a regular FCC distribution, while the overlaying dashed lines illustrate the corresponding yield surfaces.Colors distinguish the results for the different initial porosities f 0 , and the well known delay in yielding is evident for diminishing initial porosity.This feature is also represented in the classical Gurson-Tvergaard-Needleman (GTN) yield surface (Gurson, 1977;Tvergaard, 1981;Tvergaard and Needleman, 1984) which is defined by: Here, σ e is the von Mises stress, σ 0 is the matrix material flow stress, σ m = σ kk /3 is the mean stress, f is the void volume fraction, and q 1 = 1.5 and q 2 = 1 are the Tvergaardconstants (Tvergaard, 1981).As both the von Mises equivalent stress and the mean  Comparing results from the random void distributions to that of the FCC unit cells, the regular void distribution generally displays more plastic resistance, although it is not a strict upper bound.Evidently, the yield point for the FCC configuration falls below that of some of the random distributions for triaxiality values in the range T = 1 to 4, where a large dispersion in the yield points is observed for the random distributions (see Fig. 5).At this triaxiality level, the stresses transverse to the main loading axis increase the propensity to alter the localization mode as the intervoid distances vary across the unit cell.In contrast, the regular FCC configuration has the same ligament geometry between all voids and, thus, will not exhibit the same shift in localization mechanism when changing the stress state.However, the mean for the random void distributions (solid lines) consistently gives less plastic resistance than the FCC configuration (dashed lines), independently of the initial porosity.The differ-

J o u r n a l P r e -p r o o f
Journal Pre-proof ence is prominent at moderate to high triaxialities, while the yield surfaces practically coincide at low stress triaxialities.
The dispersion in the yield points for the random void distributions may be quantified by the standard deviation S of the distance to the origin (see Eq. ( 9)).Here, results are based on 15 randomizations of each combination of porosity and triaxiality.
Figure 6 is constructed from the unit cell results to show the yield loci defined by the mean surface plus minus the standard deviation (µ ± S) for the three investigated void volume fractions.It is seen that the standard deviation increases along the mean stress axis, i.e., with increasing triaxiality as discussed for the dispersion concerning Fig. 5.
That is, the values of S are small for low stress triaxialities approaching zero as T → 0, while values of S are the highest in the interval T = 4 to 5. Above this triaxiality level, it decreases slightly as the triaxiality goes to infinity (ρ → 1).The span of the yield loci in Fig. 6 is consistent with the span with the yield points quantified out in Fig. 5.
From Fig. 6, it also becomes clear that the dispersion of the yield point at intermediate levels of stress triaxiality and the narrowing as ρ → 1 is most prominent for high initial void volume fractions.This can be ascribed to the change between localization mechanisms.For high f 0 -values, the intervoid ligaments carry higher stresses due to the smaller volume fraction of matrix material, making the ligaments more susceptible to the deformation mechanism involving localization.
To quantify the variation in the dispersion of yield points, the standard deviation normalized by the yield stress S/σ 0 is shown as a function of ρ in Fig. 7. Values obtained from the unit cell calculations are circular markers, while a cubic Hermite interpolation of the results is a continuous dashed line.The results confirm the increase in standard deviation with increasing ρ until a peak is reached around ρ = 0.8 (corresponding to T = 4.3), after which the standard deviation decreases at higher triaxiality levels.The peak value and the subsequent drop in the standard deviation are largest for the two highest initial porosities (f 0 = 0.017 and 0.034), while the decrease at high triaxialities is more modest for the lowest initial porosity (f 0 = 0.0085).The variation in S/σ 0 with the prescribed stress state ρ is clearly reflected in the scatter of the results in Fig. 5.
Moreover, it is noticed from Fig. 7 that increasing f 0 leads to a larger value for S/σ 0 J o u r n a l P r e -p r o o f Journal Pre-proof for all triaxiality values.

Mechanisms leading to the dispersion of yield points
The plastic index PI introduced in Eq. ( 10) is considered in an attempt to link the dispersion of the yield points to the localization mechanism at play in unit cells with random void distributions.The plastic index is calculated at the limit load for all configurations of f 0 and T and displayed with circular markers in At high triaxialities, the dispersion in the PI is large, consistent with some RVEs having localized plastic deformation.This manifests itself in an inclined intersection of the yield surface with the mean stress axis (see Fig. 5), consistent with localization models like that of Thomason (1985).The classical GTN model cannot account for the dispersion of the yield points observed in Fig. 5 and mapped out in Fig. 7, which are consequences of randomness in the void distribution.However, as demonstrated in Fig. 6, the unit cell response is accurately represented by overlaying the mean distance µ by the standard deviation S suggesting that the yield locus is within

Introducing the effect of randomness into the GTN yield surface
where Φ µ is the mean of the yield locus and Φ S represents the spread of the yield surface.Thus, since the GTN model quite accurately models the mean yield surface (see Fig. 10), it is suggested to scale the distance to the origin of the GTN yield surface J o u r n a l P r e -p r o o f Journal Pre-proof with the standard deviation such that the distribution-enriched GTN yield locus, within a standard deviation, is expressed as Here, S is the standard deviation of the distance to the yield surface, which is a function of the porosity f and the stress state ρ according to Fig. 7.A continuous yield function is obtained for the stress states and porosity values considered by incorporating the cubic Hermite interpolation, and Fig. 11 depicts the new yield surfaces alongside the classical GTN yield surface.The GTN yield surface (dashed lines) represents the mean surface of the unit cell with random void distributions, while the dispersion of yield points is obtained through the scale factor (1 ± S/σ 0 ) 2 shown as solid lines on either side of the mean yield surface.The depicted confidence interval of ±S/σ 0 encloses 70% of the expected yield points observation due to random void distributions.It is worth noting that the characteristic dispersion of the yield surface is achieved such that there is little effect of the random distribution at low stress triaxiality, while the spread of the curves on either side of the mean increases with stress triaxiality until about T ≈ 4.
Moreover, the narrowing of the dispersion in the yield surface close to the mean stress axis is also captured (see Figs. 5,6,and 11).However, it should be noted that, in the high triaxiality regime the GTN model does not represent localization accurately and hence efforts to model yield surfaces analytically should be augmented by appropriate localization strategies.
The proposed distribution-enriched GTN yield surface in Eq. ( 13) can be employed in large-scale continuum modeling by assigning individual finite elements (or Gauss points) a new material value α, which determines the yield surface for this particular material point.That is, the distribution-enriched yield surface may be expressed as where the α-value could be assigned with a random spatial distribution of average value zero throughout the volume, such that it follows a normal distribution of the form J o u r n a l P r e -p r o o f

Journal Pre-proof
Here, S is the corresponding standard deviation (see Fig. 7).It would be expected that the interval α ∈ [±2S] would represent about 95% of the yield point observations obtained from unit cell calculations, while 70% of the observations would lie in the interval α ∈ [±S] (depicted in Fig. 6).Note that for a given value of the distance to the yield surface, appropriate derivatives of α must be included to obtain the yield surface normal.Adopting this procedure and introducing the yield surface from Eq. ( 14) into a continuum-based finite element calculation would reflect the dispersion of yield points when using a discretization where individual integration points represent about four voids.Increasing the number of voids described by a yield surface in an integration point would need appropriate scaling of the standard deviation up until the limit, where each integration point describes a very large ensemble of voids and a mean yield surface is appropriate.In this way, the resolution of the discretization will control the effect of the random void distributions.
In this work, unit cells with four voids have been chosen as representative of a microstructural configuration.This was deemed to be a good balance between the highly idealized single-void unit cell and a unit cell with a large number of voids which can be computationally costly and not representative of inhomogeneous deformation at a scale of a characteristic void distribution.The framework can readily accommodate unit cells with different number of voids.A larger number of voids would represent a larger material region and thus lead to smaller standard deviation of the unit cell properties.

Conclusions
The present work demonstrates how the dispersion of yield points, identified as the limit load, in ductile metals with random void distributions ties to the deformation • A strong dependency on the stress state exists.The dispersion due to random void distributions is practically zero at low triaxiality, while it grows to the largest value in the range of 4 < T < 5 and drops slightly for higher triaxialities (see Figs. 5 and 6).The reason is found in the deformation mechanism at play as macroscopic yielding prevails at low triaxiality (in line with Tekoğlu et al., 2015), while a complex mixture of localization modes can develop at higher triaxiality depending on the intervoid distance and location of the voids (see discussion in Section 3.3).
• The standard deviation follows the deformation mechanism.The plastic index introduced in Eq. ( 10) is adopted to demonstrate a correlation between the dispersion and the prevailing deformation mechanism.The index equals one when plastic straining occurs in the entire unit cell at the limit load, while small index values signal localization in a portion of the unit cell volume.It is found that the index displays a large spread for different randomization when T ≳ 4 and that the index standard deviation of the observations correlates with that of the dispersion of the yield points.Thus, it is concluded that the variations in the localization mechanism determine the dispersion of the yield points.
• The porosity influences the spread of the yield locus.The variation in the standard deviation of both PI and µ depends on the porosity as the intervoid ligament size diminishes with increasing void volume fraction.The variation is most significant for a large porosity such that the peak value of S attains the highest level at ρ ≈ 0.8 (corresponding to T ≈ 4.3) but also the largest relative drop at higher triaxialities (see Fig. 7).Thus, the width of the yield locus up until T ≈ 4.

Figure 5
Figure 5 presents the simulated yield points for all combinations of porosity f 0 , stress triaxiality T , and void distribution considered in the present work.The circular markers

J
o u r n a l P r e -p r o o f Journal Pre-proof stress enter explicitly into the GTN yield surface, it can readily be represented in the von Mises versus mean stress space.

Figure 5
Figure5shows a combined dependence on the void distribution and the stress triaxiality.For the lowest triaxiality values, i.e., the results closest to the von Mises stress axis, the effect of the void distribution is negligible, and all yield points practically coincide.This is in line with results presented inTekoğlu et al. (2015) andHolte et al. (2021), where T = 1, corresponding to ρ = 0.4, is found as the limit below which the onset of macroscopic yielding (homogeneous yielding according to Hure, 2021) co-occurs with void coalescence, i.e., with intervoid localization.However, the dispersion in the yield points amplifies when increasing the stress triaxiality indicating the activation of different or more complex localization mechanisms that depends on the interaction of voids and, thereby, their spatial distribution.The results for a specific porosity and stress state show yield points over a significant range of σ e and σ m combinations, indicating that the void distribution is an important microstructural property for porous metals at these triaxiality levels.Moreover, for the highest triaxiality (ρ ≈ 1), i.e. the results closest to the mean stress axis, the spread in yield points is moderate compared to the slightly lower triaxiality levels.This indicates a shift in the localization mechanism, leading to a smaller dependency on the void distribution.
Fig. 8.The mean value of the plastic index µ PI for each porosity value is shown as a function of ρ (solid lines).A plastic index of 1 corresponds to yielding in the entire unit cell, while low PI-values signal intense localization in a smaller part of the unit cell volume.It is observed from Fig.8that almost the entire unit cell volume deforms plastically for ρ ≲ 0.5 (corresponding to T ≲ 4/3) for all values of f 0 considered.This is well in line with the fact that macroscopic yielding (or homogeneous yielding) is the dominant deformation mechanism at low to moderate values of stress triaxiality, rendering the effect of the void distribution negligible.In contrast, the plastic index shows a much greater dispersion for higher values of ρ, ranging from 0.4 to 1.The low values of the plastic index reflect intense localization in a small portion of the unit cell, which is highly controlled by the location of the voids and, thereby, the void distribution.At high triaxiality levels, the limit load can be attained through macroscopic yielding, e.g., for equally distanced voids, and microscopic localization if voids are located in a favorable band.Moreover, the large stress components transverse to the main loading direction increases the likelihood of encountering a favorable voided band.It is this combined effect of different mechanisms or between different favorable localization modes at high stress triaxialities that makes void distribution essential to the plastic resistance.In this way, it is the localization mechanism that controls the dispersion in both the plastic index PI and the yield loci shown in Figs.5 and 6.This is also evident from Fig.9, showing the standard deviation of the plastic index S PI as a function of ρ for all values of f 0 considered.For increasing triaxiality, i.e., increasing ρ, the standard deviation increases, indicating a greater variation in the plastic deformation mechanism, and the variation compares well to that of S/σ 0 in Fig.7.As seen in Fig.7, J o u r n a l P r e -p r o o f Journal Pre-proof the standard deviation for the plastic index S PI also increases with increased f 0 since the size and position of the intervoid ligaments affect the deformation mechanism to a greater extent when the void volume increases.

Figure 10
Figure10shows the mean distance to the origin µ based on Eq. (8) as obtained from the unit cell investigations of random void distributions where each data point and corresponding standard error is based on 15 different randomizations of the void distributions.The small error bars in Fig.10show the standard error of the mean, indicating that a sufficient sample size is used.Moreover, the initial void volume fraction has little effect on the plastic resistance distance for T < 1, while the influence of porosity increases with triaxiality.Comparing the results to the predictions by the classical GTN model, a good agreement is obtained for all values of the initial porosity f 0 and stress states ρ.Thus, the GTN model may form a basis for a micro-mechanics based continuum model accounting for random void distributions if modified suitably.
mechanism and suggests a way to incorporate the findings into the Gurson-Tvergaard-Needleman yield surface.The study relies on a numerical investigation of a periodic microstructure represented by unit cells containing a limited number of randomly distributed voids.The unit cell setup is considered a mesoscale model of the material and J o u r n a l P r e -p r o o f Journal Pre-proof it provided insight into the statistical characteristics of the yield locus owing to the void distribution.The key findings for the dispersion of the yield points are

FiguresFigure 1 :
Figs. 6 and 11).The present work investigates statistical variations of yield surfaces for porous materials.It is shown that the classical GTN yield surface rather accurately models the

Figure 2 :Figure 3 :
Figure 2: Examples of three unit cell finite element meshes containing different randomizations of voids.

Figure 4 :
Figure 4: Modeling procedure for imposing periodic boundary conditions and loading to the unit cell, showing the faces, edges, and vertices defined by A, B, C, D and A ′ , B ′ , C, D ′ .The dummy nodes to control the imposed stress state are denoted N i , and the related master nodes are M i (i = 1, 2, 3), which are part of the finite element mesh.The dummy and master nodes are connected with spring elements, as illustrated.

Figure 5 :
Figure 5: Yield points for different combinations of porosity f 0 , stress state T , and void distribution (15 different randomizations of each f 0 -T -combination), alongside the calculated mean values (solid lines) and results for the FCC unit cell (dashed line).The stress state is limited to axisymmetry and controlled by the ρ = σ 2 /σ 1 = σ 3 /σ 1 as outlined in Section 2.2.Colors distinguish the different porosity values.

Figure 6 :
Figure 6: Yield surfaces constructed from the unit cell results based on the mean distance µ from the origin in the von Mises versus mean stress space, and the corresponding standard deviation S/σ 0 (see Section 2.3).The grey-shaded region encloses about 70% of the yield points.Colors distinguish the different porosity values.

Figure 7 :
Figure 7: Standard deviation S/σ 0 of the mean distance from the origin to the yield surface normalized by the yield stress σ 0 as a function of the applied stress ratio ρ for the three porosity-values f 0 considered.The dashed curves represent a piece-wise cubic Hermite interpolation of the discrete data points.Colors distinguish the different porosity values.

Figure 8 :
Figure 8: The plastic index PI introduced in Eq.(10) for all combinations of porosity f 0 , stress state T , and void distributions (15 different randomizations of each f 0 -T -combination) as a function of the ratio ρ.The mean plastic index µ PI is also shown as solid lines.Colors distinguish the different porosity values.

Figure 9 :
Figure 9: Standard deviation S PI of the plastic index PI (see Eq. 10) as a function of the applied stress ratio ρ for the three porosity-values f 0 considered.The dashed curves represent a piece-wise cubic Hermite interpolation of the discrete data points.Colors distinguish the different porosity values.

Figure 10 :
Figure 10: Mean distance from origin µ form the unit cell calculations as a function of the stress ratio ρ, alongside with error bars indicating the standard error of mean SEM = S/ √ n.Here is S the standard deviation, and n is the sample size (n = 15).The results are shown together with the predictions for the classical Gurson-Tvergaard-Needleman (GTN) model for the three values of f 0 considered.Colors distinguish the different porosity values.

Figure 11 :
Figure 11: Yield surfaces constructed from the new distribution-enriched Gurson-Tvergaard-Needleman (GTN) model in Eq. (13), showing the mean curve as a dashed line (repressed by the classical GTN model) and the dispersion of the yield surface represented by a confidence interval of ±S/σ 0 (enclosing 70% of the yield point observations).The yield surfaces are illustrated for the three initial porosity values considered distinguished by different colors.