- Email: [email protected]
Analysis of two dimensional fracture problems with multiple cracks under mixed boundary conditions
Engineering Fraerure Mechanics Vol. 34, No. 4, Printed in Great Britain.
pp.921-934 1989
~1~7~/89
$3.00 -t 0.00
Pergamon Press pk.
ANALYSIS OF TWO DIMENSIONAL FRACTURE PROBLEMS WITH MULTIPLE CRACKS UNDER MIXED BOUNDARY CONDITIONS WEN-HWA
CHEN and CHENQ-SHYOUNG
CHANG
Department of Power Mechanical Engineering, National Tsing Hua University, Hsinchu, Taiwan 30043, Republic of China Abstract-This paper presents an efficient finite element alternating method for the analysis of two dimensional mixed-mode fracture problems with multiple cracks under mixed boundary conditions. Based on the analytical solution derived for an unbounded crack body with a central crack under arbitrary crack-face tractions, the resultant residual displacements/stresses on external boundaries are repeatedly evaluated till the variation of stress intensity factors and resultant residual stresses on each crack are negligible. Excellent agreement between the computed results and referenced solutions is observed. The deformed profiIes and m~mum shear stress contours are also presented. The economic computer time used and simple finite element idealization needed show the advantages of the method.
IN DEALING with engineering fracture mechanics problems, many numerical methods have heen developed and employed because of the complicated boundary conditions. Among those numerical methods, the finite element method has received much attention by many researchers. In this method, simple quarter-point isoparametric element@, 21or rigorous hybrid singular elements[3-6] are usually adopted to model the region near the crack-tip. Numerical experiments indicate that the accuracy of the stress intensity factors obtained by these techniques largely depends on the mesh size taken and/or the computer time consumed. To compensate those shortcomings, especially for solving realistic fracture problems with multiple cracks, a finite element alternating method which involves the iterative superposition of finite element solution of an uncracked structure and analytical solution of an abounded crack body subjected, to arbitrary normal and shear loadings on crack surfaces is then studied. Vijayakumar and Atluri[7] derived the complete analytical solution of a three dimensional body with an elliptical crack under arbitrary crack-face tractions. The solution later was successfully implemented in conjunction with the finite element method to compute mode I stress intensity factors of an elliptical crack[8] and multiple cracks[9]. Unfortunately, since the elliptical axes taken in Vijayakumar and Athni[7] are finite, the solutions obtained cannot be directly applied to solve two dimensional problems. Besides, those studies as mentioned above are only limited to the problems under force boundary conditions. Recently, based on the solution obtained by Sneddon and Elliott[lO] for the problem with crack surface under symmetric normal pressures, the two dimensional analytical solution of an infinite plate with a crack under arbitrary normal and shear loadings together with an efficient fimte element alternating t~hnique have been successfully developed by the authors[ 11] to deal with mixed-mode fracture problems with multiple cracks. Similarly, however, only force boundary conditions are considered. Hence, for more general engineering applications, the extension of the previous work[l I] to analyse fracture problems under displacement/force boundary conditions becomes the main objective of this work. To satisfy the prescribed mixed boundary conditions, based on the analytical solution derived, the resultant residual nodal displa~m~ts and resultant residual equivalent nodal forces are repeatedly computed during the iterative sohrtion process until the variation of stress intensity factors and resultant residual stresses on each crack converge to a small value. The deformed profiles and maximum shear stress contours are then presented. The high efficiency and accuracy of the technique developed can be noted. 921
WEN-HWA CHEN and CHENQ-SHYOUNG CHANG
922
ANALYTICAL SOLUTIONS
As shown in Fig. l(a), a two dimensional infinite plate with a crack at x = 0 and - 4 d y < a, where a is half of the crack length, is considered. The crack surfaces are subjected to arbitrary normal and shear loadings which are represented by a polynomial of any orders. The complete analytical solution can be superposed by two analytical solutions of normal and shear loadings, respectively. For different loading considerations on crack surfaces, the boundary conditions are to be:
(0 stresses ((i=, cryy,r xY) and displa~ments (u,, u,,)-*Oas x2 + y*-+co (ii) for normal loading case (see Fig. lb): atx-0
(iia) rxY= 0 (iib) a,, = - i @(y/a)”
atx=Oand[yl
n-0
and (iic) U*= 0
atx=Oandlyl>a
(iii) for shear loading case (see Fig. lc): (iiia) a, = 0
atx==O
(iiib) TV,, = -
atx=Oand/y(,
and (iiic) f+ = 0
atx==Oand/y(>a.
Fig. I. An infinite plate with a crack under arbitrary crack surface tractions. (a) Arbitrary crack surface tractions. (b) Arbitrary normal crack surface tractions, (c) Arbitrary shear crack surface tractions.
Analysis of 2D fracture problems
923
In the above conditions, N is an arbitrary positive integer and C$ and C!,?) are polynomial coefficients. From the symmetry and anti-symmetry of loading cases with respect to y-axis, only the semi-infinite elastic medium for x 2 0 needs to be solved. Based on the plane strain assumptions and Sneddon and Elliott’s solution[lO] for symmetric normal loading, neglecting the effect of body forces, the analytical solutions of stresses and displacements are derived as follows: (a) for normal loading case: 0 xx=--
3 j:
* 0W@W)cos(W Oyy= - ?r s zxy
=
-
24, =
7
2
&)eetXd<
(1)
- Sx)e-[“d5
(2)
@‘,‘)(<)cos(~y) + &“(<)sin(
{:
+ &‘Wsin(WNl
{31’YOsin(tv)
-
&‘)
>WO>>&
-{*
)cos(
om ($(‘)(t
dt
(3)
1 - v) + 5x]e + $5
(4)
- 2v - tx)e-‘x$
(5)
s and uy = - y
0m{$1’)(<)sin({y) - @)(<)cos([y)}(l s (b) for shear loading case: d XX=-- ‘,” 10a #PWsin(
- &2’(t>cos(~~)}5e-tx d5 - &2)Wcos(t~)l(*
* 0m@\2)(0sin(W Oyy= - IL s z XY=--
1
24,= 7
om
s
- We-”
{$\2’(Ocos(5~) + $i2)(Wn(5y))U
Joa
(6) dt
(7)
-
(8)
{$\2)(5)sin(
(9)
and uv = - y
1: ($\2)(C;)cos(~y) + $P)(<)sin(ry))[
- 2(1 - v)]e+
7
(10)
where E is Young’s modulus and v is Poisson’s ratio. Equations (l)-(lO) have satisfied the boundary conditions (i), (iia) and (iiia) for appropriate respective cases. Further, after satisfying other boundary conditions (iib), (iic), (iiib) and (iiic), the functions $\‘)(<) and $p)([) can be determined by setting the functions #j(t) and $i2)(<) as zero and the polynomial power n in the boundary conditions (iib) and (iiib) as even. Similarly, #j(e) and $J$~)(<)are determined by setting $\‘)(r) and $p)({) as zero and the polynomial power n as odd. Define the variables
g:)(q)
= a(z/2q)“2
i
C$%f”;
It-0 wPtl)
sin@rl)
=
(wtl
W2J-
=
(7~
/2)‘/2Jl,2
+o (P)
’
!z a
= p “2F,p’(p)
I,2 (PI);
6~);
(11)
924
WEN-HWA CHEN and CHENQ-SHYOUNG
CHANG
where p is equal to 1 and 2. for normal and shear loading cases, respectively. The value q = 1 denotes the loading case when the power n is even and q = 2-represents the loading case when the power n is odd. J_,,2(plf) and J,,*(ptf) are the Bessel functions of the first kind. Two sets of dual integrals for the cases of even n and odd n are then obtained as follows:
and
From the solutions of dual integrals as shown in [12], the solutions @‘)(p) and ~~(~) of the above dual integrals can be obtained. Then, from eq. (1 I), the functions #“<<) and &p)(t) are determined to be
and
where r(n) denotes the Gamma function. Substituting the ,expressions of &p)(t) and @‘)(t) in eqs (12) and (13) into eqs (l)-(lO), the closed-form solutions for both, of normal and shear .loading cases are thus found. Through integration and appropriate superposition, the complete analytical solution of stresses and displacements can be rearranged and summarized in a matrix .form as W = WI VI
(14)
and +> =
WIWI
WI
where [IV] and [U] are the matrix of relation function; The stress vector .{a], displacement vector (u} and the coefficient vector {C> of the polynomials are
(16)
Analysis of 2D fractureproblems
From the definitions[l3],
925
the stress intensity factor vector {K} are thus determined by WI =
PI w
(17)
where [S] represents the relation matrix. Consequently, once the coefficient vector {C} in eq. (16) is found, the stress distributions and deformatians at any locations in the infinite plate and the stress intensity factors at two crack tips can be computed from eqs (14), (15) and (17), respectively.
FINITE
ELEMENT
ALTERNATING
PROCEDURES
The finite element alternating procedures in solving the fracture problem with multiple cracks in a bounded plate subjected to mixed boundary conditions can be stated as follows: (1) Solve the bounded uncrack plate under given external loadings and displacements as the problem considered and evaluate the finite element solution of residual stresses a{ at the location of fictitious crack k. The subscript k is the number of the crack (k = 1,2, . . . , K). (2) Since the crack surfaces are stress free, reverse the sign of the residual stresses u/k and fit them with polynomials of appropriate orders. The fitted stresses ai can be expressed as ai = {L}‘{C}R where (C} E is the polynomial coefficient vector of crack k and {L} is the vector composed of @/Q)” (n = 0, 1,2, . . .). The value a, is half of the length of kth crack. For better accuracy, {C}: is determined by the least-square method through the minimization of the integral, 4
(o-l;-a:)2dS
(k=1,2
,...,
K).
s Thus,
where [Hlk =
s a’
-*k
and
{R>k=
I”-ak {L} af dS.
(3) From (C)i and eq. (14), the reversed residual interaction coefficient vectors of interaction effects induced by crack k at fictitious crack I (1 = 1,2, . . . , k - 1, k + 1, . . . , K){A},k, can be evaluated. If I> k, add {A},k to (C)y and compute the updated (C)y, Thus, the coefficient vector (C)A of the residual stresses for crack k need to be released for next computation is found to be
Repeat all steps as mentioned earlier until the cumulative residual stresses on each crack are negligible. As a result, the resultant coefficient vector {C}, = (C)f + (C)i + {C}: + - * * for the current iteration is obtained. The corresponding resultant residual stresses and displacements can be computed by eq. (14) and eq. (15), respectively. (4) Substituting the resultant coefficient vectors {C}, into eq. (17), the stress intensity factors of each crack for the current iteration can be calculated and updated. (5) The stresses/displacements of the cracked plate for each iteration can be computed by the sum of finite element solutions and those solutions calculated by each resultant coefficient vectors {C}, through eqs (14) and (15). The final stresses/displacements of the cracked plate are then readily determined once the iteration is completed.
926
WEN-~A
CHEN and C~EN~S~O~G
CHANG
(6) Once the variation of stress intensity factors of current iteration for all crack tips and the ratio of resultant residual stresses/permissible stress for each crack are smaller than a small value Y (Y = O.@W, say, 1 @wmt
iteration)
_
j.pwious
iteration)
/< y
1zi&@m-m itention) J
and ” ~W’~G>,~ dS s -=k cy 2&a
forallk=1,2,...,K
the iteration is finished. Here, i is equal to I and II for mode I and mode II, respectively. (7) To satisfy the force boundary conditions on external boundaries of the fracture problem (from eq. 14), the resultant residual equivalent nodal forces {Q}: on external force boundaries at mth element are computed by the sum of the contribution of each crack (Q), (k = 1,2, . . . , AC) as
iQ>: = j, fQ3nlk where
and
PI, =
WL’IPl~ dS s %I
[G]& denotes the relation matrix and S, is the portion of the boundary of mth element where force boundary conditions applied. [N] represents the shape function, ET]&is the transfo~ation matrix, and [wk is the matrix of relation function for kth crack as defined in eq. (14). To satisfy the displacement boundary conditions on external boundaries (from eq. 15), the resultant residual nodal displacements {~}1+on external displacement boundaries at node r are calculated by the sum of the contribution of crack k (k = 1,2,. . . , K), (ujlk, and expressed as
where
and
KJh is the
matrix of relation function for kth crack as defined in eq. (15) which is estimated at node r. [Mlk is the transformation matrix taken for locating the kth crack in global coordinate system. (8) Reverse the resultant residual covalent nodal forces (Q)z and resultant residual nodal displacements (a>: calculated in procedures (7) and consider those as applied loadings and displacements acting on the external boundaries of the bounded uncrack plate. Repeat all the iteration procedures (1x7) until the variation of stress intensity factors and resultant residual stresses on each crack are negligible. RIEWLTS AND DIS~SSIONS
To demonstrate the validity of the technique developed in this work, a bounded plate with inclined cracks subjected to different force/displacement boundary conditions at upper or bottom boundaries is considered. As shown in Figs 2 and 3, the finite element meshes with 9 eight-node
927
Analysisof 2D fracture problems
Fig. 2. The rectangular plate with an inclined crack under mixed boundary conditions.
isoparametric elements and 80 degrees of freedom are adopted. All numerical computations performed by MicroVAX II computer. Example
are
I: the rectangular plate with an inclined crack
(a) Subjected to force condition at upper end and displacement condition at bottom end. As shown in Fig. 2, a rectangular plate with an inclined crack under uniform tension at y = H and fixed displacement in y-direction at y = -H is first analysed.
A’
r ---
2w
EI
ly=o
_--_--
uy==conat.
,P Pa
d
__-------
F=‘wv-
Fig. 3. The square plate with an inclined crack under displacement boundary conditions.
928
WEN-HWA
CHEN and CHEN~SHYOUN~
CHANG
Table 1. Normalized stress intensity factors for the rectangular plate with an inclined crack
4 B 30” 45” 60” 90”
41
Tada 1141
Annigeri WI
Present 9 ele,
0.73 1.488
0.386 0.706 1.01 1.446
0.363 0.685 1.066 I.406
24 ele.
Tada 1141
Annigeri WI
9 ele.
24 ele.
0.357 0.706 1.083 1.461
0.6 0.000
0.515 0.551 0.458 0.000
0.506 0.564 0.492 0.000
0.516 OS87 0.492 0.000
Present
Table 1 shows the values of normalized mode I and mode II stress intensity factors (4, I;;,) for various inclined angles 8. (Fi, F,,) = (&, ~~*)/~~. Good agreement between the computed results and referenced data by Tada et aL[14]and Annigeri and Cleary[lS] is observed. To evaluate the high efficiency of the present technique, another run with 24 elements is also performed. It is obvious that the solutions obtained have only little difference between the use of 9 elements and 24 elements. Hence, the simpler mesh of 9 elements is su~ciently good to obtain accurate stress intensity factors for this case. The excellent monotone convergence of normalized stress intensity factors for the cases of 0 = 45” and 6 = 90” is indicated in Table 2. Table 3 displays the computed normalized stress intensity factors with different orders of polynomial (N) taken for fitting the normal and shear stresses on the fictitious crack. For 8 = 90” case, since the normal stress on fictitious crack is symmetric with respect to the center of the crack, the polynomial order can be chosen as even only. It is seen that N = 2 is big enough to reach good solutions. (b) Subjected to ~~~~Iace~nt condition at both upper and bottom en&. To analyse the fracture problem under prescribed displacement boundary conditions using finite element alternating method (as shown in Fig. 31, a square plate with an inclined crack subjected to uniform ~splacements at upper and bottom boundaries is also studied. For B = 90” case, Fig. 4 displays the variations of normalized mode I stress intensity factor F$ with increasing values a/W. The results obtained by Isida[l6] are also compared. Figures 5 and 6 show the computed normalized stress intensity factors F, and F,, for the case of a/W = 0.4 vs various angles of 8 using the mesh of 9 elements and 24 elements, respectively. Both solutions are also Table 2. Convergence of normalized stress intensity factors for the rectangular plate with an inclined crack e =900
B =4S” Cycle bf iterations 1
i 3 4 5 Tada[Y4]
&
Tabie 3,
4 9 ele.
24 ele,
o.soo
o.soo
0.559 0.564 0.564 0.564 0.564
0.576 0.586 0.587 0.587 OS87
1.000 1.286 1.371 1.3% 1.404 1.406
1.000 1.322 1.420 1.450 I .459 1.461
24 ete.
9 ele.
0.500 0.623 0.664 0.679 0.683 0.685 0.73
OS00 0.626 0.676 0.695 0.703 0.706
Anniaerill51
41
24 ele.
9 ele.
1.488 1.446
0.6 0.551
0.706
Comparison of different orders of ~l~orni~s plate with an inclined crack
e=9oo
0 =4S” Order N 2 3 4 5 Tada[ 141 Annigeri[ 151
9 ele.
F,
24 eie.
9 ele.
0.706 0.706 0.706 0.706
:z
0.683 0.684 0.683 0.685 0.73 0.706
fitted for the rectangular
F,I
24 ele.
9 ele.
0.585 0.585 0.586 0.587
1.399 1.406 -
O:S63 0.564 0.6 OSSl
F,
24 ele. 1.472 1.461 -
1.488 1.446
Analysis of 2D fracture problems
LO
I
I
I
I
I
I
0.9.
0.0
I 0.8
0.4
929
I
I
0.0
0
Fig. 4. Normalized stress intensity factor F, for the square plate with a central crack.
close. The deformed profiles and maximum shear stress normalized by the average stress t3 for the cases of 8 = 45” and 8 = 90” are displayed from Figs 7-10. The average stress 6 is defined as
s W
d=
tsyudx 12W.
-W
The relative opening and sliding deformations of crack surfaces are easily seen in Figs 7 and 9. The stress concentration phenomena are also found in Figs 8 and 10.
1P
-1
0.0
0.a
Fig. 5. Normalized stress intensity factor F, for the square plate with an inclined crack (o/IV = 0.4).
WEN-HWA CHEN and CHE~Q-S~OUNG
930
-
CHANG
a elemants 24 ebment#
i
e Fig. 6. Normalized stress intensity factor Fir for the square plate with an inclined crack (a/W = 0.4).
Example 2: the rectangular plate with two inclined cracks
To examine the interaction effects among multiple cracks and displacement boundaries, the problem of a square plate with two inclined cracks subjected to uniform displacements at both ends as shown in Fig. 11 is solved. Figures 12 and 13 display the variations of normalized mode I and mode II stress intensity factors (F, , F,, ) with increasing values of inclined angles 8. Also shown are the results obtained in ref. [ 17] for an infinite plate with uniform tension Q at y = of:co, As predicted, the computed results are always larger than those of infinite plate. The deformed profiles and maximum shear stress contours for 8 = 90” and 180” are displayed clearly from Figs 14-l 7. Again, good convergence of normalized stress intensity factors and low order of ~lynomials N taken are also noted. CONCLUSIONS The finite element alternating method has been successfully extended to analyse mixed-mode fracture problems with single or multiple cracks under mixed boundary conditions. Accurate stress
s-w..-_
U-d
-lWOTEWd
Fig. 7. The deformed proties for the squam plate with an inclined crack (0 = 45”).
Fig. 8. Normalized maximum shear stress contours for the square plate with an inclined crack (0 = 45”).
Analysis of 2D fracture problems
-----_ -
931
undafonnad Daformad
Fig. 9. The deformed profiles for the square plate with an inclined crack (0 = 90”).
Fig. 10. Normalized maximum shear stress contours for the square plate with an inclined crack (8 = 90’).
J
---------
up0
1
uy=0oMt.
l.la
612
2a/W=O.6
--------m
t
ux-02~=-ooMt
_I
Fig. 11. The square plate with two inclined cracks under displacement boundary conditions.
932
WEN-HWA CHEN nrnd CHENQ-SHYOUNG
CHANG
t6
LO 0.6 a6 a4 0.0 a0
v
w
60’
so’
190’
l60’
l60”
e Fig. 12. Normalized stress intensity factor F, for the square plate with two inclined cracks.
all
Eh 0.6
0.4
a6
a0
v
60.
60’
80.
l60*
isO*
I#)’
e Fig. 13. Normalized mode II stress intensity factor F,, for the square plate with two inclined cracks.
Analysis of 2D fracture problems
--___-
undeformed
-
Ddomd
Fig. 14. The deformed profiles for the square plate with two inclined cracks (6 = 90’).
_____-
933
Fig. 15. Normalized maximum shear stress contours for the square plate with two inclined cracks (6 = 90”).
Undeformed Deformed
Fig. 16. The deformed profiles for the square plate with two inclined cracks (0 = 180”).
Fig. 17. Normalized maximum shear stress contours for the square plate with two inclined cracks (6 = 180”).
intensity factors can be evaluated using only a very limited number of regular elements and degrees of freedom. The deformed profiles and stress contours are also presented. Based on the simple finite element model adopted, the computer cost of the present computation is much less than other finite element models with quarter-point isoparametric elements or rigorous hybrid singular elements, especially for the cases with multiple cracks. REFERENCES 111R. S. Rarsoum, Application of quadratic isoparametric finite elements in linear fracture mechanics. Int. J. Fracture 10, 603-605 (1974). PI R. D. Henshell and K. G. Shaw, Crack tip finite elements are unnecessary. Inr. J. Numer. Merho& Engng 9,495-507 (1975).
934
WEN-HWA CHEN and CHENQ-SHYOUNG
CHANG
[3] S. N. Atluri, A. S. Kobayashi and M. Nakagaki, An assumed displacement hybrid finite element model for linear fracture mechanics. Int. J. Fracture 11, 257-271 (1975). [S] W. H. Chen and H. C. Lin, Flutter analysis of thin cracked panels using finite element method. AIAA J. 21.795-801 (1985). [5] W. H. Chen and K. Ting, Finite element analysis of mixed-mode thermoelastic fracture problems. Nucl. Engng and Design 90, 55-66 (1986). [6] W. H. Chen and C. B. Hwu, Hybrid singular element design for the bending analysis of bimaterial thin cracked plates. AIAA J. 24, 1399-1402 (1986). [7] K. Vijayakumar and S. N. Atluri, An embedded elliptical flaw, in an infinite solid, subject to arbitrary crack-face tractions. J. uppl. Mech. 48, 88-96 (1981). [8] T. Nishioka and S. N. Atluri, Analytical solution for embedded elliptical cracks, and finite element alternating method for elliptical surface cracks, subjected to arbitrary loadings. Engng Fracture Mech. 17, 247-268 (1983). [9] P. E. G’Donoghue, T. Nishioka and S. N. Atluri, Multiple surface cracks in pressure vessels. Engng Fracture Mech. 20, 545-560 (1984). [IO] I. N. Sneddon and H. A. Elliott, The opening of a Griffith crack under internal pressure. Q. appl. Math. 4, 262-267 (1946). [1 I] W. H. Chen and C. S. Chang, Analysis of two dimensional mixed-mode crack problems by finite element alternating method. To appear in Comput. Struct. [12] I. N. Sneddon, Fourier Transforms. McGraw-Hill, New York (1951). [13] F. Erdogan, Stress intensity factors. J. uppl. Mech. SO, Trans. ASME, 992-1002 (1983). [14] H. Tada, P. C. Paris and G. R. Irwin, The Stress Analysb ofcracks Handbook. Del Research Corporation, Hellertown, Pennsylvania (1973). [ 151 B. S. Annigeri and M. P. Cleary, Surface integral finite element hybrid (SIFEH) method for fracture mechanics. Int. J. Numer. Metho& Engng 20, 869-885 (1984). [16] M. Isida, Effect of width and length on stress intensity factors of internally cracked plates under various boundary conditions. Int. J. Fracture 7, 301-316 (1971). [17] G. C. Sih, Handbook of Stress Intensity Factors. Lehigh University, Bethlehem, Pennsylvania (1973). (Received 30 November 1988)
pp.921-934 1989
~1~7~/89
$3.00 -t 0.00
Pergamon Press pk.
ANALYSIS OF TWO DIMENSIONAL FRACTURE PROBLEMS WITH MULTIPLE CRACKS UNDER MIXED BOUNDARY CONDITIONS WEN-HWA
CHEN and CHENQ-SHYOUNG
CHANG
Department of Power Mechanical Engineering, National Tsing Hua University, Hsinchu, Taiwan 30043, Republic of China Abstract-This paper presents an efficient finite element alternating method for the analysis of two dimensional mixed-mode fracture problems with multiple cracks under mixed boundary conditions. Based on the analytical solution derived for an unbounded crack body with a central crack under arbitrary crack-face tractions, the resultant residual displacements/stresses on external boundaries are repeatedly evaluated till the variation of stress intensity factors and resultant residual stresses on each crack are negligible. Excellent agreement between the computed results and referenced solutions is observed. The deformed profiIes and m~mum shear stress contours are also presented. The economic computer time used and simple finite element idealization needed show the advantages of the method.
IN DEALING with engineering fracture mechanics problems, many numerical methods have heen developed and employed because of the complicated boundary conditions. Among those numerical methods, the finite element method has received much attention by many researchers. In this method, simple quarter-point isoparametric element@, 21or rigorous hybrid singular elements[3-6] are usually adopted to model the region near the crack-tip. Numerical experiments indicate that the accuracy of the stress intensity factors obtained by these techniques largely depends on the mesh size taken and/or the computer time consumed. To compensate those shortcomings, especially for solving realistic fracture problems with multiple cracks, a finite element alternating method which involves the iterative superposition of finite element solution of an uncracked structure and analytical solution of an abounded crack body subjected, to arbitrary normal and shear loadings on crack surfaces is then studied. Vijayakumar and Atluri[7] derived the complete analytical solution of a three dimensional body with an elliptical crack under arbitrary crack-face tractions. The solution later was successfully implemented in conjunction with the finite element method to compute mode I stress intensity factors of an elliptical crack[8] and multiple cracks[9]. Unfortunately, since the elliptical axes taken in Vijayakumar and Athni[7] are finite, the solutions obtained cannot be directly applied to solve two dimensional problems. Besides, those studies as mentioned above are only limited to the problems under force boundary conditions. Recently, based on the solution obtained by Sneddon and Elliott[lO] for the problem with crack surface under symmetric normal pressures, the two dimensional analytical solution of an infinite plate with a crack under arbitrary normal and shear loadings together with an efficient fimte element alternating t~hnique have been successfully developed by the authors[ 11] to deal with mixed-mode fracture problems with multiple cracks. Similarly, however, only force boundary conditions are considered. Hence, for more general engineering applications, the extension of the previous work[l I] to analyse fracture problems under displacement/force boundary conditions becomes the main objective of this work. To satisfy the prescribed mixed boundary conditions, based on the analytical solution derived, the resultant residual nodal displa~m~ts and resultant residual equivalent nodal forces are repeatedly computed during the iterative sohrtion process until the variation of stress intensity factors and resultant residual stresses on each crack converge to a small value. The deformed profiles and maximum shear stress contours are then presented. The high efficiency and accuracy of the technique developed can be noted. 921
WEN-HWA CHEN and CHENQ-SHYOUNG CHANG
922
ANALYTICAL SOLUTIONS
As shown in Fig. l(a), a two dimensional infinite plate with a crack at x = 0 and - 4 d y < a, where a is half of the crack length, is considered. The crack surfaces are subjected to arbitrary normal and shear loadings which are represented by a polynomial of any orders. The complete analytical solution can be superposed by two analytical solutions of normal and shear loadings, respectively. For different loading considerations on crack surfaces, the boundary conditions are to be:
(0 stresses ((i=, cryy,r xY) and displa~ments (u,, u,,)-*Oas x2 + y*-+co (ii) for normal loading case (see Fig. lb): atx-0
(iia) rxY= 0 (iib) a,, = - i @(y/a)”
atx=Oand[yl
n-0
and (iic) U*= 0
atx=Oandlyl>a
(iii) for shear loading case (see Fig. lc): (iiia) a, = 0
atx==O
(iiib) TV,, = -
atx=Oand/y(,
and (iiic) f+ = 0
atx==Oand/y(>a.
Fig. I. An infinite plate with a crack under arbitrary crack surface tractions. (a) Arbitrary crack surface tractions. (b) Arbitrary normal crack surface tractions, (c) Arbitrary shear crack surface tractions.
Analysis of 2D fracture problems
923
In the above conditions, N is an arbitrary positive integer and C$ and C!,?) are polynomial coefficients. From the symmetry and anti-symmetry of loading cases with respect to y-axis, only the semi-infinite elastic medium for x 2 0 needs to be solved. Based on the plane strain assumptions and Sneddon and Elliott’s solution[lO] for symmetric normal loading, neglecting the effect of body forces, the analytical solutions of stresses and displacements are derived as follows: (a) for normal loading case: 0 xx=--
3 j:
* 0W@W)cos(W Oyy= - ?r s zxy
=
-
24, =
7
2
&)eetXd<
(1)
- Sx)e-[“d5
(2)
@‘,‘)(<)cos(~y) + &“(<)sin(
{:
+ &‘Wsin(WNl
{31’YOsin(tv)
-
&‘)
>WO>>&
-{*
)cos(
om ($(‘)(t
dt
(3)
1 - v) + 5x]e + $5
(4)
- 2v - tx)e-‘x$
(5)
s and uy = - y
0m{$1’)(<)sin({y) - @)(<)cos([y)}(l s (b) for shear loading case: d XX=-- ‘,” 10a #PWsin(
- &2’(t>cos(~~)}5e-tx d5 - &2)Wcos(t~)l(*
* 0m@\2)(0sin(W Oyy= - IL s z XY=--
1
24,= 7
om
s
- We-”
{$\2’(Ocos(5~) + $i2)(Wn(5y))U
Joa
(6) dt
(7)
-
(8)
{$\2)(5)sin(
(9)
and uv = - y
1: ($\2)(C;)cos(~y) + $P)(<)sin(ry))[
- 2(1 - v)]e+
7
(10)
where E is Young’s modulus and v is Poisson’s ratio. Equations (l)-(lO) have satisfied the boundary conditions (i), (iia) and (iiia) for appropriate respective cases. Further, after satisfying other boundary conditions (iib), (iic), (iiib) and (iiic), the functions $\‘)(<) and $p)([) can be determined by setting the functions #j(t) and $i2)(<) as zero and the polynomial power n in the boundary conditions (iib) and (iiib) as even. Similarly, #j(e) and $J$~)(<)are determined by setting $\‘)(r) and $p)({) as zero and the polynomial power n as odd. Define the variables
g:)(q)
= a(z/2q)“2
i
C$%f”;
It-0 wPtl)
sin@rl)
=
(wtl
W2J-
=
(7~
/2)‘/2Jl,2
+o (P)
’
!z a
= p “2F,p’(p)
I,2 (PI);
6~);
(11)
924
WEN-HWA CHEN and CHENQ-SHYOUNG
CHANG
where p is equal to 1 and 2. for normal and shear loading cases, respectively. The value q = 1 denotes the loading case when the power n is even and q = 2-represents the loading case when the power n is odd. J_,,2(plf) and J,,*(ptf) are the Bessel functions of the first kind. Two sets of dual integrals for the cases of even n and odd n are then obtained as follows:
and
From the solutions of dual integrals as shown in [12], the solutions @‘)(p) and ~~(~) of the above dual integrals can be obtained. Then, from eq. (1 I), the functions #“<<) and &p)(t) are determined to be
and
where r(n) denotes the Gamma function. Substituting the ,expressions of &p)(t) and @‘)(t) in eqs (12) and (13) into eqs (l)-(lO), the closed-form solutions for both, of normal and shear .loading cases are thus found. Through integration and appropriate superposition, the complete analytical solution of stresses and displacements can be rearranged and summarized in a matrix .form as W = WI VI
(14)
and +> =
WIWI
WI
where [IV] and [U] are the matrix of relation function; The stress vector .{a], displacement vector (u} and the coefficient vector {C> of the polynomials are
(16)
Analysis of 2D fractureproblems
From the definitions[l3],
925
the stress intensity factor vector {K} are thus determined by WI =
PI w
(17)
where [S] represents the relation matrix. Consequently, once the coefficient vector {C} in eq. (16) is found, the stress distributions and deformatians at any locations in the infinite plate and the stress intensity factors at two crack tips can be computed from eqs (14), (15) and (17), respectively.
FINITE
ELEMENT
ALTERNATING
PROCEDURES
The finite element alternating procedures in solving the fracture problem with multiple cracks in a bounded plate subjected to mixed boundary conditions can be stated as follows: (1) Solve the bounded uncrack plate under given external loadings and displacements as the problem considered and evaluate the finite element solution of residual stresses a{ at the location of fictitious crack k. The subscript k is the number of the crack (k = 1,2, . . . , K). (2) Since the crack surfaces are stress free, reverse the sign of the residual stresses u/k and fit them with polynomials of appropriate orders. The fitted stresses ai can be expressed as ai = {L}‘{C}R where (C} E is the polynomial coefficient vector of crack k and {L} is the vector composed of @/Q)” (n = 0, 1,2, . . .). The value a, is half of the length of kth crack. For better accuracy, {C}: is determined by the least-square method through the minimization of the integral, 4
(o-l;-a:)2dS
(k=1,2
,...,
K).
s Thus,
where [Hlk =
s a’
-*k
and
{R>k=
I”-ak {L} af dS.
(3) From (C)i and eq. (14), the reversed residual interaction coefficient vectors of interaction effects induced by crack k at fictitious crack I (1 = 1,2, . . . , k - 1, k + 1, . . . , K){A},k, can be evaluated. If I> k, add {A},k to (C)y and compute the updated (C)y, Thus, the coefficient vector (C)A of the residual stresses for crack k need to be released for next computation is found to be
Repeat all steps as mentioned earlier until the cumulative residual stresses on each crack are negligible. As a result, the resultant coefficient vector {C}, = (C)f + (C)i + {C}: + - * * for the current iteration is obtained. The corresponding resultant residual stresses and displacements can be computed by eq. (14) and eq. (15), respectively. (4) Substituting the resultant coefficient vectors {C}, into eq. (17), the stress intensity factors of each crack for the current iteration can be calculated and updated. (5) The stresses/displacements of the cracked plate for each iteration can be computed by the sum of finite element solutions and those solutions calculated by each resultant coefficient vectors {C}, through eqs (14) and (15). The final stresses/displacements of the cracked plate are then readily determined once the iteration is completed.
926
WEN-~A
CHEN and C~EN~S~O~G
CHANG
(6) Once the variation of stress intensity factors of current iteration for all crack tips and the ratio of resultant residual stresses/permissible stress for each crack are smaller than a small value Y (Y = O.@W, say, 1 @wmt
iteration)
_
j.pwious
iteration)
/< y
1zi&@m-m itention) J
and ” ~W’~G>,~ dS s -=k cy 2&a
forallk=1,2,...,K
the iteration is finished. Here, i is equal to I and II for mode I and mode II, respectively. (7) To satisfy the force boundary conditions on external boundaries of the fracture problem (from eq. 14), the resultant residual equivalent nodal forces {Q}: on external force boundaries at mth element are computed by the sum of the contribution of each crack (Q), (k = 1,2, . . . , AC) as
iQ>: = j, fQ3nlk where
and
PI, =
WL’IPl~ dS s %I
[G]& denotes the relation matrix and S, is the portion of the boundary of mth element where force boundary conditions applied. [N] represents the shape function, ET]&is the transfo~ation matrix, and [wk is the matrix of relation function for kth crack as defined in eq. (14). To satisfy the displacement boundary conditions on external boundaries (from eq. 15), the resultant residual nodal displacements {~}1+on external displacement boundaries at node r are calculated by the sum of the contribution of crack k (k = 1,2,. . . , K), (ujlk, and expressed as
where
and
KJh is the
matrix of relation function for kth crack as defined in eq. (15) which is estimated at node r. [Mlk is the transformation matrix taken for locating the kth crack in global coordinate system. (8) Reverse the resultant residual covalent nodal forces (Q)z and resultant residual nodal displacements (a>: calculated in procedures (7) and consider those as applied loadings and displacements acting on the external boundaries of the bounded uncrack plate. Repeat all the iteration procedures (1x7) until the variation of stress intensity factors and resultant residual stresses on each crack are negligible. RIEWLTS AND DIS~SSIONS
To demonstrate the validity of the technique developed in this work, a bounded plate with inclined cracks subjected to different force/displacement boundary conditions at upper or bottom boundaries is considered. As shown in Figs 2 and 3, the finite element meshes with 9 eight-node
927
Analysisof 2D fracture problems
Fig. 2. The rectangular plate with an inclined crack under mixed boundary conditions.
isoparametric elements and 80 degrees of freedom are adopted. All numerical computations performed by MicroVAX II computer. Example
are
I: the rectangular plate with an inclined crack
(a) Subjected to force condition at upper end and displacement condition at bottom end. As shown in Fig. 2, a rectangular plate with an inclined crack under uniform tension at y = H and fixed displacement in y-direction at y = -H is first analysed.
A’
r ---
2w
EI
ly=o
_--_--
uy==conat.
,P Pa
d
__-------
F=‘wv-
Fig. 3. The square plate with an inclined crack under displacement boundary conditions.
928
WEN-HWA
CHEN and CHEN~SHYOUN~
CHANG
Table 1. Normalized stress intensity factors for the rectangular plate with an inclined crack
4 B 30” 45” 60” 90”
41
Tada 1141
Annigeri WI
Present 9 ele,
0.73 1.488
0.386 0.706 1.01 1.446
0.363 0.685 1.066 I.406
24 ele.
Tada 1141
Annigeri WI
9 ele.
24 ele.
0.357 0.706 1.083 1.461
0.6 0.000
0.515 0.551 0.458 0.000
0.506 0.564 0.492 0.000
0.516 OS87 0.492 0.000
Present
Table 1 shows the values of normalized mode I and mode II stress intensity factors (4, I;;,) for various inclined angles 8. (Fi, F,,) = (&, ~~*)/~~. Good agreement between the computed results and referenced data by Tada et aL[14]and Annigeri and Cleary[lS] is observed. To evaluate the high efficiency of the present technique, another run with 24 elements is also performed. It is obvious that the solutions obtained have only little difference between the use of 9 elements and 24 elements. Hence, the simpler mesh of 9 elements is su~ciently good to obtain accurate stress intensity factors for this case. The excellent monotone convergence of normalized stress intensity factors for the cases of 0 = 45” and 6 = 90” is indicated in Table 2. Table 3 displays the computed normalized stress intensity factors with different orders of polynomial (N) taken for fitting the normal and shear stresses on the fictitious crack. For 8 = 90” case, since the normal stress on fictitious crack is symmetric with respect to the center of the crack, the polynomial order can be chosen as even only. It is seen that N = 2 is big enough to reach good solutions. (b) Subjected to ~~~~Iace~nt condition at both upper and bottom en&. To analyse the fracture problem under prescribed displacement boundary conditions using finite element alternating method (as shown in Fig. 31, a square plate with an inclined crack subjected to uniform ~splacements at upper and bottom boundaries is also studied. For B = 90” case, Fig. 4 displays the variations of normalized mode I stress intensity factor F$ with increasing values a/W. The results obtained by Isida[l6] are also compared. Figures 5 and 6 show the computed normalized stress intensity factors F, and F,, for the case of a/W = 0.4 vs various angles of 8 using the mesh of 9 elements and 24 elements, respectively. Both solutions are also Table 2. Convergence of normalized stress intensity factors for the rectangular plate with an inclined crack e =900
B =4S” Cycle bf iterations 1
i 3 4 5 Tada[Y4]
&
Tabie 3,
4 9 ele.
24 ele,
o.soo
o.soo
0.559 0.564 0.564 0.564 0.564
0.576 0.586 0.587 0.587 OS87
1.000 1.286 1.371 1.3% 1.404 1.406
1.000 1.322 1.420 1.450 I .459 1.461
24 ete.
9 ele.
0.500 0.623 0.664 0.679 0.683 0.685 0.73
OS00 0.626 0.676 0.695 0.703 0.706
Anniaerill51
41
24 ele.
9 ele.
1.488 1.446
0.6 0.551
0.706
Comparison of different orders of ~l~orni~s plate with an inclined crack
e=9oo
0 =4S” Order N 2 3 4 5 Tada[ 141 Annigeri[ 151
9 ele.
F,
24 eie.
9 ele.
0.706 0.706 0.706 0.706
:z
0.683 0.684 0.683 0.685 0.73 0.706
fitted for the rectangular
F,I
24 ele.
9 ele.
0.585 0.585 0.586 0.587
1.399 1.406 -
O:S63 0.564 0.6 OSSl
F,
24 ele. 1.472 1.461 -
1.488 1.446
Analysis of 2D fracture problems
LO
I
I
I
I
I
I
0.9.
0.0
I 0.8
0.4
929
I
I
0.0
0
Fig. 4. Normalized stress intensity factor F, for the square plate with a central crack.
close. The deformed profiles and maximum shear stress normalized by the average stress t3 for the cases of 8 = 45” and 8 = 90” are displayed from Figs 7-10. The average stress 6 is defined as
s W
d=
tsyudx 12W.
-W
The relative opening and sliding deformations of crack surfaces are easily seen in Figs 7 and 9. The stress concentration phenomena are also found in Figs 8 and 10.
1P
-1
0.0
0.a
Fig. 5. Normalized stress intensity factor F, for the square plate with an inclined crack (o/IV = 0.4).
WEN-HWA CHEN and CHE~Q-S~OUNG
930
-
CHANG
a elemants 24 ebment#
i
e Fig. 6. Normalized stress intensity factor Fir for the square plate with an inclined crack (a/W = 0.4).
Example 2: the rectangular plate with two inclined cracks
To examine the interaction effects among multiple cracks and displacement boundaries, the problem of a square plate with two inclined cracks subjected to uniform displacements at both ends as shown in Fig. 11 is solved. Figures 12 and 13 display the variations of normalized mode I and mode II stress intensity factors (F, , F,, ) with increasing values of inclined angles 8. Also shown are the results obtained in ref. [ 17] for an infinite plate with uniform tension Q at y = of:co, As predicted, the computed results are always larger than those of infinite plate. The deformed profiles and maximum shear stress contours for 8 = 90” and 180” are displayed clearly from Figs 14-l 7. Again, good convergence of normalized stress intensity factors and low order of ~lynomials N taken are also noted. CONCLUSIONS The finite element alternating method has been successfully extended to analyse mixed-mode fracture problems with single or multiple cracks under mixed boundary conditions. Accurate stress
s-w..-_
U-d
-lWOTEWd
Fig. 7. The deformed proties for the squam plate with an inclined crack (0 = 45”).
Fig. 8. Normalized maximum shear stress contours for the square plate with an inclined crack (0 = 45”).
Analysis of 2D fracture problems
-----_ -
931
undafonnad Daformad
Fig. 9. The deformed profiles for the square plate with an inclined crack (0 = 90”).
Fig. 10. Normalized maximum shear stress contours for the square plate with an inclined crack (8 = 90’).
J
---------
up0
1
uy=0oMt.
l.la
612
2a/W=O.6
--------m
t
ux-02~=-ooMt
_I
Fig. 11. The square plate with two inclined cracks under displacement boundary conditions.
932
WEN-HWA CHEN nrnd CHENQ-SHYOUNG
CHANG
t6
LO 0.6 a6 a4 0.0 a0
v
w
60’
so’
190’
l60’
l60”
e Fig. 12. Normalized stress intensity factor F, for the square plate with two inclined cracks.
all
Eh 0.6
0.4
a6
a0
v
60.
60’
80.
l60*
isO*
I#)’
e Fig. 13. Normalized mode II stress intensity factor F,, for the square plate with two inclined cracks.
Analysis of 2D fracture problems
--___-
undeformed
-
Ddomd
Fig. 14. The deformed profiles for the square plate with two inclined cracks (6 = 90’).
_____-
933
Fig. 15. Normalized maximum shear stress contours for the square plate with two inclined cracks (6 = 90”).
Undeformed Deformed
Fig. 16. The deformed profiles for the square plate with two inclined cracks (0 = 180”).
Fig. 17. Normalized maximum shear stress contours for the square plate with two inclined cracks (6 = 180”).
intensity factors can be evaluated using only a very limited number of regular elements and degrees of freedom. The deformed profiles and stress contours are also presented. Based on the simple finite element model adopted, the computer cost of the present computation is much less than other finite element models with quarter-point isoparametric elements or rigorous hybrid singular elements, especially for the cases with multiple cracks. REFERENCES 111R. S. Rarsoum, Application of quadratic isoparametric finite elements in linear fracture mechanics. Int. J. Fracture 10, 603-605 (1974). PI R. D. Henshell and K. G. Shaw, Crack tip finite elements are unnecessary. Inr. J. Numer. Merho& Engng 9,495-507 (1975).
934
WEN-HWA CHEN and CHENQ-SHYOUNG
CHANG
[3] S. N. Atluri, A. S. Kobayashi and M. Nakagaki, An assumed displacement hybrid finite element model for linear fracture mechanics. Int. J. Fracture 11, 257-271 (1975). [S] W. H. Chen and H. C. Lin, Flutter analysis of thin cracked panels using finite element method. AIAA J. 21.795-801 (1985). [5] W. H. Chen and K. Ting, Finite element analysis of mixed-mode thermoelastic fracture problems. Nucl. Engng and Design 90, 55-66 (1986). [6] W. H. Chen and C. B. Hwu, Hybrid singular element design for the bending analysis of bimaterial thin cracked plates. AIAA J. 24, 1399-1402 (1986). [7] K. Vijayakumar and S. N. Atluri, An embedded elliptical flaw, in an infinite solid, subject to arbitrary crack-face tractions. J. uppl. Mech. 48, 88-96 (1981). [8] T. Nishioka and S. N. Atluri, Analytical solution for embedded elliptical cracks, and finite element alternating method for elliptical surface cracks, subjected to arbitrary loadings. Engng Fracture Mech. 17, 247-268 (1983). [9] P. E. G’Donoghue, T. Nishioka and S. N. Atluri, Multiple surface cracks in pressure vessels. Engng Fracture Mech. 20, 545-560 (1984). [IO] I. N. Sneddon and H. A. Elliott, The opening of a Griffith crack under internal pressure. Q. appl. Math. 4, 262-267 (1946). [1 I] W. H. Chen and C. S. Chang, Analysis of two dimensional mixed-mode crack problems by finite element alternating method. To appear in Comput. Struct. [12] I. N. Sneddon, Fourier Transforms. McGraw-Hill, New York (1951). [13] F. Erdogan, Stress intensity factors. J. uppl. Mech. SO, Trans. ASME, 992-1002 (1983). [14] H. Tada, P. C. Paris and G. R. Irwin, The Stress Analysb ofcracks Handbook. Del Research Corporation, Hellertown, Pennsylvania (1973). [ 151 B. S. Annigeri and M. P. Cleary, Surface integral finite element hybrid (SIFEH) method for fracture mechanics. Int. J. Numer. Metho& Engng 20, 869-885 (1984). [16] M. Isida, Effect of width and length on stress intensity factors of internally cracked plates under various boundary conditions. Int. J. Fracture 7, 301-316 (1971). [17] G. C. Sih, Handbook of Stress Intensity Factors. Lehigh University, Bethlehem, Pennsylvania (1973). (Received 30 November 1988)