Continuum interpolation of stress fields in brittle fracture

Theoretical and Applied Fracture Mechanics 11 (1989) 225-234

225

C O N T I N U U M I N T E R P O L A T I O N OF S T R E S S F I E L D S IN B R r I T L E F R A C T U R E N. A R I Department of Engineering Science, Lafayette College, Easton, Pennsylvania 18042, USA

The locally nonlinear behavior of crack stresses is first represented by adjustable boundary conditions near the crack tip. The actual boundary conditions are then determined by continuity conditions imposed on a general parameterized solution. The thus modified formulation yields finite crack tip stresses consistent with atomistic and other hybrid models of brittle fracture.

1. Introduction Distinct modes of microscopic failure mechanisms preclude a unified treatment of fracture phenomena. Specific models are needed to explain diverse fracture morphologies and the rupture phases of crack nucleation and final separation [1]. Physical geometry of a crack tip is highly irregular. Nevertheless continuum mechanics models a smooth, sharp crack. This formulation focuses on one essential common thread of fracture processes - the steep increase in stress concentrations near defect structures and crack tips. As such it provides a first order explanation of how the crack geometry interacts with applied loads and causes stress intensification. The idealized crack geometry however screens out the effects of the discrete substratum and prevents continuum models from incorporating microstructural parameters directly. Further the sharp edge assumption leads to infinite crack-tip stresses. Hence the maximum stress hypothesis becomes unusable and it is to be substituted by other failure criteria. In order to avoid physically unacceptable singularities one has to either resort to atomistic models with nonlinear interactions or modify the continuum formulation to include microstructural elements. Atomistic models yield finite crack stresses and provide realistic simulations of bond rupture, slip formation and other microscopic failure behavior. The computational resources they require however are substantial. In adtlition detailed information about interatomic forces is difficult to obtain [2, Chapter 7]. The less expensive alternative - modification of continuum results - is done either by excluding a core region near the tip from linear elastic analysis or by importing elements from solid-state physics (nonlocal theories) which in effect permit stress averaging of discrete fields, [e.g. 3-5]. Difficulties related to the selection of a core radius (not a material parameter) however remain in force. The motivation for the present paper is also a reconciliation of the continuum picture with the discrete microstructure. Implicit discontinuities associated with the classical mixed boundary conditions of a slit geometry are simply not compatible with a continuum interpolation of atomic displacements. Thus we propose to modify the boundary conditions of classical elasticity so as to reflect the discreteness of bond rupture. Our strategy is then (1) to acknowledge the absence of a crack tip boundary in the sense of the classical elasticity, (2) consequently to leave the displacement boundary conditions undetermined within an atomic distance of the tip, and (3) finally to determine a regularized finite stress solution by means of continuity conditions and physical considerations. This approach recasts the crack-tip problem as an interpolation task in hybrid (discrete/continuous) displacement space and requires the investigation of the sensitivity of the results with respect to 0167-8442/89/$3.50 © 1989, Elsevier Science Publishers B.V.

226

N. A ri / Continuum interpolation of stress fields in brittle fracture

interpolation parameters. To illustrate the method we will focus on brittle fracture with a cusp like crack tip opening and utilize the Griffith crack as an example. To this end first in Section 2 we review the classical mixed boundary conditions and contrast them with the modified ones. A parametric solution to the generalized problem is given in Section 3. In the final section we discuss the effects of the adjustable boundary conditions on the interpolated deformation fields. The current method is similar in spirit to Barenblatt cohesive zone model. It differs from his formulation, however, by not introducing any elements extraneous to the theory of elasticity. The present approach yields a simple and natural interpolation of finite crack tip stresses that compare well with other atomistic a n d / o r hybrid models of brittle fracture.

2. Modified boundary conditions for a Griffith crack An isotropic, homogeneous, linearly elastic plate in the (x, y ) plane contains a traction-free line crack of length 2l in its center along the x-axis. The plate is subject to uniform tension 7~y = % at y = __+o¢, (Fig. 1). The classical elasticity treatment of the thus defined Griffith crack is first to solve the following boundary value problem

(x + ~ ) ( , , , ~ + V,xy) + ~(,,,xx + u,,,y) = o,

(1)

(x + ~)(u,x~ + v,~y) + ~,(v,xx + v,~,) = 0,

(2)

ox~(x, o) = ~ ( u , y + v,x) = o,

(3)

% ( x , o) = (x + 2~,) ~,,, + x~,~ = -Oo

for Ixl < l ;

v ( x , o) = o

for Ixl > l ;

u, v--*0

(4) (s) (6)

as(x 2+y2) 1/2~,

and then to add a uniform stress field Oyy = 0 0 to the solution obtained. In (1-6), X and /z are the Lame constants, and u and v are the displacements along the x and y-axes, respectively. Here and throughout, an index following a comma represents a gradient, e.g. u, x = ~u/ax.

IlL

~-

-

r

X

E

rFig. 1. Griffith crack.

-

-12

Fig. 2. Atomistic model of a crack.

N. A ri / Continuum interpolation of stress fields in brittle fracture

227

The general solution to (1-2) is given by

u(x, y) =

(2rr ) - , / 2 f_~( i / k ) [

IklA(k) + ( I k l y - (X + 3tt)/(X + t Q ) n ( k ) ]

dk

×exp(-klyl-ikx)

v(x, y)=(2rr)-l/2f~_¢ [A(k)+yB(k)] e x p ( - k l y l - i k x ) dk,

(7) (8)

where k is the transform variable for the Fourier transform along the x-axis, i.e. (9)

~ ( k , y ) = ( 2 ~ r ) ' / 2 f _ ~ g ( x , y ) exp(ikx) dx.

A(k) and B(k)

are determined by imposing the boundary conditions (3-5). First (3) yields

(10)

B(k)-- [(X +tt)/(X + 2/*)] IklA(k), and then the ensuing pair of integral equations

f0

o,y(x, 0) = - A 0 ~°kA(k) cos kx d k = - Oo for I xl < l;

(11)

v(x, O) = (2/¢r)'/2fo~A(k)

(12)

cos

kx

dk = 0

for Ix >~ l,

are solved to obtain [6]

A ( k ) =/2ooAolJ,(t:t)/(/d ),

(13)

A o = (2/¢r)l/Z2tt (~, + tt)/(}, + 2/*).

(14)

where

Substitution of (13) into (7) and (8) leads to infinite crack tip strains and stresses. This result is clearly not consistent with a continuum picture of a crack and linear elasticity. A closer look at the classical boundary conditions imposed on a lattice model of the crack however reveals the source of the incompatibility (Fig. 2). From the figure, first we observe that the location of the crack tip is not determined precisely. It is somewhere between the last broken and the first unbroken bond. Secondly, if we assume that the crack tip is at B, then the displacement condition v = 0 requires that the corresponding "tip atom" does not move even when it is subjected to a force. Atomistic theories avoid this incompatibility by prescribing forces at A and displacements at B. In classical continuum theories however the distance measure between interior and exterior crack tip points (i.e., A and B) becomes zero. In order to restore the compatibility, in the present work we introduce a buffer zone (region III in Fig. 2) where the boundary conditions are unspecified. In the outlying zones I and II we retain classical boundary conditions. This flexibility modifies the boundary conditions (4-5) into

o~,,(x, 0) = 0 Vx,

Oyy(X, 0 )

=

--O"0 for x < l

13(x' O) = ( ~ (¢2- x2)q

(I) forf°rl1l + a = c,

(15) (II)(III)

where a is the atomic distance along the crack line and b, q are the interpolation parameters to be determined by the analysis. The form of the interpolating curve can be chosen from a large class of mathematically convenient curves. Here the particular choice (c 2 - x2) q satisfies at x -- c the continuity conditions for the displacement and displacement gradient identically for q > 1, (Cusp-like behavior). The parameter q is related to the bending force constants of the "lateral" bonds. In general higher q values

228

N. Ari / Continuum interpolation of stress fields in brittle fracture

correspond to greater shear rigidity and somewhat narrower crack tip profiles. In the next section we present a generalized solution of (15) for q > 1.

3. General solution

The modified boundary conditions (15) again lead to a standard pair of dual integral equations f0~KA(K)cosKXdK=F(X)=I fo°~A(K) cos

KXdK=G(X)

for I X I < 1 , for

(16)

IXl > 1,

(17)

where

A(K) = (Ao/t%)A(k),

G(X)={O(C2-X2)q

for for

0

X = x/l,

K = kl,

X<~C, X > C,

C = c/l,

B = hA 0(w/2)1/2OO 112q-1.

(18)

The general form of the solution is ([6], p. 86) A(K)=(2/~')

dttJo(tK)f0 dz(t 2-z2)-l/2F(z) - fl d t J o ( t K ) ~ f

dzzG(z) .

(19) Substitution of (19) into (8) yields for the displacements along the crack line (detailed derivation is given in the Appendix, eq. (A10)) v1(x,O)=Do(12-x2) +b(c2-x2)

Do

00(~. -I- 2 ~ ) 2#(2~ + #) '

[

2B F(q+ I) ] rr1/2 i.,(q+ l/2) (C2_l)q-1/2F(1, ½_q; 3; y)

1/2 1 q

for px I < l ;

(20)

l 2 -- X 2

Y=

C 2 --

X 2 ;

U I I I ( X , O ) = b ( c 2 - - X 2 ) q for

l < Ix[ < c ;

o n ( x , 0) = 0

Ix] >~c;

for

where F(x) and F(a, r, y, z) denote the G a m m a function and the Gauss' Hypergeometric function respectively. From (20) we observe that the continuity requirements on v(x, 0) at x = l and x = l + a = c are satisfied. The additional continuity constraints we have imposed in Section 2 however now restrict the values of interpolation parameters. Specifically the continuity of 3v/3x(x, O) at x = c and x = l demands q > 1,

(21) + 1).

B = ('1rl/2/2)(C2 - 1)'/2-q_F(q + 1/2)/r(q

(22)

With B and consequently b determined, the displacements and stresses depend now only on one interpolation parameter q > 1. The final expressions are

Ol(X,O)=Oo(12-x2)l/2[1-g(1, ½-q;

3;

-(12-x2)/(ca-x2)] q-b(c2-x2) q

(23)

N. Ari / Continuum interpolation of stress fields in brittle fracture

229

0. I

5 > 0.05

0 0.9985

0.999

0.9995

I

1.0005

1.061

x/l Fig. 3. Crack tip opening ( a f t ~ 0.0001).

and from (A15, A16)

Oyy(X,0)/%-- -1 +

2)1/2F( , ½; q+ X

---1+

½" (c 2- 12)/(x 2- 12)) for Ixl >c;

x(xZ-12)'/2 ( 3 ~cc2~/-~

3 (x2-12))

F 1, ~ - q ; ~; ( c2 12)

forl
3

f 0

(24)

__

I

I

0.5

1.0

1.5

2.0

(x- l)/~ x 10.2 Fig, 4. Crack tip stresses ( a f t = 0.001), P ( x ) = ( a y y ( X ,

0)//o0 4-1)(C a - 1) 1/2.

(25)

N. Ari / Continuum interpolation of stressfields in brittlefracture

230

In Figs. 3 a n d 4, we plot the crack tip displacements a n d stresses. In the following section we discuss the variation in the results with respect to the b e n d i n g p a r a m e t e r q a n d how the stress profiles relate to atomistic results.

4. Discussion F r o m Fig. 4, we observe that the m a x i m u m stress occurs between the last b r o k e n b o n d a n d the first u n b r o k e n one, i.e. l < Xmax < l + a. The n u m e r i c a l values can be s u m m a r i z e d in the form Omax//O0 =

o~(q)(C

2 -

C = (l+ a)/l,

1) -1/2,

(26)

a( q )( l / 2 a ) l/z, where a(q) varies slightly with q for 1.5 < q < 3.5, (see T a b l e 1 a n d Fig. 5). The n o r m a l stress at the tip x = c = 1 + a is given b y (from (24) or (25)) c

a y y ( X = c , O ) = o 0( c 2 - / 2 ) ' / 2

(q - 1 / 2 ) (q-l)

(27)

a n d for large q it reaches rapidly its asymptotic value

oyy ( x = c, O) ~- a o ( l / Z a ) 1/2.

(28)

q ~

Both stress expressions (26) a n d (27) are similar in form to atornistic results with their inverse square root d e p e n d e n c e o n a characteristic length, i.e.

O/Oo = ,8( q ) ( 2 l / a )1/2,

(29)

where j3(q) takes into a c c o u n t the relative b e n d i n g strengths of the bonds. We should also note the role of

Table 1 Maximum stress values as a function of q

l)/l a

q

(Xm~

1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90

0.001000 0.001000 0.001000 0.001000 0.000960 0.000880 0.000820 0.000760 0.000700 0.000640 0.000600 0.000560 0.000540 0.000500 0.000480 0.000440 0.000420 0.000400

--

a Location of the maximum stress (see Fig. 4).

a(q) 3.047688 2.561157 2.230955 2.002000 1.877427 1.824577 1.802222 1.796357 1.800723 1.811398 1.826522 1.844565 1.864599 1.886561 1.909264 1.932750 1.957034 1.981602

N. Ari / Continuum interpolation of stress fields in brittle fracture

231

Table 2 q

A a

1.10 2.20 3.30 4.40 5.50

0.002333 0.00254 0.00263 0.00268 0.00271

a A r e a u n d e r the stress curve for / ~< x ~< 1 + 2 a .

a core radius is in a natural way assumed by the lattice parameter. Further combining (29) with the Griffith relation Oo=[2Ey//crl(1-1,2)]

1/2,

(30)

where E is the modulus of elasticity, I, is the Poisson's ratio and y is the surface energy we obtain a cohesive strength estimate ire = f l ( q ) [ 4 / r r (1 - p2)] 1/2( Ey/a)'/2

(31)

o~a =

(32)

or

material constant = M.

These results are consistent with other hybrid models of highly brittle fracture (e.g. Orowan [7], Barenblatt [8]). They are also equivalent to the results of nonlocal continuum models (Eringen and his co-workers [5,9]) where the interatomic interactions are taken into account by modifying constitutive and field equations. The advantage of the present approach is that through a minimum modification of the boundary conditions and without introducing assumptions external to linear elasticity theory we can obtain comparable realistic results. The present model has also a proxy parameter for bond bending and hence distinguishes between different types of bonds. A comparison of Fig. 3 with crack profiles obtained from atomistic simulations [10] suggests a relatively small q range 1.5-5. Within this range, however, a(q) is not much sensitive to changes in q. It is interesting to note that the area under the stress curves in Fig. 3 does not vary rapidly with q and is almost a constant (see Table 2). This integral corresponds to work done by tip stresses per unit area and is related to the surface energy of crack opening. Indeed, a rough estimate of et(q) and consequently q can be obtained (31, 29 and 26), as can be seen from Table 3. The a ( q ) values do not differ much for diamond and silicon even though diamond has considerably greater shear rigidity. The corresponding q values can be interpolated from Table 1 yielding 2.4 and 2.6 for silicon and diamond respectively. With a ( q ) thus estimated, (26) or (27) can now be utilized in formulating a fracture criterion. It can be either based on m a x i m u m stress or on a threshold stress acting over a critical distance.

Table 3 E s t i m a t e s for a(q) Material

E a) (101° N m - 2 )

a a) (10 - 9 m)

oc b) (101° N m - 2 )

Diamond Silicon

100

0.15

13

0.24

a [1, p. 77], b [2, p. 5], c [2, A p p e n d i x C], d [10, p. 89].

p c)

y a) (Nm -1)

a(q)

20.5

0.1

3.2

0.28

5.35 1.46

1.91 1.86

N, Ari / Continuuminterpolationof stressfieldsin brittlefracture

232

Appendix In this appendix we provide a detailed evaluation of the normal displacement and stress expressions. To this end we first simplify the solution to the dual integral equations (i.e. A(K) in eq. (19) and then substitute it into relevant formulas. Since F(z) = 1 and G(z) vanishes for z > C, eq. (19) reduces to

A( K) = J1( K ) / K -

(2/~r)B f c dt Jo(tK)H(t), Jl

(A1)

where d Ic

H(t) = - ~

Qt(C 2 t2

(C 2 - z 2) dz ( z2 _ t 2)1/2

_

(A2)

)q-l/2

O = (q + 1/2)r(1/2)r(q + 1)/r(q + 3/2) and F(x) is the Gamma function. Next we combine eq. (8) and eq. (18) to compute 09

v(X,O)=Dolfo A ( K ) c o s K X d K , (A3)

D o = o0(h + 2/z)/2/~(X + 2/*). Upon substituting (A1) into (A3) and using the integration formula [13; 6.671.8]

fo°~Jo(Kt) cos K X d K = { 0 ( t 2 - X 2 ) - 1 / 2

for for

X < t, X > t,

(A4)

we obtain for [ X I < C, (A5)

v( X, O) = Dol[ L( X) + (2/¢r)BQK( X)], where cos KX= ( 0(1 - X2) -1/2

L ( X ) = fo °° d K ( . l l ( K ) / K )

for for

[ X[ < 1, [ X I > 1,

K ( X ) = f;:xO,x) dt (C2-t2)q-l/2(t 2 - S 2 ) -1/2. For

(A6) (A7)

I Xl < 1 (A7) becomes [13; 3.196.1] [ ( C 2 - 1 ) - s ] q-1/2 K ( X)

(X/2)JoC2-1f ds

=

[s + (1 - X2)] 1/2

(C2 -- 1)q+1/2

F(1, 1/2" 3/2 + q; - Z ) ,

(A8)

2(q + 1)(1 - X2) '/z s = t 2 - 1,

Z = (C 2 - 1 ) / ( 1 - X2),

where F(a, /3; y; z) denotes the Gauss' Hypergeometric function. A more convenient form of (A8) can be found by using the transformation formulas [13; 9.132.2, 9.121.1]

r(1, 1/2; 3/2+q; - Z ) =

r(q + 3 / 2 ) F ( - 1 / 2 ) ,-.- 1,-,~, r(~2-~q+l/-~ z rtx, 1 / 2 - q ; 3/2; - z -1) + r(q + 1/2)r(1/2) Z-1/2(1 + Z_~) q.

r(a)r(q + 1)

(A9)

N. Ari / Continuuminterpolationof stressfields in brittlefracture

233

A final c o m p a c t f o r m for v(X, 0) can now be given as

[

2B

v ( X , 0) = D01(1 - X2)1/2[1

r ( q + l)

~¢/2 r ( q +

1 / 2 ) (C2 - 1 ) q F ( 1 ; 1 / 2 -

q; 3 / 2 ; - Y)

] d

+ b l - Z q ( c 2 --

X2)

q

(A10)

for I XI < 1 ,

Y = (1 - X 2 ) / ( C 2 - 1). The normal stress is given by eqs. (11), (14) and (18)

o~y/Oo=

- f0°°KA ( K ) cos K X d K .

(A11)

I n ( A l l ) we substitute for A ( K ) from (A1) and utilize the following integration by parts.

f l % o ( K t ) t ( C 2 - t 2 ) q-l~2 d t = J l ( K t ) ( K t ) ( C 2 - t 2 ) q - 1 / 2 / K 2 ] C _ f c d t J o ( K t ) ( _ 2 t ) ( q _ 1 / 2 ) ( C 2 _ t2)q_3/2/K2.

(A12)

Further we note that by (22)

( 2 / ~ r ) B Q ( C 2 - 1) q-'/2 = 1, and obtain = - 2 ( q - 1 / 2 ) ( C 2 - 1)l/2-qfCt2(C2 - t2)q-3(2p(t, X ) dt,

(A13)

where [13; 6.671.2]

P(t, X)=

Jl(Kt) cosKXdK=t

-1

1

1 -X/(X2-t2)

1/2

for for

t>X, t
(A14)

(A13) can now be evaluated easily by means of Hypergeometric functions, [13; 3.196.1]

Oyy/OO=

-1 + [x/(g

2 - 1) 1/2] F ( 1 , 1 / 2 ; q + 1 / 2 ; ( C 2 - 1 ) / ( g 2 - 1))

for I g l > C, (A15)

Oy,/a o = - 1 + [ X ( X 2 - 1)1/2] (2q - 1 ) F ( 1 , 3 / 2 - q; 3 / 2 ; ( X 2 - 1 ) / ( C 2 - 1)) for 1 < I X I < C.

(A16)

References [1] R.O. Ritchie and A.W. Thompson, "On macroscopic and microscopic analyses for crack initiation and crack growth toughness in ductile alloys", MetallurgicalTransactions 16A, 233-248 (1985). [2] B.R. Lawn and T.R. Wilshaw, Fractureof Brittle Solids, Cambridge University Press, London (1975). [3] J.R. Rice, "A path-independent integral and the approximate analysis of strain concentration by notches and cracks", J. Appl. Mech. 35, 379-386 (1968). [4] G.C. Sih, "A special theory of crack propagation", in: Methods of Analysis and Solutions to Crack Problems, ed. by G.C. Sih, Noordhoff International, Leiden, pp. 21-45 (1973). [5] A.C. Eringen, C.G. Speziale and B.S. Kim, "Crack-tip problem in nonlocal elasticity", J. Mech. Phys. Solids25, 339-355 (1977). [6] J.N. Sneddon, Mixed Boundary ValueProblems in PotentialTheory, North Holland, Amsterdam (1966). [7] E. Orowan, "Fracture and strength of solids", Rep. Progr. Phys. 12, p. 48 (1949).

234

N. Ari / Continuum interpolation of stress fields in brittle fracture

[8] G.I. Barenblatt, "The mathematical theory of equilibrium cracks in brittle fracture", Adv. App. Mech. 7, p. 55 (1962). [9] N. Ari and A.C. Eringen, Cryst. Lattice Defects&Amorph. Mat. 10, 33-38 (1983). [10] J.E. Sinclair and B.R. Lawn, "An atomistic study of cracks in diamond-structure crystals", Proc. Roy. Soc. Lond. A329, 83-103 (1972). [12] A. Kelly, Strong Solids, Oxford University Press, Oxford (1956). [13] I.S. Gradshyteyn and I.M. Ryzhik, Table oflntegrals, Series and Products, Academic Press, New York (1965).