- Email: [email protected]
Finite element based computation of strain energy release rate by modified crack closure integral
Enginrcrins
Fmchue Mechanics
Printed in Great
Britain.
Vol. 30,
No. 2, pp. 227-231,
1988
0013-7944/88 $3.00+ .OO 1988 Pertqmon Press plc.
@
FINITE ELEMENT BASED COMPUTATION OF STRAIN ENERGY RELEASE RATE BY MODIFIED CRACK CLOSURE INTEGRAL R. SETHURAMAN and S.K.MAITI Mechanical Engineering Department, Indian Institute of Technology, Powai, Bombay 400 076, India modified crack closure integral method with square-root stress singularity elements is given for calculation of strain energy release rate for an in-plane extension of a crack. Case studies are presented to illustrate the improvement in accuracy.
Abstract-A
INTRODUCTION THE STRAIN energy release rate for an in-plane extension of a crack can be computed
through J integral [ l] and stiffness derivative procedure [2]. It can also be evaluated using Irwin’s [3] crack closure integral, whose mathematical form is as follows:
(2) where Gr, Gn are the energy release rates for mode I and mode II deformations, c,.,,, ~7~~are the normal and shear stress distributions ahead of the crack-tip, U”, U, are the relative ‘opening’ displacements of the crack faces, Aa is the virtual extension of the crack, and the crack-tip is assumed to be located at x = 0. The virtual extension Aa is equal to the size of the crack-tip element. A significant advantage of the crack closure integral method is that it allows a separation of the energies associated with the different modes in a mixed mode problem. This method was first adopted by Rybicki and Kanninen[4] in the finite element computation with 4-noded elements around the crack-tip and they termed the integral as the modified crack closure integral. Krishnamurthy er a1.[5] and Ramamurthy et a1.[6] have extended the method to include S-noded quadrilaterals and singularity elements. We present here an analysis, which is very similar to that of refs[5] and [6] but differs more in respect of assumed variation of the crack closure forces. We have assumed that the variation of stresses, hence the closure forces, is one order less than that of the displacement. This has helped to attain simple expressions for energy calculations and more accurate results using the singularity elements around the crack-tip. The method is discussed below for cases where the quarter point elements, 8-noded quadrilaterals and 4-noded quadrilaterals are used to approximate the stress and displacement field around the crack-tip. MATHEMATICAL
FORMULATION
When the quarter point singularity elements are used, the displacement along the edge 145 (Fig. la) varies as the square-root of the distance 1x1from the crack-tip. The displacement variation can be written as
U&f) = w1+ 227
57,
(3)
228
R. SETHURAMAN
and S. K. MAITI
t-’
i
.!
5
FYI
= 1 2.
k-t-t 3
4 12
FY~
uJpT>=p--& uxs
ux4
(b)
Fig. 1. (a) Quarter point element, (b) 8-noded quadrilateral and (c) 4-noded quadrilateral.
where x is suitably transformed
to 5’. The stress along 123 varies as l/A.
q(5) = Bo+-
We can therefore write
BI
1+5’
where B. and B1 are the two arbitrary constants and 4 is an appropriate natural coordinate. B. and B1 are obtained by equating the work done by ~~(5) with that by the concentrated nodal forces. This finally gives Bo= &(F,‘-4F,I), Bl=-
6&l
Aa ’
where Fyi and Fy2 are the crack closure forces at nodes 1 and 2 respectively. Using (l), (3) and (4) the following expression for Gi is obtained 1 Gi = lim Aa2Aa
q(5) W5’) dx,
where U,,, is the opening displacement at node 4. Similarly the expression for Gi, can be obtained as
(6) where F,, and FX2are the crack closure forces for mode II displacements and U,, is the relative sliding displacement at node 4. The corresponding expressions obtained by Ramamurthy et a1.[6], who made different assumptions relating to the variations of displacement and stress, are
Computation
+
GI
=
(c2l&l+
y&
+
229
of strain energy release rate
c22Fy2
KGIF,,
+
(C21&+
+
c23Fy3)
G2E2
+
C22E2
+
G3E3)
c23F,3)
(7)
&I
U4
u,,]
(8)
where FY3and Fx3 are the closure forces at node 3, UxSand U,, are the opening displacements node 5, C1i = 33~12 - 52,
Cl2 = 17 - 217r/4,
C2, = 14 - 331rl8,
C22=21~/16-7/2
at
Cl3 = 211~12- 32, and C23=8-211~18.
The condition Au+0 can be met by taking the size of the crack-tip element as small as possible. In the case of 8-noded quadrilaterals the displacement varies quadratically along 145 (Fig. lb) and the stress has a linear variation along 123. That is,
(9)
The constants B. and B1 are obtained in the same manner as stated earlier. FY1= $
[B. - BJ, and
FY2= 5 AaBo. The strain energy release rates GI and GII are then obtained as GI = &
@I
u,,
+
Fy2
(11)
uy4),
These results are the same as those obtained by Krishnamurthy ef a1.[5]. In the case of 4-noded quadrilaterals displacement has a linear variation along 13 (Fig. lc) and stress is constant along 12. Proceeding as before we obtain,
GI =-!-Fyi 2Aa Gu=-
1 F,, 2Aa
Uy3,
(13)
v,,.
(14)
These relations are again the same as those obtained Krishnamurthy et al. [5].
by Rybicki and Kanninen[4]
and
230
R. SETHURAMAN
and S. K. MAITI
Fig. 2 (a) Centre crack, (b) edge crack and (c) kinked crack.
CASE STUDIES AND DISCUSSION Results of case studies on a centre crack (Fig. 2a), an edge crack (Fig. 2b) and a kinked crack (Fig. 2c) are now given to illustrate the effectiveness of the scheme. In all the cases the quarter point elements have been used around the crack-tip to approximate the singularity. The crack-tip element size is 1% of the crack length a in the first two cases and 0.2% of the crack length a in the last case. Material is isotropic. The results in respect of SIF correction factor are presented in Tables 1 and 2 for the first two cases. In the last example the kink length is 4% of the main crack length a and the strain energy release rate for an in-plane extension of the kink has been computed by the J integral, stiffness derivative procedure and the two schemes (relations 5 and 7) of the crack closure integral procedure. These results and those due to an analytical solution[8,9] are presented in Table 3 up to a kink angle of 75”. In the first example (Table l), the reference solution differs by about 1% and less over the considered range of a/w. In the edge crack problem (Table 2), the present method (relation 5) of calculation gives an error less than 0.6% in the studied range of a/w. In these two cases the scheme (relation 7) proposed by Ramamurthy et al. gives more error. In the last example up to 75” kink angle, the results based on the present scheme differs by less than 7% from the analytical solution[8,9]. In this case the accuracy is expected to be influenced by the modelling of the singularity at the knee, which depends on the kink angle[ lo]. Our results are based on a discretization where the conventional 6-noded triangles are used around the knee. For this case the strain energy release rate based on
Table 1. Comparison of SIF correction factors for the centre crack problem SIF correction factor Y = K,/oo&
a/w
Analytical solution[7]
Crack closure integral Eq. (7)
Eq. (5)
%Error in column 4 with respect to analytical solution
0.3 0.4 0.5 0.6 0.7 0.8
1.0594 1.1118 1.1891 1.3043 1.4842 1.7989
1.1107 1.1651 1.2463 1.3679 1.5614 1.9048
1.0626 1.1147 1.1924 1.3087 1.4939 1.8223
0.3021 0.2608 0.2775 0.3373 0.6536 1.3007
Table 2. Comparison of SIF correction factors for the edge crack problem SIF correction factor Y = Kr/o,,&
a/w
Analytical solution[7]
Eq. (7)
Eq. (5)
%Error in column 4 with respect to analytical solution
0.3 0.4 0.5 0.6 0.7 0.8
1.6629 2.1066 2.8297 4.0299 6.3610 11.9890
1.7414 2.2155 2.9630 4.2148 6.6481 12.4763
1.6658 2.1191 2.8339 4.0406 6.3553 11.9191
0.1744 0.5934 0.1484 0.2655 0.0897 0.5791
Crack closure integral
Computation
231
of strain energy release rate
Table 3. Comparison of strain energy release rates for the kinked crack problem Strain energy release rate
Inclination of the kink 30” 45” 60 75”
.Z integral 0.1845 0.1508 0.1124 0.7574
x x x x
10-s 1O-3 1O-3 1O-4
Stiffness derivative procedure 0.1799 0.1471 0.1095 0.7355
x x x x
Crack closure integral Analytical solution[8,9]
1O-3 1O-3 1O-3 lo-“
0.1789 0.1497 0.1156 0.8142
x x x x
1O-3 1O-3 1O-3 lo+
Eq. (7) 0.1921 0.1569 0.1172 0.7871
x x x x
1O-3 1O-3 1O-3 lo+
Eq. (5) 0.1844 0.1513 0.1129 0.7590
x x x x
lo-) 1O-3 1O-3 1O-4
% Error in column 6 with respect to analytical solution 3.0743 1.0688 2.3356 6.7797
the stiffness derivative procedure is always less than those due to the J integral. The values based on eq. (5) are very close to those of the J integral method. REFERENCES [l] [2] [3] [4] [5] [6] [7] [8] ]9J [lo]
J. R. Rice, Trans. ASME J. appl. Mech. 35, 379-386 (1968). D. M. Parks, In?. J. Fracture 10, 487-502 (1974). G. R. Jrwin, Fracture I, Hanbuch der physik VI, (Edited by S. Fhrgge), pp. 558-590. Springer, Berlin (1958). E. F. Rybicki and M. F. Kanninen, Engng Frucrure Mech. 9, 931-938 (1977). T. Krishnamurthy, T. S. Ramamurthy, K. Vijayakumar and B. Dattaguru, Pmt. Znr. Conf. Finite Elements in Compututional Mechanics, pp. 891-900 (1985). T. S. Ramamurthy, T. Krishnamurthy, K. Badri Narayana, K. Vijayakumar and B. Dattaguru, Me&. Res. Commun. 13, 179-186 (1986). D. P. Rooke and D. J. Cartwright, Compendium of S@essZntensiry Factors. HMSO, London (1976). S. K. Main, Znt. J. Fracture 32, R33-R36 (1986). B. Cotterel and J. R. Rice, Znr. J. Fracrure 16, 155-169 (1980). M. L. Williams, Trans. ASME 1. appl. Me&. 19, 526-528 (1952). (Receioed
24 June 1987)
Fmchue Mechanics
Printed in Great
Britain.
Vol. 30,
No. 2, pp. 227-231,
1988
0013-7944/88 $3.00+ .OO 1988 Pertqmon Press plc.
@
FINITE ELEMENT BASED COMPUTATION OF STRAIN ENERGY RELEASE RATE BY MODIFIED CRACK CLOSURE INTEGRAL R. SETHURAMAN and S.K.MAITI Mechanical Engineering Department, Indian Institute of Technology, Powai, Bombay 400 076, India modified crack closure integral method with square-root stress singularity elements is given for calculation of strain energy release rate for an in-plane extension of a crack. Case studies are presented to illustrate the improvement in accuracy.
Abstract-A
INTRODUCTION THE STRAIN energy release rate for an in-plane extension of a crack can be computed
through J integral [ l] and stiffness derivative procedure [2]. It can also be evaluated using Irwin’s [3] crack closure integral, whose mathematical form is as follows:
(2) where Gr, Gn are the energy release rates for mode I and mode II deformations, c,.,,, ~7~~are the normal and shear stress distributions ahead of the crack-tip, U”, U, are the relative ‘opening’ displacements of the crack faces, Aa is the virtual extension of the crack, and the crack-tip is assumed to be located at x = 0. The virtual extension Aa is equal to the size of the crack-tip element. A significant advantage of the crack closure integral method is that it allows a separation of the energies associated with the different modes in a mixed mode problem. This method was first adopted by Rybicki and Kanninen[4] in the finite element computation with 4-noded elements around the crack-tip and they termed the integral as the modified crack closure integral. Krishnamurthy er a1.[5] and Ramamurthy et a1.[6] have extended the method to include S-noded quadrilaterals and singularity elements. We present here an analysis, which is very similar to that of refs[5] and [6] but differs more in respect of assumed variation of the crack closure forces. We have assumed that the variation of stresses, hence the closure forces, is one order less than that of the displacement. This has helped to attain simple expressions for energy calculations and more accurate results using the singularity elements around the crack-tip. The method is discussed below for cases where the quarter point elements, 8-noded quadrilaterals and 4-noded quadrilaterals are used to approximate the stress and displacement field around the crack-tip. MATHEMATICAL
FORMULATION
When the quarter point singularity elements are used, the displacement along the edge 145 (Fig. la) varies as the square-root of the distance 1x1from the crack-tip. The displacement variation can be written as
U&f) = w1+ 227
57,
(3)
228
R. SETHURAMAN
and S. K. MAITI
t-’
i
.!
5
FYI
= 1 2.
k-t-t 3
4 12
FY~
uJpT>=p--& uxs
ux4
(b)
Fig. 1. (a) Quarter point element, (b) 8-noded quadrilateral and (c) 4-noded quadrilateral.
where x is suitably transformed
to 5’. The stress along 123 varies as l/A.
q(5) = Bo+-
We can therefore write
BI
1+5’
where B. and B1 are the two arbitrary constants and 4 is an appropriate natural coordinate. B. and B1 are obtained by equating the work done by ~~(5) with that by the concentrated nodal forces. This finally gives Bo= &(F,‘-4F,I), Bl=-
6&l
Aa ’
where Fyi and Fy2 are the crack closure forces at nodes 1 and 2 respectively. Using (l), (3) and (4) the following expression for Gi is obtained 1 Gi = lim Aa2Aa
q(5) W5’) dx,
where U,,, is the opening displacement at node 4. Similarly the expression for Gi, can be obtained as
(6) where F,, and FX2are the crack closure forces for mode II displacements and U,, is the relative sliding displacement at node 4. The corresponding expressions obtained by Ramamurthy et a1.[6], who made different assumptions relating to the variations of displacement and stress, are
Computation
+
GI
=
(c2l&l+
y&
+
229
of strain energy release rate
c22Fy2
KGIF,,
+
(C21&+
+
c23Fy3)
G2E2
+
C22E2
+
G3E3)
c23F,3)
(7)
&I
U4
u,,]
(8)
where FY3and Fx3 are the closure forces at node 3, UxSand U,, are the opening displacements node 5, C1i = 33~12 - 52,
Cl2 = 17 - 217r/4,
C2, = 14 - 331rl8,
C22=21~/16-7/2
at
Cl3 = 211~12- 32, and C23=8-211~18.
The condition Au+0 can be met by taking the size of the crack-tip element as small as possible. In the case of 8-noded quadrilaterals the displacement varies quadratically along 145 (Fig. lb) and the stress has a linear variation along 123. That is,
(9)
The constants B. and B1 are obtained in the same manner as stated earlier. FY1= $
[B. - BJ, and
FY2= 5 AaBo. The strain energy release rates GI and GII are then obtained as GI = &
@I
u,,
+
Fy2
(11)
uy4),
These results are the same as those obtained by Krishnamurthy ef a1.[5]. In the case of 4-noded quadrilaterals displacement has a linear variation along 13 (Fig. lc) and stress is constant along 12. Proceeding as before we obtain,
GI =-!-Fyi 2Aa Gu=-
1 F,, 2Aa
Uy3,
(13)
v,,.
(14)
These relations are again the same as those obtained Krishnamurthy et al. [5].
by Rybicki and Kanninen[4]
and
230
R. SETHURAMAN
and S. K. MAITI
Fig. 2 (a) Centre crack, (b) edge crack and (c) kinked crack.
CASE STUDIES AND DISCUSSION Results of case studies on a centre crack (Fig. 2a), an edge crack (Fig. 2b) and a kinked crack (Fig. 2c) are now given to illustrate the effectiveness of the scheme. In all the cases the quarter point elements have been used around the crack-tip to approximate the singularity. The crack-tip element size is 1% of the crack length a in the first two cases and 0.2% of the crack length a in the last case. Material is isotropic. The results in respect of SIF correction factor are presented in Tables 1 and 2 for the first two cases. In the last example the kink length is 4% of the main crack length a and the strain energy release rate for an in-plane extension of the kink has been computed by the J integral, stiffness derivative procedure and the two schemes (relations 5 and 7) of the crack closure integral procedure. These results and those due to an analytical solution[8,9] are presented in Table 3 up to a kink angle of 75”. In the first example (Table l), the reference solution differs by about 1% and less over the considered range of a/w. In the edge crack problem (Table 2), the present method (relation 5) of calculation gives an error less than 0.6% in the studied range of a/w. In these two cases the scheme (relation 7) proposed by Ramamurthy et al. gives more error. In the last example up to 75” kink angle, the results based on the present scheme differs by less than 7% from the analytical solution[8,9]. In this case the accuracy is expected to be influenced by the modelling of the singularity at the knee, which depends on the kink angle[ lo]. Our results are based on a discretization where the conventional 6-noded triangles are used around the knee. For this case the strain energy release rate based on
Table 1. Comparison of SIF correction factors for the centre crack problem SIF correction factor Y = K,/oo&
a/w
Analytical solution[7]
Crack closure integral Eq. (7)
Eq. (5)
%Error in column 4 with respect to analytical solution
0.3 0.4 0.5 0.6 0.7 0.8
1.0594 1.1118 1.1891 1.3043 1.4842 1.7989
1.1107 1.1651 1.2463 1.3679 1.5614 1.9048
1.0626 1.1147 1.1924 1.3087 1.4939 1.8223
0.3021 0.2608 0.2775 0.3373 0.6536 1.3007
Table 2. Comparison of SIF correction factors for the edge crack problem SIF correction factor Y = Kr/o,,&
a/w
Analytical solution[7]
Eq. (7)
Eq. (5)
%Error in column 4 with respect to analytical solution
0.3 0.4 0.5 0.6 0.7 0.8
1.6629 2.1066 2.8297 4.0299 6.3610 11.9890
1.7414 2.2155 2.9630 4.2148 6.6481 12.4763
1.6658 2.1191 2.8339 4.0406 6.3553 11.9191
0.1744 0.5934 0.1484 0.2655 0.0897 0.5791
Crack closure integral
Computation
231
of strain energy release rate
Table 3. Comparison of strain energy release rates for the kinked crack problem Strain energy release rate
Inclination of the kink 30” 45” 60 75”
.Z integral 0.1845 0.1508 0.1124 0.7574
x x x x
10-s 1O-3 1O-3 1O-4
Stiffness derivative procedure 0.1799 0.1471 0.1095 0.7355
x x x x
Crack closure integral Analytical solution[8,9]
1O-3 1O-3 1O-3 lo-“
0.1789 0.1497 0.1156 0.8142
x x x x
1O-3 1O-3 1O-3 lo+
Eq. (7) 0.1921 0.1569 0.1172 0.7871
x x x x
1O-3 1O-3 1O-3 lo+
Eq. (5) 0.1844 0.1513 0.1129 0.7590
x x x x
lo-) 1O-3 1O-3 1O-4
% Error in column 6 with respect to analytical solution 3.0743 1.0688 2.3356 6.7797
the stiffness derivative procedure is always less than those due to the J integral. The values based on eq. (5) are very close to those of the J integral method. REFERENCES [l] [2] [3] [4] [5] [6] [7] [8] ]9J [lo]
J. R. Rice, Trans. ASME J. appl. Mech. 35, 379-386 (1968). D. M. Parks, In?. J. Fracture 10, 487-502 (1974). G. R. Jrwin, Fracture I, Hanbuch der physik VI, (Edited by S. Fhrgge), pp. 558-590. Springer, Berlin (1958). E. F. Rybicki and M. F. Kanninen, Engng Frucrure Mech. 9, 931-938 (1977). T. Krishnamurthy, T. S. Ramamurthy, K. Vijayakumar and B. Dattaguru, Pmt. Znr. Conf. Finite Elements in Compututional Mechanics, pp. 891-900 (1985). T. S. Ramamurthy, T. Krishnamurthy, K. Badri Narayana, K. Vijayakumar and B. Dattaguru, Me&. Res. Commun. 13, 179-186 (1986). D. P. Rooke and D. J. Cartwright, Compendium of S@essZntensiry Factors. HMSO, London (1976). S. K. Main, Znt. J. Fracture 32, R33-R36 (1986). B. Cotterel and J. R. Rice, Znr. J. Fracrure 16, 155-169 (1980). M. L. Williams, Trans. ASME 1. appl. Me&. 19, 526-528 (1952). (Receioed
24 June 1987)