Compurm & SlnrcUms Vol. 9. pp. 175482 @ Pergamon Press Ltd.. 1978. Printed in Great Britain
A MATCHING SUPERPARAMETRIC BEAM ELEMENT FOR SHELL BEAM SYSTEMS D. N. BU~GOHAI~,~ S. B. AG~WAL~ and R. S. AYYARQ Department of Civil Engineering, Indian Institute of Technology, Bombard 076, India (Received 1 Febrnary 1977; received for publication
16 September 1977)
Abstract-A three noded. 18 d.o.f. beam element is developed for beams of variable rectangular cross section, arbitrarily curved in space. The superparametric beam element is combined with the existing eight noded. 40 d.o.f. superparametric shell element for shell beam systems. To show the applicability of the formulations, some simple problems of beams and a four-beam slab are analysed.
2.2 Geometry of beam For any nodal position i the ha_lfwidth vector pz in n direction and half depth vector V, in i direction can be expressed as (Fig. I):
1. INTRODUCTION With the advent of computers, the finite element method
has become a most powerful numerical method of stress analysis. Among the large range of finite elements, the isoparametric family of elements is most efficient and versatile. For the analysis of thin and moderately thick shells, the eight noded superparametric shelf element is generally used[l, 21.To analyse shell beam systems using the superparametric shell element, the use of three dimensional elements has been suggested by Pawsey[2]. Aowever, this procedure results in a large number of unknowns; which do not improve the accuracy of the results in any way, since the displacements along the cross section of a beam can be expressed by a single node. The aim of this paper is to present a matching superparametric beam element (which is a degenerate three dimensional element) placed eccentric to the shelf element, The deveIopment of the beam element, as presented in this paper is not reported in any other available literature and is essentialiy based on the research work of AgrawaIj31. The beam element is arbitrarily curved in space and consists of three nodes along the axis of beam (Fig. 1). Each node is having six degrees of freedom, three translations u, v and w in the directions of global axes X, Y and 2 respectively and three rotations a, @ and y about the three local axes of the beam (Fig. 1). 2. FO~~~ONS 2.1 Assumptions The assumptions made for the development of beam element are: (i) The beam is of rectangular cross section, although such rectangular cross section can have variabIe dimensions along the axis of the beam which is arbitrarily curved in space. (ii) The two normais to the midsurface of beam (in the width and depth directions) remain stiff and straight during deformation, i.e. translation and rotation take place without any distortion of the cross section. In other words, warping of the cross section is neglected. However, during deformation the normals need not remain normals to the midsurface. Shear deformations are thus taken into account.
- fXiYi&lc
(1)
v3i = {&Yi&lr
- fXiYizilr.
(2)
The element geometry can be defined in terms of midsurface nodal coordinates and vectors Vzi and Vsi (Fig. 1)
The first term on the right hand side of eqn (3) defines the central &is of the beam while the fast two terms define the o~entation of the cross section. Ni are the quadratic shape functions in 4 only, given as N,=;5(1+0 N2=(l-[‘)
(4)
A+(l-8). 2.3Jispjacement field Vii, Vzi and Vsj are the vectors at node i in the directions of IocaI axes 5, n and l respectively (Fig. 1). Any line throu~_node i, ori~n~iy normal to the midsurface and along V,$ can have three independent translations and two independent rotations. They are q, vi, wi in the giobal directions X, Y and Z and ai and yj about the two IocaI axes .$ and 5; Similarly a line along v3j can have three independent translations ui, vi, wi in the globa directions X, Y, 2 and two rotations ai and j$ about the two local axes 6 and n. Thus the total number of midsurface displacements at node i are six: three translations ui, Vi and wi in the global directions X, Y and Z and three rotations LI~,pi and yi about the three local axes 4, n and y. The displacements of any point on pzi due to rotations ai and yi are: Displacement in e( tu) direction, u’ = -l/29biyi Displacement in l( V,,) direction, w’ = l/2nhiai >. (44
tAssistant Professor. *Post Doctoral Fellow. IProfessor of Civil Engineering. CAS Vol. 9, No. Z-E
17,i = {xiYi&lr
175
176
D. N.
i
/
(a)
BURAGOHAIN et al.
y,v
X,”
AND
LOCAL
lb)
PLAN
GLOBAL
VIEW
CO-ORDINATES
OF
BEAM
Fig. 1. Details of
FOR
BEAM
ELEMENT
beam element.
Similarly displacements of any point on v3,, due to rotations, ai and pi are: Displacement in [(vii) direction, U”= 1/2{didipi Displacement in q( fzi) direction, u” = - 1/2&Iiai I ’ W) The x, y, z components due to these two sets of displacements (given in eqns 4a and 4b) are
(6) 2.4 Strain displacement relation matrix [B] The procedure of calculating [B] matrix is similar to that explained by Cook[4]. However, for the sake of completeqess, it is given here in brief. The strains {E}in the global directions are expressed in terms of the derivatives in global coordinates as
where 1,_m,n are the direction cosines of the vectors VI, vz and V,. The complete 18 d.o.f. element displacement field may now be defined as
,
14 =
where Ni are shape functions defined by eqn (4), and [df] and [df*] are the matrices of direction cosines, given as
w.*1= (7)
A matchingsuperparametricbeam element for shell beam systems in which a comma denotes differentiation with respect to the subscript and
A typical block [Bi] is of the form tBi1 = I&l + lliB?l + JIB,**1
1
100000000 000010000
177
(14)
and it can be calculated after carrying out the multiplications required for eqn (13). 2.5 Element stiffness matrix [k] The element stiffness matrix is given by
The derivatives with respect to 6, q, g are related to the x, y, z Cartesian derivatives through the relation: or
(9)
{K&,KCY = VlIW.,Kz~=
where [J] is the Jacobian of the transformation of & 9, 4 derivatives to the x, y, z derivatives. Thus
(10) The complete set of derivatives of u, U, w with respect to x, y, z can be e?pressed in terms of 8, q, g derivatives using eqns (9) and (IO).
in which (D] is the stress-strain relation matrix in the global direction. The Jacobian [5] of eqn (9) contains q and g terms to the first power. For beams, where width b and depth d are small in comparison with the radii of curvature in the corresponding directions, these terms may be neglected in comparison with the terms to which they are added. Thus [5] becomes independent of 7 and 1: In that case, the matrix [&,I of eqn (14) becomes independent of q and g; $S*] becomes independent of [ and linearly dependent on 71and [[B**] becomes independent of q and linearly dependent on & In such case, explicit integration can be carried out in the width 71and depth I directions. Substituting eqn (14) into (15) and carrying out explicit integration in 7 and f directions,
where WI =
-[J*]
0
0
[J*]
-0
0
0 0
.
+:
[S**l’[D][B**]) det IJI d&
(16)
LJ*1 1 As explained by Pawsey[2] and Taylor et al. [ I] reduced
From eqn (5)
Ii +i f%.
0
0
0
0 0
0
4.4
0
0 0 0 0 0
0
0 0 0
0 0
0 0 0 0 0
0
Ni. 0 0
0
0 0
77N,e Ni 0 0 0 0
8
0
0
i=l
iii,
0
0’
0 0 1 0
0 N;: 0 0
0 0
A., ;
Substitution of eqn (12) into eqn (11) and the resultant product into eqn (7) results in {%VYu’YxvYYr’YuIT = $, [Bil~uiSwia&YiIT~ (13) The matrix IBI is made up of three 6 x 6 blocks of I&].
integration, i.e. 2-point Gauss integration in 4 is employed to evaluate the element stiffness matrix given in eqn (16). 2.6 Elastic properties matrix: [D] Since the cross section of the beam remains un-
D. N.
178
BURAGOHAINet 41,
distorted after deformation, the stresses to be considered in the local coordinates are a,, r, and rCC[B] Matrix yields strains in the global X, Y, Z directions. Therefore, the matrix [II*] is first taken in the local coordinates 5, n and 4 at the integrating point. For a beam element, the relevant stress-strain matrix in the local coordinates 6, q, t [a*] = [O*][E*] is given by ‘BEAM
000
0
oooI-vo
0
0
000
0
1
I
0
Pig. 2. Shell beam junction.
0
0
O
-57.
connecting node: (i) Transformation [TR,] of the five nodal displacements (3 translations and 2 rotations, Fig. 3) at the shell mid surface to six translations, three at the top and three at the bottom surface of the shell. (ii) Transformation [7X,] of these six translations to the corresponding beam nodal displacements: three translations and three rotations at the centre of cross section of the heam. The displacement field for the shell element is defined as [41(Fig. 3),
2G
000
I I
0
1-v
(17) where E is Young’s modulus, v is Poisson’s ratio and G = 1.2, a shape factor generally used to account for the parabolic variation of transverse shear over rectangular sections. The matrix [D] is then obtained as ID] = [TR]=[D*][TR]
{e}-~Ni{~}~tNiiil~l(;ii
(‘*) in which the transformation matrix TR is given in terms of the direction cosines of 6, 1 and 6 axes as
iI2
I
lz2
WI=
;:;,
2
2
ii:2 ::2 2$2
,;:“n,
241, 2m,m3
2n2n,
21311 2m3m,
2n3n1
km1 12mz
mlnl w2
m3n3 (, m+13m3 4md h,n2+ wb)
n,l, n24 n34 (~J2+ n,l,)
Od3 + n3M (lpm3 + 13m2) (wn3 + m3k) Grn, + l,m3) (rng, t m,n,) . hl, + n,M
where I,, m,, n, . . . are the direction cosines of the [, 3 and 5 axes at Gauss point.
(20)
1 J 7
where unr u,, w,,,, o,, pm are displaoements at shell midsurface node, ti is the shell thickness, Ni are shape functions for shell element and
2.7 Shell-beam junction A typical shell beam junction is shown in Fig. 2. In order to connect the shell element to the beam element two transformations are involved for shell element at the
[T]=
-1, -m,
[ nI
i y,v
/--x,u
Fig. 3. Local and global coordtpps
1, m,
of shell.
-6 I
(21)
A matchin
superparametric
179
beam element for shell beam systems
There l,,_m,, n, . . . are the direction cosines of vectors Vii and Vzi for shell. The transformation matrix [TR,] derived from eqn (18) is given as
3.APPLICATIONS
> r 1
0.1
i
.
3.1 Computer programming Using the proposed formulations of the superparametric beam element, a general computer program has been prepared on CDC-3600 (with a core memory of 32 IQ. In order to connect the shell and beam elements, it was necessary to develop an independent computer program for the su~~~rnet~c shell element[3] and then combine both the programs (for shell and beam) to analyse shell beam systems. Because of the limited computer facilities available (both in terms of core memory and computer time), the generalization of the program has to be restricted and only some simple problems have been analysed to show the applicability of the present formulations.
'ai
cm where uT, Y, wT and ng, iig, w, are the translations at the top and bottom of shell midsurface. The transformation [%$I is derived from beam-displacement field and is givbh as
3.2 Numerical examples The present beam element formulations have been applied to two examples of beams and one of four-beam slab. The results are compared with the available theoretical values. Solutions with refinement in element size for the same examples are also presented to study convergence. Explicit inte~ation in the width and depth directions of beam and in the thickness direction of shell is used to evaluate element stiffness.
Ui
@i
0
Pi Yi
Example
or
i&l,=
(23)
[TR21i{slb
Here I,, m,, n,. . . are direction cosines of 6, v, l axes of beam, and ug, & and qn & are 9, b coordinates of the two points of beam, *here it meets the shell (Fig. 2). The slain-displacement matrix [BJ for shell at the junction node thus becomes (24)
f_~ilshell =tBil,,,,(TR,l[TR,I*
-16
1. Straight cantilever beam This is the simplest possible example for testing the formulations. The geometrical and material properties of the beam are shown in Fig. 4 and the beam has been analysed for two load cases. Case I: Concentrated unit load at free end. Case II: Uniformly distributed load of unit intensity ~rou~out the length of beam. In each case, one, two and three elements, have been used resulting in 18, 30 and 42 unknowns respectively. For both the toad cases vertical displacement at free
A
Cl4
LOAD
CASE
I
: CONC.
LOAD
CASE
i(
: U.O.L.
E
LOAD OF
UNIT
=t.O
Fig. 4. Straight cantilever beam: Example
t.
P = 1.0 INTENSITY
180
D. N.
BU~GOHAIN et al.
end is compared with the exact solution[5] obtained with and without the effect of transverse shear in Table 1. There is no effect of refinement of etement size on vertical displacement at free end in both load cases and it is in excellent agreement with the exact solution, even when the effect of transverse shear is also considered. For both the load cases, the convergence of bending moment at fixed end is shown in Pig. 5. There is no effect of refinement of element size on moment at fixed end for load Case I (concentrated load
at free end) and it is same as that of exact moment. For load Case II (uniformly distributed load) it converges to the exact vahte with the refinement in element size. Example 2. Cantilever beam curved in plan
The geometrical and material properties of the beam are shown in Fig. 6 and is subjected to a uniformly distributed load of unit intensity. The exact solution is given by Pippard and Bakerl51. The beam has been analysed for one, two, three and four elements resulting
TabIe 1. Comparison of vertical displacement at free end for straight cantilever beam, Example f (Fig. 4) Vertical displacement Load case
Number of elements
Beam element formulation
1
2666 2066 2066 12440 12440 12440
Case I Cont. load Case II U.D.L.
2 3
1 2 3
TOTAL
No.
OF
UNKNOWNS
Exact solution without shear effect
Exact solution with shear effect
2048
2064
12288
12416
-
0.9 -
8 I
OS-
a= L
CALCULATED EXACT
MOMENT MOMENT
= LOAD
CASE
I
P = LOAD
CASE
I[
l
J
Fig. 5. Convergence of moment at fixed end: Example 1.
FIXED END
U. O.L.
OF UNIT
Fig. 6. Cantilever beam curved in plan: Example 2.
INTENSITY
181
A matching superparametric beam element for shell beam SYStemS
ing all the major difficulties encountered in any shell beam system. The solution of this problem is given by William and Scordelis[61, who used finite strip method for the analysis. The geometrical and material properties are shown in Fig. 8. The structure is simply supported at the two ends and is subjected to a uniformly distributed loading of 1000lbfsq. ft. Because of symme~y, a quarter of the structure has been analysed using 3x3 mesh resulting in 214 un~owns. The results of the present
in 18, 30, 42 and 54 unknowns respectively. The con-
vergence of bending moment and twisting moment at the fixed end with refinement in element size is shown in Fig. 7. The bending moment and twisting moment; both converge to the exact solution with refinement in element size. ~x~rnpie 3. Fair-berm stab bridge This is the simplest form of shell beam system exhibitTOTAL
NO.
OF
UNKNOWNS
CALCULATED
u=
.
M i
VALUE
THEORETICAL
0 - BENDING
Fig. 7.
EXACT
42 I
12 I
18 I
I
-
VALUE
MOMENT
= TWISTING
MOMENT
Convergence of bending and twisting moments at fixed end: Example 2.
E = 30,000
c--
3.33’---4---
t-------
LBSISQ.
FT.
3.33’ +
3.33’ 1
10.6’ -p+ 1000
L8S I SO. FT.
Fig. 8. Geometry and loading of four-beam slab bridge: Example 3. Table 2. Comparison of numerical results for four-beam slab bridge, Example 3 (Fig. 8) Beam Quantity ‘w’ displacement (ft x lo+‘) ~n~tudin~ stress a,, at the mid-sugar of slab (1080Ibfsq. ft x Id) Longitudinal stress ~~xaat the_bottom of beams f 100Iblsq. ft x 10))
Method Scordelis[6] Present (3 x 3 mesh) Scordetis Present Scordelis Present
a
b
8.69
8.69
8.69
8.69
8.73 -1.36 -1.43 4.28 4.40
8.73 -1.32 -1.43 4.30 4.48
8.73 -1.32 -1.43 4.38 4.48
8.73 -1.36 -1.43 4.28 4.48
C
d
182
D. N. BURAGOHAINet al.
analysis for vertical displacement w, longitudinal stress a,, at the mid-surface of slab and longitudinal stress CT,, at the bottom of beams at mid-span are compared with those obtained by Scordelis in Table 2. The two solutions are in excellent agreement, the present solution being slightly on the higher side. However, the small disparity in the trend can be explained by the presence of lateral bending effect. For the span to width ratio (100: 10) which is considerable in this example, such an effect is not significantly reflected in the final results based on the present formulations. 4. CONCLUSIONS
The superparametric beam element, as developed herein can be very conveniently used for the analysis of beams of variable rectangular cross-section, arbitrarily curved in space. The beam element gives very good results even for parabolic variation of bending and twisting moment without any large increase in the number of unknowns. Since effects of shear deformations are also considered in the formulations, moderately wide and deep beams can be analysed. Analysis of space grids with curved members is also possible. In such cases the six degrees of freedom at the central node of beam can be condensed out prior to assembly resulting in a reduced number of unknowns. However, to analyse shell beam systems all the degrees of freedom of beam element will be required.
With the shell-beam formulations, any plate-beam or shell-beam system with arbitrary shape and arbitrary arrangement of beams can be analysed. Acknowledgement-Theauthors acknowledge with thanks the grant of a Post Doctoral Fellowship to Dr. S. B. Agrawal by the Council of Scientific and Industrial Research, New Delhi, for carrying out research work, a part of which is reported in this paper. REFERENCES 1. 0. C. Zienkiewicz,
2.
3.
4. 5. 6.
R. L. Taylor and J. M. Too, Reduced integration technique in general analysis of plates and shells. Inl. I. Num. Meth. Engng Sci. 3, 275-290 (1971). S. F. Pawsey,The analysis of moderately thick to thin shells by the finite element method. Report No. UCSESM 70-12. Structural Engineering Laboratory, Univ. of California, Berkeley. S. B. Agrawal. Analysis of shell beam systems by discrete energy and finite element methods. Ph.D. Thesis, Dept. of Civil Engng I.I.T., Bombay (1976). R. D. Cook, Concepts and Applications of Finite Element Analysis. Wiley, New York (1974). J. S. Pippard and J. F. Baker, The analysis of Engineering Structures. Arnold, London (1957). K. J. William and A. C. Scordelis, Analysis of orthotropic folded plate with ec entric stiffeners. Report No. SESM 70-2. College of Engng., Office of Research Services, University of California, Berkeley.