- Email: [email protected]
A finite element formulation for moderately thick shells of general shape
Compurm & Srru~turcs Vol.54. No. I, pp. 49-57. 1995 Copyright 4~’ 1994 Elsevier Science Ltd
Printedin GreatBritain.Allrightsreserved cW5-7949/95 19.50+ 0.00
0045-7949(94)00304-l
A FINITE
ELEMENT FORMULATION FOR MODERATELY THICK SHELLS OF GENERAL SHAPE N. Kumbasar and T. Aksu
Department of Civil Engineering, Istanbul Technical University, 80626 Ma&k-Istanbul,
Turkey
(Received 17 October 1993) Abstract-An
deformation formulation in this study solutions of
expression for the potential energy of a shell of general shape, including thickness shear and without neglecting z/R in comparison with unity, is derived. Then a finite element based on this functional is obtained. Furthermore, a finite element shape function presented contains high degree terms to describe the surface geometry better. Accuracy of finite element moderately thick shells of general shape is shown in examples.
effects have to be considered for moderately thick expansions or basic shells, while asymptotic equations of elasticity must be used for thicker ones. There are various solution methods developed especially for computer applications, where thickness shear deformation is taken into account and the first three terms of power series of l/(1 + z/R) are considered for integration, as (1 - z/R + z */R *) [4]. The most important part of the error, that originates from the Kirchoff-Love hypothesis, is due to the neglect of thickness shear deformation. Inclusion of the shear deformation will thus add terms of order h/R. Since the neglect of the term z/R in (1 + z/R ) causes an error of the same order, the inclusion of this term may be considered as consistent and necessary. In this study, an expression of potential energy for moderately thick shells of general shape is obtained, including the effects of thickness shear deformations and considering the term z/R, generally neglected in the classical shell theory. Then a finite element formulation based on this functional is obtained.
NOTATION she11coordinates Lame parameters the principal radii of curvature of the middle surface of a shell the components of the displacement vector in the corresponding orthogonal coordinate directions tangential and normal displacements of the middle surface the rotations of tangents to the middle surface oriented along the parametric lines a, and a2 respectively the normal and shearing strains the thickness shearing strains the normal and shearing strains of the middle surface the thickness shearing strains of the middle surface the changes of curvature of the middle surface the normal stresses on two mutually perpendicular faces the shear stresses on the faces the thickness shear stresses Young’s modulus; Poisson’s ratio; shearing modulus of elasticity the thickness of shells the strain energy the strain energy of membrane terms the strain energy of bending and additional
2. STRAIN-DISPLACEMENTRELATIONS
terms the strain energy of shear terms
Strain-displacement relations obtained in threedimensional elasticity for small displacements can be expressed in terms of shell coordinates a,, a2, z and
local coordinate system local coordinates of node i shape functions of the shell element shape functions with high degree terms to describe the surface geometry
the displacements U,V, W along these coordinate directions [5]. Assuming that plane sections remain plane after deformation, the strains L, , t2, Y,~,y,_, y2._ may be written in terms of the displacements of middle surface U,U,Wand the rotations a, /I to obtain:
1. INTRODUCTION
methods are used for analysis of shells to obtain adequate solutions, with an increasing thickness-radius ratio. In contrast, the effect of thickness shear deformations and the ratio z/R in (1 + z/R) are neglected for thin shells [l, 2, 31. These 49
50
N. Kumbasar and T. Aksu
where the strains of middle surface and the curvatures are defined as 1
co -
u.1 +
’ - A,(1 +2/R,)
’
2
V-k
M’
R, +z
+zlR,)
1
dz __
=ln
mh,‘? R, + z
where
A ?.I
s fh?
-
0
=iT-Yjp co =
1
‘41.1 A,A,U
used for solutions that have the desired degree of accuracy. Some typical integrations are given below:
Similar integrations as follows:
J (-+I?:2
-,,,2
.,
s s s fh2
1 1 + z/R,
-,z,‘2
A
” = A,(1 + z/R,)
“’ + A, A,(1 :: z,R,)
B
XI
A
B.2 A,(1
+
z/R,)
1 +z/R,‘~ a
A,(1 +zlR,) ‘= A,(1 + z/R2)
4=-r,
1
1 + z/R,
A,(1 + z/R,)
1 +--r 1 +zlR,
i
+
R, +
z
+
z
R,
+
z
z dz = (R, - R,)hZ,
= R,RihZ,.
>.2
+2(R,-R,);Zc’
2 I ,jl~+R:Z,z:+~i~
+Z,
c;‘+2(R,-
’
+R~Z2~:+2vc~c~+2v-~,~, These expressions are the same as those used by Novozhilov [4], except those for the thickness shear deformations.
+;(I
-v)yY2+f(l
R,);Z,C;X~ I h2 12
-I$
-v)y;2+;(l
+);
3. STRAIN ENERGY OF MODERATELY THICK SHELLS
+(1 -v)$(R, Considering the strains obtained in Section 2, together with Hooke’s Law and inserting in the expression .(a,c,+a,~,+s,,y,,+T,,~,.-+t2_~z.-)dl’,
(3)
various types of integrals, depending on z, are encountered. It is possible to integrate them with respect to z and find exact expressions instead of expanding l/( 1 + z/R ) in power series. They may be
1 2
z,yj2
>
- R,)Z,y:r, 2
+(I b
>
;I(
z,yy2
1
I2 Cl
+i(l
v=;
(6)
The strain energy expression may be formed after performing all of these integrations. Integrating over z and rearranging all the terms of strain energy [6] in eqn (3) a field integral is obtained in the result:
1
=p
R
+h’2R + z Lz’dz m,,,>R, + z
1
Y2=A,A2(l:z,R2)a+
in terms of Z,
dz = -hZ, z
I
_h,2
1 + z/R,
i-i
R, +
+11.‘2 R
1
=p
may be obtained
(5)
-dz=[2+(2-l)Z,]h
r:
1
(4)
+;(&)o+;($J+..
;
1 -0 _ -----y:++ Y I2 - 1 + z/R,
h
z,=;(&>‘+;(&)
1
A,A2(l+z/R2)u+AZ(l+z/R2)V~2+~w 1 E-6 1 + z/R,
1 +hP-R,
~=hiZR,(l+Z,)>
+f(l
-v);(R2-
I
-v)Z,R:r:+f(l
h2 +(1 -v)7j~172+(1 ++(I
-v)(y:;+y;;)
R,)Z,y;r, -v)Z,R;r; -V)yYY! dA.
(7)
51
Finite element formulation In order to obtain an expression similar to those obtained in previous work [4], one may write,
v=
similar to Q2 given by Novozhilov, and retain the terms consistent with the accuracy of the present analysis, one may write
Eh’
Eh 2(1 -v2)
sA
QodA + 24(1 -v2)
z, Q3 dA,
Q,dA + Eh 2(1 -v2) sR
sA
(81
h2 2i ~
12R;’
(11)
Thus the following is obtained:
in which Q3 represents the effects of thickness shear stress. Q,=(er)+~;)~-2(1
-v)(&$),
(9)
which is the same expression given by Novozhilov. In order to make the expression
cy2+2(R, - R,);Z,&
Q2 = 2,
(12)
2
ry2
+R;Z,x:+Z,
Q,=;(l-v,(y;;+y;:,.
(13)
The reason for the difference of Q, from that given by Novozhilov is the independence of the functions a and j from the displacements u, v, w, due to the inclusion of shear deformations. An additional strain energy term, Q3, is also obtained for the same reason.
+2(R2-R,);Z2r;~,+R;Z,d I
4. ISOPARAMETRIC
FINITE ELEMENT
FORMULATION
z,yj2 A curved isoparametric trapezoidal finite element formulation based on the expression for the potential energy of a shell of general shape, including thickness shear deformation and without neglecting z/R in (1 + z/R), is derived. The shell element that is used has eight nodes (Fig. 1) and the parent shape function is
h2
+(I -v&R,
- R,P,yh
+(I
-V)EV,
2
+(1 -v)$(R,-
R,)Z,y;r, 1
+f(l
-v)Z,R;r;+;(l
-v)Z,R:r;
Fig. I. Curved
(10)
trapezoidal
finite element.
N. Kumbasar and T. Aksu
52 at the corner
nodes
and
three parts as bending, shear and additional improved strain energy respectively:
Ci=“,N!=f(l
-C*)(l
+rlrl!)
q,=O,N,=f(l
+‘X,)(l
-q2)
(14)
at the middle nodes [7,8]. At each node of this refined element, three displacements, U,U,W, in the directions of the local coordinates and two rotations a, fl, i.e. a total of five parameters, are defined. The displacement vector (6 ‘}, that is constituted of nodes degrees of freedom, is given below:
So there are 5 degrees of freedom at each node and 40 for each element. The finite element stiffness matrix and load vector may be obtained directly from the matrix form of the potential energy. In this study, the element stiffness matrix is calculated in two different parts, due to the occurrence of additional matrices which are obtained from improved strain energy. [kl’=[k]f+[k]r,
(16)
where [k ]‘j and [k 1; correspond to element stiffness matrices of membrane, bending, shear effect and additional element stiffness matrices of improved strain energy respectively. To avoid shear locking in thin shells, some forms of reduced integration have been suggested. Therefore in this study, the [k If element stiffness matrix can be obtained by separating into bending terms and shear terms as follows:
wl;=wl’,+~~l~.
terms of
(17)
With these definitions, the following expressions for the element stiffness matrices may be obtained in
kl’, =
(18)
PlF[~lzPl,dA. sA
In selective reduced integration, the bending term and the additional term is integrated according to the normal rule, whereas the shear term is integrated according to a lower-order rule [9, 10, I I]. In this study, for some shells, shape functions of the shell finite element with eight nodes are used to describe the surface geometry similar to the calculation of element stiffness matrices. The level of approximation obtained for the geometrical values of some shells with this element is indicated as insufficient. For that reason, a finite element shape function presented in this study contains high degree terms to describe the surface geometry better by examining publications [7,9, 121. The approach in the geometry calculation is improved with shape functions that are found by increasing the number of points on two sides of the element and making use of the Lagrange interpolation functions (Fig. 2). These shape functions are, fl,=i(l
+<<,)(l
+tl$!)]rC +(l
-9V,)(-$r12Vf-!V11,-
l)l i = 1,2,3,4
at the corner
nodes and
W,=i(l
q, = 0.
IV,= 2( 1 + <<,) (rj * - 1) (q * - l/4)
9, = f l/2,
-{*)(l
+qrj,)
i=6,8
{,=O,
i = $7
N, = -3 (1 + 5C,)V (q * - 1) (V + rl,) i = 9,10,11,12
at the middle
nodes.
Fig. 2. Element with high degree terms.
(19)
Finite element
formulation
53
Fig. 3. Shallow spherical shell.
Table 1. Deflection w (at x = 0, y = 0) of a simply-supported shallow spherical shell under uniformly distributed load (w x 10)) 0.32
1.60
3.20
6.40
Proposed finite element
330.469
51.583
I 1.482
1.929
Analytical solution, Reddy [13]
313.860
49.695
1I.265
1.976
Thickness (h )
5.
NUMERICAL EXAMPLES
written in the Fortran A computer program, programming language is developed for the analysis of moderately thick shells of general shape. In order to test the effectiveness of the refined finite element, three examples were solved and the results were
compared with other studies. Example 1. As shown in Fig. 3, a shallow spherical shell is simply supported on four sides under uniform loading. The geometric and material parameters used are 1, = l? = 32.0 in., E/p = lo’,
R = 96.0 in., v = 0.30.
This example has been solved by Reddy [ 131 for thick shell theory. The results of deflections at the middle point due to the variation in thickness h are listed in Table 1. This problem is solved using the proposedfiniteelementmodeland thecomputedresults are given for comparison. In the calculations, only a
Table
2. The maximum SAP 90
w.4
P41 0.092982
displacement
quarter of the shell is considered and divided into 64 elements due to symmetry. In this problem, geometrical values of the shell are calculated with shape functions of the shell finite element with eight nodes. It can be Seen that the results of the proposed finite element solution approach progressively the results obtained by Reddy as the thickness of the shell increases. Example 2. The performance of the shell element for thin shells is also evaluated on a standard test problem [14, 151 of a hemispherical shell with a hole, presented in Fig. 4. In this problem, geometrical values of the shell are calculated with shape functions of the shell finite element with high degree terms. The results for the maximum displacement obtained with proposed finite element are given in Table 2. The radial deflection wAis compared with the analytical solution and the various finite element solutions. Considering the various finite element solutions, the very good performance of proposed finite element is clearly evident.
w, for a hemispherical
shell with an 18” hole
Degenerated element HMSHS [15]
Isoparametric element Q4-SRI [15]
Analytical solution 114, 151
0.09494
0.09212
0.093
Proposed finite element
0.09274
N. Kumbasar and T. Aksu
R = 10.0in h = 0.04in $=160 E = 6.825x10'psi v = 0.30 P =l.Olb(on quadrant)
Fig. 4. Hemispherical shell with an 18” hole.
R = 2.5m h = l.Om E=l.Ot/m* v = 0.30 P=l.0t/m2 4R
Example 3. The shell is an axisymmetric cup with a thickness to radius ratio of 2/5 and an inside depth to radius ratio of S/5 (Fig. 5). The internal radius of the cup bend is 0.05 times the shell radius and the lip of the cup is rigidly clamped. This vessel is loaded by a constant internal pressure. In the calculations, only half of the shell is taken into consideration with one row totalling 32 elements because of axial symmetry. Geometrical values of the shell are calculated with shape functions of the shell finite element with high degree terms. Tangential, longitudinal and shear stresses, which are obtained with the proposed finite element, aregiven with the results taken from [16] for comparison. The distributions of these stresses across the shell thickness are indicated in Figs 6-l 1 for two cross-sections. The refined finite element solution of moderately thick shells has been shown to be in good agreement with the example given by Barrett and Soler [ 161. 6. CONCLUSION
I
Fig. 5. Pressure vessel.
In this study, an improved strain energy expression for moderately thick shells is developed, where the effects of thickness shear deformations and the expressions involving (1 + z/R ) are included. A curved trapezoidal finite element formulation based on this functional is derived. The shell element has eight nodes. At each node, three displacements at the directions of the local coordinates and two
0.75 0.50
P
Axiwill 399 I161 A Axihot 157 1161
l Proposed Finite Element
-
h
0.25
-
0.00
-
-0.25
-
-0.50
-
-1.60 Fig. 6. Tangential
-1.40
-1.20
-1.00
-0.80
-0.60
-0.40
(x = 2.5 h, z = 0.8h ) plotted
stress at section
-0.20
across
a1
the thickness.
0.25 1
0 Axiwill399 [I61 A Axihot 157 [161 0 Proposed Finite ElenW~t
Fig. 7. Longitudinal
stress at section
(x = 2.5
h , z = 0.8 h) plotted across the thickness.
0 Axiwill 399 [Ia] A Axihot 157 [Ia]
l Proposed Finite Element
Fig. 8. Shear stress at section
(x = 2.5 h, i = 0.8 h ) plotted
across
the thickness.
56
N. Kumbasar 1
and T. Aksu
.oo
0.75
h
0.50
‘:
0.25
7
0.00 : -0.25
4
-0.50
i
Fig. 9. Tangential
h
0.50
‘:
0.25
3 1
0.00
:
-0.25
‘:
-0.50
z
-0.75
Fig. 10. Longitudinal
0.25
0.00
(x = 2.5 h, z = 4.2 h ) plotted
across
the thickness.
0
Auwill 399 [I61 A Axihot 157 [I61 Proposed Finite Element
0
‘:
-1.oo-nY -2.50
h
stress at section
I, 1-w I, I I I I, I I I I, 1 I, I, I, -2.00 -1.30 -1 .OO -0.50 0.00 stress at section
(x = 2.5 h, z = 4.2 h ) plotted
I ‘, 0.50
across
*2
the thickness.
0 Ax&II 399 [I61 I57 1161 0 Proposed Finite Element
A Axihot
-0.25 -0.50
Fig.
11. Shear stress at section (x = 2.5 h, z = 4.2 h ) plotted across the thickness
Finite element formulation rotations are defined as a total of five unknowns. Element matrices are obtained explicitly by using this element. The approach in the geometry calculation is improved with finite element shape functions that contain high degree terms. The level of approximation obtained for the geometrical values is more accurate than those of shape functions with quadratic terms. In the solution of a problem, data are entered in respect of the finite element with 12 nodes, that contain high degree terms. When displacements which belong to last four nodes are entered as zero,
using different elements does not cause any difficulty during data input. The comparison of the results with the examples given in the literature has resulted in good agreement. The finite element is applicable for both thin and moderately thick shell analysis. The element becomes lock-free in shear actions by using reduced numerical integration when the element is applied to thin shells. The most important advantage of this formulation is the applicability of the method to problems of shells of general shape.
57
W. Fliigge, Stresses in Shells, 2nd Edn. Springer-Verlag.
7. 8. 9. 10.
Berlin (1973). V. V. Novozhilov, Thin Shell Theory, 2nd Edn. WoltersNoordhoff, Groningen (1970). H. L. Langhaar, Energy Methocis in Applied Mechanics. John Wiley and Sons, New York (1962). N. Kumbasar and T. Aksu, A refined expression of potential energy for thick shells. Bull. Technical University of Istanbul 4, 95-105 (1992). 0. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, London (1979). E. B. Becker, G. F. Carey and J. T. Oden, Finite Elements: An Introduction. Prentice-Hall, New Jersey (1981). T. J. R. Hughes, The Finite Element Method. PrenticeHall, New Jersey (1987). D. Briassoulis, On the basics of the shear locking problem of C” isoparametric plate elements. Comput. Struct. 33, 169-185 (1989).
11. D. Briassoulis, The C” structural elements with no locking mechanisms and zero energy modes. Numeta 90, 86-93
(1990).
12. D. S. Kang and T. H. H. Pian, A 20-Dof hybrid stress general shell element. Comput. Struct. 30, 789-794 (1988).
13. J. N. Reddy, Solutions of moderately thick laminated shells. J. Engng Mech. 110, 794-809 (1984). 14. A. Ibrahimbegovi9 and E. L. Wilson, Thick shell solid finite elements with independent and rotation fields. Int. J. Num. Meth. Engng 31, 1393-1414 (1991).
REFERENCES 1.
A. L. Gol’denveizer, Theory of Elastic Thin Shells. Pergamon Press, New York (1961). 2. H. Kraus, Thin Elastic Shells. John Wiley and Sons, New York (1967).
15. A. F. Saleeb, T. Y. Chang and W. Graf, A quadrilateral shell element using a mixed formulation. Comput. Struct. 26, 787-803 (1987).
16. D. J. Barrett and A. Soler, A finite element model for thick walled axisymmetric shells. J. Pressure Vessel Technol. 104, 2 15-222 (1982).
Printedin GreatBritain.Allrightsreserved cW5-7949/95 19.50+ 0.00
0045-7949(94)00304-l
A FINITE
ELEMENT FORMULATION FOR MODERATELY THICK SHELLS OF GENERAL SHAPE N. Kumbasar and T. Aksu
Department of Civil Engineering, Istanbul Technical University, 80626 Ma&k-Istanbul,
Turkey
(Received 17 October 1993) Abstract-An
deformation formulation in this study solutions of
expression for the potential energy of a shell of general shape, including thickness shear and without neglecting z/R in comparison with unity, is derived. Then a finite element based on this functional is obtained. Furthermore, a finite element shape function presented contains high degree terms to describe the surface geometry better. Accuracy of finite element moderately thick shells of general shape is shown in examples.
effects have to be considered for moderately thick expansions or basic shells, while asymptotic equations of elasticity must be used for thicker ones. There are various solution methods developed especially for computer applications, where thickness shear deformation is taken into account and the first three terms of power series of l/(1 + z/R) are considered for integration, as (1 - z/R + z */R *) [4]. The most important part of the error, that originates from the Kirchoff-Love hypothesis, is due to the neglect of thickness shear deformation. Inclusion of the shear deformation will thus add terms of order h/R. Since the neglect of the term z/R in (1 + z/R ) causes an error of the same order, the inclusion of this term may be considered as consistent and necessary. In this study, an expression of potential energy for moderately thick shells of general shape is obtained, including the effects of thickness shear deformations and considering the term z/R, generally neglected in the classical shell theory. Then a finite element formulation based on this functional is obtained.
NOTATION she11coordinates Lame parameters the principal radii of curvature of the middle surface of a shell the components of the displacement vector in the corresponding orthogonal coordinate directions tangential and normal displacements of the middle surface the rotations of tangents to the middle surface oriented along the parametric lines a, and a2 respectively the normal and shearing strains the thickness shearing strains the normal and shearing strains of the middle surface the thickness shearing strains of the middle surface the changes of curvature of the middle surface the normal stresses on two mutually perpendicular faces the shear stresses on the faces the thickness shear stresses Young’s modulus; Poisson’s ratio; shearing modulus of elasticity the thickness of shells the strain energy the strain energy of membrane terms the strain energy of bending and additional
2. STRAIN-DISPLACEMENTRELATIONS
terms the strain energy of shear terms
Strain-displacement relations obtained in threedimensional elasticity for small displacements can be expressed in terms of shell coordinates a,, a2, z and
local coordinate system local coordinates of node i shape functions of the shell element shape functions with high degree terms to describe the surface geometry
the displacements U,V, W along these coordinate directions [5]. Assuming that plane sections remain plane after deformation, the strains L, , t2, Y,~,y,_, y2._ may be written in terms of the displacements of middle surface U,U,Wand the rotations a, /I to obtain:
1. INTRODUCTION
methods are used for analysis of shells to obtain adequate solutions, with an increasing thickness-radius ratio. In contrast, the effect of thickness shear deformations and the ratio z/R in (1 + z/R) are neglected for thin shells [l, 2, 31. These 49
50
N. Kumbasar and T. Aksu
where the strains of middle surface and the curvatures are defined as 1
co -
u.1 +
’ - A,(1 +2/R,)
’
2
V-k
M’
R, +z
+zlR,)
1
dz __
=ln
mh,‘? R, + z
where
A ?.I
s fh?
-
0
=iT-Yjp co =
1
‘41.1 A,A,U
used for solutions that have the desired degree of accuracy. Some typical integrations are given below:
Similar integrations as follows:
J (-+I?:2
-,,,2
.,
s s s fh2
1 1 + z/R,
-,z,‘2
A
” = A,(1 + z/R,)
“’ + A, A,(1 :: z,R,)
B
XI
A
B.2 A,(1
+
z/R,)
1 +z/R,‘~ a
A,(1 +zlR,) ‘= A,(1 + z/R2)
4=-r,
1
1 + z/R,
A,(1 + z/R,)
1 +--r 1 +zlR,
i
+
R, +
z
+
z
R,
+
z
z dz = (R, - R,)hZ,
= R,RihZ,.
>.2
+2(R,-R,);Zc’
2 I ,jl~+R:Z,z:+~i~
+Z,
c;‘+2(R,-
’
+R~Z2~:+2vc~c~+2v-~,~, These expressions are the same as those used by Novozhilov [4], except those for the thickness shear deformations.
+;(I
-v)yY2+f(l
R,);Z,C;X~ I h2 12
-I$
-v)y;2+;(l
+);
3. STRAIN ENERGY OF MODERATELY THICK SHELLS
+(1 -v)$(R, Considering the strains obtained in Section 2, together with Hooke’s Law and inserting in the expression .(a,c,+a,~,+s,,y,,+T,,~,.-+t2_~z.-)dl’,
(3)
various types of integrals, depending on z, are encountered. It is possible to integrate them with respect to z and find exact expressions instead of expanding l/( 1 + z/R ) in power series. They may be
1 2
z,yj2
>
- R,)Z,y:r, 2
+(I b
>
;I(
z,yy2
1
I2 Cl
+i(l
v=;
(6)
The strain energy expression may be formed after performing all of these integrations. Integrating over z and rearranging all the terms of strain energy [6] in eqn (3) a field integral is obtained in the result:
1
=p
R
+h’2R + z Lz’dz m,,,>R, + z
1
Y2=A,A2(l:z,R2)a+
in terms of Z,
dz = -hZ, z
I
_h,2
1 + z/R,
i-i
R, +
+11.‘2 R
1
=p
may be obtained
(5)
-dz=[2+(2-l)Z,]h
r:
1
(4)
+;(&)o+;($J+..
;
1 -0 _ -----y:++ Y I2 - 1 + z/R,
h
z,=;(&>‘+;(&)
1
A,A2(l+z/R2)u+AZ(l+z/R2)V~2+~w 1 E-6 1 + z/R,
1 +hP-R,
~=hiZR,(l+Z,)>
+f(l
-v);(R2-
I
-v)Z,R:r:+f(l
h2 +(1 -v)7j~172+(1 ++(I
-v)(y:;+y;;)
R,)Z,y;r, -v)Z,R;r; -V)yYY! dA.
(7)
51
Finite element formulation In order to obtain an expression similar to those obtained in previous work [4], one may write,
v=
similar to Q2 given by Novozhilov, and retain the terms consistent with the accuracy of the present analysis, one may write
Eh’
Eh 2(1 -v2)
sA
QodA + 24(1 -v2)
z, Q3 dA,
Q,dA + Eh 2(1 -v2) sR
sA
(81
h2 2i ~
12R;’
(11)
Thus the following is obtained:
in which Q3 represents the effects of thickness shear stress. Q,=(er)+~;)~-2(1
-v)(&$),
(9)
which is the same expression given by Novozhilov. In order to make the expression
cy2+2(R, - R,);Z,&
Q2 = 2,
(12)
2
ry2
+R;Z,x:+Z,
Q,=;(l-v,(y;;+y;:,.
(13)
The reason for the difference of Q, from that given by Novozhilov is the independence of the functions a and j from the displacements u, v, w, due to the inclusion of shear deformations. An additional strain energy term, Q3, is also obtained for the same reason.
+2(R2-R,);Z2r;~,+R;Z,d I
4. ISOPARAMETRIC
FINITE ELEMENT
FORMULATION
z,yj2 A curved isoparametric trapezoidal finite element formulation based on the expression for the potential energy of a shell of general shape, including thickness shear deformation and without neglecting z/R in (1 + z/R), is derived. The shell element that is used has eight nodes (Fig. 1) and the parent shape function is
h2
+(I -v&R,
- R,P,yh
+(I
-V)EV,
2
+(1 -v)$(R,-
R,)Z,y;r, 1
+f(l
-v)Z,R;r;+;(l
-v)Z,R:r;
Fig. I. Curved
(10)
trapezoidal
finite element.
N. Kumbasar and T. Aksu
52 at the corner
nodes
and
three parts as bending, shear and additional improved strain energy respectively:
Ci=“,N!=f(l
-C*)(l
+rlrl!)
q,=O,N,=f(l
+‘X,)(l
-q2)
(14)
at the middle nodes [7,8]. At each node of this refined element, three displacements, U,U,W, in the directions of the local coordinates and two rotations a, fl, i.e. a total of five parameters, are defined. The displacement vector (6 ‘}, that is constituted of nodes degrees of freedom, is given below:
So there are 5 degrees of freedom at each node and 40 for each element. The finite element stiffness matrix and load vector may be obtained directly from the matrix form of the potential energy. In this study, the element stiffness matrix is calculated in two different parts, due to the occurrence of additional matrices which are obtained from improved strain energy. [kl’=[k]f+[k]r,
(16)
where [k ]‘j and [k 1; correspond to element stiffness matrices of membrane, bending, shear effect and additional element stiffness matrices of improved strain energy respectively. To avoid shear locking in thin shells, some forms of reduced integration have been suggested. Therefore in this study, the [k If element stiffness matrix can be obtained by separating into bending terms and shear terms as follows:
wl;=wl’,+~~l~.
terms of
(17)
With these definitions, the following expressions for the element stiffness matrices may be obtained in
kl’, =
(18)
PlF[~lzPl,dA. sA
In selective reduced integration, the bending term and the additional term is integrated according to the normal rule, whereas the shear term is integrated according to a lower-order rule [9, 10, I I]. In this study, for some shells, shape functions of the shell finite element with eight nodes are used to describe the surface geometry similar to the calculation of element stiffness matrices. The level of approximation obtained for the geometrical values of some shells with this element is indicated as insufficient. For that reason, a finite element shape function presented in this study contains high degree terms to describe the surface geometry better by examining publications [7,9, 121. The approach in the geometry calculation is improved with shape functions that are found by increasing the number of points on two sides of the element and making use of the Lagrange interpolation functions (Fig. 2). These shape functions are, fl,=i(l
+<<,)(l
+tl$!)]rC +(l
-9V,)(-$r12Vf-!V11,-
l)l i = 1,2,3,4
at the corner
nodes and
W,=i(l
q, = 0.
IV,= 2( 1 + <<,) (rj * - 1) (q * - l/4)
9, = f l/2,
-{*)(l
+qrj,)
i=6,8
{,=O,
i = $7
N, = -3 (1 + 5C,)V (q * - 1) (V + rl,) i = 9,10,11,12
at the middle
nodes.
Fig. 2. Element with high degree terms.
(19)
Finite element
formulation
53
Fig. 3. Shallow spherical shell.
Table 1. Deflection w (at x = 0, y = 0) of a simply-supported shallow spherical shell under uniformly distributed load (w x 10)) 0.32
1.60
3.20
6.40
Proposed finite element
330.469
51.583
I 1.482
1.929
Analytical solution, Reddy [13]
313.860
49.695
1I.265
1.976
Thickness (h )
5.
NUMERICAL EXAMPLES
written in the Fortran A computer program, programming language is developed for the analysis of moderately thick shells of general shape. In order to test the effectiveness of the refined finite element, three examples were solved and the results were
compared with other studies. Example 1. As shown in Fig. 3, a shallow spherical shell is simply supported on four sides under uniform loading. The geometric and material parameters used are 1, = l? = 32.0 in., E/p = lo’,
R = 96.0 in., v = 0.30.
This example has been solved by Reddy [ 131 for thick shell theory. The results of deflections at the middle point due to the variation in thickness h are listed in Table 1. This problem is solved using the proposedfiniteelementmodeland thecomputedresults are given for comparison. In the calculations, only a
Table
2. The maximum SAP 90
w.4
P41 0.092982
displacement
quarter of the shell is considered and divided into 64 elements due to symmetry. In this problem, geometrical values of the shell are calculated with shape functions of the shell finite element with eight nodes. It can be Seen that the results of the proposed finite element solution approach progressively the results obtained by Reddy as the thickness of the shell increases. Example 2. The performance of the shell element for thin shells is also evaluated on a standard test problem [14, 151 of a hemispherical shell with a hole, presented in Fig. 4. In this problem, geometrical values of the shell are calculated with shape functions of the shell finite element with high degree terms. The results for the maximum displacement obtained with proposed finite element are given in Table 2. The radial deflection wAis compared with the analytical solution and the various finite element solutions. Considering the various finite element solutions, the very good performance of proposed finite element is clearly evident.
w, for a hemispherical
shell with an 18” hole
Degenerated element HMSHS [15]
Isoparametric element Q4-SRI [15]
Analytical solution 114, 151
0.09494
0.09212
0.093
Proposed finite element
0.09274
N. Kumbasar and T. Aksu
R = 10.0in h = 0.04in $=160 E = 6.825x10'psi v = 0.30 P =l.Olb(on quadrant)
Fig. 4. Hemispherical shell with an 18” hole.
R = 2.5m h = l.Om E=l.Ot/m* v = 0.30 P=l.0t/m2 4R
Example 3. The shell is an axisymmetric cup with a thickness to radius ratio of 2/5 and an inside depth to radius ratio of S/5 (Fig. 5). The internal radius of the cup bend is 0.05 times the shell radius and the lip of the cup is rigidly clamped. This vessel is loaded by a constant internal pressure. In the calculations, only half of the shell is taken into consideration with one row totalling 32 elements because of axial symmetry. Geometrical values of the shell are calculated with shape functions of the shell finite element with high degree terms. Tangential, longitudinal and shear stresses, which are obtained with the proposed finite element, aregiven with the results taken from [16] for comparison. The distributions of these stresses across the shell thickness are indicated in Figs 6-l 1 for two cross-sections. The refined finite element solution of moderately thick shells has been shown to be in good agreement with the example given by Barrett and Soler [ 161. 6. CONCLUSION
I
Fig. 5. Pressure vessel.
In this study, an improved strain energy expression for moderately thick shells is developed, where the effects of thickness shear deformations and the expressions involving (1 + z/R ) are included. A curved trapezoidal finite element formulation based on this functional is derived. The shell element has eight nodes. At each node, three displacements at the directions of the local coordinates and two
0.75 0.50
P
Axiwill 399 I161 A Axihot 157 1161
l Proposed Finite Element
-
h
0.25
-
0.00
-
-0.25
-
-0.50
-
-1.60 Fig. 6. Tangential
-1.40
-1.20
-1.00
-0.80
-0.60
-0.40
(x = 2.5 h, z = 0.8h ) plotted
stress at section
-0.20
across
a1
the thickness.
0.25 1
0 Axiwill399 [I61 A Axihot 157 [161 0 Proposed Finite ElenW~t
Fig. 7. Longitudinal
stress at section
(x = 2.5
h , z = 0.8 h) plotted across the thickness.
0 Axiwill 399 [Ia] A Axihot 157 [Ia]
l Proposed Finite Element
Fig. 8. Shear stress at section
(x = 2.5 h, i = 0.8 h ) plotted
across
the thickness.
56
N. Kumbasar 1
and T. Aksu
.oo
0.75
h
0.50
‘:
0.25
7
0.00 : -0.25
4
-0.50
i
Fig. 9. Tangential
h
0.50
‘:
0.25
3 1
0.00
:
-0.25
‘:
-0.50
z
-0.75
Fig. 10. Longitudinal
0.25
0.00
(x = 2.5 h, z = 4.2 h ) plotted
across
the thickness.
0
Auwill 399 [I61 A Axihot 157 [I61 Proposed Finite Element
0
‘:
-1.oo-nY -2.50
h
stress at section
I, 1-w I, I I I I, I I I I, 1 I, I, I, -2.00 -1.30 -1 .OO -0.50 0.00 stress at section
(x = 2.5 h, z = 4.2 h ) plotted
I ‘, 0.50
across
*2
the thickness.
0 Ax&II 399 [I61 I57 1161 0 Proposed Finite Element
A Axihot
-0.25 -0.50
Fig.
11. Shear stress at section (x = 2.5 h, z = 4.2 h ) plotted across the thickness
Finite element formulation rotations are defined as a total of five unknowns. Element matrices are obtained explicitly by using this element. The approach in the geometry calculation is improved with finite element shape functions that contain high degree terms. The level of approximation obtained for the geometrical values is more accurate than those of shape functions with quadratic terms. In the solution of a problem, data are entered in respect of the finite element with 12 nodes, that contain high degree terms. When displacements which belong to last four nodes are entered as zero,
using different elements does not cause any difficulty during data input. The comparison of the results with the examples given in the literature has resulted in good agreement. The finite element is applicable for both thin and moderately thick shell analysis. The element becomes lock-free in shear actions by using reduced numerical integration when the element is applied to thin shells. The most important advantage of this formulation is the applicability of the method to problems of shells of general shape.
57
W. Fliigge, Stresses in Shells, 2nd Edn. Springer-Verlag.
7. 8. 9. 10.
Berlin (1973). V. V. Novozhilov, Thin Shell Theory, 2nd Edn. WoltersNoordhoff, Groningen (1970). H. L. Langhaar, Energy Methocis in Applied Mechanics. John Wiley and Sons, New York (1962). N. Kumbasar and T. Aksu, A refined expression of potential energy for thick shells. Bull. Technical University of Istanbul 4, 95-105 (1992). 0. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, London (1979). E. B. Becker, G. F. Carey and J. T. Oden, Finite Elements: An Introduction. Prentice-Hall, New Jersey (1981). T. J. R. Hughes, The Finite Element Method. PrenticeHall, New Jersey (1987). D. Briassoulis, On the basics of the shear locking problem of C” isoparametric plate elements. Comput. Struct. 33, 169-185 (1989).
11. D. Briassoulis, The C” structural elements with no locking mechanisms and zero energy modes. Numeta 90, 86-93
(1990).
12. D. S. Kang and T. H. H. Pian, A 20-Dof hybrid stress general shell element. Comput. Struct. 30, 789-794 (1988).
13. J. N. Reddy, Solutions of moderately thick laminated shells. J. Engng Mech. 110, 794-809 (1984). 14. A. Ibrahimbegovi9 and E. L. Wilson, Thick shell solid finite elements with independent and rotation fields. Int. J. Num. Meth. Engng 31, 1393-1414 (1991).
REFERENCES 1.
A. L. Gol’denveizer, Theory of Elastic Thin Shells. Pergamon Press, New York (1961). 2. H. Kraus, Thin Elastic Shells. John Wiley and Sons, New York (1967).
15. A. F. Saleeb, T. Y. Chang and W. Graf, A quadrilateral shell element using a mixed formulation. Comput. Struct. 26, 787-803 (1987).
16. D. J. Barrett and A. Soler, A finite element model for thick walled axisymmetric shells. J. Pressure Vessel Technol. 104, 2 15-222 (1982).