A high precision triangular laminated anisotropic cylindrical shell finite element

A HIGH PRECISION TRIANGULAR LAMINATED ANISOTROPIC CYLINDRICAL SHELL FINITE ELEMENT H. V. LAKSHM~NARAYANA and S. VISWANATH National Aeronautical Laboratory, Bangalore-560017, India (Receiud 1 February

1976)

Ah&act-The stiffness matrix for a high precision triangukr laminated anisotropic cylindrical shell finite element has been formulated and coded into a composite structural analysis program. The versatility of the element’s fomtuktion enabks its use in the analysis of multilayered composite plate and cylindrical shell type structures taking into account actual lamination parameters. The example applications presented demonstrated that accurate predictions of stresses as well as displacements are obtained with modest number of elements.

iNrRowcTroN Multilayered

composite

materials are increasingly used in the construction of plate and shell type structures for various industrial and aero-space applications, creating a parallel need for structural analysis procedures which reflect the laminated anisotropic characteristics of such materials. Although the basic equations which govern the elastic behaviour of laminated anisotropic plate and shell structures can be derived[l], analytical solutions of these complex equations, as they pertain to design and analysis, are not available. The analysis of composite structures is an area which calls for an extensive use of digital computers and numerical methods like the finite element method. The finite element method has proved to be an extremely powerful tool for the solution of problems involving complex geometries, arbitrary loadings and rather general material properties. Considerable advances have been made in the development and application of the finite element method to plate and shell type structures. Much of the work is reviewed by Gallagher[2,3]. According to these reviews, the analysis of thin shells istone of the more difficult problems that has been attempted with the finite element method. The requirements upon valid minimum potential energy solutions in thin shell finite element analysis are extremely diflicult to satisfy, and the satisfactory formulations are therefore relatively complicated. Further, for tackling problems involving description of arbitrary boundaries triangular and quadrilateral elements are obviously best suited. The most reliable and sophisticated triangular shell element formulations are those due to Dupuis[4], Cowper er al.[S], Argyris and Scharph[6] and Dawe[7]. The formulation due to Cowper et al. (51 has also been converted into a triangular cylindrical shell finite element by Lindberg and Olson[t?]. The salient features of this element are: (i) The element is fully conforming, of the displacement type and of arbitrary triangular shape, (ii) Uses higher order interpolation polynomials and has 36 degrees of freedom, (iii) Transverse displacement function contains a complete quartic polynomial plus some higher degree terms and allows cubic variation of normal slope along each edge. (iv) Tangential displacement functions are complete cubic polynomials to yield a quadratic variation of stresses within the element, (v) Gives accurate predictions of displacements as well as stresses, (vi) Satisfies sullicient conditions to guarantee rapid convergence, (vii) Permits reliable and direct 633

evaluation of stresses at nodes, (viii) The element stiffness matrix is formulated in exact closed form, (ix) Results in a smaller overall probkm size and (x) The versatility of the formulation enables its use in the analysis of Bat plates in either membrane or flexural behaviour and for the coupled membrane-flex& behaviour of cylindrical shells. The formulations of this element as given in Ref. [B] is confined to homogeneous isotropic materials. The present work is concerned with the generalization of this formulation to take into account the laminated anisotropic characteristics of multilayered composite materials. The purpose of this paper is to present the formulation and some useful applications of a high precision triangular laminated anisotropic cylindrical shell finite element. In order to perform composite structural analysis, a computer program incorporating this element has been developed. Distinct membrane and flexural stiflnesses are retained permitting analysis of composite structures subjected to inplane and bending loads. The applications presented encompassing flat plates in plane stress and in flexure and circular cylindrical shells clearly indicate the usefulness of the element in composite structural analysis. s0RMuLATIoN The geometry of an arbitrary triangular cylindrical shell element is shown in Fig. 1. The global cylindrical coordinate system is X, Y and 6, rt are taken as local coordinates for the element. R is the radius of curvature of the reference surface of the shell and I is its wall thickness. The dimesnions a, 6, c of the triangle I, 2, 3 and the rotation angle 8 are easily derived in terms of the global coordinates of the vertices[9]. The strain energy of a thin cylindrical shell in global coordinates is given by

“=;

II WITkol+WITkl~~dy

(1)

where iNIT MT ]e?’ lul’

= (NX. NY, Nxu) are the membrane forces = (Mx, MY,MXy ) are the bending stress resultants = (Ex”,l ,“, y’x~) are the reference surface strains and = (XX,XV,xm) are the bending curvatures.

For an arbitrarily laminated anisotropic shell, the con-

634

H. V.

LAKSHMINARAYANA

and S. VISWANA~H

when radially outward. Note that subscripts in eqns (4) and (5) denote differentiation. Substituting for [co] and h] from eqns (4) and (5), it is possible to express the strain energy given by eqn (3) in matrix form as

WI= IL1[dll dx dy

(6)

where

WIT= (Ux, CJ,, Vx. Vw W, Wxx, Wxv, Ww 1 and the strain energy matrix [I_.]is given in Table I. If 0, I7 denote reference surface displacements in the local coordinate system (6, n), the transformation from local to global coordinates is

VI = I&l 161

(7)

where Wall details

Fig. 1. Cylindricalshell triangular finite element -Geometry coordinatesystems.

and

and the rotation matrix [R,] is given in Table 2. Substituting eqn (7) in eqn (6). the strain energy of the element is given by

stitutive relations are expressible in the form [ 11

(2) where where [A], [B], [D] are the membrane, interaction and bending stiffness matrices of the laminate. Substituting eqns (2) in (l), we get

Ik”l’ Mlk”l + M’[~Ik”l +~~“l’~~l~~l~~~IT~~I~~l~~~dy.

[iI = VU7 WIRI with the area integration to be carried over the triangle 123. The assumed displacement functions for the element are ]51

(3)

Using Novozhilov’s thin shell theory, the straindisplacement relations are [101

[Xl= [

-wxx 1

-(WIT - Vwn) -awx, - Vx/R)

(9) (4)

(5)

where m,, a,, p,, r, and s, are just integers(51. The element stillness matrix is obtained from a calculation of its strain energy. llte displacement functions of eqns (9) are substituted into eqn (8) and the integration carried out to yield

where U, V, W are the displacement components in the global coordinate system (X, Y, 2). W being positive

u = ; [al’[kl[al

(10)

Table 1. Strain energy matrix IL.1

-m,

-42

- 2&,

-&3

-2B,,-9 _2B -20, ” R - 2Bzl/R Symmetric

2%

-2 2D

-B

a -

R

&z/R 012 243 D22

A high precisiontriangularlaminatedanisotropiccylindrical shell finite element Table 2. Rotation matrix [R.,] cc

-SC

E 2 ss SC 00 1 where

-SC -SS :: 0

SS 0 -SC 0 -SC00 cc 0 010

0 0 0 ;; ss

0 0 0

0 0 0

:

8

CC?& 2sc

-;; cc _

CC = cos*0 SS = sin’ ti SC = sin 0 cos 0

where [a] is a 40-column vector of the polynomial coefficients ai. The entries of the stiffness matrix [k] are determined in closed form and are given in the Appendix. The generalized displacements in global coordinates chosen for the element are T [WI = tu,, UX,, uu,, VI, VX,, VYl, Wl, WX,, WY,, WXX,, wJrv19WYYI,u*. . .9 u,,.. ., WY,,).

635

Boron-epoxy composite laminate with a circular hole under axial tension In Fig. 2 is shown a boron-epoxy composite panel consisting of nine plies of [Ol+ 45/O/90], laid up with a central circular hole under uniaxial tension. Symmetry of the loading, geometry and material properties make the analysis of only one quarter of the panel su5cient. The finite element idealization of the quarter panel illustrated in Fig. 3 uses 51 elements with 37 nodes resulting in 186 degrees of freedom. Typical stress distributions calculated at the nodes are shown in Figs. 4 and 5 along with the results of anisotropic elasticity solution to the problem of an infinite plate with a circular opening[ll]. The present results appear to have converged quite rapidly. The predicted tensile and compressive stress and strain concentration factors given in Table 3 together with results of anisotropic elasticity solution for the plate clearly demonstrate that the element gives excellent accuracy. An identical problem has been analysed in Ref. [12] using linear-strain triangular elements (twelve degrees of freedom). One hundred and forty one elements with 320 nodes (598 degrees of freedom) were used to describe

(11)

Further steps to be followed in the deviation of the 36 x 36 stiffness matrix relating the generalized displacements [WI and the corresponding generalized forces, and consistent load vector corresponding to uniform pressure load are very similar to those for the shallow shell element [S].

The stiffness matrix and consistent load vector for the high precision triangular laminated anisotropic cylindrical shell finite element has been coded into a composite structural analysis program. Distinct membrane and flexural stifftresses are retained permitting analysis of composite structures subjected to inplane and bending loads. The resulting force-displacement equations are solved using Choleski square root decomposition method for fixed bandwidth problems. Stress resultants and individual layer stresses are directly calculated at the nodes using the generalized displacements, constitutive relations and appropriate transformation equations. With basic lamina elastic constants and laminate details as input, the program evaluates the laminate stiffnesses for each element. It is therefore possible to analyse laminates of varying thickness taking into account the actual lamination parameters. The computer output consists of generalized displacements, stress resultants, membrane and bending strains, individual layer stresses with reference to both laminate and lamina reference axes at each node. With this information it is possible to make strength prediction of practical composite structures using appropriate failure criteria[l].

1f,., j

l~;T$ 1

Fig. 2. Boron-epoxylaminate with a circular hole under axial tension.

aR

APPLICATIONS

The versatility of the element formulation enables its use in the analysis of laminated anisotropic flat plates in either membrane or flexural behaviour and for coupled membrane-bending behaviour of circular cylindrical shells. The following example applications illustrate the usefulness of the element in composite structural andysis. AI1 calculations are carried out on an IBM 370/l% machine using double precision arithmetic.

Fig. 3. 5 x 5 Finate element grid,

H. V.

636

LAKSHMINARAYANA ad

Fig, 4. Axial stress distribution along X and Y axes. -,

Fig. 5. Tangential stress distribution along hole boundary. -,

the quarter panel. The predictions are also given in Table 3 for comparison. The fact that the present formulation yields better accuracy with about one-third the problem size demonstrate the superiority of the element. The problem has also been analysed experimentally using straingages, pbotoelastic coatings and moire techniques in Ref. [12]. Average predictions are given in Table 3. The close agreement between the present results and experimental data readily verifies the adequacy of the assumed laminated anisotropic model in the analysis as a first approximation to the actual multilayered composite material behaviour. Clamped to uniform

glass-epoxy pressure

load.

composite Fiiure

square plate subjected 6 shows a glass-epoxy

composite square plate with clamped edges subjected to uniform pressure load. Symmetry allows the analysis to be limited to one quarter of the plate. Arrangement of

S. VWANATH

Analytical solution; A, finite element results.

Analytical solution; A, finite element results.

the elements in a quadrant of the plate are also shown in Fig. 6. Predicted deflection distributions along centre lines of the plate are shown in Fii. 7 along with the results from an approximate analytical solution given in Ref. [13]. This analytical solution assumes identical deflection profiles along the centrelines of the plate. Whereas, as it should be, the present analysis predicts different detlection profiles parallel and perpendicular to the fibre direction. This prediction is in qualitative agreement with experimental observations on unidirectional composite plates in kxure using holographic techniques(l41. Bending moment distributions along the centrelines of the plate shown in Fig. 8 are fairly smooth and can be considered to have converged as well. Typical results given in Table 4 can be considered more accurate than the solution available in[13]. However, an assessment of actual accuracy of the predicted stresses

Table 3. Boroncpoxy laminate with a circular hole under axial tensiot+Compatison of results No. of UllkllOWllS

Present results Analytical solution. Ref. [12] Finite element results. Rcf [9] Experimental results. Ref. 191

Stress concentration factors at O=o”

8=!W

strain iXlncentration factors at 0-o”

O=!XP

186

- a.6957

3.4982

- 1.4725

3.449

-

- 0.6829

3.4s

- 1.46

3.45

598

-

3.60

- 1.58

3.60

-

-

3.34

- 1.5

3.34

A It&

precision tria@ar

laminated anisotropic cyihdrid

she4 finite element

637

Pressure P

X

Fig. 6. Finite element grid for clamped glass-epoxy composite square plate.

24

“Q

2.0

I

6% fQ d s

e E i

1.0

; %

0 Distance along ccntre line x/u or@

Fig. 7. Deflections of clamped&ss-cpoxy composite square plate. -, Analytical solutions; A, along X-axis; 0, along Y-axis.

and deflections is not possible due to inherent approximations in the available solution in Ref. [13]. (The s&ion in [13] when specialized to isotropic case shows an error of about 6.3% for maximum defkction and unacceptable accuracy for maximum moments.) An oflhotmpic cylindrical shell with a circular hole subjected to axial tension. The problem illustrated in Fii 9 is analysed to evaluate the coupled membrane-bending behaviour of the cylindrical shell form of the element. Taking advantage of symmetry, only one quarter of the shell needs to be analysed. Choosing the 5 x 5 finite element grid illustrated in Fig. 3.51 elements with 37 nodes (367 degrees of freedom) are used to model the quartet of the shell. The adequacy of the chosen grid to solve the present probiem was checked by solving ~II otherwise identical isotropic shell probkm. Typical results are given in Tabk 5 along with continuum sdutions[lS] for comparison to demonstrate the accuracy achieved with the rather coarse mesh used. Predicted stress distributions for both isotropic and orthotropic shells are shown in Figs. IOand II along with

0.1

0.0 Disfmc*

along

02

0.3

centre line X/o

0.4

0.5

or Y/O

Fig. 8. Moment distriitions

along centre hes of clamped 8h.w epoxy composite square plate.

the results of anisotropic elasticity solution to the problem of an infinite plate with a circular hole[ 111. The predicted stresss distriins are fairly smooth and’ appear to have converged as well. Near the hole the coupled membrane-bending behaviour appears to be

638

H. V. LAKSHMINARAYANA and S. VISWANATH Table 4. Clamped glass-epoxy composite square plate subjected to uniform pressure loadcomparison of results Maximum deflection Present results Approximate analytical solution Ref. [lo]

Bending moments atX=a/2, Y=O

WWpa4

M&a2

0.00216

- 0.0744

0.002313

-0.06517

Bending moments at X = 0, Y = a/2 Wxlpa2

M&a*

- 0.00495

- 0.00732

- 0.02938

- 0.004333

- 0.00433

- 0.01732

Myha

Table 5. Cylindrical shell with a circular hole under axial tension-Comparison of results Membrane stress ratio GJem/uOo) at#=lo” Orthotropic shell Present results Orthotropic plate Solution Ref. [12] Isotropic shell Present results Isotropic shell results from Ref. [I 11

Bending stresst ratio (Q%%o)

ate=90

at fI=O

at 8=90”

-0.301821

0.4391

- 0.503648

4.9389

- 0.4813

4.8963

- 1.26018

3.70422

- 0.67575

0.5203

- 1.25

3.6

- 0.80

0.45

t

0.0

0.0

tAt inner surface of the shell.

0 R

/

-.

.-.

y-

MATERIAL GEOMETRIC

-i_R

_ -100

PARAMETERS

CONSTANTS

ISOTROPIC (Alumlnuml E I,

ORTHOTROPIC Gloss-Epoxy COIVIPO~I~~

100X IO6 psi

6.l3XlJ

lo-oxlo6pl~

I 42X106~.n

3.64xdprl

0.53dp*t

0.3

0.27

0.3

0.0625

pri

Fig. 9. Cylindrical shell with a circular hole under axial tension. characterized adequately by the small number of elements used. Predicted membrane and bending stress concentration factors are presented in Table 5. The main difIiculty in assessing accuracy of these results arises from the fact that authors coud not find any other result for the orthotropic shell for comparison. Since the isotropic shell is a particular case of the orthotropic shell treated here it might be expected that the same order of accuracy holds good for both. coNcL.usIoNs

The stiffness matrix for a high precision triangular laminated anisotropic cylidrical shell finite elements has

been formulated, and coded into a composite structural analysis program. The versatility of the element’s formulation enables its use in the analysis of multilayered composite plate and cylindrical shell structures taking into account actual lamination parameters. The results obtained with the use of the subject finite element demonstrated good engineering accuracy for stress as well as displacements in the example applications with a modest number of elements. The ultimate objective of any analysis is as a design tool. This work implies the iterative use of the analysis procedures in a design exercise. The present finite element formulation is ideally suited for computer aided design of composite structures.

639

A high przcision triangular laminated anisotropic cylindrical shell finite element

50

4s

Q

4C

,-

3:

,-

Membronc stress Orthotropic) Membfone stress (Isotropk)

3c )-

b’ s_ e

2: i-

iE f

2( ,-

5-

,(At

inner surface)

Bending stress(lsotrqxc)

5-

@ending stress@rthotropk)

3-

__-------_ _tIl!!3 Distance, Y/a

Fig. 10. Axial stress distriiution along Y-axis. B

3

K2

z-

r

0

z = --__-_- ___---- _

c

-

00

x.1 -

Advances in Compvrarional Methods in SlrWUral Mechanics and &sign. pp. 641-678 (Edited by J. T. Oden n al.) UAH Press (1972). 4. G. Dupuis, Application of Ritz method to thin elastic shell analysis. Trans ASME, I. Appl. bkh. 38. !k. E, (4) 787 (1971). 5. G. R. Cow-per,G. M. Lindberg and M. D. Olson, A shallow shell finite element of &ngular shape. Inl. I. Solids Stmcl. 6, 1133-1156(1970). 6. 1. H. kgyh and D. W. Scbarpf. The SHEBA family Of shell elements for the matrix displacement method. AeroMKlical 1.. 72(697),873-883 (Oct. 1968). 7. D. J. Dawe, Higher-order triangular finite ekment for shell analysis. Jnl. 1. Solids SINcl. 11(IO), 1097-l 110(1975). 8. G. M. Lindberg and hf. D. Olson, A high precision triangular cylindrical shell finite element. AIM 1. 9.531-532 (1971). 9. G. R Cowper, E. Kosko, G. M. Lindberg and hf. D. Olson. A high precision triangular plate bending element. Aeronautical Report LR-514,National Research Council of Canada, Ottawa (Dec. 1968). IO. V. V. Novozhilov, The Theory of Thin Shells. 2nd Edn. Nwrdhoff, Amsterdam ( 1964). 11.S.G. L&nit&ii, Theory of EIaslicity of Anisotropic Mastic Body (Translated by P. Fern, Edited by J. I. Brandstatter) Holdenday, San Fransisco (1963). 12. I. M. Daniel, R. E. Rowlands and J. B. Whiteside, Deformation and failure of Boroncpoxy plate with circular hole. Analyses of test methods for high modulus fibers and composites, STP 521, ASTM pp. 143-164(1973). 13. J. E. Ashton and J. M. Whitney, Theory of laminated plates. In Prognss in Material scirnce scrics, Vol. IV, pp. 51-H. Techomic (1970). 14. A. Subramanian et al., Holographic analysis of some eaeineering components. NAL Technical Memorandum, TM-PR-ST201/127-76. IS. P. Van Dyke, Stresses about a circular hole in a cylindrical shell. AMA J., 3(S), 1733-1742(l%S). Stiffness Matrix [k] ki,=i,,mim~(mi+m,-2.n,+n,) j 3 i + i12(qn, + mjn,)F(m, + m, - 1, n, t n, - 1) t &kP#k

Pj- 2,nit qj)

t

t i,,mnrF(m, tPj -

1, nl+qj - 1) t L,,mS(m, t ri - 1. n, t sj)Gj t i,,m,rj(ri - I)F(m, t r, - 3. ni t s,) t 5 - 2. ni t

t i,,m,r&$ni t ilcqs,(sj-

-I c

+ L.,nflj$n, t n+, n, t n, - 2) ti,,n,pjF:(mitpi-l,n,tqi-1)

P “8

t LnAf(mi ! g E

s,- 1)

1)&n, + r, - 1, ni t s,- 2)

-2s

t

Pfi4 t qj- 2)

t l&,F(m,t rj,nit sj - l)G, t i,nirj(c - I)F(m, t r, - 2, n, t s,- I) t f$,r#(m, + rj - I, n, t s,- 2) t L,n,sj(s, - l)F(m, t r,,n,t s,- 3)

‘(

O! oc

+ &PP(P,

+ Pj - 2,4( + e)

+~w@Aj+pP~)F(P~+pj-l,q,+qj-1)

Fig. Il. Tangential stress distribution along the hole bou&ry.

+ i,,pP@,

f ‘i - 1. qi + s,,Gj

+i)6P,~(~-I)F(Pi+~-3rq,+s,) f i37PJpf(Pi

1. L. R. Calcote, The AMIYS~S of Lmninad mm. Van Nostmd. New York (1%9).

Composite Smc-

2. R H. Gallagher, Analysis of plate and sbcll structures. Proc. Conf. Appl. Fiite

Elemenr Methoa? in Civil Engng. Van-

derbilt University, U.S.A. pp. 155-206(Nov. 1969). 3. R. H. Gallagher, Application of finite ckment analysis. In

f rI - 2.Q

+ Sj - 1)

+i~~~j-~F(Pi+~-l,q~+S,-2) + &#?fqJQi + Pfi + Ll,eF@i

4r + q - 2)

+ rjv a + +.,)G,

+i,q,~(ri-,)F@,+~-2,q,+s,-l) +L,q,r~~~,+r,-l,q,+s,-2)

H. V. L~SHMIN~YANA and S. VISWANA~

640

+ i&S&-jF(Pi + ‘j3 4i

+

sf- 3).

. t L$jF(Q + rj* Si + Sj) G($j +iuq(rj-~)F(r;,+~-2,sjtsi)Gi +i,r~~(a+ri-l,sitsi-l)Gi

i,si(~i - i)F(c + rp pi f si-2)Gj +&Qrj(S- I)($- i)F(g+ rj-4*Si+sji ti6,~rj(~-l)S~(Strj-3,Si+Sj-l) t i,[r,sj(ri - I)@ - I)+ risi(ri - l)(Si - I)1

t isssis,(si - I)@ - l)F(ri + ri, Sj +Sj -4) = kii where F(m, n) = c”+‘[gm+l-(-

by+‘]

+

XF(rf+~-2,SitSj-22)

i,,resi8+q6 t 5 - 2,Sit Sj- 2) t i,~risisi(s, - l)F(ri + ‘;-- 1, Si+ Sj- 3)

t

nt,ni=O

for

pi,qi=O for

i>lO

ll>i>20

r;:,Si=O for i<21 Gi=O

for i<21,

Gi= i for ia21.

m!n!

(fntnt2)!