当前位置：文档库 › Calibration procedures for a computational model of ductile fracture

Calibration procedures for a computational model of ductile fracture

$Calibration procedures for a computational model of ductile fracture$

Calibration procedures for a computational model of ductile fracture Z.Xue a,1,M.G.Pontin b,2,F.W.Zok b ,J.W.Hutchinson a,*

School of Engineering and Applied Sciences,Harvard University,Cambridge,MA,United States b Materials Department,University of California,Santa Barbara,CA,United States

a r t i c l e i n f o Article history:Received 18August 2009Received in revised form 22October 2009Accepted 29October 2009Available online 1November 2009Keywords:Ductile fracture Computational fracture Shear fracture Damage parameters

a b s t r a c t

A recent extension of the Gurson constitutive model of damage and failure of ductile struc-

tural alloys accounts for localization and crack formation under shearing as well as tension.

When properly calibrated against a basic set of experiments,this model has the potential to

predict the emergence and propagation of cracks over a wide range of stress states.This

paper addresses procedures for calibrating the damage parameters of the extended consti-

tutive model.The procedures are demonstrated for DH36steel using data from three tests:

(i)tension of a round bar,(ii)mode I cracking in a compact tension specimen,and (iii)shear

localization and mode II cracking in a shear-off specimen.The computational model is then

used to study the emergence of the cup-cone fracture mode in the neck of a round tensile

bar.Ductility of a notched round bar provides additional validation.

Progress in computational fracture mechanics has paralleled advances in constitutive models that incorporate damage mechanisms.For many ductile structural alloys the mechanism governing failure is void nucleation,growth and coalescence.The grand challenge for these alloys is the development of a computational capability for predicting localization,crack for-mation and crack propagation under all states of stress.Capturing both tensile (mode I)and shear (mode II)fractures has been particularly challenging.When properly calibrated for a speci?c structural alloy,the Gurson model [1]and some of its close relatives,such as the Rousselier model [2],have shown considerable promise for characterizing mode I crack growth

[3–8].In addition,the models have been used to simulate transitions from mode I crack growth to mixed mode shear crack-ing in the cup-cone fracture process of round tensile bars [9,10]and in three-dimensional through-cracks in thin plates [11].Such transition problems are generally more challenging because the constitutive models have not been developed to explic-itly address damage under shear dominated conditions.

A recent extension of the Gurson model [12]speci?cally incorporates damage in shear,adding the ?exibility to address shear ruptures as well as tension dominated failures.This extension will be employed here in conjunction with a suite of three tests (round bar tension,mode I compact tension,and mode II shear-off)to calibrate the constitutive parameters for the structural steel,DH36.For veri?cation,the calibrated model is then used to study the failure details of several other problems.

To put the overall objectives of this work into some perspective,it is noted that three parameters are required to calibrate the extended Gurson model:the initial void volume fraction,f 0,a shear damage coef?cient,k x (de?ned below)and the ?nite element size,D .To accurately characterize localization and fracture,D must be on the order of the spacing between the voids 0013-7944/$-see front matter ó2009Elsevier Ltd.All rights reserved.

doi:10.1016/j.engfracmech.2009.10.007

*Corresponding author.

E-mail address:hutchinson@https://www.wendangku.net/doc/685396535.html, (J.W.Hutchinson).

1Present address:Schlumberger Reservoir Completions,Rosharon,TX,United States.

2Present address:Ceradyne,Costa Mesa,CA,United States.

Engineering Fracture Mechanics 77(2010)492–509

Contents lists available at ScienceDirect

Engineering Fracture Mechanics

j o u r n a l h o m e p a g e :w w w.elsevier.c om /loc ate/engfracmech

that dominate the fracture process,typically from tens to hundreds-of microns.With mesh requirements this ?ne,it is only possible to predict the onset and propagation of cracks in relatively small components or in larger structures where the loca-tion of the failure can be anticipated in advance.In contrast,it would not be feasible to employ a fracture model of this type to analyze fractures in large structures where the failure locations cannot be anticipated.Under such circumstances,because the ?nite element size for a large structure is necessarily orders of magnitude greater than void spacing and often larger than plate thickness,coarser criteria based on a critical effective plastic strain or a through-thickness cohesive zone must be em-ployed.These criteria must also be calibrated for each material,but against tests that make no attempt to resolve the ?ne scale fracture processes relevant for the present class of models.The two classes of fracture models complement each other.In principle,computations based on a ?ne scale model could be used to calibrate a coarse scale model.

2.The extended Gurson model

The Gurson model is an isotropic formulation that employs the mean stress,r m =r kk /3,and the effective stress,r e ???????3J 2p ??????????????????3s ij s ij =2p ,where s ij ?r ij à13r kk d ij

is the stress deviator.The extended model [12]employs,in addition,the third stress invariant

J 3?det es T?13s ij s ik s jk ?er I àr m Ter II àr m Ter III àr m Te1Twhere the expression on the right is couched in terms of principal stresses,assumed to be ordered as

r I P r II P r III .The

non-dimensional metric x er T?1à27J 3r 3e 2e2T

lies in the range,06x 61,with x ?0for all axisymmetric stress states,

r I P r II ?r III or r I ?r II P r III ;

e3T

and x ?1for all states comprised of a pure shear stress plus a hydrostatic contribution,r I ?s tr m ;r II ?r m ;r III ?às tr m es >0Te4TThe original Gurson model was formulated and calibrated based on the mechanics of void growth under axisymmetric stress states.The extension [12]does not alter the model for these states.The extension modi?es the predictions for states with non-zero x er T.In particular,a contribution to damage growth under pure shear stress states is accounted for in the extension whereas the original Gurson model predicts no change in damage for states having r m ?0.

Nomenclature

A 0;A

cross-sectional area of neck:initial,current D

characteristic element size D P ij plastic strain rate E

Young’s modulus f 0;f ;f c ;f f

void volume fraction:initial,current,onset of coalescence,failure H

plate thickness J 3

stress invariant k x

shear damage coef?cient N

strain hardening exponent q 1;q 2;q 3

?tting parameters in Gurson model R

punch radius s ij

stress deviator d punch displacement

e f

ductility—true strain in neck at failure e P M ;r M

intrinsic true plastic strain and stress in tension (damage-free)e peak T ;r peak T true strain and stress at maximum nominal stress r ij ;r e ;r m true stress,effective stress,mean stress r I P r II P r III true principal stresses

x measure of shearing relative to axisymmetric stressing

Z.Xue et al./Engineering Fracture Mechanics 77(2010)492–509493

The yield surface of the extended Gurson model is the same as the original.Including the ?tting parameters,q 1,q 2and q 3,introduced by Tvergaard [13],it is given in terms of the effective and mean stress measures by

F er e ;r m ;f T?r e r M 2t2q 1f cos h 3q 2r m r M àe1tq 3f 2Te5T

The current state is characterized by f ,the ‘‘apparent”void volume fraction,and r M ,the current effective stress governing ?ow of the damage-free matrix material.All quantities not labeled with the subscript M represent overall quantities asso-ciated with the bulk material.Normality implies that the plastic strain rate,D P

ij ,is given by

D P

ij ?

1h P ij P kl _r kl e6Twhere P ij ?@F r ij ?3s ij r M tfq 1q 2r M sin h 3q 2r m r M d ij e7T

In ?nite strain formulations,_r

ij is identi?ed with the Jaumann rate of stress.The hardening modulus,h ,is identi?ed in the Appendix A .If r m ?0,P kk ?0and the rate of plastic volume change vanishes,i.e.,D P

kk ?0;this feature persists in the exten-

sion.In the absence of nucleation,the extension of the Gurson model posits

_f ?e1àf TD p kk tk x f x er Ts ij D p ij r e e8T

The ?rst contribution is that incorporated in the original model while the second is the crux of the extension.As previ-ously noted,the modi?cation leaves the constitutive relation unaltered for axisymmetric stress states.In a state of pure

shear,however,(8)gives _f ?k x f _c

P =???3p ,where _c P is the plastic shear strain rate and k x is the shear damage coef?cient,the sole new parameter in the extended model.The inclusion of the second term in (8)rests on the notion that the volume of voids undergoing shear may not increase,but void deformation and reorientation contribute to softening and constitute an effective increase in damage [14–16].In addition,the second term can model damage generated by the nucleation in shear of tiny secondary voids in void sheets linking larger voids.Thus,in the extension,f is no longer directly tied to the plastic volume change.Instead,it must be regarded either as an effective void volume fraction or simply as a damage param-eter,as it is for example when the Gurson model is applied to materials with distinctly non-spherical voids.Further discus-sion and illustrations of the extension are given in [12],where the emphasis is on its role in shear localization.The remaining equations specifying the entire description of the model are listed in the Appendix A .Included is the speci?cation of the widely used technique [13]that accelerates damage from f ?f c to f ?f f ,at which point the material element is deleted.De-tails of the numerical algorithm used to implement the constitutive model in the ?nite element code ABAQUS Explicit [16]are also presented in the Appendix A .

3.Outline of the calibration protocol

The elastic–plastic inputs into the extended Gurson Model are the Young’s modulus,E ,the Poisson’s ratio,m ,and the intrinsic stress–strain response of the damage-free material (f 0?0).The two damage-related input parameters are the

initial Fig.1.Optical micrograph of polished and etched cross-section through DH36steel plate,showing a microstructure of ferrite (light)and pearlite (dark).

494Z.Xue et al./Engineering Fracture Mechanics 77(2010)492–509

effective void volume fraction,f0,and the shear damage coef?cient,k x.Additionally,because the constitutive model contains no material length scale,the size of the?nite element mesh,D,is calibrated through crack growth predictions,employing well-established procedures[4,7].

This paper addresses the general task of calibrating the three fracture-related parameters:f0,k x and D.The procedures are demonstrated through experiments and analyses of DH36steel(Fig.1):a high strength alloy commonly used in ship con-struction.Following extensive prior work on calibration procedures for the standard Gurson model(e.g.,[4,7]),the present study employs data from a mode I fracture test and a round bar tensile test to identify intrinsic uniaxial stress–strain behav-ior,f0and D.Additionally,a shear-off test is added to the suite of tests to determine the shear damage coef?cient,k x.The paper is organized following closely the steps in the calibration protocol:

Section4:Determination of the intrinsic stress–strain response of the undamaged material from round bar tensile tests and establishing that f0,k x and D have little in?uence on the plastic response until neck development is quite advanced.

Section5:Determination of f0and D from compact tension mode I fracture tests and establishing that k x has little in?u-ence on crack growth prediction when the crack is planar.

Section6:Determination of k x using data from shear-off tests and the previously determined f0and D.

Section7:Discussion of the applicability of the calibrated constitutive model to the cup-cone failure mode as one illus-tration and the ductility of notched round bars as another.Possible variations in the identi?cation protocol for other materials are also discussed.

The three calibration tests were conducted under quasi-static loading,while all simulations were carried out using the dynamic code ABAQUS Explicit.In order to minimize inertial effects and ef?ciently simulate the quasi-static tests in the ex-plicit code,a preliminary series of calculations with different?xed applied loading rates was performed for each test con-?guration.At some loading rate,as the rates decrease,the simulations converge to a quasi-static limit.That loading rate

was then employed in all subsequent calculations.Material strain rate dependence is ignored in the present

computations.

Fig.2.Tensile specimen geometry and?nite element mesh.

Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509495

4.Intrinsic plastic response of the undamaged material

The plastic response of the undamaged material (f 0?0)was obtained from quasi-static uniaxial tensile tests on round bars coupled with elastic–plastic ?nite element computations.The test geometry and ?nite element mesh are shown in Fig.2.The nominal axial strain e N was measured using a non-contacting laser extensometer over a central 12.7mm length within the gauge section.Prior to necking,the true (logarithmic)strain is given by e T ?ln e1te N Tand the true stress by r T ?r N e1te N T,where r N is the nominal stress (load/initial area).To ascertain the true response in the post-necking regime,computations were performed using an assumed form of the stress–strain relation (detailed below)and matching the pre-dicted nominal stress–strain curves with those obtained experimentally.To accurately capture strain localization,a ?nite strain formulation of elasto-plastic theory was employed in the ?nite element model.Four-node axisymmetric elements 0 0.1 0.2 0.3 0.4 0.5

900

800

700600500400

300

200

100

N=0.20

0.185

0.16

Experimental True strain, εT

0 0.1 0.2 0.3 0.4

700

600

500

400

300

200

100

0N=0.200.1850.16Experimental Nominal strain, εN

εσT peak T peak ,()of the true tensile stress–strain curve beyond the onset of necking and (b)the corresponding element analysis.Error bars represent the full range of experimental measurements from six tests.extensometer over a 12.7mm gauge length near the specimen center.The nominal strain,de?ned as consistently employed in both the experiments and the ?nite element calculations.The 496Z.Xue et al./Engineering Fracture Mechanics 77(2010)492–509

with reduced Gaussian integration (CAX4R in ABAQUS/Explicit [16])were used.The model was based on an axisymmetric mesh comprised of square section elements with size,D =50l m,providing more than 30elements across the gauge radius.The element size was selected to be consistent with the value emerging from the calibration of the mode I fracture data,pre-sented in the next section.Nevertheless,since the selected element size is already very much smaller than the macroscopic specimen dimensions and hence the strains are adequately resolved,further reductions in element size would have essen-tially no effect on the intrinsic (damage-free)stress–strain response.Additional computations were performed to demon-strate that f 0and k x do not affect the identi?cation of the true stress–strain curve even up to strains approaching that for rupture.

The average true stress–strain curve from ?ve tensile tests is plotted in Fig.3a.This curve was subsequently used to char-

acterize the stress–strain response for stresses below that corresponding to the load maximum,denoted r peak

.To extrapolate beyond r peak T ,a true stress–strain curve of the form r T ?r peak T ee T =e peak T

TN was assumed.A preliminary estimate of the strain hardening exponent N was obtained by a least squares ?t of the small strain data.A series of ?nite element computations was then performed to ascertain the full nominal tensile stress–strain curve,using a range of values of N ,guided by the pre-ceding curve ?tting.As shown in Fig.3b,the results for N ?0:185(and f 0?0)accurately replicate the experimental mea-surements up to the onset of rupture (at a nominal strain of e N ?0:32).In summary,the true stress–strain curve used to

characterize the damage-free material (f 0?0)is given by the experimental curve below r peak

T and the power law extrapo-lation at stresses above r peak T .

For e N <0:3,void growth has almost no effect on the tensile behavior of DH36.This result is demonstrated in Fig.4by comparing the experimental data with ?nite element computations based on a hardening exponent N ?0:185and several representative initial void volume fractions (including the Mises limit,wherein f 0?0).Other than f 0,k x and D ,the basic parameters characterizing the constitutive model that are used in all simulations in this paper are:

E ?210MPa ;m ?0:3;N ?0:185;q 1?1:5;q 2?1;q 3?2:25;f c ?0:15and f f ?0:25e9TThe comparisons show that the effects of void growth,manifested in a divergence in the stress–strain response from that of a Mises material,are important only very near the point of ?nal rupture for the DH36tensile specimen.Their effect is to accelerate the softening of the material such that the load drops more rapidly than that predicted for the damage-free mate-rial.Further details of the failure process in the neck,including formation of a cup-cone fracture surface,are presented in Section 7.

5.Determination of f 0and D from compact tension test

Compact tension tests were performed on specimens with the geometry shown in Fig.5a.Crack mouth opening displace-ment was measured using a non-contacting extensometer and a pair of ?ducial tapes mounted on the specimen edge,sep-arated by a distance of 14mm.Optical images of the broad sample surface were periodically recorded.The experimental 0 0.1 0.2 0.3 0.4

700

600

500

400

300

200

100

0f o = 0.0010.0020.003Experimental Nominal strain, εN

f o =0(Mises)k ω = 0

fraction f o on the computed nominal tensile stress–strain response.Over the pertinent range of experimental measurements up to the onset of fracture.

Z.Xue et al./Engineering Fracture Mechanics 77(2010)492–509497

measurements and observations are summarized in Figs.6and 7.Signi?cant nonlinearity due to plasticity is evident in both the load–displacement response and in the optical images at displacements above 0.5mm.Following an initial rising por-tion,the load–displacement curve reaches a maximum,at a displacement of about 3–4mm.This point corresponds to the emergence of a crack on the external surface of the sample (Fig.7d–f ).Further growth both at the surface and in the interior occurs under decreasing load.

The corresponding ?nite element model is shown in Fig.5b.In the present analysis,deformations are restricted to be symmetric with respect to the mid-plane such that a symmetry boundary condition is applied to the mid-plane.Conse-quently,the region meshed is only one half of the full specimen.Eight-node brick elements with reduced Gaussian integra-tion (C3D8R in ABAQUS/Explicit [16])were used.Iterations on element size and meshing details were made prior to arriving at the mesh used to carry out the ?nal analysis.The smallest elements at the mid-plane in the vicinity of the crack tip have dimensions 30?30?50l m with 50l m in the through-thickness direction.Near the surface of the specimen and near the tip the element dimensions are 30?30?80l m.Approximately 100elements extend from the mid-plane to the surface in the vicinity of the crack tip.The 30l m in-plane mesh at the tip allows accurate resolution of the initial tip notch.Further away from the notch tip in the region of crack propagation,the in-plane dimensions of the mesh are approximately 50?50l m.Relatively small differences in results were found from a series of computations with different meshes with ele-ment dimensions in the range from 30l m to 50l m.The mesh in Fig.5b is regarded as having a nominal (characteristic)size D =50l m.In order to improve computational ef?ciency,only the material in the region of crack propagation,which

starts

Fig.5.(a)Compact tension test geometry employed in the experimental study and (b)corresponding ?nite element model.Specimen thickness is 12.5mm.Crack mouth opening displacements were measured using a non-contacting extensometer and a pair of ?ducial tapes mounted on the specimen edge,separated by a distance of 14mm.The same de?nition was used in the subsequent ?nite element calculations.

498Z.Xue et al./Engineering Fracture Mechanics 77(2010)492–509

Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509499

from the notch tip to the left edge of the specimen and has width of7mm,was modeled using the extended Gurson model. Outside this region,the specimen was modeled using von Mises plasticity(i.e.,f0?0and k x?0).

Load–displacement predictions for four values of f0(including f0?0)and k x?2are compared with the experimental results in Fig.6.Over the range plotted,the load of the damage-free specimen increases monotonically with displacement

because there is no damage-induced softening or crack growth.In contrast,the prediction for f 0?0:001follows the exper-imental curve closely for displacements as large as 5mm.Furthermore,it predicts that cracking initiates at the center of the notch front,at a displacement of about 1mm.Thereafter,the crack grows deeper into the specimen and spreads laterally from the center (Fig.7).Upon reaching the free surface,at a displacement of 3.6mm,the load reaches a maximum and a load fall-off ensues.These results agree well with the experimental measurements.The predictions for the two larger values of f 0clearly over-predict the effect of damage and cracking at displacements below 5mm.They are particularly de?cient in predicting the displacement at the load maximum.

At displacements above 5mm,the experimental data fall below the numerical predictions for all three values of f 0.This discrepancy arises for two reasons.The symmetry imposed in the simulation precludes the transition to slant fractures that usually develop as the crack advances and the crack in the test is likely to have departed from the imposed symmetry.In addition,element deletion was used to mimic the crack propagation such that the element is deleted when f ?f f .As the crack advances,it encounters larger elements in the mesh and these dissipate more energy prior to failure than the cali-brated elements with D =50l m.It is indeed observed from Fig.8for the case of the crack month opening displacement reaching 8mm that some of the deleted elements are much larger than D =50l m.It remains for the future to verify that predictions based on the present choices of f 0and D can replicate the present experimental results for larger displacements using a computational model with no symmetry restrictions,as well as a uniform mesh with the same calibrated element size throughout the region of crack propagation.Unfortunately,this would result in a signi?cant increase in computational size that would not be feasible for the calibration procedure.3

Although the results in Fig.6b were computed with k x ?2,the shear damage coef?cient has essentially no effect on these predictions.To illustrate this,results for f 0?0:001computed with k x =2,2.5and 3are plotted in Fig.6a.The response under-goes only very slight softening with increasing k x but remains well within the range of the experimental data.The weak dependence on k x is consistent with the fact that mode I cracking occurs over the range of load–displacement data used for the

?tting.

Fig.7.Images of broad face of compact tension specimen with increasing crack mouth opening displacements.Arrows in the right column indicate the emerging near-surface crack.

More than ten days were required for each calculation based on the current mesh using a personal computer with memory requirements up to 1GB.The trade-off between ef?ciency and accuracy suggests that the present calibration strategy is a reasonable compromise.500Z.Xue et al./Engineering Fracture Mechanics 77(2010)492–509

In summary,based on the agreement between prediction and experiment for displacements below5mm,the choices

f0?0:001with D%50l m are made for DH36.

6.Determination of k x from a shear-off test

The?xture in Fig.9was designed to create a controlled test in which shear localization gives way to mode II fracture[17].

The corresponding load–displacement curve is used to infer the shear damage coef?cient,k x.In the test,a plate specimen

(3mm thick)is clamped between two thick steel platens,each with a through-hole of diameter19.2mm.Cylindrical steel

plungers,19.05mm in diameter,are inserted into each of the two holes,leaving a narrow(0.075mm)radial gap between the

plunger surface and the hole.An additional pair of plungers with slightly reduced diameter(to accommodate Te?on bear-

ings)is then inserted into the holes.The four plungers and the test specimen are then clamped together with a single bolt

passing through open holes in each of three of the plungers and the test specimen and a threaded hole in the last plunger,as

shown in Fig.9.With one side of the assembly placed on a stiff supporting base,the plunger on the opposite side is load

axially in compression.The movement of the plungers induces shear deformation within a narrow cylindrical ring in the

specimen.Failure starts as shear localizations near the upper and lower surfaces of the plate which subsequently develop

into mode II cracks as the deformation progresses into the plate.

The experimental measurements are summarized in Fig.10.The coordinate axes are the nominal applied shear stress, s P=e2p RHT(R being the plunger radius and H the plate thickness)and the normalized displacement,d=H.The resulting curves exhibit features reminiscent of those obtained in tension tests.That is,the initial linear region gives way to plasticity

at a shear stress of r O=2%240MPa(r O being the tensile yield stress,obtained from Fig.3).Following a period of strain hard-ening,the load reaches a peak,at a displacement of d=H%0.3–0.4,and subsequently diminishes with increasing displace-

ment.Scanning electron micrographs of a cross-section through a test specimen that had been interrupted following

loading to a displacement d=H%0.5are presented in Fig.11.They reveal a diffuse damage zone within the region of intense

shear as well as well-de?ned shear cracks emanating from the specimen surface in the vicinity of the plunger periphery.

A detail of the?nite element mesh is depicted in the inset of Fig.9a.Based on the prior calibrations,computations of

shear-off employ an initial void fraction f0?0:001and element size D=50l m in the region of shear localization

and

Fig.8.Evolution of plastic strain and crack growth from?nite element calculations of the compact tension test.

Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509501

cracking.As in the compact tension simulations,computational ef?ciency was enhanced by only employing the extended Gurson model and the smallest elements in the region of shear localization.Outside this region,the plate was modeled using Mises plasticity and represented by a coarser mesh.Boundary conditions were applied such that the bottom of the lower clamping plate,as well as a small section of the upper clamping plate representing the constraining effect of the clamping bolts,were restricted from all rotation and displacement.Contact between the plate and the clamps and punch was enforced with no tangential sliding.Separation was permitted when the normal traction became tensile.The punch and clamps were modeled as being elastic.Four-node axisymmetric elements with reduced Gaussian integration and hourglass control (CAX4R in ABAQUS Explicit[16])were used for all components.

The results of the simulations are plotted and compared with experimental measurements in Fig.10.Simulations for the Mises materialef0?0;k x?0Tcorrelate well with the experimental data for small displacements(d=H<0.15).However, they over-predict the stresses at larger displacements and do not reveal a load maximum.The results for the standard Gur-son modelek x?0Tyield essentially identical results up to d=H%0.35,with only small amounts of softening at larger dis-placements.Among the other simulations,the best?t of the experimental data for displacements d=H<0:5is that with k x?2:5.Fig.10gives a clear trend of the sensitivity of the predictions to the shear damage parameter:the erosion of the shear-off load is signi?cantly underestimated if k x?1and signi?cantly overestimated if k x?4.

Some details of the progression of the shear-off process–at maximum loaded=H?0:32Tand at a point just before entire ligament fracture(d=H?0:49)–are shown in Fig.12.At maximum load,shear localization and fracture has occurred at the top and bottom surfaces of the plate and some damage has occurred across the entire plate thickness.However,the level of damage in the central region is no larger than about f?0:005,well below that at which the shear stress reaches a maximum

(f%0:03)[12].The inference is that for DH36the present test leads to mode II crack propagation emanating from the

plate

Fig.9.(a)Shear-off test assembly and detail of?nite element model in the region of intense plastic deformation.(b)Macrophotograph of specimen interrupted during shear-off test and sectioned by electrodischarge machining(EDM).

502Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509

Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509503 surfaces rather than global(net-section)rupture.In the second case,for d=H?0:49,the ring cracks have extended well into the plate interior and the damage parameter in the center has almost reached the failure level,f%f F?0:25.The ligament undergoes complete fracture in the next increment and the load drops abruptly to zero.The?nal stages of the failure process are not accurately captured because the simulation does not account for friction between contacting crack surfaces.Element deletion also plays a role.In the narrow region of shear-off,the plastic strain and damage is almost uniform before the shear-off fracture occurs.Multiple elements are predicted to fail almost simultaneously predicting of a loss of load carrying

capacity that is almost certainly too rapid.Nevertheless,the main features of the initiation of damage growth in shear and its progression to a well-developed mode II crack in DH36appear to be well represented by the extended Gurson Model with k x ?2:5when f 0?0:001and D =50l m.

7.Applications of the computational model

7.1.Cup-cone fracture of a round tensile bar

In the discussion of Fig.4it was noted that damage plays almost no role in the tensile behavior of a round bar of DH36well beyond the onset of necking.Not until the load has fallen to below 60%of the maximum load does damage have notice-able effect on the overall load–elongation behavior.In this section,the computational model with the calibrated parameter values for DH36(f 0?0:001,k x ?2:5)is used to analyze the development of damage within the neck and the trajectory of the ensuing macroscopic crack.This provides a further test of the predictive capability of the extended Gurson model.The key feature of interest is the transition from normal to shear fracture that gives rise to a cup-cone appearance.

Details of the fracture surfaces of the round tensile bars are depicted in Fig.13.In addition to the classical cup-cone shape observed when the specimen is viewed at low magni?cations,three other features are evident.(i)The central (cup)region comprises equi-axed ductile dimples associated with the growth and coalescence of voids in mode I.Although a broad dis-tribution of dimple size is apparent,the average value appears to be of the order of 10l m.This dimension correlates with the spacing between pearlite colonies (Fig.1)and suggests that the pearlite serves as the principal void nucleating constit-uent.(ii)The cone region comprises highly-elongated ‘‘smeared”dimples,consistent with void coalescence by shear local-ization.The latter dimples are comparable in size to those in the cup region,suggesting that the same population of void nucleating sites is activated in both failure modes.(iii)In some regions,the cone consists of more than one shear fracture plane.For instance,on the right side of the surface in Fig.13a,the transition from cup to cone ?rst occurs by a downward deviation of the mode I crack as it grows radially from the specimen center.It subsequently deviates from this path and adopts an upward shear path,thereby yielding two distinct shear lips.Furthermore,a closer examination of the cup region near the ?rst transition further suggests that analogous processes occur at smaller length scales.That is,the mode I crack ?rst attempts to deviate onto an upward path (over a distance of about 100l m)before diving down to form the ?rst dom-inant shear lip.These deviations in crack path are likely due to the severity of the neck and the corresponding reduction in stress at even short axial distances from the neck

plane.

Fig.12.Evolution of plastic strain and damage from ?nite element calculations of a shear-off test.Inset (right)shows a shear crack developed in the vicinity of the edge of the contacting punch.

504Z.Xue et al./Engineering Fracture Mechanics 77(2010)492–509

Tvergaard and Needleman [9]carried out the ?rst detailed computational study of the failure mode in the neck of a round tensile bar based on the unmodi?ed Gurson model.Their work demonstrated that a transition from the planar mode I crack nucleated in the center of the neck to a conical mixed mode shear crack can occur for this constitutive model if a

suf?ciently Fig.13.Fracture surface of DH36tensile bar showing:(a)cup-cone failure mode;(b,c)equi-axed dimples caused by void growth and coalescence in the central region;and (c,d)elongated dimples formed by void coalescence during shear lip

formation.

Fig.14.Effects of initial element size and k x on the crack trajectory in an initially unnotched round tensile bar.

Z.Xue et al./Engineering Fracture Mechanics 77(2010)492–509505

506Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509

?ne mesh is used and if a relatively large damage level is invoked.These authors took f0?0and assumed that a4%volume fraction of voids would be nucleated under increasing strain.Thus,the total void volume fraction nucleated in their simu-lations far exceeds the void fraction considered to be representative for materials such as DH36.The present calculations suggest that,for realistic void volume fractions(of order10à3),the transition to conical shear cracking does not occur when the unmodi?ed Gurson model is employed.This?nding is borne out by an extensive study of fracture modes in round tensile bars and in-plane strain specimens by Besson et al.[10]using several damage-based constitutive models.More recently,Le-blond[18]has pursued these issues further by considering the extended Gurson Model with?ndings similar to those re-ported below.

Even when shear damage is included(k x?2:5),simulations with a square mesh in the neck(50?50l m)do not predict a transition to a conical crack.Reducing the initial element height such that the element aspect ratio at the onset of fracture is approximately unity at the onset of fracture accommodates a mixed mode conical crack propagating at roughly45°to the axis of the specimen.Although this modi?cation does not lead to a transition when the initial element width is set at

Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509507 D=50l m,a well-de?ned cup-and-cone fracture mode is predicted for slightly smaller element widths(Fig.14).The fracture patterns in Fig.14were computed using the mesh just described for deformations restricted to be axisymmetric but with no symmetry imposed with respect to the plane through the center of the neck.The mesh(40?6l m)in Fig.14c gives rise to a near-planar crack in the center of the specimen followed by the transition to a conical crack after a hesitating start in the opposite direction,broadly consistent with the experimental observations.

7.2.Ductility of straight and notched round bars

The standard de?nition of the ductility of a metal alloy is the logarithmic strain at failure of a round tensile bar as deter-mined by e f?lneA0=ATwith A0=A being the ratio of initial to?nal cross-sectional areas at the neck.The ductility predicted for the round bar of DH35with f0?0:001,k x?2:5and the40?6l m mesh is e f?1:38.This value is in close agreement with that measured experimentally:e f?1:35?0:04(from?ve specimens).The ductility prediction is not nearly as sensitive to

meshing as the prediction of the transition to the slanted fracture path.For example,the ductility predictions for the other meshes in Fig.13are e f?1:41for(a),e f?1:44for(b)and e f?1:36for(d).The fact that ductility predictions are less sen-sitive to meshing details than crack path transition is consistent with the fact that the overall load–elongation behavior is also relatively insensitive to meshing details.This can be seen in Fig.15where nominal stress–strain curves are presented corresponding to some of the same meshes used in the mode transition study in Fig.14.The cross-sectional area of the neck becomes nearly‘‘frozen”as soon as a normal localization band forms in the center of the neck much before the mode tran-sition.Thus,an accurate ductility prediction relies primarily on the ability of the constitutive model to capture the onset of a normal localization since the onset itself is not very sensitive to mesh size,assuming the mesh is adequate to accurately re-solve the stresses and strains in the neck.

As a?nal validation of the calibrated computational model,the ductility of a notched round bar of DH36has been com-puted.The specimen geometry and the mesh in the critical region are shown in Fig.16.The predicted ductility is e f?0:98. By comparison,the experimentally measured values from three test specimens fall in the range e f?0:91—0:93.Thus the model correctly predicts the reduction in ductility due to the increased stress triaxiality arising from the notch geometry.

8.Concluding remarks

This paper has demonstrated that,when properly calibrated,the extended Gurson model has considerable promise as a computational tool for predicting the formation of cracks and their subsequent propagation in ductile structural alloys under a wide range of stress states.By incorporating a parameter to characterize damage in shear,the extended model widens the scope of applications to failure modes with a heavy component of shear.The calibration protocol employed here uses three types of tests:(i)uniaxial tension of a round bar,to infer the intrinsic stress–strain behavior of the undamaged material;(ii) mode I cracking in a compact tension specimen,to calibrate the initial void volume fraction and the element size;and(iii) mode II cracking in a newly-designed shear-off test,to determine the shear damage coef?cient.For the alloy in the present study,DH36,it was established that these three calibration steps can be conducted independently,assuming that the sequential order listed above is followed.The calibration process might turn out to be more complicated for other materials, e.g.,the shear damage coef?cient might in?uence the calibration of the other two parameters in step(ii).It is worth noting that a variation on the procedure employed here in step(ii)would be to choose f0and D to?t resistance curve data in the

eD aT,extracted from a side-grooved compact tension specimen designed to sustain a form of the J-integral vs.crack growth,J

?C r Y D where C lies in the range from2to5depending on N straight crack front.The work of Xia and Shih[4]reveals that J

and f0.For DH36with N?0:185and f0?0:001,C?5such that the formula gives J IC?120kJ mà2with D=50l m.This var-eD aThas the attraction that the calibration is directly tied to the mode I toughness of the material. iation based on J

As noted in Section1,it is not feasible to use the?ne scale computational model developed here for failure analysis of large structures,except possibly when the precise location of the crack path can be anticipated.The element size in the re-gion of fracture for relatively tough structural alloys will be in the range from tens to hundreds of microns.Thus,application of damage models of the present type will usually be restricted to the study of basic aspects of crack formation and to crack-ing in structural components and in metal forming and joining processes.A method being developed[19]that is capable of analyzing the failure of large plate and shell structures is the extended?nite element method(XFEM)wherein localizations and cracks occur within large elements(compared to plate thickness,for example)and aligned in any direction.In such coarse scale formulations the fracture process is usually represented by a cohesive zone representing the overall trac-tion–separation behavior averaged through the thickness of the plate or shell.The present?ne scale computational model can be used to generate the criterion for the propagation direction and the overall traction–separation relation required for implementing the XFEM model.

The extended Gurson model can also be used to study detailed aspects of crack formation and growth as illustrated by the cup-cone failure mode of the round tensile bar.However,to properly capture the transition from mode I to shear cracking, the?nite element mesh must be designed to produce elements with nearly unit aspect ratio at failure in the rupture-critical locations.To satisfy this criterion with rectangular elements,the initial element aspect ratio(width to height)must be taken to be about e3e f=2.For DH36,with e f%1.4,the required aspect ratio is about8.This value is consistent with that used for the mesh designs that most accurately predicted the transition in failure modes(Figs.14c and d).Even more challenging are the

three-dimensional aspects of the transition of a mode I through-crack in a plate to the mixed mode slant crack that emerges when the crack advance is extensive.As the crack advances,a neck forms ahead of the current crack tip,localizing the plastic deformation and developing into a slanted shear crack in the?nal stages of separation.As noted in connection with the cup-cone simulations,the prediction of a change in direction of crack path involving a transition from a mode I to a mixed mode separation process is quite sensitive to mesh design[10,18].Further effort is needed to create more robust predictive capa-bilities.A?ne scale XFEM formulation using the extended Gurson model to generate the details of the cohesive zone behav-ior would be worth exploring.

Acknowledgment

This work was supported in part by an ONR MURI through grants to the School of Engineering and Applied Sciences at Harvard University and the Materials Department at the University of California,Santa Barbara.

Appendix A.The remaining equations governing the modi?ed Gurson model and details of the numerical algorithm The remaining equations governing increments in the modi?ed model are now listed.Void nucleation is not included but it can readily be incorporated[13,20].The consistency condition for continued plastic loading,

_F?@F

@r ij _r ijt@F

@r M

_r Mt@F

_f?0;e10T

provides the expression for the hardening modulus,

h?àe1àfTP kktk x f x

P ij s ij

@f t

h M

e1àfTr M

@r M

P ij r ij

e11THere,

@F @r M ?à

2r2e

f r m

sin h

r m

2r M

e12T

?2q1cos h 3q

r m r M

à2q3fe13T

and h M is the modulus of the matrix material de?ned in terms of the logarithmic plastic strain and true stress in uniaxial tension as

1 h M ?

d e P M

d r M

e14T

The matrix material(i.e.,the undamaged material with f?0)is de?ned by its Young’s modulus,E,Poisson’s ratio,m,and relation between logarithmic plastic strain and true stress in uniaxial tension,e P Mer MT,also considered as the relation be-tween effective plastic strain and effective stress.These are inputs to the modi?ed Gurson Model along with the new param-eter k x and the initial value of f.As in the original model,plastic work in the matrix is equated to macroscopic plastic work according to

e1àfTr M_e P M?r ij D P ij;e15Tsuch that increments in matrix?ow stress can be computed from

_r M?

h M r ij D P ij

e1àfTr Me16T

The?nal step is to identify the stress rate for?nite strain applications and to combine the elastic and plastic strain incre-ments.The stress increments,_r ij,in the above development are identi?ed with objective Jaumann increments,whose Carte-sian components coincide with true stress increments for straining in axes parallel to principal stress axes.Void damage diminishes the overall elastic moduli of the material.However,this is a small effect compared to void in?uence on plastic behavior and the effect on elasticity is neglected,as usually done in this type of model.Isotropic elastic behavior is assumed. Combining elastic strain rates,D e

,and plastic strain rates from(6)gives the total strain rate as

D ij?M ijkl_r kle17Twith instantaneous compliances

M ijkl?1tm

ed ik d jltd il d jkTà

d ij d klt

P ij P kl

508Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509

The inverse is

_r ij?L ijkl D kle18Twith instantaneous moduli

L ijkl?L e

ijkl à

L e

ijmn

P mn P rs L e

rskl htP rs L e

rsmn

P mn

where the elastic moduli are

L e ijkl ?

1tm

ed ik d jltd il d jkTt

1à2m

d ij d kl

Plastic loading has been assumed in writing both(17)and(18);if the increment is elastic,only the elastic moduli and compliances are used.The effective plastic strain rate is de?ned in terms of the logarithmic strain rates in the usual way as

_e P e ?

???????????????????

2D P

D P

e19T

The?nal failure process beginning with the onset of coalescence and terminated by element deletion is modeled in the manner that has been commonly adopted[13]wherein the growth of the effective void volume fraction is accelerated when f>f c according to

f?efT?

f;f6f c

f ct1=q1àf c

f c

efàf cTf c

(

e20T

As detailed in[13],f is replaced by f?in the yield function(5)and in all the other equations except that_f in(8)remains unchanged.The material fails when f?f f.

A variety of numerical algorithms for the integration of elasto-plastic constitutive equations have been proposed in the literature[21–23].A class of backward Euler method has proven to lead to accurate and stable results[24]and is now widely used.Aravas[25]established the backward Euler scheme for pressure-dependent plasticity.Within the same framework,the integration algorithm for the present extended Gurson model has been derived.An alternative integration method recently developed in[26]employed all six stress components as independent variables to be solved simultaneously using Newton’s method.In contrast,the algorithm was re?ned in this work such that it is more ef?cient by employing only two independent variables,the increment of mean strain and the increment of effective strain.The algorithm was implemented into ABAQUS/ Explicit[16]through its user material subroutine interface(VUMAT).Several benchmark tests described in[26]have been performed to verify the code.

References

[1]Gurson A.Continuum theory of ductile rupture by void nucleation and growth.Part I—Yield criteria and?ow rules for porous ductile media.J Engng

Mater Technol1977;99:2–15.

[2]Rousselier G.Ductile fracture models and their potential in local approach of fracture.Nucl Engng Des1987;105:97–111.

[3]Howard IC,Li ZH,Bilby BA.Ductile crack growth predictions for large center cracked panels by damage modeling using3-D?nite element analysis.

Fatigue Fract Engng Mater Struct1994;17:959–69.

[4]Xia L,Shih CF.Ductile crack growth-I a numerical study using computational cells with microstructurally-based length scales.J Mech Phys Solids

1995;43:233–59.

[5]Xia L,Shih CF,Hutchinson JW.A computational approach to ductile crack growth under large scale yielding conditions.J Mech Phys Solids

1995;43(3):389–413.

[6]Steglich D,Brocks W.Micromechanical modelling of damage and fracture of ductile materials.Fatigue Fract Engng Mater Struct1998;21(10):1175–88.

[7]Gullerud AS,Xiaosheng G,Dodds RH,Haj-Ali R.Simulation of ductile crack growth using computational cells:numerical aspects.Engng Fract Mech

2000;66:65–92.

[8]Xue L,Wierzbicki T.Ductile fracture initiation and propagation modeling using damage plasticity theory.Engng Fract Mech2008;75:3276–93.

[9]Tvergaard V,Needleman A.Analysis of the cup-cone fracture in a round bar tensile bar.Acta Metall1984;32:157–69.

[10]Besson J,Steglich D,Brocks W.Modeling of crack growth in round bars and plane strain specimens.Int J Solids Struct2001;38:8259–84.

[11]Bron F,Besson J.Simulation of the ductile tearing for two grades of2024aluminum alloy thin sheets.Engng Fract Mech2006;73:1531–52.

[12]Nahshon K,Hutchinson JW.Modi?cation of the Gurson model for shear failure.Eur J Mech A/Solids2008;27:1–17.

[13]Tvergaard V.Material failure by void growth to coalescence.Adv Appl Mech1990;27:83–151.

[14]Bao Y,Wierzbicki T.On fracture locus in the equivalent strain and stress triaxiality space.Int J Mech Sci2004;46:81–98.

[15]Barsoum I,Faleskog J.Rupture in combined tension and shear:experiments.Int J Sol Struct2007;44:1768–86.

[16]Tvergaard V.Shear deformation of voids with contact modeled by internal pressure.Int J Mech Sci2008;50:1459–65.

[17]Nahshon K,Pontin MG,Evans AG,Hutchinson JW,Zok FW.Dynamic shear rupture of plates.J Mech Mater Struct2007;2:2049–66.

[18]Leblond J-B.Studies of the cup-cone fracture mode.Work in progress.Private communication;2009.

[19]Moes N,Dolbow J,Belytschko T.A?nite element method for crack growth without remeshing.Int J Numer Meth Engng1999;46:131–50.

[20]Chu CC,Needleman A.Void nucleation in biaxially stretched sheets.J Engng Mater Technol1980;102:249–56.

[21]Cris?eld MA.Non-linear?nite element analysis of solids and structures,vol.1.West Sussex(UK):Wiley;1991.

[22]Simo JC,Hughes https://www.wendangku.net/doc/685396535.html,putational inelasticity.New York(USA):Prentice-Hall;1998.

[23]Belytschko T,Liu WK,Moran B.Nonlinear?nite elements for continua and structures.West Sussex(UK):John Wiley and Sons Ltd.;2000.

[24]Ortiz M,Popov EP.Accuracy and stability of integration algorithms for elastoplastic constitutive relations.Int J Numer Meth Engng1985;21:1561–76.

[25]Aravas N.On the numerical integration of a class of pressure-dependent plasticity models.Int J Numer Meth Engng1987;24:1395–416.

[26]Nahshon K,Xue Z.A modi?ed Gurson model and its application to punch-out experiments.Engng Fract Mech2009;76–8:997–1009.

Z.Xue et al./Engineering Fracture Mechanics77(2010)492–509509