Stress analysis around circular holes in FRP laminates under transverse load

Compurrrs & Srrummu Printed in Great Lhtain.

Vol.33.

No. 4. pp. 929-937,

lM45-7949/89 %3.00 + 0.00 \:‘ 1989 Pergamon Pressplc

1989

STRESS ANALYSIS AROUND CIRCULAR HOLES IN FRP LAMINATES UNDER TRANSVERSE LOAD T. K. PAUL? and K. M. RAO Department of Mechanical Engineering, Indian Institute of Technology, Kharagpur 721 302, India (Receiwd

17 October 1988)

Abstract-The finite element method along with LooChristensen-Wu higher order bending theory is used for evaluating stresses and stress concentration factor due to a circular hole in fibre reinforced composite laminates subjected to transverse loads. The higher order theory used here is shown to give a stress concentration factor of reliable accuracy compared with the exact elasticity solution. The variation of stress concentration factor with plate thickness, hole size and nature of load distribution is presented.

models may be viewed as power series in the thickness direction Z. The level of truncation of series in eqn (1) A significant increase in the use of fibre reinforced is consistent in the sense that the transverse shear composite laminates in aerospace and mechanical strains due to in-plane displacements U and Y are of engineering structures necessitates the rational analythe same order in z as those determined by the sis of such structures. The analysis of these structures transverse displacement W. The accuracy of the has been extensively studied by means of classical and theory is established in [I] through its application to shear deformation bending theories. These theories the problems of a bidirectional and an angle ply fail to yield results accurate enough in situations laminate subjected to sinusoidal surface loading, and where (i) the laminates are thick, (ii) the stress compared with the results obtained by elasticity gradients are high around discontinuities and contheory. The chosen displacement model characterizes centrated forces, (iii) the ratio of transverse shear both in-plane deformation and out of plane deformodulus to inplane modulus is low, (iv) the loading mation together, the feature that exists in a general gradients are very steep and (v) the material aspect laminate. The in-plane stresses o,_ o,., tO. and disratio (longitudinal modulus to transverse modulus) is placements U, I’, W predicted by higher order theory significant. These situations are usually found [eqn (l)] are in good agreement with the elasticity in composite laminates because of their inherent solution [l]. Out of plane stresses 6,. t,:, tJI, however, construction features. The analysis of laminates for differ from the elasticity solution but the magnitude results of acceptable accuracy by means of elasticity of these stresses is small compared to the former. One theory is impossible in practice. Therefore, a great can try to modify eqn (1) to improve the accuracy of deal of effort has been spent in recent years on the uZ, rII and 7jr. but the extra effort needed costs more development of theories which will serve the purpose for the purpose of evaluating in-plane stresses a,, av and at the same time be capable of being implemented and 7,).. in practice. Lo et al. [l] developed a higher order The laminated composite plates used in the conbending theory to predict the realistic flexural behavstruction of aircraft, automobiles, ships and chemical iour of composite laminated plates based on the vessels are sometimes provided with circular holes to displacement model meet functional or design requirements. In the course of their service, laminates get fractured and develop U(x,~,=)=u,(x,y)+2l.~(x,y)+z*m,(x,y) cracks. These cracks and holes act as stress raisers. Analysis of stresses around these geometric dis+ Z3%(X,Y) continuities for laminates subjected to in-plane loads has been done extensively but the same problem under transverse bending loads has received the attention of few investigators because of its analytical W(x,y,z)=M’O(X,y)+ZI:(X,)‘)+Z2~,(X,y) (1) complexity. From the point of view of fracture mechanics, three-dimensional finite element analysis in which U, V and W are total displacements along is applied by Alwar and Ramchandran [2] for studythe x, J’ and z axes respectively. I+,, uO, M’~,I,, l,, I,, ing the stress intensity factor in cracked thin plates m,, mrr m,, n,yand ny are the generalized displaceunder bending. They also studied the influence of ments in Cartesian coordinates. The displacement crack closure on stress intensity factor for plates subjected to bending [3]. Ramachandran and Alwar [4] also concluded that for a three-layered tPresent address: Central Engineering Research Institute, Durgapur 713 209, India. symmetric laminate with the central layer having a INTRODUCTION

CAS

13&B

929

930

T.

K. PAUL and K. M. RAO

through crack, the value of the stress intensity factor decreases with the increase of modulus ratio and thickness ratio of outer layers to inner layer. Analysis of stresses around a circular hole in an infinite laminated plate consisting of two bonded isotropic layers and subjected to axisymmetric bending has been done by Delale et al. [5]. Fraser.[6] used where L indicates the number of layers. By substivon Karman nonlinear bending theory for bending tuting eqns (1) and (2) in (3) and carrying on the analysis of a highly stretched plate containing an integration over the plate thickness, it is possible to eccentrically plate-reinforced circular hole and subshow that jected to a large biaxial stress system at infinity. (M3 =Plicj (4) Reissner f7] evaluated stress concentration factor in the case of pure bending of an infinite plate with a where (M} is the vector of stress resultants: circular hole using his shear defo~ation theory 181. He has shown that his results are in dose agreement (MI = IN,~,.~:N,.MJf+U?, with those of Alblas [9] obtained by elasticity theory. Lekhnitski [IO] studied a homogenous anisotropic P,P,.P.~,.~.~~~~.~.,,Q,Q,R?.R,S,.S,l’. (5a) plate with a circular hole subjected to pure bending moment along the edges. (c} is the vector of generalized strains: It may be stated that the literature available on the analysis of stresses around circular holes in a finite fe I = ]~0..~ oo., 1, uo,, -t L’~.,I,,, I,,, 2m: I,., fibre reinforced laminated plate under transverse loads is limited. Hence this needs attention for the +!s.r~r,,.,m,,,.flt,,, +M,...,-n,,‘~,.,,.n,., proper design of laminates. The present paper aims to evaluate the stress concentration factor due to circu+n?,, 1, + W”.,f, + l.(‘“..y 2nr, lar holes in cross ply finite square plates subjected to sinusoidal and uniformly distributed loads. The finite +l_ 2m, + l,,x 3n, + rn-_.?3n, + m,,,17 15b) element method along with the higher order plate and (D] is the rigidity matrix given in the Appendix. theory of eqn (1) is adopted for the analysis. The displacement components of eqn (I) can he put in matrix form FORMULATION Each layer of the laminate has fibres oriented in its plane parallel to the XY plane. As such the layers behave transversely isotropic, for which the constitutive equations are

(6a) where 1 0 0 El=

0 I r0 0

0 1

z

0

0 ?

0 0

2 0 0 ;

0 0

0

0

z3

0

2 0

0 ?

0 0

2. 0I (6b)

{ 6 ) is the generalized displacement and

{S} = [~~~~~~ol.~l,!.m~m~ml~.n,.ll and (A} is the displacement

where u,.u,,u~.r,,.,~,,~ and t,, are the stresses and f-,, $3 e:, yXr.yJ; and y._ are the engineering strains; and Q, are the components of anisotropic stiffness matrix referring to the geometric X, y and z axes. The stress resultants [NX. N,, . . . , S, and S,] corresponding to the higher order bending theory of eqn (1) are defined as

vector (6~)

vector

(A) = [U Y W]?

(6d)

Let the piate be subjected to surface tractions X, Y and Z on top surface 2 = -h, along the X. f and 2 axes, and let {q} = surface traction

vector = {X Y Z]‘.

(7)

The strain energy T, is given by T:=f

(6.4, + p,q + C.-C:+ ‘c.,,yvi + s,::r,,+

t,,;,,,)dV

(8)

Stress analysis around circular holes where V is the volume of the plate. Replacing the strains in eqn (8) by strain-displacement relations and integrating over the thickness and using the definitions of eqns (3), (4) and (5) the final form of strain energy T, is given by Ts=; SI

A

{W

1~1dA

(9a)

is the shape function matrix. Each component (i=1,2,... ,8) in [N] is a diagonal matrix [NJ =

{AK= -,z, {q} dA.

Pb)

The total potential energy T of the plate is the sum of eqns (9a) and (9b):

SI

{M}T{~} dA -

A

Substitution 7,=; ss

ss A

n,r

N,

(17)

where I is an 11 x 11 identity matrix. After substitution of eqn (15) into eqn (Sb) the vector {.G)of generalized strains becomes

where A is the area of the plate. The change in potential energy r, of external force vector {4) is Tp= -

931

ff

1= VW)

(18)

where [B] is the strain-displacement matrix of size (20 x 88) which is given in the Appendix. Using eqns @a), (I 5) and (18) in eqn (1 I), the total potential energy (‘I’<) of a typical element e of the plate turns out to be

IAX= -h, 191 dA. (10)

of eqn (4) in eqn (10) yields ,4{~)~Pl{t.) dA -

ss

@h),T_-h, (41 dA. A

Minimization of the potential energy T’ with respect to the unknown displacement 6” yields the stiffness equation of the element

(II)

The region of the plate is discretized by employing an eight noded quadratic quadrilateral isoparametric element. Let the shape functions which are established in [l l] be ni (i = 1,2, . . ,8) associated with the eight nodes. As dictated by higher order bending theory in eqn (1) the displacement vector, 6, at any point in an eiement has 11 generalized displacement components (u,,. v,. w,,, i,, I,., I:, m,, m., ml, n, and n,.), i.e. @t = [u~vowol,l,l,m,m,m,n.n,lT.

(12)

The generalized displacement vector S, of the ith node is {SiI =

[UgiUgiwoi!,if~il;im,imyimtin,in,il'

(13)

where +i, ZJ@(, . . . , nyi are values of generalized displacements at the ith node. The nodal displacement vector 6” of a typical element c is given by is}*= [S,6,6,6,6,6,6,Sg]?

(14)

In the case of an isoparametric element, the displaeement vector, 6, may be expressed by 1s 1 = where

hwv

where [Kq = element stiffness matrix:

and {F’) = element nodal force vector:

(Fe) =

Wl’[Zl:= -A,tq 1d/f.

(22)

By summing the stiffness equations of all the elements, the global stiffness equation results [K’]{S’) = {F’f

(23)

where [K] = global stiffness matrix = Z,L,, [K’], (6’) = global nodal displacement vector = Z:,L=,{A’) and {F’} = global nodal force vector = C,L=,{Fe}. L is the number of elements, in the above summations. The unknown global nodal displacement vector 6’ can be solved from eqn (23). Once the nodal displacements are known, the next stage is to evaluate the stresses a,, cry,5._ at P (the point of interest) in a typical element e. These stresses are given by

(15)

= ~Q’~PIZ’~PW~PWI

(24)

932

T. K. PAULand K. M. RAO

where [Q ‘IF =

.Q,, Qu Qu Q22 .QM Q24

Q,4

Q23

Q24

Q34

Q44I P

r1

0

0

0

iy

0

I

1 0

0

0

z

001000z0000000 0 0 0 1 0

0

0

WI, =

Q,3

Table 1 compares the result for three-layered, symmetric cross-ply (Oa/900~~}, simply supported square laminated plate obtained by three-dimensional elas-

0

0 z2 0 0 z*

0

z3 0

01

0

0

z?

0

0

z*

0

0

23 i p

0 p0

z

0

and [I& is the value of the strain-displacement matrix of eqn (A4) evaluated at the point P. NUMERICS

RESULTS

Numerical results are presented for a simply supported square, three-layered, symmetric cross-ply laminate with a circular hole of diameter D at the centre. The laminate sequence is (0°/90”/00), and the total thickness (h) is distributed among the three layers in the ratio (1:2:1) unless otherwise specified. The laminate properties referring to principal material directions are assumed to be .E_,= 25 x lo6 psi (174.6GPa), ET= 1 x 1O’psi (7GPa), G,,=OS x lo6 psi (3.5 GPa), Grr = 0.2 x IO6psi (1.4 GPa) and uLT= urr= 0.25, where subscript L refers to fibre direction and T transverse to fibre direction, E elastic modulus, G shear modulus and a Poisson’s ratio. The loadings considered are (i) uniform edge moment MO on two opposite edges parallel to y axis, (ii) uniformly distributed load of intensity p,, and (iii) sinusoidally varying load of intensity p, p = p. sin nx/a sin ny lb, where a is the plate length along the x axis and b along the y axis. The stress concentration factor variation is studied with the parameters (i) angle ((3) with respect to the x axis around the hole circumference, (ii) ratio S of in-plane dimension (a) to total plate thickness (h) and (iii) ratio of hole diameter (D) to in-plane dimension (a). Variation of cr, at 8 = 90”, and oYat @= 0” with thickness coordinate (z) is also evaluated. The results are presented in Figs 2-8 and Tables 1 and 2. The stress concentration factor (K) is stated as K=

stress at a point on the circumference of hole stress at the same point without hole ’ DISCUSSION

ticity solution and the present theory. u, is evaluated at the top and bottom surface of the laminate, and oY at the top and bottom surface of the middle layer. Reissner [7] presented the stress concentration factor for pure bending of an infinite isotropic plate with a circular hole with Poisson’s ratio u = 0.25, and the results are compared with the exact values of Alblas [9]. Table 2 compares the present results with those of [7,9]. Bending moment is applied along the edges parallel to the 4’ axis.

h

Fig. 1. Geometry.

Table 2. Stress concentration factor

AND CONCLUSION

The accuracy of the fo~ul~tion and computerization is checked against the existing results of [7, 121.

s W

0,

K (a =

b =

100, D = 2.

pure bending) Dlh

Alblas [9] (elasticity)

Reissner [7]

Present

1.0 2.0 3.0 4.0 5.0

2.052 1,938 1.865 1.841

2.038 1.922 I.875

1.830

1.835

2.150 1.850 1.840 I.835 1.825

Table 1. Comparison with 1121(a = 4 WI

6 =

1.850

20, sinusoidal load) 100 WI

(Elasticity)

Present

(Elasticity)

Present

0.248 x lO-4 +11.52 - 10.94 + 10.61 - 10.66

0.24 x lO-4 +11.1 -11.71 + 9.94 - 9.91

0.0869

0.08685

& 5390

rf:5430

&2710

h2710

Stress analysis around circular holes

933

3-

-II0

IO

20

30

40

50

60

70

60

90

2lll 0.10

8

0.15

0.20

0.25

0.30

0.35

0.40

D/a

Fig. 2. K vs 6 of CFRP (0°/90”/Oq square plate, uniform load (a = 20, b = 20, S = 10 and 100, D/a = 0.1).

Fig. 4. K vs D/a of CFRP (V/90°/Oq square plate, uniform load (a = 20, b = 20, S = 4 and 100).

In accordance with the data of Table 1, the deflection given by the present solution differs from the exact solution by about 3%. and the stresses by about 6% in the case of an extremely thick plate (S = 4); the deviation is negligibly small in the case of a thin plate (S = 100). As stated in Table 2 the stress concentration factor (K) given by the higher order theory differs from that of the exact solution by about S% in the case of a small circular hole (i.e. D/h = l), and the difference is less than 1% in the case of a large hole (i.e. D/h = 5). The comparison concludes that the higher order bending theory adopted here yields results accurate enough for engineering application. The values of the stress concentration factor (K) presented in Figs 2-7 and the stresses in Fig. 8 are obtained for a three-layered, symmetric laminate from the equations derived in formulation by particularizing them to fibre orientations of 0” or 90” in layers. The quantities KX + and KX - refer to the values of stress concentration factor (based on 0,) at

the bottom and top surfaces of the laminate, and KY + and KY - refer to the values of stress concentration factor at the bottom and top surfaces of the middle layer (based on G,,), in the case of an extremely thick plate (S = 2 and 4). For the case of moderately thick and thin plates (S = 10,20,50 and 100) the values of KX + arid KX - are identical and K (based on a,) is represented by KX. Similarly K (based on cr,) is represented by KY as the values of KY + and KY are identical. Figures 2 and 3 illustrate the variation of K around the circumference of the hole, the plate being under uniformly distributed load. It increases as the fibre direction is approached. The stress concentration factor is more than unity within the angular distance of about 25” from the fibre. The values of b, presented in Figs 4-8 are evaluated at 0 = 90” and those of aY at 0 = 0”. The stress raising ability of the hole in general decreases as the hole size increases (Figs 4 and 5) within the range of parameter values considered. However in the case of 6

K

K4

o-

_----

-‘II0

US---__

8

Fig. 3. K vs 0 of CFRP (V/90°/Oq square plate, uniform load (a = 20, b = 20, S = 2, D/a = 0.1).

0.10

0.15

0.20

0.25

0.30

0.35

0.40

D/a

Fig. 5. K vs D/a of CFRP (0”/90”/0”) square plate, sinusoidal load (a = 20, b = 20, S = 2 and 100).

T. K.

934

PAUL and

RAO

K. M.

r

loo

(a)

50 -

KY+ 2

4

I

I

I

I

IO

20

50

100

-100

I

-

-2.5

-

1

1.5

I

-0.5

s

I

I

0.5

1.5

2.5

I

Fig. 6. K vsS of CFRP (V/9O”/Oq square plate, sinusoidal

load (a = 20, b = 20, D/a = 0.1 and 0.4).

IOO

a thin plate (say S = 100) under uniformly distributed load (Fig. 4) the value of K decreases initially, but increases later. The positive and negative values of stress concentration factor visibly differ (Figs 6 and 7) in the case of thick plates (say S c 10). Figures 6 and 7 show the variation of K with the reduction of plate thickness (i.e. as S increases). The values of all KX, except that for a small hole (i.e. D/a = 0.l), increase with S in the case of a thick plate (say S < lo), and decrease with S in the case of thin plates (i.e. S > 10). The stress concentration factor KY decreases over 2 < S < 10, increases over 10 < S < 50, and then decreases in the case of small hole size (D/a = 0.1); whereas for large hole size (D/a = 0.4), KY decreases over 2 < S < 4, increases over 4 < S < 20 and then decreases over 20 < S < 100. As a whole it can also be observed from Figs 6 and 7 that the values of KX are higher than those of KY in the case of thick plates, but the reverse is true in the case of thin plates.

-1001 -2.5

K4

(b)

I

I - I.5

I

-0.5

0.5

I I.5

I 2.5

Z

2 x I04

u

6

r

0

-2

x I04

-4

I 1041 -0.1

I 0

1 0. I

Z Fig. 8. 0 vs z of CFRP (0”/9W/O’) square plate. (a) Sinusoidal load (a = 20, b = 20, S = 4, D/a = 0. I and 0.2). (b) Uniform load (a = 20, b = 20, S = 4, D/a = 0.1 and 0.2). (c) Uniform load (a = 20, b = 20, S = 100. D/a = 0.1 and 0.2).

KY+

22:o

S

Fig. 7.

K vs S

of CFRP (0°/9W/00) square plate, uniform load (a = 20, b = 20, D/a = 0.1 and 0.4).

The influence of load distribution can be studied from Figs 4-7. Uniformly distributed load (Figs 4 and 7) induces higher stress concentration factor than sinusoidally varying load (Figs 5 and 6) (see graph,

Stress analysis around circular holes

S = 100). The reduction in stress concentration factor with hole size is faster in the case of sinusoidally varying load than in the case of uniformly distributed load (Figs 4 and 5). The variation of a, and ar over plate thickness for the two types of load distribution is plotted in Fig. 8 (a)--(c). The stresses in the outer layers vary nonlinearly and those in the middle layer vary linearly in the case of thick plates [Fig. 8 (a) and (b)]. In the case of a thin plate [say S = 100, Fig. 8(c)] in the outer layers too the variation of stress is almost linear. REFERENCES

4. K. N. Ramachandran and R. S. Alwar. Laminated plate with an embedded through crack in the central layer. J. Reinforced Plas. Comp. 4, 314-323 (1985). 5. F. Delale, N. N. Kishore and A. S. D. Wang, Stress

6. 7.

8.

9.

1, K. H. Lo, R, M. Christensen and E. M. Wu, A higherorder theory of plate deformation, part 2: laminated nlates. J. annl. Mech. 44. 669-676 (1977). 2. k. S. Alwai‘and N. Ramachandran,Three dimensional finite element analysis of cracked thin plates in bending. Int. J. Numer. Meth. Engng 19, 293-303 (1983). 3. R. S. Alwar and N. Ramachandran, Influence of crack closure on stress intensity factor for plates subjected to bending. Engng Fruct. Mech. 17, 323-333 (1983).

935

10.

1I, 12.

analysis of a composite plate with a circular hole under axisymmetric bending. J. Camp. Muter. 18(3), 42&431 (1984). W. B. Fraser, Bending of a highly stretched plate containing an eccentrically plate-reinforced circular hole. Int. J. Solids Struct. 11, 501-518 (1975). E. Reissner, On transverse bending of plates including the effect of transverse shear defo~ation. Int. J. Solids Strucr. 11, 569-573 (1975). E. Reissner, The effect of transverse shear deformation on the bending of elastic plates. J. uppl. Mech. 12, A69-All (1945). J. B. Alblas, Theorie van de driedimensionale spanning stoestand in een doorboorde. Dissertation, Delft (1957). S. G. Lekhnitski, Anisotropic Plates (English translation from Russian 2nd Edn, by S. W. Tsai and J. Cheron), pp. 394419. Gordon and Breach, New York (1968). 0. C. Zienkiewin, The Finite Element Method, 3rd Edn. pp. 1555160. McGraw-Hill, New Delhi (1979). N. J. Pagan0 and S. J. Hatfield. Elastic behavior of multilayered bidirectional composites. AIAA Jnl 10, 931-933

For Appendix see over

(1972).

where

Symmetric

,I

0

(1. z, z‘. z’. +, z’, :“)Q,,dz.

4,

0 0

0

0

A,, A60

0 0 0 0

0

Fz4

Fu G,, G, GM

Fzz

Fzd G,z Gz2

F,z

FM G,,

0

0 0

4, 44 Ew Eu FM 0 0

0

0 0 0 0 0 0 0 0

0

0

D,

D,

0

0

Dx E,z 4 41 44 F,z

&,

D,,

Dz2

D,,

4, D,, 4, EIL 6, 4, F,,

D,,

DE

4,

PI = P, B, 4 B, & 4, B, 41

[A,,. 4,. C,,.D,,. 6,. F,,. G‘,,l=

APPENDIX

cs

4, 4,

0 0 0 0

0

0 0 0 0 0 0 0

0

0

B, C,, CM

4

0 0 0 0

0

0 0 0 0 0 0 0

0

0

C,, G, 4, 4, 4,

0 0 0 0

0

0 0 0 0 0 0 0

0

0

0

C, 4, D, 46 46

C,,

0 0 0 0

0

0 0 0 0 0 0 0

0

b43)

642)

Stress analysis around circular holes

~~~~~OOOO~-O

000000000000

00000000000

,;

,2

00000000

0

z-0

0

0

00000

E

0

,:

z: 0

0000

0

0

0

,:o

0

0

0

0

~~000000

0

~~0000000

0

0

ho

0

0

0

0

z-0

0

0

0

0

0

0

0

0

,;

iooooo

~-0000

2

4

%

0

0

0

0

0

0

0

0

0

0

.Y

s 0

=:

0

~0000000000000000

0

~~0000000000000000

hi L

II

s

-

h

0

0

0

0

0

=--.

,:

0000

0

,:o

z-

0

0

931