Model studies of ductile fracture—I. Formulation

Model Studiesof Ductile FractureI. Formulation*’ by s.

NEMAT-NASSER

and M. TAYA

Department of Civil Engineering, The Technological Institute, Northwestern University, Evanston, Illinois The paper is the first report on a systematic model study of ductile fracture, which has been initiated by the authors. The models include a matrix subdivided by grain boundaries, and containing hard inclusions. A systematic variational formulation for the rate problem is presented, and some of the inherent difficulties are discussed and illustrated. The effect of a large elastic compressibility is accurately accounted for. Results are compared with those of previous investigators.

ABSTRACT

I. Zntroduction Metallographical studies have revealed that the features of ductility of two-phase alloys (temperatures below 50% melting) which consist of a matrix subdivided by grain boundaries, and of hard inclusions, include: (1) formation and growth of voids around inclusions and/or at grain boundaries in slip bands; and (2) formation of intense shear bands between inclusions. The two processes are strongly coupled. Holes nucleate around inclusions, in order to relieve stress, but they are blunted as a result of intense local plastic flow of the matrix. The growth of holes is initially a stable process, but at a certain stage it may destabilize, resulting in a quick coalescence of voids. It is believed that the ductility is influenced by the shape, size, and distribution of the inclusions, by the three-dimensionality of the stress-state, and by the temperature. Recently there has been some interest in numerical simulation of this phenomenon; see, for example, Needleman (1)for a discussion and additional references. To better understand ductility in metals, a systematic numerical model study has been initiated by the authors. The model consists of a “unit cell” which is believed to represent “magnified” continuum particles. Since in most continuum theories the deformation is regarded to be “locally homogeneous”, it is reasonable to characterize the macroscopic material inelasticity by intense inhomogeneity of micro-deformations. The locally homogeneous deformation implies that the unit cell is subjected at its boundary to uniform displacements. Initial calculations are therefore performed for this model. The numerical formulation is based on a consistent and accurate description of the kinematics, the dynamics, and the constitutive relations, at finite strains and rotations, where, for the matrix, a JZ yield criterion is employed; see, for example, Hutchinson (2), and Nemat-Nasser (3). *This work has been supported GK-31352

to Northwestern

by the National University.

463

Science

Foundation

under

Grant

S. Nemat-Nasser and M. Taya In this paper, the formulation of the numerical procedure is presented, and the accuracy of the method is investigated by means of simple examples. Application to more complex problems will be reported in subsequent works. Section II summarizes basic kinematics, dynamics, and illustrative constitutive relations. Section III gives a brief account of the discretization procedure. In Sect. IV simple problems are solved numerically and analytically, providing a comparison for the assessment of the accuracy of the numerical method. This section also includes some results on the response of an elasto-plastic rectangular bar containing a circular hole at the center, and being subjected to uniaxial extension. ZZ. Continuum

Basis

The deformation can be expressed as x, = x,(X*, t), a, A = 1, 2, 3, which gives particle position x, in current (at time t) configuration c, in terms of that in the initial configuration C, the mapping being one-to-one and invertible, so that .I = det Ix,,*1 is finite and positive; comma followed by an index denotes partial differentiation with respect to the corresponding coordinate, and a fixed rectangular Cartesian coordinate system with the unit base vectors e,, is used. It is also useful to consider a convected coordinate embedded in the body, initially coincident with the fixed rectangular Cartesian coordinate system, and deforming with the body. The covariant base vectors of the convected system then are gA = &A&, repeated indices to be summed. Let u, = u,(x, t) be the velocity field, x = x,e,. Then the deformation rate tensor is (2.1)

91= Dabeneb = EABgAgB,

where Dab = $(u,,, + ub,?) are the components with respect to the fixed coordinate system, whereas E AB = Dabxa,A+ are covariant components with respect to the convected coordinate system; E AS is the rate of change of the Lagrangian strain, EAB = &(gA* gB - S,), 8, being the Kronecker delta. Let GJ be the rate of work per unit mass. Denote by p. and p the mass density in C and c, respectively. Then, if 3 is the Cauchy stress, 7 = (po/p)3 is the Kirchhoff stress, and !YR is the first Piola-Kirchhoff stress, we obtain

=

’

(2.2)

Tza3ia,A, PO

where 3 = Tabeaeb, 7 = ?&e&b = 7ABgAgB define the corresponding components. In (2.2), Tz,= (po/p)X &Tab are the components of the first PiolaKirchhoff stress, TR, in the fixed coordinate system. 464

Journal of The

Franklin Institute

Model Studies of Ductile Fracture-I

As an objective stress,

stress rate we use the Jaumann

derivative

of the Cauchy

aTab = ihb - WacTcb - WbcTca,

(2.3)

9t

where w& = 4(2)a.b- vb+). Then the Jaumann derivative of the Kirchhoff stress becomes (2.4)

c,c

where v~,~= -b/p. As has been pointed out by Hill (4), the rate problem concisely by the virtual work statement

JV

?&x

&,A

dV=

IV

pofh 6v, dV+

Is

can be expressed

?a 6v, dS,

(2.5)

where fh and ?a = N A?”Aa are prescribed rate Of body fOrCe (per Unit maSS) and surface traction (per unit initial area), NA being the components of the exterior unit normal in the initial configuration referred to the fixed coordinate system; V and S denote the initial volume and surface of the body, respectively. Direct calculation gives

$2,

3,

(2.6)

?‘a Sv, ds,

(2.7)

PO P Brab 6&A = - -----D,b-~TabS(2D,D,b-V,,V,b) P PO St

so that (2.5) can be replaced by

- -w

p srab

PO

1 I

aD,b -iTab ~(2DacDcb - vc,a&,b) dv =

”

pj‘, 6v, dv +

Is

where v and s are volume and surface of the body in current configuration c. This equation will be the basis of our numerical formulation. Except for the factor of p/pa, it is the same equation as that used in (5). Equation (2.7) is exact and coincides with that presented by Hill (5). For incompressible materials, or when small volume changes are involved, the factor p/p0 can be ignored. However, for large deformation of compressible materials, (2.7) must be used, otherwise significant errors may be accumulated. This is illustrated by means of examples, in Sect. IV. Finally, following Hill, we introduce the rate potential u(B) = ;CabceDabDce,

(2.8)

so that 9r,&Qt = aulaD,b, where C&e has different values for elastic unloading and plastic loading. For example, as has been discussed by Hutchinson (2), Vol. 302, Nos. 5 & 6, November/December

1976

465

S. Nemat-Nasser

and M. Taya

Needleman (l), and McMeeking and Rice (6), for Prandtl-Reuss mation one may set

type approxi-

where prime denotes the deviatoric part of the Cauchy stress, T2=$TLbThb, and (Y is unity if T)abDab2 0 and T= T,,,, otherwise (Yis zero, where T,,, is the maximum value of the equivalent Cauchy stress attained during the deformation history preceding the current deformation state. In (2.9), h = (l/E, - l/E)-l, where E is Young’s modulus and E, is the tangent modulus. Various aspects of this and related approximations are discussed in (6). With (2.8), (2.7) can be written as SJ(Eb)= 0, where J(B) =

I[c

p u(B) - Ta/,DacDcb+ ~T~&,,V~,b PO

- &?I,

1 Is dv -

f&,

ds,

(2.10)

and where the variation is on kinematically admissible velocity fields. For the numerical calculation we also will need an explicit one-dimensional relation between the true stress a, and the logarithmic strain E = In (l/lo), where lo and I are the initial and the final lengths, respectively. For illustration, we shall assume that

(2.11)

where uy is the yield stress. Hence we obtain E, = (E/n)(ay/u)n-l modulus.

for the tangent

HZ. Solution Procedure

We confine attention to two-dimensional problems. For numerical calculations we use (2.10) and triangular elements with a linear velocity field in each element. The procedure is standard, and will not be repeated here; see, for example, Desai and Abel (7). Denote the nodal velocities by t. Then, for arbitrary variation of these velocities, Eq. (2.7) reduces to (K+K+=Q,

(2.12)

where K is the matrix corresponding to the first term in the integrand of the left-hand side of (2.7), K, is the matrix corresponding to the second term, and Q is the nodal value of the rate of loading corresponding to the right-hand of Eq. (2.7). In (2.12) we have separated terms which are related to the stress field Tab. This will allow an easy estimate of the possible bifurcation which may take place. At a given loading state one may consider a proportional loading, ATab, use the value of the matrix C&e for that state, and estimate the value of 466

Journal of The Franklin Institute

Model Studies of Ductile Fracture-I A which leads to bifurcation, approximately from det jK+hKTl = 0; the corresponding eigenmode is then estimated from (K+hK=)f=O, where V must be normalized, e.g. V * f = 1. Let 8 be a monotonically increasing time-like parameter defined in such a manner that 8 = 0 corresponds to the initial configuration. We assume that all applied loads and prescribed displacements are defined as functions of 0. Hence all material time derivatives are interpreted as derivatives with respect to 8. Suppose that all field quantities are known at the current configuration which pertains to 8 = &,. Consider an increment, A@ resulting in an incremental change of all field quantities. In particular, we have (8) (2.13)

PO

If we set ApIp =-Au,,,, Au, = A&,, and retain terms linear in the incremental values (wherever convenient), we obtain Ae+Aw,,T,b+AWbcTca = Tdeo)+

fab

A&

(2.14)

where the quantity in the brackets is evaluated for 8 = eo, and p. is the initial mass-density.

ZV. Examples To check the accuracy of the numerical procedure, we now consider some simple examples. Similar examples have been worked out by Osias (9), using different constitutive relations which do not relate to the rate potential (2.8). Whereas (2.8) implies that the Jaumann derivative of the Kirchhoff stress is proportional to the deformation rate tensor, the proportionality factor being C &e, Osias assumes that the Jaumann derivative of the Cauchy stress is given by such a relation; this assumption leads to a non-symmetric stiffness matrix, and may have other less desirable consequences. Example 1. We consider a uniaxial extension of a rectangular block in plane strain or plane stress, and use constitutive Eq. (2.9) with (Y= 0. We set x1 = x, x2 = y, x3 = z, and consider velocity field: vl = cx ; u2 = key ; and o3 = 0 (for plane strain), kcz (for plane stress). Then, direct calculation gives: In A, = co; In A, = kc8; and In A, = 0 (plane strain), kc0 (plane streess), where A denotes stretch (final length divided by initial length) in the direction of the corresponding subscript, and 0 is the “load parameter”. Since the Jaumann derivative of the Kirchhoff stress in this case coincides with the corresponding material derivative, for plane strain we have: i 11=Ec[l+v(k-l)]/(l+v)(l-2~); iZ2= Ec[k + v(l- k)]/(l+ v)(l-2~); and i33 = Ecv(l+ k)/(l+ v)(l-2~). For plane stress, on the other hand, we obtain: ii1 = Ec(l+ vk)/(l- v’); iZ2= EC(Y + k)l(l- v’); and OS3 = -vc(l+ k)/(l- u). Hence we obtain, for plane Vol. 302, Nos. 5 & 6, Novemkr/December

1976

467

S. Nemat-Nasser and M. Taya

6

0

1.0

I.5

2.0

2.5

3.0

3.5

4.0

FIG.1. *., = A;/‘“-‘);

Tll = Eh~(l-zy)‘(l-y) In AJ(1 - v2), and P, = TiiA, = strain: EA;l In AJ(l- v’), where A, is the current cross-sectional area normal to the x-axis, and P, is the resultant force in the x-direction. For plane stress we obtain: A, = A, = A;“, T1* = EA?‘-l In A,, and P, = Eh;’ In A,. For both plane stress and plane strain, the maximum resultant force corresponds to A, = A:= e = 2.7183 . . . , which is independent of v. Figures 1 and 2 depict analytical and numerical results. In these figures the solid lines are the analytical results, the ,{px

(IO34 1 cm -6

-6

OV

1.0

I.5

2.0

2.5

3.0

3.5

4.0 A*=

+ 0

FIG.

468

2.

Journal of The Franklin Institute

Model Studies of Ductile Fracture-I TABLE

I

Numerical and Analytical Values of A: and P:, for Simple Extension

Plane strain Analytical Numerical solution solution *=f A: Y= 0.4 I v=f p: z-

v = 0.4

Plane stress Analytical Numerical solution solution

2.71 (2.71)

2.72

2.71 (2.71)

2.72

2.71 ((3.34))

2.72

2.71 ((3.01))

2.72

0.491 (0.493)

0.491

0.369 (0.370)

0.368

0.439 ((0.525))

0.438

0.369 ((0.407))

0.368

dots are the numerical results using Eq. (2.7), with two elements as shown, and the crosses correspond to numerical results ignoring the factor p/p0 in Eq. (2.7); i.e. the formulation in (6). As is seen, this has led to substantial numerical errors. We also point out that our results differ from those given by Osias. For example, his A: equals el”” for plane stress, and e(‘-“)‘” for plane strain. For v = f (incompressible material) these reduce to our results, since the Kirchhoff and Cauchy stresses coincide in such a case. Table I compares these results. In this table, numbers in parentheses are from Osias’ work, whereas numbers in double parentheses are obtained by neglecting p/p0 in (2.7). Example 2. We consider simple tension of elastoplastic material, and for simplicity assume that the slope of the V--E curve, i.e. h, is constant. For illustration we use the same material parameters as used by Osias (9), in order to be able to compare the results. The analytical solution is possible for the plane stress case. In the elastic range the solution is the one given in our previous example. At the yield limit, T i1 = cry, we calculate the corresponding 1, arriving at: my = extensions, using the solution of Example E(A$y’)2”-’ In Aj;Y’; and A,(y) = hey) z = (AiV)-“. For T, 12 uy, the constitutive relations for plane stress with constant h are

‘11=

(l-

i22 = (I-

Eu2)F

C(l+gT;22)Dll+(v_gT;1T;2)Dzzl,

E v2)F [(v - gTLTL)Dn+

D33=-;(i~+im),

9E g = 4hT;,’

(I+ gT:$Dzl, F= 1+E(5-4u) 4h (1-O

3 presents numerical and analytical results for the indicated values of the parameters. Since a very small volume change occurs in this case, the factor p/p0 does not have much effect on the accuracy of the results. Figure

Vol.302.Nos.5& 6,NovemberlDecember 1976

469

S. Nemat-Nasser

and M. Taya

.-6

1.0

I.1

1.2

1.3

FIG.

I.4

I.5 .Xx=+

0

3.

Other examples. We have also worked out two other examples, namely the simple shear of a cube, and simultaneous simple extension and rotation of a block. In both cases the numerical results are very close to the analytical ones, as in the other two previous examples. In order to save space, we shall therefore not report these results here, but simply note that our numerical procedure is such that it properly accounts for large rigid-body rotations.

V. Extension of a Rectangular Bar Hauing an Initially Circular Hole at the Center As our final example we consider a uniaxial extension of a rectangular bar containing an initially circular cylindrical hole (plane strain case) at its center. We confine our attention to the symmetric mode of deformation, and seek to analyze one quarter of the bar, as shown in Fig. 4. For the material properties

IO cm 1 FIG.

470

4.

Journal

of TheFranklin

Institute

Model Studies of Ductile Fracture-l u*= 0.00217

u*= 0.00157

El

El

FIG.5. we choose: E=SxlO’kg/cm*, aY=2x102kg/cm2, v=& and n=4. The model is divided into 110 elements containing 72 nodes, as shown in Fig. 4. The displacement is applied incrementally. Figures 5(a-d) show the development of the plastic zone around the hole, as the applied displacement at the end of the bar is increased; 2u* denotes the total displacement of the two ends of the bar, so that u*/lO is the average axial extension, i.e., change of length per unit initial length. Although the results for only one value of the Poisson ratio are shown, we note that this ratio has some effect on the development of the plastic zone. For example, a narrow plastic band at a 45” angle develops more clearly for v = 4, than for u = 4. We finally note that the pattern of the plastic zone and its evolution obtained in this work are similar to those reported by Needleman (l), and Miura and Asaoka (10). Discussion The incremental method used here is an up-to-date Lagrangian one which has some simplicity. But since the reference configuration is changed after each incremental step, some numerical errors accumulate. This is particularly pronounced, if the Cauchy stress and the Kirchhoff stress are equated on the ground that the current configuration is used as the reference one. For compressible materials the error may be unacceptable for large deformations. To avoid this difficulty, we refer the Kirchhoff stress, 7, to the initial configuration, while all the calculations are done on the basis of the current configuration. Our numerical results show that this procedure tends to improve the accuracy considerably. Finally we point out that the use of a rate constitutive relation in terms of the Kirchhoff stress, which admits a rate potential, is not only theoretically appealing and more sound, but also leads to a simpler numerical procedure. Finally we wish to point out that the conventional finite element procedure used here is not effective enough for the considered class of problems. An improved procedure will be reported in a subsequent paper.

Vol. 302, Nos. 5 & 6, November/December

1976

471

S. Nemat-Nasser

and M. Taya

References

(1)A. Needleman,

“Void growth in an elastic-plastic medium”, J. uppl. Me&., Vol. 39. p. 964, 1972. (2) J. W. Hutchinson, “Finite strain analysis of elastic-plastic solids and structures”, Numerical Solution of Nonlinear Structural Problems, Edited by R. F. Hartung, ASME, p. 17, 1973. (3) S. Nemat-Nasser, “Continuum bases for consistent numerical formulations of finite strains in elastic and inelastic structures”, Finite Element Analysis of Transient Nonlinear Structural Behavior, Edited by T. Belytschko et al. ASME, Vol. AMD-14, p. 85, 1975. (4) R. Hill, “Some basic principles in the mechanics of solids without a natural time”, J. Mech. Phys. Solids, Vol. 7, p. 209, 1959. (5) R. Hill, “Eigenmodel deformations in elastic/plastic continua”, J. Mech. Phys. Solids, Vol. 15, p. 371, 1967. (6) R. M. McMeeking and J. R. Rice, “Finite-element formulations for problems of large elastic-plastic deformation”, Inf. J. Solids Strut. Vol. 11, p. 601, 1975. (7) C. S. Desai and J. F. Abel, “Introduction to the Finite-Element Method”, Van Nostrand Reinhold, New York, 1972. (8) S. Nemat-Nasser and H. D. Shatoff, “A consistent numerical method for the solution of non-linear elasticity problems at finite strains”, SIAM J. appl. Math., Vol. 20, p. 462, 1971. solids”, NASA CR-2199, (9) J. R. Osias, “Finite deformation of elastic-plastic March, 1973. (10)I. Miura and K. Asaoka, “Simulation on the initiation and growth of ductile fracture voids”, J. Japan Inst. Metals, Vol. 39, p. 1025, 1975.

472

Journal of The Franklin

Institute