Geometrically nonlinear analysis of planar circular arches based on rigid element concept — A structural approach

Engineering Structures 30 (2008) 955–964 www.elsevier.com/locate/engstruct

Geometrically nonlinear analysis of planar circular arches based on rigid element concept — A structural approach J.D. Yau a , Y.B. Yang b,∗ a Department of Architecture and Building Technology, Tamkang University, Taipei, 10620, Taiwan b Department of Civil Engineering, National Taiwan University, Taipei, 10617, Taiwan

Received 27 February 2007; received in revised form 22 June 2007; accepted 26 June 2007 Available online 30 July 2007

Abstract To overcome the difficulty involved in selecting proper shape functions for simulating the bending-tension coupling of a curved beam, a nonconventional “structural” approach is presented in this paper. For curved-beam elements with small subtended angles, the elastic stiffness matrix is derived as the composition of two chordwise straight beam elements used to represent the curved beam. In contrast, the geometrical stiffness matrix of the curved beam is derived by the rigid element concept, through transformation of the geometrical stiffness matrix of the rigid straight beam spanning the two ends of the curved bean from the rectangular to the curvilinear coordinates. Compared with the conventional finite element procedures relying strongly on numerical integrations, the present approach has the advantage of being simple in formulation, but also explicit in expressions. The numerical studies indicate that the derived curved beam element has good convergence characteristics upon mesh refinement for the linear problems studied, and is capable of solving the stability and nonlinear problems involving large-displacement postbuckling response. c 2007 Elsevier Ltd. All rights reserved.

Keywords: Arch; Buckling; Curved beam; Geometric nonlinearity; Rigid body rule

1. Introduction Masonry arches were often used by ancient Romans as the supporting foundation of their bridges or buildings when the Roman civilization pervaded the European continent. Along with the expansion of their empire, semicircular arches became a typical symbol of the Roman architecture [1]. The basic feature of an arch is that it can sustain the self-weight and external loads in a compressive manner, i.e. without introducing any tensile force on the structure. Such a structural form is particularly suitable for building materials that are strong in resisting compressive, rather than tensile forces, including stone and concrete. With mechanical properties similar to the arch structures, a vertical curved beam can transfer the gravity loads through the bending-tension coupling effect, by which the material can be used in a more efficient manner than a straight beam. However, the coupling effect has made the deformation behaviors of curved beam structures much more complicated than structures composed of straight beams. ∗ Corresponding author. Tel.: +886 2 3366 4245; fax: +886 2 2362 2975.

E-mail address: [email protected] (Y.B. Yang). c 2007 Elsevier Ltd. All rights reserved. 0141-0296/$ - see front matter doi:10.1016/j.engstruct.2007.06.003

Concerning the finite element formulation of a curved beam element, the use of low-order independent polynomial interpolations for the normal (e.g. cubic polynomials) and tangential (e.g. linear polynomials) displacements may bring about significant errors in the computational results. Such a phenomenon, commonly known as membrane locking, has been reviewed by Prathap [2]. To overcome such a problem, efforts have been undertaken by researchers using approaches such as the strain element technique, reduced integration method, and isoparametric elements [3–5]. In the literature, a number of curved-beam elements were developed for the free vibration, buckling and postbuckling analyses of curved beams [6–20], based on various theories. It was pointed out by Pi et al. [19, 20] that some general purpose finite element codes, such as ABAQUS and ANSYS, were unable to accurately predict the elasto-plastic behaviour of members curved in space. Though some progress has been made in the past in the study of curved beam problems, the coupling behaviour (bendingtension and bending-torsion) of curved beams remains a mathematical hindrance in the derivation of a consistent displacement (strain) field for the curved beam element aimed

956

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964

at avoiding the membrane locking problem. In the out-of-plane buckling analysis of curved beams, Yang and co-workers [21, 22] have demonstrated that a curved beam can be simulated by a number of straight-beam elements through consideration of the equilibrium for structural joints connecting noncollinear members in the deformed configuration. As can be seen from the above review, there is an apparent lack of a simple and straightforward approach for formulating a curved beam element that is free of membrane locking in the elastic stiffness for linear analysis, and can duly take into account the effect of curvature on the geometrical stiffness for buckling analysis. On the other hand, planar curved beams constitute a special class of structures in engineering applications, for which both the in-plane deformation and outof-plane buckling behaviours have always been the concerns of structural designers. For the reasons stated above, a nonconventional structural approach will be proposed herein for deriving the planar curved beam element. For curved beams with small subtended angles, the elastic stiffness matrix of the curved beam will be derived as the composition of two chordwise straight-beam elements used to represent the curved beam. Based on the concept of rigid displacements [23,24], the geometric stiffness matrix of the curved beam will be derived by transforming the geometric stiffness matrix of the straight beam with identical nodal degrees from the rectangular coordinates to the curvilinear coordinates. To examine the applicability and accuracy of the curved-beam element presented herein, four examples on linear static, buckling and geometrically nonlinear analyses of curvedbeam structures will be studied.

Fig. 1. Finite element modelling of a circular curved beam.

Fig. 2. Modelling of a curved-beam element by two chordwise straight-beam elements.

3. Elastic stiffness matrix of planar curved-beam element

the curved beam into a number of curved-beam elements, as shown in Fig. 1. To circumvent the problem associated with the selection of proper shape functions for treating the coupling effect of bending-extension deformations, each planar curved-beam element with a small subtended angle (=2ϕ) will be replaced by two chordwise straight-beam elements, as schematically depicted in Fig. 1. The curved beam element of concern is shown in Fig. 2 with a radius R and subtended angle 2ϕ. The z-axis in Fig. 2 represents one of the principal axes of the cross-section, and the x-axis is tangent to the curvilinear axis of the beam. Concerning derivation of the elastic stiffness matrix for describing the linear behaviour, the curved-beam element acb in Fig. 2 will be approximated by the two chordwise straight beam elements ac and cb, also named as beam elements i and j, respectively, which share a common auxiliary node c at the midpoint of the curved- beam element. In addition, θi and θ j represent the included angles of segments ac and cb, respectively, with respect to the chord ab of the curved element in the X –Z coordinate system shown in Fig. 2. With reference to Fig. 2, the elastic stiffness matrix of the straight-beam element along the chord ac or cb can be expressed in a partitioned form in the respective local coordinates as # " [k11,n ] [k12,n ]   (1) ke,n n=i, j = [k21,n ] [k22,n ]

For analysis of a circular curved beam subjected to external loadings by the finite element method, it is reasonable to divide

in which the subscript n (=i, j) indicates the straight beam element of concern. The submatrices in Eq. (1) are [22,25]

2. Basic assumptions for planar curved beams In this section, a simple, nonconventional structural approach will be presented for formulation of the planar curved-beam element, which can be regarded primarily as extensions of the theories for the straight-beam element concerning derivation of the elastic stiffness matrix and geometrical stiffness matrix. The following are the assumptions adopted for the planar curved beam: (1) The material is elastic and homogeneous; (2) The cross-section of the curved beam is uniform and doubly symmetrical; (3) Every cross-section remains rigid, i.e. undistorted, during deformation; (4) The length and radius of the curved beam are large in comparison with the cross-sectional dimensions of the beam; (5) The shearing deformation on the curved beam is negligible; and (6) Only concentrated loadings are allowed to act at the two ends of the curved-beam element.

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964

EA  Ln    0  k11,n =     0 

0

0



 6E I  12E I  − 2 , L 3n Ln   6E I 4E I  − 2 Ln Ln  EA − 0  Ln  12E I T  0    − 3 k12,n = k21,n =   Ln   6E I 0 L 2n   EA 0 0  Ln   12E I 6E I     0   k22,n =   L 3n L 2n  ,    4E I  6E I 0 Ln L 2n

(2a)

0



 6E I   − 2 , Ln   2E I  Ln

Fig. 3. Nodal forces of a curved beam element.

(2b)

(2c)

where E = elastic modulus, A = cross sectional area, I = moment of inertia, and L n |n=i, j = length of the element i or j, as shown in Fig. 2. To assemble the stiffness matrices of the two straight-beam elements connected at the auxiliary node c, the stiffness matrices in each local coordinate system should be transformed into the common X –Z coordinate system [25], that is,  i  i   [kaa ] [kac ] T , [ki ] = [T (θi )] ke,i [T (θi )] = i T i [kac ] [kcc ] # " (3a,b) j j    T    [kcc ] [kcb ] k j = T (θ j ) ke, j T (θ j ) , j j [kcb ]T [kbb ] where the submatrices in Eqs. (3a,b) are given by i ] = [Ti ]T [k11,i ][Ti ], [kaa i [kac ] = [Ti ]T [k12,i ][Ti ], i [kcc ] = [Ti ]T [k22,i ][Ti ],

(4a–c)

j

[kcc ] = [T j ]T [k11, j ][T j ], j [kcb ] = [T j ]T [k12, j ][T j ], j [kbb ] = [T j ]T [k22, j ][T j ],

957

(4d–f)

and the transformation matrices [Ti ] and [T j ] can be expressed as   [Tn ] [0] [T (θn )]n=i, j = [0] [Tn ]    cos θn sin θn 0 − sin θn cos θn 0  [0]     0 0 1     =  . (5) cos θ sin θ 0 n n    − sin θn cos θn 0  [0] 0 0 1 For the curved beam element acb shown in Fig. 2, we shall consider the special case where the connecting node c is located at the midpoint of the curved element, i.e. the

two chordwise straight elements have identical length. In this study, a clockwise rotation θ y is taken as positive. Based on the hypothesis of shallowness, the subtended angle 2ϕ of the curved beam element is assumed to be small. As a consequence, the tangential angle ϕ at the two ends a and b with respect to the chords ac and cb can be approximated by the included angles θi and −θ j , that is, sin θi = sin(ϕ − θi ), or, θi = −θ j = ϕ/2. By introducing this relation into Eqs. (4) and (5) and assembling the two element stiffness matrices by the direct stiffness method [25], one can obtain the following approximate elastic stiffness matrix for the curved-beam element acb:   i i [kaa ] [kac ] [0] n= Xj  i T j j  i (6) [ks ]9x9 = [kn ] = [kac ] [kcc ] + [kcc ] [kcb ] . j j n=i [0] [kcb ]T [kbb ] Then, the matrix equation of force equilibrium for the curvedbeam element with nodes a, b, and c, as shown in Fig. 3, can be given as follows:     i  i [kaa ] [kac ] [0] {da } {da }  i T j j  i [ks ] {dc } = [kac ] [kcc ] + [kcc ] [kcb ] {dc }     j j {db } {db } [0] [kcb ]T [kbb ]   { pa } = {0} , (7)   { pb } where {dn }n=a,b,c = nodal displacement vectors, and { pa } and { pb } are the nodal force vectors at nodes a and b, respectively, namely,      un   Fxn  {dn } = vn { pn } = Fzn , . (8a, b)     θ yn n=a,b,c M yn n=a,b Considering only the force equilibrium of the auxiliary node c, one can write the following equation of equilibrium based on Eq. (7):   j j i T i [kac ] {da } + [kcc ] + [kcc ] {dc } + [kcb ]{db } = {0}, (9) from which the nodal displacement vector {dc } can be solved and expressed in terms of the nodal displacement vectors {da } and {db } as  −1   j j i i T {dc } = − [kcc ] + [kcc ] [kac ] {da } + [kcb ]{db } . (10) Then, substituting Eq. (10) into Eq. (7) yields the following condensed stiffness equation for the curved-beam element with

958

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964

 

i i [kaa ] − [kac ]

 k¯e,c = 

 −1 j i i T [kcc ] + [kcc ] [kac ] symm.

  −1 j j i [kcc ] + [kcc ] [kcb ]   −1 . j j T j j i [kbb ] − [kcb ] [kcc ] + [kcc ] [kcb ] i −[kac ]

Box I.

the two nodes a and b:       {da } { pa } k¯e,c = {db } { pb }

(11)

where [k¯e,c ] is the elastic stiffness matrix of the curved-beam element in the X –Z coordinate system, given in Box I. To obtain the elastic stiffness matrix of the curved-beam element in the curvilinear coordinates x–z as indicated in Fig. 2, a further transformation of the coordinate system is necessary, that is,     ke,c = [Tc ]T k¯e,c [Tc ] , (12) where the matrix [Tc ] represents the transformation between the X –Z system and the local x–z curvilinear system for the two nodes a and b of the element with rotations ϕ and −ϕ, respectively, that is, cos ϕ  sin ϕ   0 [Tc (ϕ)] =   0  0 0

− sin ϕ cos ϕ 0 0 0 0

0 0 1 0 0 0

0 0 0 cos ϕ − sin ϕ 0

0 0 0 sin ϕ cos ϕ 0

0 0  0 . 0 0 1

(13) Fig. 4. A two-dimensional beam undergoing rigid body rotation: (a) C1 configuration, (b) C2 configuration.

The above procedure for derivation of the elastic stiffness matrix for the planar curved beam element with two nodes a and b based on the straight-beam element stiffness matrices is characterized by the fact that it is based purely on the consideration of structural geometry. Of course, the accuracy of the curved-beam element can be improved by using more chordwise straight elements, rather than only two elements, in the modelling. 4. Geometric stiffness matrix of planar curved beam element In an incremental nonlinear analysis based on the updated Lagrangian formulation, the incremental element stiffness equation for a planar beam element from the last calculated configuration C1 (see Fig. 4(a)) to the current deformed configuration C2 (see Fig. 4(b)) can be written as [22] ([ke ] + [k g ]){u} + {1 F} = {2 F},

(14)

where [ke ] = elastic stiffness matrix, [k g ] = geometrical stiffness matrix, {1 F} = initial force vector acting at the C1 configuration, and {2 F} = current force vector at the deformed configuration C2 . The nodal displacement vector {u} for the incremental step from C1 to C2 in terms of the x 1 –y 1 coordinates (see Fig. 4) is

 {u} = u a va θa u b vb θb . (15)

Correspondingly, the initial force vector {1 F} acting on the element is (see Fig. 4(a))

T {1 F} = Fxa Fza M ya Fxb Fzb M yb . (16) Or, by equilibrium,

{1 F} = −Fxb −Fzb

M ya

Fxb

Fzb

M yb

T

(17)

where Fxb = axial force, Fzb = shear force, M yb = bending moment, and Fzb = (M ya + M yb )/L. Before we proceed with derivation of the geometrical stiffness matrix for the curved beam element, two important concepts must be mentioned here. First, as far as the rigid behaviour is concerned, the behaviour of a two-node curved beam depends exclusively on the behaviour of the two end points a and b, but not on the shape or elastic properties of the element, and thus is identical to the behaviour of the straight beam consisting of the same end points a and b. Second, the geometrical stiffness matrix for a two-dimensional straight beam can be derived in an approximate, but rather accurate manner, by considering only the rigid body rotations [23,24,26], due to the fact that that the rigid body displacements constitute a great portion of the buckling displacement of the beam elements, as schematically shown in Fig. 5. In this section, the geometric stiffness matrix of a two-dimensional straight beam element based on the rigid concept will be first summarized.

959

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964

{2 F} ∼ = −Fxb + Fzb θr

−Fxb θr − Fzb

M ya

Fxb − Fzb θr

Fzb + Fxb θr

M yb

T

.

Box II.

Substitution of Eqs. (17) and Box II into Eq. (18) yields the following expression for the geometrical stiffness matrix [k g ]:   k g {u r } = {2 F} − {1 F}

T = Fzb θr −Fxb θr 0 −Fzb θr Fxb θr 0 . (22) For a beam element under a small rigid rotation θr , the following approximation can be adopted: va − vb θr ∼ . (23) = L By the preceding relation, one can manipulate the right-hand side of Eq. (22) to obtain the geometrical stiffness matrix for the straight-beam element as follows [24,28]: 0    1  kg =  L 

−Fzb Fxb

0 0 0

symm.

Fig. 5. Schematic representation of buckling shape for a cantilever.

As shown in Fig. 4(b), when a straight-beam element is subjected to a small rigid rotation θr , no forces will be generated by the elastic stiffness matrix [ke ]. Thus, the incremental element stiffness equation in Eq. (14) reduces to [k g ] {u r } + {1 F} = {2 F}

(18)

where {u r } denotes the rigid displacement vector,  {u r }T = 0 0 θr 0 −Lθr θr ,

(19)

and {2 F} denotes the resulting forces acting on the element at C2 after the rigid rotation, as indicated in Fig. 4(b). According to the rigid body rule [27], the initial nodal forces acting on an element should rotate following the rigid rotation, while their magnitudes remain unchanged. Consequently, the nodal forces {2 F} acting on the element after the rigid rotation should be equal in magnitude to the initial nodal forces {1 F} acting at C1 , but rotate by the angle θr , as shown in Fig. 4(b), i.e. {2 F} = [T (θr )]{1 F},

(20)

where for a small rigid rotation θr , the transformation matrix [T (θr )] can be approximated as follows: 1 θr  0 [T (θr )] ∼ = 0 0 0

−θr 1 0 0 0 0

0 0 0 0 1 0 0 1 0 θr 0 0

0 0 0 −θr 1 0

0 0  0 . 0 0 1

(21)

As a result, the nodal forces {2 F} reduce to the equation in Box II.

0 Fzb 0 0

Fzb −Fxb 0 −Fzb Fxb

0 0  0 . 0 0 0

(24)

As can be seen from Eq. (24), the geometrical stiffness matrix for the straight beam depends only on the nodal forces and element length, but not on the material and sectional properties of the beam, since only rigid displacements have been considered. Such a matrix is exactly the stiffness matrix for the rigid beam element [23,24]. For two-point rigid solid elements, whether they are straight or curved, or whether they are thick or thin, they will exhibit the same rigid body behaviour as long as they have identical degrees of freedom at the two end points. Thus, the geometrical stiffness matrix as derived above for the rigid straight beam element is exactly the same as that for the rigid curved beam element, except that different coordinate systems have been used, namely, a rectangular coordinate system is used for the straight beam, and a curvilinear coordinate system for the curved beam. Thus, a transformation of the coordinate systems is required to obtain the geometrical stiffness matrix for the curved beam from the straight beam. With reference to Fig. 6, the nodal forces (Fxa , Fza , Fxb , Fzb , M ya , M yb ) of the straight beam can be related to the nodal forces ( f xa , f za , f xb , f zb , M ya , M yb ) of the curved beam as      cos ϕ − sin ϕ 0  f xa   Fxa  Fza =  sin ϕ cos ϕ 0 f za ,     M ya 0 0 1 M ya      (25a, b) cos ϕ sin ϕ 0  f xb   Fxb  Fzb = − sin ϕ cos ϕ 0 f zb .     M yb 0 0 1 M yb From the first two equations of Eq. (25b), the following equations of equilibrium can be written: Fxb = f xb cos ϕ + f zb sin ϕ.

(26a)

960

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964

0     1  k¯ g =  L 

f xb sin ϕ − f zb cos ϕ f xb cos ϕ + f zb sin ϕ

0 0 0

0 − f xb sin ϕ + f zb cos ϕ 0 0

symm.

− f xb sin ϕ + f zb cos ϕ −( f xb cos ϕ + f zb sin ϕ) 0 f xb sin ϕ − f zb cos ϕ f xb cos ϕ + f zb sin ϕ

0 0  0 . 0 0 0

Box III.

0           1  k g,c = [Tc ]T k¯ g [Tc ] ∼ = k g,c0 + 1k g,c0 =  2Rϕ   3ϕ f − 2 f xb zb  1   +  2R  

f xb + 2 f zb ϕ 2( f zb − f xb ϕ)

0 0 0

− f zb f xb

0 f zb 0 0

symm.

f xb ϕ f xb + 2 f zb ϕ 0 − f xb + 2 f zb

symm.

0 0 0

0 0  0  0 0 0 0

f zb − f xb 0 − f zb f xb

− f xb 0 0  0 0  − f xb 0 2( f xb ϕ − f zb ) 0 0

Box IV.

1 ϕ  0 [Tc (ϕ)] '  0 0 0

Fig. 6. A two-node planar curved beam element.

Fzb =

M ya + M yb = − f xb sin ϕ + f zb cos ϕ, L

(26b)

Substituting Eqs. (26) into Eq. (24) yields the geometrical stiffness matrix of the straight beam in terms of the nodal forces of the curved beam as given in Box III. The preceding geometrical stiffness matrix can be transformed to the one for the curved beam in the curvilinear coordinate system as follows: 

   k g,c = [Tc ]T k¯ g [Tc ] ,

(27)

where the transformation matrix [Tc ] has been given in Eq. (13). In the finite element modelling, we shall use only curved beam elements with small subtended angles. Under such a condition, the equations of equilibrium in Eq. (26) and the transformation matrix in Eq. (13) reduce to Fxb = f xb + f zb ϕ, Fzb = − f xb ϕ + f zb ,

(28a, b)

−ϕ 1 0 0 0 0

0 0 1 0 0 0

0 0 0 1 −ϕ 0

0 0 0 ϕ 1 0

0 0  0 . 0 0 1

(28c)

By the use of the preceding relations, and by neglecting the terms involving ϕ 2 and higher orders, the geometrical stiffness matrix for the curved beam in Eq. (27) can be rederived and expressed in an explicit form as given in Box IV. Here, one can observe that the geometrical stiffness matrix for a rigid curved beam with a small subtended angle consists of two parts. The first part [k g,c0 ] represents the uncoupled effect of the nodal forces f xb and f zb , and the second part [1k g,c0 ] the coupling effect between the nodal forces f xb and f zb due to the curvature effect of the curved beam. Moreover, for the limit case when (R → ∞, ϕ → 0) but 2Rϕ → L, then f xb ∼ = Fxb , f zb ∼ = Fzb and the geometrical stiffness matrix for the curved beam given in Box IV reduces to that in Eq. (24) for the straight beam. From the equation given in Box IV, one observes that as the angle ϕ approaches zero (ϕ → 0), the accuracy of the matrix [k g,c ] will be affected, since ϕ appears in the denominator. This forms a limit for reducing the angle ϕ or for subdivision of the beam. For the sake of computational accuracy, it is suggested that the denominator 2Rϕ be kept in a order close to the length L of the curved beam, i.e. 2Rϕ ≈ L or ϕ ≈ L/(2R). At this point, some remarks will be made concerning the properties of the stiffness matrices derived. First of all, the assumption of shallowness for the curved beam, i.e. with small

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964

961

of convergence or the number of iterations. In this study, the curved beam element matrix derived herein, including the elastic stiffness matrix [ke,c ] in Eq. (12) and the geometric stiffness matrices [k g,c ] in Eq. (27), will be employed to construct the structural stiffness matrix [K ]. The corrector phase deals with recovery of the element forces at C2 from the existing nodal forces and the element displacement increments {1u} made available from the structural displacement increments {1U } via the predictor. This phase determines primarily the accuracy of the solution. In this study, the element forces {22 f } at C2 with reference to C2 will be computed as follows: {22 f } = {11 f } + {1 f } , Fig. 7. Motion of body in three-dimensional space.

subtended angles 2ϕ, has been adopted in derivation of the elastic stiffness matrix [ke,c ] in Section 3. Such an assumption was not definitely needed in derivation of the geometrical stiffness matrix [k g,c ] in Section 4. For instance, the full geometrical stiffness matrix [k g,c ] as given in Eq. (27) (for all ranges of ϕ) can be directly used. However, it has been numerically demonstrated that basically no difference can be made between the results obtained using the full matrix in Eq. (27) or the simplified matrix given in Box IV (with small angles ϕ). As such, we prefer to use the simplified version given in Box IV for its explicitness. 5. Incremental nonlinear analysis Three typical configurations can be identified for the incremental nonlinear analysis of structures, as shown in Fig. 7. The configuration C0 denotes the undeformed configuration, C1 the last calculated configuration, and C2 the current configuration of the structure. In an incremental–iterative analysis based on the updated Lagrangian formulation, the element stiffness equations can be assembled to form the incremental equations of equilibrium of the structure with respect to the last known or initial configuration C1 , given the applied load increments {2 P} − {1 P}: [K ] {1U } = {2 P} − {1 P}

(29)

where [K ] denotes the structural stiffness matrix assembled from the element stiffness matrices, {1U } the displacement increments of the structure generated during the incremental step from C1 to C2 , and {1 P} and {2 P} respectively denote the external loads applied on the structure at C1 and C2 . The solution of the incremental nonlinear equations in Eq. (29) will be attempted by a combination of the incremental and iterative procedures. Three key phases are essential to an incremental–iterative nonlinear analysis, i.e. the predictor, corrector and equilibrium checking phases [22,29]. The predictor phase is concerned with solution of the structural displacement increments {1U } under given load increments {1P} = {2 P} − {1 P} from the structural equations in Eq. (29), which may affect the speed

(30)

where {11 f } denotes the initial nodal forces acting on the curved beam element at C1 with reference to C1 , and {1 f } the force increments generated during the incremental step from C1 to C2 . It should be noted that the nodal forces {11 f } that were initially acting at C1 with reference to C1 have been treated as the forces acting at C2 with reference to C2 , according to the rigid body rule [22,27]. The elastic force increments {1 f } of the element will be computed as    
{1 f } = ke,c + k g,c {1u} , (31) where {1u} denotes the element displacement increments. For each iterative step, the displacement increments {1u} can be decomposed into two parts as the natural deformations {1u}n and rigid displacements {1u}r . Since the present curved beam element is rigid body qualified, no actions will be induced by the stiffness matrices [ke,c ] + [k g,c ] when undergoing the rigid displacements {1u}r . Thus, the force increments {1 f } can also be computed by considering merely the elastic stiffness matrix, that is, {1 f } = [ke,c ]{1u}n . After all the element forces are computed and expressed in the deformed configuration of the structure, one can compare them with the total applied loads {2 P} to obtain the unbalanced forces of the structure at C2 , as is typical for the equilibrium-checking phase. Whenever the unbalanced forces are greater than preset tolerances, the unbalanced forces should be regarded as the applied loads and another iteration involving the three phases should be repeated. Details concerning the incremental–iterative procedure for the geometrically nonlinear analysis of structures are available in Ref. [22]. For the purpose of tracing the load–deflection response of the curved beam structure, the generalized displacement control (GDC) method proposed by Yang and Shieh [30] will be adopted for its general stability in dealing with multi loops in the postbuckling response. By this method, the load increment size is determined as a function of the general stiffness parameter (GSP) [30], which serves as a reliable indicator for reversing the direction of loading once a limit point is detected (see Fig. 8). A further description of this method is available in Ref. [30]. 6. Numerical verifications In order to examine the validity and accuracy of the curved beam element proposed herein, linear, buckling and

962

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964

Table 1 Convergence study of the curved beam element on deflection analysis No. of elements

4

10

20

Analytical (1)

β = 90◦

0.075 m (2.7%a )

0.078 m (1.2%) 0.468 m (0.6%)

0.079 m (0.0%) 0.470 m (0.2%)

0.079 m 0.471 m

β = 180◦

0.451 m (4.2%)

a Percentage error compared with analytical solution.

Fig. 8. General characteristics of nonlinear system.

geometrically nonlinear analyses for different curved beam structures will be carried out and the results obtained will be compared with those existing in the literature. 6.1. Deflection of cantilevered curved beams with tip load In this example, the linear deflection analysis of a cantilevered curved beam subjected to a vertical load at the free end will be used to illustrate the applicability of the elastic stiffness matrix derived for the curved beam element. For the cantilevered circular ring with a subtended angle β subjected to a vertical tip load P in Fig. 9, the following properties are adopted: P = 1 kN, R = 10 m, E A = 106 kN, and E I = 104 kN m2 . For the two subtended angles π/2 and π considered, the analytical solutions for the vertical deflection at the tip are [31]  π P R3 π   , for β = , 2 1 = 4 EI 3 (32)   3π P R , for β = π. 2 EI With the curved beam represented by different numbers of curved beam elements, the numerical results obtained by the present approach for the vertical tip deflection are compared with the analytical ones in Table 1. As can be seen, larger errors exist for curved beams with larger subtended angles. However, the error decreases as the finite element mesh is refined. Overall, the present solutions converge in an asymptotic manner to the analytical solutions upon mesh refinement. Thus, it is confirmed the elastic stiffness matrix derived in this paper for the curved beam element is suitable for the linear analysis of planar circular arches. 6.2. Buckling analysis To demonstrate the applicability of the geometrical stiffness matrix derived for the curved beam element given in Box IV

Fig. 9. Cantilevered circular beam under vertical load.

to general buckling analysis, the in-plane buckling analysis of a pinned circular arch under radial compression q distributed along the axis of the beam will be investigated, as shown in Fig. 10. The distributed radial load q is simulated as statically equivalent concentrated loads acting on the two nodes of each element used to represent the arch. The following properties are adopted for the curved beam: R = 24.257 m, E = 200 GPa, A = 20.26 × 10−4 m2 , and I = 32.67 × 10−8 m4 . For a uniformly compressed circular arch, the radial load q will produce a constant axial force F0 on the beam, which is equal to q R, and the critical load can be obtained from Ref. [32] as   E I π2 F0,cr = qcr R = 2 − 1 . (33) R ϕ2 For the arch represented by different numbers of curved-beam elements, the critical loads computed for different subtended angles 2ϕ have been listed in Table 2. Obviously, larger errors occur for arches with larger subtended angles using the same number of elements, due to violation of the shallowness hypothesis for larger subtended angles. However, the error decreases drastically as more elements are used for the mesh. In fact, the level of errors is an indicator of the shallowness applicable. It is confirmed that the present curved beam element can be reliably used in modelling the buckling loads of pinned circular arches with various subtended angles, if a sufficient number of elements is used in the representation. 6.3. Postbuckling analysis of shallow arch The geometrical nonlinear behaviour of the shallow arch shown in Fig. 11 was analysed by Yang and Kuo [22], Chapter 7, by the GDC method using a 25 straight-beam element representation. The following data are adopted in this study: L = 100 in., h = 5 in., E = 2000 psi, I = 1 in.4 , and A = 1 in.2 . In this finite element modelling, the arch was first discretized into 14 curved beam elements. Two loading cases

963

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964 Table 2 Convergence study of the present curved beam element on buckling loads No. of elements

4

2ϕ = 30◦ 2ϕ = 60◦ 2ϕ = 90◦

18.68 kN (18%a )

2ϕ = 150◦

4.66 kN (19.8%) 2.06 kN (23.3%) 0.71 kN (34%)

10

20

Analytical (F0,cr )

16.22 kN (2.1%) 4.10 kN (5.4%) 1.81 kN (8.5%) 0.61 kN (15%)

16.17 kN (1.8%) 4.02 kN (3.3%) 1.77 kN (6.3%) 0.59 kN (11%)

15.89 kN 3.89 kN 1.67 kN 0.53 kN

a Percentage error compared with analytical solution.

Fig. 13. Hinged circular arch. Fig. 10. Pinned-end circular arch under radial loads.

Fig. 11. Shallow arch.

Fig. 14. Load–deflection curves for arch under symmetrical loading.

6.4. Postbuckling analysis of hinged semicircular arch

Fig. 12. Load–deflection curves for shallow arch.

are considered. For the perfect loading case, the vertical load is applied at the apex of the arch. For the imperfect loading case, the element nearest to the apex on one side was further divided into two elements to provide an additional node for application of loading to produce the effect of imperfection. From Fig. 12, it is clear that the results obtained by the present approach agree very well with those by Yang and Kuo [22]. This example confirms the validation and accuracy of the present curved beam element for use in geometrical nonlinear analysis.

Fig. 13 shows a hinged semicircular arch subjected to a central vertical load. The following data are considered: L = 100 in., E = 2000 psi, I = 1 in.4 , and A = 1 in.2 . Two loading cases are studied. One is called the symmetrical loading in which the vertical load is applied at the apex of the arch. The second case is the asymmetric loading in which the imperfection is introduced by shifting the vertical load to the node nearest to the centre of the arch, as indicated in Fig. 13. By using 26 curved-beam elements in the representation, the load deflection curves computed for the hinged circular arch subjected to the symmetrical and asymmetrical loadings have been plotted in Figs. 14 and 15, respectively. As can be seen, the results obtained by the present curved beam approach agree generally well with those by Yang and Kuo [22], even though the circular arch has experienced a very large deformation, such that the vertical deflection almost reaches the span dimension of the arch.

964

J.D. Yau, Y.B. Yang / Engineering Structures 30 (2008) 955–964

Fig. 15. Load–deflection curves for arch under asymmetrical loading.

7. Concluding remarks To circumvent the mathematical hindrance encountered in formulation of the curved beam element using the conventional shape function approach, a nonconventional structural approach is proposed in this paper. The elastic stiffness and geometrical stiffness matrices of the curved beam element are treated in different ways. For the curved beam with small subtended angles, the elastic stiffness matrix is derived as the composition of two chordwise straight beam elements used to represent the curved beam. In contrast, the derivation of the geometrical stiffness matrix is based on the concept of rigid element, that is, by transformation of the geometrical stiffness matrix derived for the rigid straight beam sharing the same end nodes as the curved beam from the rectangular to the curvilinear coordinates. The advantage of the present approach is that there is no need to deal with the bending-tension coupling effect of the curved beam, thereby avoiding problems such as membrane locking, and that all the stiffness matrices are presented in explicit form. The applicability and accuracy of the derived matrices for the curved beam element have been demonstrated in the solution of several examples involving the linear, stability and postbuckling nonlinear behaviours of curved beam structures. Compared with the use of curved beam element with no restriction on the subtended angles, e.g. Chapter 7 of Ref. [22], slightly more elements should be used by the present approach. On the other hand, the present approach outperforms the conventional straight beam approach, as the curvature and instability properties have been considered in the formulation, though in an approximate manner. Acknowledgements The authors would like to thank the financial support through a grant (NSC 95-2211-E-032-053) from National Science Council of Taiwan. References [1] Heyman J. The science of structural engineering. Imperial College Press; 1999.

[2] Prathap G. The curved beam/deep arch/finite ring element revisited. Int J Numer Methods Eng 1985;21:389–407. [3] Choi JK, Lim JK. Simple curved shear beam elements. Comm Numer Methods Engrg 1993;9:659–69. [4] Ashwell DG, Sabir AB. Limitations of certain curved finite elements when applied to arches. Int J Mech Sci 1971;13:133–9. [5] Ashwell DG, Sabir AB, Roberts TM. Further studies in the application of curved finite elements to circular arches. Int J Mech Sci 1971;13:507–17. [6] Sabir AB, Ashwell DG. A comparison of curved beam elements when used in vibration problems. J Sound Vib 1971;18:555–63. [7] Veletsos AS, Austin WJ, Lopes Pereira CA, Wung SJ. Free inplane vibration of circular arches. ASCE J Eng Mech Div 1972;98:311–29. [8] Sabir AB, Lock AC. Large deflexion geometrically nonlinear finite element analysis of circular arches. Int J Mech Sci 1973;15:37–47. [9] Dawe DJ. Curved finite elements for the analysis of shallow and deep arches. Comput Struct 1974;4:559–80. [10] Dawe DJ. Numerical studies using circular arch finite elements. Comput Struct 1974;4:729–40. [11] Meck HR. An accurate polynomial displacement function for finite ring elements. Comput Struct 1980;11:265–9. [12] Prathap G, Babu CR. Field consistency and violent oscillations in the finite element method. Int J Numer Methods Eng 1987;24:2017–33. [13] Balasubramanian TS, Prathap G. A field consistent higher order curved beam element for static and dynamic analysis of stepped arches. Comput Struct 1989;33:281–8. [14] Pandian N, Appa Rao TVSR, Chandra S. Studies on performance of curved beam finite elements for analysis of thin arches. Comput Struct 1989;31:997–1002. [15] Yang YB, Kuo SR, Cherng YD. Curved beam elements for nonlinear analysis. ASCE J Eng Mech 1989;115(4):840–55. [16] Krishnan A, Dharmaraj S, Suresh YJ. Free vibration studies on arches. J Sound Vib 1995;186:856–63. [17] Krishnan A, Suresh YJ. A simple cubic linear element for static and free vibration analyses of curved beams. Comput Struct 1998;68:473–89. [18] Pi YL, Trahair NS. Non-linear buckling and postbuckling of elastic arches. Eng Struct 1998;20(7):571–9. [19] Pi YL, Bradford M, Uy B. A spatially curved-beam element with warping and Wagner effects. Int J Numer Methods Eng 2005;63:1342–69. [20] Pi YL, Bradford M, Uy B. A rational elasto-plastic spatially curved thinwalled beam element. Int J Numer Methods Eng 2007;70:253–90. [21] Yang YB, Kuo SR, Yau JD. Use of straight-beam approach to study buckling of curved beams. ASCE J Struct Eng 1991;117(7):1963–78. [22] Yang YB, Kuo SR. Theory and analysis of nonlinear framed structures. Singapore: Prentice Hall; 1994. [23] Yang YB, Lin SP, Chen CS. Rigid body concept for geometric nonlinear analysis of 3d frames, plates and shells based on the updated Lagrangian formulation. Comp Methods Appl Mech Eng 2007;196:1178–92. [24] Yang YB, Lin SP, Leu LJ. Solution strategy and rigid element for nonlinear analysis of elastically structures based on updated Lagrangian formulation. Eng Struct 2007;29(6):1189–200. [25] McGuire W, Gallagher RH. Matrix structural analysis. New York: John Wiley & Sons, Inc.; 1979. [26] Yang YB, Chang JT, Yau JD. A simple nonlinear triangular plate element and strategies of computation for nonlinear analysis. Comp Meth Appl Mech Eng 1999;178:307–21. [27] Yang YB, Chiou HT. Rigid body motion test for nonlinear analysis with beam elements. ASCE J Eng Mech 1987;9:1404–19. [28] Yang YB, Kuo SR, Wu YS. Incrementally small-deformation theory for nonlinear analysis of structural frames. Eng Struct 2002;24:783–93. [29] Yang YB, Yau JD, Leu LJ. Recent development on nonlinear and postbuckling analysis of framed structures. Appl Mech Rev 2003;56(4): 431–49. [30] Yang YB, Shieh MS. Solution method for nonlinear problems with multiple critical points. AIAA J 1990;28(12):2110–6. [31] Timoshenko SP. Strength of material: Part I. 3rd ed. NJ: D. Van Nostrand Company, Inc.; 1955. [32] Timoshenko SP, Gere JM. Theory of elastic stability. 2nd ed. NY: McGraw-Hill; 1961.