A special purpose element for shell-beam systems

Com,wers & Structures Vol. 29. No. 2, Printed in Great Britain.

pp.301-308.

00457949/a $3.00 + 0.00 0 1988 Pergamon Press plc

1988

A SPECIAL PURPOSE ELEMENT SHELL-BEAM SYSTEMS

FOR

GABOR M. V~R&

Department of Mechanical Engineering, Technical University, H-l 111 Budapest, Hungary (Received 17 July 1987) Abstract-This paper presents an isoparametric finite element for reinforced shells and plates. The formulation is based on general beam theory and takes into account both transverse shear deformation and torsional warping. The element exhibits complete two-line compatibility. Numerical examples are presented in order to demonstrate the validity of the formulation and the possibilities of application.

I. INTRODUCTION In recent years there has been a tendency to develop

special purpose finite element techniques for the analysis of reinforced plates and shells. The traditional procedure for analysis is the orthotropic model, where the stiffeners are not considered as discrete members but their effect is averaged or smeared out. The case of relatively widely spaced reinforcing bars requires another local approach. Most finite element models utilise shell and beam elements to represent reinforced structure. The nodes of the beam element, through constraint conditions, are forced to undergo the displacements and rotations prescribed by the corresponding shell nodes. The constraint condition can be introduced either through manipulation of the discretized equations as in [1,2] or through the governing variational principle employing Lagrange multipliers [3]. Beam elements, compatible with the thick shell elements, may be obtained either as a special form of the general three dimensional isoparametric elements or as an application of the one dimensional isoparametric concept and an appropriate beam theory. Among others Buragohain et al. [l] and Ferguson and Clark [4] developed curved beam elements of variable rectangular cross section, which have been obtained as a special form of the general three dimensional isoparametric element. Although the shear deformations are taken into account, the distortion of the cross section and the torsional properties are neglected. Thus the degenerate three dimensional elements, placed eccentric to the shell, fail to conveniently represent the torsional response of the beam. A different approach was introduced by Jirousek [S] in which the formulation of the element was based on an appropriate beam theory which properly accounts for the reinforcement eccentricity with respect to the middle surface of the shell, the transverse shear and the effect of the shear center location. Although the torsional response for any compact cross section was taken into account in the generalized

stress-strain relationship, this family neglects the torsional warning in displacement compatibility conditions. The aforementioned concepts of modelling the reinforcing beams use the requirement of displacement compatibility along a so called joint line. In most practical cases the line compatibility gives acceptable results, just as the area of the beam-shell joint surface is negligible as compared to the total surface of the structure. In certain cases more realistic results can be expected if two joint nodes are used in

each nodal beam section and the shell between the midsurface normals of joint points is considered as part of the stiffener section. The perfect displacement continuity can be achieved via the nodal point continuity only if the same shape functions, say Lagrangian, and the usual six degrees of freedom are used in the shell and beam nodes respectively. If the influence of the torsional beam warping on deformation is included an additional beam degrees of freedom is needed. Taking this seventh parameter as the rate of twist-the change of torsional rotation-a Hermitian shape function must be used which violates the full displacement continuity. This problem is surmounted by taking the seventh displacement parameter corresponding to the bimoment as an independent variable. The same idea was used for the slope in the Timoshenko beam and Mindlin shell theories. The goal of the study undertaken at the Chair of Mechanics of TU Budapest may be summarized as the development of an approach to the joint problem in which the displacement compatibility is assured along the real joint surface. To do this proper beam and shell elements are needed with suitable nodal degrees of freedom and shape functions. The advantage of this concept is that the cross effect between the inplane shell displacements and the torsional warping of the reinforcement can be taken into consideration. An additional advantage is that some engineering oriented results, such as the internal forces of the beam and the shear and normal force distribution between the beam and the shell or plate, may be implemented. In Sec. 3 a new stiffened shell element is derived with 301

302

GABORM. WR&

Fig. 1. Geometry of general quadratic beam element.

different nodal degrees of freedom by combining a thick shell element with a beam element. Although the process described hereafter may be applied to any number of beam element nodes only the quadratic case is detailed here. 2. DISCRIPTION

OF THE BEAM ELEMENT

The development of the beam element as presented here is essentially based on the research work of Kiss [6]. Figure 1 shows a typical slender beam element of general section. Let G be a point on the centroidal axis, the centroid of the section, and x, y, z a local Cartesian frame in which z is normal to the cross section while x and y are its principal axes of inertia. This leads to the matrix of orthogonal transformation T=[e,,q,e,l,

(1)

nates of a point in the cross section. The vector e3 is the unit normal of the cross section which has been assumed tangential to the centroidal line. The last term in (2) represents the tangential displacement where cp(x, y) is the warping function for free torsion of initially straight beams. Its application to curved beams is merely an approximation [7j. The seven independent generalized displacement parameters are u, v, w, a, fl, y, and 9, functions of the arc length z only, and the shape functions are 1, x, y and rp which fulfill the orthogonality conditions IxdA=j”ydA

I”xcidR

=[“qd,4

=j”yqdA

=O.

(4b)

24= UC- yy L,=u,+yx

(5)

w=w,+ay-Bx+lkp xR+@e,,

(4a)

In the local frame the strains following from the general displacements

in which e, , e2, and e3 stand for unit vectors lying in the x, y and z local directions respectively. Neglecting the transverse deformations of the cross section the assumed displacement field is u=u,+$

=IxydA

(2)

can be written as where U&Z) = u e, + 0 e2 + w e3 = uli, + u2i2+ u3i3 ~(z)=ae,+Be,+re,=~,i,+~zi,+~3i, R=xe,+ye, &YZ

e, = l,i, + m,i, + n,i,

are the displacement,

=vb+y’x

+a

aq

+-9. ay

(3)

rotation and the local coordi-

The prime stands for derivation with respect to arc length z. Rearranging the shear strain expressions, so

A special purpose element for shell-beam systems

303

stresses-or resultants-an be calculated. The components by use of orthogonality conditions (4) are: -shear forces F, =

GE,dA

= GA&:*

sA Fy =

GcyrdA = GA+ sA

+& -B - SYS) -normal

force F, =

Ec,dA

= EAw;;;

sA

+(ub + a + 9x,),

(84

where the first terms according to St Venant’s theorem are the strain distribution from free torsion, the second terms with parameter (9 - y’) are an approximative contribution from constrained torsion while the third terms correspond to shear forces. The St Venant torsion problem determines the stress and strain distributions only [S], while the displacement field requires the additional specification of the local centre of axial rotation. This point called the torsion (shear) centre S is defined by the orthogonality conditions (4b). Here we note that different-but identical with eqn (8a)-decompositions are given in El: E xz =

acp z -(Y

,

Y

-bending

moments

M, =

Ec,y dA = EI,a’ sA

My =

E&,x dA = EI,/V; sA

-torsion

moment

M, =

G(E~~x - eXzy) dA sA

=GJy’-GZ,($-y’)=M,-M,;

-bimoment

-Ys)

[

1

+@

-Y’g

B=

Ee,qdA=Er9’.

(9)

sA

+(4-B ayz=r’

-Y’Y,)

arp(x -ay

The cross sectional properties arising here are

-xs)

[

+(s 1

+(&+a

y$

+y’x,)

@b) I,=

and in [9]:

Exz=Q [

[

l-=

c cp2dA

$Y-Yd 1 +@-Y’)(Y-Ys)

+wi -B Ey:‘9

c x2dA

g-(x-x,)

+(&+a

J = Z, + Z, - Z, .

- Y’YYS) -(I9-y’)(x-x,)

1

w

fy’x,).

Both decompositions fail to conveniently represent the transverse shear forces of the beam. This is due to the fact that not only the third but the second terms too are equivalent to inplane internal forces. Moreover in (8~) instead of ‘, the rate of twist, the warping parameter 9 appears in the first term. For an elastic homogeneous material with modulus of elasticity E and shear modulus G the generalized

y

(10)

The resultants listed in (9) justify the seemingly arbitrary decomposition of the shear stresses @a). Namely, there is no shear force from the first two terms and no torque from the third term of (8a). The strain energy associated with a cross section is

=O.S sA

[EE~+G(E~~+E~~) z Xi Y:

+G(E,“;~ + E;*) + G(E:: + E;:)] dA.

(11)

GABOR

304

The integrals of products .s~$$ Ed&, a:~&, . . . are vanishing. A better representation of the shear strain energy contribution can be obtained by replacing the last term in (11) by

M. VTIROS

The generalized local strains in accordance with (7) and (g) are:

where D, and D,, are the shear deformation coefficients. Finally, the strain energy (per unit length) expression is obtained by the integration of (11):

+Er9’2

+ GJy’* 4- GZJ9 - Y’)~

-k G&(u;; - B - +.8

t CD,&

+ CL+ ~xs)~],

(13)

according to [q. Even though the isoparametric finite element formulation can be found in the literature, the basic matrices are presented here for the sake of completeness. The beam element is arbitrarily curved in space and consists of three nodes along the axis line. Each node has seven degrees of freedom. The centriodal axis will be defined in terms of dimensionless arc length 1; as

By using eqns (15) and (f6), a typical submatrix BG, is expressed as

B, =

0

0

Ui

0

0

0

0 0

0

0

0

aj

0

0

0

0

0

0

ilj

0

0

0

0

0

0

0

a,

0

0

0

0

0

0

0

a,

0

0

0

0

0

-f+

Ni

a,

0

0

0

-N,

0

YSN,

0

a,

0

Nj

0

0

M%

where

where N,(c) are the quadratic Lagrange interpolating ~~ynomia~s and X, YGItZ, are the global Cartesian coordinates of the nodes. For any cross section ofthe element the (10) cross sectional properties and the shear deformation coefficients specified in nodes are interpolated in the same manner. At any cross section let AC be a vector of seven independent parameters comprising the global displacements, rotations and the warping parameter as defined in (3):

z’ being the modulus of tangent vector. The element strain energy expresion is obtained from (13) as

where where AGi is the vector of the nodal degrees of freedoms A~=[u,,,u,,u,,,~t,rl/z,4k~,91i.

(16)

The nodal displacements and rotations associated with the local frame x, y and z are related to the global components by the orthogonal transformation (1):

and B,$ is the 21 x 21 stiffness matrix of the beam element. The vector U, contains the nodal variables (16) as

When the beam element is used as a thick shell stiffener any nodal cross section will be associated with two joint line points A, and & Along the joint fines the requirement of displacement compatibility between the beam and shell elements

305

A special purpose element for shell-beam systems

Fig. 2. Definition of joint lines of reinforced shell element.

RAi= [XGA,Y,,, .&Ii and as shown in Fig. 2, will be used to specify the position of joint points with respect to Gi. Let AAi and A, be the vectors of Ai and Bi point displacement parameters respectively: exists.

R,=

The

vectors

100

[Xcs, YcB, Z,&,

A,i = [UAI,u,,,

~,a,

$1,

$2,

0

@ABi

=

A,91

Asi=[uBI,UBZI1183,~1,~2,~j,Ql,

(20)

1 0

-

o

z”B - yAB f3q"B

z

o

AB

xAB

m3(PAB

001

YAB

-x”B

0

0 0 0

1

0

0

n3qAB 0

0 0 0

0

1

0

0

0 0 0

0

0

1

0

0 0 0

0

0

0

1 (23)

and from the beam displacement field definition, (2) and (3), one obtains the transformations AGi= 8,,AAi

and

AGi= eBiABir

(21)

By using transformation (21) to beam element strain energy (19), the element stiffness matrix at joint points becomes

where

100 0

1 0

0 0 8, =

0 z,,

1 -Y,,

- Z GB yGB

-l,cp,

and the strain energy in terms of B joint line variables

0 - XGB -m3vB x GB 0 +wB

000

1

0

0

0

000

0

1

0

0

000

0

0

1

0

000

0

0

0

1

LJB= 0.5 UTKBU B B Bv

(25)

u; = [A&, A&,

(26)

where

(22) and the same holds for eAi. The vector AAi is related to vector A, by a simple transformation

If the one joint line concept is used then RAi = R, and for membrane stiffener RAi = R, = 0. The derivation of the nodal forces due to loads follows the standard methodology and transformations outlined before.

3. THICK STIFFENED

A,i = @ii’@BiABi = eABiABi,

where

A&].

SHELL

ELEMENT

The nine noded Lagrangian isoparametric shell element is used in order to analyse shell-beam systems [lo]. The points of the shell element are

GABORM. ViS~ijs

306

0.30 m Case b

,x

Bear!Y!Y*,zB,

Fig. 3. Excentric concentrated load F = 1000kg on a double-T section cantilever: E = 210,OOO.O kg/cm*, v = 0.15 (Example 1).

interpolated from global coordinates of nodes on the midsurface by Lagrangian shape functions as UT = [A:, A;, . . . , A;].

(29)

The formulation of the strain energy follows the standard methodology which finally can be simply written as where t is the shell thickness and I,, m,, n, are the direction cosines of the e3 vector pointing in the midsurface normal direction. The local coordinate system required for the use of shell theory can be derived from (27). After taking [ = 0, the calculation of unit base vectors is as follows:

e3 =

(a x g2)/k1x g21

e, = e3 x i,

e, = e3 x e,.

The displacement field is defined as

Us = 0.5 UTKSU,

(30)

where KS is the 54 x 54 stiffness matrix. The three rotational degrees of freedoms may cause singularities in the stiffness matrix which are removed by introducing small fictious stiffnesses. This is done by using the average rotational stiffnesses multiplied by a properly small constant. In order to illustrate the process evaluating the stiffened shell element, suppose that the joint line is defined by shell nodes 3, 6, 9 as it is shown on Fig. 2. The seventh displacement parameter of the beam node creates difficulties as it does not appear in shell node vector (29). Thus the shell matrix must be expanded to a 57 x 57 matrix by inserting zero rows and columns at the appropriate places, and defining new nodal vectors in shell joint nodes as A,i = In,,,

uz,2,

u,,3,

ti,,

$2,

$3,

i =

SIT,

3,699,

and a new extended element vector as U,T= [A:, AT, A:,

where uli, uzi, nji are the global displacements at midsurface node and tili, JIZi, J/3i are three global rotations. Let Ai and U denote the vectors containing six nodal and 54 element degrees of freedom respectively,

A:,

A;,

A&,

A:,

A:,

&I.

(31)

After joining the beam and shell elements the external nodes of the new element on the joint face are B,,$ and B9. In order to eliminate the Ai internal joint nodes the transformation U, = 8 U, is used, where U; = [AT, A:,

Ai39

and according to (23)

A:,

A;,

AL,

A;.

Ai-,

A&l,

A special purpose element for shell-beam systems Table 1. Deflection and stresses for the double-T section cantilever of Fig. 3 (Example 1) K&m) -0.3043 a -0.2888 b -0.2722 -0.2780 t; Ref. [S] -0.2813

@=

a,,

adktd~~z~

-7.74 -8.35 -8.06 -8.22

-4.02 -3.47 -3.73 -3.67

Q&V, =G 2.55 2.62 2.37 2.48

0.97 0.83 1.09 0.86

rk I, @Ala,,RI, @,,, 1, I, @m J.

(32)

Here I denotes the 6 x 6 unit matrix. Note that t~nsfo~ation 6 includes the joint nodal values of the cp warping function as well. This fact guarantees that the relative displacements of joint nodes depend on the warping ability of the stiffener. In order words, the displacements of the stiffened shell element are constrained by the eqn (2) beam displacement field between the A and B joint lines. The eqn (30) shell strain energy becomes

307

warping of the solid rectangle is small enough to take r = 0 in (10) and rpAB= 0 in transformation (23). Therefore the number of degrees of freedom is equal to 254 in cases (b) and (c). (d) The entire structure is modelled exclusively by 15 quadratic thick shell elements resulting in 372 degrees of freedom. Table 1 presents a compa~son of results for vertical displacement under the concentrated load, longitudinal membrane stress at the midsurface points of the slab and longitudinal stress at the bottom of the beams at the fixed end of the cantilever. The aA deflection is close to the result obtained in [5]. The stiffer behaviour observed in case (c) is obviously due to the full displacement compatibility between the slab and the beams. Example 2 This example was chosen to study both the influence of reiatively wide beams and the beam warping on the whole structure. The problem is solved by using three finite element representations:

2Us = U;K;U, = U;@-K$@ U, = U;K;UB, (33) where Ki is the extended 57 x 57 stiffness matrix. This is the final form of shell strain energy to which the eqn (25) beam strain energy can be added. 4. APPLICATIONS

In application of the reinforced shell element the first problem is the accurate specification of the cross sectional properties. Recently a more complete formulation including the effect of transverse contraction appeared in [I 13. Two numerical examples are presented in order to show the accuracy and the additional benefits of the two-line compatibility formulation presented in this paper. Ex~~Ie

I

This example was taken from [5] where the author solved the problem by using three quadratic-cubic and 12 cubic-cubic thick shell elements resulting in 540 degrees of freedom. Here the excentrically loaded double-T section cantilever is analysed by using four different repre~ntations: (a) The structure is considered as a straight cantilever beam with thin walled section under bending and constrained torsion [7]. (b) In the second case, the structure is represented by an assembly of nine quadratic thick shell elements and six quadratic reinforced shell elements. In this case the one line compatibility is assured and the reinforcement section is a 0.3 x 1.13 m* rectangle. (c) As in case (b), but the two-line compatibility is assured and the reinforcement section is taken as a 0.3 x 1.27 m* rectangle. Figure 3 shows the details of element layout at one typical junction. The torsional

(a) In the first case the structure is represented by nine quadratic thick shell and six quadratic reinforced shell elements where the external nodes coincide with internal joint nodes. This fact ensures the one-line compatibility of displacements. (b) AS in case (a) but the two-line compatibility is assured. In cases (a) and (b) the torsional warping is neglected. (c) As in case (b), but the torsional warping is included. Details of element division are shown on Fig. 4. The reinforcement cross sectional properties, u,~ deflection, GAl and +,,,* global rotations under the load and lon~tudinal membrane stress at the fixed end are summarised in Table 2. The essential difference observed between cases (a) and (b) is

Table 2. Cross sectional properties, deflection, rotations and stresses for the cantilever of Fig. 4 (Example 2)

A (mm2) wm”) &(mm’) &(mm4) GS (mm’) ~~~(mm2) &(mm*) qnmZ) l- (nM aA3(mm) rLAlx 104 +A2x t@ u&W a,(MPa)

a

b

c

2800 4,493,333 2,830,476 7,135,887 71.44 0 1100 397 4,227,164,540 - 1.2536 23.937 10.326 5.008 0.297

3800 5326,666 4,147,456 8,977,101, 61.47 0 1500 724 5,749,211,100 - 1.0272 f 1.143 8.454 4.189 1.259

3800 $326,666 4,147,456 8,977,lOl 61.47 3439.8 1.500 724 5,‘749,211,100 - 1.0126 9.000 8.429 4.084 1.287

308

GABOR M. VijRijs

//GbZ++e

’

CasFr

IFig. 4. Excentric concentrated load F = 1000 N on a reinforced cantilever: E = 210,OOO.O MPa, Y = 0.3 (Example 2).

obviously due to a wide overlap of the shell and the beam. Results of cases (b) and (c) are much closer, indicating the small influence of the beam warping on the whole structure.

6.

7. REFERENCES

1. D. N. Buragohain, S. B. Agrawal and R. S. Ayyar, A matching superparametric beam element for shell beam systems. Comput. Sfruct. 9, 175-182 (1978). 2. D. Bushnell, Computerized analysis of shell-governing equations. Comput. Sfruct. 18, 471-536 (1984). 3. J. R. G’Lery and I. Harari, Finite element analysis of stiffened plates. Compuf. Strucf. 21, 973-984 (1985). 4. G. H. Ferguson and R. D. Clark, A variable thickness curved beam and shell stiffening element with shear deformations. Int. J. Numer. Meth. Engng 14, 581-592 (1979). 5. J. Jirousek, A family of variable section curved beam

8.

9.

10. 11.

and thick-shell or membrane-stiffening isoparametric elements. Inr. .I. Numer. h4eth. Engng 17, 171-186 (1981). F. Kiss and J. Kiss, The new family of beam elements for ASKA. Proc. Int. FEM Congress, Baden-Baden (Edited by IKOSS Gmbh, Stuttgart), pp. 381-446 (1986). G. Wempner, Mechanics of Solids with Application to Thin Bodies. Sijthoff & Noordhoff, Rockville, MD (1981). G. M. Voriis, A variational principle in torsion problems of composite rods. Periodica Polytechnica 23, 367-376 (1979). S. Krenk and 0. Gunneskov, A triangulation procedjure for elastic cross sections with moderate wall thickness. Comput. Struct. 24, 1-12 (1986). D. Jangblad, QUAD9, a general shell element implemented in ASKA. SAAB-SCANIA report TKHR-3420 (1984). K. S. Surana, Isoparametric elements for cross sectional properties and stress analysis of beams. Inf. J. Numer. Meth. Engng 14, 475497 (1979).