- Email: [email protected]
Crack-tip fields in a viscoplastic material
hr. 1. Solids S~ruc~ura Vol. 23, No. 7. pp. 1035-1052. 1987 Printed in Great Britain.
CRACK-TIP
0
FIELDS IN A VISCOPLASTIC
ooze-7683/87 s3.00+ .a0 1987 Pcrgamon Journals Ltd.
MATERIAL
J. D. ACHENBACH,N. NISHIMURA~and J. C. SUNG Department of Civil Engineering, Northwestern University, Evanston, IL 60201, U.S.A. (Received 26 May 1986 ; in revised form 13 October 1986)
Abstract-An integral representation for the particle velocity in terms of a Green’s function and certain linear combinations of the inelastic strain rates has been used in this paper, both for a numerical method to compute full field solutions and to develop unequivocal asymptotic expressions for the fields near a stationary crack tip. Specific rest&s have been obtained for a solid whose constitutive behavior is represented by the Rodner-Partom model. It is shown that the leading term of the near-tip particle velocity is of order r”*, and the higher-order terms are of the forms r log r and r. Expressions have been derived for the angular variations and for the multiplying timedependent intensity factors. The r log r term is absent for the Mode III case. Two questions have been addressed in further detail : the dependence of the intensity factors on time and the importance of the higher-order terms. The numerical results show a stress intensity factor which decays with time. At a small distance from the crack tip the numerically computed normal&d opening stress r22(xI.0, ~)/T~~(x,,O,O)has been compared with a one-term asymptotic representation, i.e. with &(f)/&(O). The two curves diverge at very small values of time. The inclusion of a second term in the asymptotic expression for the stress gives very good agreement over a larger time range.
1. INTRODUCTION
A clear understanding of the mathematical structure of the fields of stress and deformation near a stationary crack tip is an essential prerequisite to an investigation of conditions for the onset of crack propagation. In a viscoplastic solid the fields in the immediate vicinity of a crack tip can be very intricate. In this paper near-tip viscoplastic fields have been investigated in detail for a class of constitutive models for which the inelastic strain rate approaches a bounded value as the stresses increase beyond bounds. The Bodner-Partom model[l], which has been shown to describe the mechanical response of a wide class of materials, displays this particular behavior of the inelastic strain rate. Fields near a crack tip are usually analyzed asymptotically by assuming a near-tip field of the general form r’f(O). When this expression is substituted into the governing equations and into appropriate boundary conditions, and terms of equal order in r are collected, there results a non-linear eigenvalue problem, which yields both p and S(O). For viscoelastic solids this method was used by Riedel[2] and Riedel and Rice[3], who investigated the neartip fields of stationary cracks for Mode III and Mode I, respectively. They considered a viscoelastic constitutive equation of the power law type. For the Mode I case their asymptotic solution displays the well-known HRR singularity if the power is greater than one. When the power is less than one, the solution becomes square-root singular. Hui and Riedel[4] considered quasistatically moving cracks for both Mode I and Mode III conditions. When the power is greater than 3, their singularity is different from both the HRR and the square-root singular forms, while for a power less than 3 it reduces to a square-root singularity. Lo[5] followed the Hui-Riedel approach, but he took into account the effect of inertia. The implications of square-root singular crack-tip fields for viscoplastic solids have been explored by Hart[6] and Freund and Hutchinson[‘l]. For a Mode III quasi-static steady-state analysis, Hart[6] related the stress intensity factor at the tip to the corresponding far-field quantity. He also discussed transient cases on the basis of certain assumptions. Freund and Hutchinson[7] developed a similar relation between the tip and the far-field stress intensity factors for the dynamic steady-state Mode I case, using an essentially linear relation between the stress and the inelastic strain rate. In this paper a more fundamental approach towards the computation of near-tip fields is presented, in the sense that asymptotic forms are not assumed a priori. Instead, the t Present address : School of Civil Engineering, Kyoto University, Kyoto 606, Japan. 1035
1036
J.
D.
ACHl%MCli
cl
d.
‘5 0 Fig. 1. Crack-tip geometry. particle velocity field is expressed by a representation integral over the complete domain of the cracked body. This integral involves an elastic Green’s function and linear combinations of the inelastic strain rates. Since the inelastic strain rates are bounded at the crack tip, the integral can be expanded in terms of the distance to the crack tip, to yield an unequivocal expression for the near-tip field. Details have been worked out for a stationary crack tip. It is shown that the leading term of the particle velocity is of order r”‘, and the next terms are of the forms r log r and r. Expressions have been derived for the angular variations and for the multiplying time-dependent intensity factors. The r log r term is absent for the Mode III case. Two points are of interest : the dependence of the multiplying terms on time and the importance of the higher-order terms. The representation of the particle velocity in terms of an appropriate Green’s function and the rates of change of time of the inelastic strains has not only been used to obtain near-tip asymptotic expansions of relevant field variables, but also to develop numerical solutions by the use of the boundary integral equation method, in conjunction with an iteration technique. The numerical calculations have been carried out under the assumption of small-scale yielding. Thus, the equations of the Bodner-Partom model have been used in a region around the crack tip. The conditions on the boundary of this region are provided by the singular term of a corresponding elastic solution. The results show a stress intensity factor which decays with time. At a small distance from the tip the numerically computed normalized opening stress t22(~,,0, 1)/r22(~,r0,0) has been compared with a one-term asymptotic representation, i.e. with K,(r)/K,(O). The two curves diverge at very small values of time. The inclusion of a second term in the asymptotic expression for the stress gives very acceptable agreement over a longer period of time. As early as 1971 the boundary integral equation method (BIEM), also called the boundary element method (BEM), was applied by Swedlow and Cruse[8] for elastoplastic material behavior. In recent years the method has been used extensively for the numerical analyses of viscoelastic/plastic solids of various shapes including bodies containing cracks. Recent work was summarized by Mukherjee[9]. 2. FORMULATION
The states of deformation that will be considered are two-dimensional, either plane or antiplane strain. The body contains a straight crack. At one of the crack tips, say A ,, the applied loads give rise to elastic stress intensity factors KF(r), K:(t) and K:,(t), if the material were linearly elastic. We will, however, consider a material that displays inelastic behavior as the stress level increases. Hence there will be a domain of inelastic deformation near the crack tip. It is, however, assumed that this domain is small as compared to the crack length, and that linear elasticity applies outside the near-tip domain of inelasticity. To analyze the field right near the tip A i, we follow the arguments that have been introduced for small-scale yielding. Thus, we magnify the geometrical scale near A ,, so that the geometry becomes one of a semi-infinite crack in a full space, as shown in Fig. 1. A Cartesian coordinate system is centered at the crack tip. For the plane strain case, possible boundary conditions on the bounding surfaces of the body are replaced by the asymptotic condition
u --, g
[K~(t)f,(e)+K~(r)f,,(e)l
as r -+ co, where p is the shear modulus, and (r, 0) are polar coordinates centered at the tip
Crack-tip fields in a viscoplastic material
1037
in the magnified configuration. Note that the definition of the stress intensity factors differs from the conventional one by the factor (27r)1’2. The functions f, and f,, give the angular variations of the Cartesian components of the elastic near-tip displacements for Modes I and II, i.e. -
-
cos
cos
3012 4) 012 4 f,W = (K (K+ 4) sin 8/2 - 4 sin 38/2
fl,(Q = where
K
(K+ {) sin O/2+ 4 sin 30/2 (-K+;) cos 8/2-t cos 38/2
(2a)
(2b)
is a constant defined in terms of Lame’s constants 1 and p 1+3/A
(3)
K=F. We now proceed to the governing equation for the problem. The strain rate i =
“,[Vti+ (Vti)‘],
where (‘) = d/df
(4)
is usually decomposed into the elastic part Ccand the inelastic part d’ as c = $+i’.
(5)
The elastic part is related to the stress rate + by Hooke’s law +=I1
trdc+2@
(6)
where 1 is the unit tensor and “tr” is the 3-D trace. Since the crack is stationary, the effect of inertia is not considered. Hence, the equilibrium equation is
v-i = 0.
(7)
An equivalent form of eqn (7) in terms of the displacement rate ti follows from eqns (4)(7) as A.*i = V.,j’
(8)
A* = pV2+(A+p)VV-
(9)
where the operator A* is defined as
and for convenience we have introduced the quantity ti.‘,defined by t+’= 11 tr i’+2$.
(10)
Equation (8) is the governing equation for our problem. This equation, together with the boundary conditions i:=+*n=O
on
x2=0*,
x,x0
(11)
and an asymptotic condition of the form of eqn (1) with u and @iI replaced by ti and k:,,, defines ri. When k:,, = 0, the behavior of the elastic solution away from the crack tip becomes
J. D.
1038
ACHENBACH
Ii *
O(r-
al.
et
‘i2).
(12)
To complete the formulation, we supplement these equations with the regularity requirements at the crack tip ir - bounded as r J 0
(13)
+ - 0(1/r)
(14)
as
rl0.
In the present problem, we will also discuss the Mode III case, where T,~ = r33 = u, = 0 (tl, /I = 1,2). We prefer to use the stress function 4 because it makes the stress computation easier. The stress function defines the stress components r‘. : = 7.3
(a = 1,2)
(15)
as r. = e,fl Q#J
($4 : = WW
(16)
where em8is the permutation symbol, eaB= e,,,. The compatibility condition e,&& = 0
(E, : = &,j)
(17)
together with eqns (5), (6) and (16), yields the governing equation
The boundary condition on the crack faces and the regularity condition analogous to requirement (14) yield i = const
(19)
on the crack faces with the constants for both faces being the same. Also, on the crack line, we have 5, = a,4 = 0.
(20)
Finally, the asymptotic condition is 4 --) -2r”*k,E,,(t)
cos 46 as
r + co.
(21)
When X:,(t) = 0, condition (21) is replaced by d; - O(r-‘I*)
as
r + co.
(22)
In this case, the constant in eqn (19) is zero because d(O) = 0. This follows from m <*dx,=-
o= s0
Oq$4 dx, = &o>-&co> I0
where eqn (22) has been used. In the sequel, we will consider only the cases where
= d(O)
(23)
1039
Crack-tip fields in ;Lviscoplastic material kfi,,.,*,w
3. REPRESENTATION
(24)
= 0.
OF THE VELOCITY FIELD
The term on the right-hand side of the displacement-rate equation, eqn (8) is of the nature of a body-force distribution. This implies that the solution to cqn (8) can bc cxprcsscd by the use of a Green’s function. Representations in terms of Green’s functions are at the basis of the boundary integral equation method (BIEM), and they have therefore become quite well known. For rate-type inelastic analysis the appropriate representations were first developed by Swedlow and Cruse[8]. This section gives a brief derivation of such representations for the purpose of making the paper self-contained. For further detail we refer to Mukherjee[9] and the papers cited there. For an unbounded body containing a semi-infinite crack, the Green’s function C(x, y) is defined by A,+G(x,y) = -18(x-y)
(25)
where x = (xi, xJ, y = (y ,, yJ. The operator A* is defined by eqn (9), and the subscript indicates differentiations with respect to xi, i = 1,2. On the faces of the crack (xi < 0, x2 = O*) we have U+LY)
=
(26)
0
where T,G stands for the elastic tractions computed from G. We also have
W, y) -
Wag
1x1) as as
V,W,Y)
- WQ/lxl)
Ix1+
(27)
00
Ix]+ 0
(28)
as 1x1-+ 0.
(29)
A standard analysis based on Muskhelishvili’s method yields C(x, y) as
-K
log
(,/Z-&
($--$I+-~
log (&+&,)-log
)+&(z-Z”
+
(++&o)} +=)I] (30)
where
K
is defined by eqn (3), and
i-/(-l),
z=xl+ix2,
z. =y,+iy2,
5=x,-ix2,
f. =y,-iy,.
(31)
Equation (30) agrees with the expressions given by Erdogan[lO]. A suitable representation of the solution to eqn (8) is obtained by the use of Green’s identity
1040
J. D. ACHENBACH er a/.
Fig. 2. Domain
(A:v) * w dA, -
da,
for use of integral
*
(A:w) v dA, = s d&z
identity.
(T,v)* w dc(T,w)*v dc sad.?, s ddn,,
(32)
which holds for an arbitrary pair (v, w), and whcrc dQR,ris the domain bounded by circles of radii E and R (0 < E < R) centered at the tip, and by the crack faces, as shown in Fig. 2. We now substitute ii(x) for v(x) and G(x,y) for w(x) in eqn (32). Subsequent use of the divergence theorem and the regularity conditions (13), (14), (28) and (29), followed by letting R + co, E + 0, and using conditions (12) and (27), gives
it(x, t) =
s
V$(Y,
R2
4 *b'(y, t) dA,
(33)
or in component form
tij(X, t) =
VykG,j(YvX)~h(Yt 1)dAy.
(34)
Here V, = i(d/8y,)+j(d/&), and dA, = dy, dy2. In deriving these equations we have assumed that the integrals in eqn (33) or eqn (34) converge. For this to be true, it is sufficient to assume that 6’ satisfies (35)
lti’l - O(Pl)
as
r + co
(36)
with al > -1,
a2 < -1.
(37a, b)
Although inequality (37a) appears unnecessarily restrictive, it turns out that these inequalities guarantee the satisfaction of conditions (13) and (14) automatically. It is noted, however, that inequalities (37) are just sufficient conditions. For the Mode III case a similar procedure gives
d(X> 4= where
l2a,,
G(y, x)e,&~,
0 d4
(38)
Crack-tip fields in a viscoplastic matcrial
W,Y)
(& - &0) ($-
= - $ log
(Jz+J&l)
&I)
1041
(39)
(J~+JzlJ).
Hence we have r$(x, 1) = epqax,
a,,WY,
4. ASYMPTOTIC
x)e,d(y,
0 dA,.
(40)
EXPANSIONS
In the following sections, we employ the potential representation for i given in Section 3. The results of this section are analytical, while numerical results are presented in Section 7. For a certain class of constitutive models the quantity b’, given by eqn (lo), is finite. For example, a viscoplastic model proposed by Bodner and Partom[l] shows this property. In this section we will obtain the asymptotic expression for ti for the case that h,$i&(r, 0,.t) exists for all 8. In a conventional asymptotic analysis one seeks near-tip fields of the forms
where ~j, tii (j = 1,2 , . . .) are functions which are increasingly less singular near r = 0 with increasing subscript. Since near r = 0 ii =
d’(0,8, t)+o(l)
(43)
is bounded, and
we find by substituting expressions (41) and (42) into eqn (8) A*u, = 0 A*irz-V~t$,
(45) = 0.
(46)
The operator A*, which is defined by eqn (9), is here used in polar coordinates. Equation (45) implies that u, is of the same form as for the linearly elastic case. The appropriate separation of variables form of eqn (46) can only be guessed at this stage. One might try til =
tqe,t).
(47)
This expression for I(& r) does, however, not necessarily lead to a solution of eqn (46) because the corresponding homogeneous problem has non-trivial solutions. For the Mode I case consider, for example, the homogeneous equation A*[rf(e, t)] = 0 with boundary conditions
(48)
J. D. ACHKNBACH et (11.
1042
0 =
T[rl(e, f)] = 0
7[:
i, * T[rl(e, r)] = 0
0 = 0:
and
(49)
i, * [rr(6, t)] = 0
(jOa, b)
where the operator T has the same meaning as in eqn (26), but in polar coordinates, and i, are unit vectors. Equation (48) has a solution corresponding to a uniaxial tensile field in the xi-direction. Hence the approach based on eqn (47) becomes ambiguous at this point. An unequivocal form of the term it2 can, however, be obtained directly from the representation given by eqn (34). Using polar coordinates (r, 6) and (ro, 0,) for y and x, we split cqn (33) into two inlcgrals
ss n
Co, 00,O=
‘0
--I 0
rVG(r, 8; ro, 8,) -S’(r, 0, t) dr de
rVG(r, 19;ro, 0,) *ci.‘(r, 0, I) dr df?.
+
(51)
Here the gradient operator is taken with respect to the polar coordinates r and 8. Next VG is expanded into a power series of r. We can write r < r. :
Mnc I (4 00)
(52)
(n+ I)/2
r > r. :
rVG=
f n--l
0 y
N,+ I (6 eo>.
(53)
The coefficients in eqns (52) and (53) can be obtained from eqn (30). It is found that No is independent of do. Substitution of eqns (52) and (53) into eqn (51) yields after some rearrangements Co,
00)
=
f0 +tl(e0)rY2
+i,
(r0,0,)+i,
(ro, 0,)
(54)
where
ss n
i,
=
m
--II 0
N,(O) -d’(r, 0) dr de
(55)
r- ‘/‘N, (0,&J - ci’(r, 0) dr dt7
(56)
*tf’(r, 0) dr de
(57)
The dependence on time t, shown in eqn (51) is implied; but not explicitly shown in eqns (54)-(58) and in the sequel. In the next step the orders in r. of the integrals (57) and (58) are determined.
Crack-tip fields in a viscoplastic material
1043
Since ti’(r, 0) has an estimate given by eqn (43) it follows that eqn (57) is of order ro, i.e.
where
A@,) =
1 -No(e)-2N1(8, 80)
M,, #,&)
(60)
A discussion of the order of the integral in eqn (58) is given in the Appendix. The result is
L (ro,0,) = i,(e,)r, log To +roWo)+ 4ro)
(61)
where
t2(eo) = -
s’
N2(e, eopii(o,
e) de.
(62)
--I
The function B(O,) is defined by eqn (A6). By combining eqns (54), (59) and (61) it is concluded that ri(r,, 6,) is of the form
ii(ro,eO) = ~o+i,(~o)r~~2+~2(~0)r0 log ro+~3(Oo)ro+o(r0)
(63)
where
i,(e,)
=
A(e,)+B(e,).
(64)
The constant to in eqn (63) corresponds to a rigid-body motion ; it will be left out from further considerations. Some additional manipulations reduce i,(6,) and i,(O,) to the forms f,
(0,) = [&f,(0,)+ li;,f,,(e,w~
(65)
where f,(e,) and f,,(S,) are defined by eqns (2a) and (2b), and
r
- ‘/2D,,,I- d’(r, 8, t) dr de.
(67)
E,,r, * c+(O, 8, t) de.
w-0
Also
In eqns (67) and (68)
J. D. ACHENBACH EI ~1
1044
70 30 3 3 - - sm T + - sm 2 ‘2 2
D, = i 3
38
2
2
- -SK-
3
78
2
2
+ -sm-
1 38 3 78 K - - sin- - -sm2 222
( 1 1
-cos2222 E =
’
(
36
3 78 + -cos-
(i-K)C0S28-220s48
E,, = (1+x)
--
-
(I-K)
sin 48
-2
cos-
~(COS 2e-cos
4e)
30
3
.3e
21
sm-+-an
3.
~(COS 2e-cos (K+
1)
(69b)
22
-2 sin 48 cos 20+2 cos 40>
sin 2e+2 sin 48
(~-3)
3 70 + -cos 2 2 2
(69a)
sin 28-2
(70a)
48)
sin 48
G’Ob)
It is shown in the Appendix that i,(e) can be written in the form
where
P(e, e,) eciyo,e) =
4np(~ +Kpl(e0)(2~
w (e-O,)Hl(e)+H,(0)).l’(O,u)
+g,um2n w (e-O,)H,(e)+H,(u)).~‘(o,a) -t-g,ted (2~sgn(0- hdE, I (0))- ~i(o,0)
+h(ed (2n w (6-edE, (eww4 0)
+gs(eo)(E,(~).~i(O,~))+g,(~O)(El,(~)*~i(O,e))l
(72)
and
(73)
In eqn (72)
1
g,(e) = ___ 2(1+~)
cos e [ -sin e
e1 1' gm=[co: (74a-f)
Crack-tip
C
-KCOS 2u-cos
fields in a viscoplastic
40
HA@= -((h.+ I) sin 20)/2-sin
-((K+ 40
cos 20-cos
H4(e) _(K+ +(K+ i)o
2 cos 40 sin 40 (I -K) sin 20-2
1)
sin *O-sin 48
-(K-3)
2
1045
I) sin *Q/*--sin cos 2O+cos40
40
1
1
(I -K) sin 20+2 2 cos 48 sin 48
- cos 28 +
cos
48
1
1
48
-sin 28 + sin 48
--OS 2e+c0s 48
2
material
COS w-2
COS 48
2 sin (20 -sin 48)
*(sin 28 -sin 48)
-(K-t-i)
COS 28-e
COS 48
1 ’
(75a-d)
In the expression for c,,,,, we have denoted by pf the finite part of the integral defined in the Appendix. Equation (63) implies that the dominant singularity of the strain rate has exactly the same form as in linear elasticity. The rates of the stress intensity factors &;.,i are given by eqn (67). The second singular term is of the order of log r near the tip. The magnitude of the multiplying term is determined only by the distribution of di at the tip. Because of this last property, expression (68) holds true for any crack configuration, whereas eqn (67) is true only in the sense of the “small-scale yielding”. The logarithmic singularity in the strain rate does not give rise to a logarithmic stress singularity except in T, , for Mode I. For the Mode III case, we can use the same general approach to show that
(76)
where I” K,, =
4x s n
d0
R’ S[0
-ri’,(r,
30
30
0) sin -* + h\(r,O) cos y
1
r- I” dr
(77)
(‘0
li(0,) =
OS ;
0
[C’*$(O,O) cos 0-C.C?(O,O) -
sin 81 d0
psi(o,e) sin e+e+(o,e)
cos e] de-[e7ii(o,e,)]e
(78)
where
(7% b) Note that the stress rate and hence the strain rate do not have logarithmic singularities. As for the Mode I case, we can relate the absence of the logarithmic term either to the nonexistence of a non-trivial solution to a certain boundary value problem or to the lack of a certain term in the far field expansion of the Green’s function.
J. D. ACHENBACH
1046
er al.
In the published literature, the choice of terms for asymptotic expansions is generally made on a trial and error basis. The analysis presented in this section gives a more rational foundation to the investigation of the asymptotic near-tip structure.
5. THE BODNER-PARTOM
MODEL
In Ref. [l], Bodner and Partom proposed a set of constitutive equations for an elastic-viscoplastic material. These equations have the convenient property that no separate specification of a yield criterion is required, nor is it necessary to consider loading and unloading separately. Within the context of the Bodner-Partom model both elastic and inelastic deformations are present at all stages of loading and unloading, but the inelastic deformations are very small at low stress levels. The Bodner-Partom equations have been used to model the mechanical behavior of several metals over a wide range of temperatures and strain rates. Very satisfactory agreement between predicted and experimental results has generally been obtained. For a recent discussion of the Bodner-Partom model we refer to Ref. [ 121. The constitutive equations follow by supplementing eqns (5) and (6) with an expression for the rate of inelastic strain. In a recent version of the Bodner-Partom model this strain rate is given by, see Ref. [ 121
ii =
DoL exp [ - (Z2/3J2)“/2]
JJ~
(80)
where T’ is the stress deviator f’ =
tr
f-f1
f
(81)
and
Also Z = Z, +(Z,-Z,)
eemWp.
(83)
In eqns (80) and (83), Zo, Z,, m, n and Do are material constants, while W, is the work over the inelastic strains
W, =
J
t-8
dt.
(84)
For a uniaxial tension of magnitude cr, it is simple to plot the inelastic strain rate vs a/Z. For various values of n the curves are shown in Fig. 3. The linear plot of Fig. 3 shows that d’ remains very small at small stresses with a rapid change over a small range of a/Z (the “yield” region), and that an asymptotic limit is approached as Q becomes large. It is of interest to note that the components of ii, as defined by eqn (80), are always bounded. Even if r is singular (which is the case at a crack tip), the singularity of r as r + 0 will be cancelled by the singularity of ,/J2 in the term r’/ JJ2. The term exp [ - (z2/3J2)“/2] is bounded whether J2 is singular or not. Hence eqn (80) implies that the inelastic strain rate is bounded at a crack tip. It then follows from eqn (10) that the components of rY are also bounded at a crack tip, and it can thus be concluded that the asymptotic results of Section 4 apply to the Bodner-Partom model.
1047
Crack-tip fields in a viscoplastic material
UNIAXIAL
CASE
01 0
.2
.4
.6
.8
1.0
1.2
1.4
1.6
1.8
2.0
Fig. 3. Inelastic strain rates for uniaxial tension, u, and various values of II.
6. NUMERICAL
PROCEDURE
A numerical procedure has been developed to compute the stress intensity factor and the opening stress r 22 at a position on the crack line. As pointed out by Mukherjee and Kumar[l3], a numerical evaluation based on eqn (33) is very efficient for viscoplastic problems. Here we use this approach (generally referred to as the BIE method) for cracktip fields governed by the Bodner-Partom equations. The procedure is as follows. (1) Compute the initial field. Since the inelastic strain is zero at t = 0, the initial field is simply the linearly elastic one. In the present context this field is
Hence, the solution is known for t = nAt with n = 0, where At is the time increment. (2) Subdivide the near-tip domain into cells, and compute d’ at the centroid of each cell by using eqn (80). Then assuming that ii may be taken as constant over each cell, compute li by using (33) and WIby differentiating eqn (33) (3) Compute + from eqns (5) and (6). Update r by using r[(n+ l)At] = r(nAt)+At+(nAt).
(86)
(4) Use the updated value oft to compute Pi as in (2), and repeat (2) and (3). When the edges of the cells are straight, the integrals involved in eqn (33) can be evaluated analytically, as discussed by Cruse and Polch[ 111. By picking up the square-root singular terms in these integrals and summing up their coefficients, k,;.,, can be computed. Alternatively, the computation of g ,,,! may be carried out by using eqn (67) with the piecewise constant approximation for i’.
7. NUMERICAL
RESULTS
Numerical computations have been carried out for the non-hardening case of the Bodner-Partom model. With reference to eqns (80) and (83), this case corresponds to
1048
J. D.
ACHENBACH
er
al.
Fig. 4. Cell configuration [or application of boundary integral equation method
zo/z, =
1.
(87)
It is then convenient to introduce non-dimensional quantities defined as u = u(WJ3)/(W> f = X(Z2/3)/(K:)r,
T = r/V/,/3) f = Wo)WJ3)
(gga, b) Wa, b)
(90) Here the dot notation represents the derivative with respect to f. In terms of these dimensionless quantities, eqns (80), (6) and (8) convert into 8’ = [T’/($T’*T’)“2]
exp [-i(iT’*T’)-“1
(93) where 0 and v2 are the de1 operator and the Laplacian with respect to the dimensionless spatial coordinates fi, and T’ is the dimensionless stress deviator. Since n/p and (n+p)>/p can be expressed in terms of the elastic Poisson’s ratio, v, it now follows that numerical values have to be chosen only for v and n. Here we choose n = 1 and
v = 0.2.
(94)
For the numerical procedure based on the boundary integral equation method described in Section 6, the cell arrangement shown in Fig. 4 has been used. The expression for the rate of change of time of the stress intensity factor is given by eqn (67). Once the inelastic strain rate has been computed, the term tij can be determined and eqn (67) can be evaluated numerically. A subsequent numerical integration with respect to time produces
Crack-tip
0
Fig. 5. Normalized
2
4
fields in a viscoplastic
6
8
10
1049
material
12
opening stress at f = 0.0014 and ratio E,(r) dimensionless time, I.
14
76
= K,(f)/K,(O).
78
20
both as function
of
&(Q Figure 5 shows Ki(Q/Xr vs the dimensionless time r. It is noted that Ki@)/il~ decreases with time. Also shown in Fig. 5 is the normalized stress T&,, 0, T)/T&i, O,O), evaluated at R , = 0.0014. The two curves diverge rapidly as the dimensionless time increases. It is instructive to compute the actual physical distance to the crack tip which would correspond to di = 0.0014. Suppose we consider a crack in an unbounded medium under far-field uniform tension 222 = co. Then Kf = ~*~(~/2). For titanium, Ref. f 1) suggests the following value for the yield stress, by N 250 N mm-*, while Z - 1200 N mm-*. Let cro N ~rv/2. Equation (89) then yields xl N 1.63 x IO- * ~2,. Hence 2, = 0.0014 corresponds to x, N 2.279 x lo- 5u, which is a very short distance to the crack tip. The two curves in Fig. 5 thus imply that even very close to the crack tip, the singular term by itself would give a very poor approximation to the near-tip stress T**, with an error which increases with time f. The relation between physical time f and dimensionless time fcan be obtained from eqn (89) as I = 0.635 x IO?, where the following values from Ref. [I] have been used DON lo4 s- ’ and p= 0.44 x lo5 N mm-*. As predicted by the constitutive model used in this paper, the relaxation of the stress intensity factor takes place very rapidly. An improved asymptotic representation of the near-tip stresses can be obtained by including higher-order terms. The calculation of these terms requires the evaluation of the integrals for i,, eqn (68), and f,(t),), eqn (71). The r log r term in eqn (63) contributes only a constant term to fz2. A constant contribution is also obtained from the term i,(O,)ro. After the calculations have been completed we find on the crack line (0 = 0)
T 22
-
4r,
- 3 376
Js
-
or
(96) This two-term expansion gives excellent agreement with the numerical result up to I= 4, as can be seen from the numerical values listed in Table 1. The agreement deteriorates rapidly for f > 4. Figure 6 shows ~~~,(~)~~~~vs the dimensionless time f. Also shown is the ratio T2(Z1, 0, f)/T2(Z1, 0,O) evaluated at 2, = 0.0047. As for the Mode I case, the two curves diverge rapidly as f increases. It is relatively simple to write out explicitly the asymptotic results for the Mode III case. At the crack tip we obtain from eqn (91)
J. D.
1050
ACHESBACH er al.
Table I, Comparison of asymptotic rcprcscnfatlons with numerical results calculated by the Mode 1 722WJ3)
Mode III r,/(Gy’3)
WD, Z/J3
Two terms
I
13.04
1 L
II.05
BIEM
0
I
0
2
4
Two terms
Error (%)
13.10 II.17 9.12 7.12 6.73 6.35 5.97 5.60 5.235 4.87 4.52 4.17 3.84 3.504
8.93 6.84 6.43 6.02 5.62 5.21 4.8 1 4.42 4.03 3.63 3.25 2.87
3 4 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0
BIEM
22.18 18.78 14.47 10.08 9. I94 8.32 7.444 6.555 5.68
3.93 4.5 5.2 5.9 6.96 8.0 9.24 10.84 12.95 15.36 18.0
BIEM
Error (%)
22.28 18.82 14.35 10.42 9.745 9.104 8.502 7.94 7.41
3.26 5.65 8.6 12.4 17.4 23.3
I
6
8
10
12
14
16
18
20
E Fig. 6. Normalized anti-plane shear stress at f = 0.0047 and ratio R,,,(r) = K,,,(I)/K,,,(O), both as functions of dimensionless time, I.
Ei(o,q
=
(-sin:e4’). cos
(97)
Use of this result in eqns (76)-(78), and subsequent integration with respect to time yields
The values of T2 at 2, = 0.0047, Z2 = 0 computed according to eqn (98) and by the BIEM have been listed in Table 1. Good agreement is observed up to I = 4. The results presented in this section show that the square root singular term presents a very poor approximation to the crack-tip field even at very small distances from the crack tip. The approximation becomes worse as time increases. The use of additional terms in the asymptotic expression for the particle velocity increases the length of time for agreement between the asymptotic representation and the numerically computed stress component at a crack-line point just ahead of the crack tip. Acknowledgemenr-This work was carried out in the course of research sponsored by the U.S. Of&x of Naval Research (Contract No. NOO014-76-C-0063).
1051
Crack-tip fields in a viscoplastic material REFERENCES
1. S. R. Bodner and Y. Partom, Constitutive equations for elastic-viscoplastic strain-hardening materials. J. Appl. Mech. 42,385-389 (1975). 2. H. Riedel, Cracks loaded in anti-plane shear under creep conditions. Z. Mefall. 69(12), 755-760 (1978). 3. H. Riedel and J. R. Rice, Tensile cracks in creeping solids, Frucfure Mechanics (Praceedings oj the 12th National Symposium on Fracture Mechanics, St. Lous. May 1979) (Edited by P. C. Paris), ASTM SIT 700, pp. 112-130. American Society for Testing and Materials, Philadelphia (1979). 4. C. Y. Hui and H. Riedel, The asymptotic stress and strain field near the tip of a growing crack under creep conditions. Int. /. Fracture 17(4), 409425 (1981). 5. K. K. Lo, Dynamic crack-tip fields in rate-sensitive solids. J. Mech. Phys. Solidc 31(4), 287-305 (1983). 6. E. W. Hart, A theory for stable crack extension rates in ductile materials. Inf. /. Solidr Sfrucfures 16, 807823 (1980). 7. L. B. Freund and J. W. Hutchinson, High strain-rate crack growth in ratedependent plastic solids. /. Mech. Phys. Solids 33(2), 169-191 (1985). 8. J. L. Swedlow and T. A. Cruse, Formulation of boundary-integral equations for three-dimensional elastoplastic flow. Inr. J. Solids Strucfures 7, 1673-1683 (1971). 9. S. Mukhejee, Boundary Element Merhoa!s in Creep and Fracture. Applied Science, London (1982). 10. F. Erdogan, On the stress distribution in plates with collinear cuts under arbitrary loads, Proceedings o/lfh U.S. National Congress of Applied Mechanics (Edited by R. M. Rosenberg), pp. 547-553, American Society of Mechanical Engineers, New York (1962). I 1. T. A. Cruse and E. Z. Polch, Elastoplastic BIE analysis of cracked plates and related problems, Inl. 1. Numer. Me/h. Engng X3,429-452 (1986). 12. S. R. Bodner, Review of a unified elasto-viscoplastic
theory, Interim Scientific Report, prepared for the Air Force Office of Scientific Research, Boiling AFB, D.C. 20332 (October 1984). 13. S. Mukhejee and V. Kumar, Numerical analysis of time-dependent inelastic deformation in metallic media using the boundary-integral equation method. 1. Appl. Mech. 45,785-790 (1978). 14. R. P. Kanwal, Generalized Functions: Theory and Technique, pp. 75-104. Academic Press, New York (1983).
APPENDIX:
FINITE PARTS OF INTEGRALS
The “finite part” of an integral, denoted by pf, is defined as (see Ref. [14], pp. 75-104)
part as R + ~O,E-+ 0
’ /(x) k-divergent
1.
(Al)
A typical example is log ro,
1 I-
pf ;“f
=
n= 1
r:-”
I-n’
642) n> I.
For small r,, it then follows that (log
ro)1'(0,6)+o(ro), n =
Pf
d’(O,fI)+o(r~-“),
I n >
1
(A3)
where we have used that b’(r, 0) is bounded at r = 0, see eqn (43). It is convenient to rewrite eqn (58) in the form
(A4) By using eqn (A3), eqn (A4) can be rewritten as I, (O,&) = f2(OO)r0log
r0+roW0)+4r0)
(AS)
where f,(O,) has been defined in eqn (62) and
(A6) Equation (A5) shows that the order of I, as r. + 0 is defined by terms r. log r. and rr,.
1. D.
1052
ACHENHACH cl 01.
It is too complicated to compute all the terms in eqns (57) and (A6) by using the expansions given by eqns (52) and (53). Another way of computing the relevant terms follows by considering the following integral : 0
Pf
rVG dr = r,,
3
c 2M.+,(M,)+ ._,n+3
2 2N,+,G’JQ .__,n-1
n+I
1
-N&Wo
Jog ro
(A7)
where we have used eqns (52), (53). (Al) and (A2). We define 0:
PC&0,) = ; Equations
(64), (60) and (A6)-(A8)
~,(4!)=
Pf
s Cl
rVG dr+N,(O, 0,) log rO.
(A8)
then yield
P(0.0,)*~+‘(0,0) dO+C
(A9)
The terms with the constants c, and C?i,come from the first integral of (A6). The constant C?,,(c,) vanishes ford’ having Mode I (II) symmetry. The manipulations that are required to evaluate eqn (A8) by the use of eqn (30) are lengthy and complicated. The result is given by eqn (72). Now that the third term in expansion (63) is known explicitly, namely il* = i,(O)r log f
(AlO)
Ii, = f,(fI)r
(All)
and it is known that ti, is of the form
an alternative way of computing i,(O) is provided by returning to the substitution method outlined by cqns (41)(46). By collecting terms of order r WCfind A*[r+,(U)] = A*u,+V*o;.
(A12)
A direct calculation shows that rA*& is a function of 0 only. Hence the terms on the right-hand side of eqn (A12) are r- ’ times a function of 0. The standard Fredholm argument shows that for boundary conditions appropriate to Mode I, eqn (A12) does have a solution which is unique to within a term which would give rise to a uniform contribution to T, ,. The solution to eqn (A12) can now be determined by well-established methods. The result is again given by cqn (72) if the term H, and the last term which includes g6. which are for Mode II, arc omitted.
CRACK-TIP
0
FIELDS IN A VISCOPLASTIC
ooze-7683/87 s3.00+ .a0 1987 Pcrgamon Journals Ltd.
MATERIAL
J. D. ACHENBACH,N. NISHIMURA~and J. C. SUNG Department of Civil Engineering, Northwestern University, Evanston, IL 60201, U.S.A. (Received 26 May 1986 ; in revised form 13 October 1986)
Abstract-An integral representation for the particle velocity in terms of a Green’s function and certain linear combinations of the inelastic strain rates has been used in this paper, both for a numerical method to compute full field solutions and to develop unequivocal asymptotic expressions for the fields near a stationary crack tip. Specific rest&s have been obtained for a solid whose constitutive behavior is represented by the Rodner-Partom model. It is shown that the leading term of the near-tip particle velocity is of order r”*, and the higher-order terms are of the forms r log r and r. Expressions have been derived for the angular variations and for the multiplying timedependent intensity factors. The r log r term is absent for the Mode III case. Two questions have been addressed in further detail : the dependence of the intensity factors on time and the importance of the higher-order terms. The numerical results show a stress intensity factor which decays with time. At a small distance from the crack tip the numerically computed normal&d opening stress r22(xI.0, ~)/T~~(x,,O,O)has been compared with a one-term asymptotic representation, i.e. with &(f)/&(O). The two curves diverge at very small values of time. The inclusion of a second term in the asymptotic expression for the stress gives very good agreement over a larger time range.
1. INTRODUCTION
A clear understanding of the mathematical structure of the fields of stress and deformation near a stationary crack tip is an essential prerequisite to an investigation of conditions for the onset of crack propagation. In a viscoplastic solid the fields in the immediate vicinity of a crack tip can be very intricate. In this paper near-tip viscoplastic fields have been investigated in detail for a class of constitutive models for which the inelastic strain rate approaches a bounded value as the stresses increase beyond bounds. The Bodner-Partom model[l], which has been shown to describe the mechanical response of a wide class of materials, displays this particular behavior of the inelastic strain rate. Fields near a crack tip are usually analyzed asymptotically by assuming a near-tip field of the general form r’f(O). When this expression is substituted into the governing equations and into appropriate boundary conditions, and terms of equal order in r are collected, there results a non-linear eigenvalue problem, which yields both p and S(O). For viscoelastic solids this method was used by Riedel[2] and Riedel and Rice[3], who investigated the neartip fields of stationary cracks for Mode III and Mode I, respectively. They considered a viscoelastic constitutive equation of the power law type. For the Mode I case their asymptotic solution displays the well-known HRR singularity if the power is greater than one. When the power is less than one, the solution becomes square-root singular. Hui and Riedel[4] considered quasistatically moving cracks for both Mode I and Mode III conditions. When the power is greater than 3, their singularity is different from both the HRR and the square-root singular forms, while for a power less than 3 it reduces to a square-root singularity. Lo[5] followed the Hui-Riedel approach, but he took into account the effect of inertia. The implications of square-root singular crack-tip fields for viscoplastic solids have been explored by Hart[6] and Freund and Hutchinson[‘l]. For a Mode III quasi-static steady-state analysis, Hart[6] related the stress intensity factor at the tip to the corresponding far-field quantity. He also discussed transient cases on the basis of certain assumptions. Freund and Hutchinson[7] developed a similar relation between the tip and the far-field stress intensity factors for the dynamic steady-state Mode I case, using an essentially linear relation between the stress and the inelastic strain rate. In this paper a more fundamental approach towards the computation of near-tip fields is presented, in the sense that asymptotic forms are not assumed a priori. Instead, the t Present address : School of Civil Engineering, Kyoto University, Kyoto 606, Japan. 1035
1036
J.
D.
ACHl%MCli
cl
d.
‘5 0 Fig. 1. Crack-tip geometry. particle velocity field is expressed by a representation integral over the complete domain of the cracked body. This integral involves an elastic Green’s function and linear combinations of the inelastic strain rates. Since the inelastic strain rates are bounded at the crack tip, the integral can be expanded in terms of the distance to the crack tip, to yield an unequivocal expression for the near-tip field. Details have been worked out for a stationary crack tip. It is shown that the leading term of the particle velocity is of order r”‘, and the next terms are of the forms r log r and r. Expressions have been derived for the angular variations and for the multiplying time-dependent intensity factors. The r log r term is absent for the Mode III case. Two points are of interest : the dependence of the multiplying terms on time and the importance of the higher-order terms. The representation of the particle velocity in terms of an appropriate Green’s function and the rates of change of time of the inelastic strains has not only been used to obtain near-tip asymptotic expansions of relevant field variables, but also to develop numerical solutions by the use of the boundary integral equation method, in conjunction with an iteration technique. The numerical calculations have been carried out under the assumption of small-scale yielding. Thus, the equations of the Bodner-Partom model have been used in a region around the crack tip. The conditions on the boundary of this region are provided by the singular term of a corresponding elastic solution. The results show a stress intensity factor which decays with time. At a small distance from the tip the numerically computed normalized opening stress t22(~,,0, 1)/r22(~,r0,0) has been compared with a one-term asymptotic representation, i.e. with K,(r)/K,(O). The two curves diverge at very small values of time. The inclusion of a second term in the asymptotic expression for the stress gives very acceptable agreement over a longer period of time. As early as 1971 the boundary integral equation method (BIEM), also called the boundary element method (BEM), was applied by Swedlow and Cruse[8] for elastoplastic material behavior. In recent years the method has been used extensively for the numerical analyses of viscoelastic/plastic solids of various shapes including bodies containing cracks. Recent work was summarized by Mukherjee[9]. 2. FORMULATION
The states of deformation that will be considered are two-dimensional, either plane or antiplane strain. The body contains a straight crack. At one of the crack tips, say A ,, the applied loads give rise to elastic stress intensity factors KF(r), K:(t) and K:,(t), if the material were linearly elastic. We will, however, consider a material that displays inelastic behavior as the stress level increases. Hence there will be a domain of inelastic deformation near the crack tip. It is, however, assumed that this domain is small as compared to the crack length, and that linear elasticity applies outside the near-tip domain of inelasticity. To analyze the field right near the tip A i, we follow the arguments that have been introduced for small-scale yielding. Thus, we magnify the geometrical scale near A ,, so that the geometry becomes one of a semi-infinite crack in a full space, as shown in Fig. 1. A Cartesian coordinate system is centered at the crack tip. For the plane strain case, possible boundary conditions on the bounding surfaces of the body are replaced by the asymptotic condition
u --, g
[K~(t)f,(e)+K~(r)f,,(e)l
as r -+ co, where p is the shear modulus, and (r, 0) are polar coordinates centered at the tip
Crack-tip fields in a viscoplastic material
1037
in the magnified configuration. Note that the definition of the stress intensity factors differs from the conventional one by the factor (27r)1’2. The functions f, and f,, give the angular variations of the Cartesian components of the elastic near-tip displacements for Modes I and II, i.e. -
-
cos
cos
3012 4) 012 4 f,W = (K (K+ 4) sin 8/2 - 4 sin 38/2
fl,(Q = where
K
(K+ {) sin O/2+ 4 sin 30/2 (-K+;) cos 8/2-t cos 38/2
(2a)
(2b)
is a constant defined in terms of Lame’s constants 1 and p 1+3/A
(3)
K=F. We now proceed to the governing equation for the problem. The strain rate i =
“,[Vti+ (Vti)‘],
where (‘) = d/df
(4)
is usually decomposed into the elastic part Ccand the inelastic part d’ as c = $+i’.
(5)
The elastic part is related to the stress rate + by Hooke’s law +=I1
trdc+2@
(6)
where 1 is the unit tensor and “tr” is the 3-D trace. Since the crack is stationary, the effect of inertia is not considered. Hence, the equilibrium equation is
v-i = 0.
(7)
An equivalent form of eqn (7) in terms of the displacement rate ti follows from eqns (4)(7) as A.*i = V.,j’
(8)
A* = pV2+(A+p)VV-
(9)
where the operator A* is defined as
and for convenience we have introduced the quantity ti.‘,defined by t+’= 11 tr i’+2$.
(10)
Equation (8) is the governing equation for our problem. This equation, together with the boundary conditions i:=+*n=O
on
x2=0*,
x,x0
(11)
and an asymptotic condition of the form of eqn (1) with u and @iI replaced by ti and k:,,, defines ri. When k:,, = 0, the behavior of the elastic solution away from the crack tip becomes
J. D.
1038
ACHENBACH
Ii *
O(r-
al.
et
‘i2).
(12)
To complete the formulation, we supplement these equations with the regularity requirements at the crack tip ir - bounded as r J 0
(13)
+ - 0(1/r)
(14)
as
rl0.
In the present problem, we will also discuss the Mode III case, where T,~ = r33 = u, = 0 (tl, /I = 1,2). We prefer to use the stress function 4 because it makes the stress computation easier. The stress function defines the stress components r‘. : = 7.3
(a = 1,2)
(15)
as r. = e,fl Q#J
($4 : = WW
(16)
where em8is the permutation symbol, eaB= e,,,. The compatibility condition e,&& = 0
(E, : = &,j)
(17)
together with eqns (5), (6) and (16), yields the governing equation
The boundary condition on the crack faces and the regularity condition analogous to requirement (14) yield i = const
(19)
on the crack faces with the constants for both faces being the same. Also, on the crack line, we have 5, = a,4 = 0.
(20)
Finally, the asymptotic condition is 4 --) -2r”*k,E,,(t)
cos 46 as
r + co.
(21)
When X:,(t) = 0, condition (21) is replaced by d; - O(r-‘I*)
as
r + co.
(22)
In this case, the constant in eqn (19) is zero because d(O) = 0. This follows from m <*dx,=-
o= s0
Oq$4 dx, = &o>-&co> I0
where eqn (22) has been used. In the sequel, we will consider only the cases where
= d(O)
(23)
1039
Crack-tip fields in ;Lviscoplastic material kfi,,.,*,w
3. REPRESENTATION
(24)
= 0.
OF THE VELOCITY FIELD
The term on the right-hand side of the displacement-rate equation, eqn (8) is of the nature of a body-force distribution. This implies that the solution to cqn (8) can bc cxprcsscd by the use of a Green’s function. Representations in terms of Green’s functions are at the basis of the boundary integral equation method (BIEM), and they have therefore become quite well known. For rate-type inelastic analysis the appropriate representations were first developed by Swedlow and Cruse[8]. This section gives a brief derivation of such representations for the purpose of making the paper self-contained. For further detail we refer to Mukherjee[9] and the papers cited there. For an unbounded body containing a semi-infinite crack, the Green’s function C(x, y) is defined by A,+G(x,y) = -18(x-y)
(25)
where x = (xi, xJ, y = (y ,, yJ. The operator A* is defined by eqn (9), and the subscript indicates differentiations with respect to xi, i = 1,2. On the faces of the crack (xi < 0, x2 = O*) we have U+LY)
=
(26)
0
where T,G stands for the elastic tractions computed from G. We also have
W, y) -
Wag
1x1) as as
V,W,Y)
- WQ/lxl)
Ix1+
(27)
00
Ix]+ 0
(28)
as 1x1-+ 0.
(29)
A standard analysis based on Muskhelishvili’s method yields C(x, y) as
-K
log
(,/Z-&
($--$I+-~
log (&+&,)-log
)+&(z-Z”
+
(++&o)} +=)I] (30)
where
K
is defined by eqn (3), and
i-/(-l),
z=xl+ix2,
z. =y,+iy2,
5=x,-ix2,
f. =y,-iy,.
(31)
Equation (30) agrees with the expressions given by Erdogan[lO]. A suitable representation of the solution to eqn (8) is obtained by the use of Green’s identity
1040
J. D. ACHENBACH er a/.
Fig. 2. Domain
(A:v) * w dA, -
da,
for use of integral
*
(A:w) v dA, = s d&z
identity.
(T,v)* w dc(T,w)*v dc sad.?, s ddn,,
(32)
which holds for an arbitrary pair (v, w), and whcrc dQR,ris the domain bounded by circles of radii E and R (0 < E < R) centered at the tip, and by the crack faces, as shown in Fig. 2. We now substitute ii(x) for v(x) and G(x,y) for w(x) in eqn (32). Subsequent use of the divergence theorem and the regularity conditions (13), (14), (28) and (29), followed by letting R + co, E + 0, and using conditions (12) and (27), gives
it(x, t) =
s
V$(Y,
R2
4 *b'(y, t) dA,
(33)
or in component form
tij(X, t) =
VykG,j(YvX)~h(Yt 1)dAy.
(34)
Here V, = i(d/8y,)+j(d/&), and dA, = dy, dy2. In deriving these equations we have assumed that the integrals in eqn (33) or eqn (34) converge. For this to be true, it is sufficient to assume that 6’ satisfies (35)
lti’l - O(Pl)
as
r + co
(36)
with al > -1,
a2 < -1.
(37a, b)
Although inequality (37a) appears unnecessarily restrictive, it turns out that these inequalities guarantee the satisfaction of conditions (13) and (14) automatically. It is noted, however, that inequalities (37) are just sufficient conditions. For the Mode III case a similar procedure gives
d(X> 4= where
l2a,,
G(y, x)e,&~,
0 d4
(38)
Crack-tip fields in a viscoplastic matcrial
W,Y)
(& - &0) ($-
= - $ log
(Jz+J&l)
&I)
1041
(39)
(J~+JzlJ).
Hence we have r$(x, 1) = epqax,
a,,WY,
4. ASYMPTOTIC
x)e,d(y,
0 dA,.
(40)
EXPANSIONS
In the following sections, we employ the potential representation for i given in Section 3. The results of this section are analytical, while numerical results are presented in Section 7. For a certain class of constitutive models the quantity b’, given by eqn (lo), is finite. For example, a viscoplastic model proposed by Bodner and Partom[l] shows this property. In this section we will obtain the asymptotic expression for ti for the case that h,$i&(r, 0,.t) exists for all 8. In a conventional asymptotic analysis one seeks near-tip fields of the forms
where ~j, tii (j = 1,2 , . . .) are functions which are increasingly less singular near r = 0 with increasing subscript. Since near r = 0 ii =
d’(0,8, t)+o(l)
(43)
is bounded, and
we find by substituting expressions (41) and (42) into eqn (8) A*u, = 0 A*irz-V~t$,
(45) = 0.
(46)
The operator A*, which is defined by eqn (9), is here used in polar coordinates. Equation (45) implies that u, is of the same form as for the linearly elastic case. The appropriate separation of variables form of eqn (46) can only be guessed at this stage. One might try til =
tqe,t).
(47)
This expression for I(& r) does, however, not necessarily lead to a solution of eqn (46) because the corresponding homogeneous problem has non-trivial solutions. For the Mode I case consider, for example, the homogeneous equation A*[rf(e, t)] = 0 with boundary conditions
(48)
J. D. ACHKNBACH et (11.
1042
0 =
T[rl(e, f)] = 0
7[:
i, * T[rl(e, r)] = 0
0 = 0:
and
(49)
i, * [rr(6, t)] = 0
(jOa, b)
where the operator T has the same meaning as in eqn (26), but in polar coordinates, and i, are unit vectors. Equation (48) has a solution corresponding to a uniaxial tensile field in the xi-direction. Hence the approach based on eqn (47) becomes ambiguous at this point. An unequivocal form of the term it2 can, however, be obtained directly from the representation given by eqn (34). Using polar coordinates (r, 6) and (ro, 0,) for y and x, we split cqn (33) into two inlcgrals
ss n
Co, 00,O=
‘0
--I 0
rVG(r, 8; ro, 8,) -S’(r, 0, t) dr de
rVG(r, 19;ro, 0,) *ci.‘(r, 0, I) dr df?.
+
(51)
Here the gradient operator is taken with respect to the polar coordinates r and 8. Next VG is expanded into a power series of r. We can write r < r. :
Mnc I (4 00)
(52)
(n+ I)/2
r > r. :
rVG=
f n--l
0 y
N,+ I (6 eo>.
(53)
The coefficients in eqns (52) and (53) can be obtained from eqn (30). It is found that No is independent of do. Substitution of eqns (52) and (53) into eqn (51) yields after some rearrangements Co,
00)
=
f0 +tl(e0)rY2
+i,
(r0,0,)+i,
(ro, 0,)
(54)
where
ss n
i,
=
m
--II 0
N,(O) -d’(r, 0) dr de
(55)
r- ‘/‘N, (0,&J - ci’(r, 0) dr dt7
(56)
*tf’(r, 0) dr de
(57)
The dependence on time t, shown in eqn (51) is implied; but not explicitly shown in eqns (54)-(58) and in the sequel. In the next step the orders in r. of the integrals (57) and (58) are determined.
Crack-tip fields in a viscoplastic material
1043
Since ti’(r, 0) has an estimate given by eqn (43) it follows that eqn (57) is of order ro, i.e.
where
A@,) =
1 -No(e)-2N1(8, 80)
M,, #,&)
(60)
A discussion of the order of the integral in eqn (58) is given in the Appendix. The result is
L (ro,0,) = i,(e,)r, log To +roWo)+ 4ro)
(61)
where
t2(eo) = -
s’
N2(e, eopii(o,
e) de.
(62)
--I
The function B(O,) is defined by eqn (A6). By combining eqns (54), (59) and (61) it is concluded that ri(r,, 6,) is of the form
ii(ro,eO) = ~o+i,(~o)r~~2+~2(~0)r0 log ro+~3(Oo)ro+o(r0)
(63)
where
i,(e,)
=
A(e,)+B(e,).
(64)
The constant to in eqn (63) corresponds to a rigid-body motion ; it will be left out from further considerations. Some additional manipulations reduce i,(6,) and i,(O,) to the forms f,
(0,) = [&f,(0,)+ li;,f,,(e,w~
(65)
where f,(e,) and f,,(S,) are defined by eqns (2a) and (2b), and
r
- ‘/2D,,,I- d’(r, 8, t) dr de.
(67)
E,,r, * c+(O, 8, t) de.
w-0
Also
In eqns (67) and (68)
J. D. ACHENBACH EI ~1
1044
70 30 3 3 - - sm T + - sm 2 ‘2 2
D, = i 3
38
2
2
- -SK-
3
78
2
2
+ -sm-
1 38 3 78 K - - sin- - -sm2 222
( 1 1
-cos2222 E =
’
(
36
3 78 + -cos-
(i-K)C0S28-220s48
E,, = (1+x)
--
-
(I-K)
sin 48
-2
cos-
~(COS 2e-cos
4e)
30
3
.3e
21
sm-+-an
3.
~(COS 2e-cos (K+
1)
(69b)
22
-2 sin 48 cos 20+2 cos 40>
sin 2e+2 sin 48
(~-3)
3 70 + -cos 2 2 2
(69a)
sin 28-2
(70a)
48)
sin 48
G’Ob)
It is shown in the Appendix that i,(e) can be written in the form
where
P(e, e,) eciyo,e) =
4np(~ +Kpl(e0)(2~
w (e-O,)Hl(e)+H,(0)).l’(O,u)
+g,um2n w (e-O,)H,(e)+H,(u)).~‘(o,a) -t-g,ted (2~sgn(0- hdE, I (0))- ~i(o,0)
+h(ed (2n w (6-edE, (eww4 0)
+gs(eo)(E,(~).~i(O,~))+g,(~O)(El,(~)*~i(O,e))l
(72)
and
(73)
In eqn (72)
1
g,(e) = ___ 2(1+~)
cos e [ -sin e
e1 1' gm=[co: (74a-f)
Crack-tip
C
-KCOS 2u-cos
fields in a viscoplastic
40
HA@= -((h.+ I) sin 20)/2-sin
-((K+ 40
cos 20-cos
H4(e) _(K+ +(K+ i)o
2 cos 40 sin 40 (I -K) sin 20-2
1)
sin *O-sin 48
-(K-3)
2
1045
I) sin *Q/*--sin cos 2O+cos40
40
1
1
(I -K) sin 20+2 2 cos 48 sin 48
- cos 28 +
cos
48
1
1
48
-sin 28 + sin 48
--OS 2e+c0s 48
2
material
COS w-2
COS 48
2 sin (20 -sin 48)
*(sin 28 -sin 48)
-(K-t-i)
COS 28-e
COS 48
1 ’
(75a-d)
In the expression for c,,,,, we have denoted by pf the finite part of the integral defined in the Appendix. Equation (63) implies that the dominant singularity of the strain rate has exactly the same form as in linear elasticity. The rates of the stress intensity factors &;.,i are given by eqn (67). The second singular term is of the order of log r near the tip. The magnitude of the multiplying term is determined only by the distribution of di at the tip. Because of this last property, expression (68) holds true for any crack configuration, whereas eqn (67) is true only in the sense of the “small-scale yielding”. The logarithmic singularity in the strain rate does not give rise to a logarithmic stress singularity except in T, , for Mode I. For the Mode III case, we can use the same general approach to show that
(76)
where I” K,, =
4x s n
d0
R’ S[0
-ri’,(r,
30
30
0) sin -* + h\(r,O) cos y
1
r- I” dr
(77)
(‘0
li(0,) =
OS ;
0
[C’*$(O,O) cos 0-C.C?(O,O) -
sin 81 d0
psi(o,e) sin e+e+(o,e)
cos e] de-[e7ii(o,e,)]e
(78)
where
(7% b) Note that the stress rate and hence the strain rate do not have logarithmic singularities. As for the Mode I case, we can relate the absence of the logarithmic term either to the nonexistence of a non-trivial solution to a certain boundary value problem or to the lack of a certain term in the far field expansion of the Green’s function.
J. D. ACHENBACH
1046
er al.
In the published literature, the choice of terms for asymptotic expansions is generally made on a trial and error basis. The analysis presented in this section gives a more rational foundation to the investigation of the asymptotic near-tip structure.
5. THE BODNER-PARTOM
MODEL
In Ref. [l], Bodner and Partom proposed a set of constitutive equations for an elastic-viscoplastic material. These equations have the convenient property that no separate specification of a yield criterion is required, nor is it necessary to consider loading and unloading separately. Within the context of the Bodner-Partom model both elastic and inelastic deformations are present at all stages of loading and unloading, but the inelastic deformations are very small at low stress levels. The Bodner-Partom equations have been used to model the mechanical behavior of several metals over a wide range of temperatures and strain rates. Very satisfactory agreement between predicted and experimental results has generally been obtained. For a recent discussion of the Bodner-Partom model we refer to Ref. [ 121. The constitutive equations follow by supplementing eqns (5) and (6) with an expression for the rate of inelastic strain. In a recent version of the Bodner-Partom model this strain rate is given by, see Ref. [ 121
ii =
DoL exp [ - (Z2/3J2)“/2]
JJ~
(80)
where T’ is the stress deviator f’ =
tr
f-f1
f
(81)
and
Also Z = Z, +(Z,-Z,)
eemWp.
(83)
In eqns (80) and (83), Zo, Z,, m, n and Do are material constants, while W, is the work over the inelastic strains
W, =
J
t-8
dt.
(84)
For a uniaxial tension of magnitude cr, it is simple to plot the inelastic strain rate vs a/Z. For various values of n the curves are shown in Fig. 3. The linear plot of Fig. 3 shows that d’ remains very small at small stresses with a rapid change over a small range of a/Z (the “yield” region), and that an asymptotic limit is approached as Q becomes large. It is of interest to note that the components of ii, as defined by eqn (80), are always bounded. Even if r is singular (which is the case at a crack tip), the singularity of r as r + 0 will be cancelled by the singularity of ,/J2 in the term r’/ JJ2. The term exp [ - (z2/3J2)“/2] is bounded whether J2 is singular or not. Hence eqn (80) implies that the inelastic strain rate is bounded at a crack tip. It then follows from eqn (10) that the components of rY are also bounded at a crack tip, and it can thus be concluded that the asymptotic results of Section 4 apply to the Bodner-Partom model.
1047
Crack-tip fields in a viscoplastic material
UNIAXIAL
CASE
01 0
.2
.4
.6
.8
1.0
1.2
1.4
1.6
1.8
2.0
Fig. 3. Inelastic strain rates for uniaxial tension, u, and various values of II.
6. NUMERICAL
PROCEDURE
A numerical procedure has been developed to compute the stress intensity factor and the opening stress r 22 at a position on the crack line. As pointed out by Mukherjee and Kumar[l3], a numerical evaluation based on eqn (33) is very efficient for viscoplastic problems. Here we use this approach (generally referred to as the BIE method) for cracktip fields governed by the Bodner-Partom equations. The procedure is as follows. (1) Compute the initial field. Since the inelastic strain is zero at t = 0, the initial field is simply the linearly elastic one. In the present context this field is
Hence, the solution is known for t = nAt with n = 0, where At is the time increment. (2) Subdivide the near-tip domain into cells, and compute d’ at the centroid of each cell by using eqn (80). Then assuming that ii may be taken as constant over each cell, compute li by using (33) and WIby differentiating eqn (33) (3) Compute + from eqns (5) and (6). Update r by using r[(n+ l)At] = r(nAt)+At+(nAt).
(86)
(4) Use the updated value oft to compute Pi as in (2), and repeat (2) and (3). When the edges of the cells are straight, the integrals involved in eqn (33) can be evaluated analytically, as discussed by Cruse and Polch[ 111. By picking up the square-root singular terms in these integrals and summing up their coefficients, k,;.,, can be computed. Alternatively, the computation of g ,,,! may be carried out by using eqn (67) with the piecewise constant approximation for i’.
7. NUMERICAL
RESULTS
Numerical computations have been carried out for the non-hardening case of the Bodner-Partom model. With reference to eqns (80) and (83), this case corresponds to
1048
J. D.
ACHENBACH
er
al.
Fig. 4. Cell configuration [or application of boundary integral equation method
zo/z, =
1.
(87)
It is then convenient to introduce non-dimensional quantities defined as u = u(WJ3)/(W> f = X(Z2/3)/(K:)r,
T = r/V/,/3) f = Wo)WJ3)
(gga, b) Wa, b)
(90) Here the dot notation represents the derivative with respect to f. In terms of these dimensionless quantities, eqns (80), (6) and (8) convert into 8’ = [T’/($T’*T’)“2]
exp [-i(iT’*T’)-“1
(93) where 0 and v2 are the de1 operator and the Laplacian with respect to the dimensionless spatial coordinates fi, and T’ is the dimensionless stress deviator. Since n/p and (n+p)>/p can be expressed in terms of the elastic Poisson’s ratio, v, it now follows that numerical values have to be chosen only for v and n. Here we choose n = 1 and
v = 0.2.
(94)
For the numerical procedure based on the boundary integral equation method described in Section 6, the cell arrangement shown in Fig. 4 has been used. The expression for the rate of change of time of the stress intensity factor is given by eqn (67). Once the inelastic strain rate has been computed, the term tij can be determined and eqn (67) can be evaluated numerically. A subsequent numerical integration with respect to time produces
Crack-tip
0
Fig. 5. Normalized
2
4
fields in a viscoplastic
6
8
10
1049
material
12
opening stress at f = 0.0014 and ratio E,(r) dimensionless time, I.
14
76
= K,(f)/K,(O).
78
20
both as function
of
&(Q Figure 5 shows Ki(Q/Xr vs the dimensionless time r. It is noted that Ki@)/il~ decreases with time. Also shown in Fig. 5 is the normalized stress T&,, 0, T)/T&i, O,O), evaluated at R , = 0.0014. The two curves diverge rapidly as the dimensionless time increases. It is instructive to compute the actual physical distance to the crack tip which would correspond to di = 0.0014. Suppose we consider a crack in an unbounded medium under far-field uniform tension 222 = co. Then Kf = ~*~(~/2). For titanium, Ref. f 1) suggests the following value for the yield stress, by N 250 N mm-*, while Z - 1200 N mm-*. Let cro N ~rv/2. Equation (89) then yields xl N 1.63 x IO- * ~2,. Hence 2, = 0.0014 corresponds to x, N 2.279 x lo- 5u, which is a very short distance to the crack tip. The two curves in Fig. 5 thus imply that even very close to the crack tip, the singular term by itself would give a very poor approximation to the near-tip stress T**, with an error which increases with time f. The relation between physical time f and dimensionless time fcan be obtained from eqn (89) as I = 0.635 x IO?, where the following values from Ref. [I] have been used DON lo4 s- ’ and p= 0.44 x lo5 N mm-*. As predicted by the constitutive model used in this paper, the relaxation of the stress intensity factor takes place very rapidly. An improved asymptotic representation of the near-tip stresses can be obtained by including higher-order terms. The calculation of these terms requires the evaluation of the integrals for i,, eqn (68), and f,(t),), eqn (71). The r log r term in eqn (63) contributes only a constant term to fz2. A constant contribution is also obtained from the term i,(O,)ro. After the calculations have been completed we find on the crack line (0 = 0)
T 22
-
4r,
- 3 376
Js
-
or
(96) This two-term expansion gives excellent agreement with the numerical result up to I= 4, as can be seen from the numerical values listed in Table 1. The agreement deteriorates rapidly for f > 4. Figure 6 shows ~~~,(~)~~~~vs the dimensionless time f. Also shown is the ratio T2(Z1, 0, f)/T2(Z1, 0,O) evaluated at 2, = 0.0047. As for the Mode I case, the two curves diverge rapidly as f increases. It is relatively simple to write out explicitly the asymptotic results for the Mode III case. At the crack tip we obtain from eqn (91)
J. D.
1050
ACHESBACH er al.
Table I, Comparison of asymptotic rcprcscnfatlons with numerical results calculated by the Mode 1 722WJ3)
Mode III r,/(Gy’3)
WD, Z/J3
Two terms
I
13.04
1 L
II.05
BIEM
0
I
0
2
4
Two terms
Error (%)
13.10 II.17 9.12 7.12 6.73 6.35 5.97 5.60 5.235 4.87 4.52 4.17 3.84 3.504
8.93 6.84 6.43 6.02 5.62 5.21 4.8 1 4.42 4.03 3.63 3.25 2.87
3 4 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0
BIEM
22.18 18.78 14.47 10.08 9. I94 8.32 7.444 6.555 5.68
3.93 4.5 5.2 5.9 6.96 8.0 9.24 10.84 12.95 15.36 18.0
BIEM
Error (%)
22.28 18.82 14.35 10.42 9.745 9.104 8.502 7.94 7.41
3.26 5.65 8.6 12.4 17.4 23.3
I
6
8
10
12
14
16
18
20
E Fig. 6. Normalized anti-plane shear stress at f = 0.0047 and ratio R,,,(r) = K,,,(I)/K,,,(O), both as functions of dimensionless time, I.
Ei(o,q
=
(-sin:e4’). cos
(97)
Use of this result in eqns (76)-(78), and subsequent integration with respect to time yields
The values of T2 at 2, = 0.0047, Z2 = 0 computed according to eqn (98) and by the BIEM have been listed in Table 1. Good agreement is observed up to I = 4. The results presented in this section show that the square root singular term presents a very poor approximation to the crack-tip field even at very small distances from the crack tip. The approximation becomes worse as time increases. The use of additional terms in the asymptotic expression for the particle velocity increases the length of time for agreement between the asymptotic representation and the numerically computed stress component at a crack-line point just ahead of the crack tip. Acknowledgemenr-This work was carried out in the course of research sponsored by the U.S. Of&x of Naval Research (Contract No. NOO014-76-C-0063).
1051
Crack-tip fields in a viscoplastic material REFERENCES
1. S. R. Bodner and Y. Partom, Constitutive equations for elastic-viscoplastic strain-hardening materials. J. Appl. Mech. 42,385-389 (1975). 2. H. Riedel, Cracks loaded in anti-plane shear under creep conditions. Z. Mefall. 69(12), 755-760 (1978). 3. H. Riedel and J. R. Rice, Tensile cracks in creeping solids, Frucfure Mechanics (Praceedings oj the 12th National Symposium on Fracture Mechanics, St. Lous. May 1979) (Edited by P. C. Paris), ASTM SIT 700, pp. 112-130. American Society for Testing and Materials, Philadelphia (1979). 4. C. Y. Hui and H. Riedel, The asymptotic stress and strain field near the tip of a growing crack under creep conditions. Int. /. Fracture 17(4), 409425 (1981). 5. K. K. Lo, Dynamic crack-tip fields in rate-sensitive solids. J. Mech. Phys. Solidc 31(4), 287-305 (1983). 6. E. W. Hart, A theory for stable crack extension rates in ductile materials. Inf. /. Solidr Sfrucfures 16, 807823 (1980). 7. L. B. Freund and J. W. Hutchinson, High strain-rate crack growth in ratedependent plastic solids. /. Mech. Phys. Solids 33(2), 169-191 (1985). 8. J. L. Swedlow and T. A. Cruse, Formulation of boundary-integral equations for three-dimensional elastoplastic flow. Inr. J. Solids Strucfures 7, 1673-1683 (1971). 9. S. Mukhejee, Boundary Element Merhoa!s in Creep and Fracture. Applied Science, London (1982). 10. F. Erdogan, On the stress distribution in plates with collinear cuts under arbitrary loads, Proceedings o/lfh U.S. National Congress of Applied Mechanics (Edited by R. M. Rosenberg), pp. 547-553, American Society of Mechanical Engineers, New York (1962). I 1. T. A. Cruse and E. Z. Polch, Elastoplastic BIE analysis of cracked plates and related problems, Inl. 1. Numer. Me/h. Engng X3,429-452 (1986). 12. S. R. Bodner, Review of a unified elasto-viscoplastic
theory, Interim Scientific Report, prepared for the Air Force Office of Scientific Research, Boiling AFB, D.C. 20332 (October 1984). 13. S. Mukhejee and V. Kumar, Numerical analysis of time-dependent inelastic deformation in metallic media using the boundary-integral equation method. 1. Appl. Mech. 45,785-790 (1978). 14. R. P. Kanwal, Generalized Functions: Theory and Technique, pp. 75-104. Academic Press, New York (1983).
APPENDIX:
FINITE PARTS OF INTEGRALS
The “finite part” of an integral, denoted by pf, is defined as (see Ref. [14], pp. 75-104)
part as R + ~O,E-+ 0
’ /(x) k-divergent
1.
(Al)
A typical example is log ro,
1 I-
pf ;“f
=
n= 1
r:-”
I-n’
642) n> I.
For small r,, it then follows that (log
ro)1'(0,6)+o(ro), n =
Pf
d’(O,fI)+o(r~-“),
I n >
1
(A3)
where we have used that b’(r, 0) is bounded at r = 0, see eqn (43). It is convenient to rewrite eqn (58) in the form
(A4) By using eqn (A3), eqn (A4) can be rewritten as I, (O,&) = f2(OO)r0log
r0+roW0)+4r0)
(AS)
where f,(O,) has been defined in eqn (62) and
(A6) Equation (A5) shows that the order of I, as r. + 0 is defined by terms r. log r. and rr,.
1. D.
1052
ACHENHACH cl 01.
It is too complicated to compute all the terms in eqns (57) and (A6) by using the expansions given by eqns (52) and (53). Another way of computing the relevant terms follows by considering the following integral : 0
Pf
rVG dr = r,,
3
c 2M.+,(M,)+ ._,n+3
2 2N,+,G’JQ .__,n-1
n+I
1
-N&Wo
Jog ro
(A7)
where we have used eqns (52), (53). (Al) and (A2). We define 0:
PC&0,) = ; Equations
(64), (60) and (A6)-(A8)
~,(4!)=
Pf
s Cl
rVG dr+N,(O, 0,) log rO.
(A8)
then yield
P(0.0,)*~+‘(0,0) dO+C
(A9)
The terms with the constants c, and C?i,come from the first integral of (A6). The constant C?,,(c,) vanishes ford’ having Mode I (II) symmetry. The manipulations that are required to evaluate eqn (A8) by the use of eqn (30) are lengthy and complicated. The result is given by eqn (72). Now that the third term in expansion (63) is known explicitly, namely il* = i,(O)r log f
(AlO)
Ii, = f,(fI)r
(All)
and it is known that ti, is of the form
an alternative way of computing i,(O) is provided by returning to the substitution method outlined by cqns (41)(46). By collecting terms of order r WCfind A*[r+,(U)] = A*u,+V*o;.
(A12)
A direct calculation shows that rA*& is a function of 0 only. Hence the terms on the right-hand side of eqn (A12) are r- ’ times a function of 0. The standard Fredholm argument shows that for boundary conditions appropriate to Mode I, eqn (A12) does have a solution which is unique to within a term which would give rise to a uniform contribution to T, ,. The solution to eqn (A12) can now be determined by well-established methods. The result is again given by cqn (72) if the term H, and the last term which includes g6. which are for Mode II, arc omitted.