Free vibration response of a three-node shear flexible finite element for moderately deep shells

Finite Elements in Analysis and Design 39 (2003) 257 – 281 www.elsevier.com/locate/"nel

Free vibration response of a three-node shear (exible "nite element for moderately deep shells Humayun R.H. Kabir1 ∗ Department of Civil Engineering, Kuwait University, 13060 Safat, Kuwait Received 28 August 2001; accepted 7 January 2002

Abstract A free vibration response of a shear-locking free isoparametric shear (exible three-node triangular "nite element suitable for moderately deep shell con"gurations is presented. The moderately deep shell con"guration based on the Sanders’ shell theory is incorporated into the element formulation. A shear correction term is introduced in transverse shear strain components to alleviate the shear-locking phenomenon. The element is developed with a full integration scheme. Natural frequencies and mode shapes are obtained and compared for doubly curved shells with the available analytical and "nite element solutions. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Three-node element; Shells; Shear locking; Vibration

1. Introduction In industries such as the space shuttle program of NASA, international space center, etc., where reusable hardware require certi"cation for its safe structural use, the threshold of such certi"cation is based on the margin of safety with respect to yield or ultimate strength of materials. The structural safety analysis is performed in most of the cases using linear=nonlinear "nite element models. An improvement or accuracy in the "nite element methodology would save very expensive reusable hardware, and provide an e:cient design to new hardware as well. In modeling shell-like structures, triangular elements are required along with quadrilateral elements to adopt the proper geometry of complex bodies. Most popular elements in industries are three- and four-node shell elements. Threeor four-node shear (exible isoparametric elements show shear-locking phenomenon in thin shell situations when elements are formulated in a conventional way with a full integration scheme. These ∗

Tel.: +965-481-7240. E-mail address: naim [email protected] (H.R.H. Kabir). 1 Present address: USA Solid Rocket Booster Element, NASA Kennedy Space Center, FL, USA.

0168-874X/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 2 ) 0 0 0 6 7 - 7

258

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 1. A typical doubly curved panel.

locking eGects have been circumvented using various schemes as can be found in the works of Dvorkin and Bathe [1], Bathe and Dvorkin [2], Zienkiewicz and Lefebvre [3], Papadopoulus and Taylor [4], Verwood and Kok [5], Zienkiewicz et al. [6], Tessler [7], Kabir [8,9], Onate et al. [10], Sydenstricker and Landau [11], etc. among others. However, to take advantage of shell curvature eGects in the element formulation it is necessary to include shell theories in C 0 element formulation. The main objective of the present eGort is to develop and study frequency response and mode shapes of three-node element for doubly curved shells using moderately deep shell theory of Sanders [12]. The second objective of the eGort is to compare the results with the available analytical and "nite element solutions. Fig. 1 shows the middle surface of a doubly curved panel of rectangular planform of a total thickness of h. x1 and x2 axes represent the directions of the lines of curvature of the middle surface, while the x3 -axis is a straight line perpendicular to the middle surface. a and b are spans along the curved length of axes x1 and x2 , respectively. R1 and R2 are radii of curvature along axes x1 and x2 , respectively. The equations of motion based on Sanders’ shell theory are as follows [12–14]: N1; 1 + N6; 2 + cM6; 2 +

Q1 = C1

; R1

(1a)

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

259

Fig. 2. A triangular element with local and natural coordinate system.

N6; 1 + N2; 2 − cM6; 1 + Q1; 1 + Q2; 2 −

Q2 = C2

; R2

(1b)

N1 N2 − = C3

; R1 R2

(1c)

M1; 1 + M6; 2 − Q1 = C4

;

(1d)

M6; 1 + M2; 2 − Q2 = C5

;

(1e)

where




2P2 Ci = P 1 + Ri








ui; tt

P3 + P2 + Ri

 u(i+3); tt

C3

= P1 u3; tt ; 
2P3

ui; tt + (P3 )u(i+3); tt Ci+3 = P2 + Ri

(i = 1; 2);

(2a) (2b)

(i = 1; 2)

(2c)

in which surface parallel and rotational inertias are included. De"nitions of Pi (i = 1; 2; 3) are as follows:  h=2 (P1 ; P2 ; P3 ) = (1; x3 ; x32 ) d x3 ; (2d) −h=2

where  is the mass density of the material. N1 ; N2 ; N6 are surface parallel stress resultants, while M1 ; M2 ; M6 are moment resultants, and Q1 and Q2 are the transverse shear stress resultants, all per unit length. c is de"ned as 12 (1=R1 − 1=R2 ). The stress resultants related to the strains through the constitutive law are as follows: {N } = {N1 ; N2 ; N6 }T = [A]{}; ˜

(3a)

{M } = {M1 ; M2 ; M6 }T = [D]{};

(3b)

P {Q} = {Q1 ; Q2 }T = [G]{};

(3c)

260

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 3. Finite element representation of a doubly curved panel with triangular elements (32 × 32 elements).

Fig. 4. Variations of normalized !1∗ with respect to radius/span for a panel with a=h = 10.

P are the matrices of the membrane, bending and transverse shear stiGnesses, where [A], [D], and [G] respectively. For an isotropic shell with constant thickness of h, the de"nitions of constitutive matrices are as follows:   1  0 Eh   1 0 [A] = (4a)  ; (1 − 2 ) 1 (1 − ) Sym 2 [D] =

1 12

h2 [A];

P = k 2 hG [G]



1 0

(4b)  0 ; 1

(4c)

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

261

Fig. 5. Variations of normalized !1∗ with respect to radius/span for a panel with a=h = 100.

where E and G are Young’s and transverse shear moduli, respectively.  denotes Poisson’s ratio. k 2 is taken as 56 [15]. {}; ˜ {}, and {}, ˜ respectively, de"ned in terms of the displacements and their derivatives are as follows:  u3      u1; 1 +    R1     ˜1      u 3 ˜2 = u2; 2 + ; (5a)        R2    ˜6     u2; 1 + u1; 2     u4; 1        1     u5; 2 2 = ; (5b)           3 (u5; 1 + u4; 2 ) − c(u2; 1 + u1; 2 )      u3; 2 + u5 − u2   4 R2  = : (5c) {} ˜ = u1    5  u3; 1 + u4 − R1 2. Finite element formulation The "nite element formulation based on strain–displacement relationships as mentioned in Eq. (5) with isoparametric interpolation functions leads to a shear-locking phenomenon using a full integration scheme. In order to alleviate the locking eGects, a correction in transverse shear strain

262

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 6. Variations of normalized !1∗ with respect to radius/span for a panel with a=h = 1000.

Fig. 7. Variations of normalized !2∗ with respect to radius/span for a panel with a=h = 10.

components is considered [9,8]. The corrected strain–displacement relationships for the case of transverse shear strains are now expressed in the following form in terms of element local orthogonal coordinate system (x1
; x2
; x3
):

∗ 4 = 4 +

(x2
)

(6a)

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

263

Fig. 8. Variations of normalized !2∗ with respect to radius/span for a panel with a=h = 100.

Fig. 9. Variations of normalized !2∗ with respect to radius/span for a panel with a=h = 1000.

or


∗ 4 = u3; 2 + u5 −

∗ 5 = 5 +

(x1
)

u2
+ Ru3; x2
; R
2

(6b) (7a)

264

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 10. Variations of normalized frequency with respect to mode shapes for a panel with a=h = 100, and R1 =a = R2 =b = 40.

or


∗ 5 = u3; 1 + u4 −

u1
+ Ru3; x1
: R
1

(7b)

Radius of curvature, in element local Cartesian coordinate system, is calculated from the surface geometry of the shell at pre-processing phase. The surface geometry of the shell is used to mesh elements upon the shell surface. (x2
) and (x1
) are shear correction terms expressed in the form of Ru3; x2
and Ru3; x1
, respectively. Ru3 is expressed in the following polynomial form [8,9]: Ru3 = a1 x1
2 + a2 x2
2 + a3 x1
x2
+ a4 x1
+ a5 x2
+ a6 :

(8)

The unknown constant coe:cients ai (i = 1 − 6) are obtained by expressing the transverse shear strain in natural coordinate system (r; s; t, Fig. 2) as shown in the following expression: ∗r = r + Rr ; where r = u3; r + ’n − Rr = Ru3; r ;

(9) ur ; Rr

(10) (11)

r is the transverse shear strain along r. ’n is the rotation of normal along the direction r, and ur is the in-plane displacement along r. Rr is the radius of curvature along the direction r. The constant coe:cients of the expression Ru3 in Eq. (8) are obtained satisfying ∗r; n = 0

(12a)

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

265

Fig. 11. Variations of normalized frequency with respect to mode shapes for a panel with a=h=1000, and R1 =a=R2 =b=40.

and Ru3 = 0

(12b)

at three nodes. The displacements in terms of nodal values with respect to local coordinate system (x1
; x2
; x3
) are {u
} = Si uj
i ;

(13)

where Si , shape functions at the reference surface of the shell element, are given as follows: S1 = (1 − r − s)ei!t ;

(14a)

S2 = (r)ei!t ;

(14b)

S3 = (s)ei!t ;

(14c)

where ! and t denote frequency and time, respectively. Eq. (13) can now be expressed as 
 u1       1  
     u   2         u  

P P P {u } = u3 = S 1 S 2 S 3 ; u2    3 
       u4   u   
 u5 where



S1 0   SP 1  =  0  0 0

0 S1 0 0 0

0 0 S1 0 0

0 0 0 S1 0

 0 0   0 ;  0 S1

(15)

(16a)

266

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 12. Mode shape 1 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (KEYA-FREQ).

Fig. 13. Mode shape 2 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (KEYA-FREQ).

Fig. 14. Mode shape 3 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (KEYA-FREQ).

Fig. 15. Mode shape 4 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (KEYA-FREQ).

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

267

Fig. 16. Mode shape 5 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (KEYA-FREQ).

Fig. 17. Mode shape 6 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (KEYA-FREQ).

Fig. 18. Mode shape 7 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (KEYA-FREQ).



S2 0  SP 2  =  0 0 0  S3 0  SP 3  =  0 0 0

0 S2 0 0 0

0 0 S2 0 0

0 0 0 S2 0

0 S3 0 0 0

0 0 S3 0 0

0 0 0 S3 0

 0 0  0 ; 0 S2  0 0  0 ; 0 S3

u 1  = { u1
1 u2
1 u3
1 u4
1 u5
1 }T ;

(16b)

(16c)

(17a)

268

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 19. Mode shape 1 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (NASTRAN).

u 2  = { u1
2 u2
2 u3
2 u4
2 u5
2 }T ;

(17b)

u 3  = { u1
3 u2
3 u3
3 u4
3 u5
3 }T ;

(17c)

where ui
j (i = 1–5) are displacements at node j. Strains in terms of local coordinate system of an element can be now expressed as   @S1 S1 1 0 x3
@S 0 @x1
R1 @x1    
@S1 @S1 S1 0 x 0

  3 @x2 @x2 R2    

@S1

@S1
@S1
@S1
 (1−c x3 ) @x
0 x3 @x
x3 @x
: : :  { } =  (1+c x3 ) @x2  1 2 1     R!∗T ()S ∗T ∗T ∗T ∗T S1 @S1 − +R! ()S +R! ()S R! ()S S +R! ()S  1 1 @x
1 1 1 1 R
2 2   ∗T ∗T ∗T ∗T ∗T S1 @S1 − R
+R! ()S1 R! ()S1 +R! ()S1 S1 +R! ()S1 R! ()S1 5×15 @x
1

1

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

269

Fig. 20. Mode shape 2 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (NASTRAN).

 1  u        2 ; u       3  u  15×1 {
} = [B]ui
j ;

(18)

(19)

where R!∗T contains respected terms as obtained from the operation, substitution, and proper transformation of Eqs. (8) – (13) to local orthogonal Cartesian coordinate system of the element.

270

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 21. Mode shape 3 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (NASTRAN).

Introduction of Galerkin’s approach provides the following expression to the generation of elemental stiGness matrix for the dynamic analysis:  v

T




[B] [C][B]{U } dv + !

2

 A

[M ]{U
} dA = 0:0;

(20)

where [C] is a constitutive matrix for materials, and  1  u       
2 {U } = u  ;     3   u 

(21)

P SP 1  SP 2  SP 3  ]; [M ] = [ SP 1  SP 2  SP 3  ]T [][

(22)

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

271

Fig. 22. Mode shape 4 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (NASTRAN).

where



P1 +

2P2 R
1

    0    0 [] P =   P3   P2 +
 R1   0

0

0

P1 +

0

0

0

P1

0

0

0

P3

0

0

P+

P2 +

P2 R2

P3 R
2

P2 R
1

0



  P2  P1 +
 R2    : 0    0     P3

(23)

A full integration scheme [8,9] is used to obtain the element stiGness, and mass matrices. The full assembled matrices are linked with IMSL [16] to solve eigenvalue and eigenvector problems. The mode shapes are plotted using SURFER [17].

272

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 23. Mode shape 5 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (NASTRAN).

3. Numerical results and discussions The following example problems, pertaining to spherical panels of square planform—which are special cases of doubly curved panels—are considered to illustrate the present "nite element procedure: a=b = 1; a=h = 10; 100 and 1000; R1 =a = R2 =b = 5; 10; 20 and 40: Boundary conditions: All translations and rotations are restrained at the edges. a=h = 1000 is considered to study the performance of the element in a very thin con"guration.

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

273

Fig. 24. Mode shape 6 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (NASTRAN).

For the sake of convenience the following nondimensionalized frequency is de"ned:  !i∗ = !i a2 =E=h; where !i , i = 1; 2; 3; : : : represents the natural frequencies in numerical orders of magnitudes. The present results are compared with the available analytical solutions [14], and commercially available "nite element packages [18,19]. The analytical solutions [14], suitable for moderately thick shallow shell panels, are obtained using assumed displacement functions as proposed for the shell theory considered in the present study. The displacement-based solution functions are assumed in the following form: ui =

∞  ∞ 

Aimn sin(*m x1 ) sin(+n x2 )ei!t

(i = 1; 2; 3; 4; 5);

m=1 n=1

where Aimn denotes Fourier constants. *m and +n are de"ned as m,=a and n,=b, respectively. Converged results are obtained with m = n ¿ 5.

274

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 25. Mode shape 7 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (NASTRAN).

The shell element Shell43—Plastic Shell of ANSYS [18] is considered to compare the present results. The three-node triangular element is based on single-point integration on in-plane, and two-point integration along the through thickness. Mixed interpolation functions as of Bathe and Dvorkin [2], and Dvorkin and Bath [1] are considered along with the vertex rotations of Allman [20] in the element formulation of ANSYS [18]. The three-node triangular shell element CTRIA3 of NASTRAN is considered to compare the present results. The following executable commands are used in NASTRAN data deck: PARAM, K6ROT, 2.0 PARAM, SNORM, 2.0 The present methodology is referred to as KEYA-FREQ. The spherical panel considered is meshed with 32 × 32 elements, as shown in Fig. 3, though few elements are su:cient for the convergence. First lowest seven eigenvalues and mode shapes are obtained and compared. Fig. 4 plots variations of fundamental frequency, obtained from the analytical solution, ANSYS, NASTRAN, and present KEYA-FREQ with respect to radius/span for a=h = 10. The present solution

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

275

Table 1 Comparisons of mode shapes for a doubly curved shell with R1 =a = R2 =b = 40, and a=h = 1000 Mode number (eigenvalue)

KEYA-FREQ

NASTRAN

ANSYS

1 (!1 )

1. Peak value at the center of the shell panel 2. Diagonal symmetry

1. Peak value at the center of the shell panel 2. Symmetry about the center of the shell

Peak values at a1 (min) and a2 (max)

2 (!2 )

Same as NASTRAN

Peaks at the center of zones A(min) and D(max)

Peaks near the location b2 (max) and at the center of zones A(min) and B(min)

3 (!3 )

Same as NASTRAN

Peaks at the center of zones B(min) and C(max)

Peaks at b1 (max) and b2 (min)

4 (!4 )

Peaks at the center of the zones A and D(max), and at the center of the shell panel (min)

Peaks at the center of zones A(max), B(min), C(min), and D(max)

Same as NASTRAN

5 (!5 )

Similar to NASTRAN

Similar to NASTRAN

6 (!6 )

Peaks near the center of zones B(max) and D(max), and at the center of shell (min) Similar to ANSYS

Peaks at a1 (max) and a2 (max), and b1 (min) and b2 (min) Peak at the center (max); and at the center of zones A, B, C, and D(min)

7 (!7 )

1. Peaks at zones B(max) and C(min) 2. Pattern similar to NASTRAN

Peaks at zones A(max) and D(min)

Similar to KEYA-FREQ

Same as NASTRAN

agrees closely to the analytical one, as they are formed from the same roots. Results obtained from NASTRAN remain little far away then ANSYS with respect to the present method. Similar trends are also observed for the case of a=h=100 and 1000 as shown in Figs. 5 and 6, respectively. Similar variations are plotted for the second lowest frequencies for the case of a=h = 10, 100, and 1000 in Figs. 7–9, respectively. Figs. 10 and 11 plot frequencies with respect to mode shapes (1–7) for a shell with a=h = 100 and 1000, respectively. In Fig. 10, one can see the closeness of values of !2∗ and !3∗ , and !5∗ and !6∗ for all the "nite element results for a=h = 100. However, for a=h = 1000 (Fig. 11) KEYA-FREQ and NASTRAN show !2∗ ≈ !3∗ , and !1∗ ≈ !2∗ ≈ !3∗ , respectively. This trend is not noticed in the results of ANSYS.

276

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 26. Mode shape 1 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (ANSYS).

Fig. 27. Mode shape 2 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (ANSYS).

Seven mode shapes for the corresponding lowest frequencies are presented from each analysis group for a shell with a=h = 1000, and R1 =a = R2 =b = 40 for an illustrative purpose. Mode shapes in Figs. 12–18, 19 –25, 26 –32 are for KEYA-FREQ, NASTRAN, and ANSYS, respectively. Plots for KEYA-FREQ are obtained using SURFER [17]. In order to compare the mode shapes, the shell

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

277

Fig. 28. Mode shape 3 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (ANSYS).

geometry is divided into four zones (Fig. 33) and locations as de"ned below: Zone A: 0:0 6 x1 =a 6 0:5; 0:0 6 x2 =b 6 0:5; Zone B: 0:5 6 x1 =a 6 1:0; 0:0 6 x2 =b 6 0:5; Zone C: 0:0 6 x1 =a 6 0:5; 0:5 6 x2 =b 6 1:0; Zone D: 0:5 6 x1 =a 6 1:0; 0:5 6 x2 =b 6 1:0; Location a1 : x1 =a = 0:25; x2 =b = 0:50; Location a2 : x1 =a = 0:75; x2 =b = 0:50; Location b1 : x1 =a = 0:50; x2 =b = 0:25; Location b2 : x1 =a = 0:50; x2 =b = 0:75: Table 1 presents a comparison of mode shapes obtained from KEYA-FREQ, NASTRAN, and ANSYS. Noticeable disparities are observed in few mode shapes among all "nite element methods for the example considered. 4. Conclusion A three-node shell element with a constant curvature and full integration scheme has been developed using Sanders’ moderately deep shell theory suitable for moderately thick and thin thickness con"gurations. The inclusion of curvature eGects in the element formulation did not show any locking phenomena in thin shell range. An improvement in frequency response has been achieved in comparison with the analytical solution for a thick shell. Seven mode shapes corresponding to the lowest frequencies have been studied and compared with NASTRAN and ANSYS for a thin shell con"guration, and noticeable disparities are observed in few mode shapes among them.

278

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 29. Mode shape 4 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (ANSYS).

Fig. 30. Mode shape 5 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (ANSYS).

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 31. Mode shape 6 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (ANSYS).

Fig. 32. Mode shape 7 for a panel with a=h = 1000, and R1 =a = R2 =b = 40 (ANSYS).

279

280

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

Fig. 33. Plan of a panel with zones and locations.

References [1] E.N. Dvorkin, K.J. Bathe, A continuum mechanics based four-node shell element for general nonlinear analysis, Eng. Comput. 1 (1984) 77–88. [2] K.J. Bathe, E.H. Dvorkin, Short communication—a four-node plate bending element based on Mindlin/Reissner plate theory and mixed interpolation, Int. J. Numer. Methods Eng. 21 (1985) 367–383. [3] O.C. Zienkiewicz, D. Lefebvre, Three-"eld mixed approximation and the plate bending problem, Comm. Appl. Numer. Methods 3 (1987) 301–309. [4] P. Papadopoulus, R.L. Taylor, A triangular element based on Reissner-Mindlin plate theory, Int. J. Numer. Methods Eng. 30 (1990) 1029–1049. [5] M.H. Verwwood, A.W.M. Kok, A shear locking free six-node Mindlin plate bending element, Comput. Struct. 36 (1990) 547–555. [6] O.C. Zienkiewicz, R.L Taylor, P. Papadopoulos, E. Onate, Plate bending elements with discrete constraints, new triangular elements, Comput. Struct. 35 (1990) 505–522. [7] A. Tessler, A C 0 anisotropic three node shallow shell element, Comput. Method Appl. Mech. Eng. 78 (1990) 89–103. [8] H.R.H. Kabir, A shear-locking free isoparametric three-node triangular element for moderately-thick and thin plates, Comput. Struct. 35 (1992) 503–519. [9] H.R.H. Kabir, A shear-locking free robust isoparametric three-node triangular element for general shells, Comput. Struct. 51 (1994) 425–436. [10] E. Onate, F. Zarate, F. Flores, Simple triangular element for thick and thin plate and shell analysis, Int. J. Numer. Methods Eng. 37 (1994) 2569–2582. [11] R.M. Sydenstricker, L. Landau, Study of some triangular discrete Reissner–Mindlin plate shell elements, Comput. Struct. 78 (2000) 21–33. [12] J.L. Sanders Jr., An improved "rst-approximation theory for thin shells, NASA Technical Report R-24, 1959. [13] J.N. Reddy, Exact solutions of moderately thick laminated shells, J. Eng. Mech. ASCE 110 (1984) 794–809.

H.R.H. Kabir / Finite Elements in Analysis and Design 39 (2003) 257 – 281

281

[14] H.R.H. Kabir, R.A. Chaudhuri, On Gibbs-phenomenon-free Fourier solution for "nite shear-(exible laminated clamped curved panels, Int. J. Eng. Sci. 32 (1994) 501–520. [15] E. Reissner, On the theory of bending of elastic plates, J. Math. Phys. 23 (1944). [16] IMSL Inc., Math/Library, Houston, TX, USA, 1991. [17] SURFER, Scienti"c Software Corp., VA, USA, 1994. [18] ANSYS Version 5.6, Canonsburg, PA, USA, 1999. [19] MSC/NASTRAN Version 70.7, Macneal–Schwendler Corporation, CA, USA, 2000. [20] D.J. Allman, A compatible triangular element including vertex rotations for plane elasticity analysis, Comput. Struct. 19 (1984) 1–8.