The objective is to identify the micromechanism(s) of ductile crack advance, and isolate the key microstructural and material parameters that a ect these micromechanisms and fracture toughness of ductile structural materials. Three dimensional, nite element, nite deformation, small scale yielding calculations of mode I crack growth are carried out for ductile material matrix containing two populations of void nucleating particles using an elasto-viscoplastic constitutive framework for progressively cavitating solid. The larger par- ticles or inclusions that result in void nucleation at an early stage are modeled discretely while smaller particles that require large strains to nucleate voids are homogeneously dis- tributed. The size, spacing and volume fraction of inclusions introduce microstructure-based length-scales. In the calculations, ductile crack growth is computed and fracture toughness is characterized. Several features of crack growth behavior and dependence of fracture tough- ness on microstructural and material parameters observed in experiments, naturally emerge in our calculations. The extent to which the microstructural and material parameters a ect the micromechanisms of ductile crack advance and, hence, the macroscopic fracture tough- ness of the material is discussed. The results presented provide guidelines for microstructural engineering to increase ductile fracture toughness, for example, the results show that for a material with small inclusions, increasing the mean inclusion spacing has a greater e ect on fracture toughness than for a material with large inclusions. Keywords: Fracture mechanism, Fracture toughness, Voids and inclusions, Finite elements, Ductile fracture 1. Introduction The objective of this work is to correlate the macroscopic fracture toughness of heteroge- neous ductile materials with the fracture mechanisms operating at the microscale. The focus is con ned to ductile structural materials undergoing fracture due to nucleation, growth and coalescence of microscale voids (Tipper, 1949; Goods and Brown, 1979; Tvergaard, 1989;  Email address: ankit.sri@tamu.edu (Ankit Srivastava) Preprint submitted to JMPS March 23, 2020 Benzerga and Leblond, 2010; Benzerga et al., 2016). The microstructure of such materi- als can be idealized as a ductile material matrix with randomly distributed inclusions or second phase particles, as is encountered in a variety of ferrous and aluminum alloys. The second phase particles in such materials are either added intentionally to the material's microstructure, for precipitation-hardening or dispersion-strengthening (Martin, 2012), or their formation is an unavoidable consequence of processing (Birat, 2016). There is an important di erence between the characterization of fracture in a deforma- tion eld that is more or less uniform and the characterization of fracture in a deformation eld that is heterogeneous (Rice and Johnson, 1970; Rice, 1976; McMeeking, 1977; Aravas and McMeeking, 1985a; Needleman and Tvergaard, 1987; Gao et al., 1998; Tang et al., 2013; Srivastava et al., 2014). In the former, an unstructured continuum description of duc- tile fracture, in principle, can be based on material parameters such as strength and strain hardening exponent. On the other hand, in the latter, the macroscopic fracture toughness must involve characteristic length-scales associated with the material microstructure. This expected dependence of fracture toughness on material length-scales is consistent with the experimental results of ductile fracture due to nucleation, growth and coalescence of mi- croscale voids in ductile structural materials (Edelson and Baldwin Jr, 1959; Kra t, 1964; Birkle et al., 1966; Spitzig, 1969; Hahn et al., 1972; Passoja and Hill, 1974; Green and Knott, 1976; Lautridou and Pineau, 1981; Garrison and Moody, 1987; Dubensky and Koss, 1987; Becker, 1987; Garrison et al., 1997; Benzerga et al., 1999; Jie et al., 2007; Garrison Jr and Wojcieszynski, 2009). Several attempts have been made in the past to correlate material microstructural param- eters such as volume fraction, size and spacing of the inclusions/voids with the macroscopic fracture toughness of the material. The mode I plane strain fracture toughness, KIC, of several ductile materials has been shown to follow KIC=EN(2L)1 2(Kra t, 1964) or KICp20E(Hahn and Rosen eld, 1975), where Eis Young's modulus, Nis the work hardening coecient, 0is the ow strength, Lis the size of the process zone and is the spacing of the nearest void/inclusion from the initial crack tip. Often the process zone size has been hypothesized to correlate with the microstructural length-scale set by the pre-existing inclusions in the material (Birkle et al., 1966; Spitzig, 1969; Passoja and Hill, 1974). These suggest that the fracture toughness of the material predominantly depends on a single length-scale introduced by the mean (or extreme) spacing of the inclusions/voids. On the other hand, the fracture toughness of the material has been suggested to depend not only on the mean spacing but also on the size of the inclusions/voids. For example, the material's fracture toughness de ned in terms of JIChas been suggested to correlate with two microstructural length-scales as, JIC/0 2=Ni1=3 ln (rc=r0), whereNiis the number of inclusions per unit volume (1 =N1=3 iis the mean inclusion spacing, l0),r0is the initial size of the inclusions and rcis the radius of the cavities at failure nucleated from these inclu- sions (Amar and Pineau, 1985). Dubensky and Koss (1987) carried out an experimental study to distinguish the in uence of overall volume fraction, spacing and size of pre-existing voids/holes on material's ductility. Their experimental observation suggests that the spacing and size of the pre-existing voids has a larger e ect on ductility compared to their volume 2 fraction. The dominant e ect of the spacing and/or size of the pre-existing inclusions/voids over their volume fraction on macroscopic fracture toughness has also been suggested by Benzerga et al. (1999); Van Stone and Psioda (1975); Shabrov and Needleman (2002). The extent to which fracture toughness of a ductile material depends on the overall microstructural parameter i.e. volume fraction, or spacing and size of pre-existing voids has been associated with the micromechanisms of ductile crack advance. In general, two micromechanisms of ductile crack advance in a ductile material with preexisting distribution of voids have been proposed: crack advance by a void-by-void micromechanism and crack advance by multiple-void interaction (Tvergaard and Hutchinson, 2002; Kim et al., 2003; Pardoen and Hutchinson, 2003; Gao et al., 2005; Kim et al., 2007; Tvergaard, 2007; Teko glu et al., 2015). It has been suggested that for void-by-void micromechanism, JIC=(0l0) or c=l0(cbeing the critical crack-tip opening displacement) depends on the initial void volume fraction whereas for multiple-void interaction micromechanism, JIC=(0l0) orc=l0does not depend on the initial void volume fraction. At the microscale, crack advance by the void- by-void micromechanism involves coalescence of the crack tip with the closest void (Rice and Johnson, 1970; Aravas and McMeeking, 1985a), whereas crack advance by multiple- void interaction involves simultaneous interaction of multiple voids ahead of the crack tip (Tvergaard and Hutchinson, 2002). The transition from the crack advance by the void-by- void micromechanism to multiple-void interaction has been shown to depend on the initial void volume fraction (Tvergaard and Hutchinson, 2002; Kim et al., 2003), wherein void- by-void micromechanism is operative at low initial void volume fraction and multiple-void interaction micromechanism is operative at high initial void volume fraction. However, Gao et al. (2005), suggested that other factors such as initial distribution of the voids can also a ect the transition from void-by-void to multiple-void interaction micromechanism of ductile crack advance. There is a distinct advantage of modeling discrete voids in a ductile matrix in the nite element calculations, e.g. Tvergaard and Hutchinson (2002); Gao et al. (2005), as this ap- proach directly models the void growth behavior in the vicinity of a crack tip. However, it doesn't not naturally account for ductile crack advance. In addition, the class of ductile structural materials considered here, often contains a three-dimensional distribution of two populations of void nucleating particles: larger particles or inclusions (e.g. MnS inclusions in steels) that nucleate voids at relatively small strains and smaller particles (e.g. carbides in steels) that nucleate voids at larger strains. The voids nucleate either by debonding or cracking of the particles (Puttick, 1959; Rogers, 1960; Gurland and Plateau, 1963). Further- more, the three-dimensional distribution of inclusions in a ductile matrix can promote crack meandering (Needleman and Tvergaard, 1987; Tang et al., 2013; Srivastava et al., 2014) and microcracking (Tinet et al., 2004; Srivastava et al., 2014) ahead of the initial crack tip well before the overall crack growth. Nevertheless, in Srivastava et al. (2014) it was shown that for a xed initial inclusion size, JIC=(0l0) depends on the inclusion volume fraction for small inclusion volume fractions whereas for larger inclusion volume fractions JIC=(0l0) does not depend on the inclusion volume fraction. This raises several fundamental questions associated with inclusion driven ductile fracture: (i) Is the dependence of JIC=(0l0) on inclusion volume fractions is due to the transition from void-by-void to multiple-void inter- 3 action micromechanism of ductile crack advance as has been suggested for ductile materials with pre-existing voids? (ii) How does the interlacing of two microstructural length-scales, inclusion spacing and size, a ect the micromechanisms of ductile crack advance? (iii) How does the matrix material properties such as, strain hardening exponent, susceptibility to secondary void nucleation and energy dissipated in the growth of the nucleated voids prior to crack advance a ect the micromechanisms of ductile crack advance? To this end, nite element, nite deformation calculations are carried out using a consti- tutive framework for progressively cavitating ductile solid. The matrix material is modeled as an isotropic hardening elasto-viscoplastic solid together with two population of void nu- cleating particles. Larger particles or inclusions that result in void nucleation at an early stage are modeled discretely while smaller particles that require large strains to nucleate voids are homogeneously distributed. The size, spacing and volume fraction of inclusions introduce microstructure-based characteristic length-scales into the formulation. The cal- culations are carried out for small scale yielding conditions under remote Mode I loading. In the calculations, ductile crack growth are computed to the sucient extent to be able to characterize the fracture toughness of the material in terms of JICusing a procedure mimicking the ASTM:E1820-11 (2011). The extent to which the microstructural parame- ters such as volume fraction, size and spacing of the inclusions, and the properties of the matrix material a ect the micromechanisms of ductile crack advance and hence the macro- scopic fracture toughness of the material is discussed. The results presented also provide guidelines for microstructural engineering to increase fracture toughness. 2. Problem formulation and numerical method 2.1. Constitutive relation The constitutive framework used here is the modi ed Gurson elastic-viscoplastic consti- tutive relation for a progressively cavitating solid (Gurson, 1975; Tvergaard and Needleman, 1984; Chu and Needleman, 1980) with the ow potential having the form  =2 e 2+ 2q1fcosh3q2h 2 1(q1f)2= 0 (1) whereq1= 1:25,q2= 1:0 are parameters introduced in Tvergaard (1981, 1982b), fis the e ective void volume fraction,  is the matrix ow strength, and e2=3 20:0; h=1 3:I;0=hI (2) The function f, introduced in Tvergaard and Needleman (1984), is given by f=f; f <f c fc+ (1=q1fc)(ffc)=(fffc); ffc(3) where the values of fc= 0:12 andff= 0:25 are used. Calculations are also carried out to explore the e ect of the parameter fcon the fracture toughness of the material for xed microstructures. 4 The rate of deformation tensor is written as the sum of an elastic part, de=L1:^, a viscoplastic part, dp, and a part due to thermal straining, d= _I, so that d=L1:^+ _I+dp(4) Here, ^is the Jaumann rate of Cauchy stress,  is the temperature, = 1105/K is the thermal expansion coecient and Lis the tensor of isotropic elastic moduli. The plastic part of the strain rate, dp, is given by Pan et al. (1983) dp=" (1f)_ :@ @# @ @(5) The matrix plastic strain rate, _, is given by _= _0 g(;)1=m ; g(;) =0G() [1 + =0]N(6) with =R_dtand0=0=E. In the calculations the values of Young's modulus, E= 70GPa, Poisson's ratio, = 0:3, initial ow strength, 0= 300MPa, ( 0=0=E= 0:00429), strain hardening exponent, N= 0:1, strain rate sensitivity exponent, m= 0:01, and refer- ence strain rate, _ 0= 103s1, are used. The value of _ 0provides a reference strain rate so that for _= _0, =gin Eq. (6). The value of _ 0= 103s1is used here since it and other material parameters have the same values as employed in Srivastava et al. (2014); Osovski et al. (2015a); Srivastava et al. (2017). Parametric studies are also carried out to explore the e ect of the matrix strain hardening exponent, N, on the fracture toughness of the material for xed microstructures. Adiabatic conditions are assumed so that cp@ @t= :dp(7) withbeing the current density of the material (the initial density, 0= 7600kg/m3), cp= 465J/(kgK), = 0:9, and the temperature-dependence of the ow strength, Eq. (6), is given by G() = 1 +bGexp(c[0273]) [exp(c[0])1] (8) withbG= 0:1406 andc= 0:00793/K. In Eq. (8),  and  0are in K and  0= 293K. Also, the initial temperature is taken to be uniform and 293K. The initial void volume fraction is taken to be zero and the evolution of the void volume fraction is governed by _f= (1f)dp:I+_fnucl (9) where the rst term on the right hand side of Eq. (9) accounts for void growth and the second term accounts for void nucleation. The constitutive framework employed in this work includes many of the hardening and softening mechanisms shown by ductile metallic materials. The values of the constitutive 5 and material parameters used here follows the work of Srivastava et al. (2014). Most ma- terial parameters, such as elastic constants and reference yield stress, are representative of aluminum alloys. The initial density, however, is taken to be greater than that for aluminum to increase the stable time increment in the dynamic calculations. 2.2. Numerical implementation A mode I small scale yielding boundary value problem is analyzed for a slice of material with an initial through thickness crack, as shown in Fig. 1. The slice of material analyzed has the dimensions hxhyhz, where,hx=hy= 0:4m andhz=hx= 0:0125. The nite element mesh consists of 428 ;256 twenty node brick elements giving 1 ;868;230 nodes and 5;604;690 degrees of freedom. Ten uniformly spaced elements are used through the thicknesshz. A uniform 20864 in-plane ( xyplane) mesh is used in a 0 :02m0:006m region immediately in front of the initial crack front with in-plane elements of dimension 9:62105m by 9:38105m. The element dimension ex= 9:62105m serves as a normalization length. The initial crack front with an initial opening b0= 1:95exlies along (0;0;z) and the crack faces remain traction free. In the xyplane, velocity boundary conditions corresponding to linear isotropic mode I crack tip eld are prescribed, _ux(R;;z;t ) =2 (1 +)_KI Er R 2cos 2 12+ sin2 2 (10) _uy(R;;z;t ) =2 (1 +)_KI Er R 2sin 2 22cos2 2 (11) Here,R2=x2+y2and= tan1(y=x) for points on the boundary of the region analyzed, _KIis the prescribed rate of increase of the mode I stress intensity factor. Also, symmetry boundary conditions are imposed on z= 0,hzplanes. The value of _KI= _00pex = 1359:4 is prescribed for it to be same as in Srivastava et al. (2014, 2017). The nite element calculations are based on the dynamic principle of virtual work using a nite deformation Lagrangian convected coordinate formulation. The calculations are based on the dynamic principle of virtual work for numerical convenience, hence to minimize the wave e ects, initial velocity elds at t= 0 consistent with Eqs. (10) and (11) are prescribed throughout the region analyzed. Eight point Gaussian integration is used in each twenty-node element for integrating the internal force contributions and twenty-seven point Gaussian integration is used for the element mass matrix. Lumped masses are used so that the mass matrix is diagonal. The discretized equations are integrated using the explicit Newmark -method ( = 0) (Belytschko et al., 1976). The constitutive updating is based on the rate tangent modulus method proposed in Peirce et al. (1984), while material failure is implemented via the element vanish technique proposed in Tvergaard (1982a). When the value of the void volume fraction fat an integration point reaches 0 :9ff, the value of fis kept xed so that the material deforms with a very low ow strength. The entire element is taken to vanish when three of the eight integration points in the element have reached this stage. 6 Figure 1: Sketch of the initially cracked block of material analyzed together with the magni ed view of the three-dimensional random distributions of inclusions in front the initial crack. For the two distributions shown: (top right) radius of the inclusions, r0= 1:5exand volume fraction, n= 0:024; and (bottom right) radius of the inclusions, r0= 3:0exand volume fraction, n= 0:024. The parameter exis the normalization length parameter. 3. Microstructure generation In the calculations the material microstructure is characterized by two populations of void nucleating particles: (i) uniformly distributed small particles that are modeled by plastic strain controlled void nucleation; and (ii) large particles or inclusions that are modeled as stress controlled nucleation sites. In each case, void nucleation is assumed to be described by a normal distribution following Chu and Needleman (1980). For plastic strain nucleation _fnucl=k_ ; k=fN sNp 2exp" 1 2N sN2# (12) withfN= 0:04,N= 0:3 andsN= 0:1. Parametric studies are also carried out to explore the e ect of the mean equivalent plastic strain to void nucleation in the matrix material, N, on the fracture toughness for xed microstructures. For stress controlled nucleation _fnucl=k[_+ _h]; k=fN sNp 2exp" 1 2+hN sN2# (13) if (+h)(+h)max, where the maximum is taken over the previous mechanical history, and@(+h)=@t> 0. Otherwise, k= 0. To generate the random distribution of inclusions (modeled discretely as stress controlled nucleation sites, Eq. (13)) within the uniform mesh region of volume, Vu= 200ex60ex50ex, in front of the initial crack front, we rst randomly generate predetermined number of inclusion, Nincl, centers (x0;y0;z0) with the restriction that the center to center distance of two neighboring inclusions is at least twice 7 the radius, r0. The value of fNin Eq. (13) at a material (Gaussian) point ( x;y;z ) in the initial undeformed con guration, is then assigned as, fN=8 < :fNforp (xx0)2+ (yy0)2+ (zz0)2r0 0 forp (xx0)2+ (yy0)2+ (zz0)2>r 0(14) The values fN= 0:04,sN=0= 0:2 andN=0= 1:5 are used in the calculations. The inclusion volume fraction n, and the mean inclusion spacing l0are estimated as, n= (Nincl 4r3 0=3)=Vuandl0= (Vu=Nincl)1=3, respectively. The inclusions introduce microstructure- based characteristic length-scales into the formulation, one associated with their spacing and the other with their size. Two random distributions of inclusions with n= 0:024, and r0= 1:5exand 3:0exin the volume Vu, ahead of the initial crack front are shown in Fig. 1. 4. Results and discussion The calculations are carried out for three radii, r0= 1:5ex, 3:0exand 4:5ex(whereexis the normalization length parameter) of inclusions, and for a xed r0the calculations are carried out for six inclusion volume fractions, n= 0:012, 0:024, 0:048, 0:071, 0:095 and 0:19. In addition, for a xed r0andn, the calculations are carried out for ve random distributions of inclusion centers i.e. ve realizations of the same overall microstructure. Parametric studies are also carried out to explore the e ect of the matrix material properties on the micromechanisms of ductile crack advance and fracture toughness of the material. 4.1. Crack growth resistance Crack growth resistance ( JR) curves, plots of Jversus a, are extracted from all the calculations following the procedure adopted in Srivastava et al. (2014); Osovski et al. (2015a,b); Srivastava et al. (2017). The value of Jis calculated from the applied value of KIusing the relation (Rice, 1968), J=K2 I(12) E(15) and the crack growth extension,  a, is obtained by averaging the maximum projected length of thef= 0:1 contours through the thickness of the block of the material analyzed, Fig. 1. The normalized JRcurves for two inclusion radii, r0= 1:5exand 4:5exand for each r0, two inclusion volume fractions, n= 0:024 and 0:095 are shown in Fig. 2. For a xed r0andn, the crack growth resistance curves are shown for three realizations of the same overall microstructure. In Fig. 2, the value of Jis normalized by 0l0, while ais normalized byl0. As seen in Fig. 2, for a xed inclusion volume fraction, the level of the normalized JRcurve is greater for the smaller inclusion radius, and for a xed inclusion radius, the level of the normalized JRcurve is greater for the smaller inclusion volume fraction. The variation in the JRcurves between di erent realizations of the microstructure with a xed r0andn, mainly stems from the initial stage of crack growth. The initiation of crack growth 8 (a) (b) Figure 2: Applied normalized J,J=(0l0), versus normalized crack extension,  a=l0, curves for two inclusion radii,r0= 1:5exand 4:5exwith overall inclusion volume fractions, (a) n= 0:024 and (b) n= 0:095. For a xedr0andn, results are shown for three random distributions of the inclusions. depends on the interaction of the initial crack front with the nearest inclusions. Thus, for smaller inclusion volume fraction, Fig. 2a, the variation among the various realizations is greater than the larger inclusion volume fraction, Fig. 2b. The fracture initiation toughness, JIC, is computed based on a procedure mimicking the ASTM:E1820-11 (2011), where a power law of the form J=A(a)Bis t to the portion of theJRcurve in between two exclusion lines, J= 20(a1:5ex) andJ= 20(a5ex). The value of JICis then de ned as the intersection of the curve J=A(a)Band the lineJ= 20(a2ex). The variation of normalized fracture toughness, JIC=(0l0), with inclusion volume fraction, n, for three inclusion radii are shown in Fig. 3. The error bars in Fig. 3 show the standard errors (standard deviation divided by the square root of the number of realizations) for realizations of inclusion distributions having the same value of r0andn. The dash-dot lines in Fig. 3 for r0= 1:5exshow a t to JIC=(0l0) =C1nC2+C3r0+C4, and forr0= 3:0exand 4:5exshow a t to JIC=(0l0) =C3r0+C4, whereC1= 0:01,C2=0:8, C3=0:1exandC4= 0:85. The results presented in Fig. 3 show that for a xed set of material parameters, for the smallest inclusion radius considered here, r0= 1:5ex, the value ofJIC=(0l0) decreases rapidly with increasing nforn0:071 and then tends to saturate. But for the two larger inclusion radii, r0= 3:0exand 4:5ex, the value of JIC=(0l0) is rather independent of n, for the range of nconsidered. In all the previous works, on the interaction of: (i) explicitly modeled two-dimensional distribution of (cylindrical plane-strain) voids or three-dimensional distribution of (spher- ical) voids with a blunting crack tip (Rice and Johnson, 1970; McMeeking, 1977; Aravas and McMeeking, 1985a,b; Hom and McMeeking, 1989; Kim et al., 2003), or (ii) explicitly 9 Figure 3: Variation of normalized fracture toughness, JIC,JIC=(0l0), with inclusion volume fraction, n, for three inclusion radii, r0= 1:5ex, 3:0exand 4:5ex. The lines are the tted curves and the data labeled, Exp, are the experimental results summarized by McMeeking (1977). In the experiments, the crack-tip opening displacement at initiation of crack growth, c, were measured that are converted to JIC, following the relation, JIC=m0c(Shih, 1981), where m= 1. modeled two-dimensional or three-dimensional distribution of (void nucleating) inclusions with a growing crack tip (Needleman and Tvergaard, 1987; Srivastava et al., 2014), the focus was con ned on correlating macroscopic fracture initiation toughness de ned using c=l0orJIC=(0l0) with the microstructural parameter de ned either using the ratio of mean void/inclusion spacing and size, l0=2r0or void/inclusion volume fraction, n(where n/(r0=l0)3). These calculations predict that for a xed set of material parameters, the value ofc=l0orJIC=(0l0) is uniquely related to l0=2r0orn, and the value of c=l0or JIC=(0l0) initially increases rapidly with increasing l0=2r0(decreasing n) and then tends to saturate. These calculations did not predict the variation in the values of JIC=(0l0) for a xedn. Our results, Fig. 3, show that for high inclusion volume fraction JIC=(0l0) depends on the inclusion size whereas for low inclusion volume fraction JIC=(0l0) depends on both the inclusion size and the inclusion volume fraction. The normalized experimental results plotted in Fig. 3 follows the same trend even though the experimental results are for a wide variety of engineering metals and alloys. The experimental results shown in Fig. 3 were summarized by McMeeking (1977), further details of the experimental results can be found in the