Rigid Body Concept for Geometric Nonlinear Analysis of 3D Frames, Plates and Shells Based on the Updated Lagrangian Formulation

(1)

Rigid body concept for geometric nonlinear analysis of 3D frames,

plates and shells based on the updated Lagrangian formulation

Y.B. Yang

*

_{, S.P. Lin, C.S. Chen}

Department of Civil Engineering, National Taiwan University, Taipei 10617, Taiwan Received 26 October 2005; received in revised form 18 July 2006; accepted 26 July 2006

Abstract

In the nonlinear analysis of elastic structures, the displacement increments generated at each incremental step can be decomposed into two components as the rigid displacements and natural deformations. Based on the updated Lagrangian (UL) formulation, the geometric

stiﬀness matrix [kg] is derived for a 3D rigid beam element from the virtual work equation using a rigid displacement ﬁeld. Further, by

treating the three-node triangular plate element (TPE) as the composition of three rigid beams lying along the three sides, the [kg] matrix

for the TPE can be assembled from those of the rigid beams. The idea for the UL-type incremental-iterative nonlinear analysis is that if the rigid rotation effects are fully taken into account at each stage of analysis, then the remaining effects of natural deformations can be treated using the small-deformation linearized theory. The present approach is featured by the fact that the formulation is simple, the expressions are explicit, and all kinds of actions are considered in the stiffness matrices. The robustness of the procedure is demonstrated in the solution of several benchmark problems involving the postbuckling response.

Keywords: Beam; Geometric nonlinear analysis; Plate; Postbuckling; Rigid element; Shell

1. Introduction

The nonlinear analysis of elastic structures is usually conducted in an incremental-iterative way based on the three conﬁgurations: the initial conﬁguration C0, last

calcu-lated conﬁguration C1, and current deformed conﬁguration

C2, as indicated inFig. 1. In a step-by-step nonlinear

anal-ysis, we are interested in the behavior of the structure dur-ing the incremental step from C1to C2. The deformations

occurring within each incremental step are assumed to be small, but the displacements accumulated for all incremen-tal steps can be very large. The concept to be presented herein is based on the updated Lagrangian (UL) formula-tion, in that all quantities of the structure are expressed with reference to the last conﬁguration C1.

The displacement increments generated at each incre-mental step of an elastic nonlinear analysis can be com-posed into two components as the rigid displacements and natural deformations[1,2]. For most structures encountered in practice, the rigid component constitutes a much larger portion of the displacement increments at each incremental step with respect to the deformational component. For a UL-type incremental-iterative analysis, the idea is that if the rigid rotation eﬀects for elements with initial forces (or stresses) are fully taken into account at each stage of analysis, then the remaining eﬀects of natural deformations can be treated using the small-deformation linearized theory.

Concerning the incremental-iterative procedure, distinc-tion should be made between the predictor and corrector stages [3,4]. The predictor relates to solution of the dis-placement increments {U} of the structure for given load increments {P} based on the structural equation [K]{U} = {P}. This stage determines the trial direction of iteration of the structure in the load–deﬂection space and 0045-7825/$ - see front matter 2006 Elsevier B.V. All rights reserved.

doi:10.1016/j.cma.2006.07.013

*

Corresponding author. Tel.: +886 2 3366 4245; fax: +886 2 2362 2975. E-mail address:ybyang@ntu.edu.tw(Y.B. Yang).

(2)

thus affects the number of iterations or speed of conver-gence. For this reason, the stiffness matrix [K] used in the structural equation need not be exact, but should be kept rigid-body qualified to avoid convergence to incorrect directions. In the UL formulation, the corrector refers to recovery of the force increments {2f} at C2from the

dis-placement increments {u} made available through the structural displacement increments {U}, and the superim-position of these force increments with the initial nodal forcesf1

1fg following the rigid body rule[5,6]for obtaining

the total element forcesf2

2fg at C2.

In this paper, a rigid-body qualiﬁed geometric stiﬀness matrix [kg] will be derived for the 3D beam element from

the virtual work equation by assuming the displacement ﬁeld to be of the rigid type. Such an element is referred to as the rigid element. To the knowledge of the authors, no similar elements were presented by other scholars to explicitly accommodate the rigid behaviors of structures. For the 3D beam, the initial surface tractions may generate some moment terms upon 3D rotations during the incre-mental step from C1 to C2, commonly known as the

moments induced by the semitangential torques and quasi-tangential bending moments[1,7]. Naturally, all such terms should be included in the virtual work formulation for the 3D beam. The other issue to be considered for the space frames is the equilibrium of angled joints in the rotated conﬁguration C2, rather than in C1. Based on such

a consideration, only the symmetric part of the geometric stiﬀness matrix of each element has to be retained in the structural stiﬀness matrix, as the antisymmetric parts of all the elements meeting at the same joint cancel out with each other[6,7].

As for the analysis of plate/shell problems, a triangular plate element (TPE) with three translational and three rota-tional degrees of freedom (DOFs) at each of the three tip nodes will be considered, for its compatibility with the 12-DOF beam element derived above. Since the rigid body behavior of each finite element is solely determined by its external shape or nodal DOFs, the geometric stiffness matrix for the TPE is derived by treating the TPE as the composition of three rigid beams lying along the three sides. The geometric stiffness matrix so derived is explicit and capable of dealing with all kinds of in-plane and out-plane actions.

For a review of related works on geometric nonlinear analysis of structures, Ref.[8]may be consulted, in which a total of 122 papers were cited. The purpose of this paper is not to review any related works. Rather, efforts will be focused on application of the rigid body concept and deri-vation of rigid-body qualified geometric stiffness matrices [kg] for the 3D beam element and TPE. The elastic stiffness

matrices [ke] adopted are those readily available in the

literature, namely, the elastic stiﬀness matrix [ke] adopted

for the 3D beam element is the one commonly used [6,9], and the elastic stiﬀness matrix [ke] adopted for the TPE is

constructed as the composition of Cook’s plane hybrid ele-ment for membrane actions[10]and the hybrid stress model (HSM) of Batoz et al. for bending actions[11]. For the sake of brevity, repetition of relevant derivations is kept to the minimum.

2. Theory of three-dimensional beams

Before we proceed to derive the rigid element for the 3D beam, a summary of the theory for the 3D beam with bisymmetric solid cross-sections is ﬁrst given. The beam element considered has a total of 12 DOFs as shown in

Fig. 2, with x denoting the centroidal axis and (y, z) the two principal axes of the cross-section. Based on the UL formulation, the virtual work equation for a 3D beam at C2, but with reference to C1, can be expressed in a

linear-ized form as[6]: Z V ðE1exxd1exxþ 4G1exyd1exyþ 4G1exzd1exzÞ dV þ Z V ð1_s xxd1gxxþ 21sxyd1gxyþ 21sxzd1gxzÞ dV ¼2 1R 1 1R; ð1Þ

where E and G denote the elastic and shear modulus, respectively, V is the volume of the element, and the fac-tors of 4 and 2 are added to account for the symmetry of shear strains, i.e., 1exy=1eyx, 1exz=1ezx, 1gxy=1gyx, 1gxz=1gzx, (1sxx, 1sxy, 1sxz) are the initial (Cauchy) axial

and shear stresses, and d denotes the variation of the quantity following. The linear and nonlinear components of the strain increments can be expressed with reference to C1as 1 ya F 1 xa F 1 za F 1 xa M 1 ya M 1 za M 1 yb F 1 xb F 1 zb F 1 xb M 1 yb M 1 zb M y

z

x L a _b

Fig. 2. Three-dimensional beam element. P Q P Q P Q 0 C 1 C C2 1 X 3 X 2 X

(3)

1exx¼ oux ox; 1exy¼ 1 2 oux oy þ ouy ox ; 1exz¼ 1 2 oux oz þ ouz ox ; ð2Þ 1gxx¼ 1 2 oux ox 2 þ ouy ox 2 þ ouz ox 2 " # ; 1gxy¼ 1 2 oux oy oux ox þ ouy oy ouy ox þ ouz oy ouz ox ; 1gxz¼ 1 2 oux oz _ou x ox þ ouy oz _ou y ox þ ouz oz _ou z ox ; ð3Þ where (ux, uy, uz) are the displacements of a generic point

N(y, z) at cross-section x of the beam. By allowing the sur-face tractions 1

1tk and 2

1tk to exist only at the two end

sec-tions, the external virtual works 1

1R and

2

1R at C1and C2,

respectively, can be written as

1 1R¼ Z S ð1 1txduxþ11tyduyþ11tzduzÞ dS; 2 1R¼ Z S ð2 1txduxþ21tyduyþ21tzduzÞ dS; ð4a; bÞ

where S is the surface area of the element at C1.

Based on the Bernoulli–Euler hypothesis of plane sec-tions remaining plane and normal to the centroidal axis after deformation, the displacements (ux, uy, uz) of the

gen-eric point N can be related to the displacements (u, v, w) of the centroid of the same cross-section as

ux¼ u yv0 zw0; uy¼ v zhx; uz¼ w þ yhx; ð5Þ

where hx is the angle of twist. For the 3D beam, the

dis-placement increments {u} are

fug ¼ uh a va wa hxa hya hza ub vb wb hxb hyb hzbiT:

ð6Þ Correspondingly, the nodal forcesf1

1fg and f21fg acting on

the element at C1and C2are

f1 1fg ¼ h1_F xa 1Fya 1Fza 1Mxa 1Mya 1Mza 1_F xb 1Fyb 1Fzb 1Mxb 1Myb 1Mzbi T ; ð7Þ f2 1fg ¼ h2_F xa 2Fya 2Fza 2Mxa 2Mya 2Mza 2_F xb 2Fyb 2Fzb 2Mxb 2Myb 2Mzbi T ; ð8Þ

where the left subscript ‘‘1’’ is omitted from each term on the right hand side, i.e., 1_F

xa11Fxa; 1_M xa11Mxa; . . . ; 2_F xa21Fxa; 2_M

xa21Mxa; . . ., etc. All the nodal forces and

moments should be interpreted as the stress resultants, as will be given below.

2.1. Stress resultant deﬁnitions

Let us cut a beam at section x and consider its left por-tion as inFig. 3a. Also, let us attach two sets of coordinates to the centroid C of the cross-section. The ﬁrst set is the g–f coordinates embedded at the cross-section, with g and f denoting the two principal axes. The second set is assumed to have an origin ﬁxed at the centroid C, referred to as the

1_x–1_y–1_z_{axes at C}

1, withð1y;1zÞ coincident with (g, f) (see

Fig. 3a), and as the2_x–2_y–2_z_{axes at C}

2, withð2y;2zÞ parallel

toð1_y _and 1_{zÞ (see}_{Fig. 3}_b).

Based on the conditions of equilibrium, the forces acting on section x at Ck, with k = 1, 2, can be interpreted as

stress resultants over the cross-section of area A[6,7]:

k_F x¼ Z A k 1SxxdA; k_F y ¼ Z A k 1SxydA; k_F z¼ Z A k 1SxzdA; ð9Þ k_M x ¼ Z A ðk 1Sxz k_yk 1Sxy k_{zÞ dA;} k_M y ¼ Z A ðk 1Sxx k_zk 1Sxz k_{xÞ dA;} k_M z¼ Z A ðk 1Sxy k_xk 1Sxx k_{yÞ dA;} ð10Þ whereðk 1Sxx; k 1Sxy; k

1SxzÞ are the second Piola–Kirchhoﬀ

stres-ses acting at Ckwith reference to C1. For the element at C1,

(4)

these stresses reduce to the Cauchy stresses (1sxx,1sxy,1sxz),

and the embedded coordinates (g, f) coincide with the ð1_y;1_{zÞ axes, i.e.,}

1_x_{¼ 0;} 1_y_{g ¼ y;} 1_z_{f ¼ z:} _ð11Þ

It follows that the forces and moments in Eqs.(9) and (10)

reduce to those commonly known for the beam at C1. In

particular, the moments acting at C1are

1_M x¼ Z A ð1 sxzy1sxyzÞ dA; 1_M y¼ Z A 1_s xxzdA; 1Mz¼ Z A 1_s xxydA: ð12Þ

For the element to deform from C1 to C2, the

displace-ments (ux, uy, uz) of the generic point N at C1 were given

in Eq. (5). Consequently, the coordinates of point N at C2are (Fig. 3b)

2_x_{¼ yh}

zþ zhy; 2y¼ y zhx; 2z¼ z þ yhx: ð13Þ

Meanwhile, the second Piola–Kirchhoﬀ stresses can be expressed in an incremental form,

2 1Si¼

1

siþ Si ð14Þ

for i = xx, xy, xz, where Siare the stress increments. Since

the meanings for the forces (2Fx,2Fy,2Fz) at C2 are

self-explanatory, we shall focus on derivation of the moments (2Mx,2My,2Mz) at C2. By substituting Eqs. (13) and (14)

into Eq.(10), along with the stress resultant deﬁnitions in Eq. (12), the following can be obtained for the moments at C2: 2_M x¼ Z A ð2 1Sxz y 2 1Sxy zÞ dA; 2_M y¼ Z A 2 1Sxx z dA 1_M zhxþ 1 2 1_M xhz; 2_M z¼ Z A 2 1Sxx y dA þ 1_M yhx 1 2 1_M xhy; ð15Þ

where products of displacement and stress increments are neglected and the factor 1/2 is adopted for bisymmetric so-lid cross-sections. The expressions in Eq.(15) are derived based merely on the hypothesis of planar sections in Eq.

(5) and the stress resultant deﬁnitions in Eq. (10). Evi-dently, the torque1Mxshould be interpreted as the

semi-tangential moment, as indicated by the terms 1/21Mxhz

and 1/21

Mxhy, and the bending moments 1My and 1Mz

as the quasi-tangential moments, as indicated by the terms

1

Myhxand1Mzhx[1,6,7,12].

2.2. Strain energy term

With the linear strain components in Eq.(2)and the dis-placements in Eq.(5), one can obtain the strain energy dU in variational form from the ﬁrst integral of Eq.(1)as fol-lows[6]: dU ¼ Z V ðE1exxd1exxþ 4G1exyd1exyþ 4G1exzd1exzÞ dV ¼ Z L 0 ðEAu0_du0_{þ EI} yw00dw00þ EIzv00dv00þ GJ h0xdh 0 xÞ dx; ð16Þ where A is the cross-sectional area, L is the length, Iyand Iz

are the moments of inertia, and J is the torsional constant. From the strain energy term dU, a 12· 12 elastic matrix [ke] can be derived:

dU ¼ fdugT½kefug ð17Þ

as available elsewhere[6,9]. However, for the case of rigid displacements to be considered in the following, the strain energy term dU simply vanishes.

2.3. Virtual work of initial stresses

By substituting the nonlinear strain components in Eq.

(3) into the second integral of Eq. (1), along with the dis-placements in Eq.(5) and the deﬁnitions for forces in Eq.

(9) (for k = 1) and moments in Eq. (12), the virtual work dV of the initial stresses can be derived as[6]

dV ¼ Z V ð1_s xxd1gxxþ 21sxyd1gxyþ 21sxzd1gxzÞ dV ¼1 2 ZL 0 1_F xdðu02þ v02þ w02Þ þ1Fx Iy Adw 002_þIz Adv 002 þ1_F x Iyþ Iz A dh 02 x dx þ ZL 0 ½1_M zdðw0h 0 xÞ 1_M ydðv0h 0 xÞ 1_M ydðu0w00Þ þ1Mzdðu0v00Þ dx þ ZL 0 ½1_F ydðw0hx u0v0Þ 1Fzdðv0hxþ u0w0Þ dx þ1 2 ZL 0 ½1_M xdðv00w0Þ 1Mxdðw0v0Þ dx; ð18Þ

where higher order terms have been neglected.

2.4. External virtual work terms

By treating the surface tractions at C1as Cauchy

stres-ses, i.e.,1 1tx 1_s xx,11ty 1_s xy,11tz 1_s

xz, and by using the

dis-placements in Eq.(5) and the deﬁnitions for forces in Eq.

(9)(for k = 1) and moments in Eq.(12), the external virtual work1

1R at C1, Eq.(4a), can be written as 1

1R¼ fdug T

f1

1fg; ð19Þ

where the displacements {u} and initial forces f1

1fg were

deﬁned in Eqs. (6) and (7).

Similarly, by treating the surface tractions at C2as the

second Piola–Kirchhoﬀ stresses, i.e., 2 1tx 2 1Sxx, 2 1ty 2 1Sxy, 2 1tz 2

1Sxz, expressing the stresses in incremental form as

in Eq. (14), and using the displacement ﬁeld in Eq. (5)

and the deﬁnitions for forces in Eq. (9) (for k = 2) and moments in Eq.(12), one derives from Eq.(4b)the external virtual work2

(5)

2 1R¼ fdug T f2 1fg þ 1_M zhx 1 2 1_M xhz dhyþ 1Myhxþ 1 2 1_M xhy dhz L 0 ; ð20Þ where f2

1fg was deﬁned in Eq. (8). Clearly, the semi- and

quasi-tangential properties of the torque 1Mxand bending

moments 1My, 1Mz, respectively, are maintained in this

expression.

3. Geometric stiﬀness matrix for the rigid beam element Now, we shall derive the geometric stiﬀness matrix for the 3D rigid beam. For an element subjected only to nodal loads, the bending moments can be related to those at the two ends as 1_M y¼ 1Myað1 x=LÞ þ1Mybðx=LÞ; 1_M z¼ 1Mzað1 x=LÞ þ1Mzbðx=LÞ: ð21Þ

Further, based on the conditions of equilibrium, the fol-lowing can be written:

1_F x¼ 1Fxa¼1Fxb; 1Fy¼ 1Fya¼1Fyb¼ 1 Lð 1_M zaþ1MzbÞ; ð22Þ 1_F z¼ 1Fza¼1Fzb¼ 1 Lð 1_M yaþ1MybÞ; 1Mx¼ 1Mxa¼1Mxb: ð23Þ For a rigid rotation, the axial and angular displacements can be expressed as

u¼ ð1 x=LÞuaþ ðx=LÞub; hx ¼ ð1 x=LÞhxaþ ðx=LÞhxb

ð24Þ subject to the constraints:

ua¼ ub; hxa¼ hxb: ð25a; bÞ

The lateral and rotational displacements can be expressed as

v¼ ð1 x=LÞvaþ ðx=LÞvb; w¼ ð1 x=LÞwaþ ðx=LÞwb;

ð26Þ hy¼ ð1 x=LÞhyaþ ðx=LÞhyb; hz¼ ð1 x=LÞhzaþ ðx=LÞhzb

ð27Þ subject to the constraints:

hya¼ hyb; hza¼ hzb: ð28Þ

For the rigid beam, the strain energy vanishes, i.e., dU = 0. Consequently, the virtual work equation in Eq.(1), along with the aid of Eqs. (16), (18)–(20), reduces to the following: 1 2 Z L 0 1_F xdðu02þ v02þ w02Þ þ1Fx Iy Adw 002_þIz Adv 002 þ1_F x Iyþ Iz A dh 02 x dx þ Z L 0 ½1_M zdðw0h0xÞ 1_M ydðv0h0xÞ 1_M ydðu0w00Þ þ1Mzdðu0v00Þdx þ Z L 0 ½1_F ydðw0hx u0v0Þ 1Fzdðv0hxþ u0w0Þdx þ1 2 Z L 0 ½1_M xdðv00w0Þ 1Mxdðw00v0Þdx ¼ fdugTðf2 1fg f 1 1fgÞ þ 1_M zhx 1 2 1_M xhz dhyþ 1Myhxþ 1 2 1_M xhy dhz L 0 : ð29Þ

Since this is a problem with constraints, the following points must be considered in deriving each of the terms involved in Eq.(29) using the rigid displacement ﬁeld: (1) The following variational terms are equal to zero: du02₌

0, dh02_x ¼ 0, dv002_{= 0, dw}002_{= 0, d(u}0_w00_{) = 0, d(u}0_v00_{) = 0, as}

a direct consequence of the rigid displacement ﬁeld. (2) The variation of a quantity that is zero in the rigid displace-ment ﬁeld need not be zero. For instance, u0_{= (u}

b ua)/

L = 0 according to Eq. (25a), but du0_{= (du}

b dua)/L is

not equal to zero since dub5dua. (3) Even though some

of the variational terms in Eq.(29)are known to be zero, such as u0_dv0_{, u}0_dw0_{, h}0

xdv0, h 0

xdw0, they are kept in the

der-ivation for the sake of symmetry of the related entries in the stiﬀness matrix. (4) The second derivative terms v00

and w00 _{should be treated as h}0

z and h

0

y, respectively, and

interpolated by Eq. (27) in order to capture the property of rotations.

With the above considerations in mind, we may substi-tute the force expressions in Eqs. (21)–(23) and the dis-placement ﬁelds in Eqs. (24), (26) and (27) into the virtual work equation for the rigid beam in Eq. (29). By taking into account the arbitrary nature of the variations, we can derive the stiﬀness equation for the rigid beam ele-ment as ½kg beam fug ¼ f2 1fg f 1 1fg; ð30Þ

where the geometric stiﬀness matrix is

½kg beam ¼ ½g ½ha ½g ½hb ½ha T ½ia ½ha T ½0 ½g ½ha ½g ½hb ½hb T ½0 ½hb T ½ib 2 6 6 6 6 6 4 3 7 7 7 7 7 5 : ð31Þ

Here, [0] is a 3· 3 null matrix containing all zero entries, and ½g ¼ 0 1_F yb=L 1Fzb=L 1_F yb=L 1Fxb=L 0 1_F zb=L 0 1Fxb=L 2 6 6 4 3 7 7 5; ð32Þ

(6)

½hb ¼ 0 0 0 1_M yb=L 1Mxb=2L 0 1_M zb=L 0 1Mxb=2L 2 6 4 3 7 5; ð33Þ ½ib ¼ 0 0 0 1_M zb 0 1Mxb=2 1_M yb 1Mxb=2 0 2 6 4 3 7 5: ð34Þ

The submatrices [ha] and [ia] associated with end a can be

obtained by replacing the terms (1Mxb,1Myb,1Mzb) in the

submatrices [hb] and [ib] by the terms (1Mxa,1Mya,1Mza),

respectively. The matrix as presented in Eq.(31) is asym-metric, due to the asymmetry of the submatrices [ia] and

[ib] originating from the boundary terms in Eq. (29).

The geometric stiﬀness matrix derived in Eq.(31)is the same as that derived from the incrementally small-defor-mation theory [13]. It is qualiﬁed by the rigid body rule

[5,6], in the sense that when an element with the initial nodal forcesf1

1fg is subjected to a rigid rotation {u}r, the

total forcesf2

1fg obtained from Eq.(30)will have a

magni-tude equal to that of the initial nodal forcesf1

1fg, but with

the directions of the acting forces rotated by an angle equal to the rigid rotation. Of interest to note is that the eﬀects of all kinds of actions, i.e., the axial forces, shear forces, bend-ing moments, and torques, undergobend-ing the rigid rotations are considered in the [kg]beam matrix. In particular, this

matrix reduces to the one for the truss element acted upon by a pair of axial forces undergoing the rigid rotations. 4. Geometric stiﬀness matrix for the triangular plate element (TPE)

Let us consider the triangular plate element (TPE) shown inFig. 4, which has three translational and three rotational DOFs at each of the three tip nodes [2]. We choose this element simply because it is compatible with the 12-DOF beam element derived above.

In Section3, the geometric stiﬀness matrix for the rigid beam has been explicitly given in terms of the nodal forces and element length. Thus, as far as the rigid body behavior

of an element is concerned, only the initial forces acting on the element and the external shape of the element need to be considered. The elastic properties that are essential to the deformation of the element, such as Young’s modulus, cross-sectional area and moments of inertia, can be totally ignored. Based on such an idea, the rigid behavior of the TPE can be simulated as if it is composed of three rigid beams lying along the three sides of the element, as shown in Fig. 5. It is in this sense that the geometric stiﬀness matrix for the rigid TPE will be derived.

Fig. 6shows the nodal forces acting on each of the three beam elements, named as beam 12, beam 23 and beam 31. In this ﬁgure, 1_Fij

k and 1M

ij

k denote the nodal force and

nodal moment, respectively, acting on beam ij at the C1

conﬁguration, with the right subscript k denotes the direc-tion and the nodal point at which the force or moment is acting. In the following, we shall show how to determine the nodal forces acting on beam ij:

(1) Equations of equilibrium for each node of the TPE: There are three forces and three moments acting at each node of the TPE. The following is the vector of forcesf1

1fg

for the element at the C1conﬁguration:

f1 1fg ¼ 1_F x1 1Fy1 1Fz1 1Mx1 1My1 1Mz1 1Fx2 1Fy2 1Fz2 1_M x2 1My2 1Mz2 1Fx3 1Fy3 1Fz3 1Mx3 1My3 1Mz3 T : ð35Þ 1 x θ 1 y θ 1 z θ 2 x θ 2 y θ 2 z θ 3 x θ 3 y θ 3 z θ 1 u 2 u 3 u 1 v 1 w (X , Y1 1)

(

X , Y2 2

)

2 v 2 w

(

X , Y3 3

)

3 w 3 v 1 X 1 Y 1_Z C(0, 0) :Centriod C

Fig. 4. Three-node triangular plate element.

1 1 1 2 2 2 3 3 3 a a a b b b

Fig. 5. TPE treated as the composition of three rigid beams.

1 23 ya F 1 23 xa F 1 23 za F 1 23 xa M 1 23 ya M 1 23 za M 1 23 yb F 1 23 xb F 1 23 zb F 1 23 xb M 1 23 yb M 1 23 zb M 1 12 ya F 1 12 xa F 1 12 za F 1 12 xa M 1 12 ya M 1 12 za M 1 12 yb F 1 12 xb F 1 12 zb F 1 12 xb M 1 12 yb M 1 12 zb M 1 31 ya F 1 31 xa F 1 31 za F 1 31 xa M 1 31 ya M 1 31 za M 1 31 yb F 1 31 xb F 1 31 zb F 1 31 xb M 1 31 yb M 1 31 zb M

1 ₁

2

3

(7)

Considering the equilibrium of the forces and moments acting at each node, we can write

(a) Node 1 1_F12 xaþ1F 31 xb¼1Fx1; 1F12yaþ1F 31 yb¼1Fy1; 1F12zaþ1F 31 zb¼1Fz1; 1_M12 xaþ 1_M31 xb¼ 1_M x1; 1M12yaþ 1_M31 yb¼ 1_M y1; 1M12zaþ 1_M31 zb¼ 1_M z1: ð36Þ (b) Node 2 1_F23 xaþ1F 12 xb¼1Fx2; 1F23yaþ1F 12 yb¼1Fy2; 1F23zaþ1F 12 zb¼1Fz2; 1_M23 xaþ 1_M12 xb¼ 1_M x2; 1M23yaþ 1_M12 yb¼ 1_M y2; 1M23zaþ 1_M12 zb¼ 1_M z2: ð37Þ (c) Node 3 1_F31 xaþ1F 23 xb¼1Fx3; 1F31yaþ1F 23 yb¼1Fy3; 1F31zaþ1F 23 zb¼1Fz3; 1_M31 xaþ 1_M23 xb¼ 1_M x3; 1M31yaþ 1_M23 yb¼ 1_M y3; 1M31zaþ 1_M23 zb¼ 1_M z3: ð38Þ (2) Equations of equilibrium of beam ij:

All the forces and moments acting on beam ij (with ij = 12, 23 or 31) must satisfy the conditions of equilibrium for the beam, as given below:

(a) Equilibrium of forces:

1_Fij xaþ 1_Fij xb¼ 0; 1_Fij yaþ 1_Fij yb¼ 0; 1_Fij zaþ 1_Fij zb¼ 0: ð39Þ (b) Equilibrium of moments: Yi1Fijzaþ Yj1Fijzb þ1 Mij_xaþ1 Mij_xb¼ 0; Xi1Fijzaþ Xj1Fijzb þ1_Mij yaþ 1_Mij yb¼ 0; Yi1Fijxaþ Yj1F ij xb þ Xi1Fijyaþ Xj1F ij yb þ1 Mij_zaþ1 Mij_zb¼ 0: ð40Þ for ij = 12, 23 or 31, where (Xi, Yi) are the

coordi-nates of node i.

(3) Equations of equilibrium for the TPE:

The conditions of equilibrium must be obeyed by the TPE as a whole, as given below:

(a) Equilibrium of forces: X3 k¼1 1_F xk¼ 0; X3 k¼1 1_F yk ¼ 0; X3 k¼1 1_F zk¼ 0: ð41Þ (b) Equilibrium of moments: X3 k¼1 Yk1Fzkþ1Mxk ¼ 0; X 3 k¼1 Xk1Fzkþ1Myk ¼ 0; X3 k¼1 Yk1Fxkþ Xk1Fykþ1Mzk ¼ 0: ð42Þ

Now, we can solve Eqs.(36)–(42)to obtain the nodal forces for beam ij as 1_Fij xa¼ 1 3 1_F xijþ fx; 1Fijya¼ 1 3 1_F yijþ fy; 1Fijza¼ 1 3 1_F zijþ fz; 1_Fij xb¼ 1F ij xa; 1F ij yb¼ 1F ij ya; 1F ij zb¼ 1F ij za ð43Þ and the nodal moments as

1_Mij xa¼ 1 3 1_M xijþ 1 3Yij 1_F zjþ ðmx YifzÞ; 1_Mij ya¼ 1 3 1_M yij 1 3Xij 1_F zjþ ðmyþ XifzÞ; 1_Mij za¼ 1 3 1_M zij 1 3 Yij 1_F xj Xij1Fyj þ ½mzþ ðYifx XifyÞ; 1_Mij xb¼ 1_Mij xa Yij1Fijza; 1_Mij yb¼ 1_Mij yaþ Xij1Fijza; 1_Mij zb¼ 1_Mij zaþ Yij1Fijxa Xij1Fijya; ð44Þ where a variable with right subscript ij denotes the diﬀer-ence of two quantities, e.g., 1Mxij=1Mxi1Mxj, Xij=

Xi Xj, and fx, fy, fz, mx, myand mz are the constants to

be determined. One feature with the rigid element is that there are more unknowns than the equations that can be utilized. For instance, the total number of equations given in Eqs. (36)–(42) is 42. By deducting the 12 dependent equations, the total number of equations that can be used is 30, but the total number of unknown forces in Eqs.(43) and (44)is 36, with 12 for each of the three beam elements. This means that the force distribution of a rigid element is not unique, unlike that of an elastic element, and that six more conditions should be made in order to obtain a solu-tion. In this paper, we simply let the six variables fx, fy, fz,

mx, myand mzequal to zero.

All the nodal forces and moments solved and given in Eqs. (43) and (44) have been expressed with reference to the coordinates of the TPE, as indicated in Fig. 6. They have to be transformed into the local coordinates of each element in order to compute the geometric stiﬀness matrix [kg]beamof the element. Let 11f

ij

n o

denote the nodal forces of element ij in the global coordinates of the TPE and

1 1f

ij

n o

in the local coordinates of element ij,

1 1f ij n o ¼ h1_Fij xa 1_Fij ya 1_Fij za 1_Mij xa 1_Mij ya 1_Mij za 1_Fij xb 1_Fij yb 1_Fij zb 1_Mij xb 1_Mij yb 1_Mij zbi T ; ð45Þ 1 1f ij n o ¼ h1_Fij xa 1_Fij ya 1_Fij za 1_Mij xa 1_Mij ya 1_Mij za 1_Fij xb 1_Fij yb 1_Fij zb 1_Mij xb 1_Mij yb 1_Mij zbi T : ð46Þ The following is the transformation between the two sets of forces for element ij:

1 1f ij n o ¼ ½T 1 1f ij n o ; ð47Þ

(8)

where the transformation matrix [T] is ½T ¼ ½R ½0 ½0 ½0 ½0 ½R ½0 ½0 ½0 ½0 ½R ½0 ½0 ½0 ½0 ½R 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 ; ½R ¼ Xij Lij Yij Lij 0 Yij Lij Xij Lij 0 0 0 1 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 ð48Þ with Lij denoting the length of element ij. By substituting

the nodal forces 1 1f

ij

n o

into Eq.(31), the geometric stiﬀness matrix [kg]beam_ijfor element ij can be obtained. Finally, by

assembling the geometric stiffness matrices for the three ele-ments following the standard finite element procedure, the geometric stiffness matrix [kg]TPE for the TPE can be

obtained ½kgTPE¼

X

ij¼12;23;31

½T T½kgbeam ij½T ; ð49Þ

which has been explicitly given inAppendix.

It is easy to verify that the geometric stiﬀness matrix [kg]TPEderived above for the TPE passes the rigid body test

[5,6], in that the resulting forcesf2

1fg computed from Eq.

(30) for the TPE undergoing a rigid rotation maintain a magnitude equal to that of the initial forcesf1

1fg, but with

the directions of action rotated by an angle equal to the rigid rotation.

5. Joint equilibrium condition in the rotated position Both the geometric stiﬀness matrix [kg]beamfor the rigid

beam, as given in Eq.(31), and the [kg]TPEmatrix for the

TPE, as given in the Appendix, are asymmetric, as indi-cated by the submatrices [ik] and [Ik], respectively, which

relate to the behavior of nodal moments undergoing 3D rotations[6,7]. Such a property of asymmetry is restricted to the element level, but not the structure level, which adds little extra eﬀort to nonlinear analysis. To explain such a fact, let us decompose the geometric stiﬀness matrices for the rigid beam and TPE, as given in Eq. (31) and the

Appendix, respectively, into the symmetric and anti-sym-metric parts: ½kg beam ¼ ½kg beam sym þ ½kg beam anti-sym; ð50Þ ½kg TPE ¼ ½kg TPE sym þ ½kg TPE anti-sym; ð51Þ

where the anti-symmetric parts are

½kg beam anti-sym¼ ½0 ½0 ½0 ½0 ½0 ½Aa ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½Ab 2 6 6 6 6 6 4 3 7 7 7 7 7 5 ; ð52Þ ½kg TPE anti-sym¼ ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½A1 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½A2 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½0 ½A3 2 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 5 ð53Þ

as indicated by the submatrix [A]k, which is a function of

the nodal moments,

½Ak ¼ 1 2 0 1_M zk 1Myk 1_M zk 0 1Mxk 1_M yk 1Mxk 0 2 6 4 3 7 5 ð54Þ

with the left subscript k = a, b for the beam element, and k = 1, 2, and 3 for the TPE.

In the literature[6,7], it has been demonstrated that the anti-symmetric parts of the geometric stiﬀness matrices of all beam elements meeting at the same joint will cancel out, once the conditions of equilibrium for the joint in the rotated conﬁguration C2are enforced. As a result, the

stiﬀness matrix assembled for the structure turns out to be symmetric. The same is also true for the plates and shells represented by the TPE, as will be proved below.

Using the symbol eijk to denote a permutation symbol

[14], the anti-symmetric [A] for each node of the TPE can be represented as follows: Aij ¼ 1 2eijk 1_M k ði; j; k ¼ 1; 2; 3Þ; ð55Þ

which implies that1M1=1Mx,1M2=1My, etc. Let [U]

de-note the transformation matrix from the local coordinates of an element to the global coordinates of the structure. Then,

½UT½U ¼ ½I; ð56Þ

where [I] is an identity matrix, or

UpiUpj¼ dij ð57Þ

with dijdenoting the Kronecker delta. Based on the

tenso-rial operations, the following is shown to be valid: eijkUpiUqjUrk¼ epqr ði; j; k; p; q; r ¼ 1; 2; 3Þ: ð58Þ

The nodal moment vector {M} in the element coordinates can be transformed to the structure coordinates as

fMg ¼ ½UfMg; ð59Þ

where fMg denotes the nodal moments in the structure coordinates, or

1_M

p¼ Upi1Mi ði; p ¼ 1; 2; 3Þ: ð60Þ

Similarly, the anti-symmetric matrix [A] can be trans-formed to the global coordinates as

(9)

or

Apq¼ UpiAijUqj ði; j; p; q ¼ 1; 2; 3Þ: ð62Þ

Substituting Eq.(55)into Eq.(62)yields Apq¼ UpiUqjAij ¼ 1 2eijkUpiUqj 1_M k ði; j; k; p; q; r ¼ 1; 2; 3Þ: ð63Þ By the following relations:

eijkUpiUqj ¼ eijtUpiUqjdtk¼ eijtUpiUqjUrtUrk ¼ epqrUrk: ð64Þ

Eq.(63)can be rewritten as follows: Apq¼ 1 2epqrUrk 1_M k¼ 1 2epqr 1_M r ðk; p; q; r ¼ 1; 2; 3Þ: ð65Þ

Next, let us consider a nodal point of the plate/shell structure where n TPE elements lying along diﬀerent direc-tions are connected. For the nodal point to be in equilib-rium in the rotated conﬁguration C2, the sum of moments

exerted by all the elements meeting at the node along the three global axes must be equal to zero,

Xn s¼1 1_MðsÞ r ¼ 0: ð66Þ Consequently, we have Xn s¼1 AðsÞ_pq ¼X n s¼1 1 2epqr 1_MðsÞ r ¼ 1 2epqr Xn s¼1 1_MðsÞ r ¼ 0 ðp; q; r ¼ 1; 2; 3Þ: ð67Þ

This implies that when the anti-symmetric matrices ½A, which are expressed in the global coordinates, are summed over all the TPE elements sharing a common node, the resulting stiffness matrix associated with this node is zero. Thus, we have proved that the total stiffness matrix for the structure in the rotated configuration C2is symmetric,

although the element stiﬀness matrices are known to be asymmetric.

6. Predictor and corrector for incremental-iterative analysis Only incremental-iterative nonlinear analysis of the UL type is of concern in this study. The idea is that if the rigid rotation effects are fully taken into account at each stage of analysis, then the remaining effects of natural deformations can be treated using the small-deformation linearized the-ory. In this regard, two major stages can be identified

[3,4]. The ﬁrst is the predictor or trial stage, which relates to solution of the structural displacements {U}, given the load increments {2P} {1

P}, as indicated by the following equation:

½KfU g ¼ f2_P_{g f}1_Pg; _ð68Þ

where {1P} and {2P} denote the applied loads acting on the structure at C1 and C2, respectively. Once the structural

displacement increments {U} are made available, the dis-placements increments {u} for each element can be

com-puted accordingly. It is known that the predictor aﬀects only the number of iterations or the speed of convergence

[3,4]. Therefore, the structural stiffness matrix [K] is al-lowed to be approximate, but to the limit that the direction of iterations is not misguided. In addition to the elastic stiff-ness matrix, this will require the use of geometric stiffstiff-ness matrices that are capable of simulating the rigid rotational behavior of the initially stressed elements, such as those presented herein for the 3D beam and TPE elements.

For the 3D beam element, both the elastic stiﬀness matrix [ke] available in Refs.[6,9] and the geometric

stiﬀ-ness matrix [kg] in Eq. (31) will be used in constructing

the structural stiﬀness matrix [K]. For the TPE, the elastic stiﬀness matrix [ke] is obtained as the composition of

Cook’s plane hybrid element for membrane actions [10]

and the HSM element of Batoz et al. for bending actions

[11]. The elastic stiﬀness matrices [ke] so obtained, plus

the geometric stiﬀness matrices [kg] TPE

given in the Appen-dix, will be used to assemble the structural stiﬀness matrix [K] for the plates and shell considered.

The corrector stage is concerned with the recovery of the element forces 2

2f

at C2with reference to C2. For

incre-mental analysis of the UL type, two contributions need to be considered. The ﬁrst is the initial nodal forces 1

1f

existing at C1and expressed with reference to C1.

Accord-ing to the rigid body rule, the initial nodal forcesf1 1fg will

be rotated by an angle equal to the rigid rotation (with no limit in magnitude of rotation) generated during the incre-mental step from C1to C2, while their magnitude remains

unchanged [5,6]. Thus, one can directly treat the initial nodal forcesf1

1fg as the forces acting at C2and expressed

with respect to C2. As the rigid rotation eﬀect is already

taken into account, the force increments {2f} can be

com-puted from the displacement increments {u} using the elastic stiﬀness matrix [ke] alone based on the

small-defor-mation linearized theory, i.e.,

f2fg ¼ ½kefug: ð69Þ

Summing the above two eﬀects yields the total element forces 2 2f at C2as 2 2f ¼ f1 1fg þ f2fg; ð70Þ

where the reference conﬁguration is C2.

As the structural and element displacements, i.e., {U} and {u}, are made available from the predictor stage in Eq.(68), one may compute the nodal coordinates and ele-ment orientations for the structure at the displaced conﬁg-uration C2. Here, the nodal rotations generated at each

incremental step are assumed to be so small that the law of commutativity applies. Then, by summing the element forces 2

2f

computed from Eq. (70)for each node of the structure at C2 and by comparing them with the total

applied loads {2P}, the unbalanced forces for the structure at C2can be obtained. Another iteration involving the

pre-dictor and corrector stages should be conducted to elimi-nate the unbalanced forces, should they be greater than preset tolerances.

(10)

In a UL-type incremental-iterative nonlinear analysis, the accuracy of the solution is governed primarily by the corrector; while the predictor can only aﬀect the speed of convergence or the number of iterations[3,6]. For the pur-pose of investigating the capacity of the geometric stiﬀness matrix derived in this paper, the following combinations of predictor and corrector will be tested in the analysis of the space frames:

• Predictor: The structural stiﬀness matrix [K] in Eq.(68)

is assembled using either of the following:

(P1) [ke] + [kg]beam, where [kg]beamis the one derived in

this paper for the rigid beam.

(P2) [ke] + [kg] + [ki], where [kg] is the (conventional)

geometric stiﬀness matrix and [ki] the induced

moment matrix available on pp. 361–362 of Ref.

[6].

(P3) [ke] matrix only.

• Corrector:

(C1) {2f} = [ke]{u} as proposed in Eq.(69).

(C2) {2f} = ([ke] + [kg] + [ki]){un}, where {un} denotes

the natural deformations of the element.

As for the plates and shells, the TRIC element developed by Argyris et al.[15]will be used to yield solutions for comparison. The following are the combinations to be used:

• Predictor: The structural stiﬀness matrix [K] in Eq.(68)

is assembled using either of the following:

(P1) [ke]TPE+ [kg]TPE, where [ke]TPEis composed of the

elastic stiﬀness matrices available in Refs.[10,11]; (P2) [ke]

TRIC

+ [kg] TRIC

, where both the stiﬀness matri-ces are given in Ref.[15].

(P3) [ke] TPE

matrix only. • Corrector:

(C1) {2f} = [ke]TPE{u}, as proposed in Eq. (69) of this

paper.

(C2) {2f} = [ke]TRIC{u}.

7. Numerical examples

The generalized displacement control (GDC) method[16]

is adopted for tracing the nonlinear load–deﬂection response of the problems studied. With the aid of the

gen-eralized stiffness parameter (GSP), this method can auto-matically adjust the load increment sizes to reflect the variation in stiffness of the structure, while reversing the direction of loading when passing a limit point. It can easily deal with the various critical points encountered in the postbuckling response of a structure. Since all the analyses are allowed to self adjust their load increments by the GDC method, the total computation time consumed by each scheme may not truly reflect its efficiency of computation. Example 1. Fig. 7 shows a single beam restrained against the rotations about the x and y axes at the left end and against the rotations about the y and z axes at the other end. The following data are adopted: length L= 100 mm, cross-sectional area A = 0.18 cm2, torsional constant J = 2.16 mm4, moments of inertia Iy= 0.54 mm4,

Iz= 1350 mm4, Young’s modulus E = 71,240 N/mm2,

and shear modulus G = 27,190 N/mm2. The beam is subjected to a moment of Mza at the left end, and a

perturbation of Mxb= 0.01Mzaat the right end (to trigger

the lateral buckling). Ten elements are used for the beam. The theoretical critical moment for the beam under uniform moment is Mcr¼ p ffiffiffiffiffiffiffiffiffiffiffiffiffiEIyGJ

p

=L¼ 1493:3 N mm

[6,7].

The numerical results obtained are plotted in Fig. 8, along with the computation time of each analysis listed in

Table 1. The peak load inFig. 8is 1463.0 N mm. As can be seen, the proposed scheme P1C1 can be used to yield solutions that are as accurate as the conventional one, though at the cost of longer computation time. It is inter-esting that the same problem can be solved using the elastic stiﬀness matrix [ke] alone in all stages of analysis, as

indi-cated by P3C1 in Table 1. Of course, this is achieved at the cost of much longer computation time.

Fig. 8. Moment–displacement curves for single beam: (a) Uxa, (b) Uzb.

a _x L y z xb M za M b

(11)

Example 2. Fig. 9 shows an angled frame under uniform bending (i.e., with Mza= Mzc) and a shear load Fzb(N) of

magnitude 5· 105Mza[1,6,7]. The length L is 240 mm and

other properties are the same as in Example 1. Due to symmetry of the frame, only the left half modeled by 10 elements is analyzed. The member ab is restrained against rotations about the x, y axes and translations along the y, z axes. The translation along the x-axis at node b is also restrained. The theoretical critical loads for the frame are: Mcr¼ ppffiffiffiffiffiffiffiffiffiffiffiffiffiEIyGJ=L¼ 622:2 N mm.

From the results plotted inFig. 10, one observes that the proposed scheme P1C1 can be reliably used to trace the postbuckling behavior of the frame, though at the expense of slightly longer computation time, compared with the conventional one (P2C2). Of interest is that using the [ke]

matrix alone in both the predictor and corrector, as indi-cated by P3C1, can arrive exactly at the same solution, but with much longer computation time (seeTable 2). Example 3. The right-angled frame subjected to an in-plane load Fybat the free end inFig. 11(a) was studied[1],

for which the critical load is Fyb,cr= 1.088 N. The same

material and geometry data as those used inExample 2are adopted. Each member of the frame is modeled by 10 beam elements. The entire frame was also modeled by 38 triangular elements as in Fig. 11b. An imperfection load Fzb equal to one thousandth of Fyb is applied at the free

end.

From the results plotted in Fig. 12, it is seen that the proposed scheme P1C1 using either the beam or TPE approach is capable of tracing the postbuckling response of the angled frame. Table 3 indicates that the time con-sumed by the present beam approach is slightly longer than the conventional approach. The same is also true for the TPE approach, compared with the TRIC approach. It should be noted that the number of iterations at each incre-mental step is about the same for both the TPE and TRIC

approaches. However, less time is consumed by the TRIC approach in computing the stiﬀness matrices, since it con-siders only the in-plane actions, while the TPE concon-siders all kinds of in-plane and out-of-plane actions.

No convergent solution was obtained by the P3C1 scheme using only the [ke] matrix, due to the fact that the

trial directions produced by the predictor based on the [ke] matrix alone deviates too much from the equilibrium

path in the region near the bifurcation point.

Example 4. Fig. 13shows a spherical shell subjected to the central point load P with all edges hinged and immovable

[17]. The data adopted herein are: half of side length a = 784.9 mm, thickness t = 99.45 mm, Young’s modulus radius of curvature R = 2540 mm, E = 68.95 N/mm2, and Poisson’s ratio v = 0.3. Because of symmetry, only one quarter of the shell is modeled as an 8· 8 mesh and analyzed.

The results obtained for the load–deﬂection response of the central point using diﬀerent combinations are plotted in

Fig. 14, along with the running time listed inTable 4. Basi-cally no distinction can be made between the results obtained by diﬀerent approaches. Compared with the Fig. 10. Moment–displacement curves for hinged angled frame: (a) Uxa,

(b) hza.

Table 2

Running time for analysis ofExample 2

Combination P1C1 P2C2 P3C1

CPU time (s) 116.91 109.25 1081.13

Table 1

Running time for analysis ofExample 1

Combination P1C1 P2C2 P3C1 CPU time (s) 55.99 38.02 474.33 a x L y z za M b zc M c zb F

(12)

TRIC approach, slightly longer computation time is required by the TPE approach due to inclusion of all kinds of in-plane and out-of-plane actions in the stiﬀness matri-ces by the TPE approach, though the number of iterations at each incremental step is about the same. Of interest is

yb F b zb F y x z 240 mm 240 mm a yb F b zb F y x z 30 mm h= 0.6 mm t= 225 mm 255 mm 255 mm

Fig. 11. Right angled frame with ﬁxed support: (a) beam elements, (b) plate elements.

Fig. 12. Load–displacement curves for angled frame: (a) beam approach, (b) plate approach.

Table 3

Running time for analysis ofExample 3 Approach Combination P1C1 P2C2 P3C1 Beam 37.59 36.16 Divergent Plate 1492.22 1211.25 Divergent P 2a 2a R 2 2540 mm 784.90 mm 99.45 mm 68.95 N/mm 0.3 R a t E v = = = = = t

(13)

that using the [ke]TPEstiﬀness matrix alone, as indicated by

P3C1, is good enough for obtaining the solution, though at the cost of extra computation time.

Example 5. A cylindrical shell subjected to a central load P on the top surface is shown inFig. 15 [18]. The longitudinal boundaries of the shell are hinged and immovable, whereas the curved edges are completely free. The following data are adopted: E = 3.10275 kN/mm2, R = 2540 mm, L = 254 mm, h = 12.7 mm, h = 0.1 rad, v = 0.3. Because of symmetry, only one quarter of the shell is considered and modeled by an 8· 8 mesh.

The results obtained for the central deﬂection of the shell are plotted in Fig. 16, along with the computation time of each analysis in Table 5. Again, no distinction can be made between the results obtained by diﬀerent approaches. The computation time consumed by the pres-ent TPE approach is slightly longer than the TRIC approach for the reasons stated in Example 3. The P3C1 approach with the [ke] matrix alone is able to trace the

entire load–deﬂection response, though at the cost of extra computation time.

Example 6. This example is identical toExample 5except that the thickness of the shell is halved to h = 6.35 mm

[18]. The results obtained using an 8· 8 mesh have been shown in Fig. 17, together with the computation time in Table 6. A comparison of Fig. 17 with Fig. 16

indicates that reducing the thickness of the shell by half has resulted not only in a drastic decline of the limit load capacity, but also in a shift of the response curve pattern to embrace a snap-back region. As can be seen from Table 6, slightly longer computation time is con-sumed by the present TPE approach, compared with the TRIC approach, for which the reason was given in

Example 3. Again, this problem can be solved using exclusively the [ke] matrix in each stage of analysis, as

Table 6

Running time (s) for analysis ofExample 6

Mesh Combination P1C1 P2C2 P3C1 8· 8 197.08 177.78 559.73 P 2L h R θ

Fig. 15. Hinged cylindrical shell.

Fig. 16. Central deﬂection of hinged cylindrical shell with h = 12.7 mm.

Table 5

Mesh Combination P1C1 P2C2 P3C1 8· 8 397.42 388.57 668.48 -0.4 -0.2 0 0.2 0.4 0.6 0.8 0 5 10 15 20 25 30 Central deflection (mm) Load P (kN) P1C1 (Present) P2C2 (TRIC) P3C1

Fig. 17. Central deﬂection of hinged cylindrical shell with h = 6.35 mm. Fig. 14. Load–deﬂection curve for hinged spherical shell.

Table 4

Mesh Combination

P1C1 P2C2 P3C1

(14)

indicated by P3C1, though much longer computation time is required.

8. Concluding remarks

The procedure presented in this paper is based on the idea that if the rigid rotation effects are fully taken into account at each stage of the incremental-iterative nonlinear analysis, then the remaining effects of natural deformations can be treated using the small-deformation linearized the-ory. Based on the UL formulation, two stages are consid-ered essential, i.e., the predictor stage for computing the structural displacements given the load increments, and the corrector stage for recovering the element forces. The former affects the number of iterations, but the latter deter-mines the accuracy of solution.

As for use in the predictor, the geometric stiﬀness matrix [kg] is derived for a 3D rigid beam element from the virtual

work equation by using a rigid displacement ﬁeld, and the geometric stiﬀness matrix [kg] for the rigid triangular plate

element (TPE) is assembled from those for the three rigid beams lying along the three sides of the element. One advantage with the present approach is that the geometric stiﬀness matrices [kg] for both the 3D rigid beam and rigid

TPE are explicitly given and all kinds of in-plane and out-of-plane actions are taken into account. Using the rigid ele-ment concept, the work of derivation is greatly reduced for both elements, as only rigid displacement ﬁelds are required.

As for the corrector, the rigid body rule is used to update the initial nodal forcesf1

1fg existing at C1with no limit on

the magnitude of rigid rotations. With this, the force incre-ments generated at each incremental step are computed using only the elastic stiﬀness matrix [ke] based on the

small-deformation linearized theory.

The robustness of the proposed combination of predic-tor and correcpredic-tor has been demonstrated in the solution of a number of frame and shell problems involving post-buckling responses. A slightly longer computation time is required for the TPE approach due to consideration of all kinds of actions, compared with the TRIC approach that considers only in-plane actions, based on the frame-work of the GDC method for self adjustment of load incre-ments. For problems with no abrupt change in the slope of the load–deﬂection curves, the entire postbuckling response can be traced by using only the elastic stiﬀness matrix [ke]

in each stage of analysis, though at the cost of extra com-putation time.

Acknowledgements

The research works reported herein have been con-ducted as the result of a series of research projects granted by the ROC National Science Council to the senior author, including NSC 80-0410-E002-18, NSC 80-0410-E002-59, and NSC 81-0410-E002-05.

Appendix. Geometric stiﬀness matrix for the TPE

The geometric stiﬀness matrix derived in this paper for the TPE is

½kg TPE

¼

½E1 ½F1 ½E12 ½G12 ½E31 ½H31

½F1 T ½I1 ½H12 T ½0 ½G31 T ½0 ½E12 T ½H12 ½E2 ½F2 ½E23 ½G23 ½G12 T ½0 ½F2 T ½I2 ½H23 T ½0 ½E31 T ½G31 ½E23 T ½H23 ½E3 ½F3 ½H31 T ½0 ½G23 T ½0 ½F3 T ½I3 2 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 5 ; where ½Eij ¼ aij _cij _dij cij _bij _eij dij eij _fij 2 6 4 3 7 5; ½Gij ¼ nij _pij ₀ qij _oij ₀ rij _sij _mij 2 6 4 3 7 5; ½Hij ¼ gij _iij ₀ jij _hij 0 kij lij mij 2 6 4 3 7 5; ½I1 ¼ t12 1 þ t312 u121 þ u312 w121 þ w312 v12 1 þ v 31 2 t 12 1 þ t 31 2 y12 1 þ y 31 2 x12 1 þ x312 z121 þ z312 0 2 6 4 3 7 5; ½I2 ¼ t23 1 þ t 12 2 u 23 1 þ u 12 2 w 23 1 þ w 12 2 v23 1 þ v122 t231 þ t122 y23 1 þ y122 x23 1 þ x 12 2 z 23 1 þ z 12 2 0 2 6 4 3 7 5; ½I3 ¼ t31 1 þ t232 u311 þ u232 w311 þ w232 v31 1 þ v 23 2 t 31 1 þ t 23 2 y31 1 þ y 23 2 x31 1 þ x232 z311 þ z232 0 2 6 4 3 7 5; ½E1 ¼ ½E12 þ ½E31; ½E2 ¼ ½E23 þ ½E12;

½E3 ¼ ½E31 þ ½E23;

½F1 ¼ ½H12 ½G31; ½F2 ¼ ½H23 ½G12;

½F3 ¼ ½H31 ½G23:

In the above expressions, the parameters are given as follows: aij ¼ 1 L4_ij XijY 2 ijF ij xa X 2 ijYijFijya L 2 ijYijFijya ; bij¼ 1 L4 ij XijY2ijF ij xa X 2 ijYijFijyaþ L 2 ijXijFijxa ; cij_¼ 1 L4 ij Y3_ijFij xaþ X 3 ijF ij ya ; dij¼ 1 L2 ij XijFijza; eij ¼ 1 L2_ijYijF ij za; fij¼ 1 L2 ij XijFijxaþ YijFijya ; gij_¼ 1 2L4 ij XijY2ijM ij xa X 2 ijYijMijya L 2 ijYijMijya

(15)

hij ¼ 1 2L4 ij XijY2ijM ij xa X 2 ijYijMijyaþ L 2 ijXijMijxa ; iij¼ 1 2L4 ij Y3_ijMij_xa XijY2ijM ij yaþ L 2 ijYijMijxa ; jij_¼ 1 2L4 ij X2_ijYijMijxa X 3 ijM ij ya L 2 ijXijMijya ; kij ¼ 1 L2_ijXijM ij za; l ij ¼ 1 L2_ijYijM ij za; m ij ¼ 1 2L2 ij XijMijxaþ YijMijya ; nij¼1 L2 ij XijYijFijzaþ 1 2L4 ij XijY2ijM ij xa X 2 ijYijMijya L 2 ijYijMijya ; oij_¼1 L2 ij XijYijFijza 1 2L4 ij XijY2ijM ij xa X 2 ijYijMijyaþ L 2 ijXijMijxa ; pij ¼1 L2_ijY 2 ijF ij zaþ 1 2L4 ij Y3 ijM ij xa XijY2ijM ij yaþ L 2 ijYijMijxa ; qij¼ 1 L2 ij X2_ijFij_za 1 2L4 ij X2_ijYijMijxa X 3 ijM ij ya L 2 ijXijMijya ; rij_¼1 L2 ij XijMijza 1 L2 ij XijYijFijxa X 2 ijF ij ya ; sij¼ 1 L2_ijYijM ij za 1 L2_ij Y 2 ijF ij xa XijYijFijya ; tij₁¼ 1 L2 ij XijYijMijza; u ij 1¼ 1 L2 ij Y2_ijMij za; v ij 1¼ 1 L2 ij X2_ijMij za; wij₁¼ 1 2L2 ij XijYijMijxaþ Y 2 ijM ij ya ; xij1¼ 1 2L2ij XijYijMijxa X 2 ijM ij ya L 2 ijM ij ya ; yij₁¼ 1 2L2 ij X2_ijMij xaþ XijYijMijya ; zij₁¼ 1 2L2 ij Y2 ijM ij xa XijYijMijyaþ L 2 ijM ij xa ; tij₂¼ 1 L2 ij XijYijMijzaþ 1 L2 ij XijY2ijF ij xa X 2 ijYijFijya ; uij₂¼ 1 L2 ij Y2 ijM ij zaþ 1 L2 ij Y3 ijF ij xa XijY2ijF ij ya ; vij₂¼ 1 L2_ijX 2 ijM ij za 1 L2_ij X 2 ijYijFijxa X 3 ijF ij ya ; wij₂¼ 1 2L2 ij XijYijMijxaþ Y 2 ijM ij ya ; xij₂¼ XijFijzaþ 1 2L2 ij XijYijMijxa X 2 ijM ij ya L 2 ijM ij ya ; yij2¼ 1 2L2 ij X2_ijMij_xaþ XijYijMijya ; zij₂¼ YijFijzaþ 1 2L2 ij Y2_ijMij xa XijYijMijyaþ L 2 ijM ij xa :

In the above, the nodal forces for element ij, that is, Fij_xa; Fij_ya; . . . ; Mij_xa; Mij_ya; . . . ; with the left superscript ‘‘1’’ omitted, have been given in Eqs. (43) and (44), Xij=

Xi Xj, Yij= Yi Yj, and Lijis the length of element ij.

References

[1] J.H. Argyris, O. Hilpert, G.A. Malejannakis, D.W. Scharpf, On the geometrical stiﬀness matrix of a beam in space – a consistent v.w. approach, Comp. Methods Appl. Mech. Engrg. 20 (1979) 105–131. [2] Y.B. Yang, J.T. Chang, J.D. Yau, A simple nonlinear triangular plate

element and strategies of computation for nonlinear analysis, Comp. Methods Appl. Mech. Engrg. 178 (1999) 307–321.

[3] Y.B. Yang, L.J. Leu, Force recovery procedures in nonlinear analysis, Comput. Struct. 41 (6) (1991) 1255–1261.

[4] G. Turkalj, J. Brnic, J. Prpic-Orsic, ESA formulation for large displacement analysis of framed structures with elastic-plasticity, Comput. Struct. 82 (2004) 2001–2013.

[5] Y.B. Yang, H.T. Chiou, Rigid body motion test for nonlinear analysis with beam elements, J. Engrg. Mech., ASCE 113 (9) (1987) 1404–1419.

[6] Y.B. Yang, S.R. Kuo, Theory and Analysis of Nonlinear Framed Structures, Prentice Hall, Englewood Cliﬀs, NJ, 1994.

[7] Y.B. Yang, S.R. Kuo, Frame buckling analysis with full consider-ation of joint compatibilities, J. Struct. Engrg., ASCE 118 (1992) 871– 889.

[8] Y.B. Yang, J.D. Yau, L.J. Leu, Recent developments on geometri-cally nonlinear and postbuckling analysis of framed structures, Appl. Mech. Rev. 56 (4) (2003) 431–449.

[9] W. McGuire, R.H. Gallagher, R.D. Ziemian, Matrix Structural Analysis, second ed., John Wiley, NY, 2000.

[10] R.D. Cook, A plate hybrid element with rotational d.o.f. and adjustable stiﬀness, Int. J. Numer. Methods Engrg. 24 (1987) 1499– 1508.

[11] J.L. Batoz, K.J. Bathe, L.W. Ho, A study of three-node triangular plate bending elements, Int. J. Numer. Methods Engrg. 15 (1980) 1771–1812.

[12] H. Ziegler, Principles of Structural Stability, second ed., Birkha¨user, Stuttgart, Germany, 1977.

[13] Y.B. Yang, S.R. Kuo, Y.S. Wu, Incrementally small-deformation theory for nonlinear analysis of structural frames, Engrg. Struct. 24 (2002) 783–798.

[14] L.E. Malvern, Introduction to the Mechanics of a Continuous Medium, Prentice Hall, Englewood Cliﬀs, NJ, 1969.

[15] J. Argyris, L. Tenek, L. Olofsson, TRIC: a simple but sophisticated 3-node triangular element based on 6 rigid-body and 12 straining modes for fast computational simulations of arbitrary isotropic and lami-nated composite shells, Comp. Methods Appl. Mech. Engrg. 145 (1997) 11–85.

[16] Y.B. Yang, M.S. Shieh, Solution method for nonlinear problems with multiple critical points, AIAA J. 28 (12) (1990) 2110–2116. [17] G. Horrigmoe, P.G. Bergan, Numerical analysis of free-form shells by

ﬂat ﬁnite elements, Comp. Methods Appl. Mech. Engrg. 16 (1978) 11–35.

[18] K.S. Surana, Geometrically non-linear formulation for the three dimensional solid-shell transition ﬁnite elements, Comput. Struct. 15 (5) (1982) 549–566.