- Email: [email protected]
p-Version finite element approximations of stress intensity factors for cracked plates including shear deformation
Pe,gamon
p-VERSION
Engineering Fracture Mechanics Vol. 52, No. 3, pp. 493-502, 1995 Copyright ~? 1995 ElsevierScience Ltd 0013-7944(95)00019--4 Printed in Great Britain. All rights reserved 0013-7944/95 $9.50+ 0.00
FINITE ELEMENT
APPROXIMATIONS
OF
STRESS INTENSITY FACTORS FOR CRACKED PLATES INCLUDING SHEAR DEFORMATION KWANG SUNG WOO Department of Civil Engineering, Yeungnam University, Gyongsan 712-749, Korea and CHAE GUE LEE Department of Civil Engineering, Chosun University, Kwangju 501-759, Korea Abstract--The new p-version crack model is proposed to estimate the bending stress intensity factors of the thick cracked plate under flexure. The proposed model is based on high order theory and hierarchical C°-plate element including shear deformation. The displacement fields are defined by integrals of Legendre polynomials which can be classified into three groups, e.g. basic, side and internal mode. The computer implementation allows arbitrary variations of p-level up to a maximum value of 10. The bending stress intensity factors are computed by the virtual crack extension approach. The effects of ratios of thickness to crack length (h/a), crack length to width (a/W) and boundary conditions are investigated. Very good agreement with the existing solution in the literature is shown for the uncracked plate as well as the cracked plate.
1. INTRODUCTION ONE OFTHE BASICrequirements for the strength analysis of plates containing flaws or cracks is the knowledge of the singular character of the stress field in the neighborhood of the point at the crack tip. The SIF (stress intensity factor), which is in fact the strength of the stress singularities at crack tips, is applied to predict the static strength of cracked bodies by Irwin [1] for plane extension, symmetric with respect to the crack. The stress intensity factors have been shown to control the rate of crack propagation under cyclic loading in such situations. Several investigators have discussed the nature of the local stresses around a sharp crack in a thin plate subjected to out-of-plane bending loads. Based on the Poisson-Kirchhoff theory of thin plate and the technique of Fadle eigen function expansion, Williams [2] found that the elastic bending stresses near the tip of a semi-infinite crack vary as the inverse square root of the radial distance from the crack front. His results were not complete in that the strength or magnitude of the local stresses was left undetermined. Sih et al. [3] cleared the way for finding the coefficients in the eigen function expansions by application of the theory of complex functions. The results in refs [2, 3] were obtained from the classical fourth-order theory of thin plates, however, the edge conditions at the crack surfaces are satisfied only in an approximate manner in that the three physically natural boundary conditions of prescribing bending moment, twisting moment and transverse shear stress are replaced by two conditions. Owing to such a replacement, the stress distribution in the immediate neighborhood of the crack edges will naturally be affected and will not be accurate. To overcome the aforementioned shortcoming, Knowles and Wang [4] used Reissner's sixth-order plate theory to satisfy all the physically natural boundary conditions along the crack edges with the help of singular integral equations. However, the extension of the method to cracked plates with finite plate thickness does not appear to be tractable. For the examination of the effect of plate thickness, Hartranft and Sih [5] put forward higher order theories that included the effect of transverse shear deformations and essentially solved the problem of infinite plates subjected to uniform moment. Murthy et al. [6] used Reissner's plate theory and applied collocation, and successive integration methods to satisfy the boundary conditions for a finite plate with a through crack. The associated numerical difficulties led them to conclude the finite element method is better suited for such problems. Wilson and Thompson [7] used the h-version of the finite element method 493
494
KWANG SUNG WOO and CHAE GUE LEE
to compute the SIF for the through-thickness crack in a plate under flexure. They used a large number of three-node triangular C~-elements (actually 257 and 305 node models) and computed the value of SIF using the direct method. Hilton [8] used 48 20-noded three-dimensional isoparametric elements including the special singularity element (279 node model) to compute a SIF. Yagawa and Nishioka[9] used a superposition method of analytical and finite element solutions on the basis of 36 eight-noded isoparametric elements (133 node model) together with the reduced integration technique. Alwar and Nambissan [10] used 110 three-dimensional quadratic isoparametric elements, including degenerate triangular quarter-point prism elements at the crack tip region, to model one-eighth of the plate and computed the SIF by COD method. The results presented by earlier workers [8, 9] for the finite plate under uniform moments differ significantly from the results by Alwar and Nambissan. They indicated that the difference may be due to the coarseness of the mesh used by those. Chen and Chen [11] presented the hybrid~lisplacement singular element model based on the Kirchhoff plate theory for the bending analysis of thin cracked infinite plate and transverse pressure analysis of thin finite cracked plate with various boundary conditions. Kang and De Saxce [12] used 35 hybrid mongrel elements, which is the special crack tip element, for infinite cracked plate under bending and finite thin plate under transverse pressure. Since the singularity at the through crack corner of finite plate subjected to various loads and boundary conditions is an unresolved issue, the computation of the SIF of plate is a very interesting problem. In this paper, the LEFM analysis of plates subjected to out-of-plane loading is considered using the p-version of the finite element method. The superiority of the p-version of the finite element method in LEFM computations was demonstrated by Mehta [13], Basu et al. [14], and Woo and Basu [15]. It was shown by numerical experimentations and analytical proof [16] that the rate of convergence of p-extension is twice the rate of h-extension in the finite element modeling of plane stress/strain problems with sharp cracks when the number of degrees of freedom is increased by uniform or quasi-uniform mesh refinement. The objective of this study is to determine the SIF for the finite plate with different ratios of thickness to crack length (h/a), crack length to width (a/W) and boundary conditions under distributed edge moments and transverse pressures by the p-version crack model based on the virtual crack extension method and integrals of Legendre polynomials.
2. STRESS FIELD NEAR CRACK TIP AND STRESS INTENSITY FACTORS Irwin [1] obtained the form of the elastic stress distribution in the vicinity of a crack tip in extensional problems. Williams [2] extended his analysis to thin plates subjected to bending out of the plane. In each case it is shown that the significant stresses in the vicinity of the crack tip are those associated with the singularity of stress of the order r- ~/2,where r is the radial distance from the crack tip. Moreover, the distribution of stress in each, extension or bending, is unique, i.e. its functional form in terms of coordinates measured from the crack tip is always the same. Hence, for bending, Williams' results will be modified to define moment intensity factors in bending in a manner consistent with Irwin's definitions.
12Mrz ~F
O"0
=
30
3 + 5v
< LC°S
<>'r -
h3
12M,oz
r,0- ~ h
7+v
07 3KIGz cosmic+
3o ,+3v cos
-cos 2
{[ - s i n
+
7+ v
2J
,_,,
sin~ ~ +
~
F 30 L-sin 2
3 -q- 5v . 07 3KnGz ] 5 + 3v s m 2 ] - - - ~ - ;
[ 3o .o]3. - -oz}
+ sin ~- + sin
o}3,<,oz [ - c o s 302 l_,, + 3v cos ~
,
(1)
where z is a coordinate perpendicular to the middle plane of the plate, G and h are associated with the shear modulus of elasticity and thickness of the plate, respectively. When the crack is tilted at an angle /7 with reference to the plane about which a bending moment of magnitude Mo is applied, Sih et al. [3], and Hartranft and Sih [5] have provided the moment intensity factors shown in eqs (2)-(4).
495
p-Version finite element approximations of S1Fs
I ' e 0x~ X x ~ l . ~ l" ~
Consistent Range of Energy Release Rate
,.m p-level=7
~ I,~ I
.~0
a/W=0.1
\
h/a=0.1
\
\
,
IrV:-5
,
,
,
I(~i-O le~-7
6a/a
Fig. I. Sensitivity test of
,
W~:-8
\
I e'~-g
G(s) and 6a/a.
K~ = q) (1)M ox/~ sin 2 fl
(2)
Kn = ~(l)Mox//a sin fl cos fl
(3)
Km=
(1x+/ ~v)h D(1)Mox/a sin fl cos fl,
(4)
where the functions ~(1), ~(1) and f~(1) are non-dimensional bending stress intensity factors computed numerically from integral equations. Their values are a function of h/ax/1-O for different Poisson's ratios. In the case of mode I, the moment intensity factor K~ can be modified to the bending stress intensity factor k~ under the plane strain condition such that k, = *(l)abx/a.
(5)
where 6
(6)
ab=~SMo. 3. HIERARCHICAL C°-PLATE E L E M E N T S
3.1. Integrals of Legendre shape functions The general form of the integral of Legendre polynomials is defined over the standard domain such that
N[ 2
, ej(t) dt,
Table 1. Convergence of the maximum transverse displacement; Simply supported (S.S.)
(7)
w, x 10-2
Clamped condition (C.S.)
p-Level
Point load
Bending moment
Uniform load
Point load
4
0.48053 (5.22%)
0.04356 (--8.41%)
0.02728 (80.25%)
0.04830 (82.13%)
5
0.49163 (3.03%)
0.03979 (-0.9%)
0.09510 (31.16%)
0.15580 (36.44%)
6
0.49470 (1.79%)
0.04045 (-0.67%)
0.14036 ( - 1.60%)
0.23716 (3.25%)
7
0.50255 (0.87%)
0.04013 (0.12%)
0.13808 ( - 0.04%)
0.23995 (2.11%)
8
0.50482 (0.43%) 0.50699
0.04027 (-0.22%) 0.04018
0.13856 (0.30%) 0.13814
0.24241 (I.11%) 0.24513
Theory [17]
496
KWANG SUNG WOO and CHAE GUE LEE
x
@
®
®
®
2t.. Ld
,J
r"
2W
"1 w
Fig. 2. p-Version finite element model. where 1
Pj( t ) = ~
~ut( t 2- 1)i,
(8)
where i = 0, 1, 2 . . . . . For a standard quadrilateral domain, with 3, r / ( - 1 ~< ~ ~< 1 and - 1 ~< r/~< 1) as the coordinate of a point, the shape functions are built from the one-dimensional shape functions Fi(~) as follows [14, 15]: (1) The four vertex modes are F~(~) • Fj(r/) with i,j = 1, 2. (2) The side modes can be obtained by multiplying Fp(~) by (q + 1) and (q - 1) for sides, r / = + 1 and Fp(~/) by (3 + 1) and (3 - 1) for the sides ~ = + 1. Here p t> 2 is the order of the shape function, i.e. p = 2 indicates the quadratic shape functions. (3) The internal modes are valid for p / > 4 only and can be obtained by taking the product F~(~). Fj(rl) so that i + j = p , and both i and j are greater than or equal to 2. 3.2. Formulation of C°-plate element For the description of the Reissner-Mindlin plate problem, let E be Young's modulus, v the Poisson's ratio, G the shear modulus and h the thickness of the plate. We denote the deflection at a point (x, y) by co, the rotation about the y-axis by 0x and the rotation about the x-axis by 0y. The generalized strains E = (Z x, @~) of the plate are composed of a strain vector Z=
[
ax'
By'
\ay +ax]j
(9)
t.2
Alwar
1.0
k•0.6 •
0.0
p-version(4Elements, p=8)
0.4
0,Z
0"8.0
I
i
i
i
0.6
1.0
1.6
Z.O
2.5
Fig. 3. Variation of @(1) of infinite plate when a/W = 0.03.
p-Version finite element approximations of SIFs
i.o
~
1.4
497
• •
/d2a ~ 4 h/'2a=2
•
hP..a= 1.5 h/2a= 1
k.~ I,t h,"2a=0.5
g
h#.a= 0.3
•:.,eo.o
hC2a=0.1 0,0
0"00.0
i
I
i
0.1
0.2
0.3
I
0.4 a/W
i
I
i
0.iS
D.O
0.7
Fig. 4. Variation of ~ ( l ) with different ratios of
h/a
0.9
a/W.
and
and a rotation vector cI)= - ~ x -
-Oy
(lO)
.
The generalized stress resultants a = (mx, my, mxy, qx, qy)r are grouped into two parts, Q = (qx, qy)~ with the form
the
bending
M
moments
(11)
M = (m,., m~., m.~:.) r
=Dr'z
and shear forces (12)
Q =D,'~,
(13)
where
Dr=D
.
v
1
0 0 -
°1
0 1-v
Ds=
S "
g-
['0
(14)
The elasticity matrix denoted by D and S can be expressed by Eh 3
D-
S-
12(1 - v a)
Eh
(15)
2(1 + v)c~'
where c~ is the shear factor.
1.4
• p-version(4 Elements, p=8)
1,2 ..................."
•.....
....................\
.a
1.O
0.8 I a/W=0.5
0.6
0.4 0.0
!
I
i
I
0.5
! .0
1.5
2,0
h/7_a Fig. 5. Comparison o f ~ ( l ) with finite clement solutions when a/W = 0.5.
498
K W A N G S U N G W O O and C H A E G U E LEE
Effect of Poisson's 1.1
Ratio
(a/W=O.l)
v=0.1 a v=0.2 a ,v = 0.3
0
=........~=------'''---=
= ~=0.4
0.7
r~1
" ~
a
,
ff
R
l
I
l
I
i
0"0.0
0.5
1.0
I.§
2.0
2.5
h/a4i'6 Fig. 6. Effect of Poisson's ratio with different thickness to crack length ratios when
4. C O M P U T A T I O N
a/W = 0.1.
OF S T R E S S I N T E N S I T Y F A C T O R S
The finite element method has been used by a number of investigators to determine elastic stress intensity factors for cracked bodies. The characteristic elastic square root singularity has been represented by the use of the virtual crack extension method in this work. For the virtual crack extension method, it may be shown in the absence of body forces.
fo
6 A ( s ) ds =
LG(s)
{u},
(16)
where G(s) is Griffith's energy and a function of position s along the crack front, a is the length of the crack front and {u } is the vector of nodal point displacements found from the finite element computation. The change in the stiffness matrix A[K], for a given virtual crack extension may be written as a forward difference, i.e. A[K] ~ [K]a+,o -- [K],.
(17)
Consider a crack of length a which advances by an incremental amount 6a, thereby causing a release of strain energy of amount 6U. Thus the incremental crack surface 6A (s) for axisymmetric bodies is defined by: (18)
6 A ( s ) = h(a + 6a) - ha,
where h is the thickness of the plate. I
1.51
• . . 0 . . . . . . . 0 . . . . . . . . . . . . . . . . . 0 a / W = 0.6 ..0""'"
1.3 ...."
.. O'
O'
o _,"
O.g
.....
.-"
~o,~,o o
'"
...0"
•
~.... .."'
..o'
.,
"'"
.0""
.................
0 . . . . O .......
. .0""
...0 .......
. ....o. .... ,.0"" _ ...o 0"""
0
a/W= 0.I
0 ................. 0 a/W=oa
^ .......... ....... v-.
....o-""i o- ....... o ....... o'ii I
o adW=o.3 - a/W=o.2
.. • ,.,,a
........ oa/W=0.z
.... . ~.::..... g.;.......o~ • v
i:~::?:g:f ~0.7 r
o
Hartranft p-version(4 Element,p=8)
0.5
O'u.O
i
i
i
!
i
0.5
1.0
1.6
20.
E.6
3.0
Na4T6 Fig. 7. Variation of q~(1) associated with
h/ax/~
and
a/W.
p-Version finite element approximations of SIFs
=
1.4
• p-version(4 Element,p=8)
1.2
o Chen A Kang
499
Simply Sul~g}orted
0.6
0.8
•
Clamped Condition 0., o.2 0.0 0.0
' 0.1
' 0.2
' 0.3
' 0.4
' 0.5
' 0.0
' 0.7
' 0.8
• 0.9
1.0
a/W
Fig. 8. Normalized bending stress intensity factors vs a/W ratio with different boundary conditions. Therefore, the strain energy release rate G (s) for axisymmetric cracked bodies can be expressed as: 6U
G(s) = ~A (s) '
(19)
then the stress intensity factors are directly related to the value of G(s) caused by a crack extension in the appropriate mode. Hartranft and Sih [5] have provided k t expressions under plane strain condition including the effect of the plate thickness at the upper surface of the cracked plate (z = h/2). The relationship between k, and G for plane strain condition is expressed by: kl =
~f GE
(20/
l__V2,
in which G is the strain energy release rate under mode I action, v is Poisson's ratio and E is Young's modulus. However, since the quantity G(s) is very sensitive to the crack length increment 6a, the sensitivity test was investigated to find the consistent range of energy release rate between G(s) and 6a (Fig. 1). F r o m this, cSa was adapted by 10-7a in this study. 5. N U M E R I C A L R E S U L T S
5.1. Uncracked plates under transverse loads and bending moments The first example is an analysis of a square uncracked plate with length W and thickness h which is adopted as an example to verify the validity of C°-hierarchical plate element proposed in 1.4
~ 1.2
•
v:0.1
o v=0.2
A V=0.3
Simply Supported
~" 1.0 =
0.8
0.0 Con
on
0.8
0.7
0.4 0.2 0.a 0.0
0.!
0.2
0.3
0.4
0.5
0.8
0.0
!.0
a/W Fig. 9. Variation of normalized bending stress intensity factors with a/W and Poisson's ratio.
500
KWANG SUNG WOO and CHAE GUE LEE
Po
m
2W
Fig. 10. Configurationof cracked square plate with two opposite edges simply supported and the other edges free. this work before proceeding with the cracked plate problems. Due to symmetry, a quarter plate is modeled with only the single element. The geometric and mechanical data of the problem are: W = 10, h = 0.1, E = 107 and v = 0.3. Three different types of loading are applied, e.g. transversly uniform pressures with unitary intensity, a concentrated load with P = 4 at the center of the plate and constant distributed bending moment Mo per unit length. Various boundary conditions are considered in order to analyze the sensitivity of the finite solutions for simply supported (S.S.) and clamped conditions (C.S.). Table 1 shows the results from the convergence test of the center transverse displacement of a thin plate (W/h = 100) compared with the thin plate solution by Timoshenko and Woinowsky-Krieger, for the four combinations of loading-boundary conditions. In all the cases, the results by the p-version model exhibit good convergence behavior within 1.0% error when the p-level is 8. 5.2. A centrally cracked square plate under uniform edge moments The second example is a centrally cracked square plate subjected to uniform edge normal bending moments. Only a quarter of the structure is discretized into 2 × 2 mesh, because of symmetry (Fig. 2). The geometric and mechanical data of the problem are given by the length of the square plate with W = 10, E = 1 × 106, V ----- 0.3 and constant normal bending moment Mo per unit length. Since the virtual crack extension method is applied to calculate the stress intensity factor in this study, it is necessary to determine the crack length increment 6a. As mentioned earlier, the energy release rate G(s) does not show consistent values whatever crack length increment 6a is used. Thus, the crack length increment 6a has been adapted by 10-Ta from the sensitivity test associated with G(s) and 6a/a in Fig. 1. The p-version results of the non-dimensional bending stress intensity factor q~(1) defined by eq. (5) have been compared with three-dimensional finite element solutions by Alwar and Nambissan [10], and theoretical solutions of the infinite cracked plates obtained by Hartranft and Sih [5] when a/W = 0.03 which are plotted in Fig. 3. It is noted that the p-version four-element model with p = 8, in the case of the infinite plate, which is almost equal to Alwar's analysis, shows a good comparison with Hartranft's theoretical solution within 4% relative errors. The non-dimensional bending stress intensity factors q~(1) for the finite plate with different ratios of a/W have been shown in Fig. 4 as the thickness of crack length ratio h/2a varies. The effects of shear deformation of the finite cracked plate have been tested with the same p-version model in Fig. 5. The p-version finite element model is very close to the three-dimensional analysis by Alwar and Nambissan[10], however, two-dimensional analysis by Yagawa and Nishioka [9] gives large difference in reference to both the p-version model and Alwar's solution when a~ W = 0.5. Figure 6 represents the influence of the Poisson's ratio varied from 0.1 to 0.4 on the determination of non-dimensional bending stress intensity factors as the thickness is increased. As we are aware of it from this figure, the non-dimensional bending stress intensity factors have been affected by the variation of Poisson's ratio. Unfortunately, there are no published papers in the literature related to the effect of Poisson's ratio up to this stage. The non-dimensional bending stress intensity factors q~(1) for the finite plate with respect to the ratio h / a x / ~ have been shown in Fig. 7 as the ratio of a/W varies. The p-version four-
p-Version finite element approximations of SIFs
501
1.3 1.2 1.1 1.0
W/h= 10
/
• p-version(4 Element, p=8) o SOSA
0.9
~
0.7 0.8
o/
ooj"
O.fi 0.4 0.3
0.0
i
0.l
I
I
I
I
I
I
I
|
0.Z
0.3
0.4
0.5
0.O
0.7
0.8
0.g
.0
o/w Fig. I |. G r o w t h o f normalized bending stress intensity factors with
a/W ratio.
element model shows an excellent agreement with the infinite plate solution by Hartranft when a~ W is zero. 5.3. A centrally cracked square plate under uniform pressures The third example has been tested to show the effect of bending stress intensity factors with different crack lengths and boundary conditions in the case of a centrally cracked square plate under uniform transverse loading Po. The finite element mesh and geometric/mechanical data are the same as a centrally cracked square plate problem under distributed edge moments. The bending stress intensity factors with respect to crack length for simply supported (S.S.) and clamped condition (C.S.) have been investigated in Fig. 8 where the stress intensity factors were normalized by the reference value of the bending stress intensity factor of a / W = 0.5 with simply supported condition and this problem was also solved under the condition of thin cracked plate, e.g. W/h = 30. The p-version finite element model shows a good comparison with the hybrid~lisplacement singular element approach by Chen and Chen [11] and the hybrid mongrel element model by Kang and De Saxce [12], respectively. It is also noted that the bending stress intensity factor increases to a maximum value at a certain ratio of a~ W and begins to decrease even when the crack length is growing. As we expected, the effect of Poisson's ratio is very large as well as the centrally cracked plate under uniform edge moments. However, it is observed that Poisson's ratio varied from 0.1 to 0.3 has an effect on the normalized bending stress intensity factor in S.S. boundary condition rather than in C.S. boundary condition (Fig. 9). The configuration of the square cracked plate with two opposite edges simply supported and the other two edges free is shown in Fig. 10 where the loading condition and mechanical properties are unchanged like the previous problem. The solutions obtained by the four-element p-version model with p = 8 are in good accordance with those by Sosa and Eichen [18]. In this case, the non-dimensional bending stress intensity factor grows rapidly as the crack length increases (Fig. ll).
6. SUMMARY AND CONCLUSIONS The new p-version crack model which has been proposed is based on the Reissner-Mindlin theory, integral of Legendre polynomials and virtual crack extension method. Also, this model is under the scope of theory of small scale yielding. Throughout treating the examples, the cracked square plates have been investigated with respect to the effects of the size of finite plate, Poisson's ratio, loading conditions of uniform edge moments and transverse pressures, various boundary conditions and variation of thickness. The physical domains are modeled by only four C°-hierar chical plate elements with p-level = 8 to get the same levels of accuracy from the special singular finite element solutions in the literature. From this study, it is concluded that the proposed crack
502
KWANG SUNG WOO and CHAE GUE LEE
m o d e l s h o w s s u p e r i o r p e r f o r m a n c e to the existing c r a c k m o d e l o n the basis o f h y b r i d s i n g u l a r e l e m e n t s in the sense o f m o d e l i n g simplicity, a c c u r a c y a n d a t r e m e n d o u s s a v i n g in C P U a n d u s e r ' s time. A l s o , it is a p p a r e n t t h a t p - v e r s i o n c r a c k m o d e l is v e r y s u i t a b l e f o r L E F M analysis i r r e s p e c t i v e of loadings and geometric/mechanical properties.
REFERENCES [1] G. R. Irwin, Analysis of stresses and strains near the end of a crack traversing a plate. J. appl. Mech. 79, 361-364 (1957). [2] M. L. Williams, The bending stress distribution at the base of a stationary crack. J. appl. Mech. 28, 78-82 (1961). [3] G. C. Sih, P. C. Paris and F. Erdogan, Crack-tip, stress-intensity factors for plane extension and plate bending problems. J. appl. Mech. (Trans. of ASME), 306-312 (1962). [4] J. K. Knowles and N. M. Wang, On the bending of an elastic plate containing a crack. J. Math. Phys. 39, 223-236 (1960). [5] R. J. Hartranft and G. C. Sih, Effect of plate thickness on the bending stress distribution around through cracks. J. Math. Phys. 47, 276-291 (1968). 16] M. V. V. Murthy, K. M. Raju and S. Viswanath, On the bending stress distribution at the tip of a stationary crack for Reissner's theory. Int. J. Fracture 17, 537-552 (1981). [7] W. K. Wilson and D. G. Thompson, On the finite element method for calculating stress intensity factors for cracked plates in bending. Engng Fracture Mech. 3, 97-102 (1971). [8] P. D. Hilton, A special finite element approach for three-dimensional crack problems, in Mechanics of Fracture (Edited by G. C. Sih), Vol. 3. Noordhoff, Leyden (1977). [9] G. Yagawa and T. Nishioka, Finite element analysis of stress intensity factor for plate extension and plate bending problems. Int. J. numer. Meth. Engng 14, 727-740 (1979). [10] R. S. Alwar and K. N. R. Nambissan, Three-dimensional finite element analysis of cracked thick plates in bending. Int. J. numer. Meth. Engng 10, 551-564 (1983). [11] W. H. Chen and P. Y. Chen, A hybrid~lisplacement finite element model for the bending analysis of thin cracked plates. Int. J. Fracture 24, 83 106 (1986). [12] C. H. Kang and G. De Saxce, Computation of stress intensity factors for plate bending problem in fracture mechanics by hybrid mongrel finite element. Comput. Structures 42, 581-589 (1992). [13] A. K. Mehta, p-Convergent finite element approximations in linear elastic fracture mechanics. Doctoral Dissertation, Washington University, St. Louis, U.S.A. (1978). [14] P. K. Basu, N. U. Ahmed and K. S. Woo, Plates and shells with crack-like flaws. Computer Utilization in Structural Engineering (Edited by J. K. Nelson), Proc. 7th Structures Congress, ASCE, 286-295 (1989). [I 5] K. S. Woo and P. K. Basu, Analysis of singular cylindrical shells by P-version of finite element method. Int. J. Solids Structures 25, 151-165 (1989). [16] B. A. Szabo, Estimation and control of error based on p-convergence, in Accuracy Estimates and Adaptive Refinements in Finite Element Computation. Wiley (1986). [17] S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells. 2nd edn. McGraw-Hill (1969). [18] H. A. Sosa and J. W. Eichen, Computation of stress intensity factors for plate bending via a path-independent integral. Engng Fracture Mech. 25, 451-462 (1986). (Received 11 May 1994)
p-VERSION
Engineering Fracture Mechanics Vol. 52, No. 3, pp. 493-502, 1995 Copyright ~? 1995 ElsevierScience Ltd 0013-7944(95)00019--4 Printed in Great Britain. All rights reserved 0013-7944/95 $9.50+ 0.00
FINITE ELEMENT
APPROXIMATIONS
OF
STRESS INTENSITY FACTORS FOR CRACKED PLATES INCLUDING SHEAR DEFORMATION KWANG SUNG WOO Department of Civil Engineering, Yeungnam University, Gyongsan 712-749, Korea and CHAE GUE LEE Department of Civil Engineering, Chosun University, Kwangju 501-759, Korea Abstract--The new p-version crack model is proposed to estimate the bending stress intensity factors of the thick cracked plate under flexure. The proposed model is based on high order theory and hierarchical C°-plate element including shear deformation. The displacement fields are defined by integrals of Legendre polynomials which can be classified into three groups, e.g. basic, side and internal mode. The computer implementation allows arbitrary variations of p-level up to a maximum value of 10. The bending stress intensity factors are computed by the virtual crack extension approach. The effects of ratios of thickness to crack length (h/a), crack length to width (a/W) and boundary conditions are investigated. Very good agreement with the existing solution in the literature is shown for the uncracked plate as well as the cracked plate.
1. INTRODUCTION ONE OFTHE BASICrequirements for the strength analysis of plates containing flaws or cracks is the knowledge of the singular character of the stress field in the neighborhood of the point at the crack tip. The SIF (stress intensity factor), which is in fact the strength of the stress singularities at crack tips, is applied to predict the static strength of cracked bodies by Irwin [1] for plane extension, symmetric with respect to the crack. The stress intensity factors have been shown to control the rate of crack propagation under cyclic loading in such situations. Several investigators have discussed the nature of the local stresses around a sharp crack in a thin plate subjected to out-of-plane bending loads. Based on the Poisson-Kirchhoff theory of thin plate and the technique of Fadle eigen function expansion, Williams [2] found that the elastic bending stresses near the tip of a semi-infinite crack vary as the inverse square root of the radial distance from the crack front. His results were not complete in that the strength or magnitude of the local stresses was left undetermined. Sih et al. [3] cleared the way for finding the coefficients in the eigen function expansions by application of the theory of complex functions. The results in refs [2, 3] were obtained from the classical fourth-order theory of thin plates, however, the edge conditions at the crack surfaces are satisfied only in an approximate manner in that the three physically natural boundary conditions of prescribing bending moment, twisting moment and transverse shear stress are replaced by two conditions. Owing to such a replacement, the stress distribution in the immediate neighborhood of the crack edges will naturally be affected and will not be accurate. To overcome the aforementioned shortcoming, Knowles and Wang [4] used Reissner's sixth-order plate theory to satisfy all the physically natural boundary conditions along the crack edges with the help of singular integral equations. However, the extension of the method to cracked plates with finite plate thickness does not appear to be tractable. For the examination of the effect of plate thickness, Hartranft and Sih [5] put forward higher order theories that included the effect of transverse shear deformations and essentially solved the problem of infinite plates subjected to uniform moment. Murthy et al. [6] used Reissner's plate theory and applied collocation, and successive integration methods to satisfy the boundary conditions for a finite plate with a through crack. The associated numerical difficulties led them to conclude the finite element method is better suited for such problems. Wilson and Thompson [7] used the h-version of the finite element method 493
494
KWANG SUNG WOO and CHAE GUE LEE
to compute the SIF for the through-thickness crack in a plate under flexure. They used a large number of three-node triangular C~-elements (actually 257 and 305 node models) and computed the value of SIF using the direct method. Hilton [8] used 48 20-noded three-dimensional isoparametric elements including the special singularity element (279 node model) to compute a SIF. Yagawa and Nishioka[9] used a superposition method of analytical and finite element solutions on the basis of 36 eight-noded isoparametric elements (133 node model) together with the reduced integration technique. Alwar and Nambissan [10] used 110 three-dimensional quadratic isoparametric elements, including degenerate triangular quarter-point prism elements at the crack tip region, to model one-eighth of the plate and computed the SIF by COD method. The results presented by earlier workers [8, 9] for the finite plate under uniform moments differ significantly from the results by Alwar and Nambissan. They indicated that the difference may be due to the coarseness of the mesh used by those. Chen and Chen [11] presented the hybrid~lisplacement singular element model based on the Kirchhoff plate theory for the bending analysis of thin cracked infinite plate and transverse pressure analysis of thin finite cracked plate with various boundary conditions. Kang and De Saxce [12] used 35 hybrid mongrel elements, which is the special crack tip element, for infinite cracked plate under bending and finite thin plate under transverse pressure. Since the singularity at the through crack corner of finite plate subjected to various loads and boundary conditions is an unresolved issue, the computation of the SIF of plate is a very interesting problem. In this paper, the LEFM analysis of plates subjected to out-of-plane loading is considered using the p-version of the finite element method. The superiority of the p-version of the finite element method in LEFM computations was demonstrated by Mehta [13], Basu et al. [14], and Woo and Basu [15]. It was shown by numerical experimentations and analytical proof [16] that the rate of convergence of p-extension is twice the rate of h-extension in the finite element modeling of plane stress/strain problems with sharp cracks when the number of degrees of freedom is increased by uniform or quasi-uniform mesh refinement. The objective of this study is to determine the SIF for the finite plate with different ratios of thickness to crack length (h/a), crack length to width (a/W) and boundary conditions under distributed edge moments and transverse pressures by the p-version crack model based on the virtual crack extension method and integrals of Legendre polynomials.
2. STRESS FIELD NEAR CRACK TIP AND STRESS INTENSITY FACTORS Irwin [1] obtained the form of the elastic stress distribution in the vicinity of a crack tip in extensional problems. Williams [2] extended his analysis to thin plates subjected to bending out of the plane. In each case it is shown that the significant stresses in the vicinity of the crack tip are those associated with the singularity of stress of the order r- ~/2,where r is the radial distance from the crack tip. Moreover, the distribution of stress in each, extension or bending, is unique, i.e. its functional form in terms of coordinates measured from the crack tip is always the same. Hence, for bending, Williams' results will be modified to define moment intensity factors in bending in a manner consistent with Irwin's definitions.
12Mrz ~F
O"0
=
30
3 + 5v
< LC°S
<>'r -
h3
12M,oz
r,0- ~ h
7+v
07 3KIGz cosmic+
3o ,+3v cos
-cos 2
{[ - s i n
+
7+ v
2J
,_,,
sin~ ~ +
~
F 30 L-sin 2
3 -q- 5v . 07 3KnGz ] 5 + 3v s m 2 ] - - - ~ - ;
[ 3o .o]3. - -oz}
+ sin ~- + sin
o}3,<,oz [ - c o s 302 l_,, + 3v cos ~
,
(1)
where z is a coordinate perpendicular to the middle plane of the plate, G and h are associated with the shear modulus of elasticity and thickness of the plate, respectively. When the crack is tilted at an angle /7 with reference to the plane about which a bending moment of magnitude Mo is applied, Sih et al. [3], and Hartranft and Sih [5] have provided the moment intensity factors shown in eqs (2)-(4).
495
p-Version finite element approximations of S1Fs
I ' e 0x~ X x ~ l . ~ l" ~
Consistent Range of Energy Release Rate
,.m p-level=7
~ I,~ I
.~0
a/W=0.1
\
h/a=0.1
\
\
,
IrV:-5
,
,
,
I(~i-O le~-7
6a/a
Fig. I. Sensitivity test of
,
W~:-8
\
I e'~-g
G(s) and 6a/a.
K~ = q) (1)M ox/~ sin 2 fl
(2)
Kn = ~(l)Mox//a sin fl cos fl
(3)
Km=
(1x+/ ~v)h D(1)Mox/a sin fl cos fl,
(4)
where the functions ~(1), ~(1) and f~(1) are non-dimensional bending stress intensity factors computed numerically from integral equations. Their values are a function of h/ax/1-O for different Poisson's ratios. In the case of mode I, the moment intensity factor K~ can be modified to the bending stress intensity factor k~ under the plane strain condition such that k, = *(l)abx/a.
(5)
where 6
(6)
ab=~SMo. 3. HIERARCHICAL C°-PLATE E L E M E N T S
3.1. Integrals of Legendre shape functions The general form of the integral of Legendre polynomials is defined over the standard domain such that
N[ 2
, ej(t) dt,
Table 1. Convergence of the maximum transverse displacement; Simply supported (S.S.)
(7)
w, x 10-2
Clamped condition (C.S.)
p-Level
Point load
Bending moment
Uniform load
Point load
4
0.48053 (5.22%)
0.04356 (--8.41%)
0.02728 (80.25%)
0.04830 (82.13%)
5
0.49163 (3.03%)
0.03979 (-0.9%)
0.09510 (31.16%)
0.15580 (36.44%)
6
0.49470 (1.79%)
0.04045 (-0.67%)
0.14036 ( - 1.60%)
0.23716 (3.25%)
7
0.50255 (0.87%)
0.04013 (0.12%)
0.13808 ( - 0.04%)
0.23995 (2.11%)
8
0.50482 (0.43%) 0.50699
0.04027 (-0.22%) 0.04018
0.13856 (0.30%) 0.13814
0.24241 (I.11%) 0.24513
Theory [17]
496
KWANG SUNG WOO and CHAE GUE LEE
x
@
®
®
®
2t.. Ld
,J
r"
2W
"1 w
Fig. 2. p-Version finite element model. where 1
Pj( t ) = ~
~ut( t 2- 1)i,
(8)
where i = 0, 1, 2 . . . . . For a standard quadrilateral domain, with 3, r / ( - 1 ~< ~ ~< 1 and - 1 ~< r/~< 1) as the coordinate of a point, the shape functions are built from the one-dimensional shape functions Fi(~) as follows [14, 15]: (1) The four vertex modes are F~(~) • Fj(r/) with i,j = 1, 2. (2) The side modes can be obtained by multiplying Fp(~) by (q + 1) and (q - 1) for sides, r / = + 1 and Fp(~/) by (3 + 1) and (3 - 1) for the sides ~ = + 1. Here p t> 2 is the order of the shape function, i.e. p = 2 indicates the quadratic shape functions. (3) The internal modes are valid for p / > 4 only and can be obtained by taking the product F~(~). Fj(rl) so that i + j = p , and both i and j are greater than or equal to 2. 3.2. Formulation of C°-plate element For the description of the Reissner-Mindlin plate problem, let E be Young's modulus, v the Poisson's ratio, G the shear modulus and h the thickness of the plate. We denote the deflection at a point (x, y) by co, the rotation about the y-axis by 0x and the rotation about the x-axis by 0y. The generalized strains E = (Z x, @~) of the plate are composed of a strain vector Z=
[
ax'
By'
\ay +ax]j
(9)
t.2
Alwar
1.0
k•0.6 •
0.0
p-version(4Elements, p=8)
0.4
0,Z
0"8.0
I
i
i
i
0.6
1.0
1.6
Z.O
2.5
Fig. 3. Variation of @(1) of infinite plate when a/W = 0.03.
p-Version finite element approximations of SIFs
i.o
~
1.4
497
• •
/d2a ~ 4 h/'2a=2
•
hP..a= 1.5 h/2a= 1
k.~ I,t h,"2a=0.5
g
h#.a= 0.3
•:.,eo.o
hC2a=0.1 0,0
0"00.0
i
I
i
0.1
0.2
0.3
I
0.4 a/W
i
I
i
0.iS
D.O
0.7
Fig. 4. Variation of ~ ( l ) with different ratios of
h/a
0.9
a/W.
and
and a rotation vector cI)= - ~ x -
-Oy
(lO)
.
The generalized stress resultants a = (mx, my, mxy, qx, qy)r are grouped into two parts, Q = (qx, qy)~ with the form
the
bending
M
moments
(11)
M = (m,., m~., m.~:.) r
=Dr'z
and shear forces (12)
Q =D,'~,
(13)
where
Dr=D
.
v
1
0 0 -
°1
0 1-v
Ds=
S "
g-
['0
(14)
The elasticity matrix denoted by D and S can be expressed by Eh 3
D-
S-
12(1 - v a)
Eh
(15)
2(1 + v)c~'
where c~ is the shear factor.
1.4
• p-version(4 Elements, p=8)
1,2 ..................."
•.....
....................\
.a
1.O
0.8 I a/W=0.5
0.6
0.4 0.0
!
I
i
I
0.5
! .0
1.5
2,0
h/7_a Fig. 5. Comparison o f ~ ( l ) with finite clement solutions when a/W = 0.5.
498
K W A N G S U N G W O O and C H A E G U E LEE
Effect of Poisson's 1.1
Ratio
(a/W=O.l)
v=0.1 a v=0.2 a ,v = 0.3
0
=........~=------'''---=
= ~=0.4
0.7
r~1
" ~
a
,
ff
R
l
I
l
I
i
0"0.0
0.5
1.0
I.§
2.0
2.5
h/a4i'6 Fig. 6. Effect of Poisson's ratio with different thickness to crack length ratios when
4. C O M P U T A T I O N
a/W = 0.1.
OF S T R E S S I N T E N S I T Y F A C T O R S
The finite element method has been used by a number of investigators to determine elastic stress intensity factors for cracked bodies. The characteristic elastic square root singularity has been represented by the use of the virtual crack extension method in this work. For the virtual crack extension method, it may be shown in the absence of body forces.
fo
6 A ( s ) ds =
LG(s)
{u},
(16)
where G(s) is Griffith's energy and a function of position s along the crack front, a is the length of the crack front and {u } is the vector of nodal point displacements found from the finite element computation. The change in the stiffness matrix A[K], for a given virtual crack extension may be written as a forward difference, i.e. A[K] ~ [K]a+,o -- [K],.
(17)
Consider a crack of length a which advances by an incremental amount 6a, thereby causing a release of strain energy of amount 6U. Thus the incremental crack surface 6A (s) for axisymmetric bodies is defined by: (18)
6 A ( s ) = h(a + 6a) - ha,
where h is the thickness of the plate. I
1.51
• . . 0 . . . . . . . 0 . . . . . . . . . . . . . . . . . 0 a / W = 0.6 ..0""'"
1.3 ...."
.. O'
O'
o _,"
O.g
.....
.-"
~o,~,o o
'"
...0"
•
~.... .."'
..o'
.,
"'"
.0""
.................
0 . . . . O .......
. .0""
...0 .......
. ....o. .... ,.0"" _ ...o 0"""
0
a/W= 0.I
0 ................. 0 a/W=oa
^ .......... ....... v-.
....o-""i o- ....... o ....... o'ii I
o adW=o.3 - a/W=o.2
.. • ,.,,a
........ oa/W=0.z
.... . ~.::..... g.;.......o~ • v
i:~::?:g:f ~0.7 r
o
Hartranft p-version(4 Element,p=8)
0.5
O'u.O
i
i
i
!
i
0.5
1.0
1.6
20.
E.6
3.0
Na4T6 Fig. 7. Variation of q~(1) associated with
h/ax/~
and
a/W.
p-Version finite element approximations of SIFs
=
1.4
• p-version(4 Element,p=8)
1.2
o Chen A Kang
499
Simply Sul~g}orted
0.6
0.8
•
Clamped Condition 0., o.2 0.0 0.0
' 0.1
' 0.2
' 0.3
' 0.4
' 0.5
' 0.0
' 0.7
' 0.8
• 0.9
1.0
a/W
Fig. 8. Normalized bending stress intensity factors vs a/W ratio with different boundary conditions. Therefore, the strain energy release rate G (s) for axisymmetric cracked bodies can be expressed as: 6U
G(s) = ~A (s) '
(19)
then the stress intensity factors are directly related to the value of G(s) caused by a crack extension in the appropriate mode. Hartranft and Sih [5] have provided k t expressions under plane strain condition including the effect of the plate thickness at the upper surface of the cracked plate (z = h/2). The relationship between k, and G for plane strain condition is expressed by: kl =
~f GE
(20/
l__V2,
in which G is the strain energy release rate under mode I action, v is Poisson's ratio and E is Young's modulus. However, since the quantity G(s) is very sensitive to the crack length increment 6a, the sensitivity test was investigated to find the consistent range of energy release rate between G(s) and 6a (Fig. 1). F r o m this, cSa was adapted by 10-7a in this study. 5. N U M E R I C A L R E S U L T S
5.1. Uncracked plates under transverse loads and bending moments The first example is an analysis of a square uncracked plate with length W and thickness h which is adopted as an example to verify the validity of C°-hierarchical plate element proposed in 1.4
~ 1.2
•
v:0.1
o v=0.2
A V=0.3
Simply Supported
~" 1.0 =
0.8
0.0 Con
on
0.8
0.7
0.4 0.2 0.a 0.0
0.!
0.2
0.3
0.4
0.5
0.8
0.0
!.0
a/W Fig. 9. Variation of normalized bending stress intensity factors with a/W and Poisson's ratio.
500
KWANG SUNG WOO and CHAE GUE LEE
Po
m
2W
Fig. 10. Configurationof cracked square plate with two opposite edges simply supported and the other edges free. this work before proceeding with the cracked plate problems. Due to symmetry, a quarter plate is modeled with only the single element. The geometric and mechanical data of the problem are: W = 10, h = 0.1, E = 107 and v = 0.3. Three different types of loading are applied, e.g. transversly uniform pressures with unitary intensity, a concentrated load with P = 4 at the center of the plate and constant distributed bending moment Mo per unit length. Various boundary conditions are considered in order to analyze the sensitivity of the finite solutions for simply supported (S.S.) and clamped conditions (C.S.). Table 1 shows the results from the convergence test of the center transverse displacement of a thin plate (W/h = 100) compared with the thin plate solution by Timoshenko and Woinowsky-Krieger, for the four combinations of loading-boundary conditions. In all the cases, the results by the p-version model exhibit good convergence behavior within 1.0% error when the p-level is 8. 5.2. A centrally cracked square plate under uniform edge moments The second example is a centrally cracked square plate subjected to uniform edge normal bending moments. Only a quarter of the structure is discretized into 2 × 2 mesh, because of symmetry (Fig. 2). The geometric and mechanical data of the problem are given by the length of the square plate with W = 10, E = 1 × 106, V ----- 0.3 and constant normal bending moment Mo per unit length. Since the virtual crack extension method is applied to calculate the stress intensity factor in this study, it is necessary to determine the crack length increment 6a. As mentioned earlier, the energy release rate G(s) does not show consistent values whatever crack length increment 6a is used. Thus, the crack length increment 6a has been adapted by 10-Ta from the sensitivity test associated with G(s) and 6a/a in Fig. 1. The p-version results of the non-dimensional bending stress intensity factor q~(1) defined by eq. (5) have been compared with three-dimensional finite element solutions by Alwar and Nambissan [10], and theoretical solutions of the infinite cracked plates obtained by Hartranft and Sih [5] when a/W = 0.03 which are plotted in Fig. 3. It is noted that the p-version four-element model with p = 8, in the case of the infinite plate, which is almost equal to Alwar's analysis, shows a good comparison with Hartranft's theoretical solution within 4% relative errors. The non-dimensional bending stress intensity factors q~(1) for the finite plate with different ratios of a/W have been shown in Fig. 4 as the thickness of crack length ratio h/2a varies. The effects of shear deformation of the finite cracked plate have been tested with the same p-version model in Fig. 5. The p-version finite element model is very close to the three-dimensional analysis by Alwar and Nambissan[10], however, two-dimensional analysis by Yagawa and Nishioka [9] gives large difference in reference to both the p-version model and Alwar's solution when a~ W = 0.5. Figure 6 represents the influence of the Poisson's ratio varied from 0.1 to 0.4 on the determination of non-dimensional bending stress intensity factors as the thickness is increased. As we are aware of it from this figure, the non-dimensional bending stress intensity factors have been affected by the variation of Poisson's ratio. Unfortunately, there are no published papers in the literature related to the effect of Poisson's ratio up to this stage. The non-dimensional bending stress intensity factors q~(1) for the finite plate with respect to the ratio h / a x / ~ have been shown in Fig. 7 as the ratio of a/W varies. The p-version four-
p-Version finite element approximations of SIFs
501
1.3 1.2 1.1 1.0
W/h= 10
/
• p-version(4 Element, p=8) o SOSA
0.9
~
0.7 0.8
o/
ooj"
O.fi 0.4 0.3
0.0
i
0.l
I
I
I
I
I
I
I
|
0.Z
0.3
0.4
0.5
0.O
0.7
0.8
0.g
.0
o/w Fig. I |. G r o w t h o f normalized bending stress intensity factors with
a/W ratio.
element model shows an excellent agreement with the infinite plate solution by Hartranft when a~ W is zero. 5.3. A centrally cracked square plate under uniform pressures The third example has been tested to show the effect of bending stress intensity factors with different crack lengths and boundary conditions in the case of a centrally cracked square plate under uniform transverse loading Po. The finite element mesh and geometric/mechanical data are the same as a centrally cracked square plate problem under distributed edge moments. The bending stress intensity factors with respect to crack length for simply supported (S.S.) and clamped condition (C.S.) have been investigated in Fig. 8 where the stress intensity factors were normalized by the reference value of the bending stress intensity factor of a / W = 0.5 with simply supported condition and this problem was also solved under the condition of thin cracked plate, e.g. W/h = 30. The p-version finite element model shows a good comparison with the hybrid~lisplacement singular element approach by Chen and Chen [11] and the hybrid mongrel element model by Kang and De Saxce [12], respectively. It is also noted that the bending stress intensity factor increases to a maximum value at a certain ratio of a~ W and begins to decrease even when the crack length is growing. As we expected, the effect of Poisson's ratio is very large as well as the centrally cracked plate under uniform edge moments. However, it is observed that Poisson's ratio varied from 0.1 to 0.3 has an effect on the normalized bending stress intensity factor in S.S. boundary condition rather than in C.S. boundary condition (Fig. 9). The configuration of the square cracked plate with two opposite edges simply supported and the other two edges free is shown in Fig. 10 where the loading condition and mechanical properties are unchanged like the previous problem. The solutions obtained by the four-element p-version model with p = 8 are in good accordance with those by Sosa and Eichen [18]. In this case, the non-dimensional bending stress intensity factor grows rapidly as the crack length increases (Fig. ll).
6. SUMMARY AND CONCLUSIONS The new p-version crack model which has been proposed is based on the Reissner-Mindlin theory, integral of Legendre polynomials and virtual crack extension method. Also, this model is under the scope of theory of small scale yielding. Throughout treating the examples, the cracked square plates have been investigated with respect to the effects of the size of finite plate, Poisson's ratio, loading conditions of uniform edge moments and transverse pressures, various boundary conditions and variation of thickness. The physical domains are modeled by only four C°-hierar chical plate elements with p-level = 8 to get the same levels of accuracy from the special singular finite element solutions in the literature. From this study, it is concluded that the proposed crack
502
KWANG SUNG WOO and CHAE GUE LEE
m o d e l s h o w s s u p e r i o r p e r f o r m a n c e to the existing c r a c k m o d e l o n the basis o f h y b r i d s i n g u l a r e l e m e n t s in the sense o f m o d e l i n g simplicity, a c c u r a c y a n d a t r e m e n d o u s s a v i n g in C P U a n d u s e r ' s time. A l s o , it is a p p a r e n t t h a t p - v e r s i o n c r a c k m o d e l is v e r y s u i t a b l e f o r L E F M analysis i r r e s p e c t i v e of loadings and geometric/mechanical properties.
REFERENCES [1] G. R. Irwin, Analysis of stresses and strains near the end of a crack traversing a plate. J. appl. Mech. 79, 361-364 (1957). [2] M. L. Williams, The bending stress distribution at the base of a stationary crack. J. appl. Mech. 28, 78-82 (1961). [3] G. C. Sih, P. C. Paris and F. Erdogan, Crack-tip, stress-intensity factors for plane extension and plate bending problems. J. appl. Mech. (Trans. of ASME), 306-312 (1962). [4] J. K. Knowles and N. M. Wang, On the bending of an elastic plate containing a crack. J. Math. Phys. 39, 223-236 (1960). [5] R. J. Hartranft and G. C. Sih, Effect of plate thickness on the bending stress distribution around through cracks. J. Math. Phys. 47, 276-291 (1968). 16] M. V. V. Murthy, K. M. Raju and S. Viswanath, On the bending stress distribution at the tip of a stationary crack for Reissner's theory. Int. J. Fracture 17, 537-552 (1981). [7] W. K. Wilson and D. G. Thompson, On the finite element method for calculating stress intensity factors for cracked plates in bending. Engng Fracture Mech. 3, 97-102 (1971). [8] P. D. Hilton, A special finite element approach for three-dimensional crack problems, in Mechanics of Fracture (Edited by G. C. Sih), Vol. 3. Noordhoff, Leyden (1977). [9] G. Yagawa and T. Nishioka, Finite element analysis of stress intensity factor for plate extension and plate bending problems. Int. J. numer. Meth. Engng 14, 727-740 (1979). [10] R. S. Alwar and K. N. R. Nambissan, Three-dimensional finite element analysis of cracked thick plates in bending. Int. J. numer. Meth. Engng 10, 551-564 (1983). [11] W. H. Chen and P. Y. Chen, A hybrid~lisplacement finite element model for the bending analysis of thin cracked plates. Int. J. Fracture 24, 83 106 (1986). [12] C. H. Kang and G. De Saxce, Computation of stress intensity factors for plate bending problem in fracture mechanics by hybrid mongrel finite element. Comput. Structures 42, 581-589 (1992). [13] A. K. Mehta, p-Convergent finite element approximations in linear elastic fracture mechanics. Doctoral Dissertation, Washington University, St. Louis, U.S.A. (1978). [14] P. K. Basu, N. U. Ahmed and K. S. Woo, Plates and shells with crack-like flaws. Computer Utilization in Structural Engineering (Edited by J. K. Nelson), Proc. 7th Structures Congress, ASCE, 286-295 (1989). [I 5] K. S. Woo and P. K. Basu, Analysis of singular cylindrical shells by P-version of finite element method. Int. J. Solids Structures 25, 151-165 (1989). [16] B. A. Szabo, Estimation and control of error based on p-convergence, in Accuracy Estimates and Adaptive Refinements in Finite Element Computation. Wiley (1986). [17] S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells. 2nd edn. McGraw-Hill (1969). [18] H. A. Sosa and J. W. Eichen, Computation of stress intensity factors for plate bending via a path-independent integral. Engng Fracture Mech. 25, 451-462 (1986). (Received 11 May 1994)