- Email: [email protected]
High-order triangular finite element for shell analysis
HIGH-ORDER
TRIANGULAR FINITE FOR SHELL ANALYSIS
ELEMENT
D. J. D:\\rr Department
of Civil Engineering.
Ilniversity
of Krmir&m.
Hlrmingh.lnr.
I’r~gI:~ntl
Abstract-A curved-shell finite element of triangular shape i\ d~\cl-ii~itl !t hich 1%h;lwd (II) ~w~vcntional \hell theory expressed in terms of surface coordinate\ and dirpl;~cemcnr\. kwh ~rl’ IhL> ~hrw wt’.lcc’ dispkemrnt components is independently reprewnted I~> ;I twt~-Jill~c~~\i~ln;II pal! II~~~I.II IA ~~~~~~lr.;iiI1~~.J-r~uinli~ illdcr giving the element a total of 54 degrees of frwdom. T\to p;lrlicul,w ~c~~mc~ric l’~wm~ Ilf the clement XC considered, viz. doubly-curved shallow and circular cylindrkll. The h+h lc\cI t)f ~LL’LWI;I<~ 1%hiih can br achieved using few elements is Jemonstrakd in a range of prl,hlcm\ u hw L‘~w~p;lrwm i\ m;~de \vith preview finite element solutions.
I. INTR~)DL’CTION
Recent studies by the author[ I. 21 of the comparative have demonstrated displacement
the characteristics
fields and on assumed-strain
interpolated
displacement
components
efficient in a range of applications nearly-inextensional
behaviour.
eticirn
of v;u-ious arch finite elements
of models based both on rrlativcly
high-order
of quintic order h:~ve been shou,n to be consistently shallow :ind Jeep gccxnctrir\
embracing This efficiency
and
is \uch that in mc&!lling
in the calculated
applications[2]; polynomial
maximum
calculated
displacement
displacement
component
forces and bending moments arc‘ ;_Iso accunlrr.
on problem geometry all applications;
is restricted
force
components
the
Elements based on
i\ rcduccd to cubic order are
as QXII ;I\ on ;m clement basis. In particular,
to ;L cubic v;lriation
and the calculated
and
arches with a
i\ JCYSthan 04’~; in an) of the considered
fields in which one or both
very much less efficient on a degree-of-freedom the tangential
very
extensional
circular
single element based on surface (i.e. tangential and normal) quintic displacement error
polynomial
fields. Arch moc!cI , I~;~\ell on (two) independently-
di~trihutiow;
this waviness becomes very pronounced
tierzen
cfficirncy
I’ot, \tIch ;:iAcl\
if
Jcpcnds markedly exhibit
in nc~lrl!,-in~utensional
waviness in
problems where
its order can be very many times greater than the [rue mugrtitudc of the forcr’. This ia so whether the associated normal component of displacement ih of iubis 01’ quiiitic or&r and, in general. there is no improvement in increasing the order of the :l~!rm;lI ccrmponent I’rom cubic to quintic unless the order of the tangential a significant
improvement
component
incrc;i XXI 1though ~\ccptionally
is \imil;trly
where the arch geometry
it; sh;~llo\,~ .~nd r~lati\.cly
The arch studies point to the use of quintic polynomi:li for doubly-curved
shell elements.
elements of triangular
planform
lllspl;l~~mcnl
In fact there aIrcad> c&t
;r\~un;plic:n k
with quintic polynomial
component5
iii ihe literature I’IW
there is
thick). as the basis
two conforming
AI three displacement
components[3,4] but these components are carlesiun ones. C‘itrtrAn disp!;tcement components were chosen since they allow the precise rrprrscnt;ltion of lhc ri$-bcA!~ motions of an element ibf the when the geometry of the undeformccl shell is dctinccl ill linc,ir cc~mllination\ 111 :hc‘ norh of Argyris same set of basic functions used to define the ~li\pl;~ccmrnt~. and Scharpf[3] the development of the shell theory WI \\hich the c’lemcnt is based uses the natural strain concept and its relationship of Dupuis and Goel[4]
to clahsic;tl thin \hcll Ihcory iq not obvious.
is based on shell theory \\rittcn
in term\ of :I
The work
cuordinatr:
system
in place of the usual curvilinear system and the equation\ ;IIY crprc\\cd in relation to the height of the curved middle surface above a reference plane. Po>\iblc difficultlcs ;Ih\ociatrd with such schemes as these have been pointed out by Morris!!]. In ;I further stud! Dupuis[h] has described a triangular shell element based on a modified repre\ent;ltion 01 the bhcll geometry in which the element is straight-sided in the plane of shell parameter\. Kigid-hod! motions are represented precisely only by making an appropriate approximation IO the shell middle surface and although the derivation is, in principle, quite general it is spe&list’d to the LYISC\vhcrc the displacement components
are represented
by third-order
polynomials.
motion
states
(whose
inextensional
mode\.
states can be closely shell theory
nature
i\ dihcusacd
ik not neccwry approximated.
and the ahsumption
This C;II~ hc achieved
polynomi;rl wrface of “sensitive” easily
allows
thr: exilct
paper det;lil\
element
which helong\
surface
component
polynomial here
the
~huracteristic~ radius of twkt
limited
twent!
to t\vo
iii ~l~wd
to \tuJ!,
and bending
The :Ipproach
of st;md;trlI
IYIWXI on it)tl\~~litl~~l\:l;
ti.c.
-one term\
~~‘cwlpon~lll~~
classes
of \hcIi.
h;i\iciill\
.hc!i
\hcll
i L~,~k.!!
iI ~oiiiplctc ionstr;lint<
The philowph!
of
rri;lnylt,l!
par~;~gr;~pl:. l‘h~~ or&l
to M hich thrw
of freedom.
specific
\haIw
of :I ;onlornling
in the preceding
of
Hith
form
;imi numeric;il
with :I nltrncric;ll
in detail the preciw
momc‘nt~
de~crihed
and Olson19-121 assumption
,I\ -IIL’~I
i\ \utticientlv
i4 con~tr;rinecl-quintic having
small error a\wciated
displacement
wch
q:ri!it~~
;lrc ;ipIliiclI i
i\ iI general
gcomctr\
\\ hich
OIW IV!
exhibit
1111:
of general \hellh. \ iz. the shalltw shell with twj principal ratlii of cur\;ltuir’ I~III~. .: antI the circular c.ylindric;li khcll. ‘l‘hi\ rc\triction !\ ~mpcwd MI that the cIi-ritc’l>:
Gnce it is proposed forces
\tatc\
polynt~mi~~l surfaa
and application
ha\ ;L total of 54 degree\ i\
may he cibtaincd
any possible
in ;m clcmcnt
rrprewitation
IO the r~atcgor~ dcscrilwl polynomi;il
analyG3
3titfnesF
geometric
the tlewlopmun~
in two dimension\
and thus the element
components
h;l\ic
high. I I I> !‘llrthr*r notl:cI ix:’ l.lll! displncement fields arc p;irticuMy appropriiltc iii dcaliii g \%lth ;!I1 11111’1l~1 prohlcm~ in shell worklS.Sj antI thxt ~111anal!\~ \ 11;,,1’,1 ,111 i’,lll! ~‘!!!lrlll.~i
claw
The prewnt
of other
and rapid C~U~VCI~C’IICC’pro\‘iditg
of independently-intcrpol;lted
\o long a\ the order oi the polynomial
theory
in 171) or indeed
for high accuracy
hcrc ik tlircctlh
.;\ 0iAQ
schcmt.
of c:llcul;M
nwn:il
I’trta ~limiii;ll!tibl:
I\ ~~~~I~~,IcI~cJ JC,II .lhil
1114t’la~‘c’r,l~iit\. mL*mhi Lint’
stlrf;lcca
~~*‘/atcd to that xdnptcd
cliffawcc
:I con~trainril-qtiiiitic
components
wriation
over the shell middle
the fund:unentnl
intcgtxticiri
intqqticrn
that whilst
I’:II.~IC~ IT\ (‘,I\\
the txlrlicr
dicpl;icemtnt
\\trrk
component.
arc taken to hr cubic pol~nomial4.
plea’, I IINII~~I~~
i\ h;isctl 1111tht,
the
I\\ 11 t:lngcnti;ll
Lsfd. 1I?! arntciid
t’~~.pcr.
th;11. III!
;I wt numhcr of clement\. inzrca4ng the 11rdc1 of the t;lngcnti;ll ~‘~~iiipoii~iit~ IO L‘lw4!-:lin~~i quintic will lead to only :I mwginal impro\~ment in ;lccur;lc! achit~\~~~tl ;II the C.XIWIIW ,li .)I. increase ;i
in the degree< of frcrdom
T;iylor’~
limiting
of CW’I, II%\ contention
series hut i\ certainly
not correct
C~IW of the dtxp circulw
etticient
i\ has.~l on error ~l\nsider;ttic~n~
:ii 211 \ituaticjn\
arch ;I motlcl with quintic
on :I de~ree-rlf-f’r~~c(l~~nl
h;l\i\
th:ln
\in
:I\ m~nliciiicd
tangenti;li
J awrc\pondirig
IIVII;.
c;lrlic:.
cli~placc‘mcnt
II’ lhc
i\ nlu~b WOI~
I~~IJc:! u iti)
;;rnecnti;ii
displacement. Any its\essment
of the merit\
must take account geometry
and loading hehaviour.
hending.
etc.: rigid-body
structure
with
detailed
error
arch [ 1X-151).
or the shell might
prediction large.
Accordingly.
numerical
application earlier
which
include\
ma!’
function
clement\.
ha\c
to \peciIic
nil the c;ltc~oriec
hc tlominatctl primarily
;IIMI\
difficult
prewnthd
for
h!
plohlem~.
Hwcver.
prc\cntcd
noted ;Iho\ C. In the prcscntation
\\ ith Ioc:li /I~~!~‘\ 1): circurn\t;~ll~c.~
I\ noted
prohlcm
of
;I LI iLk
.!
wpplcment~cl I’:III~C
it’
th~~ilgh I~;!I thv
her‘. i\ cwmincd
wch ;lpplic;ltiotl\.
cast* of the ;Irch 1I. 21.rmhracc
:!!’ !j!
of ;I general illr: c:i ,hcll III
\implitied
II~)!.
hehaviorr-.
In thcw
modellin~
\IIIIL’III~~~-
I)ependill;
memhranc
to awtruct.
the
of the clement
XI\ ot \hell
;I\ ;I mcmhr;Inc
of :I linite clrmcnt i\ \~r\’
heen
the ilccur;lc!
for the limiting
for the
, of the ~lcmcnl
may or rn;~! not he of importnnc‘c.
of the ilccurilc\
unequal
estimilte\
presented
motions
propowd
of p~~~sihlc application
the shell rehponw
hending
v;llid theoretical
of ;I finite element
of the di\ercit!.
4lcII
..ll.i‘ll!.ll
,~~lci~~ !:I lb: ~.,*\uItll,:j!.:.
$11 ,lct~~ilotl iiumCri~;::
illlll
.:, .irlt~
p:~rticular attenticjn i\ paid tcr ;I
1L
lesser improvement in \hallo\r :lpplicaticrll~ \t;Igc one of the wfercl*~ ,)f this papc~’ h.1..
;I
1Ih] in cl:hich ;I farnil!, of \hcll
\I
)‘SlS
n1:1.
\lI
c!cmcnts
:.IILUII
rt4;ltcl.: ;I~ thut
s
the many shell theorie\ availahle in the literature that of Koitcrl 171 is :Idopted hcrc. In thih theory (unlike wnc others I the homogenecw\ strain-displacement equation\ ;irc’ \;lti\Ii(.(i in ;I general rigid-hod! moticjn. ‘I hc theor! ;14\unw an orthogonal curvilinwr ic~cwclin;ltc \\ \I;‘ITI !>\I: . .. ii,, thk \!‘\tern ntx~l I~OI cc)inciJl: \rith the /!I){: Y ,li. tirincip;ll ~~trr\:~tirr~~ -II !Ib.~ ,hcil bliit‘;i, Of
I099
High-order triangular finite element for shell analysis
notation
of Koiter
the strain energy
the sum of the extensional
density.
and flexural
dV,,, of a shell of thickness
energies
h may be expressed
as
in the form
or
Here
and the (6 x 6) matrix along orthogonal
[B] is easily constructed.
parametric
and T are the physical Also the quantity
where
curves
components
C ic defined
E and u are Young’s
As mentioned
earlier,
types of shell-those evaluated
in closed
conveniently
form.
expressed
strains
these curves:
and twist referred
xl, e
to these curves.
as
modulus
of shallow
strain-displacement
lI and lI are the extensional
of the changes of curvature
and Poisson’s
the present
simplified
The quantities
(Y, p and $ is the shear strain between
paper
geometry
specifically
and those of circular
equations For both
ratio respectively.
is concerned
with
two representative
cylindrical
geometry-so
can be used and thus the element
types
of shell the strain-displacement
that
stiffness
can be
equations
can be
as IE] =
[Wd)
13)
where (4)
For the shallow parallel
shell
to suitable
component.
global
Following
II and
V are the tangential
Cartesian
coordinates
the usual assumptions
ID1=
of shallow 0
0
l/R,
0
0
0
0
I
l/R:
0
0
0
0
I
I
0
2/T
0
0
0
0
I
0
0
0
0
0
I
0
0
I
0
0000
directions
orthogonal
to X:
[D]
and twist of the middle surface
runs parallel
U. V and
Y coordinate
to the longitudinal
W are displacement
and the shell normal respectively.
is
(5)
related to
axis and Y is the components The matrix
in the
[D] in
to be
0
I
0 I 0
I 0 0
l/R 0 0
0 0 I
0 0 0
0
-l/R
0
0
0
I
$R 0
-;R 0
00
0 0
0 0
0I
00I
00
[Dl=
0 0 00
0 -1 where
0
radii of curvature
shell the X coordinate
of the X coordinate,
this case is taken
0
the matrix
in directions displacement
1
0 0
where R,, R, and T are the (constant) the X. Y coordinates. coordinate
shell theory.
0
0000
For the cylindrical
components
-1
L(10
curvilinear
displacement
X and Y. and W is the normal
R is the mean radius of the cylinder.
0 0 a
(6)
I). .t.
I Ino
lhar:
For the tinitc clement ;tnulysis the domain of the shell in the XY plane (representing plane in the case of the shallow shell or the developed
plane for the cylindrical
the base
shell) is divided
into triangular elements :I\ shown in Fig. I. Within an element it is required that each of the three displacement components varies basically as a quintic polynomial in the two coordinates. As is well known. the use of ;I quintic polynomial displacement field originated in the development of refined conforming triangular plate bending elements by a number of independent investigators (see [Is]).
Corresponding
to ;I complete
degrees of freedom comprising together
quintic
the deflection
with the normal deric;ltivc
polynomial
such an element has twenty-one
with its first and second derivatives
at the mid-point
of each side. However,
at each vertex
a more convenient
eighteen dof (constr;linrtl-quinli~) tli+lc.cmcnt field is obtained when the mid-side degrees of freedom are rlimin;~tcJ 11~ re\tricting the \.uriution of normal slope along an edge to a cubic polynomial and it i\ thkb crln~tr;linecl-qtlintic polynomial displacement component. Cunneition of itll nodal quantities boundaries
of each displacement
component
of the normal derivative of the tangential with the potential energ!’ principle.
I’
/
which is selected here for each ensures continuity at inter-element
and of its normal derivative, components
is not strictly
although the continuity
necessary
for compliance
x,u
/
/”
/’
i”’ _.-_.
I
-.-____--..L_
-
Having dccidcd in hro;rd term\ the hiI& element
stiffness c:m he implcmcnted
extend and modify general closed-form
x,u
_.___
in
of the shell element the detailed calculation
:I
number of ways, but it is most convenient
of the here to
the nppro:lch described by Cowper et a/.[9. lo]. This approach leads to a integration formulas for the evaluation of stiffness terms via the introduction
of local coor-din;& ;IYC'Y I :~ntl y (we Fig. I). The corresponding local displacement components arc the tungcnti;d cmcs II ;~d I’ tllP;I\tIrCd in the directionr of s and y respectively, together with the component
norm;ll to illc &cll \urf;lce.
The component\
with like cuprccrions components
of tht.. conforming
IV. which coincides with the global component, displacement
for I‘ ;Ind I(’ giving 60 independent
;lr’e more \uit;lhl!,
W.
field are assumed to be of the form
coefficients
Ai in all. The displacement
c~:presd
where mi . . si are integer< runnln,c hetween t) and 5. The above polynomial
does not contain an
s4~ term so th;lt the norm:~l hlopc \;lriation along edge l-2 is automatically restricted to a cubic polynomial: two further cclnstr;lint conditions per component ;lre applied later to similarly restrict the normal slope variation
;dong the remaining
thus reduced to IX per component.
5-l in dl.
edges. The number of independent
coefficients
is
1101
High-order triangular finite element for shell analysis
Following the example of [lo] the strain energy for the shallow or the cylindrical shell element can then be expressed (by substituting eqn (5) or (6) into eqn (1) via eqn (3) and transforming to local displacements and coordinates) in the form of an integral over local coordinates as
[d][L]{d} dx d?,
(81
where (9) and IL1 = w41T[Dl’vww?41.
(10)
Matrices [I&] and [t] have the same general meaning here as in [lo]; the rotation matrix [I&], which links the column vectors {d} and {d} (i.e. {d} = [R,](d)) is in fact unchanged whilst matrix [L] does, of course, differ in detail both between the shallow and cylindrical shell analyses described here and the earlier analysis. The strain energy can be further expressed via eqn (7) in the form Vo= ;[Al[klIAl.
(11)
Here {A} is the column vector of the sixty coefficient Ai and [k] is a (60 X 60) “stiffness” matrix whose elements kii are explicitly defined in terms of the elements of the (8 x 8) matrix [i] by a modified form of the equation given in the Appendix of [ 101.Since this modified equation is a long one it will not be presented in full here; rather the necessary changes from the original equation will be noted. For the shallow shell element described here the expression for kij is that of the original equation but rows 8 to 13 inclusive as printed may be deleted since the elements of matrix [L] which are involved (L (1,6) through to L (4,B)) are in this case zero. For the cylindrical shell element presented here the expression for kri is the complete one of [IO] with the addition of the two terms tL(6,7)[n(ri
- l)r,q + rj(ri - t)ris]F(ri
t rj -3, si + Sj- I)
(12) + L(7, 8)[risisi(si
- l)+ rjsjsi(si - l)]F(ri + rj - I, si + SI- 3)
where F(m, n) is defined in [lo]. The “stiffness” matrix [k] corresponds to the quadratic form (eqn 11) of the strain energy in terms of the polynomial coefficients {A}. The energy can be further expressed in a quadratic form of the global degrees of freedom as v, = T”[ WZIP-GI~ W:I
(13)
[&I = [RlT[~~l’[kl~~4Rl
(14)
where,
is the required global stiffness matrix. Here the column vector of global degrees of freedom {Wz} of length 54 comprises the quantities au au a2u
u --
9ax’ay’ax”axau’p9
a2u
a2u
v
“..
a2v -@w**--ayz
a2w
IlO2
I). .I, IhW’F
at each of nodes I. 2 and 3 in turn. The matrices I I’, 1and I R I here. of course, take a different form
from that used in the work of Cowper et a/./9. IO]. Matrix (T,I is a transformation matrix relating the coefficients Ai to the lucid degrees of freedom. II. ti~r/i~v I.. ~ir’~~jdy’) at each node. using the houndary component
conditions
and incorporeting
the two constraint
which are required, ;IS mentioned
equetinn5
per dispiaccmcnt
ithove. to restrict the normal 4ope variatic~u alang
edges 9-3 and 3-l. Matrix 1T,] i\ of order (60 ‘. !-II itnd is the inverse of it (60 x WI matrix ( ‘1’1with the last six columns of this inverse deleted: tl~tail~ of mntriu 1?‘I ;trc given in ‘l‘uhlc 1. ;Matri\ 1K 1 ic ;I rotation
matrix u~d to transform
global degrcrs Ilf frcetlom
to I[Ic:I~ ont’~ and is dct;~ilcd ill
Table 2.
S2 0
.
0’
Sj
0
0
0
[PI
0
s2
s4
0
0
0
34
0
0
0 0
0
si
0
Sti ,s
%
0
St
0
C 0
% 0
%
T:thle 2. Matrix I R I 71
IRI
=
0
0
,;
0
0
0
rL
0
C
0
0
'2
0
0
c
.:
0
0
.:
0
r..
0
0
0
0
0
0 0
9
0
0
':
0
0
0
0
0
0
8
0
G
0
0
0
0
cc2
SC
0
G
0
0
SC
3-
0
ir
IJ
c2
0
C
0
0
SC
0
J
c
,;’
^3c2
6
-3c
0
L’
$j
0
0
G
0
c
0
3
0
0
0
_f,>
_.2
i:
;.
0
:;.’
-nc
0
c.
i’
2
0
0
c
0
3
C
C
u
0
0
0
0
t
-3
[r*l
=
0
0
_,,2 $4, cl% -200’ c
-9
High-order triangular finite element for shell analysis
1103
Finally, the stress resultants (N,, N2 and S) and couples (MI, M2 and W) are obtained from the calculated strains using the relationships N, = C(E, + VEZ),
MI
iv2
MZ=$&+q,),
=
C(E2
s+1-
+
YE,),
= $
(XI
+
vx2),
(15)
Ch’
V)C$,
w=T(l-Y)T.
3. NUMERICAL
STUDIES
The two forms of the shell element described in the last section have been used in convergence studies of some well-documented problems. The elements are incorporated into the BERSAFE finite element system [ 191and calculations have been performed on an IBM 360 series computer using double precision arithmetic throughout. The distributions of displacements, stress resultants and couples are obtained by calculating these quantities at points along the element edges at intervals of one tenth of the side length. In representing a distributed loading the load vector is derived in a consistent manner through virtual work considerations, the necessary integrals being calculated in closed form. In considering the specification of boundary conditions only the necessary kinematic conditions have been applied; no attempt has been made to satisfy any force boundary conditions in addition. Infinitely -long, clamped circular cylinder
This singly-curved problem, illustrated in Fig. 2, was chosen as a first check on programme validity. With Poisson’s ratio Y = O-0 and all longitudinal displacements zero the problem is equivalent to the centrally-loaded, clamped arch. Four geometries are considered corresponding to the two values of the angle p and the two values of thickness shown in the figure. The shallow form of the shell element is used in analysing the problem corresponding to the smallest value of /3. Two gridworks are used in this symmetric problem, one being that shown in Fig. 2, and the other a coarse gridwork of similar form but with only 2 elements in the Y direction. Uniform
line load
/
Finite element
edge
/3= 900 or 7; h= l/l6 or I
u= 0.0
Fig. 2. Details of infinitely-long,clamped cylinder.
The percentage error in the calculated normal deflection under the load for the four geometries is given in Table 3. It is seen that the error increases with increase in angle /3 and with reduction in thickness. Although the magnitude of the error is small it is generally somewhat larger than that for the corresponding quintic-quintic arch[l, 21 since here in the shell analysis constrained-quintic displacement components are used. It is noted that the numerical studies
!)
tlsing arch elements would
be
very
1’ direction normal
I)\\!1
show that I!] Ihis :qqAi(::~II!)n the &ickn<) much
nearly-inextensional increased
I
lc\\
Ihan
geometry
are in excehs
the error\
of WY
for the corresponding displacement
rhat ;mtl
cil
of’ the \hell
111~ prcccni
c‘kmrnt~.
in u\itlg I\~O ;tnd four quintic-~uhic liYr
shell t’lenicn!\
I c+Avel!,
,md Ihc\c
\inct*
Ihcy ;lrt’ h,14
elements
For
the
of 19. IO)
deep.
thin.
arch element5
in the
LWCW\ will be somewhar on 2
component.
This doubly-curved
shalloh
shell i\ shown iI? Fig. 3. The loading is id uniform
acting over the whole shell surface
and the hhell edges :Ire freely
normal
(i.e. Wand
pressure V zero on
edge AD, W and I; zero on edge AH). Det;rileJ results for this problem are available t’ol comparison. both exact 1201 and hased on the 3, ~.lof’clcmen~ of (‘owper t;t d. 191.Two geometries are considered membrane primarily
corresponding
to \lalueh of il shell parameter
and bending behaviour as
ii
membrane
are more localised
are of imporr;lncc
Kkil.
in this diticult
wilh zone% ol’ hcnJil,_ 2 near the supported
for the shell N ith the ~l;~licr
equal to 042 and 0405. problem.
Both
the shell functioning
edge\:
these bending
1.alum I)I‘ K/r ii. ’ !i.c. for the thinner
zone>
shell if K and
1. are constant).
Uniform
normal pressure,
Fc,
3 x 3R mesh
The calculations mesh of elements problem
using the 54 Jot’ shallo\\, in :I symmetric
quadrant
cicment
are largely
and some use i$ also made of the rctincd clement
of the total strain energy calcltlatcd
hnscd on the use of
though this is not the most efficient
for holh shell gromctrie5
;I
uniform
arrangement
in thi\
meshe> shown in Fig. 3. Convergence u4ng hoth the present elementa
the 36 dof elements is shown in Fig. 4. It i\ clc‘;~r that in thia particular >hallow convergence of the strain energy kng the prcscnt clement i\ very littlc degree-of-freedom hasi\ th:m thaI of the carlicr L~Icmcnl. Some typical distribution\ of di\placemrnt \. force\ :Ind moment\ calculated element are illustrated in Fig\. 51 ~mclh \vhtw comp;lricon ih mude with the exact
and
application the different on in using the 5-I tlo~ solution. Quite
clearly the behaviour of the thinner shell is the more difficult to accommodate in the finite elrmcnt bvith it. This difficulty i5 largely overcome calculation because of the steeper gradient , ;I\\ociatcd by suitable refinement of the mesh and considering the \mall number of elements that are used in analysing this difficult problem the nccuraq of the results is high. The precision of the calculated bending moment distributions presented here is little different from that obtained by Cowper et 0/.[9] but there is a considerable improvement in the calculated force distribution. with much of the oscillation of the force eliminntcd.
1105
High-order triangularfinite elementfor she11 analysis I l Idmerh l
w
tie l
2x2
klQ .2x2 # 3x3
2x2 o
n 3x3 l
2x2
2X2R
0 .3X3R 3x30 .3X3R 3x3 0
03X3R
I 3
Rh/?=O~OZ
Final D&t 0 Rcf: 191 0 Present
of Freedom
’
Rh/p
=O-005
i
Fig. 4. Convergenceof strainenergy.sphericalcap probtem. C
B
3 10 -- % $Rh 0.05
(a) Rh/I%l.02
(blRh/f=OOoS Fig. 5. DistributionsalongBC for sphericalcap problem.
-.I’ ----
x
Exact solution 1x1 FiniiO oktall! 2x2 3x3 solutions 3dA
The shell roof is unalyxxl element
meshes.
normal
The
hrrc using the 5-l dof deep cylindrical
\elf-vhcight
loads of Irigoncn~it‘tric
vertical
as it Taylor’s
displacement.
series
loading
form
itnd c\prc\\ing
of terms
force and moment
i% xxotnmodated
up IO
at the central
section
seen that the finite element
r’rsult4 ;trc \‘cr! a~~ttrntt
generally
with
(0,
good agreement
Vert~col
Oeflect,on
and
the “CKICI”
Long,tudmal
the sino
the fourth
hhrll clement
hy resolving and
power
co\inc
ol’ the ;mgle frottt
of the angle.
BC arc shown
and unticlrttt
into tangentktl Distributions
in Fig. 7. where
with even the coarsesl
;IIKI
IIIC 01
ir c;tn IV
possible mesh givtng
t~wtl~\.
Forrt
4.
b)Bendmg momenta
Details of the rapid concrrgencc (A’ point v;tlue\ of displacement. force and moment :rnd 01 the total strain energy are recorded in Table 4. C‘omparative results for an "CXilCt" analysis 191 arc‘ listed but since these are based largely on shallo\\-shell theory convergence is to slightly differenl levels. It is noted that the rate of ccmvrrgcncc of the present finite element results (including particularly that of the strain energy) ia much more rapid, on a degree-of-freedom as well a\ on :tn
High-order
triangular
finite element
for shell anal!;sk
Table 4. Shell roof numerical
result<
element basis, than that of the 36 dof element results [9]. This is illustrated for one quantity-the vertical deflection at the centre of the straight edge-in Fig. 8. where results of various other finite element studies[23-251 are also included. The performance of the present element compares well with that of any other available element in this application. TOTAL 0
DEGREES
OF FREEDOM 400
200
Q
600
h31 Isaparametric.full mtegratlon
A 1231Isoparametnc, reduced integrat Ion
Fig. 8. Shell roof deflection
comparison.
Pinched cylinder with free ends In this problem the ends of the cylinder shown in Fig. 9 are completely free and the relevant data is L = IO.35in, R = 4.953 in. E = IO.5x IO”psi and I’ = 0.3125. This pinched cylinder problem has been studied in a number of finite element investigations. Most of these are concerned with a shell thickness of 0.094 in. (R/h = 52+3)[26-301but there also exist solutions corresponding to a thickness of O~OlS4Xin.(R/h = 3201[79.30]. The problem is almost inextensional and is sensitive to the representation of rigid-body motions: an inextensional solution by Timoshenko exists[31] though this gives a somewhat low value for the deflection under load. The problem is complicated by the very steep bending moment gradients local to the loading points and it is desirable to take account of this by employing some refinement of the finite element mesh in this region: such refinement is easy to apply where the elements are of triangular shape. Since accurate comparative values are available only for the magnitude of the deflection under the load, convergence of this quantity alone is examined for both thicknesses of shell in Table 5. The differences in the quoted comparative values of displacement reflect to some extent at least the small differences in the shell theories used. It is seen that convergence of the finite element results is rapid and that very few elements in the symmetric octant are required for a
Fig. Y. The pinched cylinder. ‘l‘;thk 5. Deflection of free pinched cylinder
sh.ite
pinal de@es of freedom
o1ezlent mesh
result of adequate engineering problems and particularly
accuracy.
DBflection under load (in.j R/I = 52.8 P = 100 lb.
R/h = 320 P = 0.1 lb.
It is noted that results for neurlv-inextensionnl
arch
for the pinched ring problem [?I indicate that $hcll elements based on
independently-interpolated restricted to a cubic variation
polynomial
displacements
in which the membrane
will not deal very efficiently
components
with the free pinched cylinder
;trc
prohlem.
Pinched cylinder with supported ends For this problem reference may again he made to Fig. 9 where in this case L/K 1. R/h = 100, v = 0.3. Here. though, the cylinder ends are supported by a diaphragm ti.c. the displacement components W and V are zero around the curved edge AD) which increase> the problem difficulty. Detailed solutions of thi\ prohlem hased both on ;I finite element study using their 36 dof element and on ;I double Fourier series solution are given hy Lindherg et rd. I I I I. The finite element results for the 54 dof cylindrical typical
distributions
and compared
of displacement
components
shell element are given in the form of
(Fig. 10) and of forces and moments (Fig. I I )
with the solutions of [I I]. In considering
uniform
meshes there is clearly
a very
significant improvement on an element basis in using the present element and there is also an improvement in the calculated deflection under load (and hence the strain energy) on ;I degree-of-freedom basis. The very considerable oscillation in the membrane force distribution calculated using the 36 dof element is drastically reduced in using the higher-order element. With regard to the force distributions
it is noted that the calculated
finite element values of membrane
force local to the applied load are converging on rather higher levels than those given by the Fourier series solution. Of course. here again. to make effective use of the refined 54 dof element some grading of the mesh should be used in the region of steep force and moment gradients. It can be seen from Figs. 10 and I I that when this is done considerable improvement in accuracy results, to the extent that the calculated strain energy for the 4 x 4A mesh solution differs only by, 0.31% compared with the Fourier series solution. (Note that the 2 x 2A mesh is basically a uniform 2 x 2 mesh but has the same kind of refinement in the rectangle local to the applied load as is shown for the 4x4A mesh in Fig. 9).
High-order
EhW P
triangular
finite element
1109
for shell analysis
’ -50
-100 -
-150
Sties solution
EhU P
REF.
PRESENT
Ill1
Fig. IO. Supported
pinched
cylinder.
displacement
distributions.
-10 -_._&2
Series sduhon
---3x3 ---_&XL
-20 0.2 !.k P
Finite dem~nt solutions
0.1
Fig.
I I.
Supported
pinched cylinder.
distribution;
along BC.
1. CONCI.l!SIONS
The numerical results indicate the high accuracy which can be obtained in a wide range of shell problems using a conforming, triangular finite element formulated in terms of conventional shell theory and based on polynomial surface displacement components which are each of constrained-quintic order. From the results presented here and from those of the limiting case of the arch it is concluded that the 54 dof clement is generally more efficient than are elements which have some or all polynomial displacement components of lower order. In particular whilst the
I I IO
efficiency
I). .I fh\if
of the 54 dof element
that of a related
36 dof element
gains in efficiency waviness
in applications
of calculated
of comparatively
will perhaps involving
displacement.
high-order.
differ
little in ~mc
(with cubic membrane
curved.
shells of deep.
force or moment displacement
>hi~llo~-\hell
components).
thin geometry.
distrihution~
element\
applications
from
there can he very cignifkrnr The oscillation
OI
which appears charactcrl\tlc
is Ggnificanth
reduced
in the 5-i dirl
element. .As with
other
high-order
connection
of displacement
variational
formulation.
over-compatibility
does
displacement derivatives
elements
From the results obtained not
circumstances
the connection
accommodate
some specific
appear
to
of particular physical
the
else
of
the present
at the node5 whose compatihilit~ lead
here and those obtained to an!
displacement
condition
significant derivative\
esccs~
element
involve\
i\ not required earlier
h! the
for ;uxhe\
stiffness.
III
t hi\ sot1ita
will need to he relaxed
10
of the khell.
4~~,lo#~/edgemenfs--The numerical r.e\uIb prr\cnled in thi\ paper wcrc oht;nned using the computing facilille> 611 lhc Central Electricity Generating Board at Berkeley Nucbr I.ahoratorie\. Berkclcy. Gloucestershire. The support of thi\ uork hy the Board is gratefully acknowledged and particular thank . :Lre due IO Wr. T K. Hcllcn f~)r hi\ help in incclrpclKlting the \hell elements into the BERSAFE finite clement *)‘stcm.
I, I). J. Dawe. Curved finite element3 for the analyris of \h;tlhn\ and deep srchch. i’orrrprrf~~\ ~rnd SrrrrcYwcs 4, ,559 I lY7Jl studies using circular arch fmi~r elements. Qvnprrr~n trrtd Slrrrrtures 1. 719 (19731. 2. D. J. Dawe. Numerical 3. J. H. Argyris and D. W. Scharpf. The SHERA farnil! of \hcll elements for the matrix tli\placcmcnt method. .41,rt, .I. 71. x73 (19681. 4. (;. Dupuis and .l. J. GoBI. .A curved tinitc clcmcn~ lor thm sI;L\tic *hell\. /III. J. So/id\ .SIUW. 6. 1113 I lWl1 5. ,A J. Morris. A deficient! in current tinite ~CIWIII~ is11-thin rhcll ;tipplic;~liunu. /!I/. .1. S~~litlv Srrm.!. 9.
t I%nI. G. E. Strickland
and W. .A. Loden. A doubly-curved trinngular \hcll element. /VII(. 3111 C‘r~r~l’ !\lrirri~ .M’th~dv or V~IOY. Wright-Patter~nn Air Force Base. Ohill t IWI \hcll di\crIIlc clcmcnt. ,\I \ \ .1. 5. -2’ ! I’hT! 3. k. R. Bopner. R. I.. Fo\ and I :\. Stl.:lin functilm\. I,ir .I .Mrt+i. Sci. 14. ITI llY??l. 31. R. D. Henshell. T. J. Bond and J. 0. Makoju, Ring linnc clcments for ;I\i\!mmctric thin \hrll analysi+. Prol,. ( tw,. C’nriotional Method: in Engitfeerinfi. l!niversity of Southampton 11972). and S. Woinow\ky-Krieger. ‘Fhcc1r.v c!fPfrrtc.v mu~She/ls 2nd Edn. McGraw-Hill. New York I 19591. 31. S. P. Timoshenko 2.5.
:MK/I.
TRIANGULAR FINITE FOR SHELL ANALYSIS
ELEMENT
D. J. D:\\rr Department
of Civil Engineering.
Ilniversity
of Krmir&m.
Hlrmingh.lnr.
I’r~gI:~ntl
Abstract-A curved-shell finite element of triangular shape i\ d~\cl-ii~itl !t hich 1%h;lwd (II) ~w~vcntional \hell theory expressed in terms of surface coordinate\ and dirpl;~cemcnr\. kwh ~rl’ IhL> ~hrw wt’.lcc’ dispkemrnt components is independently reprewnted I~> ;I twt~-Jill~c~~\i~ln;II pal! II~~~I.II IA ~~~~~~lr.;iiI1~~.J-r~uinli~ illdcr giving the element a total of 54 degrees of frwdom. T\to p;lrlicul,w ~c~~mc~ric l’~wm~ Ilf the clement XC considered, viz. doubly-curved shallow and circular cylindrkll. The h+h lc\cI t)f ~LL’LWI;I<~ 1%hiih can br achieved using few elements is Jemonstrakd in a range of prl,hlcm\ u hw L‘~w~p;lrwm i\ m;~de \vith preview finite element solutions.
I. INTR~)DL’CTION
Recent studies by the author[ I. 21 of the comparative have demonstrated displacement
the characteristics
fields and on assumed-strain
interpolated
displacement
components
efficient in a range of applications nearly-inextensional
behaviour.
eticirn
of v;u-ious arch finite elements
of models based both on rrlativcly
high-order
of quintic order h:~ve been shou,n to be consistently shallow :ind Jeep gccxnctrir\
embracing This efficiency
and
is \uch that in mc&!lling
in the calculated
applications[2]; polynomial
maximum
calculated
displacement
displacement
component
forces and bending moments arc‘ ;_Iso accunlrr.
on problem geometry all applications;
is restricted
force
components
the
Elements based on
i\ rcduccd to cubic order are
as QXII ;I\ on ;m clement basis. In particular,
to ;L cubic v;lriation
and the calculated
and
arches with a
i\ JCYSthan 04’~; in an) of the considered
fields in which one or both
very much less efficient on a degree-of-freedom the tangential
very
extensional
circular
single element based on surface (i.e. tangential and normal) quintic displacement error
polynomial
fields. Arch moc!cI , I~;~\ell on (two) independently-
di~trihutiow;
this waviness becomes very pronounced
tierzen
cfficirncy
I’ot, \tIch ;:iAcl\
if
Jcpcnds markedly exhibit
in nc~lrl!,-in~utensional
waviness in
problems where
its order can be very many times greater than the [rue mugrtitudc of the forcr’. This ia so whether the associated normal component of displacement ih of iubis 01’ quiiitic or&r and, in general. there is no improvement in increasing the order of the :l~!rm;lI ccrmponent I’rom cubic to quintic unless the order of the tangential a significant
improvement
component
incrc;i XXI 1though ~\ccptionally
is \imil;trly
where the arch geometry
it; sh;~llo\,~ .~nd r~lati\.cly
The arch studies point to the use of quintic polynomi:li for doubly-curved
shell elements.
elements of triangular
planform
lllspl;l~~mcnl
In fact there aIrcad> c&t
;r\~un;plic:n k
with quintic polynomial
component5
iii ihe literature I’IW
there is
thick). as the basis
two conforming
AI three displacement
components[3,4] but these components are carlesiun ones. C‘itrtrAn disp!;tcement components were chosen since they allow the precise rrprrscnt;ltion of lhc ri$-bcA!~ motions of an element ibf the when the geometry of the undeformccl shell is dctinccl ill linc,ir cc~mllination\ 111 :hc‘ norh of Argyris same set of basic functions used to define the ~li\pl;~ccmrnt~. and Scharpf[3] the development of the shell theory WI \\hich the c’lemcnt is based uses the natural strain concept and its relationship of Dupuis and Goel[4]
to clahsic;tl thin \hcll Ihcory iq not obvious.
is based on shell theory \\rittcn
in term\ of :I
The work
cuordinatr:
system
in place of the usual curvilinear system and the equation\ ;IIY crprc\\cd in relation to the height of the curved middle surface above a reference plane. Po>\iblc difficultlcs ;Ih\ociatrd with such schemes as these have been pointed out by Morris!!]. In ;I further stud! Dupuis[h] has described a triangular shell element based on a modified repre\ent;ltion 01 the bhcll geometry in which the element is straight-sided in the plane of shell parameter\. Kigid-hod! motions are represented precisely only by making an appropriate approximation IO the shell middle surface and although the derivation is, in principle, quite general it is spe&list’d to the LYISC\vhcrc the displacement components
are represented
by third-order
polynomials.
motion
states
(whose
inextensional
mode\.
states can be closely shell theory
nature
i\ dihcusacd
ik not neccwry approximated.
and the ahsumption
This C;II~ hc achieved
polynomi;rl wrface of “sensitive” easily
allows
thr: exilct
paper det;lil\
element
which helong\
surface
component
polynomial here
the
~huracteristic~ radius of twkt
limited
twent!
to t\vo
iii ~l~wd
to \tuJ!,
and bending
The :Ipproach
of st;md;trlI
IYIWXI on it)tl\~~litl~~l\:l;
ti.c.
-one term\
~~‘cwlpon~lll~~
classes
of \hcIi.
h;i\iciill\
.hc!i
\hcll
i L~,~k.!!
iI ~oiiiplctc ionstr;lint<
The philowph!
of
rri;lnylt,l!
par~;~gr;~pl:. l‘h~~ or&l
to M hich thrw
of freedom.
specific
\haIw
of :I ;onlornling
in the preceding
of
Hith
form
;imi numeric;il
with :I nltrncric;ll
in detail the preciw
momc‘nt~
de~crihed
and Olson19-121 assumption
,I\ -IIL’~I
i\ \utticientlv
i4 con~tr;rinecl-quintic having
small error a\wciated
displacement
wch
q:ri!it~~
;lrc ;ipIliiclI i
i\ iI general
gcomctr\
\\ hich
OIW IV!
exhibit
1111:
of general \hellh. \ iz. the shalltw shell with twj principal ratlii of cur\;ltuir’ I~III~. .: antI the circular c.ylindric;li khcll. ‘l‘hi\ rc\triction !\ ~mpcwd MI that the cIi-ritc’l>:
Gnce it is proposed forces
\tatc\
polynt~mi~~l surfaa
and application
ha\ ;L total of 54 degree\ i\
may he cibtaincd
any possible
in ;m clcmcnt
rrprewitation
IO the r~atcgor~ dcscrilwl polynomi;il
analyG3
3titfnesF
geometric
the tlewlopmun~
in two dimension\
and thus the element
components
h;l\ic
high. I I I> !‘llrthr*r notl:cI ix:’ l.lll! displncement fields arc p;irticuMy appropriiltc iii dcaliii g \%lth ;!I1 11111’1l~1 prohlcm~ in shell worklS.Sj antI thxt ~111anal!\~ \ 11;,,1’,1 ,111 i’,lll! ~‘!!!lrlll.~i
claw
The prewnt
of other
and rapid C~U~VCI~C’IICC’pro\‘iditg
of independently-intcrpol;lted
\o long a\ the order oi the polynomial
theory
in 171) or indeed
for high accuracy
hcrc ik tlircctlh
.;\ 0iAQ
schcmt.
of c:llcul;M
nwn:il
I’trta ~limiii;ll!tibl:
I\ ~~~~I~~,IcI~cJ JC,II .lhil
1114t’la~‘c’r,l~iit\. mL*mhi Lint’
stlrf;lcca
~~*‘/atcd to that xdnptcd
cliffawcc
:I con~trainril-qtiiiitic
components
wriation
over the shell middle
the fund:unentnl
intcgtxticiri
intqqticrn
that whilst
I’:II.~IC~ IT\ (‘,I\\
the txlrlicr
dicpl;icemtnt
\\trrk
component.
arc taken to hr cubic pol~nomial4.
plea’, I IINII~~I~~
i\ h;isctl 1111tht,
the
I\\ 11 t:lngcnti;ll
Lsfd. 1I?! arntciid
t’~~.pcr.
th;11. III!
;I wt numhcr of clement\. inzrca4ng the 11rdc1 of the t;lngcnti;ll ~‘~~iiipoii~iit~ IO L‘lw4!-:lin~~i quintic will lead to only :I mwginal impro\~ment in ;lccur;lc! achit~\~~~tl ;II the C.XIWIIW ,li .)I. increase ;i
in the degree< of frcrdom
T;iylor’~
limiting
of CW’I, II%\ contention
series hut i\ certainly
not correct
C~IW of the dtxp circulw
etticient
i\ has.~l on error ~l\nsider;ttic~n~
:ii 211 \ituaticjn\
arch ;I motlcl with quintic
on :I de~ree-rlf-f’r~~c(l~~nl
h;l\i\
th:ln
\in
:I\ m~nliciiicd
tangenti;li
J awrc\pondirig
IIVII;.
c;lrlic:.
cli~placc‘mcnt
II’ lhc
i\ nlu~b WOI~
I~~IJc:! u iti)
;;rnecnti;ii
displacement. Any its\essment
of the merit\
must take account geometry
and loading hehaviour.
hending.
etc.: rigid-body
structure
with
detailed
error
arch [ 1X-151).
or the shell might
prediction large.
Accordingly.
numerical
application earlier
which
include\
ma!’
function
clement\.
ha\c
to \peciIic
nil the c;ltc~oriec
hc tlominatctl primarily
;IIMI\
difficult
prewnthd
for
h!
plohlem~.
Hwcver.
prc\cntcd
noted ;Iho\ C. In the prcscntation
\\ ith Ioc:li /I~~!~‘\ 1): circurn\t;~ll~c.~
I\ noted
prohlcm
of
;I LI iLk
.!
wpplcment~cl I’:III~C
it’
th~~ilgh I~;!I thv
her‘. i\ cwmincd
wch ;lpplic;ltiotl\.
cast* of the ;Irch 1I. 21.rmhracc
:!!’ !j!
of ;I general illr: c:i ,hcll III
\implitied
II~)!.
hehaviorr-.
In thcw
modellin~
\IIIIL’III~~~-
I)ependill;
memhranc
to awtruct.
the
of the clement
XI\ ot \hell
;I\ ;I mcmhr;Inc
of :I linite clrmcnt i\ \~r\’
heen
the ilccur;lc!
for the limiting
for the
, of the ~lcmcnl
may or rn;~! not he of importnnc‘c.
of the ilccurilc\
unequal
estimilte\
presented
motions
propowd
of p~~~sihlc application
the shell rehponw
hending
v;llid theoretical
of ;I finite element
of the di\ercit!.
4lcII
..ll.i‘ll!.ll
,~~lci~~ !:I lb: ~.,*\uItll,:j!.:.
$11 ,lct~~ilotl iiumCri~;::
illlll
.:, .irlt~
p:~rticular attenticjn i\ paid tcr ;I
1L
lesser improvement in \hallo\r :lpplicaticrll~ \t;Igc one of the wfercl*~ ,)f this papc~’ h.1..
;I
1Ih] in cl:hich ;I farnil!, of \hcll
\I
)‘SlS
n1:1.
\lI
c!cmcnts
:.IILUII
rt4;ltcl.: ;I~ thut
s
the many shell theorie\ availahle in the literature that of Koitcrl 171 is :Idopted hcrc. In thih theory (unlike wnc others I the homogenecw\ strain-displacement equation\ ;irc’ \;lti\Ii(.(i in ;I general rigid-hod! moticjn. ‘I hc theor! ;14\unw an orthogonal curvilinwr ic~cwclin;ltc \\ \I;‘ITI !>\I: . .. ii,, thk \!‘\tern ntx~l I~OI cc)inciJl: \rith the /!I){: Y ,li. tirincip;ll ~~trr\:~tirr~~ -II !Ib.~ ,hcil bliit‘;i, Of
I099
High-order triangular finite element for shell analysis
notation
of Koiter
the strain energy
the sum of the extensional
density.
and flexural
dV,,, of a shell of thickness
energies
h may be expressed
as
in the form
or
Here
and the (6 x 6) matrix along orthogonal
[B] is easily constructed.
parametric
and T are the physical Also the quantity
where
curves
components
C ic defined
E and u are Young’s
As mentioned
earlier,
types of shell-those evaluated
in closed
conveniently
form.
expressed
strains
these curves:
and twist referred
xl, e
to these curves.
as
modulus
of shallow
strain-displacement
lI and lI are the extensional
of the changes of curvature
and Poisson’s
the present
simplified
The quantities
(Y, p and $ is the shear strain between
paper
geometry
specifically
and those of circular
equations For both
ratio respectively.
is concerned
with
two representative
cylindrical
geometry-so
can be used and thus the element
types
of shell the strain-displacement
that
stiffness
can be
equations
can be
as IE] =
[Wd)
13)
where (4)
For the shallow parallel
shell
to suitable
component.
global
Following
II and
V are the tangential
Cartesian
coordinates
the usual assumptions
ID1=
of shallow 0
0
l/R,
0
0
0
0
I
l/R:
0
0
0
0
I
I
0
2/T
0
0
0
0
I
0
0
0
0
0
I
0
0
I
0
0000
directions
orthogonal
to X:
[D]
and twist of the middle surface
runs parallel
U. V and
Y coordinate
to the longitudinal
W are displacement
and the shell normal respectively.
is
(5)
related to
axis and Y is the components The matrix
in the
[D] in
to be
0
I
0 I 0
I 0 0
l/R 0 0
0 0 I
0 0 0
0
-l/R
0
0
0
I
$R 0
-;R 0
00
0 0
0 0
0I
00I
00
[Dl=
0 0 00
0 -1 where
0
radii of curvature
shell the X coordinate
of the X coordinate,
this case is taken
0
the matrix
in directions displacement
1
0 0
where R,, R, and T are the (constant) the X. Y coordinates. coordinate
shell theory.
0
0000
For the cylindrical
components
-1
L(10
curvilinear
displacement
X and Y. and W is the normal
R is the mean radius of the cylinder.
0 0 a
(6)
I). .t.
I Ino
lhar:
For the tinitc clement ;tnulysis the domain of the shell in the XY plane (representing plane in the case of the shallow shell or the developed
plane for the cylindrical
the base
shell) is divided
into triangular elements :I\ shown in Fig. I. Within an element it is required that each of the three displacement components varies basically as a quintic polynomial in the two coordinates. As is well known. the use of ;I quintic polynomial displacement field originated in the development of refined conforming triangular plate bending elements by a number of independent investigators (see [Is]).
Corresponding
to ;I complete
degrees of freedom comprising together
quintic
the deflection
with the normal deric;ltivc
polynomial
such an element has twenty-one
with its first and second derivatives
at the mid-point
of each side. However,
at each vertex
a more convenient
eighteen dof (constr;linrtl-quinli~) tli+lc.cmcnt field is obtained when the mid-side degrees of freedom are rlimin;~tcJ 11~ re\tricting the \.uriution of normal slope along an edge to a cubic polynomial and it i\ thkb crln~tr;linecl-qtlintic polynomial displacement component. Cunneition of itll nodal quantities boundaries
of each displacement
component
of the normal derivative of the tangential with the potential energ!’ principle.
I’
/
which is selected here for each ensures continuity at inter-element
and of its normal derivative, components
is not strictly
although the continuity
necessary
for compliance
x,u
/
/”
/’
i”’ _.-_.
I
-.-____--..L_
-
Having dccidcd in hro;rd term\ the hiI& element
stiffness c:m he implcmcnted
extend and modify general closed-form
x,u
_.___
in
of the shell element the detailed calculation
:I
number of ways, but it is most convenient
of the here to
the nppro:lch described by Cowper et a/.[9. lo]. This approach leads to a integration formulas for the evaluation of stiffness terms via the introduction
of local coor-din;& ;IYC'Y I :~ntl y (we Fig. I). The corresponding local displacement components arc the tungcnti;d cmcs II ;~d I’ tllP;I\tIrCd in the directionr of s and y respectively, together with the component
norm;ll to illc &cll \urf;lce.
The component\
with like cuprccrions components
of tht.. conforming
IV. which coincides with the global component, displacement
for I‘ ;Ind I(’ giving 60 independent
;lr’e more \uit;lhl!,
W.
field are assumed to be of the form
coefficients
Ai in all. The displacement
c~:presd
where mi . . si are integer< runnln,c hetween t) and 5. The above polynomial
does not contain an
s4~ term so th;lt the norm:~l hlopc \;lriation along edge l-2 is automatically restricted to a cubic polynomial: two further cclnstr;lint conditions per component ;lre applied later to similarly restrict the normal slope variation
;dong the remaining
thus reduced to IX per component.
5-l in dl.
edges. The number of independent
coefficients
is
1101
High-order triangular finite element for shell analysis
Following the example of [lo] the strain energy for the shallow or the cylindrical shell element can then be expressed (by substituting eqn (5) or (6) into eqn (1) via eqn (3) and transforming to local displacements and coordinates) in the form of an integral over local coordinates as
[d][L]{d} dx d?,
(81
where (9) and IL1 = w41T[Dl’vww?41.
(10)
Matrices [I&] and [t] have the same general meaning here as in [lo]; the rotation matrix [I&], which links the column vectors {d} and {d} (i.e. {d} = [R,](d)) is in fact unchanged whilst matrix [L] does, of course, differ in detail both between the shallow and cylindrical shell analyses described here and the earlier analysis. The strain energy can be further expressed via eqn (7) in the form Vo= ;[Al[klIAl.
(11)
Here {A} is the column vector of the sixty coefficient Ai and [k] is a (60 X 60) “stiffness” matrix whose elements kii are explicitly defined in terms of the elements of the (8 x 8) matrix [i] by a modified form of the equation given in the Appendix of [ 101.Since this modified equation is a long one it will not be presented in full here; rather the necessary changes from the original equation will be noted. For the shallow shell element described here the expression for kij is that of the original equation but rows 8 to 13 inclusive as printed may be deleted since the elements of matrix [L] which are involved (L (1,6) through to L (4,B)) are in this case zero. For the cylindrical shell element presented here the expression for kri is the complete one of [IO] with the addition of the two terms tL(6,7)[n(ri
- l)r,q + rj(ri - t)ris]F(ri
t rj -3, si + Sj- I)
(12) + L(7, 8)[risisi(si
- l)+ rjsjsi(si - l)]F(ri + rj - I, si + SI- 3)
where F(m, n) is defined in [lo]. The “stiffness” matrix [k] corresponds to the quadratic form (eqn 11) of the strain energy in terms of the polynomial coefficients {A}. The energy can be further expressed in a quadratic form of the global degrees of freedom as v, = T”[ WZIP-GI~ W:I
(13)
[&I = [RlT[~~l’[kl~~4Rl
(14)
where,
is the required global stiffness matrix. Here the column vector of global degrees of freedom {Wz} of length 54 comprises the quantities au au a2u
u --
9ax’ay’ax”axau’p9
a2u
a2u
v
“..
a2v -@w**--ayz
a2w
IlO2
I). .I, IhW’F
at each of nodes I. 2 and 3 in turn. The matrices I I’, 1and I R I here. of course, take a different form
from that used in the work of Cowper et a/./9. IO]. Matrix (T,I is a transformation matrix relating the coefficients Ai to the lucid degrees of freedom. II. ti~r/i~v I.. ~ir’~~jdy’) at each node. using the houndary component
conditions
and incorporeting
the two constraint
which are required, ;IS mentioned
equetinn5
per dispiaccmcnt
ithove. to restrict the normal 4ope variatic~u alang
edges 9-3 and 3-l. Matrix 1T,] i\ of order (60 ‘. !-II itnd is the inverse of it (60 x WI matrix ( ‘1’1with the last six columns of this inverse deleted: tl~tail~ of mntriu 1?‘I ;trc given in ‘l‘uhlc 1. ;Matri\ 1K 1 ic ;I rotation
matrix u~d to transform
global degrcrs Ilf frcetlom
to I[Ic:I~ ont’~ and is dct;~ilcd ill
Table 2.
S2 0
.
0’
Sj
0
0
0
[PI
0
s2
s4
0
0
0
34
0
0
0 0
0
si
0
Sti ,s
%
0
St
0
C 0
% 0
%
T:thle 2. Matrix I R I 71
IRI
=
0
0
,;
0
0
0
rL
0
C
0
0
'2
0
0
c
.:
0
0
.:
0
r..
0
0
0
0
0
0 0
9
0
0
':
0
0
0
0
0
0
8
0
G
0
0
0
0
cc2
SC
0
G
0
0
SC
3-
0
ir
IJ
c2
0
C
0
0
SC
0
J
c
,;’
^3c2
6
-3c
0
L’
$j
0
0
G
0
c
0
3
0
0
0
_f,>
_.2
i:
;.
0
:;.’
-nc
0
c.
i’
2
0
0
c
0
3
C
C
u
0
0
0
0
t
-3
[r*l
=
0
0
_,,2 $4, cl% -200’ c
-9
High-order triangular finite element for shell analysis
1103
Finally, the stress resultants (N,, N2 and S) and couples (MI, M2 and W) are obtained from the calculated strains using the relationships N, = C(E, + VEZ),
MI
iv2
MZ=$&+q,),
=
C(E2
s+1-
+
YE,),
= $
(XI
+
vx2),
(15)
Ch’
V)C$,
w=T(l-Y)T.
3. NUMERICAL
STUDIES
The two forms of the shell element described in the last section have been used in convergence studies of some well-documented problems. The elements are incorporated into the BERSAFE finite element system [ 191and calculations have been performed on an IBM 360 series computer using double precision arithmetic throughout. The distributions of displacements, stress resultants and couples are obtained by calculating these quantities at points along the element edges at intervals of one tenth of the side length. In representing a distributed loading the load vector is derived in a consistent manner through virtual work considerations, the necessary integrals being calculated in closed form. In considering the specification of boundary conditions only the necessary kinematic conditions have been applied; no attempt has been made to satisfy any force boundary conditions in addition. Infinitely -long, clamped circular cylinder
This singly-curved problem, illustrated in Fig. 2, was chosen as a first check on programme validity. With Poisson’s ratio Y = O-0 and all longitudinal displacements zero the problem is equivalent to the centrally-loaded, clamped arch. Four geometries are considered corresponding to the two values of the angle p and the two values of thickness shown in the figure. The shallow form of the shell element is used in analysing the problem corresponding to the smallest value of /3. Two gridworks are used in this symmetric problem, one being that shown in Fig. 2, and the other a coarse gridwork of similar form but with only 2 elements in the Y direction. Uniform
line load
/
Finite element
edge
/3= 900 or 7; h= l/l6 or I
u= 0.0
Fig. 2. Details of infinitely-long,clamped cylinder.
The percentage error in the calculated normal deflection under the load for the four geometries is given in Table 3. It is seen that the error increases with increase in angle /3 and with reduction in thickness. Although the magnitude of the error is small it is generally somewhat larger than that for the corresponding quintic-quintic arch[l, 21 since here in the shell analysis constrained-quintic displacement components are used. It is noted that the numerical studies
!)
tlsing arch elements would
be
very
1’ direction normal
I)\\!1
show that I!] Ihis :qqAi(::~II!)n the &ickn<) much
nearly-inextensional increased
I
lc\\
Ihan
geometry
are in excehs
the error\
of WY
for the corresponding displacement
rhat ;mtl
cil
of’ the \hell
111~ prcccni
c‘kmrnt~.
in u\itlg I\~O ;tnd four quintic-~uhic liYr
shell t’lenicn!\
I c+Avel!,
,md Ihc\c
\inct*
Ihcy ;lrt’ h,14
elements
For
the
of 19. IO)
deep.
thin.
arch element5
in the
LWCW\ will be somewhar on 2
component.
This doubly-curved
shalloh
shell i\ shown iI? Fig. 3. The loading is id uniform
acting over the whole shell surface
and the hhell edges :Ire freely
normal
(i.e. Wand
pressure V zero on
edge AD, W and I; zero on edge AH). Det;rileJ results for this problem are available t’ol comparison. both exact 1201 and hased on the 3, ~.lof’clcmen~ of (‘owper t;t d. 191.Two geometries are considered membrane primarily
corresponding
to \lalueh of il shell parameter
and bending behaviour as
ii
membrane
are more localised
are of imporr;lncc
Kkil.
in this diticult
wilh zone% ol’ hcnJil,_ 2 near the supported
for the shell N ith the ~l;~licr
equal to 042 and 0405. problem.
Both
the shell functioning
edge\:
these bending
1.alum I)I‘ K/r ii. ’ !i.c. for the thinner
zone>
shell if K and
1. are constant).
Uniform
normal pressure,
Fc,
3 x 3R mesh
The calculations mesh of elements problem
using the 54 Jot’ shallo\\, in :I symmetric
quadrant
cicment
are largely
and some use i$ also made of the rctincd clement
of the total strain energy calcltlatcd
hnscd on the use of
though this is not the most efficient
for holh shell gromctrie5
;I
uniform
arrangement
in thi\
meshe> shown in Fig. 3. Convergence u4ng hoth the present elementa
the 36 dof elements is shown in Fig. 4. It i\ clc‘;~r that in thia particular >hallow convergence of the strain energy kng the prcscnt clement i\ very littlc degree-of-freedom hasi\ th:m thaI of the carlicr L~Icmcnl. Some typical distribution\ of di\placemrnt \. force\ :Ind moment\ calculated element are illustrated in Fig\. 51 ~mclh \vhtw comp;lricon ih mude with the exact
and
application the different on in using the 5-I tlo~ solution. Quite
clearly the behaviour of the thinner shell is the more difficult to accommodate in the finite elrmcnt bvith it. This difficulty i5 largely overcome calculation because of the steeper gradient , ;I\\ociatcd by suitable refinement of the mesh and considering the \mall number of elements that are used in analysing this difficult problem the nccuraq of the results is high. The precision of the calculated bending moment distributions presented here is little different from that obtained by Cowper et 0/.[9] but there is a considerable improvement in the calculated force distribution. with much of the oscillation of the force eliminntcd.
1105
High-order triangularfinite elementfor she11 analysis I l Idmerh l
w
tie l
2x2
klQ .2x2 # 3x3
2x2 o
n 3x3 l
2x2
2X2R
0 .3X3R 3x30 .3X3R 3x3 0
03X3R
I 3
Rh/?=O~OZ
Final D&t 0 Rcf: 191 0 Present
of Freedom
’
Rh/p
=O-005
i
Fig. 4. Convergenceof strainenergy.sphericalcap probtem. C
B
3 10 -- % $Rh 0.05
(a) Rh/I%l.02
(blRh/f=OOoS Fig. 5. DistributionsalongBC for sphericalcap problem.
-.I’ ----
x
Exact solution 1x1 FiniiO oktall! 2x2 3x3 solutions 3dA
The shell roof is unalyxxl element
meshes.
normal
The
hrrc using the 5-l dof deep cylindrical
\elf-vhcight
loads of Irigoncn~it‘tric
vertical
as it Taylor’s
displacement.
series
loading
form
itnd c\prc\\ing
of terms
force and moment
i% xxotnmodated
up IO
at the central
section
seen that the finite element
r’rsult4 ;trc \‘cr! a~~ttrntt
generally
with
(0,
good agreement
Vert~col
Oeflect,on
and
the “CKICI”
Long,tudmal
the sino
the fourth
hhrll clement
hy resolving and
power
co\inc
ol’ the ;mgle frottt
of the angle.
BC arc shown
and unticlrttt
into tangentktl Distributions
in Fig. 7. where
with even the coarsesl
;IIKI
IIIC 01
ir c;tn IV
possible mesh givtng
t~wtl~\.
Forrt
4.
b)Bendmg momenta
Details of the rapid concrrgencc (A’ point v;tlue\ of displacement. force and moment :rnd 01 the total strain energy are recorded in Table 4. C‘omparative results for an "CXilCt" analysis 191 arc‘ listed but since these are based largely on shallo\\-shell theory convergence is to slightly differenl levels. It is noted that the rate of ccmvrrgcncc of the present finite element results (including particularly that of the strain energy) ia much more rapid, on a degree-of-freedom as well a\ on :tn
High-order
triangular
finite element
for shell anal!;sk
Table 4. Shell roof numerical
result<
element basis, than that of the 36 dof element results [9]. This is illustrated for one quantity-the vertical deflection at the centre of the straight edge-in Fig. 8. where results of various other finite element studies[23-251 are also included. The performance of the present element compares well with that of any other available element in this application. TOTAL 0
DEGREES
OF FREEDOM 400
200
Q
600
h31 Isaparametric.full mtegratlon
A 1231Isoparametnc, reduced integrat Ion
Fig. 8. Shell roof deflection
comparison.
Pinched cylinder with free ends In this problem the ends of the cylinder shown in Fig. 9 are completely free and the relevant data is L = IO.35in, R = 4.953 in. E = IO.5x IO”psi and I’ = 0.3125. This pinched cylinder problem has been studied in a number of finite element investigations. Most of these are concerned with a shell thickness of 0.094 in. (R/h = 52+3)[26-301but there also exist solutions corresponding to a thickness of O~OlS4Xin.(R/h = 3201[79.30]. The problem is almost inextensional and is sensitive to the representation of rigid-body motions: an inextensional solution by Timoshenko exists[31] though this gives a somewhat low value for the deflection under load. The problem is complicated by the very steep bending moment gradients local to the loading points and it is desirable to take account of this by employing some refinement of the finite element mesh in this region: such refinement is easy to apply where the elements are of triangular shape. Since accurate comparative values are available only for the magnitude of the deflection under the load, convergence of this quantity alone is examined for both thicknesses of shell in Table 5. The differences in the quoted comparative values of displacement reflect to some extent at least the small differences in the shell theories used. It is seen that convergence of the finite element results is rapid and that very few elements in the symmetric octant are required for a
Fig. Y. The pinched cylinder. ‘l‘;thk 5. Deflection of free pinched cylinder
sh.ite
pinal de@es of freedom
o1ezlent mesh
result of adequate engineering problems and particularly
accuracy.
DBflection under load (in.j R/I = 52.8 P = 100 lb.
R/h = 320 P = 0.1 lb.
It is noted that results for neurlv-inextensionnl
arch
for the pinched ring problem [?I indicate that $hcll elements based on
independently-interpolated restricted to a cubic variation
polynomial
displacements
in which the membrane
will not deal very efficiently
components
with the free pinched cylinder
;trc
prohlem.
Pinched cylinder with supported ends For this problem reference may again he made to Fig. 9 where in this case L/K 1. R/h = 100, v = 0.3. Here. though, the cylinder ends are supported by a diaphragm ti.c. the displacement components W and V are zero around the curved edge AD) which increase> the problem difficulty. Detailed solutions of thi\ prohlem hased both on ;I finite element study using their 36 dof element and on ;I double Fourier series solution are given hy Lindherg et rd. I I I I. The finite element results for the 54 dof cylindrical typical
distributions
and compared
of displacement
components
shell element are given in the form of
(Fig. 10) and of forces and moments (Fig. I I )
with the solutions of [I I]. In considering
uniform
meshes there is clearly
a very
significant improvement on an element basis in using the present element and there is also an improvement in the calculated deflection under load (and hence the strain energy) on ;I degree-of-freedom basis. The very considerable oscillation in the membrane force distribution calculated using the 36 dof element is drastically reduced in using the higher-order element. With regard to the force distributions
it is noted that the calculated
finite element values of membrane
force local to the applied load are converging on rather higher levels than those given by the Fourier series solution. Of course. here again. to make effective use of the refined 54 dof element some grading of the mesh should be used in the region of steep force and moment gradients. It can be seen from Figs. 10 and I I that when this is done considerable improvement in accuracy results, to the extent that the calculated strain energy for the 4 x 4A mesh solution differs only by, 0.31% compared with the Fourier series solution. (Note that the 2 x 2A mesh is basically a uniform 2 x 2 mesh but has the same kind of refinement in the rectangle local to the applied load as is shown for the 4x4A mesh in Fig. 9).
High-order
EhW P
triangular
finite element
1109
for shell analysis
’ -50
-100 -
-150
Sties solution
EhU P
REF.
PRESENT
Ill1
Fig. IO. Supported
pinched
cylinder.
displacement
distributions.
-10 -_._&2
Series sduhon
---3x3 ---_&XL
-20 0.2 !.k P
Finite dem~nt solutions
0.1
Fig.
I I.
Supported
pinched cylinder.
distribution;
along BC.
1. CONCI.l!SIONS
The numerical results indicate the high accuracy which can be obtained in a wide range of shell problems using a conforming, triangular finite element formulated in terms of conventional shell theory and based on polynomial surface displacement components which are each of constrained-quintic order. From the results presented here and from those of the limiting case of the arch it is concluded that the 54 dof clement is generally more efficient than are elements which have some or all polynomial displacement components of lower order. In particular whilst the
I I IO
efficiency
I). .I fh\if
of the 54 dof element
that of a related
36 dof element
gains in efficiency waviness
in applications
of calculated
of comparatively
will perhaps involving
displacement.
high-order.
differ
little in ~mc
(with cubic membrane
curved.
shells of deep.
force or moment displacement
>hi~llo~-\hell
components).
thin geometry.
distrihution~
element\
applications
from
there can he very cignifkrnr The oscillation
OI
which appears charactcrl\tlc
is Ggnificanth
reduced
in the 5-i dirl
element. .As with
other
high-order
connection
of displacement
variational
formulation.
over-compatibility
does
displacement derivatives
elements
From the results obtained not
circumstances
the connection
accommodate
some specific
appear
to
of particular physical
the
else
of
the present
at the node5 whose compatihilit~ lead
here and those obtained to an!
displacement
condition
significant derivative\
esccs~
element
involve\
i\ not required earlier
h! the
for ;uxhe\
stiffness.
III
t hi\ sot1ita
will need to he relaxed
10
of the khell.
4~~,lo#~/edgemenfs--The numerical r.e\uIb prr\cnled in thi\ paper wcrc oht;nned using the computing facilille> 611 lhc Central Electricity Generating Board at Berkeley Nucbr I.ahoratorie\. Berkclcy. Gloucestershire. The support of thi\ uork hy the Board is gratefully acknowledged and particular thank . :Lre due IO Wr. T K. Hcllcn f~)r hi\ help in incclrpclKlting the \hell elements into the BERSAFE finite clement *)‘stcm.
I, I). J. Dawe. Curved finite element3 for the analyris of \h;tlhn\ and deep srchch. i’orrrprrf~~\ ~rnd SrrrrcYwcs 4, ,559 I lY7Jl studies using circular arch fmi~r elements. Qvnprrr~n trrtd Slrrrrtures 1. 719 (19731. 2. D. J. Dawe. Numerical 3. J. H. Argyris and D. W. Scharpf. The SHERA farnil! of \hcll elements for the matrix tli\placcmcnt method. .41,rt, .I. 71. x73 (19681. 4. (;. Dupuis and .l. J. GoBI. .A curved tinitc clcmcn~ lor thm sI;L\tic *hell\. /III. J. So/id\ .SIUW. 6. 1113 I lWl1 5. ,A J. Morris. A deficient! in current tinite ~CIWIII~ is11-thin rhcll ;tipplic;~liunu. /!I/. .1. S~~litlv Srrm.!. 9.
t I%nI. G. E. Strickland
and W. .A. Loden. A doubly-curved trinngular \hcll element. /VII(. 3111 C‘r~r~l’ !\lrirri~ .M’th~dv or V~IOY. Wright-Patter~nn Air Force Base. Ohill t IWI \hcll di\crIIlc clcmcnt. ,\I \ \ .1. 5. -2’ ! I’hT! 3. k. R. Bopner. R. I.. Fo\ and I :\. S
:MK/I.