- Email: [email protected]
Calculation of stress intensity factors of cracked reissner plates by the sub-region mixed finite element method
Compufws & Srrucrurrs Vol. 30. No. 4. pp. 837-840. Printed in Great Britain.
~5-7~91~8 $3.00 + 0.00 Pergamon Press ptc
I988
CALCULATION OF STRESS INTENSITY FACTORS OF CRACKED REXSSNER PLATES BY THE SUB-REGION MIXED FINITE ELEMENT METHOD MINFENG HUANG and YUQIW LONG
Tsinghua University, Beijing, The People’s Republic of China Abstract-The stress intensity factors of Reissner plates are calculated by the sub-region mixed finite element method. Stress intensity factors of modes I, II and III for some finite and infinite plates are calculated. The relative errors are less than 1%.
I. INTRODU~ION
2. STRESS FUNCTIONS AROUND THE CRACK-TiP
The early analysis of plate fracture problems was based on the Kirchhoff plate theory. References [l-3], although based on the Reissner plate theory, have only solved some plate fracture problems of mode I. In this paper, the stress intensity factors of modes I, II and III for the Reissner plate are calculated by the sub-region mixed finite element method [4$ The vicinity of the crack tip is treated as the complementary energy region (C-region), in which stress parameters (fif are unknowns, and the other part of the plate as the potential energy region (P-region), in which unknowns are nodal displacements (sj. The present method gives results with high accuracy. For a Reissner plate, the strain~isplacement and stress-strain relations are
In the polar coordinate system shown in Fig. 1, the stress field for symmetric bending can be expressed as follows [3]: --MrjD = OSA,,r -i/T(cos 3%/2 - 5 cos 0/2) x (1 - v) i- A,,(1 - v)
x (1 +cosz~)+(c~2~) x D,,r(cos 38 - cos %) + l.SA,,r’iZ(l - v)(cos58/2 -I- 3cos%/2) f 2A,,r(i - v) x (COS3% f kos e 1-b opj 3
Pa)
f&,(1
GM
-(C/2D~~,2r(cos3~
with
-v)(l
-cos2%) -cos8)
+ l.SA,,r’~“(v- 1) Iv
IBd=D v
I 1 0 [
0 0 (1 -v)/2
W
-2M,,/D(i
where w is the deflection, (+ 1 = [$,, Il/JT are the cross sectional rotations, fx > = [~X,xY,K+,.]~are the curvatures, {r f = [yr, rJf are the shearing strains, {M} = [M,, M,., M.J’ are the moments, (Q > = [QX,QJ are the transverse shear forces, D = Eh3/12(1 - G), C = Ch/1.2,h is the thickness, E is the modulus of elasticity, G is the shear modulus and v is the Poisson’s ratio.
-v)
= -A,,r-“2(sin
%/2 + sin 3%/2)
- 2A,,sin 2% + 2D12rk2 x (sin %- sm 3%) + 3&r
112
x (sin S/2 - sin Se/Z> + 4A,s x (-sin 837
3%+ i sin %)+ O(r3’*)
838
MINFENGHUANG
--MO/D = OSA,,r -‘:2(v - 1) (sin 3ej2
Y
e3 0
L
and YUQIULQNG
+ sin e/2) - (C/4D)D,,r”*
x
x (sin Se/2 - sin e/2)
+ 1.5A23r’S2(v- l)(sin58/2
Fig. 1.
- sin e/2) + 2A,,r(v - 1)
QrK = P/(1 + VIAIF I;2
x (sin 38 - 3 sine) + 0(r3i2)
x (cos 0/2+3cos36/2)
-2M,/D(l
- v) = A,,r -"2(COS3e/2 +fCOSo/z)
+ D,,[cosB+ o.wr*
+ D2,r’,2k2(cos 5012 - cos e/2)
x (COSB -cos3e)]+[2/(1 + v)]
+ 3A2,r “‘(~0s 50/2 - 0.2 cos e/2)
x A,,r312(cos8/2 - 5cos5e/2)
+ 4A,,r(cos 30 - cos 0) + 0(r3’2) Q,/C = [2/(3 + 3v)]A2,r”2
+ [2/(3 + 3v)]A,g 2
x (sin e/2 + sin 30/Z)
x (cos e - cos 38) + 0(r512)
+ D2,[(0.5r -“2 sin e/2 + (k2r312/3)
Qo/C = -[6/(1+ v)]A,,r”*
x (sin e/2 - sin 56/z)]
x (sin e/2 + sin 3e/2)
+ [2/(5 + Sv)]A,,r””
+&[-sine
- 1.5k2r2 x (sin e/2 - sin 5e/2)
x(sine
-jsin3e)]
+[lO/(l
+v)] + D,,[1.5 sin 30/2 - (3k2rs12/5)
x A,,r3’2(sin 5e/2 - sin e/2)
x (sin 7e/2 - sin 3e/2)]
+ [2/( 1 + v)]A,g2( - sin e
+ [2/( 1 + v)]A2,r2(sin 8 - 3 sin 38)
+ f sin 38) + O(r5,‘2). For skew symmetric pure resultants may be written as - M,/D = 0.5Ag
twisting,
(4) the
stress
-“*( 1 - v)
+ D,g3j2 sin Se/2 + 0(r5’2) Qo/C = [2/(1 + v)]A2,r”2 x (cos0/2+fcos30/2) + D2,[1.5r m”2cos0 + (5k2r3’2/3)
x (sin 38/2 - $sin e/2) x (cOse -
+ (C/4D)D,,r”’ x (sin se/2 - sin e/2)
0.2 cos 5e/2)1
+ [2/( 1 + v)]A23r3~2 x (C~S e/2 -
0.2 cos 5e/2)
+ 1.5A2,r”*(1 -v)
+ D23[1.5r’~2~~~30/2 + (7k2r512/5)
x (sin 5e/2 + 0.6 sin e/2)
x (cos 30/2 - $0~ 76/2)]
+ 2A,,r (1 - v) (sin 38 + sin e)
+ [6/( I + v)]A2,r2(cos 0 - cos 38)
+ 0(r”2)
+ 2.5D2,r3’*cos 50/2 + O(r’!*)
(5)
839
Calculation of stress intensity factors of cracked Reissner plates in which k2 = C/2D(l - v) = 5/2h2. Denote {P} = [A,,, A,,,
42,
Al39
AMIT for the symmetric
ifi} = [LA,>A,,,
A,,,
D23r
case
A243 D2,1T
b
2w=200
-I
Fig. 2.
for the skew symmetric case and the above stress resultants can be written as IN
=
From the stationary conditions of the energy functional K, (6) and {b} can be obtained as
PEtI(8) (6)
IQ) = F’sl~B~~
3. ENERGY FUNCTIONAL OF THE SUEREGION MIXED FINITE ELEMENT METHOD
The energy functional of the sub-region variational principle [4] is n =n,-n,+H
PC
(7)
with
‘Psl)dV. (10) energy on the interline r,,‘
(M:$:: + M:,ll/r: + Q$@)ds TP(
(15)
K, = 2@(1
(9)
s
IPI = w’m{~).
mixed
where [K] is the global stiffness matrix for the P-region and {p} is the equivalent nodal load vector. rr,. is the complementary energy in the C-region
Hp,=
(14)
Finally the stress intensity factors can be determined from the following relations
in which n, is the potential energy in the P-region
H,,, is the additional
is>=(I~1+ lw~cl-‘vl-‘b~
(11)
- v)A,,
K,, = -2@(1
- v)A,,/3
K,,, = -fiCD2,/2.
4. NUMERICAL
EXAMPLES
Example 1. Stress intensity factor K,” of the infinite cracked plate subjected to uniform moment A4 (Fig. 2)
The length of crack is 2a. In order to approximate the infinite plate, 2L = 2W = 20a is taken. The mesh for l/4 plate is shown in Fig. 3, in which the C-region is a rectangular region with the crack tip at the center. In the P-region the eight-node isoparametric elements are used. The computed results of K,” for this infinite plate are shown in Table 1, in which a comparison is made with the analytic solution[5] and other finite element results [3]. The relative errors of the present results are all below 1%.
where M;, M;, and Q:, are the boundary forces of the C-region on T,,‘ and I/J:, $I! and wr are the boundary displacements of the P-region on rpr. Equation (11) can also be written as f-f,,.= Substituting
{P)‘[Fl{S).
(16)
(12)
eqns (8) (9) and (12) into (7) we have b_;’
Fig. 3.
840
MINFENGHUANGand YUQIULONG Table 1. K; and K, for crack plate (/M$) K;
for infinite plate
hla
Anal. [5]
Present
0.2 0.5 1.0 1.5 2.0
0.647 0.693 0.741 0.781 0.816
0.6510 0.6970 0.7450 0.7848 0.8186
K, for finite plate
Ref. [3] 0.6726 0.7352 0.8144
Present
Ref. [3]
K;OIK,
0.7376 0.8108 0.8985 1.0059 1.0830
0.8036 0.8985 -
1.141 1.170 1.208 1.288 1.327
Table 2. K,, and K,,, for infinite plate K,,,l[fiHlh(l
K,,IHJ
+
hla
Present
Ref. [6]
Error
Present
Ref. [6]
0.2 0.5 1.0 1.5 2.0
0.1942 0.3974 0.5135 0.6467 0.7010
0.193 0.395 0.510 0.643 0.697
0.63% 0.62% 0.59% 0.58% 0.59%
0.1386 0.1437 0.1219 0.0818 0.0605
0.140 0.145 0.123 0.0825 0.061
Example 2. Stress intensity factor KI ofjinite cracked plate with 2L = 2 W = 4a under uniform moment M The results of K, are also shown
~11 Error - 1.OO% -0.92% -0.88% -0.89% -0.90%
O.OSa to 0.2a, the computation significant difference.
results
make
no
in Table 1.
Example 3. Stress intensity factors K,[ and KI,, of infinite cracked plate subjected to pure twisting moments H
The results of K,, and K,,, are shown in Table 2. 5. CONCLUSIONS It can be concluded that the sub-region mixed finite element method is an efficient method for the determination of stress intensity factors in cracked Reissner plates. The numerical results show that the present method can be used to calculate K,, K,, and K,,, of finite and infinite plates with various thickness and gives higher accuracy. It is suggested that the best size of the C-region is r,, = 0. la. In practice, however, when r0 varies from
REFERENCES 1.
R. S. Barsoum, A degenerate solid element for linear fracture analysis of plate bending and general shells. Inr.
J. Numer. Meth. Engng 10, 551-564 (1976). 2. G. Yagawa, Finite element analysis of stress intensity factors for plate extension and bending problems. ht. J. Numer. Meth. Engng 14, 727-740 (1979). 3. Li Yingzhi and Chuntu Liu, Analysis of Reissner’s plate bending fracture problem. Acfa Mech. Sin. 4, 366-375 (1983). 4. Yuqiu Long, Bingchen Zhi, Wenqi Kuang and Jian
Shari,, Subregion mixed finite element method for the calculation of stress intensity factor. Proc. Inf. Coflf. on FEM, Shanghai, pp. 738-740 (1982). 5. F. J. Hartanft and G. C. Sih, Effect of plate thickness on the bending stress distribution around through cracks. J. Math. Whys. 47, 276291 (1968). 6. N. M. Wang, Twisting of an elastic plate containing a crack. Int. J. Fracture Mech. 6, 367-378 (1970).
~5-7~91~8 $3.00 + 0.00 Pergamon Press ptc
I988
CALCULATION OF STRESS INTENSITY FACTORS OF CRACKED REXSSNER PLATES BY THE SUB-REGION MIXED FINITE ELEMENT METHOD MINFENG HUANG and YUQIW LONG
Tsinghua University, Beijing, The People’s Republic of China Abstract-The stress intensity factors of Reissner plates are calculated by the sub-region mixed finite element method. Stress intensity factors of modes I, II and III for some finite and infinite plates are calculated. The relative errors are less than 1%.
I. INTRODU~ION
2. STRESS FUNCTIONS AROUND THE CRACK-TiP
The early analysis of plate fracture problems was based on the Kirchhoff plate theory. References [l-3], although based on the Reissner plate theory, have only solved some plate fracture problems of mode I. In this paper, the stress intensity factors of modes I, II and III for the Reissner plate are calculated by the sub-region mixed finite element method [4$ The vicinity of the crack tip is treated as the complementary energy region (C-region), in which stress parameters (fif are unknowns, and the other part of the plate as the potential energy region (P-region), in which unknowns are nodal displacements (sj. The present method gives results with high accuracy. For a Reissner plate, the strain~isplacement and stress-strain relations are
In the polar coordinate system shown in Fig. 1, the stress field for symmetric bending can be expressed as follows [3]: --MrjD = OSA,,r -i/T(cos 3%/2 - 5 cos 0/2) x (1 - v) i- A,,(1 - v)
x (1 +cosz~)+(c~2~) x D,,r(cos 38 - cos %) + l.SA,,r’iZ(l - v)(cos58/2 -I- 3cos%/2) f 2A,,r(i - v) x (COS3% f kos e 1-b opj 3
Pa)
f&,(1
GM
-(C/2D~~,2r(cos3~
with
-v)(l
-cos2%) -cos8)
+ l.SA,,r’~“(v- 1) Iv
IBd=D v
I 1 0 [
0 0 (1 -v)/2
W
-2M,,/D(i
where w is the deflection, (+ 1 = [$,, Il/JT are the cross sectional rotations, fx > = [~X,xY,K+,.]~are the curvatures, {r f = [yr, rJf are the shearing strains, {M} = [M,, M,., M.J’ are the moments, (Q > = [QX,QJ are the transverse shear forces, D = Eh3/12(1 - G), C = Ch/1.2,h is the thickness, E is the modulus of elasticity, G is the shear modulus and v is the Poisson’s ratio.
-v)
= -A,,r-“2(sin
%/2 + sin 3%/2)
- 2A,,sin 2% + 2D12rk2 x (sin %- sm 3%) + 3&r
112
x (sin S/2 - sin Se/Z> + 4A,s x (-sin 837
3%+ i sin %)+ O(r3’*)
838
MINFENGHUANG
--MO/D = OSA,,r -‘:2(v - 1) (sin 3ej2
Y
e3 0
L
and YUQIULQNG
+ sin e/2) - (C/4D)D,,r”*
x
x (sin Se/2 - sin e/2)
+ 1.5A23r’S2(v- l)(sin58/2
Fig. 1.
- sin e/2) + 2A,,r(v - 1)
QrK = P/(1 + VIAIF I;2
x (sin 38 - 3 sine) + 0(r3i2)
x (cos 0/2+3cos36/2)
-2M,/D(l
- v) = A,,r -"2(COS3e/2 +fCOSo/z)
+ D,,[cosB+ o.wr*
+ D2,r’,2k2(cos 5012 - cos e/2)
x (COSB -cos3e)]+[2/(1 + v)]
+ 3A2,r “‘(~0s 50/2 - 0.2 cos e/2)
x A,,r312(cos8/2 - 5cos5e/2)
+ 4A,,r(cos 30 - cos 0) + 0(r3’2) Q,/C = [2/(3 + 3v)]A2,r”2
+ [2/(3 + 3v)]A,g 2
x (sin e/2 + sin 30/Z)
x (cos e - cos 38) + 0(r512)
+ D2,[(0.5r -“2 sin e/2 + (k2r312/3)
Qo/C = -[6/(1+ v)]A,,r”*
x (sin e/2 - sin 56/z)]
x (sin e/2 + sin 3e/2)
+ [2/(5 + Sv)]A,,r””
+&[-sine
- 1.5k2r2 x (sin e/2 - sin 5e/2)
x(sine
-jsin3e)]
+[lO/(l
+v)] + D,,[1.5 sin 30/2 - (3k2rs12/5)
x A,,r3’2(sin 5e/2 - sin e/2)
x (sin 7e/2 - sin 3e/2)]
+ [2/( 1 + v)]A,g2( - sin e
+ [2/( 1 + v)]A2,r2(sin 8 - 3 sin 38)
+ f sin 38) + O(r5,‘2). For skew symmetric pure resultants may be written as - M,/D = 0.5Ag
twisting,
(4) the
stress
-“*( 1 - v)
+ D,g3j2 sin Se/2 + 0(r5’2) Qo/C = [2/(1 + v)]A2,r”2 x (cos0/2+fcos30/2) + D2,[1.5r m”2cos0 + (5k2r3’2/3)
x (sin 38/2 - $sin e/2) x (cOse -
+ (C/4D)D,,r”’ x (sin se/2 - sin e/2)
0.2 cos 5e/2)1
+ [2/( 1 + v)]A23r3~2 x (C~S e/2 -
0.2 cos 5e/2)
+ 1.5A2,r”*(1 -v)
+ D23[1.5r’~2~~~30/2 + (7k2r512/5)
x (sin 5e/2 + 0.6 sin e/2)
x (cos 30/2 - $0~ 76/2)]
+ 2A,,r (1 - v) (sin 38 + sin e)
+ [6/( I + v)]A2,r2(cos 0 - cos 38)
+ 0(r”2)
+ 2.5D2,r3’*cos 50/2 + O(r’!*)
(5)
839
Calculation of stress intensity factors of cracked Reissner plates in which k2 = C/2D(l - v) = 5/2h2. Denote {P} = [A,,, A,,,
42,
Al39
AMIT for the symmetric
ifi} = [LA,>A,,,
A,,,
D23r
case
A243 D2,1T
b
2w=200
-I
Fig. 2.
for the skew symmetric case and the above stress resultants can be written as IN
=
From the stationary conditions of the energy functional K, (6) and {b} can be obtained as
PEtI(8) (6)
IQ) = F’sl~B~~
3. ENERGY FUNCTIONAL OF THE SUEREGION MIXED FINITE ELEMENT METHOD
The energy functional of the sub-region variational principle [4] is n =n,-n,+H
PC
(7)
with
‘Psl)dV. (10) energy on the interline r,,‘
(M:$:: + M:,ll/r: + Q$@)ds TP(
(15)
K, = 2@(1
(9)
s
IPI = w’m{~).
mixed
where [K] is the global stiffness matrix for the P-region and {p} is the equivalent nodal load vector. rr,. is the complementary energy in the C-region
Hp,=
(14)
Finally the stress intensity factors can be determined from the following relations
in which n, is the potential energy in the P-region
H,,, is the additional
is>=(I~1+ lw~cl-‘vl-‘b~
(11)
- v)A,,
K,, = -2@(1
- v)A,,/3
K,,, = -fiCD2,/2.
4. NUMERICAL
EXAMPLES
Example 1. Stress intensity factor K,” of the infinite cracked plate subjected to uniform moment A4 (Fig. 2)
The length of crack is 2a. In order to approximate the infinite plate, 2L = 2W = 20a is taken. The mesh for l/4 plate is shown in Fig. 3, in which the C-region is a rectangular region with the crack tip at the center. In the P-region the eight-node isoparametric elements are used. The computed results of K,” for this infinite plate are shown in Table 1, in which a comparison is made with the analytic solution[5] and other finite element results [3]. The relative errors of the present results are all below 1%.
where M;, M;, and Q:, are the boundary forces of the C-region on T,,‘ and I/J:, $I! and wr are the boundary displacements of the P-region on rpr. Equation (11) can also be written as f-f,,.= Substituting
{P)‘[Fl{S).
(16)
(12)
eqns (8) (9) and (12) into (7) we have b_;’
Fig. 3.
840
MINFENGHUANGand YUQIULONG Table 1. K; and K, for crack plate (/M$) K;
for infinite plate
hla
Anal. [5]
Present
0.2 0.5 1.0 1.5 2.0
0.647 0.693 0.741 0.781 0.816
0.6510 0.6970 0.7450 0.7848 0.8186
K, for finite plate
Ref. [3] 0.6726 0.7352 0.8144
Present
Ref. [3]
K;OIK,
0.7376 0.8108 0.8985 1.0059 1.0830
0.8036 0.8985 -
1.141 1.170 1.208 1.288 1.327
Table 2. K,, and K,,, for infinite plate K,,,l[fiHlh(l
K,,IHJ
+
hla
Present
Ref. [6]
Error
Present
Ref. [6]
0.2 0.5 1.0 1.5 2.0
0.1942 0.3974 0.5135 0.6467 0.7010
0.193 0.395 0.510 0.643 0.697
0.63% 0.62% 0.59% 0.58% 0.59%
0.1386 0.1437 0.1219 0.0818 0.0605
0.140 0.145 0.123 0.0825 0.061
Example 2. Stress intensity factor KI ofjinite cracked plate with 2L = 2 W = 4a under uniform moment M The results of K, are also shown
~11 Error - 1.OO% -0.92% -0.88% -0.89% -0.90%
O.OSa to 0.2a, the computation significant difference.
results
make
no
in Table 1.
Example 3. Stress intensity factors K,[ and KI,, of infinite cracked plate subjected to pure twisting moments H
The results of K,, and K,,, are shown in Table 2. 5. CONCLUSIONS It can be concluded that the sub-region mixed finite element method is an efficient method for the determination of stress intensity factors in cracked Reissner plates. The numerical results show that the present method can be used to calculate K,, K,, and K,,, of finite and infinite plates with various thickness and gives higher accuracy. It is suggested that the best size of the C-region is r,, = 0. la. In practice, however, when r0 varies from
REFERENCES 1.
R. S. Barsoum, A degenerate solid element for linear fracture analysis of plate bending and general shells. Inr.
J. Numer. Meth. Engng 10, 551-564 (1976). 2. G. Yagawa, Finite element analysis of stress intensity factors for plate extension and bending problems. ht. J. Numer. Meth. Engng 14, 727-740 (1979). 3. Li Yingzhi and Chuntu Liu, Analysis of Reissner’s plate bending fracture problem. Acfa Mech. Sin. 4, 366-375 (1983). 4. Yuqiu Long, Bingchen Zhi, Wenqi Kuang and Jian
Shari,, Subregion mixed finite element method for the calculation of stress intensity factor. Proc. Inf. Coflf. on FEM, Shanghai, pp. 738-740 (1982). 5. F. J. Hartanft and G. C. Sih, Effect of plate thickness on the bending stress distribution around through cracks. J. Math. Whys. 47, 276291 (1968). 6. N. M. Wang, Twisting of an elastic plate containing a crack. Int. J. Fracture Mech. 6, 367-378 (1970).