- Email: [email protected]
A shear deformable curved shell element of quadrilateral shape
Compurrrs & Srrucrurrs Vol Printed in Great Britain.
33. No. 4. pp. 987-1001.
1989
Oc45-7949/89 s3.00 + 0 cnl CT‘ 1989 Perpmon Press plc
A SHEAR DEFORMABLE CURVED SHELL ELEMENT QUADRILATERAL SHAPE H. V. LAKSHMINARAYANAand K. KA~LASH Structural Sciences Division, National Aeronautical Laboratory, Bangalore-
OF
017, India
(Received 20 October 1988) Abstract-Formulation and numerical evaluation of QUAD8**, a shear deformable curved shell element of quadrilateral shape, is presented. This element. which has eight nodes and six engineering degrees of freedom per node, is free from locking problems, has no spurious zero energy modes and provides numerical results of acceptable accuracy in a number of test and example problems. The standing of QUAD8** compared with alternatives should consider (1) the simplified input data, (2) the ease of use, (3) improved performance and (4) the simplicity of programming and computer implementation.
INTRODUCTION of simple to use and reliable curved shell elements has been and continues to be an active area of finite element research. Recent advances in this area are contained in a monograph [l] and a review paper [2]. Of particular relavance to the present work is a quadrilateral curved shell element existing in MSC/NASTRAN [3], whose performance was found to be disappointing [4]. This has motivated us to re-examine this element, identify and correct its weaknesses and offer a new element as an alterative. This report presents the formulation and numerical evaluation of a shear deformable curved shell element of quadrilateral shape designated QUAD8**. The assumed strain method is adopted here to alleviate the membrane and shear locking phenomena. Previous contributors to assumed strain curved shell elements include Park and Stanley [S], Huang and Hinton [6], Bathe and Dvorkin [7], Jang and Pinsky [8] and MacNeal [9]. The new element was demonstrated to be reliable and to lead to numerical results of acceptable accuracy in a proposed set of standard test problems [4]. The standing of QUADS**, when compared with alternatives, must take into account (1) the simplified input data, (2) the ease of use, (3) improved performance and (4) simplicity of programming and computer implementation. The report is organized as follows: the geometry of the subject element is defined first, followed by the co-ordinate systems used in the formulation. Succeeding sections describe the Jacobian matrix, displacement interpolations, definition of element strain matrices, constitutive relations, numerical integration formulae, the elastic stiffness matrix calculation, and alleviation of locking phenomena. Numerical solutions to a standard set of test problems and comparison of QUAD8** performance with alternatives are presented next. Of the many problems analysed with this element, four examples The development
987
are presented here: Cook’s plane stress problem, Morley’s skew plate problem, a cantilevered curved beam, and a spherical shell under internal pressure. We conclude by identifying directions for further work. FORMULATION Geometry and coordinate systems
Consider the typical shell element of Fig. l(a). The outer and inner surfaces of the element are curved, whereas the sections across the thickness are generated by straight lines. The geometric description of the element requires the specification of two vectors at each of the eight middle-surface nodes. One is the position vector Ri of the node, i, with the three global Cartesian components Xi, Yi, Z,, where the subscript I’ identifies the node number; the other is the unit normal vector V: along with the shell wall thickness ti of the same nodes. This geometric representation does not ensure continuity of the slope between adjacent elements. Therefore, a shell surface idealized with these elements cannot be exact. However, the accuracy of the approximation increases as the number of elements is increased. The element formulated for the analysis of shells with continuous as well as discontinuous curvatures carries six engineering degress of freedom (three translations and three rotations) at each of its eight mid-surface nodes. The nodal degrees of freedom are illustrated in Fig. l(b). Stiffness for rotation about the normal to the element middle surface is zero at all the nodes. Therefore, a small fictitious rotational spring is introduced at each node, corresponding to the rotation about the normal. Referring to Fig. l(b), let 5, q be the two curvilinear co-ordinates in the element middle surface and t a linear co-ordinate in the thickness direction. If we assumed that 5, q and r vary between - 1 and 1 on the respective faces of the element we can write
H. V. LAKSHMINARAYANA and K. KAILASH
988
(a)
1-4 4
i
7
3
-EL!-
2
8
6
5
1
l-
2
t
2
1
Fig. 1. A curved shell element of quadrilateral shape. (a) Element geometry, (b) coordinate systems and nodal degrees of freedom, (c) parent element and node numbering scheme.
a relationship between the global Cartesian coordinates (X, Y, Z) of any point and the natural co-ordinates (4, I, T) in the form
(1)
Here N,(<, q) is a shape function taking a value of unity at the nodes i and zero at all other nodes. They
are derived as shape functions of a parent element shown in Fig. l(c) and are given in Table 1. In the present formulation, the vector V; is not input by the user but is computed using the co-ordinates of the eight mid-surface nodes as described later. The shell thickness tin eqn (1) is assumed to be constant as its actual value is unimportant for calculating element strain matrices. Actual thickness and its variation within the element is introduced into the constitutive relations.
Table 1. Shape functions and their derivatives Edge nodes (5,7)
Edge nodes (6, 8)
a(1 + &)(I + %)(50+ ‘lo- 1)
f(l - ?)(l
f(1 - a2)(1 + 50)
N:,
i; (1 + %)(250 + %)
-5(l
NV
$(I +50)(2%+&I)
;(l-,F2)
Function
Corner nodes (1, 2, 3, 4)
N,
+ rlo)
+ ‘lo)
;(l -rl(l
-tq) + 5”)
lo = &; ‘lo= qq,; values of 5,. 1,. are defined with reference to the nodes in the diagram given below:
e
1
(-1.-l)
l
5 (O.-l)
2 (1.-l)
989
A shear deformable curved shell element of quadrilateral shape
Let us also define a local Cartesian co-ordinate - system (X, Y, Z) with 2 being normal to the element middle surface and 8, p tangent to it, as shown in Fig. I(b). The orientations of 1, y,z with respect to X, Y, Z are established uniquely using a procedure given in [3]. At any point (5, q) on the element mid-surface (7 = 0) we can define the following two vectors which are tangent to the mid-surface:
.
surface, i.e. IV,>
=
iv,
1 x
(6)
F21.
It is evident that Vi in eqn (1) can be computed using eqn (6). Finally, the local Cartesian to global Cartesian co-ordinate transformation matrix [Eg] is defined as
and This matrix is also used to transform the global displacements (U, V, W) to local displacements --(K K w)
The unit vectors of the local Cartesian co-ordinate system will be determined from the tangent vectors in an ‘unbiased’ fashion, i.e. in such a way that the orientation of the vectors is independent of the sequence of the grid point numbers, except for multiples of 90”. Define:
(8)
Jacobian matrix The transpose of the Jacobian matrix is defined as
-ax ax
ax-
Z&Z dx
[@ c=
axax away --+--+--
(
azaz
x aq at au at al?>i
I$,,
= [J] = 1
(2)
ay ay ay ~~~~ az az az -at; Zy aT_
(9)
By applying the derivatives to eqn (1) it is seen that
where I;, I,, are lengths corresponding to unit increments in 5 and q, and C is the cosine of the angle between the two tangent vectors. The first two unit vectors of the local Cartesian co-ordinate system are
rax
axi
x&w,Iv21=
ay
ay
x
F
44 [ -S/l,
-BP, 41, 1
(3)
az dz _x qi
+
Q and p are chosen to ensure that V, and V, have unit length and are perpendicular. 1 x =2
[
1 (I+(1 1
1 fi=,
[
(1--.
1 +(-)I/2
1
(1 + C)“2
1 1
(5)
We note that if the tangent vectors were perpendicular (i.e. if C = 0) then GL= 1, p = 0 (V, is tangent to the r direction and V, is tangent to the q direction). The third unit vector V, is normal to the middle
7;{3~N~~ N7i OJ).
(10)
In eqn (lo), Nti and N+ are the derivatives of shape functions N,(5, 9) with respect to 5, r~and the vector (ek ejy, e,), is the local normal at one of the eight nodes. In the present formulation, we will need the following two particular forms:
(11)
(12)
H.
990
Combining
V. ~AKSHMINA~AYANA and K. KAILASH
eqns (10)-(12),
matrix:
Definition of membrane strains, bending curvatures and transverse shear strains
Dispiucemmtt iflterpafution
The global displacements (U, V, W) are interpolated from the nodal point displacements and nodal point rotations using the same shape functions IV,{&q) as in eqn (1):
Strains and curvatures are defined in the local --- -. Cartesian co-ordinate system (A’, Y, Z). with A’, Y m the mid-surface and z being normal to it. The normal strain in the 2 direction 62, is neglected. for consistency with the zero normal stress assumption, The membrane strains ate
= -1im r-n
The bending curvatures are where at each mid-surface node i:
The parametric derivatives of the displacements U, V, W are given bY
The transverse shear strains are
au au ?j-zqlt 1 au
du
1
[ &f=
aw -i?w iiw
In the local Cartesian co-ordinate system (2, Y, .?), derivatives of displacements 0, p, k? with respect to X, r, 2 are given by the terms in the following
The evaluation of the membrane strains and the transverse shear strains requires the selection of appropriate terms in eqn (18), evaluated for 5 = 0. However, the evaluation of bending curvatures requires the differentiation with respect to r of terms in eqn (18). Detailed expressions for ail the shell strain components in terms of nodal degrees of freedom may be derived and are given by [3]:
A shear deformable curved shell element of quadrilateral shape
where
[CT”] =B,001 I 0
Bz
Bz
0 [E’]
B,
0
(24)
991
is the nodal dispia~ement vector and [BE], [BB] and [BS] are the so called element strain matrices. They are constructed as
(25)
(26) Element stress resultants and constitutit>e relations (27)
The stress resultants corresponding to the strains and curvatures are defined and related by constitutive relations in the usual manner. The general form of the constitutive relation used in the present formation is
(28)
(29) where (30)
(Nj=
Nxx N,, 1!NXP
(31)
are the membrane
forces per unit length,
(MJ= (32)
Mxx My, ! MKY1
are the bending moments per unit length, and
(33)
In evaIuating the above matrices, t/2 is set equal to unity in the calculation of [q and [,?I because the factor t/2 cancels everywhere. It is to be noted that the membrane strains [cm] do not depend on rotations, which indicates that the present formulation can be used in a shell analysis where bending stiffness is neglected. It is convenient to rewrite eqns (22)-(24) as follows: (34) (35) (36) where
and
are the transverse shear forces per unit length. We note that the element stresses have been defined in the local Cartesian co-ordinate system. Therefore it is also necessary to specify the elasticity matrices [A], [B], [D] and [E] in the same co-ordinate system. A,, B,, D, (i, J = 1, 2, 3) and Ei_,(i, J = 1, 2) are also called membrane, bending, membrane-bending coupling and transverse shear stiffness coefficients. The general form of eqn (40) is included in this element formulation to enable the analyst to model laminated anisotropic materials and sandwich construction in addition to homogeneous isotropic solids. Actual thickness and its variation within the element is accounted for in the calculation of [A], [B], [D], and [E] matrices. It is also a standard practice to include shear correction factors in evaluating the [E] matrix. When thermal loads are considered, eqn (40) has to be modified to include thermal forces and thermal moments. Numerical integration ~ormuI~e The stiffness matrix, and indeed all other element matrices, involve integrals over the area of element
H. V. LAKSHMINARAYANA
992
and K.
KAILASH
Table 2. Gauss point locations and the weighting factors (a) 3 x 3 Network of Gauss points t’l e Node points x Gauss points
Gauss point, g
CR
%
W,
JO.6
1 2
-
3 4 5 6 7 8
,‘,0.6 - JO.6 0.0 ,:‘M -JO.6 0.0
SXS
0.0
2.0 x s x (1.0 -S) -v’O.6 0.0 0.0 0.0
SXS 2.0 x s
X(1.0-S)
4.0 x (I .o - Sl’ 2.0 x s x (I .o -
S)
SXS 3.0 x s
x (I.0
-S)
(b) 2 x 2 Network of Gauss points
I
i
I
/
I
x Gausspoints
Gauss point, g 1 2 3 4
mid-surface,
which are generally
Three (two for invoking the reduced integration technique) Gauss points are required in the [, g directions for parabolic shape functions used in the present etement. The general formula for numerical integration is
of the form
where J* = l&l
- C”)’ :,
(42)
Numerical integration within the appropriate - 1, + 1 limits is carried out using Gauss quadrature.
where W’, is a weighing factor at integration point g whose co-ordinates are (,, ‘la and J: is the value of J* at the integration point. The location of integration points and the weighting factors for 3 x 3 and 2 x 2 networks of Gauss points are given in Table 2.
A shear deformable curved shell element of quadrilateral shape
993
Elastic stiffness matrix
Resolving membrane and shear locking phenomena
The expression for evaluation of the element stiffness matrix [K] is obtained from the strain energy U@)of the element.
To resolve the locking problems, the assumed strain method [2] is used. This approach is felt to be quite general and the one currently gaining favour [I]. To alleviate shear and membrane locking we interpolate the transverse shear strains and the membrane strains independently with appropriately chosen interpolation functions-and tie the coefficients in these interpolations to the strain components evaluated directly by differentiating the displacement fteld at certain points within the element. While studying various possible assumed strain interpolation schemes analytically and numerically it was found that the scheme described below offers much promise. To eliminate the shear locking problem, the procedure of inte~lating transverse shear strain components yEr and ysi, based on natural co-ordinates, and transforming to components based on Cartesian co-ordinates yij, ypi, has been found to play a key role in ensuring improved performance when the mesh is distorted or curved. For the quadratic isoparametric element, the transverse shear strains ySr, yEIshould be assumed as polynomials of at least the following form [6]:
1 IJ”’ = - areawmu~m~ 2s
+ wYmw
+ vYPlwI
+ [XblTPl(Xb~
+ WITWlW)>dfaW,
(44)
where the integration is carried out over the element mid-surface area. Using eqns (34)-(36), the expression for strain energy of the element reduces to the standard form:
with
+
[BB]r[B][BE] + [BB]qD][BB]
+ [BS]r[E] [KS]) d(area).
(46)
When the numerical integration formuia is applied to eqn (46), the equation for the element stiffness matrix computation becomes
+ fm?]f[B][BE]+ [BB]qD][BB]
It is implied that the element strain matrices have to be evaluated at each Gauss point and the elasticity matrices [A], [B], [D] and [El have to be specified with reference to the local Cartesian co-ordinates at each Gauss point. Use of reduced integration (2 x 2 network of Gauss points) in the evaluation of eqn (47) generally results in a more flexible element which has superior convergence characteristics over the standard approach of using full integration (3 x 3 network of Gauss points). The present fo~ulation with full integration exhibits both shear locking and membrane locking phenomena. Use of reduced integration is thus once again imperative for improved performance, and the resulting element is called QUAD8*. Numerical experiments with QUADS* on a set of standard test problems proposed by MacNeal and Harder [4] revealed persistence of locking problems. The task now before us is to modify this formulation to resolve the locking phenomena. This is not easy, particularly when combined with the goal of not permitting any spurious zero energy modes in the element. CA.S
33 CF
Y,,=C,+C*~+~,rl+C,5~~C,5’.
(48)
Consequently, five sampling points are necessary for construction of the assumed transverse shear strain fields. The location of sampling points for ycr and Y,,~are different, and are shown in Fig. 2. The transverse shear strain field in the element is therefore expressible in the form YET = R,t&O&$ + TRY&&
+ R,(~v)Y& +
Wsh:,I;
Y?Jr = &(Slr)Y?, I$ + %@&,,I~ +
+ R&v)Y& (49) + ~~(~~)~~~I~
s,(~?)r,,l~ + m)Y,,I&
(501
where Ri(& 7) and S,(<, q) are the interpolation functions given in Table 3 and y:,(g, y,, jg, . _. , are the strain components at sampling points (I), (2), etc., directly evaluated from displacement interpolations in the usual manner. In the sequel, the following matrix is used to transform the local Cartesian coordinate transverse shear strains 722 and ypz to the natural co-ordinate strains ycr and yrn:
The use of a consistent sequence of inte~olations
and
H, V. LAKSHMINARAYANA and K. KAILASH
994
9 (-l,l) 0
2
1
7 I
T
I I \
0
0 (-1,-l 1
7
I
3x
I 5 -c-)
l&Iv3
i-1,1)
(l,l)
0
I I I I I I I I II I
4 l&Et
(1,-l)
Y m ~
0 (1,-l)
:~4_-_-_______-+:,
I_,32 t,
x3
c
5 _---_-___-____-~:
0 (-1 ,-I 1
(a) For r[r
0
3 0 (1,-l)
(b) For yn7
Fig. 2. Location of points used for interpolation of transverse shear strains.
transformations from the natural co-ordinate system to the local Cartesian co-ordinate system allows the element to be accurate even in the shape of an arbitrary quadrilateral. To resolve the membrane locking problem, the membrane strain components trx, err and cxr have to be interpolated independently such that exx is linear in x and quadratic in y, whereas trr is linear in 7 and quadratic in 8 and txP is quadratic in both f and 7. Therefore, for cxx, err and 6Kpthe interpolation is as follows:
trated in Figure 3; and the superscript DI signifies that these strains are directly evaluated from the displacement interpolations in the usual way. In the computation, a transformation of strains between the local Cartesian co-ordinate system (x, 7, .!?) and an orthogonal curvilinear co-ordinate system (7, S, t) is established and used. The origin of the (7, S, t) system is located at the centre of the element. The axis t is oriented normal to the middle surface and the (y, S, 0) surface coincides with the (5, q, 0) surface. Symmetry between the natural co-ordinate system Table 3. Interpolation
functions for
tranverse shear strains
c*p = +ys + g;.
(52)
In eqn (52), Ri(<, q) and Si(r, q) are the interpolation functions listed in Table 4; subscripts (1), (2), (3), etc., refer to sampling points whose locations are illus-
A shear deformable curved shell element of quadrilateral shape (5, q, T) and the orthogonal curvilinear co-ordinate system (y, S, t) is maintained to ensure that the (y, S, t) system is uniquely defined. In practice, programming and computer implementation of the assumed strain method described above simply involves replacing the original strain matrices - [BE] and [BS] by the modified matrices [BE] and [BS ] everywhere. This will also enable us to employ full integration (3 x 3 network of Gauss points) in the evaluation of the element stiffness matrix, thereby avoiding the occurrence of spurious zero energy modes. The new element is designated [QUADS**].
995
Table 4. Interpolation functions for membrane strains
Pressure loads
The work equivalent nodal forces in the X, Y, Z directions resulting from normal pressure P acting in the z direction over the element mid-surface is
where the subscripts g and i refer to Gauss point and node number respectively. We note that the nodal point moments are null. The normal pressure P may be uniform or may vary over the element mid-surface. In the latter case, the normal pressure at the Gauss point P, is interpolated from corner node point values specified by the user as
a = l/J3.
where h,(<, q) is the bilinear interpolation
function (55)
with
Stress recovery
Element stresses will be recovered at the Gauss points. Stress resultants will be computed in the local Cartesian co-ordinate system in planes parallel to the mid-surface. The stress resultants are calculated using eqns (40). At each Gauss point, global co-ordinates,
1)
(--i.1) 0
1
2
T
I‘ I / I I I 1
(1,l) 0
i 1
*
2 I I I
I I I
I
I ;
6 0
(-1,-l)
(a)
ForezR
t
0
1143 If73
c-)and;
E~,F
(1,-l)
(b)ForepT
and;
~27
Fig. 3. Location of points used for interpolation of membrane strains.
H. V. LAKSHMINARAYANA and K.
996
KAILASH
LBEBQUDJLAMSTF Ipiiq 1 SHLOTR
1
1 BEQUDI
(
BEQUAD BSQUAD BBQUAD SHLGEM
y--JSKOUADJ i JACOB 3 1
s is
LOAD
0
SHLOTR JACOB3 SFRTN GAUSSQ
1 BBQUAD
1
) BQQUAD ]
Fig. 4. Program organization.
directions of the local Cartesian co-ordinate axes and shell thickness are also printed out along with the values Nxx. NFY. NXP, MXX, MPP, MXY, QXZ, QPZ, in the same order. Such components are indeed directly of interest, and by using them we can compute shell surface stress o~,~. flpr, T,?~, 5xZ. ~~2. We note that the shear stress components 7.pzand tp2 are in fact zero at the shell surfaces. However, the values directly obtained for these components are the average values across the shell wall thickness. The maximum value occurs on the mid-surface and is equal to 1.5 times the average value. If nodal point stresses are to be reported, there is a need for a rational and consistent procedure to ensure that the stresses at common nodes of adjoining elements are computed in the same direction before averaging.
of
COMPUTER
IMPLEMENTATION
We have discussed the element formulation and derived the specific expressions needed in the calculation of element matrices. An important advantage of our formulation is the simplicity of programming. In essence, an assemblage of a number of subroutines computes all the relevant matrices for the QUAD8** and QUADS* elements. The specific subroutines are identified in Fig. 4 and their functions are briefly described here. SKQUAD BQQUAD
STQUAD
stiffness Calculation of element matrix Calculation of kinematically consistent load vector due to pressure and volume loading Calculation of element stress resultants
BSBQUD BEBQUD BEQUDI BSQUD 1 SFRTN6 GAUSS6 LMAT GAUSS5 SFRTNS
Calculation of elasticity matrices [A], [Bl. IDI. El Calculation of element strain matrices (BE]. [BB]. and [BS] Calculation of nodal normal vectors V;(i=1,8) Calculation of the transformation matrix [Egl Calculation of Jacobian matrices [J] and [J”] Calculation of shape functions N, (i = 1, 8) and their derivatives Setting up of positions and weights of Gauss points
Calculationof modified matrices [BE] and [BS].
strain
In Fig. 4, the data input subroutine is named INPUT. The subroutine which accepts all load data is called LOAD. Assembly of the element matrices and solution of the resulting system of algebraic equations is carried out in SOLVE. NUMERICAL
EVALUATION
Test problems
MacNeal and Harder [4] have proposed a standard set of problems to test finite element accuracy. A description of each of these problems in terms of geometry, boundary conditions, applied loads, material properties, and finite element meshes is available in [4]. Theoretical results for the test problems are also tabulated in [4] and these have been used to present the results in a normalized form, i.e. the reported value is the ratio of the value obtained here to the theoretical value. The test results are also reduced to a letter grade, using the following rules: Rule
Grade
2% > error 10% > error > 2% 20% > error > 10% 50% > error > 20% error > 50%
A B C D F
For developers of new finite elements these test problems help to find and correct element errors and weaknesses. Almost every problem is capable of evoking results ranging from excellent to poor. Therefore a low grade in one or more problems should not disqualify the element under test for practical applications.
A shear deformable curved shell element of quadrilateral shape
997
Table 5. Results for straight cantilever beam Normalized tip_ displacement in direction of load _ Tip loading direction
QUAD8 [4]
(a) Rectangular elements Extension In-plane shear Out-of-plane shear Twist (b) Trapezoidal elements Extension In-plane shear Out-of-plane shear Twist (c) Parallelogram elements Extension In-plane shear Out-of-plane shear Twist
QUAD8’
QUAD8**
0.999 0.987 0.991 0.950
0.999 0.987 0.992 0.953
0.998 0.985 0.996 0.944
0.999 0.946 0.998 0.943
0.999 0.967 0.986 0.960
0.998 0.906 0.995 0.944
0.999 0.995 0.985 0.965
0.999 0.966 0.977 0.960
0.998 0.965 0.991 0.891
Table 6. Results for curved beam Normalized tip displacement in direction of load Tip loading direction -
QUAD8 [4]
QUAD8*
QUAD8**
In-plane Out-of-plane
1.007 0.971
0.968 0.940
0.955 0.962
The test results are given in Tables 5-l 1 and summarized in Table 12. The significant thing to note is that QUAD8* fails the rectangular plate tests, indicating shear locking, whereas QUADS** reaches
a B grade (an exception is a rectangular plate with clamped boundaries subjected to a point load at its centre which has a C grade for uniform mesh; however, a graded mesh is found to improve the
Table 7. Result for twisted beam Normalized tip displacement in direction of load Tip loading direction
QUAD8 [4]
QUAD8’
QUAD8** -
In-plane Out-of-plane
0.998 0.998
0.993 0.985
0.998 0.998
Table 8. Results for rectangular plate-simple supports; uniform load Normalized lateral deflection at centre No. of node spaces per edge of model
QUAD8 [4]
QUAD8’
QUAD8**
1.OOo
fails
1.016
1.OOo
fails
1.010
(a) Aspect ratio = 1.0 (b) Aspect rat”0 = 5.0 8
Table 9. Results for rectangular plate-clamped supports; concentrated load Normalized lateral deflection at centre No. of node spaces per edge of model
QUAD8 141
OUADB’
0.997
fails
1.090
0.975
fails
0.867
OUAD8**
(a) Aspect ratio = 1.0 (b) Aspect raze = 5.0 8
H. V. LAKSHMINARAYANA and K. KAILASH
998
Table 10. Results for Scordelis-Lo
roof
Normalized vertical deflection at mid-point of free edge No. of node spaces per edge of model 4 8 12
QUAD8 [4]
QUADS*
QUAD8**
0.984 0.997
0.9912 1.006 1.ooo
0.8739 0.9762 0.9787
Table 11. Results for hemispherical shell Normalized radial deflection at load point No. of node spaces
per edge of model
QUAD8 141
4 8 12 16
OUAD8*
OUAD8**
0.178 0.805 0.975 0.998
0.034 0.383 0.765 0.940
0.121 0.823 0.992
grade). QUAD** passes the patch tests, and receives good grades on all the tests, including the spherical shell problem for which QUAD8 achieves a C grade.
Some examples
Of the many problems analysed with the QUAD8** element, four examples are presented here to demonstrate the performance of the QUAD8** element in the analysis of complex problems. Example 1: Cook’s membrane problem. This problem, whose details are given in Fig. 5, is used to evaluate the membrane capabilities of a shell element in the shape of arbitrary quadrilaterals. The vertical defection at point C is compared with the converged finite element solution reported by Cook [lo].
._~
QUAD8** converges quickly and gives an accuracy of 0.6% even with a coarse mesh (4 x 4. as shown in Fig. 5). This example demonstrates the insensitivity of the element to mesh distortion. Example 2: Morley’s skew plate problem. This problem, shown in Fig. 6, is the most severe test for the bending capabilities of a shell element. It is to be noted that the finite element mesh features distorted elements with edges not aligned with the global axes. The convergence of central deflection is shown in Fig. 6. Although a rather large number of elements is required with the chosen mesh, improved results can be obtained by using a different mesh. Example 3: cantilevered curved beam. This problem, illustrated in Fig. 7. has been proposed by
Table 12. Summary of test results for shell elements Element loading
Element shape?
QUAD4 [41
QUAD8 [41
X
Irregular
A
C
Pass
PdSS
X
Irregular All
A A
D A
Pass A
Pass A
X
Regular
B
A
A
A
X
Irregular
F
B
B
B
X
Regular
A
A
A
A
X
Irregular
B
A
A
A
All
B
B
B
B
X
X X X
Regular Regular Regular Regular
C B A B
A B A B
B B A F
B B A B
X
X
Regular
B
A
A
A
X
X
Regular
A
C
B
B
In-plane 1. Patch test
2. Patch test 3. Straight beam, extension 4. Straight beam, bending 5. Straight beam, bending 6. Straight beam, bending 7. Straight beam, bending 8. Straight beam, 9. 10. 11. 12. 13. 14.
twist Curved beam Curved beam Twisted beam Rectangular plate (N = 4) Scordelis-Lo roof (N = 4) Spherical shell (N = 8)
Out-of-plane
X
t Regular means that element shape has not been intentionally distorted.
QUAD8*
QUAD8**
A
999
shear deformable curved shell element of quadrilateral shape R = 1.6m h =0.5m
r =O.Olm P=l.O
h (thickness) = 1 .O Target: V, = 23.91
Fig. 5. Tapered skew cantilever plate under uniform shear at the free end. Fig. 7. Cantilevered curved beam. Robinson [12] for the evaluation of shell elements. The finite element model used divides the beam into four equal-sized elements circumferentially. Predicted free-end displacement in the x-direction is 3.0045 x 10m4 (3.0733 x 10m4), in the y-direction is 1.9133 x 10m4(1.9565 x 10e4) and rotation about the z-axis is 2.5467 x 10m4 (2.6087 x 10m4), where the
numbers in the brackets are the theoretical results. Figure 8 shows that predicted bending moments are in good agreement with theory. Example 4: spherical shell under internal pressure.
The analysis of a thin-walled spherical shell under internal pressure, using the type of mesh shown in
4 X 4 mesh Simply
supported
E=30
x 106
Boundary condition on four edges
edges
u = 0.3 b=
w~
1.0
Thickness Uniform
01 5
h = 0.01
(EXACT)
D = E h3/12
W = 0
=408X .
1O-4
Xpb4 D
(1 -uz)
Pressure = P
I
I
I
I
10
15
I
20
25
30
No. of nodes per edge of model
Fig. 6. Morley skew plate under uniform pressure load.
H. V. LAKSHMINARAYANA and K. KAILASH
1000
1.8 +iii 5 '.62
1.4-
g
1.2-
is l.Oi! 5;
0.8:Zi & 0.6 =
0.4-
10
0
1
I
I
I
I
I
I
20
30
I
40
60
60
70
80
QO
Angle 0 (degrees) Fig. 8. Bending moment distribution for the cantilevered curved beam.
Thickness
t = 0.2
Modulus of elasticity
= E
Poisson’s ratio u = 0.3 Pressure P * 1 .O
Fig.
9. Spherical shell under internal pressure.
Fig. 9, appears
to be one of the most severe (and meaningful) validation tests for a curved sheli element. Because of symmetry, only one octant of the shell extending between the equator and the pole was discretized using QUAD8** elements. The exact solution to the problem is a uniform radial displacement and a pure membrane stress state [I l]. Predicted displacements and stress resultants are given in Table 13. Here NOBand Nb4 denote membrane stress resultants in the hoop and meridional directions respectively. The finite element sofutions are seen to be in excellent agreement with the exact solution. CONCLUSIONS AND DIRECTIONS FOR FURTHER WORK
A shear deformable quadrilateral shell structures
of arbitrary
finite element for
and critically evaluated using a standard set of test problems. The assumed strain method was used in the formulation of the QUADg** element, which is free of membrane and shear locking problems, has no spurious zero energy modes and yields numerical results which are quite insensitive to mesh distortion. Good accuracy and fast convergence toward theoretical solutions were observed for both displaeements and stress resultants in a number of example problems. Future work will address extension of this fo~ulation to thermal stress, buckling and vibration analyses. The present formulation is applicable to homogeneous isotropic as well as laminated anisotropic materials. However, a critical evaluation of the convergence, accuracy and computational etliciency of [QUADS**] when applied to static, dynamic and stability analysis of composite plate and shell structures is essential. Assumed strain fields for curved shell elements must satisfy certain requirements to be convergent [8] and there is a pre-determined variationally correct form that satisfies an exact equivalence of the assumed strain approach to the Hellinger-Reissner or Hu-Washizu mixed variational principles. If the orthogonality condition is also satisfied in designing the assumed strain field interpolations, the resulting element will have faster convergence and will be free from violent stress oscillations. This is a topic for future research. The
shape has been developed
QUAD8**
element
is ideally
suitable
Table 13. Comparison of displacement and stress resultants for spherical shell under internal pressure Stress resultants Normal displacement (~~~/~~Z)
Near pole
At pole
At equator
Exact QUADS*
0.3500 0.3500
0.3500 0.3500
NoelPR 0.500 0.500
QUAD8**
0.360
0.3548
0.5021
N##/pR 0.500 0.500 0.51 t3
Near equator NeelPR 0.500 0.500 0.5027
N,, IPR 0.500 0.500 0.5025
for
A shear deformable curved shell element of quadrilateral shape
implementation in a general purpose structural analysis software. Clearly, the range of application of this type of element is very wide. Applications to practical problems featuring variable and discontinuous curvatures, geometric discontinuities (holes and cracks) and intersections can be confidently undertaken. Acknowledgement-We
thank Mr B. L. Keshavakumar for
help in computations. REFERENCES
T. J. R Hughes and E. Hinton (Editors), Finite Element Methods for Plates and Shell Structures. Pineridge Press, Swansea (1986). T. Belytschko, A review of recent developments in plate and shell elements (197551985). Comput. Mech. Adv. Trends, AMD, ASME 75, 217-231 (1986). R. H. MacNeal, Specifications for the QUAD8 quadrilateral curved shell element. MacNeal-Schwendler, Los Angeles, Memo RHM-46B (1980). R. H. MacNeal and R. L. Harder, A proposed standard set of problems to test finite element accuracy. Finite Elements Anal. Des. 1, 3-20 (1985).
1001
5. K. C. Park and G. M. Stanley, A curved C” shell element based on assumed natural coordinate strains. J. a@. Mech. 53, 278-290 (1986). 6. H. C. Huang and E. Hinton, A new nine node degenerated shell element with enhanced membrane anld shear interpolation. Ini. J. Numer. Merh. Engng 22, 73-92 (1986). 7. K. J. Bathe and E. N. Dvorkin,
A formulation of general shell elements-the use of mixed interpolation of tensorial components. Int. J. Numer. Meth. Engng 22,
697-722 (1986). 8. J. Jang and P. M. Pinsky, An assumed covariant strain based 9-node shell element. Znf. J. Numer. Meth. Engng 24, 2389-2411 (1987). 9. R. H. MacNeal, Derivation of element stiffness matrices by assumed strain distributions. Nucl. Engng Des. 70, 3-12 (1982).
10. R. D. Cook, Improved two-dimensional finite elements. J. Strucr. Div.. AXE 100, 1851-1863 (1976). 11. M. D. Olson, Analysis of arbitrary shells using shallow shell finite elements. In Thin-shell Sfrucfures (Edited bv Y. C. Fune and E. E. Sechler). Prentice-Hall. Englewood Cliffs. NT (1974). 12. J. H. Robinson, New user project: cantilevered curved beam problem for evaluating shell elements. Finife Element News No. 2 (April 1987).
33. No. 4. pp. 987-1001.
1989
Oc45-7949/89 s3.00 + 0 cnl CT‘ 1989 Perpmon Press plc
A SHEAR DEFORMABLE CURVED SHELL ELEMENT QUADRILATERAL SHAPE H. V. LAKSHMINARAYANAand K. KA~LASH Structural Sciences Division, National Aeronautical Laboratory, Bangalore-
OF
017, India
(Received 20 October 1988) Abstract-Formulation and numerical evaluation of QUAD8**, a shear deformable curved shell element of quadrilateral shape, is presented. This element. which has eight nodes and six engineering degrees of freedom per node, is free from locking problems, has no spurious zero energy modes and provides numerical results of acceptable accuracy in a number of test and example problems. The standing of QUAD8** compared with alternatives should consider (1) the simplified input data, (2) the ease of use, (3) improved performance and (4) the simplicity of programming and computer implementation.
INTRODUCTION of simple to use and reliable curved shell elements has been and continues to be an active area of finite element research. Recent advances in this area are contained in a monograph [l] and a review paper [2]. Of particular relavance to the present work is a quadrilateral curved shell element existing in MSC/NASTRAN [3], whose performance was found to be disappointing [4]. This has motivated us to re-examine this element, identify and correct its weaknesses and offer a new element as an alterative. This report presents the formulation and numerical evaluation of a shear deformable curved shell element of quadrilateral shape designated QUAD8**. The assumed strain method is adopted here to alleviate the membrane and shear locking phenomena. Previous contributors to assumed strain curved shell elements include Park and Stanley [S], Huang and Hinton [6], Bathe and Dvorkin [7], Jang and Pinsky [8] and MacNeal [9]. The new element was demonstrated to be reliable and to lead to numerical results of acceptable accuracy in a proposed set of standard test problems [4]. The standing of QUADS**, when compared with alternatives, must take into account (1) the simplified input data, (2) the ease of use, (3) improved performance and (4) simplicity of programming and computer implementation. The report is organized as follows: the geometry of the subject element is defined first, followed by the co-ordinate systems used in the formulation. Succeeding sections describe the Jacobian matrix, displacement interpolations, definition of element strain matrices, constitutive relations, numerical integration formulae, the elastic stiffness matrix calculation, and alleviation of locking phenomena. Numerical solutions to a standard set of test problems and comparison of QUAD8** performance with alternatives are presented next. Of the many problems analysed with this element, four examples The development
987
are presented here: Cook’s plane stress problem, Morley’s skew plate problem, a cantilevered curved beam, and a spherical shell under internal pressure. We conclude by identifying directions for further work. FORMULATION Geometry and coordinate systems
Consider the typical shell element of Fig. l(a). The outer and inner surfaces of the element are curved, whereas the sections across the thickness are generated by straight lines. The geometric description of the element requires the specification of two vectors at each of the eight middle-surface nodes. One is the position vector Ri of the node, i, with the three global Cartesian components Xi, Yi, Z,, where the subscript I’ identifies the node number; the other is the unit normal vector V: along with the shell wall thickness ti of the same nodes. This geometric representation does not ensure continuity of the slope between adjacent elements. Therefore, a shell surface idealized with these elements cannot be exact. However, the accuracy of the approximation increases as the number of elements is increased. The element formulated for the analysis of shells with continuous as well as discontinuous curvatures carries six engineering degress of freedom (three translations and three rotations) at each of its eight mid-surface nodes. The nodal degrees of freedom are illustrated in Fig. l(b). Stiffness for rotation about the normal to the element middle surface is zero at all the nodes. Therefore, a small fictitious rotational spring is introduced at each node, corresponding to the rotation about the normal. Referring to Fig. l(b), let 5, q be the two curvilinear co-ordinates in the element middle surface and t a linear co-ordinate in the thickness direction. If we assumed that 5, q and r vary between - 1 and 1 on the respective faces of the element we can write
H. V. LAKSHMINARAYANA and K. KAILASH
988
(a)
1-4 4
i
7
3
-EL!-
2
8
6
5
1
l-
2
t
2
1
Fig. 1. A curved shell element of quadrilateral shape. (a) Element geometry, (b) coordinate systems and nodal degrees of freedom, (c) parent element and node numbering scheme.
a relationship between the global Cartesian coordinates (X, Y, Z) of any point and the natural co-ordinates (4, I, T) in the form
(1)
Here N,(<, q) is a shape function taking a value of unity at the nodes i and zero at all other nodes. They
are derived as shape functions of a parent element shown in Fig. l(c) and are given in Table 1. In the present formulation, the vector V; is not input by the user but is computed using the co-ordinates of the eight mid-surface nodes as described later. The shell thickness tin eqn (1) is assumed to be constant as its actual value is unimportant for calculating element strain matrices. Actual thickness and its variation within the element is introduced into the constitutive relations.
Table 1. Shape functions and their derivatives Edge nodes (5,7)
Edge nodes (6, 8)
a(1 + &)(I + %)(50+ ‘lo- 1)
f(l - ?)(l
f(1 - a2)(1 + 50)
N:,
i; (1 + %)(250 + %)
-5(l
NV
$(I +50)(2%+&I)
;(l-,F2)
Function
Corner nodes (1, 2, 3, 4)
N,
+ rlo)
+ ‘lo)
;(l -rl(l
-tq) + 5”)
lo = &; ‘lo= qq,; values of 5,. 1,. are defined with reference to the nodes in the diagram given below:
e
1
(-1.-l)
l
5 (O.-l)
2 (1.-l)
989
A shear deformable curved shell element of quadrilateral shape
Let us also define a local Cartesian co-ordinate - system (X, Y, Z) with 2 being normal to the element middle surface and 8, p tangent to it, as shown in Fig. I(b). The orientations of 1, y,z with respect to X, Y, Z are established uniquely using a procedure given in [3]. At any point (5, q) on the element mid-surface (7 = 0) we can define the following two vectors which are tangent to the mid-surface:
.
surface, i.e. IV,>
=
iv,
1 x
(6)
F21.
It is evident that Vi in eqn (1) can be computed using eqn (6). Finally, the local Cartesian to global Cartesian co-ordinate transformation matrix [Eg] is defined as
and This matrix is also used to transform the global displacements (U, V, W) to local displacements --(K K w)
The unit vectors of the local Cartesian co-ordinate system will be determined from the tangent vectors in an ‘unbiased’ fashion, i.e. in such a way that the orientation of the vectors is independent of the sequence of the grid point numbers, except for multiples of 90”. Define:
(8)
Jacobian matrix The transpose of the Jacobian matrix is defined as
-ax ax
ax-
Z&Z dx
[@ c=
axax away --+--+--
(
azaz
x aq at au at al?>i
I$,,
= [J] = 1
(2)
ay ay ay ~~~~ az az az -at; Zy aT_
(9)
By applying the derivatives to eqn (1) it is seen that
where I;, I,, are lengths corresponding to unit increments in 5 and q, and C is the cosine of the angle between the two tangent vectors. The first two unit vectors of the local Cartesian co-ordinate system are
rax
axi
x&w,Iv21=
ay
ay
x
F
44 [ -S/l,
-BP, 41, 1
(3)
az dz _x qi
+
Q and p are chosen to ensure that V, and V, have unit length and are perpendicular. 1 x =2
[
1 (I+(1 1
1 fi=,
[
(1--.
1 +(-)I/2
1
(1 + C)“2
1 1
(5)
We note that if the tangent vectors were perpendicular (i.e. if C = 0) then GL= 1, p = 0 (V, is tangent to the r direction and V, is tangent to the q direction). The third unit vector V, is normal to the middle
7;{3~N~~ N7i OJ).
(10)
In eqn (lo), Nti and N+ are the derivatives of shape functions N,(5, 9) with respect to 5, r~and the vector (ek ejy, e,), is the local normal at one of the eight nodes. In the present formulation, we will need the following two particular forms:
(11)
(12)
H.
990
Combining
V. ~AKSHMINA~AYANA and K. KAILASH
eqns (10)-(12),
matrix:
Definition of membrane strains, bending curvatures and transverse shear strains
Dispiucemmtt iflterpafution
The global displacements (U, V, W) are interpolated from the nodal point displacements and nodal point rotations using the same shape functions IV,{&q) as in eqn (1):
Strains and curvatures are defined in the local --- -. Cartesian co-ordinate system (A’, Y, Z). with A’, Y m the mid-surface and z being normal to it. The normal strain in the 2 direction 62, is neglected. for consistency with the zero normal stress assumption, The membrane strains ate
= -1im r-n
The bending curvatures are where at each mid-surface node i:
The parametric derivatives of the displacements U, V, W are given bY
The transverse shear strains are
au au ?j-zqlt 1 au
du
1
[ &f=
aw -i?w iiw
In the local Cartesian co-ordinate system (2, Y, .?), derivatives of displacements 0, p, k? with respect to X, r, 2 are given by the terms in the following
The evaluation of the membrane strains and the transverse shear strains requires the selection of appropriate terms in eqn (18), evaluated for 5 = 0. However, the evaluation of bending curvatures requires the differentiation with respect to r of terms in eqn (18). Detailed expressions for ail the shell strain components in terms of nodal degrees of freedom may be derived and are given by [3]:
A shear deformable curved shell element of quadrilateral shape
where
[CT”] =B,001 I 0
Bz
Bz
0 [E’]
B,
0
(24)
991
is the nodal dispia~ement vector and [BE], [BB] and [BS] are the so called element strain matrices. They are constructed as
(25)
(26) Element stress resultants and constitutit>e relations (27)
The stress resultants corresponding to the strains and curvatures are defined and related by constitutive relations in the usual manner. The general form of the constitutive relation used in the present formation is
(28)
(29) where (30)
(Nj=
Nxx N,, 1!NXP
(31)
are the membrane
forces per unit length,
(MJ= (32)
Mxx My, ! MKY1
are the bending moments per unit length, and
(33)
In evaIuating the above matrices, t/2 is set equal to unity in the calculation of [q and [,?I because the factor t/2 cancels everywhere. It is to be noted that the membrane strains [cm] do not depend on rotations, which indicates that the present formulation can be used in a shell analysis where bending stiffness is neglected. It is convenient to rewrite eqns (22)-(24) as follows: (34) (35) (36) where
and
are the transverse shear forces per unit length. We note that the element stresses have been defined in the local Cartesian co-ordinate system. Therefore it is also necessary to specify the elasticity matrices [A], [B], [D] and [E] in the same co-ordinate system. A,, B,, D, (i, J = 1, 2, 3) and Ei_,(i, J = 1, 2) are also called membrane, bending, membrane-bending coupling and transverse shear stiffness coefficients. The general form of eqn (40) is included in this element formulation to enable the analyst to model laminated anisotropic materials and sandwich construction in addition to homogeneous isotropic solids. Actual thickness and its variation within the element is accounted for in the calculation of [A], [B], [D], and [E] matrices. It is also a standard practice to include shear correction factors in evaluating the [E] matrix. When thermal loads are considered, eqn (40) has to be modified to include thermal forces and thermal moments. Numerical integration ~ormuI~e The stiffness matrix, and indeed all other element matrices, involve integrals over the area of element
H. V. LAKSHMINARAYANA
992
and K.
KAILASH
Table 2. Gauss point locations and the weighting factors (a) 3 x 3 Network of Gauss points t’l e Node points x Gauss points
Gauss point, g
CR
%
W,
JO.6
1 2
-
3 4 5 6 7 8
,‘,0.6 - JO.6 0.0 ,:‘M -JO.6 0.0
SXS
0.0
2.0 x s x (1.0 -S) -v’O.6 0.0 0.0 0.0
SXS 2.0 x s
X(1.0-S)
4.0 x (I .o - Sl’ 2.0 x s x (I .o -
S)
SXS 3.0 x s
x (I.0
-S)
(b) 2 x 2 Network of Gauss points
I
i
I
/
I
x Gausspoints
Gauss point, g 1 2 3 4
mid-surface,
which are generally
Three (two for invoking the reduced integration technique) Gauss points are required in the [, g directions for parabolic shape functions used in the present etement. The general formula for numerical integration is
of the form
where J* = l&l
- C”)’ :,
(42)
Numerical integration within the appropriate - 1, + 1 limits is carried out using Gauss quadrature.
where W’, is a weighing factor at integration point g whose co-ordinates are (,, ‘la and J: is the value of J* at the integration point. The location of integration points and the weighting factors for 3 x 3 and 2 x 2 networks of Gauss points are given in Table 2.
A shear deformable curved shell element of quadrilateral shape
993
Elastic stiffness matrix
Resolving membrane and shear locking phenomena
The expression for evaluation of the element stiffness matrix [K] is obtained from the strain energy U@)of the element.
To resolve the locking problems, the assumed strain method [2] is used. This approach is felt to be quite general and the one currently gaining favour [I]. To alleviate shear and membrane locking we interpolate the transverse shear strains and the membrane strains independently with appropriately chosen interpolation functions-and tie the coefficients in these interpolations to the strain components evaluated directly by differentiating the displacement fteld at certain points within the element. While studying various possible assumed strain interpolation schemes analytically and numerically it was found that the scheme described below offers much promise. To eliminate the shear locking problem, the procedure of inte~lating transverse shear strain components yEr and ysi, based on natural co-ordinates, and transforming to components based on Cartesian co-ordinates yij, ypi, has been found to play a key role in ensuring improved performance when the mesh is distorted or curved. For the quadratic isoparametric element, the transverse shear strains ySr, yEIshould be assumed as polynomials of at least the following form [6]:
1 IJ”’ = - areawmu~m~ 2s
+ wYmw
+ vYPlwI
+ [XblTPl(Xb~
+ WITWlW)>dfaW,
(44)
where the integration is carried out over the element mid-surface area. Using eqns (34)-(36), the expression for strain energy of the element reduces to the standard form:
with
+
[BB]r[B][BE] + [BB]qD][BB]
+ [BS]r[E] [KS]) d(area).
(46)
When the numerical integration formuia is applied to eqn (46), the equation for the element stiffness matrix computation becomes
+ fm?]f[B][BE]+ [BB]qD][BB]
It is implied that the element strain matrices have to be evaluated at each Gauss point and the elasticity matrices [A], [B], [D] and [El have to be specified with reference to the local Cartesian co-ordinates at each Gauss point. Use of reduced integration (2 x 2 network of Gauss points) in the evaluation of eqn (47) generally results in a more flexible element which has superior convergence characteristics over the standard approach of using full integration (3 x 3 network of Gauss points). The present fo~ulation with full integration exhibits both shear locking and membrane locking phenomena. Use of reduced integration is thus once again imperative for improved performance, and the resulting element is called QUAD8*. Numerical experiments with QUADS* on a set of standard test problems proposed by MacNeal and Harder [4] revealed persistence of locking problems. The task now before us is to modify this formulation to resolve the locking phenomena. This is not easy, particularly when combined with the goal of not permitting any spurious zero energy modes in the element. CA.S
33 CF
Y,,=C,+C*~+~,rl+C,5~~C,5’.
(48)
Consequently, five sampling points are necessary for construction of the assumed transverse shear strain fields. The location of sampling points for ycr and Y,,~are different, and are shown in Fig. 2. The transverse shear strain field in the element is therefore expressible in the form YET = R,t&O&$ + TRY&&
+ R,(~v)Y& +
Wsh:,I;
Y?Jr = &(Slr)Y?, I$ + %@&,,I~ +
+ R&v)Y& (49) + ~~(~~)~~~I~
s,(~?)r,,l~ + m)Y,,I&
(501
where Ri(& 7) and S,(<, q) are the interpolation functions given in Table 3 and y:,(g, y,, jg, . _. , are the strain components at sampling points (I), (2), etc., directly evaluated from displacement interpolations in the usual manner. In the sequel, the following matrix is used to transform the local Cartesian coordinate transverse shear strains 722 and ypz to the natural co-ordinate strains ycr and yrn:
The use of a consistent sequence of inte~olations
and
H, V. LAKSHMINARAYANA and K. KAILASH
994
9 (-l,l) 0
2
1
7 I
T
I I \
0
0 (-1,-l 1
7
I
3x
I 5 -c-)
l&Iv3
i-1,1)
(l,l)
0
I I I I I I I I II I
4 l&Et
(1,-l)
Y m ~
0 (1,-l)
:~4_-_-_______-+:,
I_,32 t,
x3
c
5 _---_-___-____-~:
0 (-1 ,-I 1
(a) For r[r
0
3 0 (1,-l)
(b) For yn7
Fig. 2. Location of points used for interpolation of transverse shear strains.
transformations from the natural co-ordinate system to the local Cartesian co-ordinate system allows the element to be accurate even in the shape of an arbitrary quadrilateral. To resolve the membrane locking problem, the membrane strain components trx, err and cxr have to be interpolated independently such that exx is linear in x and quadratic in y, whereas trr is linear in 7 and quadratic in 8 and txP is quadratic in both f and 7. Therefore, for cxx, err and 6Kpthe interpolation is as follows:
trated in Figure 3; and the superscript DI signifies that these strains are directly evaluated from the displacement interpolations in the usual way. In the computation, a transformation of strains between the local Cartesian co-ordinate system (x, 7, .!?) and an orthogonal curvilinear co-ordinate system (7, S, t) is established and used. The origin of the (7, S, t) system is located at the centre of the element. The axis t is oriented normal to the middle surface and the (y, S, 0) surface coincides with the (5, q, 0) surface. Symmetry between the natural co-ordinate system Table 3. Interpolation
functions for
tranverse shear strains
c*p = +ys + g;.
(52)
In eqn (52), Ri(<, q) and Si(r, q) are the interpolation functions listed in Table 4; subscripts (1), (2), (3), etc., refer to sampling points whose locations are illus-
A shear deformable curved shell element of quadrilateral shape (5, q, T) and the orthogonal curvilinear co-ordinate system (y, S, t) is maintained to ensure that the (y, S, t) system is uniquely defined. In practice, programming and computer implementation of the assumed strain method described above simply involves replacing the original strain matrices - [BE] and [BS] by the modified matrices [BE] and [BS ] everywhere. This will also enable us to employ full integration (3 x 3 network of Gauss points) in the evaluation of the element stiffness matrix, thereby avoiding the occurrence of spurious zero energy modes. The new element is designated [QUADS**].
995
Table 4. Interpolation functions for membrane strains
Pressure loads
The work equivalent nodal forces in the X, Y, Z directions resulting from normal pressure P acting in the z direction over the element mid-surface is
where the subscripts g and i refer to Gauss point and node number respectively. We note that the nodal point moments are null. The normal pressure P may be uniform or may vary over the element mid-surface. In the latter case, the normal pressure at the Gauss point P, is interpolated from corner node point values specified by the user as
a = l/J3.
where h,(<, q) is the bilinear interpolation
function (55)
with
Stress recovery
Element stresses will be recovered at the Gauss points. Stress resultants will be computed in the local Cartesian co-ordinate system in planes parallel to the mid-surface. The stress resultants are calculated using eqns (40). At each Gauss point, global co-ordinates,
1)
(--i.1) 0
1
2
T
I‘ I / I I I 1
(1,l) 0
i 1
*
2 I I I
I I I
I
I ;
6 0
(-1,-l)
(a)
ForezR
t
0
1143 If73
c-)and;
E~,F
(1,-l)
(b)ForepT
and;
~27
Fig. 3. Location of points used for interpolation of membrane strains.
H. V. LAKSHMINARAYANA and K.
996
KAILASH
LBEBQUDJLAMSTF Ipiiq 1 SHLOTR
1
1 BEQUDI
(
BEQUAD BSQUAD BBQUAD SHLGEM
y--JSKOUADJ i JACOB 3 1
s is
LOAD
0
SHLOTR JACOB3 SFRTN GAUSSQ
1 BBQUAD
1
) BQQUAD ]
Fig. 4. Program organization.
directions of the local Cartesian co-ordinate axes and shell thickness are also printed out along with the values Nxx. NFY. NXP, MXX, MPP, MXY, QXZ, QPZ, in the same order. Such components are indeed directly of interest, and by using them we can compute shell surface stress o~,~. flpr, T,?~, 5xZ. ~~2. We note that the shear stress components 7.pzand tp2 are in fact zero at the shell surfaces. However, the values directly obtained for these components are the average values across the shell wall thickness. The maximum value occurs on the mid-surface and is equal to 1.5 times the average value. If nodal point stresses are to be reported, there is a need for a rational and consistent procedure to ensure that the stresses at common nodes of adjoining elements are computed in the same direction before averaging.
of
COMPUTER
IMPLEMENTATION
We have discussed the element formulation and derived the specific expressions needed in the calculation of element matrices. An important advantage of our formulation is the simplicity of programming. In essence, an assemblage of a number of subroutines computes all the relevant matrices for the QUAD8** and QUADS* elements. The specific subroutines are identified in Fig. 4 and their functions are briefly described here. SKQUAD BQQUAD
STQUAD
stiffness Calculation of element matrix Calculation of kinematically consistent load vector due to pressure and volume loading Calculation of element stress resultants
BSBQUD BEBQUD BEQUDI BSQUD 1 SFRTN6 GAUSS6 LMAT GAUSS5 SFRTNS
Calculation of elasticity matrices [A], [Bl. IDI. El Calculation of element strain matrices (BE]. [BB]. and [BS] Calculation of nodal normal vectors V;(i=1,8) Calculation of the transformation matrix [Egl Calculation of Jacobian matrices [J] and [J”] Calculation of shape functions N, (i = 1, 8) and their derivatives Setting up of positions and weights of Gauss points
Calculationof modified matrices [BE] and [BS].
strain
In Fig. 4, the data input subroutine is named INPUT. The subroutine which accepts all load data is called LOAD. Assembly of the element matrices and solution of the resulting system of algebraic equations is carried out in SOLVE. NUMERICAL
EVALUATION
Test problems
MacNeal and Harder [4] have proposed a standard set of problems to test finite element accuracy. A description of each of these problems in terms of geometry, boundary conditions, applied loads, material properties, and finite element meshes is available in [4]. Theoretical results for the test problems are also tabulated in [4] and these have been used to present the results in a normalized form, i.e. the reported value is the ratio of the value obtained here to the theoretical value. The test results are also reduced to a letter grade, using the following rules: Rule
Grade
2% > error 10% > error > 2% 20% > error > 10% 50% > error > 20% error > 50%
A B C D F
For developers of new finite elements these test problems help to find and correct element errors and weaknesses. Almost every problem is capable of evoking results ranging from excellent to poor. Therefore a low grade in one or more problems should not disqualify the element under test for practical applications.
A shear deformable curved shell element of quadrilateral shape
997
Table 5. Results for straight cantilever beam Normalized tip_ displacement in direction of load _ Tip loading direction
QUAD8 [4]
(a) Rectangular elements Extension In-plane shear Out-of-plane shear Twist (b) Trapezoidal elements Extension In-plane shear Out-of-plane shear Twist (c) Parallelogram elements Extension In-plane shear Out-of-plane shear Twist
QUAD8’
QUAD8**
0.999 0.987 0.991 0.950
0.999 0.987 0.992 0.953
0.998 0.985 0.996 0.944
0.999 0.946 0.998 0.943
0.999 0.967 0.986 0.960
0.998 0.906 0.995 0.944
0.999 0.995 0.985 0.965
0.999 0.966 0.977 0.960
0.998 0.965 0.991 0.891
Table 6. Results for curved beam Normalized tip displacement in direction of load Tip loading direction -
QUAD8 [4]
QUAD8*
QUAD8**
In-plane Out-of-plane
1.007 0.971
0.968 0.940
0.955 0.962
The test results are given in Tables 5-l 1 and summarized in Table 12. The significant thing to note is that QUAD8* fails the rectangular plate tests, indicating shear locking, whereas QUADS** reaches
a B grade (an exception is a rectangular plate with clamped boundaries subjected to a point load at its centre which has a C grade for uniform mesh; however, a graded mesh is found to improve the
Table 7. Result for twisted beam Normalized tip displacement in direction of load Tip loading direction
QUAD8 [4]
QUAD8’
QUAD8** -
In-plane Out-of-plane
0.998 0.998
0.993 0.985
0.998 0.998
Table 8. Results for rectangular plate-simple supports; uniform load Normalized lateral deflection at centre No. of node spaces per edge of model
QUAD8 [4]
QUAD8’
QUAD8**
1.OOo
fails
1.016
1.OOo
fails
1.010
(a) Aspect ratio = 1.0 (b) Aspect rat”0 = 5.0 8
Table 9. Results for rectangular plate-clamped supports; concentrated load Normalized lateral deflection at centre No. of node spaces per edge of model
QUAD8 141
OUADB’
0.997
fails
1.090
0.975
fails
0.867
OUAD8**
(a) Aspect ratio = 1.0 (b) Aspect raze = 5.0 8
H. V. LAKSHMINARAYANA and K. KAILASH
998
Table 10. Results for Scordelis-Lo
roof
Normalized vertical deflection at mid-point of free edge No. of node spaces per edge of model 4 8 12
QUAD8 [4]
QUADS*
QUAD8**
0.984 0.997
0.9912 1.006 1.ooo
0.8739 0.9762 0.9787
Table 11. Results for hemispherical shell Normalized radial deflection at load point No. of node spaces
per edge of model
QUAD8 141
4 8 12 16
OUAD8*
OUAD8**
0.178 0.805 0.975 0.998
0.034 0.383 0.765 0.940
0.121 0.823 0.992
grade). QUAD** passes the patch tests, and receives good grades on all the tests, including the spherical shell problem for which QUAD8 achieves a C grade.
Some examples
Of the many problems analysed with the QUAD8** element, four examples are presented here to demonstrate the performance of the QUAD8** element in the analysis of complex problems. Example 1: Cook’s membrane problem. This problem, whose details are given in Fig. 5, is used to evaluate the membrane capabilities of a shell element in the shape of arbitrary quadrilaterals. The vertical defection at point C is compared with the converged finite element solution reported by Cook [lo].
._~
QUAD8** converges quickly and gives an accuracy of 0.6% even with a coarse mesh (4 x 4. as shown in Fig. 5). This example demonstrates the insensitivity of the element to mesh distortion. Example 2: Morley’s skew plate problem. This problem, shown in Fig. 6, is the most severe test for the bending capabilities of a shell element. It is to be noted that the finite element mesh features distorted elements with edges not aligned with the global axes. The convergence of central deflection is shown in Fig. 6. Although a rather large number of elements is required with the chosen mesh, improved results can be obtained by using a different mesh. Example 3: cantilevered curved beam. This problem, illustrated in Fig. 7. has been proposed by
Table 12. Summary of test results for shell elements Element loading
Element shape?
QUAD4 [41
QUAD8 [41
X
Irregular
A
C
Pass
PdSS
X
Irregular All
A A
D A
Pass A
Pass A
X
Regular
B
A
A
A
X
Irregular
F
B
B
B
X
Regular
A
A
A
A
X
Irregular
B
A
A
A
All
B
B
B
B
X
X X X
Regular Regular Regular Regular
C B A B
A B A B
B B A F
B B A B
X
X
Regular
B
A
A
A
X
X
Regular
A
C
B
B
In-plane 1. Patch test
2. Patch test 3. Straight beam, extension 4. Straight beam, bending 5. Straight beam, bending 6. Straight beam, bending 7. Straight beam, bending 8. Straight beam, 9. 10. 11. 12. 13. 14.
twist Curved beam Curved beam Twisted beam Rectangular plate (N = 4) Scordelis-Lo roof (N = 4) Spherical shell (N = 8)
Out-of-plane
X
t Regular means that element shape has not been intentionally distorted.
QUAD8*
QUAD8**
A
999
shear deformable curved shell element of quadrilateral shape R = 1.6m h =0.5m
r =O.Olm P=l.O
h (thickness) = 1 .O Target: V, = 23.91
Fig. 5. Tapered skew cantilever plate under uniform shear at the free end. Fig. 7. Cantilevered curved beam. Robinson [12] for the evaluation of shell elements. The finite element model used divides the beam into four equal-sized elements circumferentially. Predicted free-end displacement in the x-direction is 3.0045 x 10m4 (3.0733 x 10m4), in the y-direction is 1.9133 x 10m4(1.9565 x 10e4) and rotation about the z-axis is 2.5467 x 10m4 (2.6087 x 10m4), where the
numbers in the brackets are the theoretical results. Figure 8 shows that predicted bending moments are in good agreement with theory. Example 4: spherical shell under internal pressure.
The analysis of a thin-walled spherical shell under internal pressure, using the type of mesh shown in
4 X 4 mesh Simply
supported
E=30
x 106
Boundary condition on four edges
edges
u = 0.3 b=
w~
1.0
Thickness Uniform
01 5
h = 0.01
(EXACT)
D = E h3/12
W = 0
=408X .
1O-4
Xpb4 D
(1 -uz)
Pressure = P
I
I
I
I
10
15
I
20
25
30
No. of nodes per edge of model
Fig. 6. Morley skew plate under uniform pressure load.
H. V. LAKSHMINARAYANA and K. KAILASH
1000
1.8 +iii 5 '.62
1.4-
g
1.2-
is l.Oi! 5;
0.8:Zi & 0.6 =
0.4-
10
0
1
I
I
I
I
I
I
20
30
I
40
60
60
70
80
QO
Angle 0 (degrees) Fig. 8. Bending moment distribution for the cantilevered curved beam.
Thickness
t = 0.2
Modulus of elasticity
= E
Poisson’s ratio u = 0.3 Pressure P * 1 .O
Fig.
9. Spherical shell under internal pressure.
Fig. 9, appears
to be one of the most severe (and meaningful) validation tests for a curved sheli element. Because of symmetry, only one octant of the shell extending between the equator and the pole was discretized using QUAD8** elements. The exact solution to the problem is a uniform radial displacement and a pure membrane stress state [I l]. Predicted displacements and stress resultants are given in Table 13. Here NOBand Nb4 denote membrane stress resultants in the hoop and meridional directions respectively. The finite element sofutions are seen to be in excellent agreement with the exact solution. CONCLUSIONS AND DIRECTIONS FOR FURTHER WORK
A shear deformable quadrilateral shell structures
of arbitrary
finite element for
and critically evaluated using a standard set of test problems. The assumed strain method was used in the formulation of the QUADg** element, which is free of membrane and shear locking problems, has no spurious zero energy modes and yields numerical results which are quite insensitive to mesh distortion. Good accuracy and fast convergence toward theoretical solutions were observed for both displaeements and stress resultants in a number of example problems. Future work will address extension of this fo~ulation to thermal stress, buckling and vibration analyses. The present formulation is applicable to homogeneous isotropic as well as laminated anisotropic materials. However, a critical evaluation of the convergence, accuracy and computational etliciency of [QUADS**] when applied to static, dynamic and stability analysis of composite plate and shell structures is essential. Assumed strain fields for curved shell elements must satisfy certain requirements to be convergent [8] and there is a pre-determined variationally correct form that satisfies an exact equivalence of the assumed strain approach to the Hellinger-Reissner or Hu-Washizu mixed variational principles. If the orthogonality condition is also satisfied in designing the assumed strain field interpolations, the resulting element will have faster convergence and will be free from violent stress oscillations. This is a topic for future research. The
shape has been developed
QUAD8**
element
is ideally
suitable
Table 13. Comparison of displacement and stress resultants for spherical shell under internal pressure Stress resultants Normal displacement (~~~/~~Z)
Near pole
At pole
At equator
Exact QUADS*
0.3500 0.3500
0.3500 0.3500
NoelPR 0.500 0.500
QUAD8**
0.360
0.3548
0.5021
N##/pR 0.500 0.500 0.51 t3
Near equator NeelPR 0.500 0.500 0.5027
N,, IPR 0.500 0.500 0.5025
for
A shear deformable curved shell element of quadrilateral shape
implementation in a general purpose structural analysis software. Clearly, the range of application of this type of element is very wide. Applications to practical problems featuring variable and discontinuous curvatures, geometric discontinuities (holes and cracks) and intersections can be confidently undertaken. Acknowledgement-We
thank Mr B. L. Keshavakumar for
help in computations. REFERENCES
T. J. R Hughes and E. Hinton (Editors), Finite Element Methods for Plates and Shell Structures. Pineridge Press, Swansea (1986). T. Belytschko, A review of recent developments in plate and shell elements (197551985). Comput. Mech. Adv. Trends, AMD, ASME 75, 217-231 (1986). R. H. MacNeal, Specifications for the QUAD8 quadrilateral curved shell element. MacNeal-Schwendler, Los Angeles, Memo RHM-46B (1980). R. H. MacNeal and R. L. Harder, A proposed standard set of problems to test finite element accuracy. Finite Elements Anal. Des. 1, 3-20 (1985).
1001
5. K. C. Park and G. M. Stanley, A curved C” shell element based on assumed natural coordinate strains. J. a@. Mech. 53, 278-290 (1986). 6. H. C. Huang and E. Hinton, A new nine node degenerated shell element with enhanced membrane anld shear interpolation. Ini. J. Numer. Merh. Engng 22, 73-92 (1986). 7. K. J. Bathe and E. N. Dvorkin,
A formulation of general shell elements-the use of mixed interpolation of tensorial components. Int. J. Numer. Meth. Engng 22,
697-722 (1986). 8. J. Jang and P. M. Pinsky, An assumed covariant strain based 9-node shell element. Znf. J. Numer. Meth. Engng 24, 2389-2411 (1987). 9. R. H. MacNeal, Derivation of element stiffness matrices by assumed strain distributions. Nucl. Engng Des. 70, 3-12 (1982).
10. R. D. Cook, Improved two-dimensional finite elements. J. Strucr. Div.. AXE 100, 1851-1863 (1976). 11. M. D. Olson, Analysis of arbitrary shells using shallow shell finite elements. In Thin-shell Sfrucfures (Edited bv Y. C. Fune and E. E. Sechler). Prentice-Hall. Englewood Cliffs. NT (1974). 12. J. H. Robinson, New user project: cantilevered curved beam problem for evaluating shell elements. Finife Element News No. 2 (April 1987).