https://doi.org/10.1007/s11044-021-09777-x
Three-dimensional beam element for pre- and
post-buckling analysis of thin-walled beams in multibody
systems
J.B. Jonker1
Received: 24 April 2020 / Accepted: 6 January 2021 © The Author(s) 2021
Abstract This paper presents a three-dimensional beam element for stability analysis of elastic thin-walled open-section beams in multibody systems. The beam model is based on the generalized strain beam formulation. In this formulation, a set of independent deforma-tion modes is defined which are related to dual stress resultants in a co-rotadeforma-tional frame. The deformation modes are characterized by generalized strains or deformations, expressed as analytical functions of the nodal coordinates referred to the global coordinate system. A non-linear theory of non-uniform torsion of open-section beams is adopted for the derivation of the elastic and geometric stiffness matrices. Both torsional-related warping and Wagner’s stiffening torques are taken into account. Second order approximations for the axial elonga-tion and bending curvatures are included by addielonga-tional second order terms in the expressions for the deformations. The model allows to study the buckling and post-buckling behaviour of asymmetric thin-walled beams with open cross-section that can undergo moderately large twist rotations. The inertia properties of the beam are described using both consistent and lumped mass formulations. The latter is used to model rotary and warping inertias of the beam cross-section. Some validation examples illustrate the accuracy and computational efficiency of the new beam element in the analysis of the buckling and post-buckling be-haviour of thin-walled beams under various loads and (quasi)static boundary conditions. Finally, applications to multibody problems are presented, including the stability analysis of an elementary two-flexure cross-hinge mechanism.
Keywords Thin-walled beam element· Open section · Warping · Pre- and post-buckling · Deformation modes· Geometrically nonlinear · Spacar
1 Introduction
Thin-walled members found in multibody systems are often modelled as thin-walled beams. Beam models are computationally efficient and convenient for parametric analysis and
B
J.B. Jonkerj.b.jonker@utwente.nl
shape optimization of multibody geometries. Due to their geometric characteristics, thin-walled beams with open section exhibit complex structural behaviour, including cross-sectional warping and lateral-torsional buckling. In particular, pre-buckling deformations have a predominant influence on the stability of thin-walled beams, since these beams can undergo large flexible displacements and cross-sectional twist before buckling occurs with-out exceeding their yield limit [37]. In addition, at high angles of twist, longitudinal stresses induced by the axial shortening of the longitudinal fibres cause significant nonlinear stiff-ening effects [44]. An analysis of these phenomena must include at least the warping of the beam cross-section, a proper representation of the axial shortening effect and pre-buckling deflections.
Although the development of thin-walled beam finite elements has received a lot of atten-tion in recent decades, little effort has been devoted to their implementaatten-tion in flexible multi-body systems codes. Finite element models of flexible beams in a multimulti-body system can be formulated on the basis of a one-dimensional continuum description with strain measures expressed in a global inertial frame, or based on small deflection beam theory embedded in a co-rotational frame. Simo [42] developed a nonlinear beam theory based on the inertial frame approach that is known as geometrically exact beam theory. This beam theory is ca-pable of accounting for finite rotations and large deformations. Several authors developed nonlinear geometrically exact beam finite elements incorporating warping effects [22,31, 33,43] and cross-sectional deformation [16,21,23] to mention just a few. Saravia et al. [40] presented a geometrically exact thin-walled beam element with multibody capabilities for modelling of high aspect ratio composite wind turbine blades. Since the dominant nonlin-earity of flexible multibody systems is due to large rigid-body motions, the co-rotational formulation where the rigid-body motions are separated from deformations (e.g. see [14]) could be attractive for modelling of flexible multibody systems. A few articles dealing with co-rotational formulations in nonlinear buckling and post-buckling analyses for thin-walled open section beams have been published. Hsiao at al. [12,24] presented a consistent co-rotational total Lagrangian formulation for nonlinear analysis of thin-walled beams. They investigated the effects of deformation dependent third order terms of the element nodal forces on the buckling load and on the post buckling behaviour. Kim et al. [28,29] pre-sented a beam finite element formulation for post buckling analysis of thin-walled spatial frames using finite semitangential rotations. Battini and Pacoste [9] and Alsafadie et al. [3] presented 3D co-rotational finite elements for beams with arbitrary cross-section. Garcea et al. [17] constructed a co-rotational formulation to derive geometrically exact beam and shell models for nonlinear FEM analysis. Genoese et al. [18] exploited this method to develop a geometrically exact beam model with generic section including non-uniform warping due to torsion and shear. Rinchen et al. [39] presented a thin-walled beam finite element for arbitrary open cross-sections within the co-rotational frame work of OpenSees [1]. How-ever, when modelling flexible multibody systems with elastic and rigid bodies, conventional co-rotational formulations treat rigid bodies in the same way as elastic bodies with large stiffness. Thus they are not able of to model rigid-body dynamics exactly, yielding a less efficient formulation in terms of computational time and eliminating high frequency modes of deformation. For this reason conventional co-rotational formulations are still rarely used in flexible multibody dynamic analyses.
Co-rotational beam formulations bear much resemblance to the generalized strain beam formulation proposed by Besseling [10]. This formulation refers to the use of deforma-tion modes, by Argyris at al. [5] termed ‘natural modes’, which define elastic deformations as well as implicitly rigid-body motion of the element and served as the basis for the de-velopment of a finite element based beam formulation for rigid-flexible multibody system
analysis [25]. For each beam element, a fixed number of independent deformation modes is defined which are invariant under arbitrary rigid-body motions of the element. A set of generalized strains or deformations, corresponding to each of the deformation modes is then expressed as analytical functions of the nodal coordinates referred to the global co-ordinate system. Flexible elements are modelled by allowing non-zero deformations and specifying constitutive equations relating deformations and dual stress resultants. The equi-librium equations are formulated in the deformed configuration of the beam element and are capable of accounting for arbitrary large deformations. Geometric nonlinearities aris-ing from nonlinear couplaris-ings among axial elongation, bendaris-ing and torsional deformations are accounted for by modified deformations [27,34]. The deformations can be redefined in a co-rotational frame and related to the dual stress-resultants using existing beam models at different levels of sophistication ranging from elementary thin-walled beam theory [8] to relatively advanced formulations which take into account large torsion deformation and associated shortening effects. Attard [7] developed a nonlinear theory of non-uniform tor-sion for straight prismatic thin-walled open beams. The equations were presented in finite element form [6]. Trahair [44] developed a finite beam element for analysing large twist rotations of elastic thin-walled beams of general cross-section. Mohri et al. [36,37], em-ployed similar equations to those established by Attard [7] and developed a large torsion finite element model for thin-walled Bernoulli beams including flexural torsional coupling and pre-buckling deflection effects. Dourakopoulos and Sapountzakis [15] have extended this model to beams of arbitrary cross-section by employing the boundary-element method. This article presents the development of a new finite torsion beam element, based on the gen-eralized strain beam formulation, that allows for the buckling and post-buckling behaviour of thin-walled beams in multibody systems. To this end a second order stiffness formulation is proposed which will be derived in two-steps:
1. The first step concerns the derivation of the elastic and geometric stiffness matrices using the above nonlinear theory of non-uniform torsion of open section beams [7,37,44]. Both axial shortening and nonlinear Wagner’s stiffening torques are taken into account. Since the twist rate is deformation dependent, not element size dependent [24], the angle of twist will be assumed to be moderately large. The stiffness matrices are derived using cubic (Hermitian) polynomials. The coupling among bending and torsion due to a non-coinciding centroid and shear centre is incorporated using a transformation matrix which accounts for the eccentric location of the shear centre.
2. In the second step, second order approximations for the axial elongation and bending curvatures are included in the stiffness formulation. For that, second order displacement and rotation fields are derived using a local cross-section rotation matrix of second order. The displacement and rotation fields are interpolated in the same way as in the linear case by means of linear shape functions for the axial elongation and cubic functions for the transverse displacements and twist rotations. Integrating the interpolated second order curvature and strain–displacement equations over the length of the beam element using the moment-area method [19], a set of modified deformations is obtained in which second order approximations for the axial elongation and bending curvatures are represented by additional quadratic terms in the expressions of the basic deformations.
This formulation combines the advantages of the co-rotational formulation with the consis-tency of the inertial frame approach, viz derivation of the inertia forces in terms of absolute nodal velocities and accelerations. Deformations associated with a large stiffness can be eliminated by constraining them to be zero. The inertia properties of the beam element are described using both consistent and lumped mass formulations. The latter is used to model
rotary and warping inertia of the beam cross-section. The outline of the paper is as follows: The generalized strain beam formulation framework for a thin-walled beam element is pre-sented in Sect.2, while the local beam kinematic formulation is introduced in Sect.3. The elastic and geometric stiffness matrices are presented in Sect.4, and in Sect.5a set of mod-ified deformations is derived. A lumped mass formulation describing twist rotary and warp-ing inertia of the beam’s cross-section is presented in Sect.6and the (linearized) equations of motion for the second order beam element are presented in Sect.7. Finally, numerical examples illustrating the performance of the present beam element are presented in Sect.8.
2 Generalized strain beam formulation
In the original description of the generalized strain beam formulation which was developed for stability analysis of structures [10], Euler angles were used to parameterize global nodal rotations. Van der Werff and Jonker [45] introduced a description including Euler parame-ters which is more appropriate for computations in multibody system codes and enabled an implementation in the SPACAR program [26].
2.1 Description of thin-walled beam model
The kinematic model is based on the following assumptions:
1. The beam is prismatic, initially straight and slender, i.e. the dimensions of the cross-sections are small with respect to the beam length.
2. Thin-walled open-section beams posses small torsional stiffness in comparison to their flexural rigidities. Consequently, the angle of twist will be assumed to be moderately large. However, elastic rotations due to flexure will be assumed to be small but finite. 3. The cross-section is not deformed in its own plane but is subjected to a warping
displace-ment perpendicular to the displaced cross-sectional plane. The warping distribution is assumed to be given by the product of the twist rate and the corresponding Saint Venant warping function of the cross-section.
4. The cross-section remains orthogonal to the centroid axis in the deformed configuration when out of plane warping is excluded (Euler–Bernoulli hypothesis).
Figure1shows a thin-walled beam element which is represented by an elastic line which coincides with the centroid of the cross-sections of the beam. The configuration of the el-ement is described by the position vectors rp= (Xp, Yp, Zp)Tand rq= (Xq, Yq, Zq)Tof
the centroid at the end (section) nodes p and q, respectively, in the global inertial coordinate system (OXYZ) with base vectors eX, eY, eZ. The orientations are described by orthogonal
triads of unit vectors (nx, ny, nz)which are rigidly attached to the end nodes p and q. The
vector nx is perpendicular to the cross-sectional plane of the beam and ny and nz are in
the directions of the principal axes of the cross-section. The rotation matrix from the global directions to the initially straight reference state is defined as
nx= R0eX, ny= R0eY, nz= R0eZ. (1)
The rotational part of the motion of the flexible beam is determined by the rotation of the triads, which are described by rotation matrices Rpand Rq, so
nxp= Rpnx, nyp= Rpny, nzp= Rpnz, nxq= Rqnx, nyq= Rqny, nzq= Rqnz. (2)
Fig. 1 spatial beam element,
reference and deformed state
So the complete orientations of the triads are given by the composite rotation matrices RpR
0
and RqR
0. If the beam is rigid then the rotation matrices Rpand Rqare identical and in the
initial undeformed state they are equal to the identity matrix. In the present description Euler parameters are used to parameterize the rotation matrices, but the formulation can easily be transformed if a different choice is made. If the Euler parameters are denoted by (λ0, λ)with
the scalar part λ0and the vector part λ= (λ1, λ2, λ3)T, a rotation matrix can be expressed
as [30]
R(λ0, λ)= [λ20− λ
Tλ]I + 2λ
0˜λ + 2λλT, (3)
where I is a 3× 3 unitary matrix. Furthermore use has been made of the tilde notation to denote the skew symmetric matrix associated with a vector,
˜λ = ⎡ ⎣ λ03 −λ03 −λλ21 −λ2 λ1 0 ⎤ ⎦ . (4)
By definition, the Euler parameters must satisfy the constraint equation
λ20+ λTλ= 1. (5)
This constraint condition is incorporated in the theory, using the so-called λ-element [25].
2.2 Deformation modes
The nodal coordinates of the spatial beam element are the six Cartesian coordinates rep-resenting the position vectors rp and rq, two sets of four Euler parameters (λp
0, λp)and (λq0, λq)and two warping coordinates αpand αq. A further explanation of the warping
coor-dinates follows in Sect.4. If a redundant parametrization for the rotations is used, only three of them are independent. Therefore, as the beam has six degrees of freedom as a rigid-body, eight independent deformation modes can be defined that are invariant under rigid-body motions of the element. A set of deformations εi, corresponding to each of the deformation
modes is then expressed as analytical functions of the nodal coordinates Xk εi= Di(Xk); i = 1, . . . , 8; k = 1, . . . , 16, or ε = D(X), (6) where X=rpT, (λ0, λpT), αp, rqT, (λ0, λqT), αq T . (7)
A suitable choice for the deformation functions is:
ε1= l − l0 (axial elongation), (8a)
ε2= l0arcsin (npT z nqy− np T y nqz)/2 (torsion), (8b) ε3= −l0nTln p z, ε5= l0nTln p y , ε4= l0nTln q z ε6= −l0nTln q y (bending), (8c) ε7= l20αp, ε8= l02αq (warping), (8d)
where l0is the reference length of the element, l= rq−rpis the vector from node p to node q, l= l is the distance between the nodes and nl= l/l is the unit vector directed from
node p to node q. The deformations are invariant under rigid-body displacements of the element, so constraining them to be zero enforces rigidity on the element. Furthermore, these deformations have a clear physical meaning which facilitates the description of strength and stiffness properties of the element. The first deformation ε1describes the axial elongation, ε2/ l0the twist angle, ε3− ε6represent bending deformations in the (xz)- and (xy)-planes,
respectively, and ε7/ l20, ε8/ l02describe the change of twist at the nodes p and q, respectively.
To be able to describe finite torsion deformations, the twist angle is obtained from the arcsin function. With this function twist angles−π/2 ≤ ε2/ l0≤ π/2 radians can be described. The
bending deformations are defined in terms of inner products of unit vector nland the unit
vectors ny, nzattached at the element nodes and visualized in Fig.2. Note that ε3− ε6does
not change if the beam undergoes an axial elongation with fixed orientations of the nodes and fixed direction nl. The physical dimension of all deformations is length.
2.3 Discrete stress resultants and equilibrium equations
Let us consider an equilibrium force system defined by the forces Fp, Fq, the moments Tp, Tq and the bi-moments Bp, Bq applied at the nodal points p and q of the free beam
element, which are placed in a vector of element nodal forces
F=FpT, TpT, Bp, FqT, TqT, BqT. (9) The bimoments have the physical dimension of moment multiplied by length [46]. Next we consider virtual variations of the nodal positions δrp, δrq, rotations δϕp, δϕqand warping
coordinates δαp, δαqwhich are collected in a vector of virtual nodal displacements
δu=δrpT, δϕpT, δαp, δrqT, δϕqT, δαqT. (10) The variations δϕpand δϕqdefine infinitesimally small rotations from a general
Fig. 2 Visualization of bending modes: ε3, ε4and ε5, ε6
the variations of the Euler parameters by a 3× 4 transformation matrix Λ(λ0, λ), as
δϕ= 2Λ δλ0 δλ = 2[−λ, λ0I+ ˜λ] δλ0 δλ . (11)
Because δRp= δ ˜ϕpRpand δRq= δ ˜ϕqRq, the variations of the unit vectors in Eq. (2) can
be expressed as: δnpx= δϕp× npx, δnpy= δϕp× npy, δnpz = δϕp× npz, δnqx= δϕq× nqx, δnqy= δϕq× nqy, δnqz= δϕq× nqz. (12)
The virtual variations of the deformations δε are related to the virtual nodal displacements
δuby the relationship
δε= D,XδX= D,X ∂(δX) ∂(δu) δu= D,uδu= [D,up, D,uq] δup δuq , (13)
where the components of matrix D,uare derived using relations (12). For the deformations
(δrpT) (δϕpT) (δαp) Dp ,u= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ −nT l 0 T 0 0T l0(n p z × nqy− npy× nqz)T (npyTnqy+ npzTnqz) 0 l0 l(n p z− (n T lnpz)nl)T −l0(npz× nl)T 0 −l0 l(n q z− (n T ln q z)nl)T 0T 0 −l0 l(n p y− (n T lnpy)nl)T l0(npy× nl)T 0 l0 l(n q y− (n T ln q y)nl)T 0T 0 0 0 l2 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (14a) and (δrqT) (δϕqT) (δαq) D,uq = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ nTl 0T 0 0T l0(n q y× npz − nqz× npy)T (npyTnqy+ npzTnqz) 0 −l0 l(n p z − (nTlnpz)nl)T 0T 0 l0 l(n q z− (n T ln q z)nl)T l0 l(n q z × nl)T 0 l0 l(n p y− (n T ln p y)nl)T 0T 0 −l0 l(n q y− (n T ln q y)nl)T − l0 l(n q y× nl)T 0 0 0 0 0 0 l2 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (14b)
The discrete stress resultants σiare defined to be energetically dual to the deformations εi,
that is, the virtual work supplied by the discrete stress resultants is− σTδε. According to
the principle of virtual work, the element will be in a state of equilibrium if the equality
FTδu− σTδε= 0 (15)
holds for all δu and δε depending on δu by the relationship (13). Substituting Eq. (13) into Eq. (15) yields, with the transpose matrix DT
,u,
DT,uσ= F . (16)
These are the equilibrium equations formulated in the deformed configuration of the beam element. From these equations the equilibrium nodal force systems, expressed in terms of the
Fig. 3 (a) Equilibrium nodal force system F of nodal forces, moments and bi-moments defined by
Eq. (9) and (b)–(g) discrete stress resultants σ1− σ8corresponding to Eq. (16), where θi= arcsin(εi/ l0), i= 3, . . . , 6
discrete stress-resultants σi, are calculated and visualized in Fig.3(b)–(g). In all cases,
per-fect equilibrium is obtained for arbitrary large deformations and rigid-body displacements of the element. This is a direct consequence of the invariance of the deformations under rigid-body displacements, as the discrete stress resultants have no contribution to the virtual work in Eq. (15) for virtual rigid-body displacements, so the rigid-body displacements leaving the deformations unchanged can be described. In order to identify the discrete stress resultants of a beam element, we consider the undeformed configuration of the beam element, in which the unit vectors (npx, npy, npz)and (nqx, nqy, nqz)coincide with the global coordinate axes X, Y ,
Zas shown in Fig.3(a). Subsequently, the beam is loaded by an equilibrium force system F defined in Eq. (9). The discrete stress resultants σicorresponding to the nodal point forces
and moments can then be recognized as:
σ1= −F p x = Fxq (normal force), l0σ2= −Txp= Txq (twisting moment), l0σ3= −Typ, l0σ5= −Tzp, l0σ4= Tyq l0σ6= Tzq (bending moments), l02σ7= Bp, l02σ8= Bq (bimoments). (17)
Fig. 4 Configuration of the
elastic line in a co-rotational frame (x, y, z)
3 Local beam kinematic description
In order to provide second order approximations of the deformations and explicit expres-sions for the elastic and geometric stiffness matrices, the deformations will be redefined in a local co-rotational frame that continuously translates and rotates with the element. For beams with a non-symmetrical cross-section, the bending deformations are defined with re-spect to the shear center axis whereas the remaining ones are referred to the centroid of the cross-section.
3.1 Deformation functions
The deformations (ε2, . . . , ε8) defined in Eq. (8a)–(8d) can be redefined in a single
co-rotational frame (x, y, z) with the origin at node p and base vectors ex, ey and ez
point-ing along the respective positive directions of the coordinate axes as shown in Fig.4. The
x-axis coincides with the line connecting the nodal points p and q in the current configu-ration. The y-axis is chosen in such a way that the unit vector nyp lies in the (x, y)-plane
and the z-axis completes a right-handed orthogonal system of coordinate axes. In the un-deformed reference configuration, (ex, ey, ez)coincide with the unit vectors (nx, ny, nz)in
the global frame; ex is directed along the line of centroids (elastic line) and perpendicular
to the cross-section, and ey, ezare directed along the principal axes of the cross-section. An
arbitrary point c on the elastic line is indicated by the material coordinate x, 0≤ x ≤ l0,
and its coordinates in the deformed configuration arex+ uc(x), vc(x), wc(x)
, where uc,
vc and wc are the components of the displacement vector uc of point c. An orientation is
attached to it using a local cross-section rotation matrix R(x) which describes the rotation of a cross-section relative to the local frame due to flexural and torsional deformations of the beam. The matrix R(x) is defined as the product of three successive rotations about the rotated coordinate axes parameterized by modified Euler angles ϕx(x), ϕy(x), ϕz(x), also
referred as Tait–Bryan angles, as
R= R(ϕz)R(ϕy)R(ϕx) , (18) where R(ϕx)= ⎡ ⎣10 cos ϕ0 x − sin ϕ0 x 0 sin ϕx cos ϕx ⎤ ⎦ , R(ϕy)= ⎡ ⎣ cos ϕ0 y 01 sin ϕ0 y − sin ϕy 0 cos ϕy ⎤ ⎦ ,
R(ϕz)=
⎡
⎣cos ϕsin ϕzz − sin ϕcos ϕzz 00
0 0 1
⎤ ⎦ .
In contrast with Euler angles, modified Euler angles avoid singularity problems for small rotations. Another reason for taking modified Euler angles is that for small rotations they may be uniquely defined as components of a rotation pseudovector. In accordance with the definition of the co-rotational frame, the local nodal displacements and rotations are:
uc(0)= 0, vc(0)= 0, wc(0)= 0, uc(l0)= uqc, vc(l0)= 0, wc(l0)= 0, ϕx(0)= ϕxp= 0, ϕx,x(0)= αp, ϕy(0)= ϕ p y, ϕz(0)= ϕzp, ϕx(l0)= ϕxq, ϕx,x(l0)= αq, ϕy(l0)= ϕyq, ϕz(l0)= ϕzq. (19)
In these equations uqc is the axial displacement of node q, ϕxq, ϕyp, ϕqy, ϕzp, ϕzqare the local
nodal rotations about the centroid axis x and the principal y- and z-axes, and the warping coordinates αp, αq represent the derivatives of the twist angle ϕ
x at the nodes p and q,
respectively. The condition ϕxp= 0 implies that the y-axis is chosen in such a way that the
unit vector eyplies in the (x, y)-plane. The deformed configuration of the beam is now fully
determined by the displacement components uc, vc and wc of point c and the orientation
of the cross-sectional frameRex,Rey,Rez
. By replacing the unit vectors npy, npz and nqy,
nqz in Eq. (8a)–(8d) by the respective local unit vectors epy, epz and eqy, eqz we obtain the
deformation functions in the local co-rotational frame:
ε1= uqc (axial elongation), ε2= l0arcsin (ezpTeq y − e pT y e q z)/2 (torsion), ε3= −l0eTxe p z, ε5= l0eTxe p y, ε4= l0eTxe q z ε6= −l0eTxe q y (bending), ε7= l20αp, ε8= l02αq (warping), (20) where epy= R p ey, epz = R p ez, eqy= R q ey, eqz= R q ez. (21)
The rotation matrices Rp and Rq are obtained by evaluating the elastic rotation matrix R from Eq. (18) at the nodes p and q.
3.2 Second order approximation of deformations
Expanding matrix R from Eq. (18) around the identity matrix up to second order terms yields R= ⎡ ⎢ ⎣ 1−12ϕ2 y− 1 2ϕ 2 z −ϕz+ ϕxϕy ϕy+ ϕxϕz ϕz 1−12ϕx2− 1 2ϕ 2 z −ϕx+ ϕyϕz −ϕy ϕx 1−12ϕx2− 1 2ϕ 2 y ⎤ ⎥ ⎦ . (22)
Substituting Eq. (21) with the second order rotation matrix of Eq. (22) into Eq. (20), us-ing boundary conditions (19) and disregarding third and higher order terms, we obtain the
Fig. 5 Bending deformations
due to a twist rotation
ϕqx= ε2/ l0
following second order approximations of the torsion and bending deformations:
ε2= ε 2+ 1 2l0 (ε 3− ε 4)(ε5 + ε 6) (torsion), (23a) ε3= ε 3 ε4= ε 4+ 1 l0 ε 2ε 6 ε5= ε 5 ε6= ε 6− 1 l0 ε 2ε 4 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (bending), (23b) where ε
i(i= 2, . . . , 6) are the first order deformations defined by:
ε 2= l0ϕxq (torsion), (24a) ε 3= −l0ϕyp, ε 5= −l0ϕ p z , ε 4= l0ϕyq ε 6= l0ϕ q z (bending), (24b) ε7= l02αp, ε8= l02αq (warping). (24c)
The quadratic terms in the second order expressions for the torsion and bending deforma-tions originate from the additional nonlinearity due to the relative rotadeforma-tions of the unit vec-tors epy, epz and eqy, eqz at the nodal points p and q, see Fig.4.
3.3 Non-symmetrical cross-sections
For non-symmetric cross-sections with distinct shear centre and centroid, the shear centre is the pole of the twist rotation and hence, additional apparent bending deformations will occur due to a twist rotation ϕxq= ε2/ l0if the shear centre does not coincide with the centroid. This
will be explained on the basis of Fig.5. The figure shows the elastic line of the beam element in the co-rotational frame (x, y, z) with the origin at node p. The x-axis coincides with the centroid of the cross-sections and the y-axis and z-axis are principal axes; ysand zsare the
coordinates of the shear centre with respect to the centroid. We consider the end-sections at the nodes p and q normal to the undeflected shear centre axis. The cross-section at node q rotates through an angle ϕqx= ε2/l0about the shear axis (line connecting the shear points sp
shown in Fig.5. The associated bending deformations are
ε3= − ε4= yssin ϕxq− zs(1− cos ϕxq) ,
ε5= − ε6= ys(1− cos ϕxq)+ zssin ϕqx.
(25)
For small twist rotations sin ϕqx≈ ϕqx, cos ϕxq≈ 1, we obtain
ε3= − ε4= ˆysε2, ε5= − ε6= ˆzsε2, (26)
where ˆys= ys/ l0and ˆzs= zs/ l0. If the beam experiences both bending and torsion
defor-mations, we find that
ε3s= ε3− ˆysε2, ε4s= ε4+ ˆysε2, ε5s= ε5− ˆzsε2, ε6s= ε6+ ˆzsε2,
(27)
where εs
i(i= 3, . . . , 6) represent the bending deformations with respect to the coordinate
axes parallel to the principal axes through the shear centre. From Eq. (27) it can be observed that the transformation of bending deformations from the centroid coordinate system to the system of parallel axes passing through the shear centre depends on the twist angle
ϕxq= ε2/ l0. Note that the axial elongation ε1 is not changed by this transformation. The
transformation for the deformations can be expressed in the following matrix form:
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ε1 ε2 εs 3 εs 4 ε5s εs 6 ε7 ε8 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 − ˆys 1 0 0 0 0 0 0 ˆys 0 1 0 0 0 0 0 −ˆzs 0 0 1 0 0 0 0 ˆzs 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ε1 ε2 ε3 ε4 ε5 ε6 ε7 ε8 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (28a) or εs= Esε, (28b)
where Esis an (8×8) transformation matrix which accounts for the eccentricity of the shear
centre. This matrix will be used later in Sect.4.3to express the elastic and geometric stiffness matrices of the element in terms of the local centroid coordinates of the cross-section.
4 Stiffness properties
Flexible elements are modelled by allowing non-zero deformations and specifying consti-tutive equations relating deformations and discrete stress-resultants. In order to deal with geometric nonlinearities that arise from change in configuration a second order stiffness formulation is derived in a two step approach. The first step includes the derivation of the elastic and geometric stiffness matrices. The second step (Sect.5) concerns the derivation of the modified deformations.
Fig. 6 Local coordinate system
and displacements and rotations
4.1 Second order approximation of the Green–Lagrange strains
We consider a typical cross-section of a thin-walled beam in a Cartesian co-rotational frame
(x, y, z)with base vectors ex, ey and ezpointing along the respective positive directions of
the coordinate axes as shown in Fig.6. A plane normal to the x-axis cuts the middle surface in a line called contour of the cross-section. The origin coincides with centroid c and the x-axis coincides with the centroid x-axis of the beam before deformation. The y-x-axis and z-x-axis are the initial principal axes of the cross-section; ysand zs are the coordinates of the shear
centre s; uc, vs and ws are rigid-body translations of the cross-section in the x-direction
at the centroid and the y- and z-directions at the shear centre, respectively; ϕx(x)is the
angle of twist and ϕy(x), ϕz(x)are rigid-body rotations about the y-axis and z-axis,
respec-tively. According to the second assumption (Sect.2.1) the angle of twist will be assumed to be moderately large, but elastic rotations due to flexure will be assumed to be small but finite. From the third assumption of in-plane rigid cross-sections (εyy= εzz= εyz= 0), the
beam is deformed in a way such that the cross-section remains its shape (no distortion) and neglects the in-plane Poisson’s effect. Under these conditions the transverse displacement components v, w of an arbitrary point on the cross-section are given by [11]
v(x, y, z)= vs(x)− (z − zs)sin ϕx− (y − ys)(1− cos ϕx) ,
w(x, y, z)= ws(x)+ (y − ys)sin ϕx− (z − zs)(1− cos ϕx) ,
(29a)
and according to the fourth assumption the axial displacement u is given by [7,36]
u(x, y, z)= uc(x)− y αz− z αy− ω ϕx, x, (29b)
where ucis the average axial displacement of the cross-section, αyand αzare defined by
αy= ws, xcos ϕx− vs, xsin ϕx, αz= vs, xcos ϕx+ ws, xsin ϕx, (30)
and ω(y, z) is the warping function [20], defining the warping of the cross-section with re-spect to the shear centre. The angles αyand αzapproximate the slopes between the beam
axis in the deformed state and the initial horizontal and vertical planes, respectively. The warping function with its sectorial pole at the shear centre, satisfies the following condi-tions [46]: A ω dA= 0 , A ω y dA= 0 , A ω z dA= 0 , (31)
where A represents the cross-sectional area. Assuming that the axial deformation is small, we can neglect the second order terms involving u, x and obtain the following expressions
for the nonlinear Green–Lagrange strain tensor: εxx= u, x+12 v, x2 + w, x2 , γxy= u, y+ v, x(1+ v, y)+ w, xw, y, γxz= u, z+ w, x(1+ w, z)+ v, xv, z. (32)
Substituting Eqs. (29a), (29b) into Eq. (32), the second order approximation of the Green– Lagrange strains can be written as follows:
εxx= εx− y κz+ z κy− ω ϕx, xx+12R 2ϕ2 x, x, (33a) γxy= −(z − zs+ ω, y) ϕx, x, γxz= (y − ys− ω, z) ϕx, x, (33b) where εxis defined by εx= uc, x+12 vs, x2 + ws, x2 + (zsαz− ysαy) ϕx, x, (34)
κyand κzare the bending curvatures about the main axes y and z, defined by
κy= −ws, xxcos ϕx+ vs, xxsin ϕx,
κz= vs, xxcos ϕx+ ws, xxsin ϕx
(35)
and R is the radial distance from a point on the cross-section to the shear centre, given by
R2= (y − ys)2+ (z − zs)2. (36)
Since thin-walled open-section beams posses small torsional stiffness in comparison to their flexural rigidities, nonlinear torsional curvatures can be neglected. Equation(33a) is the fun-damental expression for the axial strain of thin-walled beams that can undergo large angles of twist. The fourth term in this equation represents the warping effect. The last term is pro-portional to ϕ2
x, x and is called shortening term. Since only moderately large angles of twist
will be considered, the trigonometric functions of the angle of twist will be approximated by a truncated third order Taylor series expansion as sin ϕx= ϕx− ϕx3/6 and cos ϕx= 1 − ϕx2/2.
Substituting Eqs. (29a) and (30) into Eq. (34) and retaining all terms up to third order yields the following relationship:
εx= uc, x+12
vc, x2 + wc, x2 −12(ys2+ z2s) ϕx, x2 , (37) and through substitution in Eq. (33a) we find the following expression for the averaged axial strain: εav= 1 A A εxxdA= uc, x+12 vc, x2 + wc, x2 +12Ip A ϕ 2 x, x, (38)
where Ipis the polar moment of area about the centroid. Similarly, a second order
approxi-mation of the bending curvatures can be derived from Eq. (35) as
κy= ϕy, x+ ϕxϕz, x, κz= ϕz, x− ϕxϕy, x, (39)
where
4.2 Local equilibrium equations and constitutive laws
The local equilibrium equations can be obtained by applying the principle of virtual work, for which the work done by the discrete stress resultants σson the virtual deformations δεs
is equal to σsTδεs= l0 0 A σxxδεxx+ τxyδγxy+ τxzδγxz dA dx . (41)
After integration over the cross-section A, we obtain
σsTδεs= l0 0 N δεx+Mxδϕx, x+Myδκy+Mzδκz+B δϕx, xx+MRϕx, xδϕx, x dx , (42)
where N is the axial force acting tangentially to the deformed longitudinal axis, My, Mzare
the resulting bending moments about the y-axis and z-axis in the deformed state, respec-tively, B is the bimoment acting on the cross-section in the deformed state, Mx is de Saint
Venant twisting moment and MRis a higher order stress resultant called Wagner’s moment.
They are defined by
N= A σxxdA, My= A σxxzdA, Mz= − A σxxydA, B= A σxxωdA, Mx= A τxz(y− ys− ω, z)− τxy(z− zs+ ω, y) dA, MR= A (σxx− σx0) R2dA , (43) where σx0 is such that the axial force resultant due to cross-sectional twist is equal to
zero [44]. For a linear elastic and isotropic material with elastic modulus E and shear mod-ulus G, the stress components σxx, τxy, τxzand σx0can be expressed as follows:
σxx= Eεxx, τxy= G γxy, τxz= G γxz, σx0=12E ϕ 2 x, x 1 A A R2dA. (44)
Substituting Eqs. (33a), (33b) and (44) into the expressions of Eq. (43) and using Eq. (31), we obtain the following relationships for the stress resultants in terms of deformation com-ponents in the principal axes:
N= A EεxxdA= EAεx+12EI0ϕx, x2 , Mx= A Gγxz(y− ys− ω, z)− γxy(z− zs+ ω, y) dA= G Itϕx, x, My= A EεxxzdA= EIy κy+ βzϕx, x2 , Mz= − A EεxxydA= EIz κz− βyϕ2x, x , (45) B= A EεxxωdA= EIω ϕx, xx− βωϕ2x, x ,
MR= A Eεxx−12 I0 Aϕ 2 x, x
R2dA= EI0εx+ 2EIyβzκy− 2EIzβyκz− 2EIωβωϕx, xx
+1 2EInϕ
2
x, x.
Here, Iy and Iz are second moments of inertia about the y- and z-axes, It is the de Saint
Venant torsion constant, Iωis the warping section constant [46], I0and IRare second and
fourth moments of area about the shear centre and Inis Wagner’s section constant. These
coefficients are defined as follows:
Iy= A z2dA , Iz= A y2dA , It= A (z− zs+ ω, y)2+ (y − ys− ω, z)2 dA , Iω= A ω2dA , I 0= A R2dA , I R= A R4dA , I n= IR− I02 A. (46) Note that IRis non-zero for bisymmetric sections and has an important influence in
nonlin-ear torsion behaviour and lateral post buckling response. The geometric coefficients βy, βz
and βω are called Wagner’s coefficients; βy, βz are associated with the bending curvature
under pure torsion and βω is associated with the modification to the torsional stiffness in
the presence of a bimoment [4]. These geometric parameters are expressed by the following relationships: βy= 1 2Iz A y (y2+ z2)dA− ys, βz= 1 2Iy A z (y2+ z2)dA− zs, βω= 1 2Iω A ω(y2+ z2)dA .
A numerical method for computing the coefficients βy, βz, βω and IR is proposed in [32,
36].
Using Eqs. (37) and (39), equations (45) can be written in matrix form as
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ N Mx My Mz B MR ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ EA 0 0 0 0 EI0 0 GIt 0 0 0 0 0 0 EIy 0 0 2EIyβz 0 0 0 EIz 0 −2EIzβy 0 0 0 0 EIω −2EIωβω
EI0 0 2EIyβz −2EIzβy −2EIωβω EIn ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ εx ϕx, x κy κz ϕx, xx 1 2ϕ 2 x, x ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = D γ = D (γl+ γc) , (47)
where γl= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ uc, x−12(y2s+ z 2 s) ϕ 2 x, x ϕx, x ϕy, x ϕz, x ϕx, xx 1 2ϕ 2 x, x ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , γc= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 2(v 2 c, x+ wc, x2 ) 0 ϕxϕz, x −ϕxϕy, x 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (48)
Matrix D is the material elasticity matrix. Its components are functions of the elastic and ge-ometric parameters. The second order approximations for the axial elongation and bending curvatures which are represented by the strain vector γcwill be included in step 2 using the concept of modified deformations described in Sect.5. The strain vector γlcan be split into a linear and nonlinear part expressed in terms of a gradient vector θ as
γl=H+12A θ , (49a) where H= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , A= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 −(y2 s + z2s) ϕx, x 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ϕx, x 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , θ= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ uc, x ϕx, x ϕy, x ϕz, x ϕx, xx ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (49b)
Matrix H is a so-called sign matrix whose entries are either 1 or 0. Matrix A contains the geometric nonlinear coupling terms between torsion deformation and axial elongation. From Eq. (49a), a variation of vector γlcan be written as
δγl=
H+ A δθ . (50)
With the relationships (49a) for γland (50) for δγl, the virtual work expression in Eq. (42) becomes
δεsTσs=
l0 0
δθTHT+ AT DH+12A θdx . (51) This equation is used to derive the elastic and geometric stiffness matrices of the element.
4.3 Stiffness matrices
The derivation of the stiffness matrices is based on assumed displacement and rotation fields expressed in terms of the element deformations. We employ a linear shape function for the
axial elongation, cubic functions for the twist angle ϕxand quadratic shape functions for the
flexural rotations ϕyand ϕz, respectively. At the shear centre the effects of bending and
tor-sion are decoupled. Therefore the flexural rotations are referred to the shear centre, whereas the remaining ones are referred to the centroid of the cross-section. The displacements and rotations of an arbitrary point on the centroid axis and shear centre axis with coordinate
x= ξl0can be then expressed as: uc= ξε1, l0ϕx= −2ξ3+ 3ξ2, ξ3− 2ξ2+ ξ, ξ3− ξ2 ⎡ ⎣εε27 ε8 ⎤ ⎦ , l0ϕy= −3ξ2+ 4ξ − 1, 3ξ2− 2ξ ε3s ε4s , l0ϕz= −3ξ2+ 4ξ − 1, 3ξ2− 2ξ εs5 εs6 , (52) where the deformations (εs
3, . . . , εs6) are defined in Eq. (28a). With the interpolations in
Eq. (52), the following relationship between the vectors θ and εscan be established:
θ= B εs, (53) where B= 1 l02 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ l0 0 0 0 0 0 0 0 0 −6ξ2+ 6ξ 0 0 0 0 3ξ2− 4ξ + 1 3ξ2− 2ξ 0 0 −6ξ + 4 6ξ − 2 0 0 0 0 0 0 0 0 −6ξ + 4 6ξ − 2 0 0 0 (−12ξ + 6)/l0 0 0 0 0 (6ξ− 4)/l0 (6ξ− 2)/l0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦.
By substituting Eq. (53) into (51), we obtain the stiffness equations
σs=Ss+ Gs εs, (54) where Ss= l 0 1 0 BTHTD H B dξ , (55a) Gs(ε 2, ε7, ε8)= l0 1 0 BTATD H+1 2(H TD A+ ATD A B dξ . (55b)
The first matrix Ss is the elastic stiffness matrix. The linearized stiffness equations are
ex-pressed by the relationship
σs= SsTεs, (56)
where Ss
Tis the tangent stiffness matrix, defined by
SsT= Ss+ GsT. (57)
The second matrix is the geometric stiffness matrix Gs
Tdefined by GsTε2, ε7, ε8)= l0 1 0 BTATD H+ HTD A+3 2A TD A B dξ . (58)
The components of the (8× 8) matrices Ssand Gs
Tcan be found in Appendix A. The vectors
σand σsand their variations are related by
σ= EsTσs, σ= EsTσs, (59)
where Es is the transformation matrix defined in Eqs. (28a), (28b). By substituting
Eqs. (28a), (28b) into Eqs. (54) and (56), we obtain with Eq. (59)
σ= (S + G) ε , σ= STε , (60)
where
S= EsTSsEs, G= EsTGsEs, ST= EsTSsTE
s. (61)
These are the element stiffness equations expressed in of the centroid coordinates. Note that in the initial undeformed configuration G(0)= 0 and ST(0)= 0. Therefore, the matrices G
and STare not decisive in a linear buckling analysis.
5 Modified deformations
In this section a set of modified deformations is derived in which second order approxima-tions for the axial elongation and bending curvatures represented by the strain vector γc are accounted for by additional quadratic terms in the expressions of the basic deforma-tions. The quadratic terms are derived through integration of the interpolated second order curvature and strain–displacement equations over the length of the beam element using the moment–area method. The displacement and rotation fields are interpolated in the same way as in the linear case by means of linear shape functions for the axial elongation and cubic functions for the transverse displacements and twist rotations. With the aid of inversion of Eqs. (23a), (23b), the original deformations εi(i= 2, . . . , 6) of Eqs. (8a)–(8d) can be
sub-stituted for the first order deformations ε
i(i= 2, . . . , 6) of Eqs. (24a)–(24c). Then Eq. (52)
can be rewritten in terms of second order deformations as:
uc = ξε1, l0ϕx = −2ξ3+ 3ξ2, ξ3− 2ξ2+ ξ, ξ3− ξ2 ⎡ ⎣ε2− 1 2l0(ε3− ε4)(ε5+ ε6) ε7 ε8 ⎤ ⎦ , l0ϕy = −3ξ2+ 4ξ − 1, 3ξ2− 2ξ εs 3 εs 4− 1 l0ε2ε s 6 , l0ϕz = −3ξ2+ 4ξ − 1, 3ξ2− 2ξ ε5s ε6s+l1 0ε2ε s 4 , vc = −ξ3+ 2ξ2− ξ, ξ3− ξ2 ε5 ε6+l1 0ε2ε4 , wc= ξ3− 2ξ2+ ξ, −ξ3+ ξ2 ε3 ε4−l1 0ε2ε6 . (62)
Integrating Eq. (38) over the length of the beam using the displacement and rotation fields of Eq. (62) gives ε1= 1 0 uc,ξdξ + 1 2l0 1 0 v2c,ξ+ w2c,ξ+Ip Aϕ 2 x, ξ dξ− ε1= ε1+ 1 30l0 2 ε32+ ε3ε4+ 2 ε24+ 2 ε52+ ε5ε6+ 2 ε26 + 1 30l0 Ip Al2 0 3 ε2 6ε2− ε7− ε8 + 2ε2 7− ε7ε8+ 2ε28 − ε1, (63)
where ε1 represents the axial elongation defined in Eq. (8a). The quadratic terms in the
bending deformations account for the second order axial shortening of the beam axis due to bending. The third term represents the axial shortening due to torsion and warping and ε1
represents the over-calculated second order elongation due to additional bending deforma-tions caused by a twist rotation of the cross-section if the shear centre does not coincide with the centroid. This is illustrated in Fig.5which shows that the elastic line becomes helical due to a twist rotation ε2/ l0of the cross-section about the shear centre axis. The additional
bending deformations are defined in Eq. (26) and their second order effects on the elongation can be calculated using the quadratic expressions in the bending deformations of Eq. (63). Substituting Eq. (26) then yields
ε1= 1 30l0 2 ε32+ ε3ε4+ 2 ε24+ 2 ε 2 5+ ε5ε6+ 2 ε26 = 1 10l0 ˆy2 s+ ˆz 2 s ε22. (64)
Second order approximations for the bending deformations of the shear centre axis can be obtained using the moment–area method [19]. This method is especially suitable when the deflection at only one point, of the beam is desired, because it is possible to find this deflec-tion without first finding the complete equadeflec-tion of the deflecdeflec-tion curve. Using the nonlinear relationships of the bending curvatures in Eq. (48) and the second order approximations of the rotation fields of Eq. (62), we obtain:
εs3= l0 1 0 ϕy,ξ+ ϕxϕz,ξ (1− ξ)dξ = εs 3+ ε2 10l0 εs5+ 2ε6s + 1 30l0 3εs 5ε7− ε8(εs5+ ε6s) , εs4= l0 1 0 ϕy,ξ+ ϕxϕz,ξ ξdξ = εs 4− ε2 10l0 2εs5+ ε6s − 1 30l0 3εs6ε8− ε7(εs5+ ε s 6) , εs5= l0 1 0 ϕz,ξ− ϕxϕy,ξ (1− ξ)dξ = εs 5− ε2 10l0 εs3+ 2ε4s − 1 30l0 3εs 3ε7− ε8(εs3+ ε4s) , εs6= l0 1 0 ϕz,ξ− ϕxϕy,ξ ξdξ = εs 6+ ε2 10l0 2εs3+ ε4s + 1 30l0 3εs4ε8− ε7(εs3+ ε s 4) . (65) The quadratic terms in the expressions of (65) take into account the effect of nonlinear bending caused by torsion and warping deformation. Note that transformation (28a), (28b) also applies to the modified deformations, i.e. εs= Esε. Substituting Eq. (27) into Eq. (65)
yields, with Eqs. (63) and (64), the following set of modified deformations: ε= E(ε), (66) where ε1= ε1+ 1 30l0 2 ε32+ ε3ε4+ 2 ε42+ 2 ε 2 5+ ε5ε6+ 2 ε26 + 1 30l0 Ip Al02 3 ε2(6ε2− ε7− ε8)+ 2ε27− ε7ε8+ 2ε28 − 1 10l0 ˆy2 s+ ˆz 2 s ε22 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (elongation), (67a) ε3= ε3+ ε2 10l0 ε5+ 2ε6+ ˆzs(ε2− ε7) + 1 30l0 3ε5ε7− ε8(ε5+ ε6) ε4= ε4− ε2 10l0 2ε5+ ε6− ˆzs(ε2− ε8) − 1 30l0 3ε6ε8− ε7(ε5+ ε6) ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ bending (x–z), (67b) ε5= ε5− ε2 10l0 ε3+ 2ε4+ ˆys(ε2− ε7) − 1 30l0 3ε3ε7− ε8(ε3+ ε4) ε6= ε6+ ε2 10l0 2ε3+ ε4− ˆys(ε2− ε8) + 1 30l0 3ε4ε8− ε7(ε3+ ε4) ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ bending (x–y), ε7= ε7 ε8= ε8 (warping). (67c)
Because the rigidity against axial elongation is much larger than the flexural and torsional rigidity, the modification of ε1 is the most relevant one. The other modifications are
rele-vant when large differences in flexural rigidities occur as can be the case for thin-walled open-section beams. Substituting Eq. (6) into Eqs. (67a)–(67c) yields a set of modified de-formations expressed as functions of the global nodal coordinates Xi, namely
ε= E(D(X)) = D(X) . (68)
The modified deformations form the basis of a second order beam finite element that cap-tures non-uniform torsion and flexural–torsional coupling of thin-walled beams with un-symmetrical cross-section. In the next two sections the inertia properties and equations of motion of the second order thin-walled beam element are derived.
6 Inertia properties
The inertia properties of the beam are described using both consistent and lumped mass formulations. A consistent mass formulation has been presented in [25,34]. In this section a lumped mass formulation is derived to describe the twist, rotary and warping inertia of the beam’s cross-section. To this end, we consider a cross-section of the beam in a local frame
Fig. 7 Position and orientation
of cross-section at node p
(x, y, z)attached at node p described by the vector rp and a set of four Euler parameters
(λp0, λp)that respectively represent the position and the orientation of this reference frame with respect to the inertial frame (O, X, Y, Z), see Fig.7. The x-axis is perpendicular to the cross-section and the y- and z-axes coincide with the principal axes of the cross-section. Let x= (0, y, z)Tdenote the position vector of an arbitrary point P of the cross-section, in the
initial undeformed configuration and let αpωbe the warping displacement perpendicular to
the cross-section at node p, where ω =ω(y, z),0, 0 T, with ω(y, z) being a normalized warping function and αp a warping coordinate describing the change of twist at node p.
The global position of point P in the deformed configuration can be described as
rP= rp+ RpR0
x+ αpω , (69)
where RpR
0 is a transformation matrix from the local coordinate system to the inertial
frame. The matrices R0and Rpare defined in Eqs. (1) and (3), respectively. Differentiating
Eq. (69) with respect to time yields under the small warping deformational-displacement approximation, the velocity vector
˙rP= ˙rp+ ˙RpR
0
x+ αpω + RpR0˙αpω , (70)
where (·)denotes differentiation with respect to time. This equation can be written in terms of the time derivatives of (λp0, λp)and αpas
˙rP= ˙rp+ Bp ˙λp 0 ˙λp + RpR 0˙αpω , (71) where Bp= −2 RpR 0[ ˜x + αp˜ω] RT0Λ p , (72)
in which ˜x and ˜ω are skew-symmetric matrices defined by [41]
˜x = ⎡ ⎣ 0z −z0 y0 −y 0 0 ⎤ ⎦ and ˜ω = ⎡ ⎣00 00 −ω0 0 ω 0 ⎤ ⎦ . (73) Matrix Λpis defined by Λp= [−λp, λp0I− ˜λp] . (74)
By applying one-half of the total twist, rotary and warping inertia’s to both end points of the beam, the kinetic energy of the lumped body at node p can be written as
Tp=14ρl0
A
˙rPT˙rPdA , (75)
where ρ and A are, respectively, the mass density and the area of the cross-section. Substi-tuting Eq. (71) into Eq. (75) yields
Tp=14ρl0 A ˙rpT,[˙λp 0, ˙λ pT ], ˙αp ⎡ ⎣I B p RpR 0ω BpTBp BpTRpR 0ω symm. ωTω ⎤ ⎦ ⎡ ⎢ ⎢ ⎣ ˙rp ˙λp0 ˙λp ˙αp ⎤ ⎥ ⎥ ⎦ dA. (76) If the origin of frame (x, y, z) coincides with the centre of mass of the body, the coupling term Bpdisappears. Since the matrices Λpand Rp, R
0are not space dependent, we obtain
using Eq. (31) A BpTBpdA= 4 ΛpTJpΛp, A RpR0ω dA= 0 , A ωTω dA= Iω. (77)
Substituting Eq. (77) into Eq. (76) yields the lumped mass matrix
M =12ρl0diag AI , 4 ΛpTJpΛp, Iω, AI , 4 Λ qT JqΛq, Iω , (78)
in which I is a 3× 3 unitary matrix, Iωis the sectorial moment defined in Eq. (46) and J p , Jqare defined by Jp= R0 I+ (αp)2Iω RT0, Jq= R0 I+ (αq)2Iω RT0, (79) where I= diag(Ip, Iy, Iz)and Iω = diag(0, Iω, Iω). Using Lagrange’s equations, the
cor-responding convective vector h can be expressed as
h =12ρl0 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 8 ˙ΛpTJpΛp+ αp˙αpΛpTI ωΛ p −2 αp˙λp 0, ˙λ pT ΛpTIωΛ p ˙λp 0 ˙λp 0 8 ˙ΛqTJqΛq+ αq˙αqΛqTI ωΛ q −2 αq˙λq 0, ˙λ qT ΛqTIωΛ q ˙λq 0 ˙λq ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (80)
7 Equations of motion
The equations of motion are derived using d’Alembert’s principle, which states that the balance of virtual work,
holds for all δε and δu depending on the variations of the nodal coordinates δX, by the relations δε= D, XδX , δu= A δX , (82) where D, X= E, DD, X, A= diag I, 2Λp,1, I, 2Λq,1 , (83) in which I is a 3× 3 unitary matrix and Λ is defined in Eq. (11). The last term in Eq. (81) represents the virtual work due to the inertia forces, where M(X) is a position-dependent mass matrix and h(X, ˙X) a convective inertia term being a quadratic function of nodal velocities. Both quantities are obtained by adding the consistent and lumped inertia parts, i.e. M= Mc+ M and h= hc+ h . Substituting the compatibility relations (82) into Eq. (81)
yields, with the transposed matrices DT, Xand AT,
DT, Xσ= ATF −M ¨X+ h . (84) These are the equations of motion for the free second order beam element. With the inclusion of the second order terms in the definitions of the modified deformations εi, relationships
(60) between deformations and discrete stress resultants remain valid, i.e.
σ= (S + G) ε , σ= STε . (85)
The stress-resultants σihave now a slightly modified meaning in case of finite deformations.
In order to study the vibration properties and the elastic stability, we consider small distur-bances from an equilibrium configuration
DT, Xσo= ATFo. (86) Substituting (85) into (84), expanding the equations in terms of disturbances X, ¨Xand disregarding second and higher order terms yields the linearized equations of motion
Mo ¨X+DoT, XSTD o , X+ D oT , XXσ o− (AoTFo) , X X= 0 , (87) or in matrix form Mo ¨X+Ko+ Go X= 0 , (88)
where Kois the structural stiffness matrix and Gothe structural geometric stiffness matrix
due to the reference load Fogiving rise to reference stress resultants σo. Note that (Ko+Go)
is symmetric.
8 Numerical examples
The proposed beam element has been implemented in the SPACAR software package [2,26] under the names BEAMW and BEAMWL element. The BEAMW element uses the modi-fied deformations of Eq. (66) whereas the BEAMWL element uses the basic deformations of Eqs. (8a)–(8d) in which the nonlinear bending curvatures are not taken into account. The SPACAR programme can make computations for multibody systems with rigid and flexible
