- Email: [email protected]
The dynamic stresses induced by the propagation of skew cracks
COMPUTER MESONS IN APPLIED ME~A~CS 0 NORTH-HOLLAND PUBLISHING COMPANY
THE DYN~IC
AND ENGINNERING I5 (1978) 181-199
STRESSES INDUCED BY THE PROPAGATION OF SKEW CRACKS Jacob ABOUDI
Department of Solid Mechanics, Materialsand Structures, School of Engineering Tel-AvivUniversity,Ramt-Aviv, Israei Received 22 April I977 Revised manuscript received 30 September 1977 Skew crack propagation in an elastic medium due to the application of transient in-plane and anti-plane loads is investigated. For the anti-plane loading the motion of the crack is purely of mode III type, whereas in the case of in-plane loading a mixed mode type of fracture takes place in which both modes I and II occur. The method of solution is numerical and is based on a certain time-dependent transformation which maps the physical plane of the crack into an auxiliary plane in which the crack propagates collinearly with its propagating tip appearing always at the origin of the moving coordinate system. The transformed equations of motion are approximated by an implicit, three-level, finite difference system of equations of second-order accuracy, whose stability analysis is discussed. The reliability of the proposed method of solution is examined in several situations in which analytical results are known, and satisfactory agreement is achieved. Extension to smoothly curving cracks is discussed.
Introduction The basic problem in fracture mechanics consists of the computation of the elastic field in the vicinity of the tip of the crack and the application of a suitable fracture criterion. In recent reviews by Achenbach [ I] and Freund [ 21 analytical methods for the analysis of elastodynamic fields near propagating cracks are presented together with lists of references to related investigations. In [3] and [4] numerical methods of solution were presented for the dynamic propagation of cracks in a homogeneous elastic medium and along the interface of the two dissimilar media in mode I, II and III types of motion. In both 131 and [4] the reliab~ity and accuracy of the computed dynamic stresses and the associated stress intensity factors were checked by comparisons with available analytical solutions which can be derived in certain circumstances. In both the above two reviews and in the numerical treatments the cracks were assumed to propagate along their own plane without changing their direction. It often occurs however (e.g. when the material is under mixed mode loading) that the propagating crack deflects from its original plane and runs along a new path. The resulting skew crack propagation is mathematically much more difficult to analyze since it requires the solution of a new class of initial-boundary value problem involving both the original and the secondary crack forming the composite crack system. When the skew crack is under in-plane loading, mixed mode motion of the crack is obtained in which both the stress intensity factors in modes I and II are present. The fracture criterion for mixed mode situations is much more complicated, and two different approaches have recently evolved. The first approach is based on the energy release rate criterion for mixed mode, as was discussed recently by Nuismer 151. The second approach is that of Sih [6], who introduces the
182
J. Aboudi,
The dynamic stresses induced by the prop~gQtion ofskew
cracks
concept of the strain energy density function. Both these approaches were applied in [ 51 and [ 61 to elastostatic crack problems. Since the directions of mode I and II crack extension are not identical, it became necessary to deal with the complicated problem of a bent crack. The computation of the elastodynamic field in the vicinity of the tip of a skew crack under anti-plane loading was presented recently in [7]. The method employs the self-similarity of the field variables, which made it possible to reduce the problem to the solution of Laplace’s equation in a semi-infinite strip containing a slit. In this problem mixed modes do not occur, and the motion of the crack is purely of mode III type. However, this case is of less interest from an engineering point of view than the problem of a skew crack propagation under in-plane loading, which is more involved. Indeed, up to date no analytical information is available concerning this problem. We consider here several situations in which the boundary conditions must be imposed on surfaces which are not necessarily in the directions of the coordinates. Such situations arise, for example, when problems of elastic wedges or branching cracks in a solid are considered. It is known that such situations severely comphcate the treatment of the problem. For example, the in-plane dynamic motion of an elastic wedge whose surfaces are subjected to arbitrary loading is still unsolved, and only boundary conditions of smooth contact with a rigid plane allow an exact closed form solution (see Kostrov [8] and Zemell [9]). Similarly, the problem of a spherical cavity embedded in an elastic half-space is known to be inseparable in any orthogonal coordinate system which has both the free surface of the half-space and the loaded cavity surface as level surfaces. Even when treating such dynamic problems numerically by employing finite difference methods, one finds that essential difficulties arise when dealing with boundary conditions prescribed on surfaces which do not coincide with the grid points. We propose a numerical method of solution to the problem of a skew crack subjected to dynamic loading. The method is based on a specific time-dependent transformation which maps the physical plane of the composite propagating crack on an auxiliary plane in which the skew crack coincides with one of the axes, In the transform plane we obtain a collinear crack propagation problem, which is much easier to deal with than the original composite crack. This transformation facilitates the form of the boundary conditions at the expense of more complicated equations of motion. Another feature of this mapping is that the tip of the propagating crack appears always at the origin of the coordinate system of the transform plane. An implicit, three-level, finite difference system of equations is obtained which approximates the transformed equations of motion to second-order accuracy. The reliability of the resulting difference approximation is examined in several situations in which analytical results are known, and it is shown that satisfactory agreement is obtained. In the first situation we consider the impact of the surface of an elastic wedge by an anti-plane loading which produces propagating SH waves in the medium with singular stresses at the vertex of the wedge, as in a problem involving a crack where the stresses are singular at its tip. The region occupied by the wedge is mapped on a half-plane in the transform domain. The agreement between the numerical and analytical [lo] solutions is reasonable, and the study of the behavior of the numerical solution near the vertex is instructive for the extraction of the stress intensity factors from the crack tip stress field in the following crack problems. As another check of the present method of solution we consider the skew crack propagation caused by the incidence of a transient SH stress wave which was investigated in [ 71. In this case we obtain good correspondence by employing the asymptotic values for the stress field near the tip of the crack and the resulting analytical expression for the dynamic. stress intensity factor in order to check the numerical results for the extracted
J. Aboudi, %e dy~m~c stressesinduced by the prop~gutio~of skew cracks
183
stress intensity factors. Finally, the analytically unsolved problem of a skew crack propagation under in-plane loading is considered in which mixed modes of motion are present. The stress field and the associated stress intensity factors are exhibited graphically. The applicability of the present method of solution to other related problems is indicated (e.g. the nonuniform extension of a skew crack and the practical possibility that the original crack will not change its direction abruptly but will curve smoothly while it propagates).
1. Formulation
of the problem
Consider a homogeneous isotropic elastic medium with Lame constants X, p and density p. The two-dimensional elastodynamic equations of motion in the absence of body forces, are given by a2
a2
a2
at2
3X2
W
-u=d’-u+BB’----u+C’-----
a2
axay
21
(1)
’
where u = {u, u, w] is the displacement vector whose components respectively, t is the time, and the matrices A’, B’, C’ are
are in the x, y, z directions,
The compressional and shear wave speeds in the medium are respectively c1 = [(h f 21.r)/p]I’* and c2 = (PiP)“2. We consider situations in which the boundary conditions have to be imposed on surfaces which are not necessarily in the directions of the coordinates. Such situations arise, for example, when problems of elastic wedges or branching cracks in a solid are considered. As was mentioned before, such situations severely complicate the treatment of the problem. We shall deal with a class of problems for which the following transformation is useful:
-5= x -f(t)
3
q
=Y/fw
- 1,
z=z,
t=t,
(3)
where f(t) is a given function of time, and N(x) describes the shape of the surface on which the boundary conditions are imposed. This transfo~ation maps the curve y = H(x) on the axis n = 0 in the (t, n&plane, provided that H(x) # 0. It is obvious that in this plane the numerical treatment will be possible since the boundary conditions which were previously prescribed on a curved surface are now prescribed on the plane surface 7)= 0 of the new coordinate system. When f(t) f 0, we see that ($, q, z) is a moving coordinate system whose origin appears always at x = f(t). This situation is very useful when problems involving moving cracks are investigated since in the moving system the location of the tip of the crack appears always at the origin, provided that f(t) is chosen to represent the location of the crack tip relative to the origin of the stationary coordinates (x, y, z). Consequently, f(t) and?(t) will correspond to the velocity and acceleration of the expanding crack, respectively. In the new system (3) the elastodynamic equations (1) take the form
where A = A’ - f2(t)1, and the nonzero elements of the symmetric matrices B, C, D are given by
P&Y= (X + 2/.1)6; f
= 0 f We, PDYY
9
p/H* ,
pB,, = -(A + E.~&‘H,
PB,, = fibi + (A + 2i.0iH2,
PD,, = -(A + /.t)H’/.H2, ,oD,, = pD,, = /.d, ,
(5)
with
In the special case when H(x) is identically constant these equations reduce to those given in 131. Although after the transformation the equations of motion appear in a complicated form (4), the advantage gained is the simpli~catio~ of the form of the boundary surfaces on which the relevant boundary conditions must be satisfied. In a concrete problem the new elastodynamic equations (4) have to be solved in the transform plane together with the initial and boundary conditions. The resulting quantities shauld then be transformed to the original pIane. The exact form of the boundary conditions will depend on the problem. In general, they have in the problems considered in this paper the form
where (Y,0) is a polar coordinate system whose origin is located at (X = 0, y = H(0) 2 d), Oii are components of the stress tensor in this polar system, and gi are appropriate functions which describe the form of the applied loading, By introducing the explicit expressions of the stresses uBO,oOr, oez in terms of the displacement gradients W/at, au/aq, the equations in (6) take the form
S=O,
where
(7)
f. A&x& a,
=X+2~sin26,
,
Pie
@mmicstressesinduced by the ~rQ~g~~Q~Qf skew cracks
a,=-(a,H’+Ersin28,)1~,
a3 =-psin20,,
a4 = (h + 2p casai+
a5 = --sin 23, )
ix&= j-a&f’
a7=cos2&
a8 = -(a,
j
l&5
- a,H’)/H
>
+ cos a?,)/E) + H’aJH
,
a,, = (cos 8 I - .H’a,)/H _
ag=-sine,,
Xn the above equations 8, is the polar angle which traces the curve y = W(x), see fig. 1a.
It should be noted that the anti-plane motion in the transformed equations and boundary conditions is still uncoupled to the in-plane one.
Fig, 1. (a) A general configuration described by y = W(x), (b) a skew crack.
2. Numerical treatment The elastodynamic equations (4) are approximated by a finite difference system of equations by introducing a grid of mesh sizes A& An in the t, n directions, respectively, as well as a time increment At, and replacing ah the derivatives in (4) by their corresponding central difference expressions. Accordingly, the resulting finite difference approximation to (4) is accurate up to second order in the spatial and temporal increments, implicit (i.e. more than one grid points is involved at the advanced time level t + At, and of a three-level type (i.e. it is possible to compute the displacement vector in the region at time t + At from its values at the two previous time steps t and t - At t~o~~uut the region). Let #i,h n denote the finite difference appro~mation to u at the grid point (iA& jAn, PDF):
ui,j,,, 4 uW$, 14, nAt) = Nt,
77, t)
,
i,j=O,*l,k2
, “’ 9
n =
1,2, .. . .
(8)
186
J. Aboudi,
The dynamic stresses induced by the propagation of skew cracks
By replacing all the derivatives in (4) as described before, the following system of algebraic equations is obtained, which is of second-order accuracy (i.e. the truncation error is second order in the space and time increments): +
ui, j1
+ ~3 [Ui_l,j,n-1
-
‘3ui--1 I j t n+l
+ Gel
‘*
n+l
L”i+l, j+l,n
+j(t)flAt[Ui+l,j,n
-
f3Ui+l
ui+l,j,
n-1
+ ui-l,
, j , n+l
=
2[ 1 - AC:- BE;lUi,j,n -
‘i 1j ?n-l
1 +AE: [Ui+l,j,n+ Ui-l,j,n 1 + BE; [Ui,j+l,n+ Ui,j-l,n 1
j-l,n
-
ui-l,
j+l,n
-
ui+l,
-Ui-,,j,nI/‘+D’~“[‘i,j+l,n
j-1,n
II4
- ui,j_l,n l/2 >
(9)
where ei = At/At,
e2 = At/A7 ,
e3 = elf(t)/2
.
According to (9) it is possible to compute the displacement vector at time t + At from its values at the previous two time steps by solving at each time step for every j a tridiagonal system, for which a direct inversion algorithm can be employed and no iterations are needed. The numerical formulation of the boundary conditions (7) at n = 0 is obtained by approximating the displacement gradients by their central difference expressions for derivatives in the t direction, and by backward and forward difference expressions for derivatives in the 77direction, when approaching 17= 0 from n < 0 and n > 0, respectively. Consequently, at every time step t = n At the following system of algebraic equations in the unknowns ~~,e,~ is obtained (with eq = An/2Ai):
(10) =
‘[“6uj I ,n +“8ui >~1 >.I + AVg,/P9 Tl
E4”,[Wi+l,0,n - Wi-l,O,nl’
ulOWi,O,n = f"lOWi,Tl,n
+ 'qg3/p
*
The upper and lower signs in (10) correspond to the regions r) < 0 and n > 0, respectively, and the index i runs over all points along the &axis where the boundary conditions (6) are prescribed. In [4] the algebraic systems obtained from the numerical treatment of the boundary conditions were solved very easily by utilizing the Gauss-Seidel iterative procedure. The present case turns out to be more involved and the Gauss-Seidel method does not converge. Indeed, it can be verified that the sufficient conditions, which were shown to bp satisfied in the problem treated in [4], are not satisfied in (10). Accordingly, we solved the algebraic systems (10) by a direct method, utilizing the band structure of the matrices of coefficients of the unknowns. In order to apply the proposed difference scheme (9), it is necessary to show that it is stable (i.e. the difference between the exact and numerical solution of the difference equations (9) re-
J. Aboudi, The dynamic stresses induced by the propagation of skew cracks
187
mains bounded as the number of time steps increases); also an estimation is needed for the stability criteria. The stability analysis of (9) proceeds essentially as in [4] for the difference equations in the case when t;r(x) = cons& whose form is similar to (9). However, in the present case B and C are not constant matrices, and in (9) there exists the additional term involving the matrix D which is associated with the term au/an in (4). This term does not affect the stability criteria practically since it is known [ 111 that if D is a bounded operator, then the stability is not destroyed by the small perturbation which is represented by the term AfD in (9) (D is bounded since it was assumed that H(x) f 0). Consequently, the stability conditions appear according to [4] in the form e2 < l/max v, , r
r=l,2,3,
(11)
where E = cr = e2 (which is the usual case in practical computations), values of the real symmetric matrix M=Asin2p+Bsin2q+Csinp
and v, are the three eigen-
sinq cosp cosq,
(12)
with sin p, sin q varying independently between -1 and 1. We can estimate v, in terms of the spectral radii of the matrices A, B, C as follows: lMl
G ll All + IlBll + II CII ,
(13)
where II - II denotes the spectral norm. Since M is symmetric, its spectral radius R(M) is equal to its norm IIMII,so that R(M)dR(A)+R(B)+R(C),
(14)
since A, B, C are symmetric as well. The spectral radii R(A), R(B), Z?(C) can be easily computed from the characteristic equations of A, B and C. Since these matrices depend on the location x, conditions (11) must be satisfied for all the relevant range of variation of x. All the results given in this paper were obtained with the space increments AE = Aq = d/SO.
3. SH wave propagation in an impacted elastic wedge As a check to the proposed method of solution we apply it to the problem of an elastic wedge subjected to spatially uniform but time-dependent shear tractions which are applied to one of its surfaces while keeping the other one free. An analytical solution to this problem is available (Achenbach [lo]) that it can be employed to assess the present numerical method. Consider an elastic wedge occupying the region r > 0,O < 8 < cwwhose interior angle is CLThe wedge is at rest prior to t = 0, and at t = 0 time-dependent shear tractions are applied to one surface while keeping the other free. Thus the boundary conditions are 00
uez
=
I
0,
W) >
r>O,
e=o,
r-20,
e=(ll,
(15)
J. Aboudi,
188
7%e dynamic stresses induced by the propagation of skew cracks
with U(t) the Heaviside step function and u0 an amplitude The analytical solution to this problem is given by a0 cos 0 U(t - r sin
~ez(r, 8, t) =
L
e/c,) + C (r, 8, t)
factor.
,
0 < 0 < n/2 ) 0 2 K/2 )
(r, 6, t)
where (with k = cr/a > l/2, 0 = -cash-‘(c,t/r)) c
(r, 8, t) = -c2uo
sin(n/2k)/(nkr)U(t
- r/cz) j
FI(r/cT,
0) dr ,
r/c2
F,
(r/c* T, 0) = sinh(P/k)
sin(8/k)/{[cos(7r/2k)
- cosh(P/k)
cos(e/k)l*
For this problem the proposed numerical treatment can be employed to (3) with f(t) = 0) the boundaries of the wedge defined by
H(x)
-d,
.Y >
-d + .Y tan Q ,
XGO,
0
+ [sinh(P/k) by transforming
sin(e/k)l*}
.
(according
)
(17)
=
to the surface q = 0 of the half-plane n < 0 on which the interior region of the wedge is mapped. In fig. 2 the numerical and analytical (16) solutions for the shear stresses uBz are shown versus r/d along the ray 8 = 7r/2 when c2 t/d = 0.5 and 1 for (Y= 5x/4. The numerical solution was obtained by applying (9) with g, = g, = 0 in (lo), whereas g, was taken from (15). It is clearly seen that a reasonable agreement between the two solutions exists. Whereas the analytical solution predicts singular stresses near the vertex of the wedge (r = 0), the numerical solution yields a finite value there as can be expected. In previous numerical solutions for dynamic crack problems [ 3,4]
+-+olJ--’ 0
I
-
r/d
0
Fig. 2. Numerical solution (solid lines) and analytical solution (16) (dashed 0 = n/2. in a wedge under anti-plane loading when czt/d = 0.5 and 1.
-
0.5 r/d
-+----I
lines) for the shear stresses agZ versus r/d along
J. Aboudi, The dynamic stressesinduced by the propagation of skew cracks
189
the numerically obtained finite values near the tip of the crack were shown to correspond to the asymptotic solution near the tip of the crack and can therefore be employed to extract the stress intensity factors. Let us explore therefore these finite values for the present problem of the elastic wedge. These values are u~~/u,-,= 1.89 and 2.38, which are obtained when c2 t/d = 0.5 and 1, respectively, at the point Y= At, 8 = 7r/2 and correspond to the analytical solution (16) at Y= rat, 8 = n/2, where y = l/2. Consequently, they can be interpreted as the asymptotic solution of the wedge problem at a distance Y= yAl from the vertex. In the crack problems discussed in the sequel this type of information will be useful in extracting the associate dynamic stress intensity factors.
4. Skew crack propagation
caused by the incidence
of SH waves
Consider a homogeneous elastic medium containing a semi-infinite crack which lies along the ray 8 = (Y(fig. 1b). At time t = 0 a plane incident SH wave strikes the tip of the semi-infinite crack. The incident wave has the form
wqx,y,
= W,TU(T)
)
(18)
with
where 8, defines the orientation of the front of the plane wave, and w0 is an amplitude factor. The incident stress wave generated by (18) is given by
~$2= (TV cog8 - e,) u(7) ,
(19)
where 00 =
w&,
’
It is assumed that the crack starts to propagate at the instant t = 0 with a velocity f(t) < c2 along the ray 8 = 0. Therefore, in contrast to previous works on crack propagation in their own plane it is assumed here that the crack deviates from its original plane with an angle K - (Y.For this problem an asymptotic expression for the stress field in the vicinity of the tip of the crack - in the case of constant rate of extension f(t) = u - was given by Achenbach and Varatharajulu [ 7 I. This expression has the form u&2,$,
t) -
(t/R)“* k/W, 4) >
(20)
where (R, #> is a new set of polar coordinates centered at the crack tip (fig. lb), k, is given by k, =
(1 - u~/c;)~‘~ K/u”* ,
(21)
and @(u, $J)describes the universal dependence of the stress field on the polar angle 4 in the vicinity of the crack. The explicit form of Q(u, 4) appears in [ 71 and [ 121. When $J= 0, @(u, 0) = 2l’*. In [ 71, the function K in (21) is plotted versus the skew angle for two values of the incidence angle eo.
190
J. Aboudi,
The dynamic stresses induced by the propagation of skew cracks
In the special case when (Y= K and there is a stationary crack (u = 0), the asymptotic expression for the diffracted stress (20) reduces to
2% (1 + sin 0e)“* (c2 t/R)“* cos(@/2) .
Q,(R, 4, t) = uoz(r,8, t) - 7
(22)
We can utilize the above asymptotic expressions in order to check the proposed numerical method of solution in the present circumstances as well as the associated dynamic stress intensity factor which is extracted from this numerical solution. The skew crack problem can be handled according to the present method by transforming the surfaces of the composite crack to the planes n = TO in the plane (t, q) to which the original domain (x, y) is transformed according to (3) with f(r) representing the velocity of the crack. The function H(x) of (3) appropriate to the skew crack is given by d - x tan a ,
XGO,
d,
x20.
H(x) =
(23)
In the transform plane the tip of the porpagating crack appears always at the origin of coordinates t = 0, 7)= 0. Since the stress intensity factor is determined from the singular behavior of the stress field in the vicinity of the crack, and the latter is contributed from the diffracted part of the field (not including the incident wave), it would be sufficient to construct the solution of an initially undisturbed medium with the boundary conditions (6) in which g, = g, = 0 and
g&O
= a’e’!,
tgo,
q=o,
tao,
(24)
provided that the independent variables of ~$2 are properly transformed. It is instructive to consider first the problem of the diffraction of an incident SH plane wave by a semifinite stationary crack lying along x < 0, y = d (i.e. with a = 1~in fig. 1b) in an infinite elastic medium. This problem possesses an analytical solution due to De Hoop and reviewed in [ 131. In the case when the plane wave front is parallel to the surface of the crack (i.e. 0, = CY = K), this solution has the form (see [ lo])
qJR> $, z) = q&
0, z) = JW - R/c,) >
-n/2 < @< a/2 )
(25)
where (2R/cz)l’* 71J/20, = (t - R/cZ)l’* { [ 1 - J, tan-‘(1 /J,)] cos(#/2 - 7r/4) - 11 - J, tan-’ (l/J,)]
(26)
sin(@/2 - a/4)) ,
with Jf = R(1 - sin @)/[c,(t - R/c,)]
and.
Ji = R( 1 + sin 4)/[cl(t
- R/c,)]
.
Accordingly, the asymptotic expression for u+ near the tip of the crack has the form
J. Aboudi,
The dyrumic
191
stresses induced by the propagation of skew cracks
uaZ(R, @,Z)-2”2 o0 t [(C,t/R)“2 - n(1 - sin $~)I’~/21cos($/2 - 7r/4) - [(c,t/R)“’
- n(1 + sin @1)“~/2]sin(@/2 - 71/4)}/7r ,
(27)
which reduces to (22) when 8, = 71and R/c, t e 1. Let us define the dynamic stress intensity factors associated with the three modes of motion of the crack in the form Ti (t) = U,j(R, 0, t)R”’
(28)
)
for small R such that i = 1 when j = #, i = 2 when j = R and i = 3 when j = z. According to (28) the dynamic stress intensity factor derived from (27) has the form T, (t) = 20, (cz t)“2 /T ,
(29)
which can also be obtained as a special case of the skew crack problem with (Y= 0, = K and u = 0. In fig. 3 the numerical and analytical solutions for uez are shown when c2 t/d = 0.5 and 1 versus R/d along the ray 4 = 0. It is clearly seen that a good correspondence between the solutions is ob-
II
il
C21/d ao.5
c,t/d=l.o
01
0
(b)
0.5
cpt/d-
I
I
Fig. 3. (a) Numerical solution (solid lines) and analytical solution (25) (dashed lines) for the shear stresses oez versus R/d along @= 0 for a stationary crack with 01= n under anti-plane loading when czt/d = 0.5 and 1, (b) the dynamic stress intensity factor curves based on the numerical solution (solid lines) and analytical solution (29) (dashed lines) for the stationary crack.
192
J. Aboudi. 77~ dynamic stresses induced by the ~ro~g~ti~~
of skewcracks
tained. In the vicinity of the tip of the crack the analyticaf solution (25) predicts a singularity at R = 0, whereas the numerical solution yields a finite value for the stress at the grid point (4‘= A& 7 = 0), which is the closest grid point ahead of the tip of the crack. For example, we obtain at this point the values of ~~~/u~ = 3.23 and 4.98 when c,tld = 0.5 and I, respectively, whereas the analytical solution (25) yields at the same point the corresponding values 2.25 and 3.54 (the asymptotic expression (27) yields 2.18 and 3.50). In order to explore these values, we compared them with the analytical solution (25) (or the asymptotic expression (27)) evaluated at several locations in the vicinity of the tip of the crack at the same times and found that they correspond to the exact solution at t = ?A.$, where y = l/2. The numerical solution at the point (4:= Ak, n = 0) can be employed to get the dynamic stress intensity factor of the crack. According to (28) we compute the values of T,(t) by multiplying uaz(t = At, n = 0, t) by the square root oft = rA[ with y = l/2. However, since the difference between the numerical solution for T3(t) obtained by (28) and the analytical one given by (29) is of order of R1’2 (see (27)), we implement this when computing T3(t) by employing (28) and adding f~A$)l’~ to the result. Fig. 3 shows curves for T,(t) based on the numerical and analytical (29) solutions, and good agreement between them can be well observed. We turn presently to the main problem of skew crack propagation. In fig. 4 the numerical solu= 0 is shown for the original crack orientation tion for the stresses uot versus R/d along the ray r#~ c$/d = I.0
c,ttd 105
O_ -
0.5 R/d
i
4 I
0 (0)
/ -
0.5 R/d
I
v/cz=o 2
-
0.5 c,tld
I
(b) Fig. 4. (a) Numerical solution for the shear stresses 0~~ versus R/d along Q,= 0 for a skew crack under anti-plane loading with IX= 3n/4, propagating with the constant velocities u/c2 = 0.2 and 0.5 when czt/d = 0.5 and 1, (b) the corresponding dynamic stress intensity factor curves based on the numerical solution (solid lines) and analytical solution (dashed lines).
193
J. Aboudi, The dy~~~c stresses induced by the propQgQ?~o~ ofskew cracks
01= 3n/4 (fig. 1b) when c,t/d = 0.5 and 1 for the two cases in which the propagation velocities are constant and given by u/c, = 0.2 and 0.5. The dynamic stress intensity factor curves T3(t) shown in the figure were extracted from the numerical solution for the stress uoz obtained at the closest grid point ahead of the crack tip (4‘= Al, q = 0) as explained before. These curves are compared with stress intensity factors based on the analytical treatment in [7] and were obtained by applying (28) to (20). It is well observed that a satisfactory agreement exists between the two methods of solution. It should be mentioned that the present choice of OL= 37r/4 represents a fairly significant deviation of the crack from its original plane since in practical situations the skew angle 77- (Y does not seem to exceed the value of ~16 [ 121 so that even a better agreement between the two methods of solution would be accomplished.
5. Skew crack propagation caused by the sudden application of normal tractions on its surface We consider presently the analytically yet unsolved problem of the skew crack propagation caused by a sudden application at t = 0 of equal and opposite uniform normal tractions to the surfaces of a semi-infinite crack. The crack lies originally along the ray 0 = a! (see fig. 1b) in an initially undisturbed medium. The applied tractions produce plane waves and diffracted cylindrical waves whose center is located at the original crack tip. It is assumed that the semi-infinite crack extends at a velocity f(t) out of its plane at the instant of the application of the tractions. This problem is solved according to the proposed method of solution by using (3) to transform the surfaces of the crack to the planes Q = TO in the transform domain as described previously. The appropriate function H(x) is given by (23). The boundary conditions are given by (6), where
g,=g,=o,
g&5 f) = -u(J
[
1- U(x)1WI ,
l
q=“o,
(30)
with x = .$+f(t). According to (30) the normal tractions are applied on the surfaces of the original crack only, while keeping the newly produced surfaces of the secondary crack free of tractions. Before considering the skew crack problem it is instructive here too to treat the case of a collinearly propagating semiin~nite crack in an infinite solid which was treated by Baker [ 141 and for which some analytical results are known. The boundary conditions are given by (6) with g, = g, =Oandg, = -u,, U(t). The asymptotic expression for uGq,,which is valid near the tip of the crack, is given by [ 141: (31) where (R, 4) is the polar coordinate system centered at the crack tip. In (3 1) F+(O) is a complicated expression which depends on the material parameters and the uniform velocity u of the crack, and \k(u, $) represents the universal dependence of the singular stresses on the polar angle (bin the vicinity of the crack, and the propagation velocity u. For Qj= 0 (i.e. for observation points along the path of the crack) *(u, 0) = 1, whereas \k(O, #) = cos3 ($12). The explicit form of U($J, u) can be identified in [ 121. According to (28) the analytical expression of the dynamic stress intensity factor T, (t) is given by T,(t) = 20&l
- u/c,)CZt]1’2/aF+(0).
(32)
194
J. Aboudi, The dynamic stresses indrtced by the Fmpagation ofskew cracks
MODE II
MODE I
4-
r&c, =0.542
3
3-
c,t/d ~1.0
c,t/d = 1.0 / b” ’ 2z / b
0.5
I
-
0.3 R/d
r
0.5
I.t::::.:--‘-: 0
(a)
0.5c&/d
Fig. 5. (a) Numerical solution for the stresses ow versus R/d along # = 0 for a stationary crack with a = n under in-plane normal loading when clt/d = 0.5 and 1 (the dynamic stress intensity factor curve based on the numerical solution (solid lines) is compared with the analytical one based on (32)), (b) numerical solution for the stresses 00~ versus R/d along @J= 0 for a stationary crack with o = n under in-plane shear loading when clt/d = 0.5 and 1, and the corresponding dynamic stress intensity curve (the material constants in this case have the value ~/et = 0.542).
Unless otherwise mentioned, all the results given in this section correspond to the choice CJC, = I/2. In fig. 5 the stresses ucbe for a stationary crack (v = 0) lying along the ray a = K are shown versus R/d along Cz, = 0 when c,t/d = 0.5 and 1. The associated numerical stress intensity factor curve is extracted from the values of u*@at the closest grid point ahead of the tip of the crack, i.e. at (t = A$, n = 0), by applying (28) with R = ?A$ and y = l/2 (since here too the numerical values of o++ at this point correspond to the values of crbQ,obtained from the asymptotic expression (3 1)). This curve is compared in the figure with the corresponding one based on (32). It is seen that a very satisfactory correspondence between the two methods exists. As we need in the sequel to extract the dynamic stress intensity factor T,(t) (associated with a mode II type of motion) from the shear stresses ueR induced near the tip of the crack, let us apply the boundary conditions (6) with g, =g, = 0 and the shear impact g, = -u0 u(t) in order to produce a mode II type of opening. In this specific case we chose the elastic material with cJc, = 0.542 (steel) since for this value Sih, Embley and Ravera [ 151 present the corresponding stress intensity factor appropriate to this problem. In fig. (5) the stresses CJ+~versus R/d along 4 = 0 are shown for a stationary crack (u = 0, (Y= n) when et t/d = 0.5 and 1. The stress intensity factor r,(t) is extracted from the numerical
J. Aboudi, The dynamic stressesinduced by the propagation of skew cracks
195
solution of utiR induced by the vicinity of the crack tip according to (28) and adding the contribution (~a#‘~, y = l/2, as in the previous anti-plane case (see also [4]). The resulting values of T2(f) as exhibited in fig. 5 are in good correspondence with the curve shown in [ 151, which is based on an analytical treatment of the problem and a numerical inversion of the Laplace transforms for the evaluation of the stress intensity factor. For example, when c,t/d = 1, our numerical whereas we obtain the value 0.55 from the curve in solution yields the value 0.52 for T,(t)/u,d”2,
[151. It is still necessary to check the reliability of the numerical scheme (9)~(10) in the in-plane loading condition of a crack when the transformation (3) is involved since in the previous cases (exhibited in fig. 5) H(x) was identically constant so that this transformation was not actually active. Since no analytical results have been obtained for the skew crack problem under in-plane loading, we chose to perform such a check by considering the previous stationary crack problem with sudden normal tractions acting on its boundary but lying now along the ray cy= 3n/4 (rather than (Y= rr as in the previous cases). The stress intensity factor T,(t) which is extracted from the closest grid point ahead of the tip of the crack, i.e. from the point (t = A& 77= 0), can be compared with the analytical results given by (32). However, since the location of this point corresponds in the present situation to C#= n/4, we must implement (28) by taking into account the fact that the numerical solution for u6@obtained at this point corresponds to the asymptotic expression (3 1) in which \k(O, 4) = cos3($/2) = l/23’2, R = rag, y = l/2. Accordingly, the numerical values of T,(t) should be extracted from uGOafter dividing the latter by this value of \k (in order to recover the value of uQOat 4 = 0) and employing (28). In fig. 6 the induced stresses uGQ,by the stationary crack, are shown versus R/d along $ = a/4 when c,t/d = 0.5 and 1. In this figure the numerical values of T,(t) are shown and contrasted with the analytical values based on (32). It is clearly seen that the method of solution is reliable in the present situation too. Having established the reliability of the numerical treatment and the extraction of the dynamic stress intensity factors in mode I and II types of motion of the crack, it is possible to investigate
015 -
R/d
(a)
-
c,t /d
(b)
Fig. 6. (a) Numerical solution for the stresses om versus R/d along @= 0 for a stationary crack with Q = 3n/4 under in-plane normal loading when ct t/d = 0.5 and 1, (b) the dynamic stress intensity factor curve based on the numerical solution (solid lines) compared with the analytical one based on (32).
196
J. Aboudi, Eke dynamic stresses induced by the p~pagat~on ofskew cracks
c,i/d =0.5
c,t/d
I 8
0.5 -
0
-
I I
b
I
0
R/d
0.5 c,t/-d
~1.0
I 0.5 --R/d
-
05 c, t/d
i I
I
Fig. 7. Numerical solution for the stresses DQ,++, 06~ versus R/d along 0 = 0 for a skew crack with a = 170’ under in-plane loading (30) when ~1 t/d = 0.5 and 1 and the corresponding dynamic stress intensity factor curves in modes I and II (the secondary crack propagates with the velocities u/cl = 0.1 and 0.2).
now the main problem of a skew crack propagation under in-plane loading (30). As was mentioned before, even though the primary crack is under pure mode I loading, we obtain in general a mixed mode problem in which both modes I and II are present. In figs. 7-9 we present the normal and shear stress uc4@,uQR versus R/d computed along the ray C#= 0 for a crack which was initially oriented along (I!= 170”, 150” and 135” when et t/d = 0.5 and 1. The propagation velocities of the secondary crack are taken to be uniform with the values u/c, = 0.1 and 0.2. The associated stress intensity factors T,(t), T,(t) appropriate to modes I and II, respectively, are also shown in the figures. The stress intensity factor 1”,(t) is extracted from the values of uo9 at the grid point ahead of the tip, whereas T,(t) is extracted from uGR at that point. It is clearly observed that the stress intensity factors decrease as the propagating velocity increases. On the other hand, T,(t) increase as the skew angle ar increases, whereas T,(r) decreases as cxincreases. These results are not unexpected since in the case of a collinear crack propagation, the shearing stress uoR vanishes along the ray 4 = 0 so that r,(t) = 0 and the mixed mode propagation becomes a pure mode I type.
6. Conclusions A numerical method of solution is presented which is able to handle the problem of the dynamic curved crack propagation under anti-plane and in-plane loading. The method is based on a certain
J. Aboudi, Tire dynarni~stressesinduced by the p~pagation ofskew cracks
C, t/d
I
I 1
-
0
0.5
c,t/d +I.0
=0.5
!
I
I
R/d
0.5 -
197
I
c,t/d
Fig. 8. Same as fig. 7 but with the original crack orientation
01 =
1 SO”.
time-dependent transformation which maps the physical plane of the crack into an auxiliary plane in which the crack appears propagating col~nearly. The reliability of the method of solution is checked by comparisons with some analytical results. Several extensions are possible. For example, the method of solution is formulated such that it can handle nonuniform skew crack propagation. In this paper we assumed that the original crack changes its direction abruptly. Common experience indicates that a crack will generally follow a smooth curved path. Furthermore, it has been shown recently by Willis [ 161 that a straight crack subjected to anti-plane static loading cannot abruptly change its direction of propagation even if it is loaded asymmet~c~ly but may curve smoothly. In such cases in which the crack propagates along a smooth curved path the present numerical method of solution will yield even more accurate results since, as a rule, numerical schemes are less accurate when dealing with functions possessing piecewise continuous derivatives (as in the case of a skew crack). On the other hand, analytical treatments for cracks that curve in an arbitrary fashion are lacking at present, so that numerical methods developed for such problems seem to be valuable. In order to explain some particular phenomena in fracture mechanics whose purpose is to characterize the actual behavior of cracks by using a simple mathematical crack model, some approaches for making the crack model more precise have been proposed. In [ 171 those approaches
19%
i
$=, 8 b 0 I
-
-
0.5 R/d
c,tM
Fig. 9. Same as fig. ‘7 but with the original crack orientation o = 135”.
are divided into six different classes of problems. One of those is the class of nonlinear shaped cracks to which the skew and curved types of cracks belong. This type of crack configurations is encountered, for example, when a mixed mode propagation of a crack is involved since generally the crack does not fottow a straight path but a curved or a bent path. Several theoreticaf and experimental evidents af particular nonlinear shaped cracks are discussed in [ 1i’] and in the references cited there. Consequently, the present paper praposes a numerical method of solution applicable to several types of dynamic crack propagation which belong to the general category of nonlinear crack configurations.
Acknowledgment This research was performed under contract AFOSR-‘76-3014 of the U.S. Air Force.
References [I ] J.D. Achenbach, The dynamic effects in brittle fracture, in: S. Nemat-Nasser (ed.) Mechanics today, 1 (Pergamon Press, 1974).
J. Aboudi, The dynamic stressesinduced by the propagation of skew cracks
199
[2] L.B. Freund, The analysis of elastodynamic crack tip stress fields, in: S. Nemat-Nasser (ed.), Mechanics today, 3 (Pergamon Press, 1976). [ 31 J. Aboudi, Numerical solution of dynamic stresses induced by moving cracks, Comp. Meths. Appl. Mech. Eng. 9 (1976) 301-316. [4] J. Aboudi, The dynamic stresses induced by moving interfacial cracks, Comp. Meths. Appl. Mech. Eng. 10 (1977) 303-323. 15) R.J. Nuismer, An energy release rate criterion for mixed mode fracture, Int. J. Fracture 11 (1975) 245-250. [6] G.C. Sih, Strain-energy-density factor applied to mixed mode crack problems, Int. J. Fracture 10 (1974) 305-321. [7 ] J.D. Achenbach and V.K. Varatharajulu, Skew crack propagation at the diffraction of a transient stress wave, Q. Appl. Math. 32 (1974) 123-135. [ 81 B.V. Kostrov, Diffraction of a plane wave by a smooth rigid wedge in an unbounded elastic medium in absence of friction, J. Appl. Math. Mech. 30 (1966) 244-250. [9] S.H. Zemell, Diffraction of elastic waves by a rigidsmooth wedge, SIAM J. Appl. Math. 29 (1975) 582-596. [lo] J.D. Achenbach, Shear waves in an elastic wedge, Int. J. Solids Structs. 6 (1970) 379-388. [ 111 R.D. Richtmeyer and K.W. Morton, Difference methods for initial-value problems, 2nd ed. (Interscience, 1967). (121 J.D. Achenbach, Wave propagation, elastodynamic stress singularities, and fracture, in: W.T. Koiter (ed.) IUTAM Congress, Delft (1976) 71-87. [ 131 J.D. Achenbach, Wave propagation in elastic solids (North-Holland, 1973). (141 B.R. Baker, Dynamic stresses created by a moving crack, J. Appl. Mech. 29 (1962) 449-458. [IS] G.C. Sih, G.T. Embley and R.S. Ravera, Impact response of a finite crack in plane extension, Int. J. Solids Structs 8 (1972) 977-993. [ 161 J.R. Willis, A discussion of crack-forking in anti-plane strain deformation, Int. J. Fracture 11 (1975) 489-493. [17] H. Kitagawa, R. Yuuki and T. Ohira, Crack-morphological aspects in fracture mechanics, Eng. Fracture Mech. 7 (1975) 515-529.
THE DYN~IC
AND ENGINNERING I5 (1978) 181-199
STRESSES INDUCED BY THE PROPAGATION OF SKEW CRACKS Jacob ABOUDI
Department of Solid Mechanics, Materialsand Structures, School of Engineering Tel-AvivUniversity,Ramt-Aviv, Israei Received 22 April I977 Revised manuscript received 30 September 1977 Skew crack propagation in an elastic medium due to the application of transient in-plane and anti-plane loads is investigated. For the anti-plane loading the motion of the crack is purely of mode III type, whereas in the case of in-plane loading a mixed mode type of fracture takes place in which both modes I and II occur. The method of solution is numerical and is based on a certain time-dependent transformation which maps the physical plane of the crack into an auxiliary plane in which the crack propagates collinearly with its propagating tip appearing always at the origin of the moving coordinate system. The transformed equations of motion are approximated by an implicit, three-level, finite difference system of equations of second-order accuracy, whose stability analysis is discussed. The reliability of the proposed method of solution is examined in several situations in which analytical results are known, and satisfactory agreement is achieved. Extension to smoothly curving cracks is discussed.
Introduction The basic problem in fracture mechanics consists of the computation of the elastic field in the vicinity of the tip of the crack and the application of a suitable fracture criterion. In recent reviews by Achenbach [ I] and Freund [ 21 analytical methods for the analysis of elastodynamic fields near propagating cracks are presented together with lists of references to related investigations. In [3] and [4] numerical methods of solution were presented for the dynamic propagation of cracks in a homogeneous elastic medium and along the interface of the two dissimilar media in mode I, II and III types of motion. In both 131 and [4] the reliab~ity and accuracy of the computed dynamic stresses and the associated stress intensity factors were checked by comparisons with available analytical solutions which can be derived in certain circumstances. In both the above two reviews and in the numerical treatments the cracks were assumed to propagate along their own plane without changing their direction. It often occurs however (e.g. when the material is under mixed mode loading) that the propagating crack deflects from its original plane and runs along a new path. The resulting skew crack propagation is mathematically much more difficult to analyze since it requires the solution of a new class of initial-boundary value problem involving both the original and the secondary crack forming the composite crack system. When the skew crack is under in-plane loading, mixed mode motion of the crack is obtained in which both the stress intensity factors in modes I and II are present. The fracture criterion for mixed mode situations is much more complicated, and two different approaches have recently evolved. The first approach is based on the energy release rate criterion for mixed mode, as was discussed recently by Nuismer 151. The second approach is that of Sih [6], who introduces the
182
J. Aboudi,
The dynamic stresses induced by the prop~gQtion ofskew
cracks
concept of the strain energy density function. Both these approaches were applied in [ 51 and [ 61 to elastostatic crack problems. Since the directions of mode I and II crack extension are not identical, it became necessary to deal with the complicated problem of a bent crack. The computation of the elastodynamic field in the vicinity of the tip of a skew crack under anti-plane loading was presented recently in [7]. The method employs the self-similarity of the field variables, which made it possible to reduce the problem to the solution of Laplace’s equation in a semi-infinite strip containing a slit. In this problem mixed modes do not occur, and the motion of the crack is purely of mode III type. However, this case is of less interest from an engineering point of view than the problem of a skew crack propagation under in-plane loading, which is more involved. Indeed, up to date no analytical information is available concerning this problem. We consider here several situations in which the boundary conditions must be imposed on surfaces which are not necessarily in the directions of the coordinates. Such situations arise, for example, when problems of elastic wedges or branching cracks in a solid are considered. It is known that such situations severely comphcate the treatment of the problem. For example, the in-plane dynamic motion of an elastic wedge whose surfaces are subjected to arbitrary loading is still unsolved, and only boundary conditions of smooth contact with a rigid plane allow an exact closed form solution (see Kostrov [8] and Zemell [9]). Similarly, the problem of a spherical cavity embedded in an elastic half-space is known to be inseparable in any orthogonal coordinate system which has both the free surface of the half-space and the loaded cavity surface as level surfaces. Even when treating such dynamic problems numerically by employing finite difference methods, one finds that essential difficulties arise when dealing with boundary conditions prescribed on surfaces which do not coincide with the grid points. We propose a numerical method of solution to the problem of a skew crack subjected to dynamic loading. The method is based on a specific time-dependent transformation which maps the physical plane of the composite propagating crack on an auxiliary plane in which the skew crack coincides with one of the axes, In the transform plane we obtain a collinear crack propagation problem, which is much easier to deal with than the original composite crack. This transformation facilitates the form of the boundary conditions at the expense of more complicated equations of motion. Another feature of this mapping is that the tip of the propagating crack appears always at the origin of the coordinate system of the transform plane. An implicit, three-level, finite difference system of equations is obtained which approximates the transformed equations of motion to second-order accuracy. The reliability of the resulting difference approximation is examined in several situations in which analytical results are known, and it is shown that satisfactory agreement is obtained. In the first situation we consider the impact of the surface of an elastic wedge by an anti-plane loading which produces propagating SH waves in the medium with singular stresses at the vertex of the wedge, as in a problem involving a crack where the stresses are singular at its tip. The region occupied by the wedge is mapped on a half-plane in the transform domain. The agreement between the numerical and analytical [lo] solutions is reasonable, and the study of the behavior of the numerical solution near the vertex is instructive for the extraction of the stress intensity factors from the crack tip stress field in the following crack problems. As another check of the present method of solution we consider the skew crack propagation caused by the incidence of a transient SH stress wave which was investigated in [ 71. In this case we obtain good correspondence by employing the asymptotic values for the stress field near the tip of the crack and the resulting analytical expression for the dynamic. stress intensity factor in order to check the numerical results for the extracted
J. Aboudi, %e dy~m~c stressesinduced by the prop~gutio~of skew cracks
183
stress intensity factors. Finally, the analytically unsolved problem of a skew crack propagation under in-plane loading is considered in which mixed modes of motion are present. The stress field and the associated stress intensity factors are exhibited graphically. The applicability of the present method of solution to other related problems is indicated (e.g. the nonuniform extension of a skew crack and the practical possibility that the original crack will not change its direction abruptly but will curve smoothly while it propagates).
1. Formulation
of the problem
Consider a homogeneous isotropic elastic medium with Lame constants X, p and density p. The two-dimensional elastodynamic equations of motion in the absence of body forces, are given by a2
a2
a2
at2
3X2
W
-u=d’-u+BB’----u+C’-----
a2
axay
21
(1)
’
where u = {u, u, w] is the displacement vector whose components respectively, t is the time, and the matrices A’, B’, C’ are
are in the x, y, z directions,
The compressional and shear wave speeds in the medium are respectively c1 = [(h f 21.r)/p]I’* and c2 = (PiP)“2. We consider situations in which the boundary conditions have to be imposed on surfaces which are not necessarily in the directions of the coordinates. Such situations arise, for example, when problems of elastic wedges or branching cracks in a solid are considered. As was mentioned before, such situations severely complicate the treatment of the problem. We shall deal with a class of problems for which the following transformation is useful:
-5= x -f(t)
3
q
=Y/fw
- 1,
z=z,
t=t,
(3)
where f(t) is a given function of time, and N(x) describes the shape of the surface on which the boundary conditions are imposed. This transfo~ation maps the curve y = H(x) on the axis n = 0 in the (t, n&plane, provided that H(x) # 0. It is obvious that in this plane the numerical treatment will be possible since the boundary conditions which were previously prescribed on a curved surface are now prescribed on the plane surface 7)= 0 of the new coordinate system. When f(t) f 0, we see that ($, q, z) is a moving coordinate system whose origin appears always at x = f(t). This situation is very useful when problems involving moving cracks are investigated since in the moving system the location of the tip of the crack appears always at the origin, provided that f(t) is chosen to represent the location of the crack tip relative to the origin of the stationary coordinates (x, y, z). Consequently, f(t) and?(t) will correspond to the velocity and acceleration of the expanding crack, respectively. In the new system (3) the elastodynamic equations (1) take the form
where A = A’ - f2(t)1, and the nonzero elements of the symmetric matrices B, C, D are given by
P&Y= (X + 2/.1)6; f
= 0 f We, PDYY
9
p/H* ,
pB,, = -(A + E.~&‘H,
PB,, = fibi + (A + 2i.0iH2,
PD,, = -(A + /.t)H’/.H2, ,oD,, = pD,, = /.d, ,
(5)
with
In the special case when H(x) is identically constant these equations reduce to those given in 131. Although after the transformation the equations of motion appear in a complicated form (4), the advantage gained is the simpli~catio~ of the form of the boundary surfaces on which the relevant boundary conditions must be satisfied. In a concrete problem the new elastodynamic equations (4) have to be solved in the transform plane together with the initial and boundary conditions. The resulting quantities shauld then be transformed to the original pIane. The exact form of the boundary conditions will depend on the problem. In general, they have in the problems considered in this paper the form
where (Y,0) is a polar coordinate system whose origin is located at (X = 0, y = H(0) 2 d), Oii are components of the stress tensor in this polar system, and gi are appropriate functions which describe the form of the applied loading, By introducing the explicit expressions of the stresses uBO,oOr, oez in terms of the displacement gradients W/at, au/aq, the equations in (6) take the form
S=O,
where
(7)
f. A&x& a,
=X+2~sin26,
,
Pie
@mmicstressesinduced by the ~rQ~g~~Q~Qf skew cracks
a,=-(a,H’+Ersin28,)1~,
a3 =-psin20,,
a4 = (h + 2p casai+
a5 = --sin 23, )
ix&= j-a&f’
a7=cos2&
a8 = -(a,
j
l&5
- a,H’)/H
>
+ cos a?,)/E) + H’aJH
,
a,, = (cos 8 I - .H’a,)/H _
ag=-sine,,
Xn the above equations 8, is the polar angle which traces the curve y = W(x), see fig. 1a.
It should be noted that the anti-plane motion in the transformed equations and boundary conditions is still uncoupled to the in-plane one.
Fig, 1. (a) A general configuration described by y = W(x), (b) a skew crack.
2. Numerical treatment The elastodynamic equations (4) are approximated by a finite difference system of equations by introducing a grid of mesh sizes A& An in the t, n directions, respectively, as well as a time increment At, and replacing ah the derivatives in (4) by their corresponding central difference expressions. Accordingly, the resulting finite difference approximation to (4) is accurate up to second order in the spatial and temporal increments, implicit (i.e. more than one grid points is involved at the advanced time level t + At, and of a three-level type (i.e. it is possible to compute the displacement vector in the region at time t + At from its values at the two previous time steps t and t - At t~o~~uut the region). Let #i,h n denote the finite difference appro~mation to u at the grid point (iA& jAn, PDF):
ui,j,,, 4 uW$, 14, nAt) = Nt,
77, t)
,
i,j=O,*l,k2
, “’ 9
n =
1,2, .. . .
(8)
186
J. Aboudi,
The dynamic stresses induced by the propagation of skew cracks
By replacing all the derivatives in (4) as described before, the following system of algebraic equations is obtained, which is of second-order accuracy (i.e. the truncation error is second order in the space and time increments): +
ui, j1
+ ~3 [Ui_l,j,n-1
-
‘3ui--1 I j t n+l
+ Gel
‘*
n+l
L”i+l, j+l,n
+j(t)flAt[Ui+l,j,n
-
f3Ui+l
ui+l,j,
n-1
+ ui-l,
, j , n+l
=
2[ 1 - AC:- BE;lUi,j,n -
‘i 1j ?n-l
1 +AE: [Ui+l,j,n+ Ui-l,j,n 1 + BE; [Ui,j+l,n+ Ui,j-l,n 1
j-l,n
-
ui-l,
j+l,n
-
ui+l,
-Ui-,,j,nI/‘+D’~“[‘i,j+l,n
j-1,n
II4
- ui,j_l,n l/2 >
(9)
where ei = At/At,
e2 = At/A7 ,
e3 = elf(t)/2
.
According to (9) it is possible to compute the displacement vector at time t + At from its values at the previous two time steps by solving at each time step for every j a tridiagonal system, for which a direct inversion algorithm can be employed and no iterations are needed. The numerical formulation of the boundary conditions (7) at n = 0 is obtained by approximating the displacement gradients by their central difference expressions for derivatives in the t direction, and by backward and forward difference expressions for derivatives in the 77direction, when approaching 17= 0 from n < 0 and n > 0, respectively. Consequently, at every time step t = n At the following system of algebraic equations in the unknowns ~~,e,~ is obtained (with eq = An/2Ai):
(10) =
‘[“6uj I ,n +“8ui >~1 >.I + AVg,/P9 Tl
E4”,[Wi+l,0,n - Wi-l,O,nl’
ulOWi,O,n = f"lOWi,Tl,n
+ 'qg3/p
*
The upper and lower signs in (10) correspond to the regions r) < 0 and n > 0, respectively, and the index i runs over all points along the &axis where the boundary conditions (6) are prescribed. In [4] the algebraic systems obtained from the numerical treatment of the boundary conditions were solved very easily by utilizing the Gauss-Seidel iterative procedure. The present case turns out to be more involved and the Gauss-Seidel method does not converge. Indeed, it can be verified that the sufficient conditions, which were shown to bp satisfied in the problem treated in [4], are not satisfied in (10). Accordingly, we solved the algebraic systems (10) by a direct method, utilizing the band structure of the matrices of coefficients of the unknowns. In order to apply the proposed difference scheme (9), it is necessary to show that it is stable (i.e. the difference between the exact and numerical solution of the difference equations (9) re-
J. Aboudi, The dynamic stresses induced by the propagation of skew cracks
187
mains bounded as the number of time steps increases); also an estimation is needed for the stability criteria. The stability analysis of (9) proceeds essentially as in [4] for the difference equations in the case when t;r(x) = cons& whose form is similar to (9). However, in the present case B and C are not constant matrices, and in (9) there exists the additional term involving the matrix D which is associated with the term au/an in (4). This term does not affect the stability criteria practically since it is known [ 111 that if D is a bounded operator, then the stability is not destroyed by the small perturbation which is represented by the term AfD in (9) (D is bounded since it was assumed that H(x) f 0). Consequently, the stability conditions appear according to [4] in the form e2 < l/max v, , r
r=l,2,3,
(11)
where E = cr = e2 (which is the usual case in practical computations), values of the real symmetric matrix M=Asin2p+Bsin2q+Csinp
and v, are the three eigen-
sinq cosp cosq,
(12)
with sin p, sin q varying independently between -1 and 1. We can estimate v, in terms of the spectral radii of the matrices A, B, C as follows: lMl
G ll All + IlBll + II CII ,
(13)
where II - II denotes the spectral norm. Since M is symmetric, its spectral radius R(M) is equal to its norm IIMII,so that R(M)dR(A)+R(B)+R(C),
(14)
since A, B, C are symmetric as well. The spectral radii R(A), R(B), Z?(C) can be easily computed from the characteristic equations of A, B and C. Since these matrices depend on the location x, conditions (11) must be satisfied for all the relevant range of variation of x. All the results given in this paper were obtained with the space increments AE = Aq = d/SO.
3. SH wave propagation in an impacted elastic wedge As a check to the proposed method of solution we apply it to the problem of an elastic wedge subjected to spatially uniform but time-dependent shear tractions which are applied to one of its surfaces while keeping the other one free. An analytical solution to this problem is available (Achenbach [lo]) that it can be employed to assess the present numerical method. Consider an elastic wedge occupying the region r > 0,O < 8 < cwwhose interior angle is CLThe wedge is at rest prior to t = 0, and at t = 0 time-dependent shear tractions are applied to one surface while keeping the other free. Thus the boundary conditions are 00
uez
=
I
0,
W) >
r>O,
e=o,
r-20,
e=(ll,
(15)
J. Aboudi,
188
7%e dynamic stresses induced by the propagation of skew cracks
with U(t) the Heaviside step function and u0 an amplitude The analytical solution to this problem is given by a0 cos 0 U(t - r sin
~ez(r, 8, t) =
L
e/c,) + C (r, 8, t)
factor.
,
0 < 0 < n/2 ) 0 2 K/2 )
(r, 6, t)
where (with k = cr/a > l/2, 0 = -cash-‘(c,t/r)) c
(r, 8, t) = -c2uo
sin(n/2k)/(nkr)U(t
- r/cz) j
FI(r/cT,
0) dr ,
r/c2
F,
(r/c* T, 0) = sinh(P/k)
sin(8/k)/{[cos(7r/2k)
- cosh(P/k)
cos(e/k)l*
For this problem the proposed numerical treatment can be employed to (3) with f(t) = 0) the boundaries of the wedge defined by
H(x)
-d,
.Y >
-d + .Y tan Q ,
XGO,
0
+ [sinh(P/k) by transforming
sin(e/k)l*}
.
(according
)
(17)
=
to the surface q = 0 of the half-plane n < 0 on which the interior region of the wedge is mapped. In fig. 2 the numerical and analytical (16) solutions for the shear stresses uBz are shown versus r/d along the ray 8 = 7r/2 when c2 t/d = 0.5 and 1 for (Y= 5x/4. The numerical solution was obtained by applying (9) with g, = g, = 0 in (lo), whereas g, was taken from (15). It is clearly seen that a reasonable agreement between the two solutions exists. Whereas the analytical solution predicts singular stresses near the vertex of the wedge (r = 0), the numerical solution yields a finite value there as can be expected. In previous numerical solutions for dynamic crack problems [ 3,4]
+-+olJ--’ 0
I
-
r/d
0
Fig. 2. Numerical solution (solid lines) and analytical solution (16) (dashed 0 = n/2. in a wedge under anti-plane loading when czt/d = 0.5 and 1.
-
0.5 r/d
-+----I
lines) for the shear stresses agZ versus r/d along
J. Aboudi, The dynamic stressesinduced by the propagation of skew cracks
189
the numerically obtained finite values near the tip of the crack were shown to correspond to the asymptotic solution near the tip of the crack and can therefore be employed to extract the stress intensity factors. Let us explore therefore these finite values for the present problem of the elastic wedge. These values are u~~/u,-,= 1.89 and 2.38, which are obtained when c2 t/d = 0.5 and 1, respectively, at the point Y= At, 8 = 7r/2 and correspond to the analytical solution (16) at Y= rat, 8 = n/2, where y = l/2. Consequently, they can be interpreted as the asymptotic solution of the wedge problem at a distance Y= yAl from the vertex. In the crack problems discussed in the sequel this type of information will be useful in extracting the associate dynamic stress intensity factors.
4. Skew crack propagation
caused by the incidence
of SH waves
Consider a homogeneous elastic medium containing a semi-infinite crack which lies along the ray 8 = (Y(fig. 1b). At time t = 0 a plane incident SH wave strikes the tip of the semi-infinite crack. The incident wave has the form
wqx,y,
= W,TU(T)
)
(18)
with
where 8, defines the orientation of the front of the plane wave, and w0 is an amplitude factor. The incident stress wave generated by (18) is given by
~$2= (TV cog8 - e,) u(7) ,
(19)
where 00 =
w&,
’
It is assumed that the crack starts to propagate at the instant t = 0 with a velocity f(t) < c2 along the ray 8 = 0. Therefore, in contrast to previous works on crack propagation in their own plane it is assumed here that the crack deviates from its original plane with an angle K - (Y.For this problem an asymptotic expression for the stress field in the vicinity of the tip of the crack - in the case of constant rate of extension f(t) = u - was given by Achenbach and Varatharajulu [ 7 I. This expression has the form u&2,$,
t) -
(t/R)“* k/W, 4) >
(20)
where (R, #> is a new set of polar coordinates centered at the crack tip (fig. lb), k, is given by k, =
(1 - u~/c;)~‘~ K/u”* ,
(21)
and @(u, $J)describes the universal dependence of the stress field on the polar angle 4 in the vicinity of the crack. The explicit form of Q(u, 4) appears in [ 71 and [ 121. When $J= 0, @(u, 0) = 2l’*. In [ 71, the function K in (21) is plotted versus the skew angle for two values of the incidence angle eo.
190
J. Aboudi,
The dynamic stresses induced by the propagation of skew cracks
In the special case when (Y= K and there is a stationary crack (u = 0), the asymptotic expression for the diffracted stress (20) reduces to
2% (1 + sin 0e)“* (c2 t/R)“* cos(@/2) .
Q,(R, 4, t) = uoz(r,8, t) - 7
(22)
We can utilize the above asymptotic expressions in order to check the proposed numerical method of solution in the present circumstances as well as the associated dynamic stress intensity factor which is extracted from this numerical solution. The skew crack problem can be handled according to the present method by transforming the surfaces of the composite crack to the planes n = TO in the plane (t, q) to which the original domain (x, y) is transformed according to (3) with f(r) representing the velocity of the crack. The function H(x) of (3) appropriate to the skew crack is given by d - x tan a ,
XGO,
d,
x20.
H(x) =
(23)
In the transform plane the tip of the porpagating crack appears always at the origin of coordinates t = 0, 7)= 0. Since the stress intensity factor is determined from the singular behavior of the stress field in the vicinity of the crack, and the latter is contributed from the diffracted part of the field (not including the incident wave), it would be sufficient to construct the solution of an initially undisturbed medium with the boundary conditions (6) in which g, = g, = 0 and
g&O
= a’e’!,
tgo,
q=o,
tao,
(24)
provided that the independent variables of ~$2 are properly transformed. It is instructive to consider first the problem of the diffraction of an incident SH plane wave by a semifinite stationary crack lying along x < 0, y = d (i.e. with a = 1~in fig. 1b) in an infinite elastic medium. This problem possesses an analytical solution due to De Hoop and reviewed in [ 131. In the case when the plane wave front is parallel to the surface of the crack (i.e. 0, = CY = K), this solution has the form (see [ lo])
qJR> $, z) = q&
0, z) = JW - R/c,) >
-n/2 < @< a/2 )
(25)
where (2R/cz)l’* 71J/20, = (t - R/cZ)l’* { [ 1 - J, tan-‘(1 /J,)] cos(#/2 - 7r/4) - 11 - J, tan-’ (l/J,)]
(26)
sin(@/2 - a/4)) ,
with Jf = R(1 - sin @)/[c,(t - R/c,)]
and.
Ji = R( 1 + sin 4)/[cl(t
- R/c,)]
.
Accordingly, the asymptotic expression for u+ near the tip of the crack has the form
J. Aboudi,
The dyrumic
191
stresses induced by the propagation of skew cracks
uaZ(R, @,Z)-2”2 o0 t [(C,t/R)“2 - n(1 - sin $~)I’~/21cos($/2 - 7r/4) - [(c,t/R)“’
- n(1 + sin @1)“~/2]sin(@/2 - 71/4)}/7r ,
(27)
which reduces to (22) when 8, = 71and R/c, t e 1. Let us define the dynamic stress intensity factors associated with the three modes of motion of the crack in the form Ti (t) = U,j(R, 0, t)R”’
(28)
)
for small R such that i = 1 when j = #, i = 2 when j = R and i = 3 when j = z. According to (28) the dynamic stress intensity factor derived from (27) has the form T, (t) = 20, (cz t)“2 /T ,
(29)
which can also be obtained as a special case of the skew crack problem with (Y= 0, = K and u = 0. In fig. 3 the numerical and analytical solutions for uez are shown when c2 t/d = 0.5 and 1 versus R/d along the ray 4 = 0. It is clearly seen that a good correspondence between the solutions is ob-
II
il
C21/d ao.5
c,t/d=l.o
01
0
(b)
0.5
cpt/d-
I
I
Fig. 3. (a) Numerical solution (solid lines) and analytical solution (25) (dashed lines) for the shear stresses oez versus R/d along @= 0 for a stationary crack with 01= n under anti-plane loading when czt/d = 0.5 and 1, (b) the dynamic stress intensity factor curves based on the numerical solution (solid lines) and analytical solution (29) (dashed lines) for the stationary crack.
192
J. Aboudi. 77~ dynamic stresses induced by the ~ro~g~ti~~
of skewcracks
tained. In the vicinity of the tip of the crack the analyticaf solution (25) predicts a singularity at R = 0, whereas the numerical solution yields a finite value for the stress at the grid point (4‘= A& 7 = 0), which is the closest grid point ahead of the tip of the crack. For example, we obtain at this point the values of ~~~/u~ = 3.23 and 4.98 when c,tld = 0.5 and I, respectively, whereas the analytical solution (25) yields at the same point the corresponding values 2.25 and 3.54 (the asymptotic expression (27) yields 2.18 and 3.50). In order to explore these values, we compared them with the analytical solution (25) (or the asymptotic expression (27)) evaluated at several locations in the vicinity of the tip of the crack at the same times and found that they correspond to the exact solution at t = ?A.$, where y = l/2. The numerical solution at the point (4:= Ak, n = 0) can be employed to get the dynamic stress intensity factor of the crack. According to (28) we compute the values of T,(t) by multiplying uaz(t = At, n = 0, t) by the square root oft = rA[ with y = l/2. However, since the difference between the numerical solution for T3(t) obtained by (28) and the analytical one given by (29) is of order of R1’2 (see (27)), we implement this when computing T3(t) by employing (28) and adding f~A$)l’~ to the result. Fig. 3 shows curves for T,(t) based on the numerical and analytical (29) solutions, and good agreement between them can be well observed. We turn presently to the main problem of skew crack propagation. In fig. 4 the numerical solu= 0 is shown for the original crack orientation tion for the stresses uot versus R/d along the ray r#~ c$/d = I.0
c,ttd 105
O_ -
0.5 R/d
i
4 I
0 (0)
/ -
0.5 R/d
I
v/cz=o 2
-
0.5 c,tld
I
(b) Fig. 4. (a) Numerical solution for the shear stresses 0~~ versus R/d along Q,= 0 for a skew crack under anti-plane loading with IX= 3n/4, propagating with the constant velocities u/c2 = 0.2 and 0.5 when czt/d = 0.5 and 1, (b) the corresponding dynamic stress intensity factor curves based on the numerical solution (solid lines) and analytical solution (dashed lines).
193
J. Aboudi, The dy~~~c stresses induced by the propQgQ?~o~ ofskew cracks
01= 3n/4 (fig. 1b) when c,t/d = 0.5 and 1 for the two cases in which the propagation velocities are constant and given by u/c, = 0.2 and 0.5. The dynamic stress intensity factor curves T3(t) shown in the figure were extracted from the numerical solution for the stress uoz obtained at the closest grid point ahead of the crack tip (4‘= Al, q = 0) as explained before. These curves are compared with stress intensity factors based on the analytical treatment in [7] and were obtained by applying (28) to (20). It is well observed that a satisfactory agreement exists between the two methods of solution. It should be mentioned that the present choice of OL= 37r/4 represents a fairly significant deviation of the crack from its original plane since in practical situations the skew angle 77- (Y does not seem to exceed the value of ~16 [ 121 so that even a better agreement between the two methods of solution would be accomplished.
5. Skew crack propagation caused by the sudden application of normal tractions on its surface We consider presently the analytically yet unsolved problem of the skew crack propagation caused by a sudden application at t = 0 of equal and opposite uniform normal tractions to the surfaces of a semi-infinite crack. The crack lies originally along the ray 0 = a! (see fig. 1b) in an initially undisturbed medium. The applied tractions produce plane waves and diffracted cylindrical waves whose center is located at the original crack tip. It is assumed that the semi-infinite crack extends at a velocity f(t) out of its plane at the instant of the application of the tractions. This problem is solved according to the proposed method of solution by using (3) to transform the surfaces of the crack to the planes Q = TO in the transform domain as described previously. The appropriate function H(x) is given by (23). The boundary conditions are given by (6), where
g,=g,=o,
g&5 f) = -u(J
[
1- U(x)1WI ,
l
q=“o,
(30)
with x = .$+f(t). According to (30) the normal tractions are applied on the surfaces of the original crack only, while keeping the newly produced surfaces of the secondary crack free of tractions. Before considering the skew crack problem it is instructive here too to treat the case of a collinearly propagating semiin~nite crack in an infinite solid which was treated by Baker [ 141 and for which some analytical results are known. The boundary conditions are given by (6) with g, = g, =Oandg, = -u,, U(t). The asymptotic expression for uGq,,which is valid near the tip of the crack, is given by [ 141: (31) where (R, 4) is the polar coordinate system centered at the crack tip. In (3 1) F+(O) is a complicated expression which depends on the material parameters and the uniform velocity u of the crack, and \k(u, $) represents the universal dependence of the singular stresses on the polar angle (bin the vicinity of the crack, and the propagation velocity u. For Qj= 0 (i.e. for observation points along the path of the crack) *(u, 0) = 1, whereas \k(O, #) = cos3 ($12). The explicit form of U($J, u) can be identified in [ 121. According to (28) the analytical expression of the dynamic stress intensity factor T, (t) is given by T,(t) = 20&l
- u/c,)CZt]1’2/aF+(0).
(32)
194
J. Aboudi, The dynamic stresses indrtced by the Fmpagation ofskew cracks
MODE II
MODE I
4-
r&c, =0.542
3
3-
c,t/d ~1.0
c,t/d = 1.0 / b” ’ 2z / b
0.5
I
-
0.3 R/d
r
0.5
I.t::::.:--‘-: 0
(a)
0.5c&/d
Fig. 5. (a) Numerical solution for the stresses ow versus R/d along # = 0 for a stationary crack with a = n under in-plane normal loading when clt/d = 0.5 and 1 (the dynamic stress intensity factor curve based on the numerical solution (solid lines) is compared with the analytical one based on (32)), (b) numerical solution for the stresses 00~ versus R/d along @J= 0 for a stationary crack with o = n under in-plane shear loading when clt/d = 0.5 and 1, and the corresponding dynamic stress intensity curve (the material constants in this case have the value ~/et = 0.542).
Unless otherwise mentioned, all the results given in this section correspond to the choice CJC, = I/2. In fig. 5 the stresses ucbe for a stationary crack (v = 0) lying along the ray a = K are shown versus R/d along Cz, = 0 when c,t/d = 0.5 and 1. The associated numerical stress intensity factor curve is extracted from the values of u*@at the closest grid point ahead of the tip of the crack, i.e. at (t = A$, n = 0), by applying (28) with R = ?A$ and y = l/2 (since here too the numerical values of o++ at this point correspond to the values of crbQ,obtained from the asymptotic expression (3 1)). This curve is compared in the figure with the corresponding one based on (32). It is seen that a very satisfactory correspondence between the two methods exists. As we need in the sequel to extract the dynamic stress intensity factor T,(t) (associated with a mode II type of motion) from the shear stresses ueR induced near the tip of the crack, let us apply the boundary conditions (6) with g, =g, = 0 and the shear impact g, = -u0 u(t) in order to produce a mode II type of opening. In this specific case we chose the elastic material with cJc, = 0.542 (steel) since for this value Sih, Embley and Ravera [ 151 present the corresponding stress intensity factor appropriate to this problem. In fig. (5) the stresses CJ+~versus R/d along 4 = 0 are shown for a stationary crack (u = 0, (Y= n) when et t/d = 0.5 and 1. The stress intensity factor r,(t) is extracted from the numerical
J. Aboudi, The dynamic stressesinduced by the propagation of skew cracks
195
solution of utiR induced by the vicinity of the crack tip according to (28) and adding the contribution (~a#‘~, y = l/2, as in the previous anti-plane case (see also [4]). The resulting values of T2(f) as exhibited in fig. 5 are in good correspondence with the curve shown in [ 151, which is based on an analytical treatment of the problem and a numerical inversion of the Laplace transforms for the evaluation of the stress intensity factor. For example, when c,t/d = 1, our numerical whereas we obtain the value 0.55 from the curve in solution yields the value 0.52 for T,(t)/u,d”2,
[151. It is still necessary to check the reliability of the numerical scheme (9)~(10) in the in-plane loading condition of a crack when the transformation (3) is involved since in the previous cases (exhibited in fig. 5) H(x) was identically constant so that this transformation was not actually active. Since no analytical results have been obtained for the skew crack problem under in-plane loading, we chose to perform such a check by considering the previous stationary crack problem with sudden normal tractions acting on its boundary but lying now along the ray cy= 3n/4 (rather than (Y= rr as in the previous cases). The stress intensity factor T,(t) which is extracted from the closest grid point ahead of the tip of the crack, i.e. from the point (t = A& 77= 0), can be compared with the analytical results given by (32). However, since the location of this point corresponds in the present situation to C#= n/4, we must implement (28) by taking into account the fact that the numerical solution for u6@obtained at this point corresponds to the asymptotic expression (3 1) in which \k(O, 4) = cos3($/2) = l/23’2, R = rag, y = l/2. Accordingly, the numerical values of T,(t) should be extracted from uGOafter dividing the latter by this value of \k (in order to recover the value of uQOat 4 = 0) and employing (28). In fig. 6 the induced stresses uGQ,by the stationary crack, are shown versus R/d along $ = a/4 when c,t/d = 0.5 and 1. In this figure the numerical values of T,(t) are shown and contrasted with the analytical values based on (32). It is clearly seen that the method of solution is reliable in the present situation too. Having established the reliability of the numerical treatment and the extraction of the dynamic stress intensity factors in mode I and II types of motion of the crack, it is possible to investigate
015 -
R/d
(a)
-
c,t /d
(b)
Fig. 6. (a) Numerical solution for the stresses om versus R/d along @= 0 for a stationary crack with Q = 3n/4 under in-plane normal loading when ct t/d = 0.5 and 1, (b) the dynamic stress intensity factor curve based on the numerical solution (solid lines) compared with the analytical one based on (32).
196
J. Aboudi, Eke dynamic stresses induced by the p~pagat~on ofskew cracks
c,i/d =0.5
c,t/d
I 8
0.5 -
0
-
I I
b
I
0
R/d
0.5 c,t/-d
~1.0
I 0.5 --R/d
-
05 c, t/d
i I
I
Fig. 7. Numerical solution for the stresses DQ,++, 06~ versus R/d along 0 = 0 for a skew crack with a = 170’ under in-plane loading (30) when ~1 t/d = 0.5 and 1 and the corresponding dynamic stress intensity factor curves in modes I and II (the secondary crack propagates with the velocities u/cl = 0.1 and 0.2).
now the main problem of a skew crack propagation under in-plane loading (30). As was mentioned before, even though the primary crack is under pure mode I loading, we obtain in general a mixed mode problem in which both modes I and II are present. In figs. 7-9 we present the normal and shear stress uc4@,uQR versus R/d computed along the ray C#= 0 for a crack which was initially oriented along (I!= 170”, 150” and 135” when et t/d = 0.5 and 1. The propagation velocities of the secondary crack are taken to be uniform with the values u/c, = 0.1 and 0.2. The associated stress intensity factors T,(t), T,(t) appropriate to modes I and II, respectively, are also shown in the figures. The stress intensity factor 1”,(t) is extracted from the values of uo9 at the grid point ahead of the tip, whereas T,(t) is extracted from uGR at that point. It is clearly observed that the stress intensity factors decrease as the propagating velocity increases. On the other hand, T,(t) increase as the skew angle ar increases, whereas T,(r) decreases as cxincreases. These results are not unexpected since in the case of a collinear crack propagation, the shearing stress uoR vanishes along the ray 4 = 0 so that r,(t) = 0 and the mixed mode propagation becomes a pure mode I type.
6. Conclusions A numerical method of solution is presented which is able to handle the problem of the dynamic curved crack propagation under anti-plane and in-plane loading. The method is based on a certain
J. Aboudi, Tire dynarni~stressesinduced by the p~pagation ofskew cracks
C, t/d
I
I 1
-
0
0.5
c,t/d +I.0
=0.5
!
I
I
R/d
0.5 -
197
I
c,t/d
Fig. 8. Same as fig. 7 but with the original crack orientation
01 =
1 SO”.
time-dependent transformation which maps the physical plane of the crack into an auxiliary plane in which the crack appears propagating col~nearly. The reliability of the method of solution is checked by comparisons with some analytical results. Several extensions are possible. For example, the method of solution is formulated such that it can handle nonuniform skew crack propagation. In this paper we assumed that the original crack changes its direction abruptly. Common experience indicates that a crack will generally follow a smooth curved path. Furthermore, it has been shown recently by Willis [ 161 that a straight crack subjected to anti-plane static loading cannot abruptly change its direction of propagation even if it is loaded asymmet~c~ly but may curve smoothly. In such cases in which the crack propagates along a smooth curved path the present numerical method of solution will yield even more accurate results since, as a rule, numerical schemes are less accurate when dealing with functions possessing piecewise continuous derivatives (as in the case of a skew crack). On the other hand, analytical treatments for cracks that curve in an arbitrary fashion are lacking at present, so that numerical methods developed for such problems seem to be valuable. In order to explain some particular phenomena in fracture mechanics whose purpose is to characterize the actual behavior of cracks by using a simple mathematical crack model, some approaches for making the crack model more precise have been proposed. In [ 171 those approaches
19%
i
$=, 8 b 0 I
-
-
0.5 R/d
c,tM
Fig. 9. Same as fig. ‘7 but with the original crack orientation o = 135”.
are divided into six different classes of problems. One of those is the class of nonlinear shaped cracks to which the skew and curved types of cracks belong. This type of crack configurations is encountered, for example, when a mixed mode propagation of a crack is involved since generally the crack does not fottow a straight path but a curved or a bent path. Several theoreticaf and experimental evidents af particular nonlinear shaped cracks are discussed in [ 1i’] and in the references cited there. Consequently, the present paper praposes a numerical method of solution applicable to several types of dynamic crack propagation which belong to the general category of nonlinear crack configurations.
Acknowledgment This research was performed under contract AFOSR-‘76-3014 of the U.S. Air Force.
References [I ] J.D. Achenbach, The dynamic effects in brittle fracture, in: S. Nemat-Nasser (ed.) Mechanics today, 1 (Pergamon Press, 1974).
J. Aboudi, The dynamic stressesinduced by the propagation of skew cracks
199
[2] L.B. Freund, The analysis of elastodynamic crack tip stress fields, in: S. Nemat-Nasser (ed.), Mechanics today, 3 (Pergamon Press, 1976). [ 31 J. Aboudi, Numerical solution of dynamic stresses induced by moving cracks, Comp. Meths. Appl. Mech. Eng. 9 (1976) 301-316. [4] J. Aboudi, The dynamic stresses induced by moving interfacial cracks, Comp. Meths. Appl. Mech. Eng. 10 (1977) 303-323. 15) R.J. Nuismer, An energy release rate criterion for mixed mode fracture, Int. J. Fracture 11 (1975) 245-250. [6] G.C. Sih, Strain-energy-density factor applied to mixed mode crack problems, Int. J. Fracture 10 (1974) 305-321. [7 ] J.D. Achenbach and V.K. Varatharajulu, Skew crack propagation at the diffraction of a transient stress wave, Q. Appl. Math. 32 (1974) 123-135. [ 81 B.V. Kostrov, Diffraction of a plane wave by a smooth rigid wedge in an unbounded elastic medium in absence of friction, J. Appl. Math. Mech. 30 (1966) 244-250. [9] S.H. Zemell, Diffraction of elastic waves by a rigidsmooth wedge, SIAM J. Appl. Math. 29 (1975) 582-596. [lo] J.D. Achenbach, Shear waves in an elastic wedge, Int. J. Solids Structs. 6 (1970) 379-388. [ 111 R.D. Richtmeyer and K.W. Morton, Difference methods for initial-value problems, 2nd ed. (Interscience, 1967). (121 J.D. Achenbach, Wave propagation, elastodynamic stress singularities, and fracture, in: W.T. Koiter (ed.) IUTAM Congress, Delft (1976) 71-87. [ 131 J.D. Achenbach, Wave propagation in elastic solids (North-Holland, 1973). (141 B.R. Baker, Dynamic stresses created by a moving crack, J. Appl. Mech. 29 (1962) 449-458. [IS] G.C. Sih, G.T. Embley and R.S. Ravera, Impact response of a finite crack in plane extension, Int. J. Solids Structs 8 (1972) 977-993. [ 161 J.R. Willis, A discussion of crack-forking in anti-plane strain deformation, Int. J. Fracture 11 (1975) 489-493. [17] H. Kitagawa, R. Yuuki and T. Ohira, Crack-morphological aspects in fracture mechanics, Eng. Fracture Mech. 7 (1975) 515-529.