Co-rotational formulation for geometric nonlinear analysis of doubly symmetric thin-walled beams

Co-rotational formulation for geometric nonlinear analysis of doubly symmetric thin-walled beams

Wen Yi Lin

, Kuo Mo Hsiao

aDepartment of Mechanical Engineering, Der Lin Institute of Technology, 1 Alley 380, Ching Yun Road, Tucheng, Taiwan, ROC

bDepartment of Mechanical Engineering, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu, Taiwan, ROC Received 24 December 1999; received in revised form 14 September 2000

Abstract

A doubly symmetric thin-walled beam element with open section is derived using co-rotational (CR) total Lagrangian (TL) for- mulation. The e€ects of deformation-dependent third-order terms of element nodal forces on the buckling load and post-buckling behavior are investigated. All coupling among bending, twisting, and stretching deformations for beam element is considered by consistent second-order linearization of the fully geometrically nonlinear beam theory. However, all third-order terms of nodal forces, which are relevant to the twist rate, rate of twist rate and curvature of the beam axis are also considered. An incremental-iterative method based on the Newton±Raphson method combined with constant arc length of incremental displacement vector is employed for the solution of nonlinear equilibrium equations. The zero value of the tangent sti€ness matrix determinant of the structure is used as the criterion of the buckling state. A parabolic interpolation method of the arc length is used to ®nd the buckling load. Numerical examples are presented to demonstrate the accuracy and eciency of the proposed element and to investigate the e€ect of third-order terms of element nodal forces on the buckling load and post-buckling behavior of doubly symmetric thin-walled beams. Ó 2001 Elsevier Science B.V. All rights reserved.

Keywords: Co-rotational formulation; Thin-walled beam; Geometric nonlinearity; Buckling; Postbuckling

1. Introduction

Due to reduction of weight, material and cost, thin-walled beams with open section are extensively used in aerospace and aircraft structures, and are often designed to work under post-buckling conditions. Such

¯exible structures can undergo large displacements and rotations without exceeding their elastic limits. To understand the behaviors of such ¯exible structures and to evaluate their elastic limits many di€erent formulations and numerical procedures for the buckling and post-buckling analysis of thin-walled beams have been proposed [1±35]. The buckling of the beam structures is caused by the coupling among bending, twisting, and stretching deformations of the beam members. Thus the buckling analysis is a subtopic of nonlinear rather than linear mechanics [7]. Currently, the most popular approach for the analysis of three- dimensional beam is to develop ®nite element models. The formulations, which have been used in the literature, might be divided into three categories: total Langrangian (TL) formulation, updated Lagrangian (UL) formulation, and co-rotational (CR) formulation. In order to capture correctly all coupling among bending, twisting, and stretching deformations of the beam elements, the formulation of beam elements might be derived by the fully geometrically nonlinear beam theory [36]. The exact expressions for the el- ement nodal forces, which are required in a TL formulation for large displacement/small strain problems,

www.elsevier.com/locate/cma

*Corresponding author. Fax: +886-35-720-634.

E-mail address: [email protected] (K.M. Hsiao).

0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved.

PII: S 0 0 4 5 - 7 8 25 ( 0 1 ) 0 0 21 2- 2

are highly nonlinear functions of element nodal parameters. However, the dominant factors in the geo- metrical nonlinearities of beam structures are attributable to ®nite rotations, the strains remaining small.

For a beam structure discretized by ®nite elements, this implies that the motion of the individual elements to a large extent will consist of rigid body motion. If the rigid body motion part is eliminated from the total displacements and the element size is properly chosen, the deformational part of the motion is always small relative to the local element axes. Thus in conjunction with the CR formulation, the higher-order terms of nodal parameters in the element nodal forces may be neglected by consistent second-order linearization [32,36]. However, the values of twist rate, the rate of twist rate and curvature of the beam axis are de- formation dependent, not element size dependent. Thus their values may not always be much smaller than unity. It seems that some third-order terms of the element nodal forces, which are relevant to the twist rate, the rate of twist rate and curvature of the beam axis, may not be negligible for some cross-sections with large rotations. In [34], the e€ect of these deformation-dependent third-order terms of element nodal forces on the buckling load and post-buckling behavior was investigated for doubly symmetric beams with solid sections. It was reported in [34] that the third-order term of twist rate is the dominant third-order term of the element nodal forces and may not be negligible for the buckling and post-buckling analysis. In [35] the formulation proposed in [32] for beams with solid sections was extended to the thin-walled beam with monosymmetric open section. In [35] the third-order term of the twist rate is considered in element nodal forces. However, the e€ect of deformation-dependent third-order terms of element nodal forces on the buckling load and post-buckling behavior was not investigated in [35]. To the authors' knowledge, the e€ect of deformation-dependent third-order terms of element nodal forces on the buckling load and post- buckling behavior of thin-walled beam with open section has not been reported in the literature.

The doubly symmetric thin-walled beam has been extensively used in practice. Thus, a reliable and ef-

®cient formulation for geometric nonlinear analysis of doubly symmetrical thin-walled beam may be re- quired. The object of this paper is to present a CR TL ®nite element formulation for the geometric nonlinear analysis of doubly symmetrical thin-walled beams and to investigate the e€ects of deformation- dependent third-order terms of element nodal forces on the buckling load and post-buckling behavior of thin-walled beam.

The shear center and centroid are coincident for doubly symmetric beam but not for monosymmetric beam. Thus, the kinematics of the monosymmetric beam element is much more complicated than that of doubly symmetric beam element. The monosymmetric beam element proposed in [35] may be inecient for the analysis of doubly symmetric beams. Moreover, the beam element proposed in [35] cannot be used to investigate the e€ect of third-order terms on the buckling load and post-buckling behavior of thin-walled beam. However, even beam elements proposed in [34,35] are not suitable for the present study, it seems that the formulations given in [34,35] can be adapted for doubly symmetric thin-walled beams. Thus, the for- mulations of beam elements proposed in [34,35] are modi®ed and employed here. All third-order terms of the element nodal forces, which are relevant to the twist rate, the rate of twist rate and curvature of the beam axis, are retained here.

An incremental-iterative method based on the Newton±Raphson method combined with constant arc length of incremental displacement vector is employed for the solution of nonlinear equilibrium equations.

The zero value of the tangent sti€ness matrix determinant of the structure is used as the criterion of the buckling state. A parabolic interpolation method of the arc length is used to ®nd the buckling load. An inverse power method for the solution of the generalized eigenvalue problem is used to ®nd the corre- sponding buckling mode. In order to initiate the secondary path, at the bifurcation point a perturbation displacement proportional to the ®rst buckling mode is added. Numerical examples are presented to demonstrate the accuracy and eciency of the proposed method and to investigate the e€ect of third-order terms of the element nodal forces on the buckling load and post-buckling behavior of doubly symmetric thin-walled beams.

2. Finite element formulation

The formulations given in [34,35] are modi®ed and employed here.

2.1. Basic assumptions

The following assumptions are made in the derivation of the beam element behavior.

1. The beam is prismatic and slender, and the Euler±Bernoulli hypothesis is valid.

2. The cross-section of the beam is doubly symmetric.

3. The unit extension of the centroid axis of the beam element is uniform.

4. The cross-section of the beam element does not deform in its own plane and strains within this cross- section can be neglected.

5. The out-of-plane warping of the cross-section is the product of the twist rate of the beam element and the Saint Venant warping function for a prismatic thin-walled beam of the same cross-section.

6. The deformation displacements of the beam element are small.

In this study, Prandtl's membrane analogy and the Saint Venant torsion theory [37] are used to obtain an approximate Saint Venant warping function for a prismatic thin-walled beam.

2.2. Coordinate systems

In this paper, a CR TL formulation is adopted. In order to describe the system, we de®ne four sets of right-handed rectangular Cartesian coordinate systems:

1. A ®xed global set of coordinates, X_i^G (i ˆ 1; 2; 3) (see Fig. 1); the nodal coordinates, displacements, and rotations, and the stiffness matrix of the system are de®ned in these coordinates.

2. Element cross-section coordinates, x^S_i (i ˆ 1; 2; 3) (see Fig. 1); a set of element cross-section coordinates is associated with each cross-section of the beam element. The origin of this coordinate system is rigidly tied to the shear center of the cross-section. The x^S₁ axis is chosen to coincide with the normal of the un- warped cross-section and the x^S₂and x^S₃axes are chosen to be the principal directions of the cross-section.

3. Element coordinates; x_i(i ˆ 1; 2; 3) (see Fig. 1), a set of element coordinates is associated with each ele- ment, which is constructed at the current con®guration of the beam element. The origin of this coordi- nate system is located at node 1, and the x₁ axis is chosen to pass through two shear centers of end sections of the element; the x2and x3 axes are chosen to be the principal directions of the cross-section in the undeformed state. Note that this coordinate system is a local coordinate system not a moving co- ordinate system. The deformations, internal nodal forces and stiffness matrix of the elements are de®ned

Fig. 1. Coordinate systems.

in terms of these coordinates. In this paper the element deformations are determined by the rotation of element cross-section coordinate systems relative to this coordinate system.

4. Load base coordinates, X_i^P (i ˆ 1; 2; 3); a set of load base coordinates is associated with each con®gura- tion-dependent moment. The origin of this coordinate system is chosen to be the node where the con-

®guration-dependent moment is applied. The mechanism for generating con®guration dependent moment is described in these coordinates, and the corresponding external load and load stiffness matrix are de®ned in terms of these coordinates.

In this paper, the symbol { } denotes the column matrix. The relations among the global coordinates, element cross-section coordinates, element coordinates and load base coordinates may be expressed by

X^Gˆ AGSx^S; X^Gˆ AGEx; X^Gˆ AGPX^P; …1†

where X^Gˆ fX₁^G; X₂^G; X₃^Gg; x^S ˆ fx^S₁; x^S₂; x^S₃g, x ˆ fx1; x2; x3g, and X^Pˆ fX₁^P; X₂^P; X₃^Pg; AGS; AGE, and AGP are the matrices of direction cosines of the element cross-section coordinate system, element coordinate system, and load base coordinate system, respectively.

2.3. Rotation vector

For the convenience of later discussion, the term `rotation vector' is used to represent a ®nite rotation.

Fig. 2shows that a vector b which as a result of the application of a rotation vector /a is transported to the new position b. The relation between b and b may be expressed as [38]

b ˆ cos /b ‡ …1 cos /†…a
b†a ‡ sin /…a  b†; …2†

where / is the angle of counterclockwise rotation, and a is the unit vector along the axis of rotation.

2.4. Kinematics of beam element

The deformations of the beam element are described in the current element coordinate system. From the kinematic assumptions made in this paper, the deformations of the beam element may be determined by the displacements of the shear center axis of the beam element, orientation of the cross-section (element cross- section coordinates), and the out-of-plane warping of the cross-section. In this study only the doubly symmetric cross-section is considered. Thus the shear center and centroid of the cross-section are coinci- dent. Let Q (Fig. 1) be an arbitrary point in the beam element, and P be the point corresponding to Q on

Fig. 2. Rotation vector.

the shear center axis. The position vector of point Q in the undeformed and deformed con®gurations may be expressed as

r0ˆ xe1‡ ye2‡ ze3 …3†

and

r ˆ xc…x†e1‡ v…x†e2‡ w…x†e3‡ h1;xxe^S₁‡ ye^S₂‡ ze^S₃; …4†

where xc…x†; v…x†, and w…x† are the x1; x2 and x3 coordinates of point P, respectively, in the deformed con®guration, x ˆ x…y; z† is the Saint Venant warping function for a prismatic beam of the same cross- section, and ei and e^S_i (i ˆ 1; 2; 3) denote the unit vectors associated with the xi and x^S_i axes, respectively.

Note that e_i and e^S_i are coincident in the undeformed state. Here, the triad e^S_i in the deformed state is assumed to be achieved by the successive application of the following two rotation vectors to the triad e_i:

h_nˆ h_nn …5†

and

htˆ h1t; …6†

where

n ˆ f0; h2=…h²₂‡ h²₃†¹⁼²; h3=…h²₂‡ h²₃†¹⁼²g; …7†

t ˆ fcos hn; h3; h2g; …8†

cos hnˆ …1 h²₂ h²₃†¹⁼²; …9†

h₂ˆ dw…x†

ds ˆ dw…x†

dx dx

dsˆ w⁰

1 ‡ ec; …10†

h3ˆdv…x†

ds ˆdv…x†

dx dx dsˆ v⁰

1 ‡ e_c; …11†

ec ˆos

ox 1; …12†

in which n is the unit vector perpendicular to the vectors e1and e^S₁, and t is the tangent unit vector of the deformed shear center axis of the beam element. Note that the orientation of e^S₁coincides with that of t. h1is the rotation about vector t. hnis the angle measured from x1axis to vector t; ecis the unit extension of the shear center axis and s is the arc length of the deformed shear center axis measured from node 1 to point P.

In this paper, the symbol … †⁰ denotes … †_;xˆ o… †=ox.

Using Eqs. (2)±(8), the relation between the vectors ei and e^S_i (i ˆ 1; 2; 3) in the element coordinate system may be obtained as [32]

e^S_i ˆ Re_i; …13†

where R is the so-called rotation matrix. The rotation matrix is determined by h_i(i ˆ 1; 2; 3). Thus, h_iare called the rotation parameters in this study.

Let h ˆ fh₁; h₂; h₃g be the column matrix of rotation parameters, dh be the variation of h. The triad e^S_i (i ˆ 1; 2; 3) corresponding to h may be rotated by a rotation vector d/ ˆ fd/₁; d/₂; d/₃g to reach their new positions corresponding to h ‡ dh [32]. When h2and h3are much smaller than unity, the relationship be- tween dh and d/ may be approximated by [32]

dh ˆ 1 h3=2 h2=2

h₃ 1 0

h₂ 0 1

2 4

3

5d/ ˆ T ¹d/: …14†

The relationship among xc…x†; v…x†; w…x†, and x may be given as xc…x† ˆ u1‡

Z _x

0 ‰…1 ‡ ec†² v²_;x w²_;xŠ¹⁼²dx; …15†

where u1 is the displacement of node 1 in the x1 direction. Note that due to the de®nition of the element coordinate system, the value of u1is equal to zero. However, the variation of u1is not zero. Making use of Eq. (15), one obtains

` ˆ L ‡ u2 u1ˆ xc…L† xc…0† ˆ Z _L

0 ‰…1 ‡ ec†² v²_;x w²_;xŠ¹⁼²dx; …16†

in which ` is the current chord length of the shear center axis of the beam element, and L is the length of the undeformed beam axis, and u2 is the displacement of node 2in the x1 direction. Making use of the as- sumption of uniform unit extension, ec may be calculated using Eq. (16).

Here, the lateral de¯ections of the shear center axis, v…x† and w…x†, and the rotation about the shear center axis, h₁…x†, are assumed to be the Hermitian polynomials of x. v…x†; w…x† and h₁…x† may be expressed by

v…x† ˆ fN1; N2; N3; N4g^tfv1; v⁰₁; v2; v⁰₂g ˆ N^t_bub;

w…x† ˆ fN1; N2; N3; N4g^tfw1; w⁰₁; w2; w⁰₂g ˆ N^t_cuc; h1…x† ˆ fN1; N2; N3; N4g^tfh11; b₁; h12; b₂g ˆ N^t_dud;

…17†

where vjand wj(j ˆ 1; 2; 3) are the nodal values of v and w at nodes j, respectively, v⁰_jand w⁰_j(j ˆ 1; 2) are the nodal values of v;xand w;xat nodes j, respectively, and h1jand b_j(j ˆ 1; 2) are nodal values of h1; h1;xat nodes j, respectively. Note that, due to the de®nition of the element coordinates, the values of vj and wj(j ˆ 1; 2) are zero. However, their variations are not zero. Ni(i ˆ 1±4) are the shape functions and are given by

N1ˆ1

4…1 n†²…2 ‡ n†; N2ˆL

8…1 n²†…1 n†;

N3ˆ1

4…1 ‡ n†²…2 n†; N4ˆL

8… 1 ‡ n²†…1 ‡ n†;

…18†

where

n ˆ 1 ‡2x

L : …19†

The axial displacements of the shear center axis may be determined from the lateral de¯ections and the unit extension of the shear center axis using Eq. (15).

If x, y and z in Eq. (3) are regarded as the Lagrangian coordinates, the Green strain e₁₁; e₁₂and e₁₃are given by [39]

e₁₁ˆ1

2…r^t_;xr_;x 1†; e₁₂ˆ1

2r^t_;xr_;y; e₁₃ˆ1

2r^t_;xr_;z: …20†

Substituting Eqs. (4), (8)±(13) into Eq. (20) and retaining all terms up to the second-order and the third- order terms, which are relevant to the retained third-order terms of nodal forces, yield

e11ˆ ec yv;xx zw;xx‡ xh1;xx‡1

2e²_c‡ xech1;xx‡1

2…y²‡ z²†h²_1;x‡1 2v²_;xx‡1

2z²w²_;xx‡1 2x²h²_1;xx

yh1w;xx‡ zh1v;xx‡ yzv;xxw;xx yxv;xxh1;xx zxw;xxh1;xx yxh²_1;xw;xx zxh²_1;xv;xx; …21†

e12ˆ1

2…x;y z†h1;x‡1

4z…v;xw;xx w;xv;xx† ‡1

2x;yech1;x

‡ …x yx;y†h1;xv;xx zx;yh1;xw;xx

‡ xx;yh1;xh1;xx

; …22†

e13ˆ1

2…x;z‡ y†h1;x‡1

4y…w;xv;xx v;xw;xx† ‡1

2‰x;zech1;x‡ …x zx;z†h1;xw;xx yx;zh1;xv;xx

‡ xx;zh1;xh1;xxŠ: …23†

The underlined terms in Eq. (21) are the retained third-order terms of strains.

2.5. Nodal parameters and forces

The element employed here has two nodes with seven degrees of freedom per node. Two sets of element nodal parameters termed `explicit nodal parameters' and `implicit nodal parameters' are employed. The explicit nodal parameters of the element are used for the assembly of the system equations from the element equations. They are chosen to be u_ij (u_1jˆ u_j; u_2jˆ v_j; u_3jˆ w_j), the x_i (i ˆ 1; 2; 3) components of the translation vectors u_jat node j …j ˆ 1; 2), /_ij, the x_i (i ˆ 1; 2; 3) components of the rotation vectors /_j at node j (j ˆ 1; 2), and b_j, the twist rate of the shear center axis at node j. Here, the values of /_jare reset to zero at the current con®guration. Thus, d/_ij, the variation of /_ij, represents in®nitesimal rotations about the xiaxes [32], and the generalized nodal forces corresponding to d/_ij are mij, the conventional moments about the xiaxes. The generalized nodal forces corresponding to duij, the variations of uij, are fij, the forces in the xidirections. The generalized nodal forces corresponding to db_j, the variations of b_j, are bimoment Bj.

The implicit nodal parameters of the element are used to determine the deformation of the beam element. They are chosen to be uij, the xi (i ˆ 1; 2; 3) components of the translation vectors uj at node j (j ˆ 1; 2), h1j, b_j; v⁰_j, and w⁰_j (j ˆ 1; 2) de®ned in Eq. (17). Let h^_1j; h^_2j and h^_3j (j ˆ 1; 2; 3) denote h1j, w⁰_j and v⁰_j, respectively. The generalized nodal forces corresponding to duij; dh^_ij and db_j are fij; m^h_ij and Bj, the forces in the xi directions, the generalized moments, and bimoments, respectively. Note that m^h_ij are not conventional moments, because dh^_ij are not in®nitesimal rotations about the xi axes at deformed state.

In view of Eqs. (10) and (14), the relations between the variation of the implicit and explicit nodal parameters may be expressed as

dq_hˆ du1

dh^₁ du2

dh^₂ db 8>

>>

><

>>

>>

: 9>

>>

>=

>>

>>

;

ˆ

I₃ 0 0 0 0

T_b1 T_a1 T_b1 0 0

0 0 I₃ 0 0

Tb2 0 Tb2 Ta2 0

0^t 0^t 0^t 0^t I2

2 66 66 4

3 77 77 5

du1

d/₁ du2

d/₂ db 8>

>>

><

>>

>>

: 9>

>>

>=

>>

>>

;

ˆ Th/dq; …24†

Tbjˆ 0 0 0

h2j=L 0 0 h3j=L 0 0 2

4

3

5; Tajˆ 1 h3j=2 h2j=2 h3j 1 ‡ ec 0 h2j 0 1 ‡ ec

2 4

3

5 …j ˆ 1; 2†; …25†

where du_jˆ fdu_j; dv_j; dw_jg, dh^_j ˆ fdh_1j; dw⁰_j; dv⁰_jg, d/_jˆ fd/_1j; d/_2j; d/_3jg (j ˆ 1; 2) and db ˆ dbf ₁; db₂g;

I₂and I₃are the identity matrices of order 2  2and 3  3, respectively, and 0 and 0 are the zero matrices of order 3  3 and 3  2, respectively.

Let f ˆ ff1; m1; f2; m2; Bg; fhˆ ff1; m^h₁; f2; m^h₂; Bg, where fjˆ ff1j; f2j; f3jg; mjˆ fm1j; m2j; m3jg, m^h_j ˆ fm^h_1j; m^h_2j; m^h_3jg (j ˆ 1; 2), and B ˆ fB1; B2g, denote the internal nodal force vectors corresponding to the

variation of the explicit and implicit nodal parameters, dq and dq_h, respectively. Using the contragradient law [40] and Eq. (24), the relation between f and fh, may be given by

f ˆ T^t_h/f_h: …26†

The global nodal parameters for the structural system corresponding to the element local nodes j (j ˆ 1; 2) should be consistent with the element explicit nodal parameters. Thus, they are chosen to be U_ij, the Xi (i ˆ 1; 2; 3) components of the translation vectors Uj at node j (j ˆ 1; 2), Uij, the Xi (i ˆ 1; 2; 3) components of the rotation vectors Ujat nodes j (j ˆ 1; 2), and b_j, the twist rate of the shear center axis at node j. Here, the values of Ujare reset to zero at the current con®guration. Thus, dUij, the variations of Uij, represent in®nitesimal rotations about the Xiaxes [32], and the generalized nodal forces corresponding to dUij are the conventional moments about the Xi axes. The generalized nodal forces corresponding to dUij, the variation of Uij, are the forces in the Xidirections. The generalized nodal forces corresponding to db_j, the variation of b_j, are Bj.

2.6. Element nodal force vector

The element nodal force vector fh(Eq. (26)) corresponding to the implicit nodal parameters is obtained from the virtual work principle in the current element coordinates. It should be mentioned again that the element coordinate system is a local coordinate system not a moving coordinate system. Thus, a standard procedure is used here for the derivation of fh. For convenience, the implicit nodal parameters are divided into four generalized nodal displacement vectors u_i(i ˆ a; b; c; d), where

uaˆ fu1; u2g …27†

and ub; uc, and ud are de®ned in Eq. (17).

The generalized force vectors corresponding to dui, the variation of ui (i ˆ a; b; c; d), are faˆ ff11; f12g; fbˆ ff21; m^h₃₁; f22; m^h₃₂g;

f_cˆ ff₃₁; m^h₂₁; f₃₂; m^h₂₂g; f_dˆ fm^h₁₁; B₁; m^h₁₂; B₂g: …28†

The virtual work principle requires that du^t_afa‡ du^t_bfb‡ du^t_cfc‡ du^t_dfdˆ

Z

V…r11de11‡ 2r12de12‡ 2r13de13†dV ; …29†

where V is the volume of the undeformed beam, de1j(j ˆ 1; 2; 3) are the variations of e1jin Eqs. (21)±(23), respectively, with respect to the implicit nodal parameter. r1j (j ˆ 1; 2; 3) are the second Piola±Kirchhoff stress. For linear elastic material, the following constitutive equations are used:

r11ˆ Ee11; r12ˆ 2Ge12; and r13ˆ 2Ge13; …30†

in which E is Young's modulus and G is the shear modulus.

If the element size is chosen to be small enough, the values of the rotation parameters of the deformed element de®ned in the current element coordinate system may always be much smaller than unity. Thus the higher-order terms of rotation parameters in the element internal nodal forces may be neglected. However, in order to include the nonlinear coupling among the bending, twisting, and stretching deformations, the terms up to the second-order of rotation parameters and their spatial derivatives are retained in element internal nodal forces by consistent second-order linearization of Eq. (29). However the values of h_1;x; h_1;xx; v_;xx, and w_;xx in Eqs. (21)±(23) are deformation dependent, not element size dependent. Thus their values may not always be much smaller than unity and their third- order terms may not be negligible. Here, the third-order terms of h1;x; h1;xx; v;xx; and w;xx are also retained in Eq. (29).

From Eqs. (21)±(23), (27)±(30), we may obtain faˆ AELec 1

" 

‡3 2ec



‡1 2EIp

Z

h²_1;xdx ‡1 2EIy

Z

w²_;xxdx ‡1 2EIz

Z v²_;xxdx

‡3 2EIx

Z

h²_1;xxdx‡EIxyz

Z

h1;xxv;xxw;xxdx

#

Ga; …31†

fbˆ EIz…1 ‡ ec† Z

N⁰⁰_bv;xxdx ‡ f12Gb‡ E…Iz Iy† Z

N⁰⁰_bh1w;xxdx ‡1 2GJ

Z

…N⁰⁰_bh1;xw;x N⁰_bh1;xw;xx†dx

‡ 3EI_xyz Z

N⁰⁰_bh_1;xxw_;xxdx‡3 2EK_xz

Z

N⁰⁰_bh²_1;xxv_;xxdx ‡1 2EK_z

Z

N⁰⁰_bv³_;xxdx

‡3 2EKyz

Z

N⁰⁰_bw²_;xxv;xxdx ‡ Cb

Z

N⁰⁰_bh²_1;xv;xxdx ‡1 2EKyz

Z

N⁰⁰_bh²_1;xw;xxdx; …32†

fcˆ EIy…1 ‡ ec† Z

N⁰⁰_cw;xxdx ‡ f12Gc‡ E…Iz Iy† Z

N⁰⁰_ch1v;xxdx ‡1 2GJ

Z

…N⁰_ch1;xv;xx N⁰⁰_ch1;xv;x†dx

‡ 3EIxyz

Z

N⁰⁰_ch1;xxv;xxdx‡3 2EKxy

Z

N⁰⁰_ch²_1;xxw;xxdx ‡1 2EKy

Z

N⁰⁰_cw³_;xxdx

‡3 2EK_yz

Z

N⁰⁰_cv²_;xxw_;xxdx ‡ C_c Z

N⁰⁰_ch²_1;xw_;xxdx ‡1 2EK_yz

Z

N⁰⁰_ch²_1;xv_;xxdx; …33†

fdˆ ‰GJ ‡ EIpecŠ Z

N⁰_dh1;xdx ‡ EIx…1 ‡ 3ec† Z

N⁰⁰_dh1;xxdx ‡ E…Iz Iy† Z

Ndv;xxw;xxdx

‡1 2GJ

Z

N⁰_d…w;xv;xx v;xw;xx†dx ‡ 3EIxyz

Z

N⁰⁰_dv;xxw;xxdx‡1 2EKI

Z

N⁰_dh³_1;xdx

‡1 2EKx

Z

N⁰⁰_dh³_1;xxdx ‡ 2EIxyz

Z

Ndh1;xx…w²_;xx v²_;xx†dx ‡ Cb

Z

N⁰_dv²_;xxh1;xdx ‡ Cc

Z

N⁰_dw²_;xxh1;xdx

‡3 2EKxz

Z

N⁰⁰_dv²_;xxh1;xxdx ‡3 2EKxy

Z

N⁰⁰_dw²_;xxh1;xxdx ‡ Cd

Z

…N⁰_dh²_1;xxh1;x‡ N⁰⁰_dh²_1;xh1;xx†dx; …34†

where Gaˆ1

Lf 1; 1g; Gbˆ Z

N⁰_bv;xdx; Gcˆ Z

N⁰_cw;xdx; …35†

Iyˆ Z

z²dA; Izˆ Z

y²dA; Ipˆ Iy‡ Iz; Ixyzˆ Z

xyzdA; Ixˆ Z

x²dA;

Cbˆ1

2E…Kz‡ Kyz‡ 4Ixyz† ‡ GJz; Ccˆ1

2E…Ky‡ Kyz 4Ixyz† ‡ GJy; Cdˆ1

2E…Kxy‡ Kxz† ‡ GJx; Ky ˆ Z

z⁴dA; Kyzˆ Z

y²z²dA; Kzˆ Z

y⁴dA;

K1ˆ Ky‡ Kz‡ 2Kyz; Kxˆ Z

x⁴dA; Kxy ˆ Z

x²z²dA; Kxzˆ Z

x²y²dA;

J ˆ Z

h…

z ‡ x;y†²‡ …y ‡ x;z†²i

dA; Jy ˆ Z

…zx;z

h x†²‡ z²x²_;yi dA;

Jzˆ Z

h…x

yx;y†²‡ y²x²_;zi

dA; Jxˆ Z

x²x²_;y

‡ x²_;z

dA …36†

in which the range of integrations for the integralR

… †dx in Eqs. (31)±(36) is from 0 to L, A is the cross- sectional area, and Nj(j ˆ b; c; d) are given in Eq. (18). The underlined terms in Eqs. (31)±(34) are the third- order terms of h1;x; h1;xx; v;xx, and w;xx.

The element nodal force vector f (Eq. (26)) corresponding to the explicit nodal parameters may be obtained from Eqs. (26) and (31)±(34). Note that only the terms up to the second-order of nodal parameters and the third-order terms of h1;x; h1;xx; v;xx, and w;xx are retained in fh. Thus, the corresponding f in Eq. (26) may be rewritten as

f ˆ f_h‡ T ^t_h/

I₁₄

f¹_h; …37†

where f¹_hare the ®rst-order terms of nodal parameters of fh, and I14is the identity matrix of order 14  14.

2.7. Element tangent sti€ness matrices

The element tangent sti€ness matrix corresponding to the explicit nodal parameters (referred to as explicit tangent sti€ness matrix) may be obtained by di€erentiating the element nodal force vector f in Eq. (37) with respect to explicit nodal parameters. Using Eqs. (24) and (37), we obtain

k ˆof oqˆ of

oq_h oq_h

oq ˆ kh

h ‡ T ^t_h/

I14

k⁰_h‡ Hi

Th/; …38†

where k_hˆ of_h=oq_his the tangent sti€ness matrix corresponding to implicit nodal parameters (referred to as implicit tangent sti€ness matrix), k⁰_h are the zeroth-order terms of nodal parameters of k_h, and H is a unsymmetrical matrix and is given by

H ˆ

0 hb1 0 hb2 0

h^t_b1 h_a1 h^t_b1 0 0

0 h_b1 0 h_b2 0

h^t_b2 0 h^t_b2 ha2 0

0^t 0^t 0^t 0^t 02

2 66 66 66 4

3 77 77 77 5

; …39†

hbj ˆ 0 …1=L†m^h_2j …1=L†m^h_3j

0 0 0

0 0 0

2 4

3

5; hajˆ

0 m^h_3j m^h_2j 0 0 …1=2†m^h_1j 0 …1=2†m^h_1j 0 2

64

3

75 …j ˆ 1; 2; 3†; …40†

where 02; 0 and 0 are the zero matrices of order 2  2 , 3  3 and 3  2, respectively.

Using the direct sti€ness method, the implicit tangent sti€ness matrix kh may be assembled by the submatrices

k_ijˆof_i

ouj; …41†

where fi(i ˆ a; b; c; d) are de®ned in Eqs. (31)±(34) and ui (i ˆ a; b; c; d) are de®ned in Eqs. (17) and (27).

Note that kij are the symmetric matrices. The explicit form of kij may be expressed as kaaˆ AEL…1 ‡ 3ec†GaG^t_a;

kabˆ Ga AEG^t_b‡ EIz

Z

N⁰⁰_b^tv;xxdx ‡ EIxyz

Z

N⁰⁰_b^th1;xxw;xxdx

!

;

k_acˆ G_a AEG^t_c‡ EI_y Z

N⁰⁰_b^tw_;xxdx ‡ EI_xyz Z

N⁰⁰_b^th_1;xxv_;xxdx

!

;

kadˆ Ga EIp

Z

N⁰_d^th1;xdx‡ 3EIx

Z

N⁰⁰_d^th1;xxdx ‡ EIxyz

Z

N⁰⁰_d^tv;xxw;xxdx

!

;

k_bbˆ EI_z…1 ‡ e_c† Z

N⁰⁰_bN⁰⁰_b^tdx ‡ f₁₂ Z

N⁰_bN⁰_b^tdx

‡3 2EKxz

Z

N⁰⁰_bN⁰⁰_b^th²_1;xxdx ‡3 2EKz

Z

N⁰⁰_bN⁰⁰_b^tv²_;xxdx ‡3 2EKyz

Z

N⁰⁰_bN⁰⁰_b^tw²_;xxdx ‡ Cb

Z

N⁰⁰_bN⁰⁰_b^th²_1;xdx;

kbcˆ E…Iz Iy† Z

N⁰⁰_bN⁰⁰_c^th1dx ‡ 3EIxyz

Z

N⁰⁰_bN⁰⁰_c^th1;xxdx ‡1 2GJ

Z

…N⁰⁰_bN⁰_c^t N⁰_bN⁰⁰_c^t†h1;xdx

‡ EKyz

Z

N⁰⁰_bN⁰⁰_c^t 3v;xxw;xx‡1 2h²_1;x

 

dx;

kbdˆ E…Iz Iy† Z

N⁰⁰_bN^t_dw;xxdx ‡1 2GJ

Z

…N⁰⁰_bN⁰_d^tw;x N⁰_bN⁰_d^tw;xx†dx ‡ 3EIxyz

 Z

N⁰⁰_bN⁰⁰_d^tw;xxdx‡3EKxz

Z

N⁰⁰_bN⁰⁰_b^th1;xxv;xxdx ‡ 2Cb

Z

N⁰⁰_bN⁰_d^th1;xv;xxdx ‡ EKyz

Z

N⁰⁰_bN⁰_d^th1;xw;xxdx;

kccˆ EIy…1 ‡ ec† Z

N⁰⁰_cN⁰⁰_c^tdx ‡ f12

Z

N⁰_cN⁰_c^tdx

‡3 2EKxy

Z

N⁰⁰_cN⁰⁰_c^th²_1;xxdx ‡3 2EKy

Z

N⁰⁰_cN⁰⁰_c^tw²_;xxdx ‡3 2EKyz

Z

N⁰⁰_cN⁰⁰_c^tv²_;xxdx ‡ Cc

Z

N⁰⁰_cN⁰⁰_c^th²_1;xdx;

kcd ˆ E…Iz Iy† Z

N⁰⁰_cN^t_dv;xxdx ‡1 2GJ

Z

…N⁰_cN⁰_d^tv;xx N⁰⁰_cN⁰_dtv;x†dx ‡ 3EIxyz

 Z

N⁰⁰_cN⁰⁰_d^tv;xxdx‡3EKxy

Z

N⁰⁰_cN⁰⁰_d^th1;xxw;xxdx ‡ 2Cc

Z

N⁰⁰_cN⁰_d^th1;xw;xxdx ‡ EKyz

Z

N⁰⁰_cN⁰_d^th1;xv;xxdx;

kdd ˆ ‰GJ ‡ EIpecŠ Z

N⁰_dN⁰_d^tdx ‡ EIx…1 ‡ 3ec† Z

N⁰⁰_dN⁰⁰_d^tdx

‡3 2EKI

Z

N⁰_dN⁰_d^th²_1;xdx ‡3 2EKx

Z

N⁰⁰_dN⁰⁰_d^th²_1;xxdx ‡ 2EIxyz

Z

NdN⁰⁰_d^t…w²_;xx v²_;xx†dx

‡C_b Z

N⁰_dN⁰_d^tv²_;xxdx ‡ C_c Z

N⁰_dN⁰_dtw²_;xxdx ‡3 2EK_xz

Z

N⁰⁰_dN⁰⁰_d^tv²_;xxdx

‡3 2EKxy

Z

N⁰⁰_dN⁰⁰_dtw²_;xxdx ‡ Cd

Z

N⁰_dN⁰_d^th²_1;xx‡ N⁰⁰_dN⁰⁰_dth²_1;x‡ 2…N⁰_dN⁰⁰_d^t‡ N⁰⁰_dN⁰_d^t†h1;xh1;xx

h i

dx; …42†

where the underlined terms are the second-order terms of h1;x; h1;xx; v;xx, and w;xx.

The element tangent sti€ness matrix referred to the global coordinate system is obtained by using the standard coordinate transformation.

2.8. Load sti€ness matrix

Di€erent ways for generating con®guration-dependent moment were proposed in the literature [8,10,33,41]. Here, for simplicity, only the conservative moments generated by conservative force or

forces (with ®xed directions) are considered, and the ways for generating conservative moment proposed in [33] are employed here. For completeness, a brief description of the ways for generating conservative moment is given here. In this study, a set of load base coordinates X_i^P (i ˆ 1; 2; 3) associated with each con®guration-dependent moment is constructed at the current con®guration. The mechanism for gen- erating con®guration-dependent moment is described in these coordinates, and the corresponding ex- ternal load and load sti€ness matrix [42] are de®ned in terms of these coordinates. Unless stated otherwise, all vectors and matrices in this section are referred to these coordinates. Note that this co- ordinate system is just a local coordinate system constructed at the current con®guration, not a moving coordinate system. Thus, it is regarded as a ®xed coordinated system in the derivation of the load sti€ness matrix.

The ®rst way for generating con®guration-dependent moment may described as follows.

Consider a sphere of radius R whose centroid is rigidly connected with the structure at node O as shown in Fig. 3. Two strings wound around a great circle of the sphere are acted upon by forces of magnitude P.

Thus, the strings are always tangent to the sphere. The great circle and the forces are on the same plane at the undeformed con®guration of the structure. However, the great circle and the forces are generally not on the same plane when the structure is deformed. The origin of the load base coordinate system is chosen to be located at the node O. The X₁^Paxis is chosen to coincide with the normal of the plane of the great circle, and the X₂^P and X₃^P axes lie in the plane of the great circle.

Let A denote the contact point of the force P and the great circle. Because P is tangent to the sphere, P is perpendicular to the line OA. Let e_Adenote the unit vector in the direction of line OA. e_Amay be expressed by

e_Aˆ a=…a^ta†¹⁼²; …43†

a ˆ e^P_p n^P ˆ f0; `3; `2g; …44†

where e^P_p ˆ f`1; `2; `3g is the unit vector in the direction of P, and n^P is the unit normal of the plane of the great circle. Note that n^P coincides with e^P₁ ˆ f1; 0; 0g, the unit vector associated with the X₁^P axis.

The external moment vector at node O generated by the above-mentioned mechanism may be expressed by

M ˆ M eA e^P_p; …45†

where M ˆ 2RP is the magnitude of the moment. The corresponding load sti€ness matrix kp may be ex- pressed as

kpˆ M…`<sup>2</sup><sub>2</sub>‡ `²₃† ¹⁼²kpa‡ M…`<sup>2</sup><sub>2</sub>‡ `²₃† ³⁼²kpb; …46†

Fig. 3. Mechanism for generating con®guration-dependent moment.

where

kpaˆ

0 `1`3 `1`2

0 `2`3 `<sup>2</sup><sub>1</sub>‡ `²₃ 0 `<sup>2</sup><sub>1</sub> `²₂ `2`3

2 4

3

5; …47†

k_pbˆ

0 `1`3…`<sup>2</sup><sub>2</sub>‡ `²₃† `1`2…`<sup>2</sup><sub>2</sub>‡ `²₃† 0 `<sup>2</sup><sub>1</sub>`2`3 `²₁`<sup>2</sup><sub>2</sub> 0 `²₁`<sup>2</sup><sub>3</sub> `²₁`2`3

2 64

3

75: …48†

Three special cases shown in Fig. 4 are considered here. Following [10], they are referred to as quasitan- gential (QT) moments of the ®rst and second types, and semitangential (ST) moment. The load sti€ness matrices corresponding to QT and ST moment at the con®gurations shown in Fig. 4 may be obtained from Eqs. (46)±(48) and given by

k^QT1_p ˆ M

0 0 0

0 0 0

0 1 0

2 4

3

5; …49†

k^QT2_p ˆ M

0 0 0 0 0 1 0 0 0 2

4

3

5; …50†

k^ST_p ˆM 2

0 0 0

0 0 1

0 1 0

2 4

3

5: …51†

The second way for generating con®guration-dependent moment may described as follows.

Consider a rigid arm of length R which end is rigidly connected with the structure at node O as shown in Fig. 5. The other end of the rigid arm is acted upon by a conservative force (with a ®xed direction) of magnitude P. The origin of the load base coordinates X_i^P (i ˆ 1; 2; 3) is chosen to be located at the node O.

The X₁^P axis is chosen to coincide with the axis of the rigid arm, and the X₂^Pand X₃^Paxes are perpendicular to the rigid arm.

The external moment vector at node O generated by the above-mentioned mechanism may be expressed by

M ˆ RP t^P e^P_p; …52†

where e^P_pˆ f`1; `2; `3g is the unit vector in the direction of P, and t^P is the unit vector in the axial direction of the rigid arm. Note that t^P coincides with e^P₁ ˆ f1; 0; 0g, the unit vector associated with the X₁^P axis. The corresponding load sti€ness matrix may be expressed as

Fig. 4. Quasitangential (QT) moment and semitangential (ST) moment.

kpˆ RP 0 `2 `3

0 `1 0

0 0 `1

2 4

3

5: …53†

The load sti€ness matrix referred to the global coordinate system is obtained by using the standard coordinate transformation and may be expressed by

k^G_p ˆ AGPkpA^t_GP; …54†

where A_GP is the transformation matrix given in Eq. (1).

2.9. Equilibrium equations

The nonlinear equilibrium equations may be expressed by

W ˆ F kP ˆ 0; …55†

where W is the unbalanced force between the internal nodal force F and the external nodal force kP, where k is the loading parameter, and P is a reference loading. Note that P may require to be updated at each iteration if the applied load is con®guration dependent. F is assembled from the element nodal force vectors, which are calculated using Eqs. (31)±(34) and (37) ®rst in the current element coordinates and then transformed from current element coordinate system to global coordinate system before assemblage using standard procedure.

In this paper, a weighted Euclidean norm of the unbalanced force is employed as the error measure for the equilibrium iterations, and is given by

kWk jkj 

p 6 eN tol; …56†

where N is the number of equilibrium equations; e_tol is a prescribed value of error tolerance.

2.10. Criterion of the buckling state

Here, the zero value of the tangent sti€ness matrix determinant is used as the criterion of the buckling state. The tangent sti€ness matrix of the structure is assembled from the element sti€ness matrix and load

Fig. 5. Mechanism for generating con®guration-dependent moment by an o€-axis load.