Three-dimensional ductile fracture analysis with a hybrid multiresolution approach and microtomography

Journal of the Mechanics and Physics of Solids 61 (2013) 2108–2124

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

Three-dimensional ductile fracture analysis with a hybrid multiresolution approach and microtomography Shan Tang a,b, Adrian M. Kopacz b, Stephanie Chan O'Keeffe c, Gregory B. Olson c, Wing Kam Liu b,d,e,n a

College of Material Science and Engineering, Chongqing University, Chongqing 400044, China Department of Mechanical Engineering, Northwestern University, Evanston, IL 60208, USA Department of Material Science and Engineering, Northwestern University, Evanston, IL 60208, USA d School of Mechanical Engineering, World Class University (WCU) Program, Sungkyunkwan University, Korea e Adjunct Professor under the Distinguished Scientists Program Committee at King Abdulaziz University (KAU), Jeddah, Saudi Arabia b c

a r t i c l e in f o

abstract

Article history: Received 9 November 2012 Received in revised form 20 May 2013 Accepted 10 July 2013 Available online 20 July 2013

A modified-JIC test on CT (compact tension) specimens of an alloy (Ti-Modified 4330 steel) was carried out. The microstructure (primary and secondary inclusions) in the fracture process zone and fracture surface are reconstructed with a microtomography technique. The zig-zag fracture profile resulting from nucleation of microvoid sheets at the secondary population of inclusions is observed. Embedding the experimentally reconstructed microstructure into the fracture process zone, the ductile fracture process occurring at different length scales within the microstructure is modeled by a hybrid multiresolution approach. In combination with the large scale simulation, detailed studies and statistical analysis show that shearing of microvoids (the secondary population of voids) determines the mixed mode zig-zag fracture profile. The deformation in the macro and micro zones along with the interaction between them affects the fracture process. The observed zigzag fracture profile in the experiment is also reasonably captured. Simulations can provide a more detailed understanding of the mechanics of the fracture process than experiments which is beneficial in microstructure design to improve performance of alloys. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Multiresolution Zig-zag fracture path Void growth Ductile fracture

1. Introduction Ductile fracture of engineering alloys is a multi-stage process involving nucleation, growth and coalescence of voids across different length scales of microstructure (micron size primary particles and submicron size secondary particles) during plastic deformation. Generally, primary inclusions debond from the metal matrix and form primary voids under minimal plastic deformation. When the localized deformation between primary voids is intensified, void linkage can be triggered by nucleation of a secondary population of particles on the intervoid ligament. In these multi-stage and multiscale processes, shear coalescence of secondary voids (microvoids) plays an important role in the fracture process. However, it remains elusive. McVeigh et al. (2007) carried out simulations on a representative volume element (RVE) with randomly distributed voids under shear loading. Their work clearly demonstrated that microvoid sheeting and strong interactions between microvoids can result in the formation of shear bands and consequently failure between primary voids, as shown

n Corresponding author at: Department of Mechanical Engineering, Northwestern University, 2145 N. Sheridan Rd. B-224. Tel.: +1 847 491 7094; fax: +1 847 491 3915. E-mail address: [email protected] (W.K. Liu).

0022-5096/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jmps.2013.07.007

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2109

experimentally by Cowie et al. (1989). Benzerga et al. (2004a) and Benzerga and Leblond (2010) also characterized the crack growth in round notched bars with 2D SEM experimentally. They showed that a shear-like fracture is favored in some specific sections which leads to a zig-zag fracture path. In our experiments, the 3D zig-zag fracture surface is also reconstructed (shown in Sections 2 and 4). The shear band caused by shear softening of microvoids is clearly observed on the zig-zag fracture path. To predict the ductile fracture process in alloys, two different computational approaches were employed in the past (refer to Fig. 1). The first approach employs a constitutive model, such as the GTN (Gurson–Tvergaard–Needleman) model, to account for the void nucleation and growth process by smearing out the real microstructure. The constitutive model is then implemented into a finite element code to simulate the initiation and crack growth (Xia and Shih, 1995). In this approach, a length scale is introduced which is related to the spacing between primary voids, usually element size. Pardeon and Hutchinson (2003) improved the GTN model by incorporating the evolution of void shape and coalescence in studying the relationship between fracture toughness and microstructure successfully. However, this was done at the cost of higher complexity and additional calibration parameters. These parameters are generally difficult to extract from experiments. Benzerga et al. (2004b) and Benzerga and Leblond (2010) also proposed a complex model to account for initial anisotropy and microstructure evolution (plastic anisotropy, porosity, void shape, orientation and spacing). It is also supplemented by a micromechanical model of void-coalescence. Anisotropic crack growth in round notched bars was successfully simulated with the proposed model. Morgeneyer and Besson (2011) studied the flat to slant fracture transition with a Gurson-type model and a shear-based void nucleation term. In these works, simulations are carried out with a Gurson-type model which usually is derived through a homogenization procedure resulting in smearing of the voids. They do not distinguish the contribution of primary and secondary voids on the fracture process. Because the crude mesh used in the simulation cannot capture the zig-zag fracture path accurately, a full understanding of the mechanism for zig-zag fracture was not developed and quantification of the zig-zag fracture path was not carried out despite showing that accounting for shear-dominated fracture processes may lead to a better comprehension and prediction of slant fracture (Morgeneyer and Besson, 2011). Another important issue is that the numerical results depend on the finite element discretization because the localization caused by shear softening which is mesh dependent. To remove mesh sensitivity, a characteristic length scale is necessarily introduced for regularizing the governing equations. Then, a phenomenological nonlocal model is employed to model alloy plates of predrilled holes under uniaxial load states (Weck et al., 2008). Compared to their experiments, Weck et al. (2008) showed that the nonlocal model is not only necessary for controlling mesh sensitivity of the numerical solutions but it can also improve the agreement between simulation and experiments. In the second approach (refer to Fig. 1), finite element computations are carried out with primary voids explicitly modeled by a refined FE mesh (Tvergaard and Hutchinson, 2002; Kim et al., 2003; Petti and Dodds, 2005). A length scale related to the growth and coalescence process is introduced by considering the spacing of real voids (not the element size in the first approach). In this approach, primary voids are ideally confined to a flat crack plane. Therefore, the zig-zag behavior cannot be captured. A criterion for failure of the intervoid ligaments, thus simulating crack growth, is required to describe the shear localization within the ligament due to a secondary population of voids (microvoids are not included). It is obvious that microvoid shearing is not taken into account. The conventional plasticity model is employed and potential void size effects at the submicron scale is excluded. Although this approach, based on discrete voids, reveals basic aspects of the mechanisms of crack growth, it is generally accepted that computational complexity renders it impractical for engineering crack growth studies with realistic three-dimensional void populations or even of two-dimensional void populations where the voids do not lie on the fracture plane. Obviously, the important role of the secondary population of particles is excluded in both approaches. With the advancement of computational capability, it is possible to model the primary voids involved in the fracture process zone directly. Therefore, the role of microvoids (secondary voids) on ductile fracture can be clarified and the related fracture mechanisms can be revealed. Further, we can make a direct comparison with 3D experimental reconstruction qualitatively and quantitatively. Along this line, Tian et al. (2010) made progress on predicting the zig-zag path seen in ductile crack growth with a phenomenological material model accounting for the effects of microvoids. Although their work gave some insight into the ductile fracture process within the multiresolution framework and showed a zig-zag fracture profile similar to experiments, their work did not clarify several important issues: 1. The underlying mechanisms that cause the zig-zag fracture path are not clear because the employed material model is phenomenological.

Fig. 1. The conventional and present hybrid multiresolution approaches for ductile fracture.

2110

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2. The deformation of the macro and micro zones in the multiresolution framework, the interaction between them, and their effect on the ductile fracture process is not clearly demonstrated. 3. The fracture indicators such as crack opening distance, void growth ratio, and zig-zag wave length are not accurately captured when compared to experiments. 4. The number of primary voids involved is much less than that known to exist in the specimens used for experimentation. In the present work, a hybrid multiresolution model is employed to simulate crack growth of ductile fracture in which the advantages of each of the two approaches adopted in the past are incorporated. The proposed hybrid approach is based on advancement of experimental techniques which allow primary and secondary particles within the fracture process zone to be fully reconstructed (Chan et al., 2013). Obviously, we cannot model all of the primary and secondary particles explicitly in the fracture process zone. Based on the second modeling approach with discrete voids, we embed the primary particles into our fracture process zone directly. Generalizing the first approach with a constitutive model to consider void damage for crack growth, we apply a homogenized multiresolution model in the material at the length scale of the secondary particles to capture shear localization. In the homogenized model, we consider the third invariant of stress to have an impact on the evolution of damage which represents the effects of microvoid shearing. A length scale is also introduced which can be gleaned from cell analysis of the secondary population of voids which can then be used to regularize the softening behavior through the constitutive laws (Liu and McVeigh, 2008; McVeigh and Liu, 2008; Tang et al., 2012a; Tang et al., 2012b). Furthermore, the material models involved are physically motivated and it is easy for us to identify the related mechanisms. We will show that a combination of direct modeling of primary voids and multiresolution homogenization can simulate the mixed mode zig-zag fracture path and is comparable to our experimental observation. Acknowledging the incredible efforts of LSTC (LS-DYNA user's manual, LS-DYNA theory manual), we implemented the proposed multiresolution model into LS-DYNA (Tang et al., in press), which renders the possibility to analyze realistic three-dimensional primary void populations using parallel computations. 2. Three-dimensional multiresolution theory A detailed treatment of the multiresolution theory for alloys is given in Vernerey et al. (2007, 2008) and McVeigh and Liu (2008, 2009, 2010). We extend that work here for 3D finite element implementation using LS-DYNA (Tang et al., in press; Belytschko et al., 2000). 2.1. The weak form of multiresolution theory The variation of internal virtual power in the multiresolution theory is Z N I δP int ¼ r : δD þ ∑ fβI : ðδDI δDÞ þ β : δDI ▿gdΩ Ω

ð1Þ

I¼1

I

where r is the macroscopic stress, βI and β are micro-stress and micro-stress couples at the Ith scale. D is the rate of macro deformation and DI is the rate of deformation at the Ith scale. I ranges from 1 to N and N is the total number of length scales (in the present work, N ¼1). The variation of kinematic virtual power is Z N _ I  II Þ : δDI dΩ ð2Þ δP kin ¼ ρv_  δv þ ∑ ðD Ω

I¼1

where v is macro velocity and II is inertial tensor at the Ith scale. It is defined Z 1 ρI x⊗x dΩI II ¼ ΩI ΩI

ð3Þ

and N

ρ ¼ ρ þ ∑ ρI :

ð4Þ

I¼1

where ρ is the macroscopic density, ρI the density at the Ith scale and ΩI is the volume at the Ith scale. In the kinematic virtual power, the higher order terms have been abandoned. To further simplify the description, we assume that the microdomains do not rotate relative to the macroscopic domain. Therefore, only the effects of stretching and shearing are considered at the micro-scale and the velocity gradient at the Ith scale LI can be replaced by its symmetric counterpart DI (Vernerey et al., 2008). The variational virtual external power takes the following form: Z N t  δv þ ∑ TI : δDI dΓ ð5Þ δP ext ¼ Γt

I¼1

where t is the traction force vector and TI is the traction couple at the Ith scale.

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2111

2.2. Formulation of plasticity models In the present work, we model the primary particles directly as discussed in the introduction. That is, we embed the reconstruction of the experimentally detected primary particles into our simulation domain directly. To simplify computation, we do not consider the nucleation process of particles. That is, both primary particles and secondary particles in the simulation are modeled as voids (Particles are debonded from the matrix already). Primary particles (voids) are modeled directly with known void center location, size and shape (the same as experiments), while the secondary particles (voids) are modeled by a homogenized multiresolution model proposed in Section 2.2. We only consider the initial void volume fraction of secondary particles in this homogenized model which is measured in our experiment. Following the procedure proposed by Tang et al. (2012b), the deformation at the length scale of the secondary particles (voids) is decomposed into the macro zone and the micro zone (see Fig. 2). The macro zone can be approximated by a cell with a single void. Obviously, the macro zone cannot consider the interaction and inhomogeneous deformation between secondary population of voids (microvoids). Therefore, we design a micro zone with a length scale to mimic the highly inhomogeneous deformation and strong interactions between microstructures. The plasticity models for the macro zone and the micro zone are formulated in the following subsections. 2.2.1. Material model of macro zone The macro zone is used to describe the response resulting from a homogenous deformation state, i.e., excluding the interaction of microvoids and strain gradient effects. We apply the GTN model and derive the set of governing non-linear equations. The numerical implementation of the GTN model has been widely addressed in the literature (Aravas, 1987). In the following, we derive the equations following Aravas (1987). We modify the derivation to be consistent with the current rate formulation. The rate of deformation tensor D at the microscopic scale is decomposed into additive elastic and plastic parts Dij ¼ Deij þ Dpij

ð6Þ

where superscripts e and p represent the elastic and plastic parts, respectively. The objective rate of Cauchy stress s p s∇ ij ¼ C ijkl ðDkl Dkl Þ

where C ijkl ¼ 2G

1  2

  δik δjl þ δjk δil 13δij δkl þ Kδij δkl :

The macro shear modulus and bulk modulus are G and K, respectively. The yield function, or the flow potential, of the macro zone is given as     s 2 3q2 sm eq 2  1 þ q21 f ¼ 0; þ 2q1 f cosh Φ¼ s 2s

ð7Þ

ð8Þ

ð9Þ

where seq , sm , s and f are the effective stress, hydrostatic stress, flow stress of the matrix and current void volume fraction, respectively. The parameters q1 and q2 were introduced to improve the predictions of the original Gurson model. The flow rule gives  ∂Φ 1 ∂Φ ∂Φ ¼ λ_ δij þ nij : ð10Þ Dpij ¼ λ_ ∂sij 3 ∂sm ∂seq where λ_ is the plastic multiplier and nij ¼

3Sij : 2seq

Fig. 2. Multiresolution decomposition of deformation at the length scale of secondary particles.

2112

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

Defining ∂Φ ; Dpm ¼ λ_ ∂sm

ð11Þ

∂Φ ∂seq

ð12Þ

Dpeq ¼ λ_

Then we can write Dpij ¼ 13Dpkk δij þ Dpeq nij : Substituting this back into Eq. (7), we obtain s▿ij ¼ 2GDij ′ þ KDkk δij KDpkk δij 2GDpeq nij

ð13Þ

where Dij ′ ¼ Dij 13Dkk : The void volume fraction f is due to the expansion of the existing void f_ ¼ ð1f ÞDpkk :

ð14Þ

The plastic work in the matrix is s ε_ ¼ ð1f Þsij Dpij ; p

ð15Þ

where ε_ is the effective plastic strain rate in the matrix. We assume that the matrix response is governed by the power law relation p

s ¼ sy ð1 þ ε p =ε0 Þn

ð16Þ

where sy is the initial yielding stress of the matrix, ε0 the initial yielding strain and n is the hardening exponent. In the classical GTN model, Eq. (14) only considers the contribution of spherical void growth. However, during the fracture process, in addition to the spherical void growth, shearing of microvoids can result in void distortion, which can contribute to softening during the coalescence process of microvoids. Suppose that the material at the length scale of the microvoids is subjected to pure shear loading. Microvoid shearing can still cause the formation and propagation of shear bands as shown experimentally (Cowie et al., 1989, Benzerga et al., 2004a) and via simulation (McVeigh et al., 2007) without spherical void growth. Motivated by this point, Nahshon and Hutchinson (2008) as well as Xue (2008) modified Eq. (14) to include a contribution of shearing of voids to the damage growth rate in which states of pure shear stress, in a manner, leaves the relation unaltered for axisymmetric stress states (they apply it generally at the length scale of the primary voids). In their modification, f is no longer directly tied to the plastic volume change as in the classical GTN model. Instead, it must be regarded either as an effective void volume fraction or simply as the damage parameter, the latter of which is often the case when the GTN Model is applied to materials with non-spherical voids. However, the idea is still physically motivated and has advantages over the phenomenological material model employed in Tian et al. (2010). It can benefit us in identifying the corresponding fracture mechanisms which will be shown in our results section. Extending their idea to microvoids, an extra term which can consider the microvoid shearing should be imposed on the right side of Eq. (14): f_ ¼ ð1f ÞDpkk þ kω f ωðsÞ

Sij Dpij se

where ωðsÞ is defined   27J 3 2 ωðsÞ ¼ 1 2s3e in which the third invariant of stress J 3 ¼ 13 Sij Sik Sjk . ωðsÞ lies between 0 and 1. When ω ¼ 0, it is just the axisymmetric stress state. When ω ¼ 1, it is a stress state composed of pure shear stress plus a hydrostatic contribution. A numerical constant kω is introduced to describe the contribution of microvoid shearing to the material response. The model with the additional microvoid shearing term is called the modified GTN model throughout the paper. If we set kω ¼ 0, the modified GTN model reduces to the classical GTN model. 2.2.2. Material model of micro zone Even when a homogeneous deformation is applied to a unit cell containing microvoids as shown in Fig. 2, highly inhomogeneous deformations and localizations can arise due to the interaction of microvoids. The classical homogenization approach only takes into account the average stress and strain. However, it does not consider the fluctuation of deformation over the average deformation in the unit cell and nonuniform deformation on the boundary of the unit cell. In Eq. (1), we can see that β and β are conjugated to the difference of deformation of the microzone and the macrozone and their corresponding gradients. So, just the first-order approximation of the fluctuation of deformation over the average

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2113

deformation is considered. We then modify Fleck and Hutchinson (1993) such that it is applicable within our multiresolution framework to obtain β and β in this section. In the proposed model for the microzone, a length scale of representing the microzone is introduced. In our previous works (Vernerey et al., 2007, 2008; McVeigh and Liu, 2008, 2009, 2010), we have shown that the introduction of this length scale can regularize the macroscale softening within the constitutive relations. It also can describe the inhomogeneous deformation and strain gradient effects caused by the microstructure. In this section, stresses and strains are defined at the Ith micro zone unless otherwise stated (the index is ignored). We assume that the elastic moduli CI in the micro zone is isotropic, that is,     ð17Þ C Iijkl ¼ 2GI 12 δik δjl þ δjk δil 13δij δkl þ K I δij δkl : where GI and KI are the shear and bulk moduli at the Ith micro zone, respectively. For the micro stress at the Ith scale, the stress–strain rate relation can be written I ~ I ~ β∇ ij ¼ 2G D ij ′e þ K D kk δij : e

ð18Þ

where the relative deformation between the micro zone and the macro zone is ~ ij ¼ DI Dij ; D ij

ð19Þ

~ ij ′ ¼ D ~ ij  1 D ~ the deviatoric strain rate D 3 kk . The micro stress βij can be decomposed into deviatoric and spherical parts m βij ¼ βdev ij þ β δij

ð20Þ

where

  1 βdev ij ¼ βij 3 βll δij ;

βm ¼ 13βll :

ð21Þ

Using micro stress decomposition, the stress–strain rate relation can be simplified as follows: βdev∇ ¼ 2GI D~ ij ′e; ij

I ~ βm∇ ¼ 13β∇ ll ¼ K D kk e

ð22Þ

where ~ ij 1D ~ ~ ij ′ ¼ D D 3 kk : The micro stress couple in the micro volume is (Vernerey et al., 2007, 2008) " !e # Z ∂LImn 1 ∇  I β ijk ¼ C Iikmn yl yj dΩ: ∂xl Ω ΩI

ð23Þ

ð24Þ

Assuming that the shape of ΩI is a cubic cell with characteristic length scale lI, we can obtain ∇

β ijk ¼ and

 I2   e  e l 2 2GI DIijk þ K I  GI DIllk δij 3 12

  DIijk ¼ 12 LIij;k þ LIji;k :

ð25Þ

ð26Þ

The micro stress couple can also be decomposed dev

m

β ijk ¼ β ijk þ β k δij : Under this decomposition, the higher order stress–strain rate can be expressed as   I I ðl Þ2 1 ∇ ðl Þ2 2 dev∇ m∇ 2GI ðDIijk Þe ; β k ¼ β llk ¼ K I  GI ðDllk Þe : β ijk ¼ 3 3 12 12

ð27Þ

ð28Þ

The plastic potential is similar to the strain gradient model proposed by Fleck and Hutchinson (1993), p Φ ¼ sH e sY ðE Þ ¼ 0

where the generalized effective stress is defined sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 3 dev dev a dev dev H βij βij þ I β ijk β ijk se ¼ 2 l

ð29Þ

ð30Þ

Rt and the micro flow stress is sY ¼ sY0 þ EI Ep which is a function of the generalized effective strain Ep ¼ 0 E_p dt qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ffiffiffiffiffiffi ~ p D~ p þ ðlI2 =a2 ÞDI;p DI;p and a ¼ 18, see Vernerey et al. (2007, 2008) and EI is the Young's modulus at the I th (E_p ¼ 23D ij ij ij ij microscale.

2114

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

From the plastic flow rule, we obtain 3 dev β ~ ij ′p ¼ E_p 2 ij ; D sH e  2 a Dijk ′p ¼ E_p where sH e

l

I

ð31Þ

dev

β ijk

sH e

:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 3 dev dev a dev dev βij βij þ I β ijk β ijk ¼ 2 l

ð32Þ

ð33Þ

The numerical implementation of micro constitutive laws is discussed briefly here. The reader who feels interested in the details can refer to Tang et al. (in press). Suppose we know the information at time step n. From Eqs. (18) and (25) and the choice of the return directions, the scalar parameter Ep fully defines the updated stress state. The numerical solution must be p used to determine values for this scalar parameter such that Dij satisfies the flow rule in Eqs. (31) and (32) and the updated stresses satisfy the yield criterion in Eq. (29) at time step n þ 1=2. The full multiresolution theory, including constitutive models, has been implemented into LS-DYNA. We need to set internal force, external force, and the mass matrix for a single element following the weak forms shown in Section 2.1. 3. Finite element models Before we proceed to the finite element model to simulate the crack growth, the experiments are briefly reviewed. Details are given in Chan (2011). The material involved in the fracture experiments is Ti-modified 4330 steel. Ti-modified 4330 steel, hereafter referred to as “Mod4330”, is a high strength low martensitic steel. A modified-JIC fracture toughness test similar to the standard test described in ASTM E1820 was carried out. Miniature single-edge notched CT specimens with a narrow starter for a fatigue crack were machined from halves of broken Mod4330 Charpy bars (see Fig. 3a). The microstructure reconstructed in the fracture process (the reconstructed microstructure is in the deformed shape and size)

Fig. 3. (a) Fracture experiment, reconstructed microstructures in the fracture process zone and fracture surface; (b) the corresponding computational modeling strategy.

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2115

and the fracture profile are shown in Fig. 3a. It can be seen clearly that the fracture profile is a zig-zag connected by the shear band resulting from microvoid shear coalescence. Based on the same technique, the microstructure outside the fracture process zone in the same specimen is also reconstructed (the microstructure obviously is not the same as those involved in the real fracture process but it is closer to the initial state of the microstructure from a statistical point of view). The position, radius and shape of primary and secondary particles are included in the digital data files after the experimental reconstruction. As discussed in the introduction, primary voids nucleate at low strain levels due to decohesion of the primary particles. To distinguish primary particles from secondary particles, a cut-off radius can be obtained through statistical analysis of the experimental data (particles whose sizes are less than the cut-off radius are considered to be secondary particles). Adopting small scale yielding assumption, the experimental specimen is modeled using a circular region whose radius is 105 times larger than the size of the representative element in the fracture process zone (see Fig. 3b). A Python script was written to decompose the simulation domain into small parts. The mesh is then generated in different parts and glued together. The mesh in the outer elastic zone is the coarsest while the finest portion of the mesh is in the fracture process zone. Initially the structure is meshed with elements much smaller than the primary void size. Once this mesh is created, the elements that are in the location of the primary voids are removed based on the experimentally determined void centers and radii. We use a very refined mesh on the surfaces of the primary voids. The size of the primary void ð 4 2 μmÞ is much larger than the mesh size ð≈0:2 μmÞ. This approximation seems reasonable. In order to rationalize fracture behavior, a coordinate system is set up such that the x-axis represents the crack propagation direction, the y-axis represents the crack opening direction and the z-axis represents the thickness direction. At the outer boundary of the simulation domain, a displacement increment corresponding to mode I plane strain K-field is imposed (Kanninen and Popelar, 1985) rffiffiffiffiffiffi  h i KI R θ ux ¼ κ1 þ 2 sin 2 θ cos 2 2μ 2π rffiffiffiffiffiffi   KI R θ  uy ¼ κ þ 12 cos 2 θ ð34Þ sin 2 2μ 2π where KI represents the mode I stress intensity factor, the position of a point on the outer boundary is located by its distance R to the crack tip and the angle θ. The shear modulus is μ and κ ¼ 34ν for plane strain conditions where ν is Poisson's ratio. Boundary conditions prescribed on the plane normal to z-direction are set up to zero. All of the extra degrees of freedom are set to zero on the outer boundary (that is, inhomogeneous deformations disappear on the outer boundary). The reference stress intensity factor is K ref ¼ 102 MPa m1=2 . The material parameters used in the macro and micro zones are given in Tables 1 and 2, respectively, which are obtained from experimental data of Chan (2011) and Vernerey et al. (2007, 2008). 3.1. Fracture process (17 primary particles) To give insight into the fracture process, we embedded 17 particles (with exact shape and position give through destructive reconstruction methods) into the fracture process zone in our simulation. The finite element mesh is shown in Fig. 4. The simulation domain size, notch size and degrees of freedom are addressed in Table 3 where we summarize different simulation cases. We first consider three different levels of the initial void volume fraction: f 0 ¼ 0:06%; 0:1%; 0:2% without microvoid shearing (classical GTN model). They are refereed to as C1-1, C1-2 and C1-3, respectively, in Table 3. Subsequently, we investigate the case with an initial void volume fraction: f 0 ¼ 0:1% with the microvoid shearing (modified GTN model) (it is referred as C1-4 in Table 3). We will show that the modified GTN model can lead to a zig-zag fracture profile, larger void growth ratio in the transverse directions and relatively lower initial fracture toughness which are more consistent with experimental observations. We will also demonstrate the concurrent interaction of the macro and micro effective strains due to submicron and micron-sized void growth, interaction, and coalescence to form micro and macro minicracks which can lead to the formation of a fracture process zone. Table 1 Material properties for the macro zone. K (GPa)

G (GPa)

sy (GPa)

q1

q2

n

176.7

81.54

1.72

1.29

0.982

0.1

Table 2 Material properties for the micro zone. K 1 (GPa)

G1 (GPa)

sY0 (GPa)

l μm

17.67

8.154

0.045

1

1

2116

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

Fig. 4. Finite element meshes for case 1 with 17 primary particles.

Table 3 Summary of parameters and simulation results. The meanings for abbreviations in the table are PS: the average primary particle size; NS: notch size ðμmÞ; DM: domain size in xyz coordinate system ðμmÞ; DOF: degrees of freedom; COD: crack opening distance ðμmÞ when the full cracking takes place in the simulation domain; TG: toughness measured by stress intensity factor (MPa m1/2). It is defined when the full cracking takes place in the simulation domain; VGR: maximum void growth ratio in transverse direction (T) and loading direction (L); CL: crack length at initial full cracking ðμmÞ. In experiments, COD is measured at the unloading configuration. The fracture toughness is not measured accurately (see discussion in Section 2). Case number

PS

NS

DM

DOF

COD

TG

VGR

CL

Experiments C1-1 C1-2 C1-3 C1-4 C1-5 C1-6 C1-7 C2-1 C2-2 C3

4 8 8 8 8 4 4 4 4 4 4

3 25 25 25 25 15 15 15 15 15 3

350  250  500 3000  3000  50 3000  3000  50 3000  3000  50 3000  3000  50 3000  3000  50 3000  3000  50 3000  3000  50 3000  3000  80 3000  3000  80 3000  3000  100

1.4  107 1.4  107 1.4  107 1.4  107 5.6  107 5.6  107 5.6  107 1.9  108 1.9  108 2.2  108

15 79 77 51 48 41 43 41 43 51 23

77 219 219 130 93 117 117 115 129 122 90

T:9.0 L:5.0 T:2.7 L:3.1 T:2.8 L:3.2 T:2.9 L:3.5 T:4.4 L:2.8 T:2.5 L:3.0 T:2.5 L:3.2 T:2.4 L:3.0 T:2.6 L:3.0 T:6.0 L:3.2 T:7.4 L:3.8

350 209 222 214 224 159 171 159 251 223 312

We first set kω ¼ 0, that is, we did not consider microvoid shearing (classical GTN model). Results are shown in Fig. 5. The final fracture surface is visualized by the effective strain in the macro zone. It can be seen that the fracture path is very sensitive to the initial volume fraction of microvoids. When the initial volume fraction of microvoids is low, the fracture will advance somewhat along the x-axis and then orient itself on the plane 451 to the x-axis (Fig. 5a). It is a typical mixed mode I and II fracture. When the initial volume fraction of microvoids is higher ðf 0 ¼ 0:1%Þ, the angle of the oriented fracture path is less than that of the case with f 0 ¼ 0:06% (Fig. 5b). When the initial volume fraction of microvoids is much higher ðf 0 ¼ 0:2%Þ, the fracture will follow a planar path (Fig. 5c) which is indicative of a typical mode I fracture. The results shown in Fig. 5 suggest that larger initial void volume fraction can result in a mode I fracture process. When the initial volume fraction of microvoids is very small, the matrix material can reduce to the classical J2 Plasticity employed in Tvergaard and Hutchinson (2002) and Kim et al. (2003) in their idealized discrete void model. Our results suggest that if the primary voids are not confined to a planar plane, it may cause the mode II fracture process to develop in their simulations instead of the mode I fracture process in their simulations. We also observe that primary voids grow faster in the transverse direction (x-direction). When the primary voids go into the coalescence stage, the void growth in the loading direction (y-direction) will be larger than that in the transverse direction. Finally, the shape of most voids is prolate (the size of void in the loading direction is larger than that in the transverse direction). The same trend in void growth is shown in Tvergaard and Hutchinson (2002) and Kim et al. (2003). The detailed analysis of void growth ratio in the loading and transverse directions will be addressed quantitatively later. It should be noted that when the initial microvoid volume fraction is smaller, full fracture takes place later and exhibits higher initial fracture toughness which is consistent with Xia and Shih (1995). In the idealized discrete model in Tvergaard and Hutchinson (2002) and Kim et al. (2003), the fracture will advance along the plane with discrete primary voids. Compared to the reconstructed fracture surface shown in Fig. 1 (Chan, 2011), we can see that the fracture profiles in Tvergaard and Hutchinson (2002) and Kim et al. (2003) are greatly different from that of our experiments. Let us now consider shearing of microvoids and set kω ¼ 15. Fig. 6 shows serial snapshots of the development of the fracture surface at different levels of far-field loadings. The deformation is visualized by the effective strain of the macro zone. The ductile fracture process starts with the macro plastic deformation accumulation at the blunting crack tip. It is followed by the primary void growth (Fig. 6a). Void growth is mainly along the transverse direction. It should be noted here that void growth starts to take place first at the void ahead of the crack tip. Then shear-driven microvoid coalescence connects the primary voids and macro fracture occurs (Fig. 6b, c). Shear bands resulting from the coalescence of microvoids orient the plane around 451 to the x-axis. Repeating the process, macro fracture propagates forward in a zig-zag fashion

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2117

Fig. 5. Fracture process without microvoid shearing visualized by the effective strain in the macro zone.

(Fig. 6d). In the fracture process, growth of primary voids does not necessarily occur in the sequence of crack propagation; it is in a simultaneous and inter-competitive fashion. To illustrate the interaction between the micro and macro zones and their effects on the ductile fracture process, we further show serial snapshots of the development of the fracture surface at different levels of far-field loadings. They are visualized by the effective strain in the micro zone (shown in Fig. 7) at the same moments as in Fig. 6. When macro plastic deformation accumulates at the blunting crack tip (Fig. 6a), the micro deformation is developed within the region of macro deformation (Fig. 7a). When shear-driven microvoid coalescence connects the primary voids and forms the macro fracture surface, the region with high micro deformation will fail and be removed from the simulation domain (Fig. 7b, c). Finally, micro deformation localizes within the region with high macro deformation, leading to further propagation of the crack. Although the effective strain in the micro zone is smaller than that in the macro zone of the whole fracture process, micro deformation indicates a region with highly inhomogeneous deformation and has implications for possible crack growth directions. Through these parametric studies, we can identify three different fracture patterns (The schematic figures are also drawn in Figs. 5 and 6). Fig. 5b, c is a typical mode I fracture pattern. In this mode, primary voids grow first. The intervoid ligament is separated by the tensile stress and the fracture surface resulting from separation is flat. This fracture pattern is very similar to the idealized discrete void model (fracture is confined to a layer of elements containing voids). In Fig. 5a, the direction of crack growth deviates from the x-axis. It is a typical mixed mode I and II fracture. With the classical GTN model, the zig-zag fracture profile cannot be observed and void growth in the loading direction is larger than that in the transverse direction in the final stages of deformation. In Fig. 6 with the modified GTN model, primary voids grow first and shear coalescence of microvoids connects the primary voids, resulting in the mixed mode zig-zag fracture profile. Void growth in the loading direction is generally less than that in the transverse direction in the final stages of deformation. To examine the evolution of deformation localization, we will show the evolution of plastic strain in the macro zone as a function of the far-field loading at the midpoint between the endpoints (centers of neighboring primary voids) of the line segment (For point 1 in Fig. 8, the end points of the line segment are the crack tip and the center of void v1.) on the fracture path with the classical and modified GTN model. The local strain then represents the behavior of the intervoid region.

2118

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

Fig. 6. Fracture process with microvoid shearing visualized by the effective strain in the macro zone.

Fig. 7. Fracture process with microvoid shearing visualized by the effective strain in the micro zone.

Results are summarized in Fig. 8. Fig. 8a shows the material points (1, 2, 3, 4 and 4′) where the evolution of plastic strain will be plotted. Fig. 8b,c,d plots the evolution of the effective plastic strain vs. the far-field loading for the modified GTN model with initial void volume fraction f 0 ¼ 0:1%, the classical GTN model with initial void volume fraction f 0 ¼ 0:06% and the classical GTN model with initial void volume fraction f 0 ¼ 0:1%, respectively. As shown in Fig. 8, the plastic strain grows much faster when considering microvoid shearing in the modified GTN model. Beyond a specific far-field loading, the plastic strain has a sharp increase which is consistent with experimental measurements in Weck et al. (2008). We also see that the initial fracture takes place earlier with microvoid shearing at around 90% of Kref, which is more consistent with our experiments and the previously reported initial fracture toughness for Mod 4330 steel. The classical GTN model predicts a higher initial fracture toughness. We will now illustrate primary void growth during the fracture process which can be related to fracture toughness through some simple relations. Fig. 9a,b shows the primary void growth ratio in the transverse and loading directions with

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2119

Fig. 8. Evolution of plastic strain vs. far-field loading during fracture process: (a) Midpoints on ligaments for plotting; (b) modified GTN model with microvoid shearing and f 0 ¼ 0:1%; (c) classical GTN model with f 0 ¼ 0:06%; (d) classical GTN model with f 0 ¼ 0:1%.

Fig. 9. Final void growth ratios vs. void position (x-axis): (a) Modified GTN model considering microvoid shearing; (b) classical GTN model.

microvoid shearing and without it (initial void volume fraction f 0 ¼ 0:1%). Voids (v1–v8) are marked in Fig. 8a on the fracture path. In the simulations with the classical GTN model, results showed that transverse void growth dominates in the early stages of deformation due to void–void interactions (see Tvergaard and Hutchinson, 2002 and Kim et al., 2003). This suggests that it is easier for the voids that line up in the transverse direction to coalesce with one another. As such, a void tends to coalesce with another void located along the lateral transverse direction. As the deformation goes on, transverse void growth starts to slow down whereas the growth in the loading direction starts to increase. This reflects the start of void growth mode transition from stable void growth to uniaxial stretching in association with shear localization on intervoid ligaments. Finally, the void growth ratio in the loading direction is larger than that in the transverse direction as shown in Fig. 9b (except the primary void v3). These observations are consistent with Tian et al. (2010). However, for the modified GTN model considering microvoid shearing, the transverse void growth is dominant throughout the fracture process. As shown in Fig. 9a for all the voids examined on the fracture path, the transverse void growth ratio is larger than that in the loading direction. Stereomicroscopy is employed to measure the final void growth

2120

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

ratio on the fracture surface in our experiments. Results showed that void growth ratio in the transverse direction is larger than that in the loading direction. The average void growth ratio in the transverse and loading directions are 4.2 and 3.6, respectively. It is obvious that the modified GTN model with microvoid shearing correctly captures the characteristics of primary void growth. However, the classical GTN model, such as that used in Tian et al. (2010). cannot predict it correctly.

3.2. Fracture process (115 primary particles case) The position and size of primary particles are reconstructed from a material outside of the fracture process zone in the same specimen. In this set of simulations, the cut-off radius is 1:1 μm and the notch size is 15 μm. We first carry out simulations with and without microvoid shearing (modified GTN vs. classical GTN), referred to as C2-1 and C2-2 in Table 3. A series of snapshots on the development of the fracture surface at different levels of far-field loadings are shown in Figs. 10 (classical GTN) and 11 (modified GTN), respectively. In both simulations, the initial microvoid volume fraction is 0.1%. In Fig. 12, formation of the fracture surface is visualized by the accumulated effective plastic strain (0.2– 0.25). Similar to case 1 with 17 primary particles, the ductile fracture process starts with plastic deformation accumulation at the blunting crack tip (Fig. 10a). Subsequent primary void growth caused by multiple void interaction occurs first at the void ahead of the crack tip (Fig. 10b). Microvoid coalescence driven by the tension connects the primary voids and macro fracture occurs (Fig. 10c). The macro crack propagates forward due to subsequent primary void growth and coalescence (Fig. 10d). Finally, we see a similar fracture profile as that observed in Fig. 5a.

Fig. 10. Fracture process with 115 primary particles without considering microvoid shearing.

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2121

Fig. 11. Fracture process with 115 particles with microvoid shearing kω ¼ 15.

For case C2-2 which considers microvoid shearing (the modified GTN model), we set kω ¼ 15. As can be seen from Fig. 11, the ductile fracture process starts with plastic deformation accumulation at the blunting crack tip, similar to the case with the classical GTN model. We then can observe the primary void growth and shear band formation between the primary voids (Fig. 11a,b). The shear band orients at around 251 relative to the x-axis. The shear band grows and finally connects the primary voids and crack tip forming the macrocrack (Fig. 11c). Then, further microvoid coalescence driven by shear can connect the primary voids and propagate the crack (Fig. 11d). Similar to case 1 with 17 primary particles, growth of primary voids does not necessarily occur in the sequence of crack propagation; it is in a simultaneous and inter-competitive fashion as previously mentioned.

3.3. Comparison with experiments We now proceed to simulate our experiments. The position and size of primary particles are reconstructed from a material outside of fracture process zone in the same specimen. With a 1:1 μm cut-off radius, 132 particles are embedded into a simulation domain. The finite element mesh is shown in Fig. 12a. It is referred as C3 in Table 3. We set the parameters in the simulation the same as those used in experiments: the notch size of the fatigue pre-crack ð3 μmÞ, initial microvoid volume fraction f 0 ¼ 0:04% which is consistent with our experimental data and the thickness of domain in the simulation ð100 μmÞ (the thickness in the experimental observation is 350 μm. 100 μm is the largest thickness that our computation can afford). We observe the fracture profile at the unloading configuration. Because more voids are involved in the simulation, the computational cost is much higher in the current simulation cases than simulations shown in the previous sections. We section our simulation domain by a plane with normal direction in the z-axis. The section at z ¼ 1:5 μm is shown in Fig. 12b. We put the experiment and the simulation under the same length scale for comparison. It is obvious that the zig-zag fracture profile observed in experiments is reproduced. We can see that more zig-zag fracture patterns are observed in the experiments as compared to the simulation. However, the general non-flat fracture profile is consistent with experimental observations. In the black rectangular region, the zig-zag fracture pattern is nearly

2122

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

Fig. 12. Comparison of shear-driven fracture profiles with experiments: (a) Finite element mesh with 132 primary particles; (b) comparison of experimental fracture profile with simulation.

Fig. 13. (a) 3D fracture surface. (b) Comparison of simulation on crack opening displacement at unloading state with experiments (unit: μm).

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

2123

reproduced. In the simulation, shear bands between voids and primary voids are both observed in a similar location as that shown in the experiments. The three-dimensional fracture surface in our simulation at the unloading state is also shown in Fig. 13a. The crack path varying through the thickness showed that the zig-zag pattern is anisotropic. It is observed along the x-axis direction and not the z-axis direction. To quantitatively characterize the shear instability behavior, we plot average crack opening distance (COD) vs. x at the unloading state for both the experiment and the simulation by taking the average of all COD values within a given x-coordinate and plotting it at the respective position on the x-axis. Results are shown in Fig. 13b. It can be seen that the simulation quantitatively agrees with the experiment to some extent. However, the discrepancy still exists. The difference may be explained by the lack of a real secondary population of particles (voids), which are smeared out and represented by the proposed multiresolution model. This may also be the reason for less zig-zag patterns on the fracture path shown in Fig. 12. The other fracture indicators such as the contour of the crack opening distance, the zig-zag wave length, the void growth ratio, the void size, the void density on the fracture surfaces, and the nearest neighbor spacing reasonably agree with experiments. A detailed comparison is reported in our recent work (Chan et al., 2013).

3.4. Summary of simulations Let us summarize the 10 predictions in Table 3. Some of the prediction cases are not shown in this paper due to limitation of space. Table 3 leads to two conclusions: 1. Crack opening distance and fracture toughness when full cracking takes place in the simulation domain are closely related to notch size and primary void size. Lower notch size leads to the lower initial toughness and crack opening distance. Crack lengths during full cracking are much smaller than that seen in the experiments except for case 2. When we employ the experimental parameters in our simulation, fracture indicators in the fracture process reasonably agree with experimental results. 2. Fracture toughness and primary void growth ratio vary in a wide range depending on the fracture mechanisms. When considering microvoid shearing, the initial fracture toughness is lower, which is more consistent with experimental data. When considering microvoid shearing, the transverse void growth ratio is larger than that in the loading direction. The reverse is true when neglecting microvoid shearing.

4. Concluding remarks In this paper, we proposed a hybrid multiresolution model to simulate ductile crack growth. The proposed hybrid multiresolution model attempts to improve upon the shortcomings of previous modeling approaches. Small scale yielding in finite element simulations of 3D crack growth on Mod 4330 were carried out to elucidate the relation between fracture profile and microstructure. We showed that the hybrid multiresolution model can capture the fracture profiles in our experimental observations in combination with LS-DYNA and large scale computation. The multiresolution framework provides a more realistic description of the microstructure of alloys during the fracture process and is computationally affordable. We demonstrated the deformation in the macro zone and the micro zone as well as the interaction between them and their effects on the fracture process. We also discussed mechanisms which lead to the zig-zag fracture profile. Benzerga et al. (2004a,b), Benzerga and Leblond (2010) and Morgeneyer and Besson (2011) provided a clue (anisotropy and shear void nucleation) for the zig-zag crack growth pattern. Through separating the contribution of primary voids and secondary voids along with large scale computing, we showed that microvoid shearing can lead to zig-zag crack growth directly. Compared to previous models, it is a simpler approach which is controlled by a single parameter kω . The comparison of the zig-zag pattern in our simulation to that of experiments has also been carried out. We showed that microvoid shearing can result in the zig-zag fracture profile which is consistent with our experimental observation (Chan et al., 2012) both qualitatively and quantitatively within a reasonable range. However, the classical models with only pure spherical microvoid growth cannot capture this behavior. The reproduced fracture pattern with a detailed analysis assisted in demonstrating that the present hybrid approach is sound and that the discussed mechanism is stable. The present work is the continuation effort of our group in modeling of high strength alloy Mod4330. The aim of the present work is to discuss how to simulate the zig-zag fracture pattern observed through experimentation and compare it with those experiments both qualitatively and quantitatively. Through the calibration provided by McVeigh and Liu (2008), the length scale parameter in the multiresolution model is found to be around 1 μm for this material. This value is used in the present simulation. Larger length scale can increase the width of the shear bands between voids, resulting in higher toughness. This point is discussed in McVeigh and Liu (2009) by utilizing a similar multiresolution framework. The quantitative study on effects of the length scale parameter with the present multiresolution theory is in progress and will be reported later. Also, the present work did not consider the nucleation of the secondary particles and how this phenomenon might influence the fracture profile. It is a topic worthy of further study.

2124

S. Tang et al. / J. Mech. Phys. Solids 61 (2013) 2108–2124

Acknowledgments The support of this work by DARPA-D3D program, NSF CMMI 0823327, and IDR-CMMI-1130948 is gratefully acknowledged. WKL was partially supported by the World Class University Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology. Discussion with Tobias Erhart and Thomas Borrval at LSTC Inc is greatly acknowledged. The authors also thank Miguel Bessa and Jacob Smith for proof reading and David Chen for assistance with the adjustment of large scale computing settings on the QUEST cluster at Northwestern University. References Aravas, N., 1987. On the numerical integration of a class of pressure-dependent plasticity models. Int. J. Numer. Methods Eng. 24, 1395–1416. Belytschko, T., Liu, W.K., Moran, B., 2000. Nonlinear Finite Elements for Continua and Structures. Wiley, Chichester. Benzerga, A.A., Besson, J., Pineau, A., 2004a. Anisotropic ductile fracture part II: experiments. Acta Mater. 52, 4639–4650. Benzerga, A.A., Besson, J., Pineau, A., 2004b. Anisotropic ductile fracture part II: theory. Acta Mater. 52, 4639–4650. Benzerga, A.A., Leblond, J., 2010. Ductile fracture by void growth to coalescence. Adv. Appl. Mech. 44, 169–305. Chan, S.C., 2011. Multi-scale Tomographic Analysis of Toughness and Fatigue Strength in UHS Steels. Ph. D. Dissertation, Northwestern University. Chan, S., Tang, S., Kopacz, A., Rowenhorst, D., Spanos, G., Liu, W.K., Olson, G.B., 2013. Multiscale ductile fracture integrating tomographic characterization and 3D simulation, submitted for publication. Cowie, J.G., Azrin, M., Olson, G.B., 1989. Microvoid formation during shear deformation of ultrahigh strength steels. Metall. Trans. 20A, 143–153. Fleck, N.A., Hutchinson, J.W., 1993. A phenomenological theory for strain gradient effects in plasticity. J. Mech. Phys. Solids 41, 1825–1857. Kanninen, M.F., Popelar, C.L., 1985. Advanced Fracture Mechanics. Oxford University Press, New York. Kim, J., Gao, X., Srivatsan, T.S., 2003. Modeling of crack growth in ductile solids: a three-dimensional analysis. Int. J. Solids Struct. 40, 7357–7374. Liu, W.K., McVeigh, C., 2008. Predictive multiscale theory for design of heterogeneous materials. Comput. Mech. 42, 147–170. LS-DYNA User's Manual, LSTC Inc., 1997. LS-DYNA, Theory Manual, LSTC Inc., 1997. McVeigh, C., Vernerey, F., Liu, W.K., Moran, B., Olson, G., 2007. An interactive micro-void shear localization mechanism in high strength steels. J. Mech. Phys. Solids 55, 225–244. McVeigh, C., Liu, W.K., 2008. Linking microstructure and properties through a predictive 3 multiresolution continuum. Comput. Methods Appl. Mech. Eng. 197, 3268–3290. McVeigh, C., Liu, W.K., 2009. Multiresolution modeling of ductile reinforced brittle composites. J. Mech. Phys. Solids 57, 244–267. McVeigh, C., Liu, W.K., 2010. Multiresolution continuum modeling of micro-void assisted dynamic adiabatic shear band propagation. J. Mech. Phys. Solids 58, 187–205. Morgeneyer, T.F., Besson, J., 2011. Flat to slant ductile fracture transition: tomography examination and simulations using shear-controlled void nucleation. Scr. Mater. 65, 1002–1005. Nahshon, K., Hutchinson, J.W., 2008. Modification of the Gurson model for shear failure. Eur. J. Mech. A/Solids 27, 1–17. Pardeon, T., Hutchinson, J.W., 2003. Micromechanics-based model for trends in toughness of ductile metals. Acta Mater. 51, 133–148. Petti, J.P., Dodds, R.H., 2005. Ductile tearing and discrete void effects on cleavage fracture under small-scale yielding conditions. Int. J. Solids Struct. 42, 3655–3676. Tang, S., Greene, S., Liu, W.K., 2012a. A two-scale theory for nonlinear viscoelasticity. J. Mech. Phys. Solids 60, 199–226. Tang, S., Greene, S., Liu, W.K., 2012b. Renormalization approach to model interaction in microstructured solids: application to porous elastomer. Comput. Methods Appl. Mech. Eng. 217–220, 213–225. Tang, S., Kopacz, A.M., Chan, S., Olson, G., Liu, W.K. Concurrent multiresolution finite element: formulation and algorithmic aspects. Comput. Mech., http:// dx.doi.org/10.1007/s00466-013-0874-3, in press. Tian, R., Chan, S., Tang, S., Kopacz, A.M., Wang, J.S., Jou, H.J., Siad, L., Lindgren, L.E., Olson, G.B., Liu, W.K., 2010. A multiresolution continuum simulation of the ductile fracture process. J. Mech. Phys. Solids 58, 1681–1700. Tvergaard, V., Hutchinson, J.W., 2002. Two mechanisms of ductile fracture: void by void growth versus multiple void interaction. Int. J. Solids Struct. 39, 3581–3597. Vernerey, F., Liu, W.K., Moran, B., 2007. Multi-scale micromorphic theory for hierarchical materials. J. Mech. Phys. Solids 55, 2603–2651. Vernerey, F., Liu, W.K., Moran, B., Olson, G.B., 2008. A micromorphic model for the multiple scale failure of heterogeneous materials. J. Mech. Phys. solids 56, 1320–1347. Weck, A., Segurado, J., LLorca, J., Wilkinson, D., Bohm, H., 2008. Numerical simulations of void linkage in model materials using a nonlocal ductile damage approximation. Int. J. Fract. 148, 205–219. Xia, L., Shih, C.F., 1995. Ductile crack growth-I. A numerical study using computational cells with microstructurally-based length scales. J. Mech. Phys. Solids 43, 233–259. Xue, L., 2008. Constitutive modeling of void shearing effects in ductile fracture of porous materials. Eng. Fract. Mech. 75, 3343–3366.