A multiresolution continuum simulation of the ductile fracture process

Journal of the Mechanics and Physics of Solids 58 (2010) 1681–1700

Contents lists available at ScienceDirect

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

A multiresolution continuum simulation of the ductile fracture process Rong Tian a,1, Stephanie Chan b, Shan Tang a, Adrian M. Kopacz a, Jian-Sheng Wang c, HerngJeng Jou c, Larbi Siad d, Lars-Erik Lindgren e, Gregory B. Olson b, Wing Kam Liu a,n a

Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Rd., Evanston, IL 60208, USA Department of Materials Science and Engineering, Northwestern University, 2220 Campus Dr., Evanston, IL 60208, USA QuesTek Innovations LLC, 1820 Ridge Ave., Evanston, IL 60201, USA d Universite´ de Reims, UFR Sciences Exactes et Naturelles, 51687 Reims Cedex 2, France e Division of Material Mechanics, Lulea University of Technology, 971 87 Lulea, Sweden b c

a r t i c l e in fo

abstract

Article history: Received 16 November 2009 Received in revised form 7 June 2010 Accepted 1 July 2010

With the advancement in computational science that is stepping into the Exascale era and experimental techniques that enable rapid reconstruction of the 3D microstructure, quantitative microstructure simulations at an unprecedented fidelity level are giving rise to new possibilities for linking microstructure to property. This paper presents recent advances in 3D computational modeling of ductile fracture in high toughness steels. Ductile fracture involves several concurrent and mutually interactive mechanisms at multiple length scales of microstructure. With serial sectioning tomographic techniques, a digital dataset of microstructure features associated with the fracture process has been experimentally reconstructed. In this study, primary particles are accurately and explicitly modeled while the secondary particles are modeled by a two scale multiresolution continuum model. The present numerical simulation captures detailed characteristics of the fracture process, such as zigzag crack morphology, critical void growth ratios, local stress triaxiality variation, and intervoid ligament structure. For the first time, fracture toughness is linked to multiscale microstructures in a realistic large 3D model. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Multiresolution microstructure mechanics Multiscale constitutive equations Fracture toughness and tomography Parallel computing Finite elements

1. Introduction A severe limitation of existing modeling and CAE systems is the inability to seamlessly address the multiscale nature of the material/structural responses required to predict structural response in transient extreme environments. As depicted in Fig. 1, the damage evolution at a critical location (hot spot) of an aerospace structural system, which is made of twophase (for strength and toughness) alloys or composites, is a combination of microstructural events at multiple length scales. In this illustration, the macroscopic stress initiates the nucleation and growth of voids, creating an intense strain field between voids. Within the microstrain field, microvoids nucleate on finer secondary particles and microvoid growth

n

Corresponding author: Fax: + 1 847 491 3915. E-mail address: [email protected] (W.K. Liu). 1 Current address: Institute of Computing Technology, Chinese Academy of Sciences, No. 6 Kexueyuan South Rd., Zhongguancun, Haidian District, Beijing 100109, PR China. 0022-5096/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jmps.2010.07.002

1682

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

hot spot structure

Initial microstructures of alloy or reinforced composites

Nucleation and growth of large voids. Creation of an intense strain field between voids

Microvoid nucleation, growth and coalescence

Final fracture

Fig. 1. Damage evolution of two scales of microstructure and its interacting strain fields at a material point of the structural system.

Fig. 2. High accuracy and high resolution simulation based on experimentally reconstructed microstructure data and large-scale parallel computing.

and coalescence occurs on a finer scale. Interactions among these microstructure features eventually lead to final macroscopic ductile fracture. Currently, advanced 3D tomographic experimental tools such as the local-electron atom-probe (LEAP), focused ion beam (FIB), and scanning electron microscope (SEM) are readily available for characterizing microstructures in the ductile fracture process. The experimentally reconstructed microstructures can then be directly embedded in a numerical model to study microstructure evolution and its interaction with macroscopic properties of materials. In particular, with the help of large-scale parallel computing, the high resolution microstructure-level continuum simulation would give rise to new possibilities for linking microstructure to properties. A computational view of linking microstructure to properties is illustrated in Fig. 2.

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

1683

This work is part of the on-going ONR/DARPA Digital 3D (‘‘D3D’’) Structure program (Christodoulou, 2009), which is ultimately for the design of high performance steels and titanium alloys. Through this effort, a multiresolution continuum theory (Kadowaki and Liu, 2004, 2005; Hao et al., 2004a, 2004b, 2004c; Vernerey et al., 2006, 2007, 2008; McVeigh and Liu, 2008, 2009; Vernerey, 2006; McVeigh, 2007) has been developed. The theory overcomes one of the fundamental shortcomings of conventional continuum mechanics, that is, the inability to predict the relationship between microstructure and properties and the corresponding microstructure evolution. Based on the theory, multiresolution modeling with a finite element method has been applied to a high strength steel alloy with two scales of embedded particles (Vernerey et al., 2006, 2007, 2008; Vernerey, 2006; McVeigh, 2007) and cemented carbide (WC-Co) which undergoes three distinct scales of inhomogeneous deformation (McVeigh and Liu, 2008, 2009; Liu et al., 2009; Liu and McVeigh, 2008; McVeigh, 2007). Hao et al. (2004a, 2004b, 2004c) also incorporated quantum mechanical analysis in a multiresolution continuum simulation for the calculation of the interfacial strength between microstructure and matrix. Yin et al. (2008) applied the theory to multiscale material optimization. Tian et al. (2008, 2009) extended the theory to 3D. Based on the theory, a parallel multiscale fracture simulator was developed (Tian et al., 2008). Recently, we have extracted a digital dataset of the microstructure of the material of interest (a Ti-modified 4330 steel) and reconstruct the material microstructure on multiple length scales. On the sub-micron level, focused ion beam milling and scanning electron microscope imaging have captured sub-micron carbides and microvoids in the shear band of a steel sample. On the micron level, metallographic polishing, and optical microscope imaging have enabled reconstruction of a crack tip process zone, including the voids and inclusions in the vicinity. In this paper, we combine the experimentally reconstructed microstructure dataset, the developed 3D multiresolution continuum theory, and large-scale parallel computing in a high resolution modeling to simulate of the 3D fracture process of the Ti-modified 4330 steel. With the experimental data, we are able to accurately model the microstructure evolution and validate our simulations. The paper is arranged as follows. In Section 2, we give a brief review of the multiresolution continuum theory followed by the proposed variable length-scale theory. This two-length-scale multiresolution continuum theory for microstructured materials is employed. In Section 3, we describe the experimental reconstruction of material microstructures. In Section 4, we review various approaches to ductile fracture modeling. In Section 5, we discuss the computational framework, presenting a simulation framework for the three-dimensional ductile fracture process zone based on the experimentally reconstructed microstructures. In Section 6, we discuss the simulation results and comparisons with experiments. In Section 7, conclusions are drawn, followed by some discussion on future directions and applications in Section 8. 2. Review of the multiresolution continuum theory and development of a variable-length scale theory Traditional continuum mechanics theories focused on approximating the heterogeneous microstructure as a homogeneous continuum, which is conducive to a FEM mathematical description. Although this makes large-scale simulation of material much more efficient than modeling the detailed microstructure, the relationship between microstructure and macroscale properties becomes unclear. 2.1. Theoretical details The multiresolution theory is developed in terms of the virtual power in a body during deformation. Separate contributions to the virtual internal power are assumed to arise from homogeneous and inhomogeneous deformation. Background information can be found in the references (Hill, 1963, 1972; Nemat-Nasser and Hori, 1993; Aifantis, 1992), and in particular, those recently developed by our group (McVeigh et al., 2006; Vernerey et al., 2006, 2007, 2008; McVeigh and Liu, 2008, 2009; Vernerey, 2006; McVeigh, 2007). Continuum mechanics is based on the premise that a heterogeneous body can be approximated as a smooth continuous medium. Mathematically, this is known as homogenization: the variational internal power density dp(x) at any point in the continuum is approximated as the average power density of a superimposed representative volume element (RVE) V0 of the actual microstructure. The resulting continuum virtual power density field dp(x) is written as Z Z 1 1 dpðxÞ ¼ dpm dV0 or dpðxÞ ¼ rm : dLm dV0 ð1Þ V0 V0 V 0 V0 where V0 is the volume of the RVE superimposed at the continuum point x. The internal virtual power can be rewritten in terms of the power conjugate local stress and velocity gradient fields within the RVE, rm, and Lm respectively. Throughout this paper the subscript m refers to local fields which exist within the microstructure (RVE). Using the Hill–Mandel lemma (Hill, 1963) Z 1 dpðxÞ ¼ rm : dLm dV0 ¼ rðxÞ : dLðxÞ; dphom ðxÞ ¼ rðxÞ : dLðxÞ ð2Þ V0 V0 where r and L are the average stress and velocity gradient within the RVE defined by Z Z 1 1 rðxÞ ¼ rm dV0 ; LðxÞ ¼ Lm dV0 V0 V0 V0 V0

ð3Þ

1684

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

In many materials, deformation may become very inhomogeneous at or below the RVE scale, often resulting in terminal strain localization. The conventional homogenization approach given in Eqs. (2) and (3) is fundamentally incapable of representing deformation at these smaller scales. Hence we extend the Hill–Mandel approach to account for both homogeneous and inhomogeneous contributions such that dp at a material point (x) becomes Z Z 1 1 dpðxÞ ¼ dphom ðxÞ þ dpinh ðxÞ, dpinh ðxÞ ¼ dpm dV1  dpm dV0 ð4Þ V1 V1 V0 V0 The inhomogeneous contribution, dpinh, at a continuum point is defined as the difference between the average virtual power density at the inhomogeneous deformation scale V1 and the average virtual power density at the RVE scale V0. It can be shown as (V1 oV0 at the same x) Z 1 dpinh ðxÞ ¼ ðrm : dLm ÞdV1 rðxÞ : dLðxÞ V 1 V1 Z 1 ¼ ðrm : dLm rðxÞ : dLðxÞÞdV1 ð5Þ V1 V1 It is useful to rewrite the inhomogeneous power density dpinh directly in terms of the local inhomogeneous velocity gradient Lm  L Z 1 ðrm : dLm rðxÞ : dLðxÞÞdV1 dpinh ðxÞ ¼ V 1 V1 Z 1 ¼ b : ðdLm dLðxÞÞdV1 ð6Þ V1 V1 m where the proposed power equivalence relationship

bm : ðLm LðxÞÞ ¼ rm : Lm rðxÞ : LðxÞ

ð7Þ

has been used. This introduces a local microstress bm as a power conjugate of the local inhomogeneous velocity gradient Lm  L. To account for the locally varying nature of the inhomogeneous deformation, it is assumed here that the local velocity gradient Lm varies linearly at scale V1 (Gao et al., 1999)  Z  Z 1 1 Lm dV þ Lm r dV Uy ¼ L1 ðxÞ þG1 ðxÞUy ð8Þ Lm  V1 V1 V1 V1 where y is the local position with respect to the center of V1. The variables L1 and G1 are continuum measures corresponding to volume averages of the local velocity gradient Lm and gradient of the local velocity gradient Lm r , respectively. Substituting Eq. (8) into Eq. (6), the inhomogeneous virtual internal power density dpinh(x) can be shown to be  Z   Z  1 1 dpinh ðxÞ ¼ bm dV1 ðdL1 dLÞ þ bm  ydV1 dG1 ð9Þ V1 V V1 V This expression contains volume averages of the local microstress bm and its first moment bm  y at the scale of the inhomogeneous deformation, V1. This is clearer when the virtual internal power density is rewritten as   1 dpinh ðxÞ ¼ b1 ðxÞ : dL1 ðxÞdLðxÞ þ b ðxÞ^dG1 ðxÞ ð10Þ 1

where the continuum microstress b1 and the higher order microstress b are defined as Z Z 1 1 1 b1 ðxÞ ¼ bm dV1 ; b ðxÞ ¼ b  ydV1 V1 V1 V1 V1 m

ð11Þ

The total virtual power density is the sum of the homogeneous contribution and the inhomogeneous contribution 1

dp ¼ r : dL þ b1 : ðdL1 dLÞ þ b ^dG1

ð12Þ 1

Each of the continuum measures r, b1 , b , L, L1 , G1 which describe the virtual internal power at a continuum point are written either directly or indirectly in terms of volume averages of the known local fields rm, Lm within the superimposed RVE. The key to this theory is that several length scales, given by the characteristic sizes of the averaging n

volumes Vn, are embedded directly via the microstress couple b . Although other higher order theories introduce stress n

couples, the difference here is that a microstress couple b is averaged directly at each scale of inhomogeneous deformation. This naturally introduces length scales ln, which are equal to the width of the cube averaging volumes Vn. 2.2. A variable length scale multiresolution microstress approach The multiresolution continuum theory can be implemented in two numerical approaches. The first approach, called the N-microstress approach, introduces a microstress and a microstress couple at each scale of inhomogeneous deformation such that the inhomogeneous contribution to the internal power density is a summation over each scale (at a material point x,

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

1685

there are N nested volumes, Vn oyoV1 oV0),

dpint ¼ r : dL þ

N X

n

ðbn : ðdLn dLÞ þ b ^dGn Þ

n¼1

Z 1 b dV V n Vn m Z n 1 b ¼ b  ydV V n Vn m

bn ¼

ð13Þ

In many materials, inhomogeneous deformation may occur at successively smaller distinct scales. The alternative approach, called the variable length scale multiresolution microstresses approach, employs alternatively a single ms multiresolution microstress bms and microstress couple b so as to represent the inhomogeneous nature of material behavior at multiple scales. The averaging volumes describe the scales of the inhomogeneous deformation. Inhomogeneous deformation tends to occur at one scale at a time, first at scale V0 followed by scale V1, etc. A single averaging volume can be used which evolves with the scale of the inhomogeneous deformation. Hence only a single microstress and microstress couple are required such that ms

dpint ¼ r : dL þ bms : ðdLms dLÞ þ b ^dGms Z

1 b dV Vms Vms m Z ms 1 b ¼ b  ydV Vms Vms m

bms ¼

ð14Þ

The averaging volume Vms contracts as the scale of inhomogeneous deformation decreases. This approach reduces the number of microstresses required, and therefore the number of degrees of freedom to be solved for. The principle of virtual power can be employed along with the divergence theorem to derive the multiresolution continuum governing equations (the superscript ‘‘ms’’ is dropped) ðrbÞU r þb ¼ 0

in O

bU r b ¼ 0 in O ðrbÞUn ¼ t on Gt bUn ¼ 0 on Gt

ð15Þ

The detailed derivation can be found in Vernerey (2006) and McVeigh (2007), McVeigh and Liu (2010). A key advance is that the technique for deriving these governing equations from the standpoint of nested RVE modeling gives rise to a natural and systematic framework for deriving microstress and microstress couple constitutive relationships (Fig. 3).

Void Nucleation

Void Growth and Coalesence

Alloy + Secondary particles Interfacial Cohesive Behavior

Normal Separation(nm)

Equivalent Stress (Pa)

Normal Stress (GPa)

Microvoiding matrix + Primary particles

2.5x109 2

Φ=Φ

1.5

0.5

0 0.2 0.4 0.6 0.8 1 1.2 Engineering Strain

Q: Vector of internal state variables

Φ1 = Φ1 (β1, β1, Q))

Subatomic scale Courtesy of Prof. Freeman (Northwestern University)

Few nanometers

( σ,, Q )

1

Φ2 = Φ2 (β2, β2, Q) 10-100 nanometers

Few microns

Fig. 3. Multiresolution continuum modeling with multiscale constitutive relations.

1686

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

3. Experimental reconstructions of material microstructure inside a fracture process zone Experimental reconstructions of microstructure were obtained at the sub-micron and micron length scales. Both experiments involved a Ti-modified 4330 steel, here designated ‘‘Mod4330’’. The work explained here is an on-going effort of the ‘‘D3D’’ project and will be published shortly in more detail (Chan, 2010). The following is a brief description only of the experimental serial sectioning methods employed. For observation of microstructure on the sub-micron scale, first, a double linear shear test was performed. Specimens in the form of TEM discs were then sectioned from the shear band of the gauge section and from the undeformed bar end of the specimen. A cantilever bar set-up was employed, using a focused ion beam (FIB) to mill out a ‘‘U’’ shape to prevent redeposition onto the plane of interest. Then, serial sectioning was completed with the FEI Auto Slice & View software. This consists of an alternating cycle between material removal and imaging. The FIB is used to mill a cross-section from the cantilever, while the SEM is used to image the surface of interest. The method is explained in further detail in Uchic et al. (2006). Utilizing a combination of automatic and manual segmentation methods, various microstructure features were captured, including primary inclusions (TiN and CaS), primary oxides, primary voids/cracks, secondary carbides (TiC), and microvoids. For the micron scale data, a JIC fracture toughness test was performed in order to obtain an appropriate specimen giving insight into the ductile fracture process. Compact tension specimens were prepared and tested according to ASTM E1820; however, the test was slightly modified in that the crack was loaded only to the point of the onset of crack instability. At this point, the presence of void growth and void coalescence in the crack tip process zone is apparent. Serial sectioning was completed with mechanical polishing and light optical microscopy techniques. Automatic and manual segmentations highlighted the various microstructure features, such as the crack, including the primary voids and shear bands, and the inclusion (TiN) distribution (Chan, 2009). These reconstructed data, shown in Fig. 4, provided statistical information for microstructure generation in this study. For each feature, location of centers, corresponding sizes, and size distributions served as the input for recreating a comparable microstructure computational/mesh model at the start of the simulations. 4. Review of ductile fracture modeling Ductile fracture occurs through void nucleation, void growth, and void coalescence. Modeling of ductile fracture has been documented in the literature in a long period of time (Rice and Tracey, 1969; Rice and Johnson, 1970; Tvergaard and Needleman, 1987; Tvergaard and Hutchinson, 1992; Liu et al., 1996; Wang and Olson, 2009). Ductile fracture is a multistage process involving several interactive mechanisms at various length scales. Many researchers used effective computational models of ductile fracture, where the void containing cell elements are employed based on constitutive models such as Gurson (1977), Xia and Shih (1995), and Tvergaard and Hutchinson (2002). These simulations are capable of elucidating the trend of fracture toughness. However, these models have yet to be used to explore fundamental connections between multiscale microstructures and toughness.

Primary Inclusions Primary Voids Oxides Submicron Carbides Microvoids

crack

Shearing

Primary particles (yellow) near crack tip Reconstruction area: 633x516x259μm3

Secondary particles (grey dots) inside shear band Reconstruction area: 70x15x28μm3

Average size: 4.89μm Spacing: 19.17μm

Average size: 0.054μm Spacing: 0.66μm

Volume fraction: 0.080%

Volume fraction: 0.038% Fig. 4. Experimental reconstructions of microstructures.

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

1687

Ductile fracture

Void growth

Nucleation

Volumetric change Coalescence Multi-void interaction and void distortion (our focus)

Fig. 5. Micromechanics of ductile fracture and the recent focus.

Fig. 5 summarizes the micromechanics of ductile fracture and highlights recent research trends. Extensive studies have been focused on void nucleation and void growth, for example, Tvergaard (1981), Tvergaard and Needleman (1984), Needleman (1987), Needleman and Rice (1978), Needleman et al. (1992), Rice and Tracey (1969), McClintock et al. (1966), Lee and Mear (1992), Gologanu et al. (1993, 1994), McClintock (1968), and Thomason (1990, 1993). Modeling of multiple void interaction and void coalescence hence have received far less attention in the literature (Pardoen and Hutchinson, 2000). Progress in void coalescence study is mainly hampered by the lack of quantitative numerical and experimental results to assess the validity of theoretical models; thus void coalescence remains the least understood (Weck and Wilkinson, 2008). In contrast to void growth, void coalescence is a more complex phenomenon. (a) Void coalescence involves strong intervoid interaction, which is historically neglected by single-void analysis. While a singlevoid model may characterize void growth, it is obviously far from sufficient in studying void coalescence. Multiple-void interaction modeling is very complicated and a closed form solution is impractical. Early studies of void growth in ductile materials focused on the behavior of a single void in an infinite block of a rigid perfectly plastic material. It was found that both for a cylindrical void (McClintock, 1968) and a spherical void (Rice and Tracey, 1969), the rate of growth is strongly dependent on the level of hydrostatic tension such that ductile failure by coalescence would be promoted by a high level of stress triaxiality. The Rice–Tracey model (Rice and Tracey, 1969) predicts that the void growth rate is exponentially dependent on the level of stress triaxiality. The model was further developed by Tracey (1971) to account for material strain hardening effects on the void growth such that an increase in the strain hardening requires a larger stress triaxiality to maintain a fixed void growth rate. The Rice–Tracey-type model is based on a single void and does not take into account the change of the void shape and interactions between voids nor does it predict ultimate failure. What has been long lacking until very recently is the consideration of void shape change with low or no volume change (Lee and Mear, 1992; Gologanu et al., 1993, 1994; Hu et al., 1993; Pardoen and Hutchinson, 2000). The recent investigations of Gologanu et al. (1993, 1994) and Pardoen and Hutchinson (2000) are extensions of the Gurson–Tvergaard–Needleman (GTN) model (1977, 1984, 1987) to account for the relative void spacing and void shape effects. (b) Void coalescence is an unstable void growth stage; void distortion dominates this deformation stage. The most widely employed criterion for the onset of void coalescence states that void coalescence starts at a critical porosity, which has been regarded tentatively as a material constant (McClintock, 1968; d’Escatha and Devaux, 1979). In the widely used Gurson–Tvergaard–Needleman (GTN) model (1977, 1984, 1987), void coalescence is modeled by accelerating void growth at a phenomenological critical void volume fraction, fc. In the work of McClintock (1968), Thomason (1990, 1993), and Pardoen and Hutchinson (2000), it has been recognized that a microstructural model of void coalescence definitely requires the introduction of information related to the void/ligament dimensions and geometry. (c) Void coalescence poses challenges to both experimental measurements and computer simulations. Void coalescence can be triggered by intervoid ligament necking and/or intervoid ligament shearing, depending on void spatial distribution, which is always accompanied with strain localization and local material softening (McVeigh et al., 2007). Numerical simulation suffers from mesh dependence. Experimental measurements of void coalescence are difficult to perform. (d) Knowledge of the underlying mechanisms of softening, localization and fracture in shear is more qualitative than quantitative. The mixed mode ductile fracture is still not fully understood. In order to link the microstructure to the toughness (Olson 1997, 2000), some key questions have to be considered, such as: (a) What is the critical void growth ratio? The final void growth ratio relative to inclusion size is an experimentally measurable parameter reflecting the resistance to void coalescence.

1688

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

(b) What is the critical strain at intervoid ligaments? The critical strain is a key parameter in ductile fracture toughness, which dominates the final void growth ratio and is closely related to the strength and plasticity of the material and the size and density of microvoid-generating secondary particles. (c) How does the local stress triaxiality evolve? The local stress triaxiality is the driving force for the void growth and influences the fracture mode mixity and mode transition during the fracture process. (d) Is the scaling role of primary void 3D spacing, X0, valid? What determines zigzag wavelength caused by shear localization and coalescence? Can it be independent of X0? Based on mostly 2D experimental data available in the literature, a phenomenological toughness model, which connects the ductile fracture toughness with the microstructure and the yield strength, has been developed (Tvergaard and Hutchinson, 1992). The model answers some of these questions empirically. Next, a large-scale multiresolution continuum approach is developed to model the multiple void growth and coalescence with accuracy.

5. Computational model 5.1. Domain reduction As shown in Fig. 6, the domain of interest contains three nested subdomains: the outer K-field, the plastic zone, and the inner process zone. The dimensions of the process zone are assumed small (several hundred mm according to the experimental reconstruction) compared to the characteristic dimensions of the K-field, and any plasticity is confined to the plastic zone. Inside the numerical process zone, experimentally reconstructed microstructures are embedded directly. Under the small scale yielding assumption the displacement field along the outer boundary of the K-field can be calculated using the equations of linear elastic fracture mechanics (LEFM) (Kanninen and Popelar, 1985) rffiffiffiffiffiffi     KI R y y ux ¼ k1 þ 2sin2 cos 2m 2p 2 2 rffiffiffiffiffiffi     KI R y y uy ¼ k12 cos2 sin 2m 2p 2 2 uz ¼ 0 ð16Þ where R and y are the distance and the angle from the crack tip, respectively, m is the shear modulus, k= 3 4v for plane strain, v is the Poisson ratio.

2. A full scale model makes it possible to directly link to macroscopic properties (toughness, etc)

K-field Plastic zone

y

Process zone (PZ) x

1. Key physics of secondary particles is represented

z

Experimental reconstruction of primary particles

Fig. 6. A simulation framework of 3D ductile fracture process zone.

A homogenization model built from RVE modeling for secondary particles

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

1689

The load is implemented through increasing the stress intensity factor KI in Eq. (16) from 0 to the experimentally measured fracture toughness KIC. The maximum KI is chosen to be 1.5KIC so that the fracture process is developed and the microstructure kinematics is directly linked to fracture toughness. 5.2. Primary void coalescence and the consideration of secondary voids In the Mod4330 steel examined here, there are two potential populations of voids, which exist at two distinct scales. Based on experimental data, we consider two length scales: primary particle scale, l1 = 5 mm, and secondary particle scale, l2 =0.1 mm. The length scale l1 is resolved directly by mesh discretization (and hence large-scale parallel computing) so that primary void growth and coalescence can be modeled with accuracy. Secondary particles’ effects (a void nucleated from secondary particles is called as microvoid) are modeled by a two-scale multiresolution continuum theory, as mentioned in Eq. (15). In this present work, as the damage due to the primary void growth is modeled explicitly, the macro-stress, r, now describes the deformation of microvoiding matrix at length scale l1 (refer to Fig. 3), a damage model is employed to account for microvoiding effects; the microstress, b, describes the deformation of alloy and secondary particles (refer to Fig. 3), which reflects the resistance to inhomogeneous deformation at length scale l2, a linear elastic constitutive law is employed for the microstress; the microstress couple, b, automatically brings the length parameter l2 into play. It is this characteristic length that makes the multiresolution continuum simulation free of the known mesh dependence issue showing up in classical continuum simulation of strain-softening. The microstress and the microstress couple are a measure of the resistance of strain localization at the ligaments between secondary voids or microvoids. The saturation of the microstresses signals the loss of the resistance capability of the matrix material and implies the onset of local microvoid coalescence. The local microvoid coalescence accumulates quickly and accelerates the primary void coalescence (Fabregue´ and Pardoen, 2008). The saturation of the microstress is equivalently controlled by the damage model. Traiviratana et al. (2008) studied the growth of nanosized voids in steels through atomistic simulation, and compared the void growth with the Cocks–Ashby damage model (Cocks and Ashby, 1980) and the Gurson damage model, respectively. They concluded that significant agreement is found with the Cocks–Ashby model. The nanosized voids are the same in size as the secondary particles of the Mod4330 steel investigated in the paper. The Cocks–Ashby damage model has also been used to model void growth in metals by Bammann et al. (1993) and Konke (1997). Motivated by these previous studies, the Cocks–Ashby damage model is finally chosen. A merit of this damage model is that it contains only one parameter, i.e., the power law creep exponent. The parameter is calibrated by fitting to the macro-fracture strain (which is 1.04) for the Mod4330 steel. It is noteworthy to point out that this is the only phenomenological parameter in the computational model used in our study. 5.3. Finite element mesh A graded mesh is generated with 8-node hexahedrons, which can alleviate mesh distortion to a certain extent. The mesh in the K-field is the coarsest while the finest mesh is in the numerical process zone. We assume that all the primary particles inside the numerical process zone are debonded from the surrounded matrix material prior to loading. Based on this assumption primary particles are embedded into the mesh as voids by element deletion. While this assumption makes mesh generation simple, the interfacial effect, contact between particles and the matrix, is neglected, i.e. void nucleation is not considered in the model. The total number of dofs is 21,761,775. The mesh results in an average void size of 5 mm, which corresponds to the average particle size. 5.4. Implementations The multiresolution continuum theory requires the support of an arbitrary number of dofs per node. The theory may be implemented by meshless techniques (Liu et al., 1995; Li and Liu, 2004). A completely new and 3D finite element implementation (in C++ and MPI) of the multiresolution continuum theory is made possible by a generalized finite element framework (Tian and Yagawa, 2007; To et al., 2008a, 2008b; Tian and Yagawa, 2005; Tian et al., 2006a, 2006b; Tian, 2006; Tian and Yagawa, 2006; Nishida et al., 2007). A total Lagrangian formulation is adopted for large deformation analyses. The Jaumann rate is used for objective stress updating (Belytschko et al., 2000). The 3D simulation code is developed in a way that can be automatically reduced back to a classical continuum code by turning off the length scale parameters (Tian et al. 2008, 2009). This feature makes the comparison between multiresolution continuum modeling and the classical counterpart straightforward, and it has been used in validating our model against commercial codes (Tian et al., 2008). The 3D code has been validated by the previous 2D results (McVeigh, 2007) in an adiabatic shear band simulation (Tian et al., 2008). The central difference method is chosen to solve momentum equations because of intense local material softening and strain localization in the problem under investigation. As expected, the explicit method offers very good scalability in parallel tests (on BlueGene/L using up to 108 elements) and demonstrates robustness in dealing with strain softening and highly nonlinear material deformation.

1690

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

6. Results and discussions With the aid of the above computational model, we perform a large-scale 3D multiresolution continuum simulation of the ductile fracture process with microstructural details. We focus on the study of microstructural evolution during ductile fracture. The experimentally reconstructed microstructural dataset is directly embedded into the FE mesh model. Simulated process zone tomography is compared with experimental measurements. Microstructural deformation is linked to fracture toughness in the end. We first outline the simulation results in Fig. 7: Fig. 7a shows an optical micrograph of the experimental compacttension specimen; Fig. 7b is the 3D tomographic reconstruction of the process zone of the crack tip; Fig. 7c shows simulation results of the process zone; and Fig. 7d is the simulation results of microstructure evolution inside the process zone. For the first time we have accomplished a simulation of the fracture process zone directly based on 3D microstructures that are experimentally obtained. Detailed discussions follow. Fig. 8 shows an overview of the simulated fracture process zone from the front side (Fig. 8a), the back side (Fig. 8b), and the top side (Fig. 8c), respectively. Irregularly shaped particles show the location of TiN clusters (primary particles). The open notch in the figure is the fatigue crack tip of the compact-tension specimen. The fracture surfaces visualized are plotted with a damage parameter, D = 0.9. On the fracture surface plots, strain rate is superposed to visualize the propagating front of microvoiding. Fig. 9 shows serial snapshots of the development of the fracture process zone in time series. The typical micromechanism of ductile fracture is repeated: the ductile fracture process starts with plastic deformation accumulation at the blunting crack tip; it is followed by circumferential microvoiding triggered by accumulated plastic deformation; and, it ends with a change in void aspect ratio and coalescence of primary voids. The red strain rate contours reveal that microvoiding does not necessarily occur in the sequence of crack propagation; primary voids grow in a simultaneous and inter-competitive fashion. The key ductile fracture processes, void growth (Fig. 10a) and void coalescence (Fig. 10b), are shown in detail. Qualitative comparisons with experiments are performed and revealed physics of microstructured materials are discussed. First, as shown in Fig. 11, both experiments and simulations show a zigzag wave pattern of localization. By plotting the crack surface, we investigate the relation of void spacing and the wavelength of the zigzag localization in Fig. 12. The characteristic wavelength is measured to be around 75 mm. This length is a projected distance of void spacing on the fracture surface, which is larger than the experimentally measured 3D spacing distance 19 mm and is closer to the 2D

Fig. 7. Multiscale simulation of three-dimensional ductile fracture process. (a) Experimental specimen; (b) microstructure reconstruction at crack tip; (c) simulation results of process zone; (d) simulation results of microstructure deformation and evolution inside process zone.

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

1691

Fig. 8. Three-dimensional fracture process zone. Fracture surfaces are visualized by plotting with a damage parameter, D = 0.9. On the fracture surface plots, strain rate is superposed to visualize the front surface of micro-cracking and voiding.

Fig. 9. Three-dimensional simulation showing microstructural features during fracture process zone development. (a) Strain accumulation; (b) void growth; and (c, d) void coalescence.

spacing distance 105 mm. Fig. 12 confirms that the ductile fracture surface clearly follows the spatial distribution of the primary voids and hence the characteristic wavelength of the fracture surface should be closely related to fracture toughness. Second, a dimple structure of voiding that is observed in experiment is predicted in the simulation (Fig. 13). By referring to the original position of particles (close-ups A and B in Fig. 13b), the simulation suggests that the dimple structure originates from the microstructure. Another comparable feature between the experimental reconstruction and the simulation is the crack opening distance (COD) (Fig. 14). The experimental reconstruction area is 633 mm  259 mm in the loading plane. The numerical model size is 300 mm  50 mm accordingly. The maximum computed COD is 28 mm which is comparable to the experimental

1692

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

Fig. 10. (a) Void growth and (b) void coalescence in the fracture process.

FCT: Fatigue crack tip

FC FCT T

Process zone 318

259

a

Experimental Reconstruction (Mod4330, Serial Sectioned Image

Process zone 50 FC FCT T

205

b Simulation (Size: 205x180x50 µm3)

Size: 630x516x259 µm3) Fig. 11. Process zone: experimental reconstruction and simulation.

measurement range of 25–50 mm at full crack opening. While this value lies within the experimental range, the modeled maximum COD lies at the lower end of this region due to its limited size corresponding to the onset of crack extension. The simulation allows further investigation of details of the fracture process. First, void growth versus load step is plotted in Fig. 15. Void growth is measured up to the moment of the first void coalescence. Fig. 15a–d are respectively for voids #1, #2, #3, and #4, numbering from near the fatigue notch to the location far away from the tip. It can be observed that (a) void growth in the transverse direction is consistently larger than that in the load direction for all four voids up to the moment of the first void coalescence. (b) void growth trends are different for the transverse and loading directions; void growth in the transverse direction shows a slowing down whereas void growth in the loading direction does not. This is most marked for void #1 (Fig. 15a). (c) the numerical simulation predicts a maximum transverse void growth ratio of 3 upon the first void coalescence. (d) all void growth starts from around load step 11 and grow simultaneously in the region analyzed. Further discussions of these points are as follows.

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

1693

Fig. 12. Relation between the projected 3D void spacing and the crack surface wavelength. (a) Projected void spacing; (b) crack surface wavelength.

Fig. 13. Microstructural features. (a) Experiment; (b) simulation.

Point (a) reveals that before the first void coalescence transverse void growth dominates. The directionally rich void growth observed in the numerical simulation reveals 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 located along the lateral transverse direction (xz-plane according to Fig. 6), rather than the void ahead of it in the fracture direction. This observation well explains the fact that in both the experiment and the simulation, the void opening distance contour plot (Fig. 14) shows structural lining up in the transverse direction (from the bottom to the top of the images in Figs. 14a, b) as favored by the macroscopic plane strain constraint. Point (b) implies a void growth mode transition. Prior to the first void coalescence, i.e., the formation of the process zone, void deformation shows stable growth behavior. Approaching the point of the first void coalescence (Fig. 15), transverse void growth starts to slow down whereas the growth rate in the load direction remains the same or starts to increase. This reflects the start of void growth mode transition, from stable void growth to uniaxial stretching in association with intervoid ligament shear localization. Point (c) defines the critical void growth ratio. The numerical simulation predicts a maximum transverse void growth ratio of 3 upon the first void coalescence, which is comparable with experiments. Due to the void growth mode transition, void growth becomes an unstable uniaxial stretching in the later stage of the fracture process. At the point the crack advances, there is a rapid and unstable uniaxial stretching along the load direction whereby the void size in the loading direction can become very large. It is thus reasonable to observe a larger void growth ratio in the load direction than in the transverse direction. As the transverse void growth ratio is relatively much less sensitive to the momentum of measurement and easy to measure, consistent with standard fractographic practice, it is confirmed that the transverse void growth ratio should be taken as the measureable microstructural link to fracture toughness. Point (d) shows simultaneous void growth during the ductile fracture process, consistent with strong multiple void interactions (see also Fig. 19). Both McClintock (1968) and Thomason (1990, 1993) pointed out that any general void coalescence model requires the incorporation of at least some microstructural information related to the void/ligament dimensions and geometry. Despite

1694

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

Crack Opening Distance

127 μm Precrak

318 μm Crackfront

Process zone limit

Experiment (633x259 μm)

Opening Distance (unit: mm) -0.0133 0.000820 0.0149

Simulation (205x50 μm)

Opening distance: 28.2 μm

Fig. 14. Crack opening distance contours. Crack opening distance: 25–50 mm (experiment); 28 mm (simulation).

its importance, there is little information available for the intervoid ligament structure, however. We investigate this structure in numerical simulation. In numerical simulation, the structure of ligaments can be easily observed. Fig. 16a shows ligament walls represented by the plastic strain contour of 0.52; Fig. 16b shows the neutral surface of the process zone with zero displacement in load direction. It is seen that the constructed intervoid ligament structure is composed of all voids that contribute to the formation of the process zone. The physical implication of this intervoid ligament structure is that once this ‘‘structure’’ is broken the fracture process is complete. As such, the intervoid ligament is an appropriate region to focus on in the study of the mixed mode fracture process. At this region, we are able to investigate the critical plastic strain, the Lode parameter,2 and the stress triaxiality. First, local stress triaxiality is studied. Stress triaxiality versus plastic strain on the ligament structure is plotted in Fig. 17a. Stress triaxiality versus load is plotted in Fig. 17b. Stress triaxiality at intervoid ligaments varies during the ductile fracture process and shows a peak point of around 2.7 on average. Upon void coalescence, stress triaxiality descends rapidly. Second, the critical plastic strain at ligaments is studied. In the multiple void configuration, the progress of ductile fracture is attributed to the sequential (one after another) failure of intervoid ligaments. For example, when the stress softens at ligament #4, ligament #4 releases stress to its surrounding ligaments #3, #5, #6, and #7. Ligaments #3, #5, #6,

2 The Lode parameter is defined by m = (2s2  s1  s3)/(s1  s3), where s1 Z s2 Z s3 are the principal stress (Lode, 1925). The parameter reflects local fracture mode mixity and mode transition, whose effect is especially noticeable at low stress triaxiality (Barsoum and Faleskog, 2007).

3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0

load dir. transverse

1 11 21 31 41 51 61

Load step

Load step

4

2.5 2.0 R/R0

R/R0

3

2

1695

load dir. transverse

1 11 21 31 41 51 61

1 3.0 2.5 2.0 1.5 1.0 0.5 0.0

3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0

R/R0

R/R0

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

load dir. transverse

1.5 1.0 0.5 0.0

load dir. transverse

1 11 21 31 41 51 61

1 11 21 31 41 51 61

Load step

Load step

Fig. 15. Void growth. (a) Void #1; (b) void #2; (c) void #3; and (d) void #4.

Fig. 16. Structure of intervoid ligaments. (a) Ligament walls (represented by the plastic strain contour of 0.52); (b) crack surface of process zone (represented by the contour of zero displacement in load direction). nepsDot: plastic strain rate.

1696

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

4

4

3.5

1 2

3

3

4

5

7

3.5

8

3

2.5

Stress triaxiality

Stress triaxiality

6

2 1.5 1 0.5

2 1.5 1 0.5 0

0 -0.5

1 2 3 4 5 6 7 8

2.5

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

-0.5

-1

-1

-1.5

-1.5

0

0.1

0.2

Effective plastic strain at ligament

0.3

0.4

0.5

0.6

0.7

0.8

J/σ0/χ0

Fig. 17. The variation of the stress triaxiality at intervoid ligaments. (a) Stress triaxiality versus effective plastic strain at ligaments; (b) stress triaxiality versus remote load.

4 2 1.6

5 3

1

6

7 8

1 2

1.4

3

Normalized stres

1.2

4

1.0

5

0.8

6

0.6

7 8

0.4 0.2 0.0 0

0.5

1.5 1 Effective plastic strain

2

Fig. 18. Stress–strain curve at intervoid ligaments. Stress is normalized by initial yield stress (lines 3, 5, 7, and 8 coincide). The stresses are measured up to the first void coalescence occurrence. Stresses soften at ligaments 1, 2, 4, and 6. Ligaments 3, 5, 7, and 8 do not enter stress softening stage yet at the time of final measurement.

and #7 delay stress collapse and help to increase the stress resistance at ligament #4. This multiple void interaction, which results in the redistribution of stresses among ligaments, slows down the stress collapse at individual ligaments, e.g., at #4, and increases ligament failure strain. As seen from Fig. 18, ligament #4 fails at an equivalent plastic strain of around 1 (line 4 in Fig. 18). This interaction may also be seen from line 4 and line 6, two failed ligaments; neither show abrupt stress collapse at the local level. This observation implies that the stress redistribution amongst surrounding ligaments delays and slows down individual ligament failure and increases fracture resistance. The rapid stress collapse on the macroscopic scale is an accumulated effect of progressive (but slow) failure at intervoid ligaments. Hence, it is important to take interaction effects into account to model the critical strain of ligament failure. At the secondary void level, the critical ligament plastic strain is one of the key parameters for characterizing void coalescence. The internal structural balance of ligaments, especially in the lateral transverse direction (xz-plane in Fig. 6), cannot be simply justified by a 2D model. Single ligaments do not account for the structural balance, hence a planar two

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

ln(R/R0) (transverse direction)

2

3

10.0

4

void 1 (near crack tip) void 2 void 3 void 4 (away from tip)

1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 J/σ0/χ0

ln(R/R0) (load direction)

1

10.0

1

1697

2

3

4

void 1 (near crack tip) void 2 void 3 void 4 (away from tip)

1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 J/σ0/χ0

Fig. 19. Void growth ahead of the crack tip up to the first void coalescence. (a) Transverse growth; (b) growth in load direction. w0 is chosen to be the characteristic wavelength of crack surface as measured in Fig. 12.

void analysis may well underestimate the critical ligament plastic strain. Again, the importance of multiple void interaction effects is emphasized. Another interesting fact observed from Fig. 18 is that the critical strain at the ligament varies and is closely related to the thickness of ligament wall—the smaller the thickness is, the smaller the critical strain would be. The ligament wall thickness can be estimated from the undeformed ligament structure shown in Fig. 16. Ligament #4 has the smallest ligament thickness and correspondingly the lowest critical strain, which is followed by ligament #6, then ligament #1 and then ligament #2 in order. The ligament softens as shown in Fig. 18 also in the order of ligaments #4, #6, #1, and #2. The ligament thickness must be an important parameter in a void coalescence criterion. This observation confirms McClintock’s and Thomason’s early assertion that a general void coalescence model requires the incorporation of at least some microstructural information related to the void/ligament dimensions and geometry. In the current simulation set-up, the controlling softening mechanism of the ligament starts with high triaxiality ligament necking followed by low triaxiality ligament shearing. From Fig. 17a, the stress triaxialities at the ligaments are all larger than 1 before the ligament strain reaches 1 (the maximum is 1.3). After that, the remaining plastic strain is from ligament deformation under a stress triaxiality lower than 1. Also from line 4 in Fig. 17b, it may be seen that the shearing process is short and rapid in time, terminating in ligament failure. Finally, the relation between the normalized J-integral, J/(s0w0) and void growth ratio R/R0 is studied in Fig. 19, where s0 is the initial yield strength, w0 is chosen to be the characteristic length of crack surface wave as investigated in Fig. 12, and R/R0 is the void growth ratio in either transverse or load direction. The transverse void growth follows (Fig. 19a) a concavedown type of trend whereas the load-direction void growth (Fig. 19b) follows a concave-up type of trend. That is, the void growth in the two directions follows different trends. In a 2D context, void growth has already been reported in literature, for example, Tvergaard and Hutchinson (1992). However, to the authors’ knowledge, study of multiple void growth in a 3D context is rare. For the first time void growth in the transverse direction is investigated by the 3D model and a different void growth trend is identified. Shown in Fig. 19a, when approaching the critical point, i.e., J/(s0w0) reaches its maximum value, the void growth in the transverse direction starts to slow down, reflecting a uniaxial stretching deformation mechanism (i.e., void distortion) of void growth during crack opening along the load direction. 7. Conclusions and remarks A theoretical and computational framework has been laid out for directly linking microstructural evolution to macroscopic properties (toughness) for an ultra-high strength steel, ‘‘Mod4330’’, with particles of two different length scales. For the first time, we have accomplished a 3D computer simulation of multiple void growth and coalescence with accuracy in a realistic microstructure configuration. A preliminary understanding of the fracture process of the high strength steel has been reached. The simulation captures in detail the fracture process characteristics, such as a localized zigzag pattern, critical void growth ratio, local peak stress triaxiality, and intervoid ligament structure. Some key observations are as follows:

 Both experiments and simulations show a localized zigzag wave pattern of the fracture surface associated with microvoid softening driven shear localization. The wavelength is comparable with the projected distance of void

1698

   





R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

spacing on the fracture surface, which is greater than the experimentally measured 3D spacing distance and closer to the 2D spacing. The numerical simulation confirmed that the ductile fracture path copes with the spatial distribution of primary voids and the characteristic wavelength of the fracture surface is closely related to the ductile fracture toughness. Stress triaxiality at intervoid ligaments varies during the ductile fracture process and shows a peak point of around 2.7 on average. For the first time void growth in the transverse direction is investigated by the 3D model. Different void growth trends are identified between the load and the transverse directions. The difference reflects a uniaxial stretching deformation mechanism (i.e., void distortion) of void growth during crack opening. Void growth curves obtained from numerical simulation show simultaneous void growth during the ductile fracture process. This simultaneous void growth fracture mechanism signals a strong multiple void interaction behavior during the fracture process. The simulation reveals that during the formation of the process zone, the transverse directional void growth dominates. Due to the rapid unstable growth in the loading direction, it is difficult to justify a critical void growth ratio in this direction. In contrast, the transverse void growth ratio is less sensitive to the momentum of measurement and should be used as the parameter linking to fracture toughness. In the multiple void configuration, the internal structural balance of ligaments delays individual ligament failure and increases fracture resistance. It is therefore important to take the interaction effects into account to model the critical strain of ligament failure. A planar two void analysis may well underestimate the critical strain due to the neglect of interaction effects. The structure of the intervoid ligament, on one hand, provides an appropriate place to investigate critical plastic strain, the Lode parameter, and stress triaxiality so that we can study local critical strain, local mixed mode fracture, and mode transition during a complicated mixed mode fracture process. On the other hand, the void/ligament dimensions and geometry also provide microstructural information for developing a general void coalescence model.

Finally, we remark that though we believe that the computer simulation has provided new insights to the fracture process of a crack tip, the results given in this paper are new and preliminary. Further development of the multiresolution materials laws and calibrations with experiments are needed before we can conclude with a quantitative understanding of the fracture process mechanisms of this material. 8. Future direction and application Conventional procedures for life prediction of components subjected to fatigue are generally based on the ‘‘safe-life’’ approach. The ‘‘damage-tolerant’’ approach on the other hand is often a more economical and suitable alternative for life prediction, especially when the rate of damage accumulation is well understood and can be monitored. Unfortunately traditional crack growth models do not accurately predict crack initiation or describe crack growth at the submicrostructure scales. There is no mechanical science-based continuum theory to model various stages of damage evolution (or crack initiation) from sub-micron scale to micron scale to sub-millimeter range. The majority of existing work is focused on the prediction of macrocracks, which are detectable using non-destructive evaluation (NDE). Since the damage nucleates and progresses on a much smaller length scale than can currently be tracked using conventional NDE tools, the development of a predictive model at these micro, sub-micro scales is of utmost importance. A significant challenge to be addressed in the multiscale damage prognosis is the fact that each manufactured part or component will have its own unique material microstructural fingerprint. For example, in the future when an aircraft wing panel is made, it would be automatically fingerprinted and these characterization dataset would be associated with the part number. Treating the digital dataset as the building block, the theoretical and computational framework developed can be used to model various stages of damage evolution (or crack initiation) from sub-micron scale to micron scale to submillimeter range. A multiscale damage prognosis/NDE technique then may be further engineered by establishing a database linkage between micro-damage evolution and macroscopic material/structural responses.

Acknowledgements The support of the ONR DARPA D3D Digital Structure Consortium (award N00014-05-C-0241), NSF CMMI 0823327, and Sandia (DoE) grants are gratefully acknowledged. References Aifantis, E.C., 1992. On the role of gradients in the localization of deformation and fracture. International Journal of Engineering Science 30 (10), 1279–1299. Bammann, B.J., Chiesa, M.L., Horstemeyer, M.F., Weingarten, L.I., 1993. In: Wierzbicki, T., Jones, N. (Eds.), Failure in Ductile Materials using Finite Element Methods, Structural Crashworthiness and Failure. Elsevier Applied Science.

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

1699

Barsoum, I., Faleskog, J., 2007. Rupture mechanisms in combined tension and shear—Experiments. International Journal of Solids and Structures 44, 1768–1786. Belytschko, T., Liu, W.K., Moran, B., 2000. In: Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, Ltd., New York. Chan, S.C., 2009. Multi-scale tomographic analysis of UHS steels for toughness and fatigue strength. D 3-D Annual PI Review, Arlington, VA, July 7–9. Chan, S.C., 2010. Multi-scale tomographic analysis of UHS steels for toughness and fatigue strength. Ph.D. Thesis. Northwestern University, Evanston, IL. Christodoulou, J., 2009. Dynamic 3-dimensional digital structure: a program review. Journal of Materials 61 (10), 21. Cocks, A.C.F., Ashby, M.F., 1980. Intergranular fracture during power-law creep under multiaxial stresses. Journal of Metal Science 14 (8), 395–402. d’Escatha, Y., Devaux, J.C., 1979. Numerical study of initiation, stable crack growth and maximum load with a ductile fracture criterion based on the growth of holes. In: Landes, J.D., Begley, J.A., Clarke, G.A. (Eds.), Elastic–Plastic Fracture. American Society for Testing and Materials, Philadelphia, pp. 229–248 (ASTM STP 668). Fabregue´, D., Pardoen, T., 2008. A constitutive model for elastoplastic solids containing primary and secondary voids. Journal of the Mechanics and Physics of Solids 56, 719–741. Gao, H., Huang, Y., Nix, W.D., Hutchinson, J.W., 1999. Mechanism-based strain gradient plasticity-I. Theory. Journal of the Mechanics and Physics of Solids 47 (6), 1239–1263. Gologanu, M., Leblond, J.B., Devaux, J., 1993. Approximate models for ductile metals containing nonspherical voids—case of axisymmetric prolate ellipsoidal cavities. Journal of the Mechanics and Physics of Solids 41, 1723–1754. Gologanu, M., Leblond, J.B., Devaux, J., 1994. Approximate models for ductile metals containing nonspherical voids—case of axisymmetric oblate ellipsoidal cavities. Journal of Engineering Materials and Technology 116, 290–297. Gurson, A.L., 1977. Continuum theory of ductile rupture by void nucleation and growth: Part I—yield criteria and flow rules for porous ductile media. Journal of Engineering Materials and Technology, 2–15. Hao, S., Liu, W.K., Moran, B., Vernerey, F.J., Olson, G.B., 2004a. Multi-scale constitutive model and computational framework for the design of ultra-high strength steels. Computer Methods in Applied Mechanics and Engineering 193, 1865–1908. Hao, S., Liu, W.K., Klein, P.A., Rosakis, A.J., 2004b. Modeling and simulation of intersonic crack growth. International Journal for Solids Structure 41, 1773–1799. Hao, S., Liu, W.K., Weertman, J., 2004c. Cohesive solutions of intersonic moving dislocations. Philosophical Magazine 84, 1067–1104. Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. Journal of the Mechanics and Physics of Solids 11 (5), 357–372. Hill, R., 1972. On constitutive macro-variables for heterogeneous solids at finite strain. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 326 (1565), 131–147. Hu, K.X., Chandra, A., Huang, Y., 1993. Multiple void-crack interaction. International Journal of Solids and Structures 30, 1473–1489. Kadowaki, H., Liu, W.K., 2004. Bridging multi-scale method for localization problems. Computer Methods in Applied Mechanics and Engineering 30–32, 3267–3302. Kadowaki, H., Liu, W.K., 2005. A multiscale approach for the micropolar continuum model. Computer Modeling in Engineering & Sciences 7 (3), 269–282. Kanninen, M.F., Popelar, C.H., 1985. In: Advanced Fracture Mechanics. John Wiley and Sons, Oxford. Konke, C., 1997. Coupling of micro- and macro-damage models for the simulation of damage evolution in ductile materials. Computers & Structures 64 (1–4), 643–653. Lee, B.J., Mear, M.E., 1992. Axisymmetric deformations of power-law solids containing a dilute concentration of aligned spheroidal voids. Journal of the Mechanics and Physics of Solids 40, 1805–1836. Liu, W.K., Chen, Y.J., Belytschko, T., Luo, Y.J., 1996. Three reliability methods for fatigue crack growth. Engineering Fracture Mechanics 53 (5), 733–775. Liu, W.K., McVeigh, C., 2008. Predictive multiscale theory for design of heterogeneous materials. Computational Mechanics 42, 147–170. Liu, W.K., Siad, L., Tian, R., Lee, S., Lee, D., Yin, X., Chen, W., Chan, S., Olson, G.B., Lindgen, L.E., Horstemeyer, M.F., Chang, Y.S., Choi, J.B., Kim, Y.J., 2009. Complexity science of multiscale materials via stochastic computations. International Journal for Numerical Methods in Engineering 80 (6–7), 932–978. Liu, W.K., Jun, S., Zhang, Y.F., 1995. Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids 20, 1081–1106. Li, S., Liu, W.K., 2004. In: Meshfree Particle Methods.. Springer. ¨ Angewandte Lode, W., 1925. The influence of the intermediate principal stress on yielding and failure of iron, copper and nickel. Zeitschrift fur Mathematik und Mechanik 5, 142. McClintock, F.A., Kaplan, S.M., Berg, C.A., 1966. Ductile fracture by hole growth in shear bands. International Journal of Fracture 2 (4), 614–627. McClintock, F.A., 1968. Local criteria for ductile fracture. International Journal of Fracture 4, 103–130. McVeigh, C, Liu, W.K., 2010. Multiresolution continuum modeling of micro-void assisted dynamic adiabatic shear band propagation. Journal of the Mechanics and Physics of Solid 58 (2), 187–205. McVeigh, C., Vernerey, F., Liu, W.K., Brinson, C., 2006. Multiresolution analysis for material design. Computer Methods in Applied Mechanics and Engineering 195, 5053–5076. McVeigh, C., Vernerey, F.J., Liu, W.K., Moran, B., Olson, G.B., 2007. An Interactive microvoid shear localization mechanism in high strength steels. Journal of the Mechanics and Physics of Solids 55 (2), 224–225. McVeigh, C., 2007. Ph.D. thesis, Northwestern University. McVeigh, C., Liu, W.K., 2008. Linking microstructure and properties through a predictive multiresolution continuum. Computer Methods in Applied Mechanics and Engineering 197, 3268–3290. McVeigh, C., Liu, W.K., 2009. Multiresolution modeling of ductile reinforced brittle composites. Journal of the Mechanics and Physics of Solids 57, 244–267. Needleman, A., 1987. A continuum model for void nucleation by inclusion debonding. Journal of Applied Mechanics 54, 525–531. Needleman, A., Rice, J.R., 1978. Limits to ductility set by plastic flow localization. In: Koistinen, D.P., Wang, N.-M. (Eds.), Mechanics of Sheet Metal Forming. Plenum Publishing, New York, pp. 237–267. Needleman, A., Tvergaard, V., Hutchinson, J.W., 1992. Void growth in plastic solids. In: Argon, A.S. (Ed.), Topics in Fracture and Fatigue. Springer-Verlag, New York, pp. 145–178. Nemat-Nasser, S., Hori, M., 1993. In: Micromechanics: Overall Properties of Heterogeneous Materials. Elsevier, Amsterdam. Nishida, A., Matsubara, H., Tian, R., Hazama, O., Suzuki, Y., Araya, F., Nakajima, N., Tani, M., Kondo, M., 2007. Numerical simulation system ‘‘threedimensional virtual plant vibration simulator’’ for nuclear plants by using assembled structural analysis. Journal of Nuclear Science and Technology 6 (3), 376–382. Olson, G.B., 2000. Designing a new material world. Science 288 (5468), 993–998. Olson, G.B., 1997. Computational design of hierarchically structured materials. Science 277 (5330), 1237–1242. Pardoen, T., Hutchinson, J.W., 2000. An extended model for void growth and coalescence. Journal of the Mechanics and Physics of Solids 48, 2467–2512. Rice, J.R., Tracey, D.M., 1969. On the ductile enlargement of voids in triaxial stress fields. Journal of the Mechanics and Physics of Solids 17, 201–217. Rice, J.R., Johnson, M.L., 1970. The role of large crack tip geometry changes in plane strain fracture. In: Kanninen, M.F. (Ed.), Inelastic Behaviour of Solids. McGraw-Hill, New York, pp. 641–672. Tian, R., Moran, B., Liu, W.K., Olson, G.B., 2008. Multiscale fracture simulator. Dynamic Microstructure Design Consortium (ONR Contract: N00014-05-C0241) Base Final Report. Tian, R., Liu, W.K., Chan, S., Olson, G.B., Tang, S., Wang, J.S., Jou, H.J., Gong, J.D., Moran, B., 2009. Linking Microstructures to Fracture Toughness—predictive 3D process zone simulations. The D 3-D Annual PI Review, Evanston, IL, March 23–25.

1700

R. Tian et al. / J. Mech. Phys. Solids 58 (2010) 1681–1700

Tian, R., Yagawa, G., 2007. Non-matching mesh gluing by meshless interpolation—an alternative to Lagrange multipliers. International Journal for Numerical Methods in Engineering 71, 473–503. Tian, R., Yagawa, G., 2005. Generalized node and high-performance elements. International Journal for Numerical Methods in Engineering 64, 2039–2071. Tian, R., Yagawa, G., Terasaka, H., 2006a. Linear dependence problems of partition of unity based generalized FEMs. Computer Methods in Applied Mechanics and Engineering 195, 4768–4782. Tian, R., 2006. A PU-based 4-node quadratic tetrahedron and linear dependence elimination in three dimensions. International Journal of Computational Methods 3, 545–562. Tian, R., Matsubara, H., Yagawa, G., 2006b. Advanced 4-node tetrahedrons. International Journal for Numerical Methods in Engineering 68, 1209–1231. Tian, R., Yagawa, G., 2006. Allman’s triangle, rotational dof and partition of unity. International Journal for Numerical Methods in Engineering 69, 837–858. To, A.C., Liu, W.K., Kopacz, A., 2008a. A finite temperature continuum theory based on interatomic potential in crystalline solids. Computational Mechanics 42 (4), 531–541. To, A.C., Liu, W.K., Olson, G.B., Belytschko, T., Chen, W., Shephard, M.S., Chung, Y.W., Ghanem, R., Voorhees, P.W., Seidman, D.N., Wolverton, C., Chen, J.S., Moran, B., Freeman, A.J., Tian, R., Luo, X., Lautenschlager, E., Challoner, A.D., 2008b. Materials integrity in microsystems: a framework for a petascale predictive-science-based multiscale modeling and simulation system. Computational Mechanics 42 (4), 485–510. Thomason, P.F., 1990. In: Ductile Fracture of Metals. Pergamon Press, Oxford. Thomason, P.F., 1993. Ductile fracture by the growth and coalescence of microvoids of non-uniform size and spacing. Acta Metallurgica et Materialia 41, 2127–2134. Tvergaard, V., 1981. Influence of voids on shear band instabilities under plane strain conditions. International Journal of Fracture Mechanics 17, 389–407. Tracey, D.M., 1971. Strain hardening and interaction effects on the growth of voids in ductile fracture. Engineering Fracture Mechanics 3, 301–315. Traiviratana, S., Bringa, E.M., Benson, D.J., Meyers, M.A., 2008. Void growth in metals: atomistic calculations. Acta Materialia 56 (15), 3874–3886. Tvergaard, V., Needleman, A., 1984. Analysis of the cup-cone fracture in a round tensile bar. Acta Metallurgica 32, 157–169. Tvergaard, V., Needleman, A., 1987. An analysis of ductile rupture modes at a crack tip. Journal of the Mechanics and Physics of Solids 35 (2), 151–183. Tvergaard, V., Hutchinson, J.W., 1992. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. Journal of the Mechanics and Physics of Solids 40 (6), 1377–1397. Tvergaard, V., Hutchinson, J.W., 2002. Two mechanisms of ductile fracture: void by void growth versus multiple void interaction. International Journal of Solids and Structures 39, 3581–3597. Uchic, M., Groeber, M., Dimiduk, D., Simmons, J., 2006. 3D microstructural characterization of nickel superalloys via serial-sectioning using a dual beam FIB-SEM. Scripta Materialia 55 (1), 23–28. Vernerey, F.J., McVeigh, C., Liu, W.K., Moran, B., Tewari, D., Parks, D.M., Olson, G.B., 2006. The 3-D computational modeling of shear-dominated ductile failure in steel. Journal of The Minerals, Metals and Materials Society, 45–51. Vernerey F.J., 2006. Ph.D. thesis, Northwestern University. Vernerey, F.J., Liu, W.K., Moran, B., 2007. Multi-scale micromorphic theory for hierarchical materials. Journal of the Mechanics and Physics of Solid 55 (12), 2603–2651. Vernerey, F.J., Liu, W.K., Moran, B., Olson, G.B., 2008. A micromorphic model for the multiple scale failure of heterogeneous materials. Journal of the Mechanics and Physics of solids 56 (4), 1320–1347. Wang, J.S., Olson, G.B., 2009. Prediction of ductile fracture toughness. In: Proceedings of the 12th International Fracture Conference, Ottawa, Canada, July 12–17. Weck, A., Wilkinson, D.S., 2008. Experimental investigation of void coalescence in metallic sheets containing laser drilled holes. Acta Materialia 56, 1774–1784. Xia, L., Shih, C.F., 1995. Ductile crack growth-I. A numerical study using computational cells with microstructurally-based length scales. Journal of the Mechanics and Physics Solids 43 (2), 233–259. Yin, X., Chen, W., To, A., McVeigh, C., Liu, W.K., 2008. Statistical volume element method for predicting microstructure-constitutive property relations. Computer Methods in Applied Mechanics and Engineering 197, 3516–3529.