Three-dimensional beam element for pre- and post-buckling analysis of thin-walled beams in multibody systems

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 deformation modes is defined which are related to dual stress resultants in a co-rotational 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 nonlinear 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 elongation and bending curvatures are included by additional 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 behaviour 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.

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 shape optimization of multibody geometries. Due to their geometric characteristics, thin-walled beams with open cross-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 without 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 stiffening 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 attention in recent decades, little effort has been devoted to their implementation in flexible multibody systems codes. Finite element models of flexible beams in a multibody 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 capable 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 nonlinearity 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] presented 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]. However, 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 deformation 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 development 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 coordinate system. Flexible elements are modelled by allowing non-zero deformations and specifying constitutive equations relating deformations and dual stress resultants. The equilibrium equations are formulated in the deformed configuration of the beam element and are capable of accounting for arbitrary large deformations. Geometric nonlinearities arising from nonlinear couplings among axial elongation, bending 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 torsion 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], employed 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 generalized 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 consistency 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 presented 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. 5 a set of modified deformations is derived. A lumped mass formulation describing twist rotary and warping inertia of the beam’s cross-section is presented in Sect. 6 and 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 parameters 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 displacement 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).

Figure 1 shows 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 element is described by the position vectors $\boldsymbol{r}^{\,p}=(X^{p},Y^{p},Z^{p})^{\mathrm{T}}$ and $\boldsymbol{r}^{\,q}=(X^{q},Y^{q},Z^{q})^{\mathrm{T}}$ of the centroid at the end (section) nodes $p$ and $q$, respectively, in the global inertial coordinate system $(\mathit{OXYZ})$ with base vectors $\boldsymbol{e}_{X}$, $\boldsymbol{e}_{Y}$, $\boldsymbol{e}_{Z}$. The orientations are described by orthogonal triads of unit vectors $(\boldsymbol{n}_{x}, \boldsymbol{n}_{y}, \boldsymbol{n}_{z})$ which are rigidly attached to the end nodes $p$ and $q$. The vector $\boldsymbol{n}_{x}$ is perpendicular to the cross-sectional plane of the beam and $\boldsymbol{n}_{y}$ and $\boldsymbol{n}_{z}$ 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

$$ \boldsymbol{n}_{x}=\mathbf{R}_{0}\,\boldsymbol{e}_{X}\,, \quad \boldsymbol{n}_{y}=\mathbf{R}_{0} \,\boldsymbol{e}_{Y}\,, \quad \boldsymbol{n}_{z}=\mathbf{R}_{0}\,\boldsymbol{e}_{Z}\,. $$

(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 $\boldsymbol{R}^{p}$ and $\boldsymbol{R}^{q}$, so

$$ \textstyle\begin{array}{l} \boldsymbol{n}_{x}^{\,p} = \mathbf{R}^{p} \boldsymbol{n}_{x}, \\ \boldsymbol{n}_{y}^{\,p} = \mathbf{R}^{p} \boldsymbol{n}_{y}, \\ \boldsymbol{n}_{z}^{\,p} = \mathbf{R}^{p} \boldsymbol{n}_{z}, \end{array}\displaystyle \qquad \textstyle\begin{array}{l} \boldsymbol{n}_{x}^{\,q} = \mathbf{R}^{\,q} \boldsymbol{n}_{x}, \\ \boldsymbol{n}_{y}^{\,q} = \mathbf{R}^{\,q} \boldsymbol{n}_{y}, \\ \boldsymbol{n}_{z}^{\,q} = \mathbf{R}^{\,q} \boldsymbol{n}_{z}. \end{array} $$

(2)

So the complete orientations of the triads are given by the composite rotation matrices $\mathbf{R}^{p}\mathbf{R}_{0}$ and $\mathbf{R}^{q}\mathbf{R}_{0}$. If the beam is rigid then the rotation matrices $\mathbf{R}^{p}$ and $\mathbf{R}^{q}$ are 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 $(\lambda _{0},\boldsymbol{\lambda })$ with the scalar part $\lambda _{0}$ and the vector part $\boldsymbol{\lambda } = (\lambda _{1},\lambda _{2},\lambda _{3})^{\mathrm{T}}$, a rotation matrix can be expressed as [30]

$$ \mathbf{R}(\lambda _{0},\boldsymbol{\lambda }) = [\lambda _{0}^{2}- \boldsymbol{\lambda }^{\mathrm{T}}\boldsymbol{\lambda }]\mathbf{I} + 2\lambda _{0} \tilde{\boldsymbol{\lambda }} + 2\boldsymbol{\lambda }\boldsymbol{\lambda }^{\mathrm{T}}, $$

(3)

where $\mathbf{I}$ is a $3\times 3$ unitary matrix. Furthermore use has been made of the tilde notation to denote the skew symmetric matrix associated with a vector,

$$ \tilde{\boldsymbol{\lambda }} = \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \phantom{-}0 & -\lambda _{3} & \phantom{-} \lambda _{2} \\ \phantom{-}\lambda _{3} & \phantom{-}0 & -\lambda _{1} \\ -\lambda _{2} & \phantom{-}\lambda _{1} & \phantom{-}0 \end{array}\displaystyle \right ]. $$

(4)

By definition, the Euler parameters must satisfy the constraint equation

$$ \lambda _{0}^{2}+\boldsymbol{\lambda }^{\mathrm{T}}\boldsymbol{\lambda } = 1. $$

(5)

This constraint condition is incorporated in the theory, using the so-called $\lambda$-element [25].

Fig. 1

spatial beam element, reference and deformed state

2.2 Deformation modes

The nodal coordinates of the spatial beam element are the six Cartesian coordinates representing the position vectors $\boldsymbol{r}^{p}$ and $\boldsymbol{r}^{q}$, two sets of four Euler parameters $(\lambda _{0}^{p},\boldsymbol{\lambda }^{p})$ and $(\lambda _{0}^{q},\boldsymbol{\lambda }^{q})$ and two warping coordinates $\alpha ^{p}$ and $\alpha ^{q}$. A further explanation of the warping coordinates 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 $\varepsilon _{i}$, corresponding to each of the deformation modes is then expressed as analytical functions of the nodal coordinates $X_{k}$

$$ \textstyle\begin{array}{l} \varepsilon _{i} = D_{i}(X_{k}); \quad i=1,\ldots ,8; \; k=1,\ldots ,16, \quad \text{or}\quad \boldsymbol{\varepsilon }= \boldsymbol{D}(\boldsymbol{X}), \end{array} $$

(6)

where

$$ \boldsymbol{X} = \big[ \boldsymbol{r}^{p\mathrm{T}}, (\lambda _{0},\boldsymbol{\lambda }^{p \mathrm{T}}),\,\alpha ^{p}, \boldsymbol{r}^{q\mathrm{T}}, (\lambda _{0}, \boldsymbol{\lambda }^{q\mathrm{T}}),\,\alpha ^{q} \big]^{\mathrm{T}}\,. $$

(7)

A suitable choice for the deformation functions is:

$$\begin{aligned} \displaystyle \,\varepsilon _{1} = \, l - l_{0}\, \hspace{13.8em} &\mathrm{(axial~elongation),} \end{aligned}$$

(8a)

$$\begin{aligned} \displaystyle \,\varepsilon _{2} =l_{0}\,\mathrm{arcsin}\big((\boldsymbol{n}_{z}^{ \,p\mathrm{T}} \boldsymbol{n}_{y}^{q} - \boldsymbol{n}_{y}^{\,p\mathrm{T}} \boldsymbol{n}_{z}^{q} )/2\big)\, \hspace{4.0em} &\mathrm{(torsion),} \end{aligned}$$

(8b)

$$\begin{aligned} \textstyle\begin{array}{l} \varepsilon _{3} = -l_{0}\,\boldsymbol{n}_{l}^{\mathrm{T}}\boldsymbol{n}_{z}^{\,p}, \\ \varepsilon _{5} = \phantom{-} l_{0}\,\boldsymbol{n}_{l}^{\mathrm{T}}\boldsymbol{n}_{y}^{\,p}\,, \end{array}\displaystyle \qquad \textstyle\begin{array}{l} \varepsilon _{4} = \phantom{-} l_{0}\,\boldsymbol{n}_{l}^{\mathrm{T}}\boldsymbol{n}_{z}^{q} \hspace{3.8em} \\ \varepsilon _{6} = -l_{0}\,\boldsymbol{n}_{l}^{\mathrm{T}}\boldsymbol{n}_{y}^{q} \end{array}\displaystyle & \Biggl\} \mathrm{(bending),} \end{aligned}$$

(8c)

$$\begin{aligned} \textstyle\begin{array}{l} \varepsilon _{7} = \phantom{-} l_{0}^{2}\,{\alpha }^{p}, \hspace{3.3em} \varepsilon _{8} = \phantom{-} l_{0}^{2}\,{\alpha }^{q}\, \end{array}\displaystyle \hspace{4.3em} &\displaystyle \mathrm{(warping),} \end{aligned}$$

(8d)

where $l_{0}$ is the reference length of the element, $\boldsymbol{l} = \boldsymbol{r}^{\,q} - \boldsymbol{r}^{\,p}$ is the vector from node $p$ to node $q$, $l = \| \boldsymbol{l} \|$ is the distance between the nodes and $\boldsymbol{n}_{l}=\boldsymbol{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 $\varepsilon _{1}$ describes the axial elongation, $\varepsilon _{2}/l_{0}$ the twist angle, $\varepsilon _{3}- \varepsilon _{6}$ represent bending deformations in the ($xz$)- and ($xy$)-planes, respectively, and $\varepsilon _{7}/l_{0}^{2}$, $\varepsilon _{8}/l_{0}^{2}$ describe 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 $-\pi /2\leq \varepsilon _{2}/l_{0}\leq \pi /2$ radians can be described. The bending deformations are defined in terms of inner products of unit vector $\boldsymbol{n}_{l}$ and the unit vectors $\boldsymbol{n}_{y}$, $\boldsymbol{n}_{z}$ attached at the element nodes and visualized in Fig. 2. Note that $\varepsilon _{3}- \varepsilon _{6}$ does not change if the beam undergoes an axial elongation with fixed orientations of the nodes and fixed direction $\boldsymbol{n}_{l}$. The physical dimension of all deformations is length.

Fig. 2

Visualization of bending modes: $\varepsilon _{3}$, $\varepsilon _{4}$ and $\varepsilon _{5}$, $\varepsilon _{6}$

2.3 Discrete stress resultants and equilibrium equations

Let us consider an equilibrium force system defined by the forces $\boldsymbol{F}^{p}$, $\boldsymbol{F}^{q}$, the moments $\boldsymbol{T}^{p}$, $\boldsymbol{T}^{q}$ and the bi-moments $B^{p}$, $B^{q}$ applied at the nodal points $p$ and $q$ of the free beam element, which are placed in a vector of element nodal forces

$$ \boldsymbol{F} = \big[ \boldsymbol{F}^{p\mathrm{T}}, \boldsymbol{T}^{p\mathrm{T}},B^{\,p}, \boldsymbol{F}^{q\mathrm{T}}, \boldsymbol{T}^{q\mathrm{T}},B^{\,q} \big]^{\mathrm{T}}\,. $$

(9)

The bimoments have the physical dimension of moment multiplied by length [46]. Next we consider virtual variations of the nodal positions $\delta \boldsymbol{r}^{p}$, $\delta \boldsymbol{r}^{q}$, rotations $\delta \boldsymbol{\varphi }^{p}$, $\delta \boldsymbol{\varphi }^{q}$ and warping coordinates $\delta {\alpha }^{p}$, $\delta {\alpha }^{q}$ which are collected in a vector of virtual nodal displacements

$$ \delta \boldsymbol{u} = \big[ \delta \boldsymbol{r}^{\,p\mathrm{T}},\delta \boldsymbol{\varphi }^{\,p\mathrm{T}},\delta {\alpha }^{p}, \delta \boldsymbol{r}^{\,q \mathrm{T}},\delta \boldsymbol{\varphi }^{\,q\mathrm{T}},\delta {\alpha }^{q} \big]^{\mathrm{T}}\,. $$

(10)

The variations $\delta \boldsymbol{\varphi }^{p}$ and $\delta \boldsymbol{\varphi }^{q}$ define infinitesimally small rotations from a general configuration with components along the axes of the inertial coordinate system. These are related to the variations of the Euler parameters by a $3 \times 4$ transformation matrix $\boldsymbol{\varLambda }(\lambda _{0}, \boldsymbol{\lambda })$, as

$$ \delta \boldsymbol{\varphi } = 2\boldsymbol{\varLambda } \left [ \textstyle\begin{array}{c} \delta \lambda _{0} \\ \delta \boldsymbol{\lambda } \end{array}\displaystyle \right ] = 2[ -\boldsymbol{\lambda }, \;\; \lambda _{0} \,\boldsymbol{\mathbf{I}} + \tilde{\boldsymbol{\lambda }} ] \left [ \textstyle\begin{array}{c} \delta \lambda _{0} \\ \delta \boldsymbol{\lambda } \end{array}\displaystyle \right ] . $$

(11)

Because $\delta \mathbf{R}^{p} = \delta \tilde{\boldsymbol{\varphi }}^{p} \mathbf{R}^{p}$ and $\delta \mathbf{R}^{q} = \delta \tilde{\boldsymbol{\varphi }}^{q} \mathbf{R}^{q}$, the variations of the unit vectors in Eq. (2) can be expressed as:

$$ \textstyle\begin{array}{l} \delta \boldsymbol{n}_{x}^{p} = \delta \boldsymbol{\varphi }^{p} \times \boldsymbol{n}_{x}^{p} \, , \\ \delta \boldsymbol{n}_{y}^{p} = \delta \boldsymbol{\varphi }^{p} \times \boldsymbol{n}_{y}^{p} \, , \\ \delta \boldsymbol{n}_{z}^{p} = \delta \boldsymbol{\varphi }^{p} \times \boldsymbol{n}_{z}^{p} \, , \end{array}\displaystyle \qquad \textstyle\begin{array}{l} \delta \boldsymbol{n}_{x}^{q} = \delta \boldsymbol{\varphi }^{q} \times \boldsymbol{n}_{x}^{q} \, , \\ \delta \boldsymbol{n}_{y}^{q} = \delta \boldsymbol{\varphi }^{q} \times \boldsymbol{n}_{y}^{q} \, , \\ \delta \boldsymbol{n}_{z}^{q} = \delta \boldsymbol{\varphi }^{q} \times \boldsymbol{n}_{z}^{q} \, . \end{array} $$

(12)

The virtual variations of the deformations $\delta \boldsymbol{\varepsilon }$ are related to the virtual nodal displacements $\delta \boldsymbol{u}$ by the relationship

$$ \delta \boldsymbol{\varepsilon }=\boldsymbol{D}_{\boldsymbol{,X}}\,\delta \boldsymbol{X}=\boldsymbol{D}_{ \boldsymbol{,X}}\left [\! \frac{\partial (\delta \boldsymbol{X})}{\partial (\delta \boldsymbol{u})}\!\right ]\delta \boldsymbol{u} =\boldsymbol{D}_{\boldsymbol{,u}}\, \delta \boldsymbol{u} = [\boldsymbol{D}_{\boldsymbol{,u}}^{\,p}\,,\,\boldsymbol{D}_{\boldsymbol{,u}}^{\,q}] \left [ \textstyle\begin{array}{c} \delta \boldsymbol{u}^{p} \\ \delta \boldsymbol{u}^{q} \end{array}\displaystyle \right ] \,, $$

(13)

where the components of matrix $\boldsymbol{D}_{\boldsymbol{,u}}$ are derived using relations (12). For the deformations defined in Eq. (8a)–(8d) we obtain

and

The discrete stress resultants $\sigma _{i}$ are defined to be energetically dual to the deformations $\varepsilon _{i}$, that is, the virtual work supplied by the discrete stress resultants is $-\,\boldsymbol{\sigma }^{\mathrm{T}} \delta \boldsymbol{\varepsilon }$. According to the principle of virtual work, the element will be in a state of equilibrium if the equality

$$ \boldsymbol{F}^{\mathrm{T}} \delta \boldsymbol{u}\, - \,\boldsymbol{\sigma }^{\mathrm{T}} \delta \boldsymbol{\varepsilon } = 0 $$

(15)

holds for all $\delta \boldsymbol{u}$ and $\delta \boldsymbol{\varepsilon }$ depending on $\delta \boldsymbol{u}$ by the relationship (13). Substituting Eq. (13) into Eq. (15) yields, with the transpose matrix $\boldsymbol{D}_{\boldsymbol{,u}}^{\mathrm{T}}$,

$$ \boldsymbol{D}_{\boldsymbol{,u}}^{\mathrm{T}}\,\boldsymbol{\sigma }= \boldsymbol{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 discrete stress-resultants $\sigma _{i}$, are calculated and visualized in Fig. 3(b)–(g). In all cases, perfect 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 $(\boldsymbol{n}_{x}^{p}, \boldsymbol{n}_{y}^{p}, \boldsymbol{n}_{z}^{p})$ and $(\boldsymbol{n}_{x}^{q}, \boldsymbol{n}_{y}^{q}, \boldsymbol{n}_{z}^{q})$ coincide with the global coordinate axes $X$, $Y$, $Z$ as shown in Fig. 3(a). Subsequently, the beam is loaded by an equilibrium force system $\boldsymbol{F}$ defined in Eq. (9). The discrete stress resultants $\sigma _{i}$ corresponding to the nodal point forces and moments can then be recognized as:

$$ \textstyle\begin{array}{l@{\quad }l} \sigma _{1} = -F_{x}^{p} = \,F_{x}^{q} &\;\;\;\; \mathrm{(normal~force),} \\ l_{0}\sigma _{2} = -T_{x}^{p} = \,T_{x}^{q} &\;\;\;\; \mathrm{(twisting~moment),} \\ \textstyle\begin{array}{l} l_{0}\,\sigma _{3} = -T_{y}^{p}, \\ l_{0}\,\sigma _{5} = -T_{z}^{p}, \end{array}\displaystyle \qquad \textstyle\begin{array}{l} l_{0}\,\sigma _{4} = T_{y}^{q} \\ l_{0}\,\sigma _{6} = T_{z}^{q} \end{array}\displaystyle & \Biggl\} \mathrm{(bending~moments),} \\ \textstyle\begin{array}{l} l_{0}^{2}\,\sigma _{7} = \phantom{-} B^{p}, \hspace{2.2em} l_{0}^{2}\,\sigma _{8} = B^{q} \end{array}\displaystyle & \phantom{-}\;\mathrm{(bimoments).} \end{array} $$

(17)

Fig. 3

(a) Equilibrium nodal force system $\boldsymbol{F}$ of nodal forces, moments and bi-moments defined by Eq. (9) and (b)–(g) discrete stress resultants $\sigma _{1}-\sigma _{8}$ corresponding to Eq. (16), where $\theta _{i} = \arcsin ( \varepsilon _{i}/l_{0})$, $i=3,\dots ,6$

3 Local beam kinematic description

In order to provide second order approximations of the deformations and explicit expressions 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 respect to the shear center axis whereas the remaining ones are referred to the centroid of the cross-section.

3.1 Deformation functions

The deformations ($\varepsilon _{2},\dots ,\varepsilon _{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 $\boldsymbol{e}_{x}$, $\boldsymbol{e}_{y}$ and $\boldsymbol{e}_{z}$ pointing 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 configuration. The $y$-axis is chosen in such a way that the unit vector $\boldsymbol{n}_{y}^{\,p}$ lies in the $(x,y)$-plane and the $z$-axis completes a right-handed orthogonal system of coordinate axes. In the undeformed reference configuration, $(\boldsymbol{e}_{x},\boldsymbol{e}_{y},\boldsymbol{e}_{z})$ coincide with the unit vectors $(\boldsymbol{n}_{x},\boldsymbol{n}_{y},\boldsymbol{n}_{z})$ in the global frame; $\boldsymbol{e}_{x}$ is directed along the line of centroids (elastic line) and perpendicular to the cross-section, and $\boldsymbol{e}_{y}$, $\boldsymbol{e}_{z}$ are 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\leq x \leq l_{0}$, and its coordinates in the deformed configuration are $\big (x+u_{c}(x),\,v_{c}(x),\,w_{c}(x)\big )$, where $u_{c}$, $v_{c}$ and $w_{c}$ are the components of the displacement vector $\boldsymbol{u}_{c}$ of point $c$. An orientation is attached to it using a local cross-section rotation matrix $\overline{\mathbf{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 $\overline{\mathbf{R}}(x)$ is defined as the product of three successive rotations about the rotated coordinate axes parameterized by modified Euler angles $\varphi _{x}(x)$, $\varphi _{y}(x)$, $\varphi _{z}(x)$, also referred as Tait–Bryan angles, as

$$ \overline{\mathbf{R}} = \overline{\mathbf{R}}(\varphi _{z})\, \overline{\mathbf{R}}(\varphi _{y})\,\overline{\mathbf{R}}(\varphi _{x}) \,, $$

(18)

where

$$\begin{aligned} &\displaystyle \overline{\mathbf{R}}(\varphi _{x}) =\left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} 1 & 0 & 0 \\ 0 & \cos \varphi _{x} & -\sin \varphi _{x} \\ 0 & \sin \varphi _{x} & \phantom{-} \cos \varphi _{x} \end{array}\displaystyle \right ], \qquad \displaystyle \overline{\mathbf{R}}(\varphi _{y}) = \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \cos \varphi _{y} & 0 & \sin \varphi _{y} \\ 0 & 1 & 0 \\ -\sin \varphi _{y} & 0 & \cos \varphi _{y} \end{array}\displaystyle \right ], \\ &\displaystyle \overline{\mathbf{R}}(\varphi _{z}) = \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \cos \varphi _{z} & -\sin \varphi _{z} & 0 \\ \sin \varphi _{z} & \phantom{-}\cos \varphi _{z} & 0 \\ 0& 0 & 1 \end{array}\displaystyle \right ]. \end{aligned}$$

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:

$$ \textstyle\begin{array}{l@{\quad }l} & u_{c}(0) = 0, \\ & v_{c}(0) = 0, \\ & w_{c}(0) = 0, \\ \end{array}\displaystyle \qquad \textstyle\begin{array}{l@{\quad }l} &u_{c}(l_{0}) = u_{c}^{q}, \\ &v_{c}(l_{0}) = 0, \\ &w_{c}(l_{0}) = 0, \\ \end{array}\displaystyle \qquad \textstyle\begin{array}{l@{\quad }l} &\varphi _{x}(0) = \varphi _{x}^{p} = 0, \\ &\varphi _{x,x}(0)=\alpha ^{p}, \\ &\varphi _{y}(0) = \varphi _{y}^{p}, \\ &\varphi _{z}(0) = \varphi _{z}^{p}, \end{array}\displaystyle \qquad \textstyle\begin{array}{l@{\quad }l} &\varphi _{x}(l_{0}) = \varphi _{x}^{q}, \\ &\varphi _{x,x}(l_{0})=\alpha ^{q}, \\ &\varphi _{y}(l_{0}) = \varphi _{y}^{q}, \\ &\varphi _{z}(l_{0}) = \varphi _{z}^{q}. \end{array} $$

(19)

In these equations $u_{c}^{q}$ is the axial displacement of node $q$, $\varphi _{x}^{q}$, $\varphi _{y}^{p}$, $\varphi _{y}^{q}$, $\varphi _{z}^{p}$, $\varphi _{z}^{q}$ are the local nodal rotations about the centroid axis $x$ and the principal $y$- and $z$-axes, and the warping coordinates $\alpha ^{p}$, $\alpha ^{q}$ represent the derivatives of the twist angle $\varphi _{x}$ at the nodes $p$ and $q$, respectively. The condition $\varphi _{x}^{p} = 0$ implies that the $y$-axis is chosen in such a way that the unit vector $\boldsymbol{e}_{y}^{\,p}$ lies in the $(x,y)$-plane. The deformed configuration of the beam is now fully determined by the displacement components $u_{c}$, $v_{c}$ and $w_{c}$ of point $c$ and the orientation of the cross-sectional frame $\big (\overline{\mathbf{R}} \boldsymbol{e}_{x}\,, \overline{\mathbf{R}} \boldsymbol{e}_{y}\,,\overline{\mathbf{R}} \boldsymbol{e}_{z}\big )$. By replacing the unit vectors $\boldsymbol{n}_{y}^{p}$, $\boldsymbol{n}_{z}^{p}$ and $\boldsymbol{n}_{y}^{q}$, $\boldsymbol{n}_{z}^{q}$ in Eq. (8a)–(8d) by the respective local unit vectors $\boldsymbol{e}_{y}^{p}$, $\boldsymbol{e}_{z}^{p}$ and $\boldsymbol{e}_{y}^{q}$, $\boldsymbol{e}_{z}^{q}$ we obtain the deformation functions in the local co-rotational frame:

$$ \textstyle\begin{array}{l@{\quad }l} \displaystyle \,\varepsilon _{1} = u_{c}^{q}\, &\phantom{-}\mathrm{(axial~elongation),} \\ \displaystyle \,\varepsilon _{2} = l_{0}\,\mathrm{arcsin}\big((\boldsymbol{e}_{z}^{ \,p\mathrm{T}} \boldsymbol{e}_{y}^{q} \;-\,\; \boldsymbol{e}_{y}^{\,p\mathrm{T}} \boldsymbol{e}_{z}^{q})/2\big)\, &\phantom{-}\mathrm{(torsion),} \\ \textstyle\begin{array}{l} \varepsilon _{3} = -l_{0}\,\boldsymbol{e}_{x}^{\mathrm{T}}\boldsymbol{e}_{z}^{p}\,, \\ \varepsilon _{5} = \phantom{-} l_{0}\,\boldsymbol{e}_{x}^{\mathrm{T}}\boldsymbol{e}_{y}^{p}\,, \end{array}\displaystyle \qquad \textstyle\begin{array}{l} \varepsilon _{4} = \phantom{-} l_{0}\,\boldsymbol{e}_{x}^{\mathrm{T}}\boldsymbol{e}_{z}^{q}\, \\ \varepsilon _{6} = -l_{0}\,\boldsymbol{e}_{x}^{\mathrm{T}}\boldsymbol{e}_{y}^{q}\, \end{array}\displaystyle & \Biggl\} \,\mathrm{(bending),} \\ \textstyle\begin{array}{l} \varepsilon _{7} = \phantom{-} l_{0}^{2}\,{\alpha }^{p}, \hspace{3.0em} \varepsilon _{8} = \phantom{-} l_{0}^{2}\,{\alpha }^{q}\, \end{array}\displaystyle \hspace{5.2em} &\displaystyle \phantom{-} \mathrm{(warping),} \end{array} $$

(20)

where

$$ \textstyle\begin{array}{l} \boldsymbol{e}_{y}^{p} = \overline{\mathbf{R}}^{p} \boldsymbol{e}_{y}\,, \\ \boldsymbol{e}_{z}^{p} = \overline{\mathbf{R}}^{p} \boldsymbol{e}_{z}\,, \end{array}\displaystyle \qquad \textstyle\begin{array}{l} \boldsymbol{e}_{y}^{q} = \overline{\mathbf{R}}^{q} \boldsymbol{e}_{y}\,, \\ \boldsymbol{e}_{z}^{q} = \overline{\mathbf{R}}^{q} \boldsymbol{e}_{z}\,. \end{array} $$

(21)

The rotation matrices $\overline{\mathbf{R}}^{p}$ and $\overline{\mathbf{R}}^{q}$ are obtained by evaluating the elastic rotation matrix $\overline{\mathbf{R}}$ from Eq. (18) at the nodes $p$ and $q$.

Fig. 4

Configuration of the elastic line in a co-rotational frame $(x,y,z)$

3.2 Second order approximation of deformations

Expanding matrix $\overline{\mathbf{R}}$ from Eq. (18) around the identity matrix up to second order terms yields

$$ \displaystyle \overline{\mathbf{R}} = \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} 1-{\textstyle \frac{1}{2}}\varphi _{y}^{2} -{\textstyle \frac{1}{2}} \varphi _{z}^{2} & -\varphi _{z}+\varphi _{x}\varphi _{y} & \phantom{-}\varphi _{y}+\varphi _{x}\varphi _{z} \\ \phantom{-}\varphi _{z} & 1-{\textstyle \frac{1}{2}}\varphi _{x}^{2} -{ \textstyle \frac{1}{2}}\varphi _{z}^{2} & -\varphi _{x}+\varphi _{y} \varphi _{z} \\ -\varphi _{y} & \varphi _{x} & 1-{\textstyle \frac{1}{2}}\varphi _{x}^{2} -{\textstyle \frac{1}{2}}\varphi _{y}^{2} \\ \end{array}\displaystyle \right ]. $$

(22)

Substituting Eq. (21) with the second order rotation matrix of Eq. (22) into Eq. (20), using boundary conditions (19) and disregarding third and higher order terms, we obtain the following second order approximations of the torsion and bending deformations:

$$\begin{aligned} \displaystyle \varepsilon _{2} = \varepsilon _{2}^{\ell } + \frac{1}{2l_{0}}(\varepsilon _{3}^{\ell }-\varepsilon _{4}^{\ell })( \varepsilon _{5}^{\ell }+\varepsilon _{6}^{\ell })\, &\displaystyle \textrm{~~~(torsion),} \end{aligned}$$

(23a)

$$\begin{aligned} \left . \textstyle\begin{array}{l} \displaystyle \varepsilon _{3} = \varepsilon _{3}^{\ell }\, \\ \displaystyle \varepsilon _{4} = \varepsilon _{4}^{\ell }\, + \frac{1}{l_{0}}\varepsilon _{2}^{\ell }\,\varepsilon _{6}^{\ell }\, \\ \displaystyle \varepsilon _{5} = \varepsilon _{5}^{\ell }\, \\ \displaystyle \varepsilon _{6} = \varepsilon _{6}^{\ell }\, - \frac{1}{l_{0}}\varepsilon _{2}^{\ell }\,\varepsilon _{4}^{\ell }\, \end{array}\displaystyle \right \} \hspace{4.2em} &\displaystyle \textrm{~~(bending),} \end{aligned}$$

(23b)

where $\varepsilon ^{\ell }_{i}\, (i=2,\dots ,6)$ are the first order deformations defined by:

$$\begin{aligned} \varepsilon _{2}^{\ell } = \phantom{-} l_{0}\,\varphi _{x}^{q}\, \hspace{7.8em} &\displaystyle ~~~\textrm{(torsion),} \end{aligned}$$

(24a)

$$\begin{aligned} \textstyle\begin{array}{l} \varepsilon _{3}^{\ell } = -l_{0} \,\varphi _{y}^{p}\, ,~~~ \\ \varepsilon _{5}^{\ell } = -l_{0} \,\varphi _{z}^{p}\, ,~~~ \end{array}\displaystyle \quad \left . \textstyle\begin{array}{l} \varepsilon _{4}^{\ell } = l_{0} \,\varphi _{y}^{q}\, \\ \varepsilon _{6}^{\ell } = l_{0} \,\varphi _{z}^{q}\, \end{array}\displaystyle \right \} &\displaystyle ~~~\textrm{(bending),} \end{aligned}$$

(24b)

$$\begin{aligned} \textstyle\begin{array}{l} \varepsilon _{7} = \phantom{-} l_{0}^{2}\,{\alpha }^{p}, \hspace{2.3em} \varepsilon _{8} = l_{0}^{2}\,{\alpha }^{q}\, \phantom{-} \end{array}\displaystyle &\displaystyle ~~~\textrm{(warping).} \end{aligned}$$

(24c)

The quadratic terms in the second order expressions for the torsion and bending deformations originate from the additional nonlinearity due to the relative rotations of the unit vectors $\boldsymbol{e}_{y}^{p}$, $\boldsymbol{e}_{z}^{p}$ and $\boldsymbol{e}_{y}^{q}$, $\boldsymbol{e}_{z}^{q}$ 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 $\varphi _{x}^{q}=\varepsilon _{2}/l_{0}$ if 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; $y_{s}$ and $z_{s}$ are 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 $\varphi _{x}^{q}={\varepsilon _{2}}/{l_{0}}$ about the shear axis (line connecting the shear points $s^{p}$ and $s^{q}$). As a result of this twist rotation, the elastic line becomes helical after rotation as shown in Fig. 5. The associated bending deformations are

$$ \textstyle\begin{array}{l} \displaystyle \varepsilon _{3}=-\,\varepsilon _{4}=y_{s}\,\sin { \varphi _{x}^{q}}-{z}_{s}\,(1-\cos {\varphi _{x}^{q}}) \,, \\ \displaystyle \varepsilon _{5}=-\,\varepsilon _{6}=y_{s}\,(1-\cos { \varphi _{x}^{q}})+{z}_{s}\,\sin {\varphi _{x}^{q}}\,. \end{array} $$

(25)

For small twist rotations $\sin {\varphi _{x}^{q}}\approx \varphi _{x}^{q}$, $\cos {\varphi _{x}^{q}}\approx 1$, we obtain

$$ \displaystyle \varepsilon _{3}=-\,\varepsilon _{4}=\hat{y}_{s}\, \varepsilon _{2}\,, \qquad \displaystyle \varepsilon _{5}=-\, \varepsilon _{6}=\hat{z}_{s}\,\varepsilon _{2}\,, $$

(26)

where $\hat{y}_{s}=y_{s}/l_{0}$ and $\hat{z}_{s}=z_{s}/l_{0}$. If the beam experiences both bending and torsion deformations, we find that

$$ \textstyle\begin{array}{l} \displaystyle \varepsilon _{3}^{s}=\varepsilon _{3}-\hat{y}_{s}\, \varepsilon _{2}\,, \qquad \varepsilon _{4}^{s}=\varepsilon _{4}+ \hat{y}_{s}\,\varepsilon _{2}\,, \\ \displaystyle \varepsilon _{5}^{s}=\varepsilon _{5}-\hat{z}_{s}\, \varepsilon _{2}\,, \qquad \varepsilon _{6}^{s}=\varepsilon _{6}+ \hat{z}_{s}\,\varepsilon _{2}\,, \end{array} $$

(27)

where $\varepsilon _{i}^{s}\, (i=3,\dots ,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 $\varphi _{x}^{q}=\varepsilon _{2}/l_{0}$. Note that the axial elongation $\varepsilon _{1}$ is not changed by this transformation. The transformation for the deformations can be expressed in the following matrix form:

$$\begin{aligned} \left [ \textstyle\begin{array}{c} \varepsilon _{1} \\ \varepsilon _{2} \\ \varepsilon _{3}^{\,s} \\ \varepsilon _{4}^{\,s} \\ \varepsilon _{5}^{\,s} \\ \varepsilon _{6}^{\,s} \\ \varepsilon _{7} \\ \varepsilon _{8} \\ \end{array}\displaystyle \right ] = \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 &-\hat{y}_{s} & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \phantom{-}\hat{y}_{s} & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 &-\hat{z}_{s} & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & \phantom{-}\hat{z}_{s} & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c} \varepsilon _{1} \\ \varepsilon _{2} \\ \varepsilon _{3} \\ \varepsilon _{4} \\ \varepsilon _{5} \\ \varepsilon _{6} \\ \varepsilon _{7} \\ \varepsilon _{8} \\ \end{array}\displaystyle \right ] , \end{aligned}$$

(28a)

or

$$\begin{aligned} \boldsymbol{\varepsilon }^{s} = \mathbf{E^{s}}\,\boldsymbol{\varepsilon }, \end{aligned}$$

(28b)

where $\mathbf{E^{s}}$ is an $(8\times 8)$ transformation matrix which accounts for the eccentricity of the shear centre. This matrix will be used later in Sect. 4.3 to express the elastic and geometric stiffness matrices of the element in terms of the local centroid coordinates of the cross-section.

Fig. 5

Bending deformations due to a twist rotation $\varphi _{x}^{q}=\varepsilon _{2}/l_{0}$

4 Stiffness properties

Flexible elements are modelled by allowing non-zero deformations and specifying constitutive 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.

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 $\boldsymbol{e}_{x}$, $\boldsymbol{e}_{y}$ and $\boldsymbol{e}_{z}$ pointing 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 axis of the beam before deformation. The $y$-axis and $z$-axis are the initial principal axes of the cross-section; $y_{s}$ and $z_{s}$ are the coordinates of the shear centre $s$; $u_{c}$, $v_{s}$ and $w_{s}$ 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; $\varphi _{x}(x)$ is the angle of twist and $\varphi _{y}(x)$, $\varphi _{z}(x)$ are rigid-body rotations about the $y$-axis and $z$-axis, respectively. 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 ($\varepsilon _{yy}=\varepsilon _{zz}=\varepsilon _{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]

$$ \textstyle\begin{array}{l} v(x,y,z)=v_{s}(x)-(z-z_{s})\sin {\varphi _{x}}-(y-y_{s})(1-\cos { \varphi _{x}})\,, \\ w(x,y,z)=w_{s}(x)+(y-y_{s})\sin {\varphi _{x}}-(z-z_{s})(1-\cos { \varphi _{x}})\,, \end{array} $$

(29a)

and according to the fourth assumption the axial displacement $u$ is given by [7, 36]

$$ u(x,y,z)=u_{c}(x)-y\,\alpha _{z}-z\,\alpha _{y}-\omega \,\varphi _{x, \,x}\;, $$

(29b)

where $u_{c}$ is the average axial displacement of the cross-section, $\alpha _{y}$ and $\alpha _{z}$ are defined by

$$ \textstyle\begin{array}{l} \alpha _{y}=w_{s,\,x}\cos {\varphi _{x}}-v_{s,\,x}\sin {\varphi _{x}}\,, \qquad \alpha _{z}=v_{s,\,x}\cos {\varphi _{x}}+w_{s,\,x}\sin {\varphi _{x}} \,, \end{array} $$

(30)

and $\omega (y,z)$ is the warping function [20], defining the warping of the cross-section with respect to the shear centre. The angles $\alpha _{y}$ and $\alpha _{z}$ approximate 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 conditions [46]:

$$ \int _{A}\omega \,dA =0\,, \qquad \int _{A}\omega \, y \,dA =0\,, \qquad \int _{A}\omega \, 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:

$$ \textstyle\begin{array}{l@{\quad }l} &\displaystyle \varepsilon _{xx}= u_{,\,x} + {\textstyle \frac{1}{2}} \big(v_{,\,x}^{\,2} + w_{,\,x}^{\,2}\big)\,, \\ &\displaystyle \gamma _{xy} = u_{,\,y} + v_{,\,x}(1+v_{,\,y}) + w_{, \,x}w_{,\,y} \;, \\ &\displaystyle \gamma _{xz} = u_{,\,z} + w_{,\,x}(1+w_{,\,z}) + v_{, \,x}v_{,\,z} \;. \\ \end{array} $$

(32)

Substituting Eqs. (29a), (29b) into Eq. (32), the second order approximation of the Green–Lagrange strains can be written as follows:

$$\begin{aligned} &\varepsilon _{xx}= \varepsilon _{x}-y\,\kappa _{z}+z\,\kappa _{y}- \omega \,\varphi _{x,\,xx} +{\textstyle \frac{1}{2}}R^{2}\,\varphi _{x, \,x}^{2}\;, \end{aligned}$$

(33a)

$$\begin{aligned} & \textstyle\begin{array}{l} \gamma _{xy} = -(z-z_{s}+\omega _{,\,y})\,\varphi _{x,\,x}\;, \\ \gamma _{xz} = \phantom{-} (y-y_{s}-\omega _{,\,z})\,\varphi _{x,\,x}\;, \\ \end{array}\displaystyle \end{aligned}$$

(33b)

where $\varepsilon _{x}$ is defined by

$$ \displaystyle \varepsilon _{x} = u_{c,\,x} +{\textstyle \frac{1}{2}} \big(v_{s,\,x}^{\,2} + w_{s,\,x}^{\,2}\big) +(z_{s}\,\alpha _{z}-y_{s} \,\alpha _{y})\,\varphi _{x,\,x}\;, $$

(34)

$\kappa _{y}$ and $\kappa _{z}$ are the bending curvatures about the main axes y and z, defined by

$$ \displaystyle \textstyle\begin{array}{l} \kappa _{y}=-w_{s,\,xx}\cos {\varphi _{x}}+v_{s,\,xx}\sin {\varphi _{x}} \;, \\ \kappa _{z}=~ \phantom{-} v_{s,\,xx}\cos {\varphi _{x}}+w_{s,\,xx}\sin {\varphi _{x}} \end{array} $$

(35)

and $R$ is the radial distance from a point on the cross-section to the shear centre, given by

$$ R^{2}=(y-y_{s})^{2}+(z-z_{s})^{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 fundamental 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 proportional to $\varphi _{x,\,x}^{2}$ 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 {\varphi _{x}}=\varphi _{x}-\varphi _{x}^{3}/6$ and $\cos {\varphi _{x}}=1-\varphi _{x}^{2}/2$. Substituting Eqs. (29a) and (30) into Eq. (34) and retaining all terms up to third order yields the following relationship:

$$ \displaystyle \varepsilon _{x} = u_{c,\,x} +{\textstyle \frac{1}{2}} \big(v_{c,\,x}^{\,2} + w_{c,\,x}^{\,2}\big) -{\textstyle \frac{1}{2}}(y_{s}^{2}+z_{s}^{2}) \,\varphi _{x,\,x}^{2}\,, $$

(37)

and through substitution in Eq. (33a) we find the following expression for the averaged axial strain:

$$ \displaystyle \varepsilon _{av} =\frac{1}{A}\int _{A} \varepsilon _{xx}\,\mathrm{d}A= u_{c,\,x} +{\textstyle \frac{1}{2}} \big(v_{c,\,x}^{\,2} + w_{c,\,x}^{\,2}\big) \displaystyle +{ \textstyle \frac{1}{2}}\frac{I_{p}}{A}\,{\varphi }_{x,\,x}^{2}\;, $$

(38)

where $I_{p}$ is the polar moment of area about the centroid. Similarly, a second order approximation of the bending curvatures can be derived from Eq. (35) as

$$ \displaystyle \textstyle\begin{array}{l} \kappa _{y} = \varphi _{y,\,x} + \varphi _{x}\,\varphi _{z,\,x}\;, \qquad \kappa _{z} = \varphi _{z,\,x} - \varphi _{x}\,\varphi _{y,\,x} \;, \end{array} $$

(39)

where

$$ \varphi _{y,\,x}=-w_{s,\,xx}\;, \qquad \varphi _{z,\,x}=v_{s,\,xx}\;. $$

(40)

Fig. 6

Local coordinate system and displacements and rotations

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 $\boldsymbol{\sigma }^{s}$ on the virtual deformations $\delta \boldsymbol{\varepsilon }^{s}$ is equal to

$$ \boldsymbol{\sigma }^{s\mathrm{T}} \delta \boldsymbol{\varepsilon }^{s}=\int _{0}^{l_{0}} \int _{A}\big(\sigma _{xx}\,\delta \varepsilon _{xx} +\tau _{xy} \,\delta \gamma _{xy}+\tau _{xz}\,\delta \gamma _{xz}\big)\mathrm{d}A \,\mathrm{d}x\,. $$

(41)

After integration over the cross-section $A$, we obtain

$$ \boldsymbol{\sigma }^{s\mathrm{T}} \delta \boldsymbol{\varepsilon }^{s}=\int _{0}^{l_{0}} \Big(N\delta \varepsilon _{x} +M_{x}\,\delta \varphi _{x,\,x}+M_{y}\, \delta \kappa _{y}+M_{z}\,\delta \kappa _{z}+B\,\delta \varphi _{x,\,xx} + \textstyle M_{R}\,\varphi _{x,\,x}\delta \varphi _{x,\,x}\Big) \mathrm{d}x\,, $$

(42)

where $N$ is the axial force acting tangentially to the deformed longitudinal axis, $M_{y}$, $M_{z}$ are the resulting bending moments about the $y$-axis and $z$-axis in the deformed state, respectively, $B$ is the bimoment acting on the cross-section in the deformed state, $M_{x}$ is de Saint Venant twisting moment and $M_{R}$ is a higher order stress resultant called Wagner’s moment. They are defined by

$$ \textstyle\begin{array}{c} \displaystyle N=\int \limits _{A}\sigma _{xx}\,\mathrm{d}A, \qquad M_{y}= \int \limits _{A}\sigma _{xx}\;z\,\mathrm{d}A, \qquad M_{z}=-\int \limits _{A}\sigma _{xx}\,y\,\mathrm{d}A, \qquad B=\int \limits _{A} \sigma _{xx}\omega \,\mathrm{d}A, \\ \displaystyle M_{x}=\int \limits _{A}\big(\tau _{xz}(y-y_{s}-\omega _{, \,z})-\tau _{xy}(z-z_{s}+\omega _{,\,y})\big)\,\mathrm{d}A, \qquad M_{R}= \int \limits _{A}(\sigma _{xx}-\sigma _{x0})\,R^{2}\,\mathrm{d}A\,, \end{array} $$

(43)

where $\sigma _{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 modulus $G$, the stress components $\sigma _{xx}$, $\tau _{xy}$, $\tau _{xz}$ and $\sigma _{x0}$ can be expressed as follows:

$$ \sigma _{xx}=E\varepsilon _{xx}, \qquad \tau _{xy}=G\,\gamma _{\,xy}, \qquad \tau _{xz}=G\,\gamma _{\,xz}, \qquad \sigma _{x0}={\textstyle \frac{1}{2}}E\,\varphi _{x,\,x}^{2}\frac{1}{A}\!\int _{A}R^{2} \,\mathrm{d}A. $$

(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 components in the principal axes:

$$\begin{aligned} & \displaystyle N =\int \limits _{A} E\varepsilon _{xx}\,\mathrm{d}A = EA \varepsilon _{x} +\textstyle \frac{1}{2}EI_{0}\varphi _{x,\,x}^{2}\,, \\ &\displaystyle M_{x} =\int \limits _{A} G\big(\gamma _{xz}(y-y_{s}- \omega _{,\,z}) -\gamma _{xy}(z-z_{s}+\omega _{,\,y})\big)\, \mathrm{d}A = G\,I_{t}\varphi _{x,\,x}\,, \\ &\displaystyle M_{y} =\int \limits _{A} E\varepsilon _{xx}\,z\, \mathrm{d}A = EI_{y}\big(\kappa _{y}+\beta _{z}\varphi _{x,\,x}^{2} \big)\,, \\ &\displaystyle M_{z} =-\int \limits _{A} E\varepsilon _{xx}\,y\, \mathrm{d}A = EI_{z}\big(\kappa _{z} -\beta _{y}\varphi _{x,\,x}^{2} \big)\,, \\ &\displaystyle B = \int \limits _{A} E\varepsilon _{xx}\,\omega \, \mathrm{d}A = EI_{\omega }\big(\varphi _{x,\,xx} -\beta _{\omega } \varphi _{x,\,x}^{2}\big)\,, \\ &\displaystyle M_{R} =\int \limits _{A} E\big(\varepsilon _{xx} -{\textstyle \frac{1}{2}}\frac{I_{0}}{A}\varphi _{x,\,x}^{2} \big) \,R^{2}\,\mathrm{d}A = EI_{0}\varepsilon _{x}+2EI_{y}\beta _{z} \kappa _{y}-2EI_{z}\beta _{y}\kappa _{z} -2EI_{\omega }\beta _{\omega } \varphi _{x,\,xx} \\ &\hphantom{\displaystyle M_{R} =}{}+{\textstyle \frac{1}{2}}EI_{\mathrm{n}}\varphi _{x,\,x}^{2} \,. \end{aligned}$$

(45)

Here, $I_{y}$ and $I_{z}$ are second moments of inertia about the $y$- and $z$-axes, $I_{t}$ is the de Saint Venant torsion constant, $I_{\omega }$ is the warping section constant [46], $I_{0}$ and $I_{R}$ are second and fourth moments of area about the shear centre and $I_{n}$ is Wagner’s section constant. These coefficients are defined as follows:

$$ \textstyle\begin{array}{c} \displaystyle I_{y} = \int \limits _{A}z^{2} dA\,,\qquad I_{z} = \int \limits _{A}y^{2} dA\,,\qquad I_{t}=\int \limits _{A}\big((z-z_{s}+ \omega _{,\,y})^{2}+(y-y_{s}-\omega _{,\,z})^{2}\big)dA\,, \\ I_{\omega } = \int \limits _{A}{\omega }^{2} dA\,, \qquad I_{0}=\int \limits _{A} R^{\,2}dA\,,\qquad I_{\mathrm{R}}=\int \limits _{A} R^{\,4}dA \,, \qquad I_{\mathrm{n}}= I_{\mathrm{R}}-\textstyle \frac{I_{0}^{2}}{A} \,. \end{array} $$

(46)

Note that $I_{\mathrm{R}}$ is non-zero for bisymmetric sections and has an important influence in nonlinear torsion behaviour and lateral post buckling response. The geometric coefficients $\beta _{y}$, $\beta _{z}$ and $\beta _{\omega }$ are called Wagner’s coefficients; $\beta _{y}$, $\beta _{z}$ are associated with the bending curvature under pure torsion and $\beta _{\omega }$ 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:

$$ \begin{aligned} &\beta _{y} =\frac{1}{2I_{z}} \int _{A}y\,(y^{2}+z^{2})dA-y_{s} \,,\qquad \beta _{z} =\frac{1}{2I_{y}} \int _{A}z\,(y^{2}+z^{2})dA-z_{s} \,,\\ & \beta _{\omega }=\frac{1}{2I_{\omega }} \int _{A}\omega (y^{2}+z^{2})dA \,. \end{aligned} $$

A numerical method for computing the coefficients $\beta _{y}$, $\beta _{z}$, $\beta _{\omega }$ and $I_{\mathrm{R}}$ is proposed in [32, 36].

Using Eqs. (37) and (39), equations (45) can be written in matrix form as

$$\begin{aligned} \left [ \textstyle\begin{array}{l} N \\ M_{x} \\ M_{y} \\ M_{z} \\ B \\ M_{R} \end{array}\displaystyle \right ] &\!\! =\!\! \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} EA &~~~ 0 ~~& 0 & 0 & 0 & EI_{0} \\ 0 &~~~ GI_{t} ~~& 0 & 0 & 0 & 0 \\ 0 &~~~ 0 ~~& EI_{y} & 0 & 0 & \phantom{-}2EI_{y}\beta _{z} \\ 0 &~~~ 0 ~~& 0 & EI_{z} & 0 &-2EI_{z}\beta _{y} \\ 0 &~~~ 0 ~~& 0 & 0 & EI_{\omega } &-2EI_{\omega }\beta _{\omega } \\ EI_{0} &~~~ 0 ~~& 2EI_{y}\beta _{z} &-2EI_{z}\beta _{y} &-2EI_{\omega } \beta _{\omega } & EI_{\mathrm{n}} \\ \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{l} \varepsilon _{x} \\ \varphi _{x,\,x} \\ \kappa _{y} \\ \kappa _{z} \\ \varphi _{x,\,xx} \\ \textstyle \frac{1}{2}\varphi _{x,\,x}^{2} \end{array}\displaystyle \right ] \\ &=\,\mathbf{D}\,\boldsymbol{\gamma }\,=\,\mathbf{D}\,(\boldsymbol{\gamma }_{l}+ \boldsymbol{\gamma }_{c})\,, \end{aligned}$$

(47)

where

$$ \boldsymbol{\gamma }_{l}= \left [ \textstyle\begin{array}{l} u_{c,\,x}-{\textstyle \frac{1}{2}}(y_{s}^{2}+z_{s}^{2})\,\varphi _{x, \,x}^{2} \\ \varphi _{x,\,x} \\ \varphi _{y,\,x} \\ \varphi _{z,\,x} \\ \varphi _{x,\,xx} \\ \textstyle \frac{1}{2}\varphi _{x,\,x}^{2} \end{array}\displaystyle \right ]\,,\qquad \boldsymbol{\gamma }_{c}= \left [ \textstyle\begin{array}{c} {\textstyle \frac{1}{2}}(v_{c,\,x}^{\,2} + w_{c,\,x}^{\,2}) \\ 0 \\ \phantom{-}\varphi _{x}\,\varphi _{z,\,x} \\ -\varphi _{x}\,\varphi _{y,\,x} \\ 0 \\ 0 \end{array}\displaystyle \right ]\,. $$

(48)

Matrix $\mathbf{D}$ is the material elasticity matrix. Its components are functions of the elastic and geometric parameters. The second order approximations for the axial elongation and bending curvatures which are represented by the strain vector $\boldsymbol{\gamma }_{c}$ will be included in step 2 using the concept of modified deformations described in Sect. 5. The strain vector $\boldsymbol{\gamma }_{l}$ can be split into a linear and nonlinear part expressed in terms of a gradient vector $\boldsymbol{\theta }$ as

$$ \boldsymbol{\gamma }_{l}=\big(\mathbf{H}+\textstyle \frac{1}{2}\mathbf{A}\big) \boldsymbol{\theta }\,, $$

(49a)

where

$$ \begin{aligned} &\mathbf{H}=\left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} 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 \end{array}\displaystyle \right ]\,,\qquad \mathbf{A}=\left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} 0 &-(y_{s}^{2}+z_{s}^{2})\,\varphi _{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 & \varphi _{x,\,x} & 0 & 0 & 0 \end{array}\displaystyle \right ]\,, \\ & \boldsymbol{\theta }= \left [ \textstyle\begin{array}{l} u_{c,\,x} \\ \varphi _{x,\,x} \\ \varphi _{y,\,x} \\ \varphi _{z,\,x} \\ \varphi _{x,\,xx} \\ \end{array}\displaystyle \right ]\,. \end{aligned} $$

(49b)

Matrix $\mathbf{H}$ is a so-called sign matrix whose entries are either 1 or 0. Matrix $\mathbf{A}$ contains the geometric nonlinear coupling terms between torsion deformation and axial elongation. From Eq. (49a), a variation of vector $\boldsymbol{\gamma }_{l}$ can be written as

$$ \delta \boldsymbol{\gamma }_{l}=\big(\mathbf{H}+\mathbf{A}\big)\,\delta \boldsymbol{\theta }\,. $$

(50)

With the relationships (49a) for $\boldsymbol{\gamma }_{l}$ and (50) for $\delta \boldsymbol{\gamma }_{l}$, the virtual work expression in Eq. (42) becomes

$$ \delta \boldsymbol{\varepsilon }^{s\mathrm{T}}\boldsymbol{\sigma }^{s}=\int _{0}^{l_{0}} \delta \boldsymbol{\theta }^{\mathrm{T}} \big(\mathbf{H}^{\mathrm{T}}+\mathbf{A}^{\mathrm{T}}\big)\,\mathbf{D}\,\big(\mathbf{H}+{\textstyle \frac{1}{2}} \mathbf{A}\big)\boldsymbol{\theta }\mathrm{d}x\,. $$

(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 $\varphi _{x}$ and quadratic shape functions for the flexural rotations $\varphi _{y}$ and $\varphi _{z}$, respectively. At the shear centre the effects of bending and torsion 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=\xi l_{0}$ can be then expressed as:

$$ \textstyle\begin{array}{l@{\quad }l} \hspace{0.8em} u_{c}\hspace{-0.8em} &=\, \xi \varepsilon _{1}\,, \qquad l_{0}\varphi _{x} = \left [ -2{ \xi }^{3}+3{\xi }^{2},\; {\xi }^{3}-2{\xi }^{2}+\xi ,\; {\xi }^{3}-{\xi }^{2} \right ] \left [ \textstyle\begin{array}{c} \!\!\varepsilon _{2}\!\! \\ \!\!\varepsilon _{7}\!\! \\ \!\!\varepsilon _{8}\!\! \end{array}\displaystyle \right ]\,, \\ l_{0}\varphi _{y}\hspace{-0.8em} &= \left [ -3{\xi }^{2}+4\xi -1,\; 3{\xi }^{2}-2\xi \right ] \left [ \textstyle\begin{array}{c} \!\!\varepsilon _{3}^{s}\!\! \\ \!\!\varepsilon _{4}^{s}\!\! \end{array}\displaystyle \right ]\,, \quad l_{0}\varphi _{z} = \left [ -3{\xi }^{2}+4\xi -1,\; 3{ \xi }^{2}-2\xi \right ] \left [ \textstyle\begin{array}{c} \!\!\varepsilon _{5}^{s}\!\! \\ \!\!\varepsilon _{6}^{s}\!\! \end{array}\displaystyle \right ]\,, \vspace{0.5em} \end{array} $$

(52)

where the deformations ($\varepsilon _{3}^{s},\dots ,\varepsilon _{6}^{s}$) are defined in Eq. (28a). With the interpolations in Eq. (52), the following relationship between the vectors $\boldsymbol{\theta }$ and $\boldsymbol{\varepsilon }^{s}$ can be established:

$$ \boldsymbol{\theta }=\mathbf{B}\,\boldsymbol{\varepsilon }^{s}\,, $$

(53)

where

$$ \displaystyle \mathbf{B} = \frac{1}{l_{0}^{2}} \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} l_{0}& 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -6{\xi }^{2}+6{\xi } & 0 & 0 & 0 & 0 & 3{\xi }^{2}-4{\xi }+1 & 3{\xi }^{2}-2{\xi } \\ 0 & 0 & -6{\xi }+4 & 6{\xi }-2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -6{\xi }+4 & 6{\xi }-2 & 0 & 0 \\ 0 & (-12{\xi }+6)/l_{0} & 0 & 0 &0 & 0 & (6{\xi }-4)/l_{0} & (6{\xi }-2)/l_{0}\\ \end{array}\displaystyle \right ]. $$

By substituting Eq. (53) into (51), we obtain the stiffness equations

$$ \boldsymbol{\sigma }^{s}=\big(\mathbf{S}^{s}+\mathbf{G}^{s}\big)\, \boldsymbol{\varepsilon }^{s}\,, $$

(54)

where

$$\begin{aligned} \displaystyle \mathbf{S}^{s} &=l_{0}\int _{0}^{1}\mathbf{B}^{\mathrm{T}}\big(\mathbf{H}^{\mathrm{T}}\,\mathbf{D}\,\mathbf{H}\big) \mathbf{B}\,\mathrm{d}\xi \,, \end{aligned}$$

(55a)

$$\begin{aligned} \mathbf{G}^{s}(\varepsilon _{2},\varepsilon _{7},\varepsilon _{8}) &=l_{0} \int _{0}^{1}\mathbf{B}^{\mathrm{T}}\big(\mathbf{A}^{\mathrm{T}} \,\mathbf{D}\,\mathbf{H} +\frac{1}{2}(\mathbf{H}^{\mathrm{T}}\, \mathbf{D}\,\mathbf{A}+\mathbf{A}^{\mathrm{T}}\,\mathbf{D}\,\mathbf{A} \big)\big)\mathbf{B}\,\mathrm{d}\xi \,. \end{aligned}$$

(55b)

The first matrix $\mathbf{S}^{s}$ is the elastic stiffness matrix. The linearized stiffness equations are expressed by the relationship

$$ \Delta \boldsymbol{\sigma }^{s}=\mathbf{S}^{s}_{\mathrm{T}}\;\Delta \boldsymbol{\varepsilon }^{s}\,, $$

(56)

where $\mathbf{S}^{s}_{\mathrm{T}}$ is the tangent stiffness matrix, defined by

$$ \mathbf{S}^{s}_{\mathrm{T}}=\mathbf{S}^{s}+\mathbf{G}^{s}_{\mathrm{T}}\,. $$

(57)

The second matrix is the geometric stiffness matrix $\mathbf{G}^{s}_{\mathrm{T}}$ defined by

$$ \mathbf{G}^{s}_{\mathrm{T}}\big(\varepsilon _{2},\varepsilon _{7}, \varepsilon _{8})=l_{0}\int _{0}^{1}\mathbf{B}^{\mathrm{T}} \big(\mathbf{A}^{\mathrm{T}}\,\mathbf{D}\,\mathbf{H} +\mathbf{H}^{\mathrm{T}}\,\mathbf{D}\,\mathbf{A}+\frac{3}{2}\mathbf{A}^{\mathrm{T}}\, \mathbf{D}\,\mathbf{A}\big)\mathbf{B}\,\mathrm{d}\xi \,. \vspace{-0.2em} $$

(58)

The components of the $(8\times 8)$ matrices $\mathbf{S}^{s}$ and $\mathbf{G}^{s}_{\mathrm{T}}$ can be found in Appendix A. The vectors $\boldsymbol{\sigma }$ and $\boldsymbol{\sigma }^{s}$ and their variations are related by

$$ \boldsymbol{\sigma } = \mathbf{E}^{s\mathrm{T}}\,\boldsymbol{\sigma }{\,^{s}}\,, \qquad \Delta \boldsymbol{\sigma } = \mathbf{E}^{s\mathrm{T}}\,\Delta \boldsymbol{\sigma }{\,^{s}}\,, $$

(59)

where $\mathbf{E}^{s}$ is the transformation matrix defined in Eqs. (28a), (28b). By substituting Eqs. (28a), (28b) into Eqs. (54) and (56), we obtain with Eq. (59)

$$ \boldsymbol{\sigma } = (\mathbf{S}+\mathbf{G})\,\boldsymbol{\varepsilon }\,,\qquad \Delta \boldsymbol{\sigma } = \mathbf{S}_{\mathrm{T}}\,\Delta \boldsymbol{\varepsilon }\,, $$

(60)

where

$$ \mathbf{S} = \mathbf{E}^{s\mathrm{T}}\,\mathbf{S}{\,^{s}}\,\mathbf{E}^{s} \,,\qquad \mathbf{G} = \mathbf{E}^{s\mathrm{T}}\,\mathbf{G}{\,^{s}}\, \mathbf{E}^{s}\,,\qquad \mathbf{S}_{\mathrm{T}} = \mathbf{E}^{s \mathrm{T}}\,\mathbf{S}^{s}_{\mathrm{T}}\,\mathbf{E}^{s}\,. $$

(61)

These are the element stiffness equations expressed in of the centroid coordinates. Note that in the initial undeformed configuration $\mathbf{G}(0)=0$ and $\mathbf{S}_{\mathrm{T}}(0)=0$. Therefore, the matrices $\mathbf{G}$ and $\mathbf{S}_{\mathrm{T}}$ are not decisive in a linear buckling analysis.

5 Modified deformations

In this section a set of modified deformations is derived in which second order approximations for the axial elongation and bending curvatures represented by the strain vector $\boldsymbol{\gamma }_{c}$ are accounted for by additional quadratic terms in the expressions of the basic deformations. 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 $\varepsilon _{i}\, (i=2,\dots ,6)$ of Eqs. (8a)–(8d) can be substituted for the first order deformations $\varepsilon ^{\ell }_{i}\, (i=2,\dots ,6)$ of Eqs. (24a)–(24c). Then Eq. (52) can be rewritten in terms of second order deformations as:

$$ \textstyle\begin{array}{ll} \displaystyle \phantom{-} \overline{u}_{c} &= \xi \varepsilon _{1}\,, \\ \displaystyle l_{0}\overline{\varphi }_{x} &= \left [ -2{\xi }^{3}+3{ \xi }^{2},\; {\xi }^{3}-2{\xi }^{2}+\xi ,\; {\xi }^{3}-{\xi }^{2} \right ] \left [ \textstyle\begin{array}{l} \displaystyle \varepsilon _{2}-\textstyle \frac{1}{2l_{0}}( \varepsilon _{3}-\varepsilon _{4})(\varepsilon _{5}+\varepsilon _{6}) \\ \varepsilon _{7} \\ \varepsilon _{8} \end{array}\displaystyle \right ]\,, \\ \displaystyle l_{0}\overline{\varphi }_{y}&= \left [ -3{\xi }^{2}+4\xi -1, \; 3{\xi }^{2}-2\xi \right ] \left [ \textstyle\begin{array}{l} \varepsilon _{3}^{s} \\ \varepsilon _{4}^{s}-\frac{1}{l_{0}}\varepsilon _{2}\, \varepsilon ^{s}_{6} \end{array}\displaystyle \right ]\,, \\ \displaystyle l_{0}\overline{\varphi }_{z}&= \left [ -3{\xi }^{2}+4\xi -1, \; 3{\xi }^{2}-2\xi \right ] \left [ \textstyle\begin{array}{l} \varepsilon _{5}^{s} \\ \varepsilon _{6}^{s}+\frac{1}{l_{0}}\varepsilon _{2}\, \varepsilon ^{s}_{4} \end{array}\displaystyle \right ]\,, \\ \displaystyle \phantom{-} \overline{v}_{c} &= \left [-{\xi }^{3}+2{\xi }^{2}-\xi ,\;\; {\xi }^{3}-{ \xi }^{2} \right ] \left [ \textstyle\begin{array}{l} \varepsilon _{5} \\ \varepsilon _{6}+\frac{1}{l_{0}}\varepsilon _{2}\,\varepsilon _{4} \end{array}\displaystyle \right ]\,, \\ \displaystyle \phantom{-} \overline{w}_{c} &= \left [ {\xi }^{3}-2{\xi }^{2}+\xi ,\;\; -{\xi }^{3}+{ \xi }^{2}\right ] \left [ \textstyle\begin{array}{l} \varepsilon _{3} \\ \varepsilon _{4}-\frac{1}{l_{0}}\varepsilon _{2}\,\varepsilon _{6} \end{array}\displaystyle \right ]\,. \end{array} $$

(62)

Integrating Eq. (38) over the length of the beam using the displacement and rotation fields of Eq. (62) gives

$$ \textstyle\begin{array}{ll} \displaystyle \overline{\varepsilon }_{1}=\int \limits _{0}^{1} \overline{u}_{c,\xi }\,\mathrm{d}\xi & \displaystyle + \frac{\;1}{2l_{0}}\int \limits _{0}^{1} \big(\overline{v}_{c,\xi }^{2}+ \overline{w}_{c,\xi }^{2} +\frac{I_{p}}{A}\,\overline{\varphi }_{x,\, \xi }^{2}\big) \mathrm{d}\xi -\Delta \varepsilon _{1} = \\ \displaystyle \hspace{5em} \,\varepsilon _{1} & \displaystyle +\frac{1}{30l_{0}}\big(2\, \varepsilon _{3}^{2} + \varepsilon _{3}\,\varepsilon _{4} + 2\, \varepsilon _{4}^{2} + 2\,\varepsilon _{5}^{2} + \varepsilon _{5}\, \varepsilon _{6} + 2\,\varepsilon _{6}^{2}\big) \\ & \displaystyle +\frac{1}{30l_{0}}\,\frac{I_{p}}{Al_{0}^{2}}\Big(3\, \varepsilon _{2}\big(6\varepsilon _{2}-\varepsilon _{7}-\varepsilon _{8} \big) +2\varepsilon _{7}^{2} - \varepsilon _{7}\,\varepsilon _{8} + 2 \varepsilon _{8}^{2}\Big)-\Delta \varepsilon _{1}\,, \end{array} $$

(63)

where $\varepsilon _{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 $\Delta \varepsilon _{1}$ represents the over-calculated second order elongation due to additional bending deformations caused by a twist rotation of the cross-section if the shear centre does not coincide with the centroid. This is illustrated in Fig. 5 which shows that the elastic line becomes helical due to a twist rotation $\varepsilon _{2}/l_{0}$ of 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

$$ \Delta \varepsilon _{1} \displaystyle =\frac{1}{30l_{0}}\big(2\, \varepsilon _{3}^{2} + \varepsilon _{3}\,\varepsilon _{4} + 2\, \varepsilon _{4}^{2} + 2\,\varepsilon _{5}^{2} + \varepsilon _{5}\, \varepsilon _{6} + 2\,\varepsilon _{6}^{2}\big) =\frac{1}{10l_{0}} \big({\hat{y}_{s}}^{2}+{\hat{z}_{s}}^{2}\big)\,\varepsilon _{2}^{2}\,. $$

(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 deflection without first finding the complete equation of the deflection 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:

$$ \textstyle\begin{array}{ll} \displaystyle \overline{\varepsilon }_{3}^{s}= l_{0}\int \limits _{0}^{1} \,\big(\overline{\varphi }_{y,\xi }+\overline{\varphi }_{x} \overline{\varphi }_{z,\xi }\big)(1-\xi )\mathrm{d}\xi \,=&\, \varepsilon _{3}^{s} \displaystyle +\frac{\varepsilon _{2}}{10l_{0}} \big(\varepsilon _{5}^{s} + 2\varepsilon _{6}^{s}\big) + \frac{1}{30l_{0}} \big(3\varepsilon _{5}^{s}\,\varepsilon _{7} - \varepsilon _{8}(\varepsilon _{5}^{s}+\varepsilon _{6}^{s})\big)\,, \\ \hspace{2em} \overline{\varepsilon }_{4}^{s}= \displaystyle l_{0}\int \limits _{0}^{1} \,\big(\overline{\varphi }_{y,\xi }+\overline{\varphi }_{x} \overline{\varphi }_{z,\xi }\big)\,\xi \,\mathrm{d}\xi \,=&\, \varepsilon _{4}^{s} \displaystyle -\frac{\varepsilon _{2}}{10l_{0}} \big(2\varepsilon _{5}^{s} + \varepsilon _{6}^{s}\big) - \frac{1}{30l_{0}} \big(3\varepsilon _{6}^{s}\,\varepsilon _{8} - \varepsilon _{7}(\varepsilon _{5}^{s} + \varepsilon _{6}^{s})\big)\,, \\ \overline{\varepsilon }_{5}^{s}= \displaystyle l_{0}\int \limits _{0}^{1} \,\big(\overline{\varphi }_{z,\xi }-\overline{\varphi }_{x} \overline{\varphi }_{y,\xi }\big)(1-\xi )\mathrm{d}\xi \,=&\, \varepsilon _{5}^{s} \displaystyle -\frac{\varepsilon _{2}}{10l_{0}} \big(\varepsilon _{3}^{s} + 2\varepsilon _{4}^{s}\big) - \frac{1}{30l_{0}} \big(3\varepsilon _{3}^{s}\,\varepsilon _{7} - \varepsilon _{8}(\varepsilon _{3}^{s}+\varepsilon _{4}^{s})\big)\,, \\ \hspace{2em} \overline{\varepsilon }_{6}^{s}= \displaystyle l_{0}\int \limits _{0}^{1} \,\big(\overline{\varphi }_{z,\xi }-\overline{\varphi }_{x} \overline{\varphi }_{y,\xi }\big)\,\xi \,\mathrm{d}\xi \,=&\, \varepsilon _{6}^{s} \displaystyle +\frac{\varepsilon _{2}}{10l_{0}} \big(2\varepsilon _{3}^{s} + \varepsilon _{4}^{s}\big) + \frac{1}{30l_{0}} \big(3\varepsilon _{4}^{s}\,\varepsilon _{8} - \varepsilon _{7}(\varepsilon _{3}^{s} + \varepsilon _{4}^{s})\big)\,. \end{array} $$

(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. $\overline{\boldsymbol{\varepsilon }}^{s} = \mathbf{E^{s}}\, \overline{\boldsymbol{\varepsilon }}$. Substituting Eq. (27) into Eq. (65) yields, with Eqs. (63) and (64), the following set of modified deformations:

$$ \displaystyle \overline{\boldsymbol{\varepsilon }} = \boldsymbol{E}(\boldsymbol{\varepsilon }), $$

(66)

where

$$ \textstyle\begin{array}{l} \left.\textstyle\begin{array}{l} \displaystyle \overline{\varepsilon }_{1} = \varepsilon _{1} +\frac{1}{30l_{0}}\big(2 \,\varepsilon _{3}^{2} + \varepsilon _{3}\,\varepsilon _{4} + 2\, \varepsilon _{4}^{2} + 2\,\varepsilon _{5}^{2} + \varepsilon _{5}\, \varepsilon _{6} + 2\,\varepsilon _{6}^{2}\big) \\ \displaystyle \hspace{1.1em} +\frac{1}{30l_{0}}\,\frac{I_{p}}{Al_{0}^{2}}\Big(3\,\varepsilon _{2}(6 \varepsilon _{2}-\varepsilon _{7}-\varepsilon _{8}) +2\varepsilon _{7}^{2} - \varepsilon _{7}\,\varepsilon _{8} + 2\varepsilon _{8}^{2}\Big) \\ \displaystyle \hspace{1.1em} -\frac{1}{10l_{0}}\big(\,{\hat{y}_{s}}^{2}+{\hat{z}_{s}}^{2}\,\big)\, \varepsilon _{2}^{2}\, \\ \end{array}\displaystyle \quad \right \} \mathrm{(elongation),} \\ \end{array} $$

(67a)

$$ \textstyle\begin{array}{l} \left.\textstyle\begin{array}{l} \displaystyle \overline{\varepsilon }_{3} =\varepsilon _{3} + \frac{\varepsilon _{2}}{10l_{0}}\big(\varepsilon _{5} + 2\varepsilon _{6}+ \hat{z}_{s}(\varepsilon _{2}-\varepsilon _{7})\big) \displaystyle + \frac{1}{30l_{0}} \big(3\varepsilon _{5}\,\varepsilon _{7}- \varepsilon _{8}(\varepsilon _{5}+\varepsilon _{6})\big)\, \\ \displaystyle \overline{\varepsilon }_{4} =\varepsilon _{4}- \frac{\varepsilon _{2}}{10l_{0}}\big(2\varepsilon _{5}+\varepsilon _{6}- \hat{z}_{s}(\varepsilon _{2}-\varepsilon _{8})\big) \displaystyle - \frac{1}{30l_{0}} \big(3\varepsilon _{6}\,\varepsilon _{8}- \varepsilon _{7}(\varepsilon _{5}+\varepsilon _{6})\big)\, \\ \end{array}\displaystyle \quad \right \} \textrm{bending~($x$--$z$),} \\ \end{array} $$

(67b)

$$ \textstyle\begin{array}{l} \left.\textstyle\begin{array}{l} \displaystyle \overline{\varepsilon }_{5}=\varepsilon _{5} - \frac{\varepsilon _{2}}{10l_{0}}\big(\varepsilon _{3} + 2\varepsilon _{4}+ \hat{y}_{s}(\varepsilon _{2}-\varepsilon _{7})\big) \displaystyle - \frac{1}{30l_{0}}\big(3\varepsilon _{3}\,\varepsilon _{7}- \varepsilon _{8}(\varepsilon _{3}+\varepsilon _{4})\big)\, \\ \displaystyle \overline{\varepsilon }_{6}=\varepsilon _{6}+ \frac{\varepsilon _{2}}{10l_{0}}\big(2\varepsilon _{3}+\varepsilon _{4}- \hat{y}_{s}(\varepsilon _{2}-\varepsilon _{8})\big) \displaystyle + \frac{1}{30l_{0}} \big(3\varepsilon _{4}\,\varepsilon _{8}- \varepsilon _{7}(\varepsilon _{3}+\varepsilon _{4})\big)\, \\ \end{array}\displaystyle \quad \right \} \textrm{bending~($x$--$y$),} \\ \end{array} $$

$$ \textstyle\begin{array}{l@{\quad }l} \textstyle\begin{array}{l} \displaystyle \overline{\varepsilon }_{7} = \varepsilon _{7}\, \\ \overline{\varepsilon }_{8} = \varepsilon _{8} \end{array}\displaystyle & \biggl\} \mathrm{(warping).} \end{array} $$

(67c)

Because the rigidity against axial elongation is much larger than the flexural and torsional rigidity, the modification of $\overline{\varepsilon }_{1}$ is the most relevant one. The other modifications are relevant 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 deformations expressed as functions of the global nodal coordinates $X_{i}$, namely

$$ \overline{\boldsymbol{\varepsilon }} = \boldsymbol{E}(\boldsymbol{D}(\boldsymbol{X})) = \overline{\boldsymbol{D}}(\boldsymbol{X})\,. $$

(68)

The modified deformations form the basis of a second order beam finite element that captures non-uniform torsion and flexural–torsional coupling of thin-walled beams with unsymmetrical 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 $(x,y,z)$ attached at node $p$ described by the vector $\boldsymbol{r}^{\,p}$ and a set of four Euler parameters $(\lambda _{0}^{p},\boldsymbol{\lambda }^{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 $\boldsymbol{x} = (0,\,y,\,z)^{\mathrm{T}}$ denote the position vector of an arbitrary point P of the cross-section, in the initial undeformed configuration and let $\alpha ^{p}\boldsymbol{\omega }$ be the warping displacement perpendicular to the cross-section at node $p$, where $\boldsymbol{\omega }\,= \big (\omega (y,z),0,0\big )^{\mathrm{T}}$, with $\omega (y,z)$ being a normalized warping function and $\alpha ^{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

$$ \boldsymbol{r}^{\,\mathrm{P}}=\boldsymbol{r}^{\,p}+\mathbf{R}^{p}\boldsymbol{R}_{0}\big( \boldsymbol{x}+\alpha ^{p}\boldsymbol{\omega }\big), $$

(69)

where $\mathbf{R}^{p}\mathbf{R}_{0}$ is a transformation matrix from the local coordinate system to the inertial frame. The matrices $\mathbf{R}_{0}$ and $\mathbf{R}^{p}$ are 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

$$ \dot{\boldsymbol{r}}^{\,\mathrm{P}}=\dot{\boldsymbol{r}}^{\,p}+\dot{\mathbf{R}}^{p} \mathbf{R}_{0}\big(\boldsymbol{x} +\alpha ^{p}\boldsymbol{\omega }\big)+\mathbf{R}^{p} \mathbf{R}_{0}\,\dot{\alpha }^{p}\boldsymbol{\omega }\,, $$

(70)

where $(\overset{\cdot }{})$ denotes differentiation with respect to time. This equation can be written in terms of the time derivatives of $(\lambda _{0}^{p},\boldsymbol{\lambda }^{p})$ and $\alpha ^{p}$ as

$$ \dot{\boldsymbol{r}}^{\,\mathrm{P}}=\dot{\boldsymbol{r}}^{\,p}+\mathbf{B}^{p} \left [ \textstyle\begin{array}{c} \dot{\lambda }_{0}^{p} \\ \dot{\boldsymbol{\lambda }}^{p} \end{array}\displaystyle \right ]+ \mathbf{R}^{p}\mathbf{R}_{0}\,\dot{\alpha }^{p}\boldsymbol{\omega }\,, $$

(71)

where

$$ \mathbf{B}^{p}=-2\,\mathbf{R}^{p}\mathbf{R}_{0}\,[\tilde{\boldsymbol{x}}+ \alpha ^{p}\tilde{\boldsymbol{\omega }}]\,\mathbf{R}_{0}^{\mathrm{T}}\, \overline{\boldsymbol{\varLambda }}^{\,p}, $$

(72)

in which $\tilde{\boldsymbol{x}}$ and $\tilde{\boldsymbol{\omega }}$ are skew-symmetric matrices defined by [41]

$$ \tilde{\boldsymbol{x}} = \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \phantom{-}0 & -z & \phantom{-} y \\ \phantom{-}z & \phantom{-}0 & \phantom{-}0 \\ -y & \phantom{-}0 & \phantom{-} 0 \end{array}\displaystyle \right ] \quad \text{and}\quad \tilde{\boldsymbol{\omega }} = \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} 0 & \phantom{-} 0 & \phantom{-} 0 \\ 0 & \phantom{-} 0 &-\omega \\ 0 & \phantom{-} \omega & \phantom{-} 0 \end{array}\displaystyle \right ]. $$

(73)

Matrix $\overline{\boldsymbol{\varLambda }}^{\,p}$ is defined by

$$ \displaystyle \overline{\boldsymbol{\varLambda }}^{\,p}=[ -\boldsymbol{\lambda }^{p},\, \lambda _{0}^{p}\boldsymbol{I} - \tilde{\boldsymbol{\lambda }}^{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

$$ \displaystyle T^{p}={\textstyle \frac{1}{4}}\rho l_{0}\,\int _{A} \dot{\boldsymbol{r}}^{\,\mathrm{P}\mathrm{T}}\dot{\boldsymbol{r}}^{\,\mathrm{P}}dA\,, $$

(75)

where $\rho$ and $A$ are, respectively, the mass density and the area of the cross-section. Substituting Eq. (71) into Eq. (75) yields

$$ \displaystyle T^{p}={\textstyle \frac{1}{4}}\rho l_{0}\!\!\int _{A}\Big[\,\dot{\boldsymbol{r}}^{\,p\mathrm{T}}\!, [\dot{\lambda }_{0}^{p}, \dot{\boldsymbol{\lambda }}^{p\mathrm{T}}],\dot{\alpha }^{p} \Big] \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \mathbf{I} & \mathbf{B}^{p} & \mathbf{R}^{p}\mathbf{R}_{0}\, \boldsymbol{\omega } \\ &\mathbf{B}^{p\mathrm{T}}\mathbf{B}^{p}& \mathbf{B}^{p\mathrm{T}} \mathbf{R}^{p}\mathbf{R}_{0}\,\boldsymbol{\omega } \\ &\mathrm{symm.\hspace{1em}} & \boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\omega } \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c} \dot{\boldsymbol{r}}^{\,p} \\ \left [ \textstyle\begin{array}{c} \dot{\lambda }_{0}^{p} \\ \dot{\boldsymbol{\lambda }}^{p} \end{array}\displaystyle \right ] \\ \dot{\alpha }^{p} \end{array}\displaystyle \right ]dA. $$

(76)

If the origin of frame $(x,y,z)$ coincides with the centre of mass of the body, the coupling term $\boldsymbol{B}^{p}$ disappears. Since the matrices $\overline{\boldsymbol{\varLambda }}^{p}$ and $\mathbf{R}^{p}$, $\mathbf{R}_{0}$ are not space dependent, we obtain using Eq. (31)

$$ \int _{A}\mathbf{B}^{p\mathrm{T}}\mathbf{B}^{p}\,dA = 4\; \overline{\boldsymbol{\varLambda }}^{\,p\mathrm{T}}\overline{\mathbf{J}}^{\;p}\; \overline{\boldsymbol{\varLambda }}^{\,p}, \quad \int _{A}\mathbf{R}^{p} \mathbf{R}_{0}\,\boldsymbol{\omega }\,dA=\mathbf{0}\,, \quad \int _{A} \boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\omega }\,dA=I_{\omega }\,. $$

(77)

Substituting Eq. (77) into Eq. (76) yields the lumped mass matrix

$$ \displaystyle \mathbf{M}_{\ell }={\textstyle \frac{1}{2}}\rho l_{0}\; \mathrm{diag}\big(A\mathbf{I}\,,\;4\,\overline{\boldsymbol{\varLambda }}^{\,p \mathrm{T}}\overline{\mathbf{J}}^{\;p}\;\overline{\boldsymbol{\varLambda }}^{\,p}, \,I_{\omega }\,,\; A\mathbf{I}\,,\;4\,\overline{\boldsymbol{\varLambda }}^{\,q \mathrm{T}}\overline{\mathbf{J}}^{\;q}\;\overline{\boldsymbol{\varLambda }}^{\,q}, \,I_{\omega } \big)\,, $$

(78)

in which $\mathbf{I}$ is a $3\times 3$ unitary matrix, $I_{\omega }$ is the sectorial moment defined in Eq. (46) and $\overline{\mathbf{J}}^{\;p}\!$, $\overline{\mathbf{J}}^{\;q}$ are defined by

$$ \overline{\mathbf{J}}^{\;p}\!\! = \mathbf{R}_{0}\big[\, \overline{\mathbf{I}}+(\alpha ^{p})^{2}\,\overline{\mathbf{I}}_{ \omega }\big]\,\mathbf{R}_{0}^{\mathrm{T}}\,,\qquad \overline{\mathbf{J}}^{\;q}\!\! = \mathbf{R}_{0}\big[\, \overline{\mathbf{I}}+(\alpha ^{q})^{2}\,\overline{\mathbf{I}}_{ \omega }\big]\,\mathbf{R}_{0}^{\mathrm{T}}\,, $$

(79)

where $\overline{\mathbf{I}} = \mathrm{diag}(I_{\mathrm{p}},\;I_{y},\;I_{z})$ and $\overline{\mathbf{I}}_{\omega }\,=\,\mathrm{diag}(0,\;{I}_{\omega },\;I_{ \omega })$. Using Lagrange’s equations, the corresponding convective vector $\boldsymbol{h}_{\ell }$ can be expressed as

$$ \displaystyle \boldsymbol{h}_{\ell }={\textstyle {\textstyle \frac{1}{2}}}\rho l_{0} \left [ \textstyle\begin{array}{c} \boldsymbol{0} \\ \left [ \textstyle\begin{array}{c} 8\big(\,\dot{\overline{\boldsymbol{\varLambda }}}^{\,p\,\mathrm{T}} \overline{\mathbf{J}}^{\;p}\;{\overline{\boldsymbol{\varLambda }}}^{\,p}\! +\; \alpha ^{p}\dot{\alpha }^{p}\;\overline{\boldsymbol{\varLambda }}^{\,p\,\mathrm{T}} \,\overline{\mathbf{I}}_{\omega }\;\overline{\boldsymbol{\varLambda }}^{\,p}\, \big) \\ -2\,\alpha ^{p}\big[\dot{\lambda }_{0}^{p},\dot{\boldsymbol{\lambda }}^{\,p\, \mathrm{T}}\big]~ \overline{\boldsymbol{\varLambda }}^{\,p\,\mathrm{T}}\, \overline{\mathbf{I}}_{\omega }\;\overline{\boldsymbol{\varLambda }}^{\,p} \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c} \dot{\lambda }_{0}^{p} \\ \dot{\boldsymbol{\lambda }}^{p} \end{array}\displaystyle \right ] \\ \boldsymbol{0} \\ \left [ \textstyle\begin{array}{c} 8\big(\,\dot{\overline{\boldsymbol{\varLambda }}}^{\,q\,\mathrm{T}} \overline{\mathbf{J}}^{\;q}\;{\overline{\boldsymbol{\varLambda }}}^{\,q}\! +\; \alpha ^{q}\dot{\alpha }^{q}\;\overline{\boldsymbol{\varLambda }}^{\,q\,\mathrm{T}} \,\overline{\mathbf{I}}_{\omega }\;\overline{\boldsymbol{\varLambda }}^{\,q}\, \big) \\ -2\,\alpha ^{q}\big[\dot{\lambda }_{0}^{q},\dot{\boldsymbol{\lambda }}^{\,q\, \mathrm{T}}\big]~ \overline{\boldsymbol{\varLambda }}^{\,q\,\mathrm{T}}\, \overline{\mathbf{I}}_{\omega }\;\overline{\boldsymbol{\varLambda }}^{\,q} \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c} \dot{\lambda }_{0}^{q} \\ \dot{\boldsymbol{\lambda }}^{q} \end{array}\displaystyle \right ] \end{array}\displaystyle \right ]\,. $$

(80)

Fig. 7

Position and orientation of cross-section at node $p$

7 Equations of motion

The equations of motion are derived using d’Alembert’s principle, which states that the balance of virtual work,

$$ \displaystyle \overline{\boldsymbol{\sigma }}^{\mathrm{T}} \delta \overline{\boldsymbol{\varepsilon }} = \boldsymbol{F}^{\mathrm{T}} \delta \boldsymbol{u}-\delta \boldsymbol{X}^{\mathrm{T}}\big(\mathbf{M}\,\ddot{\boldsymbol{X}} + \boldsymbol{h}\big)\, ,$$

(81)

holds for all $\delta \boldsymbol{\overline{\varepsilon }}$ and $\delta \boldsymbol{u}$ depending on the variations of the nodal coordinates $\delta \boldsymbol{X}$, by the relations

$$ \displaystyle \delta \overline{\boldsymbol{\varepsilon }} = \overline{\boldsymbol{D}}_{ \boldsymbol{,\,X}} \delta \boldsymbol{X}\,, \quad \boldsymbol{\delta u} = \boldsymbol{A}\, \boldsymbol{\delta X}\,, $$

(82)

where

$$ \displaystyle \overline{\boldsymbol{D}}_{\boldsymbol{,\,X}} = \boldsymbol{E}_{\boldsymbol{,\,D}}\; \boldsymbol{D}_{\boldsymbol{,\,X}}\;, \qquad \displaystyle \boldsymbol{A} = \mathrm{diag} \big(\mathbf{I},\;2\boldsymbol{\varLambda }^{p},\;1,\;\mathbf{I},\;2\boldsymbol{\varLambda }^{q}, \;1\big)\,, $$

(83)

in which $\mathbf{I}$ is a $3\times 3$ unitary matrix and $\boldsymbol{\varLambda }$ is defined in Eq. (11). The last term in Eq. (81) represents the virtual work due to the inertia forces, where $\mathbf{M}(\boldsymbol{X})$ is a position-dependent mass matrix and $\boldsymbol{h}(\boldsymbol{X},\dot{\boldsymbol{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. $\mathbf{M}=\mathbf{M}_{c}+\mathbf{M}_{\ell }$ and $\boldsymbol{h}=\boldsymbol{h}_{c}+\boldsymbol{h}_{\ell }$. Substituting the compatibility relations (82) into Eq. (81) yields, with the transposed matrices $\overline{\boldsymbol{D}}_{,\,\boldsymbol{X}}^{\mathrm{T}}$ and $\boldsymbol{A}^{\mathrm{T}}$,

$$ \displaystyle \overline{\boldsymbol{D}}_{,\,\boldsymbol{X}}^{\mathrm{T}}\, \overline{\boldsymbol{\sigma }} = \boldsymbol{A}^{\mathrm{T}}\boldsymbol{F} - \big(\mathbf{M}\, \ddot{\boldsymbol{X}} + \boldsymbol{h}\big) \,. $$

(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 $\overline{\varepsilon }_{i}$, relationships (60) between deformations and discrete stress resultants remain valid, i.e.

$$ \displaystyle \overline{\boldsymbol{\sigma }} = (\mathbf{S}+\mathbf{G})\, \overline{\boldsymbol{\varepsilon }}\,,\qquad \Delta \overline{\boldsymbol{\sigma }} = \mathbf{S}_{\mathrm{T}}\,\Delta \overline{\boldsymbol{\varepsilon }}\,. $$

(85)

The stress-resultants $\overline{\sigma }_{i}$ have now a slightly modified meaning in case of finite deformations.

In order to study the vibration properties and the elastic stability, we consider small disturbances from an equilibrium configuration

$$ \displaystyle \overline{\boldsymbol{D}}_{,\,\boldsymbol{X}}^{\mathrm{T}}\, \overline{\boldsymbol{\sigma }}^{\mathrm{o}} = \boldsymbol{A}^{\mathrm{T}}\boldsymbol{F}^{\mathrm{o}} \,. $$

(86)

Substituting (85) into (84), expanding the equations in terms of disturbances $\Delta \boldsymbol{X}$, $\Delta \ddot{\boldsymbol{X}}$ and disregarding second and higher order terms yields the linearized equations of motion

$$ \displaystyle \mathbf{M}^{o}\,\Delta \ddot{\boldsymbol{X}} +\big[\, \overline{\boldsymbol{D}}_{,\,\boldsymbol{X}}^{\mathrm{oT}}\,\mathbf{S}_{\mathrm{T}}\, \overline{\boldsymbol{D}}_{,\,\boldsymbol{X}}^{\mathrm{o}} +\overline{\boldsymbol{D}}_{,\, \boldsymbol{X}\boldsymbol{X}}^{\mathrm{oT}}\,\overline{\boldsymbol{\sigma }}^{\mathrm{o}} -(\boldsymbol{A}^{\mathrm{oT}}\boldsymbol{F}^{\mathrm{o}})_{,\,\boldsymbol{X}}\big]\,\Delta \boldsymbol{X} = \mathbf{0}\,, $$

(87)

or in matrix form

$$ \displaystyle \mathbf{M}^{o}\,\Delta \ddot{\boldsymbol{X}}+\big(\mathbf{K}^{o} + \mathbf{G}^{o}\big)\,\Delta \boldsymbol{X} = \mathbf{0}\,, $$

(88)

where $\mathbf{K}^{o}$ is the structural stiffness matrix and $\mathbf{G}^{o}$ the structural geometric stiffness matrix due to the reference load $\boldsymbol{F}^{o}$ giving rise to reference stress resultants $\overline{\boldsymbol{\sigma }}^{\mathrm{o}}$. Note that $(\mathbf{K}^{o} + \mathbf{G}^{o} )$ 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 modified 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 elements. The motion of the system can be simulated for given initial conditions, the equations of motion can be linearized about an arbitrary state of motion, stationary solutions can be determined and with the linearized equations, eigenfrequencies and corresponding frequency modes as well as buckling loads can be determined. A Newton–Raphson method for calculating and continuing static solutions for flexible multibody systems [35] is employed for the solution of the quasi-static equilibrium equations (86). In order to validate the stiffness and inertia formulation, flexural–torsional (post-)buckling analyses and vibration analyses of beams with symmetric and asymmetric cross-sections are performed first. Finally, two multibody examples are presented, namely the post-buckling behaviour of a rotating cantilever beam with C-shaped cross-section and a stability analysis of a flexible cross-hinge mechanism.

8.1 Flexural–torsional buckling of a C-shaped beam under axial load

In this example we determine the global buckling load of a simply supported C-shaped beam subjected to an axial force $F$ at the centroid $c$, see Fig. 8. The cross-sectional properties of the beam are presented in Table 1. Because of symmetry only one-half of the beam need to be modelled. The half-beam is divided into 2, 4 and 16 BEAMW elements of equal length. In Table 2, the buckling loads are presented. When 16 finite elements are adopted, the error with respect to the analytic solution of Cortinez and Piovian [13] is negligible. Next, the eigenfrequencies of modes 1–6 and mode 12 have been computed and compared with the analytical solution values of Cortinez and Piovan [13]. The beam is divided into 30 BEAMW elements of equal length. The results are shown in Table 2. Rotary inertia is not included. There is a good agreement with the analytic results.

Fig. 8

Simply supported C-shaped beam with geometrical properties

Table 1 Cross-sectional properties

Table 2 Buckling loads (kN) and flexural–torsional frequencies (Hz)

8.2 Lateral-torsional post-buckling of an I-section cantilever

This example was originally proposed by Manta and Goncalves [33] and consists of an I-section cantilever subjected to a vertical force $F$ at the free end applied at the centroid $c$ of the cross-section. The cross-section and its geometrical properties are given in Fig. 9 (left) and Fig. 10 (right). The elastic constants are $E\,=\,20900$ kN/cm² and $G\,=\,8038$ kN/cm². A perturbation force $F_{Y}=\,0.01$ N is applied and kept constant, while force $F$ is increased to 16 kN. In Fig. 10 (left), the results obtained with the present beam finite element are compared with those in [33], obtained using 20 equal length two-node geometrically exact beam elements. It can be observed that using 20 equal length BEAMW elements gives an excellent match with those of [33]. The graph also shows that even with 10 equal length BEAMW elements already a result is obtained that practically match those obtained in [33]. To obtain a comparable result using BEAMWL elements, double the number of elements is needed, so 20 BEAMWL elements are needed to match the result obtained with 10 BEAMW elements. Figure 9 (right) shows the deformed configuration of the cantilever beam.

Fig. 9

(Left) Cross-section; (Right) Deformed configuration of I-section cantilever

Fig. 10

(Left) Load–displacement curves; (Right) Cross-sectional properties

8.3 C-shaped section cantilever beam loaded by an eccentric transverse force at the free end

This example was originally proposed by Gruttmann et al. [22] and concerns the lateral–torsional buckling of a C-shaped section cantilever shown in Fig. 11(a), subjected to an eccentric vertical force $F$ at the free end applied at point $D$ of the cross-section. The dimensions of the cross-section and material properties are presented in Table 3. Figure 12 shows the load–displacement curves of point $D$ for the case where force $F$ is increased to 20 kN. In Fig. 12(a) the results of the developed beam element are compared with those given in [22] and by Manta and Goncalves [33], both obtained using 30 equal-length two-node geometrically exact beam elements. The result obtained with 30 equal-length BEAMW elements shows a good agreement with these results. The graph also shows that using 10 equal-length BEAMW elements leads to a result that practically match these results. No visual difference could be observed between the results obtained with 30 BEAMW and 30 BEAMWL elements. Figure 12(b) shows that even with 6 BEAMW elements already a rather accurate result is obtained. From this figure a significant difference can be seen between the results obtained with 6 BEAMW and 6 BEAMWL elements. The latter corresponds fairly well with the 6 FE solution in [33]. This indicates that the second order beam formulation proposed in this paper requires a less number of elements to achieve the same accuracy. Figure 11(c) shows the deformed configuration. From this figure it can be observed that torsional buckling occurs locally. In this area element twist angles can be as high as 36^∘ when 6 BEAMW elements are used.

Fig. 11

Lateral-torsional buckling of C-shaped section cantilever

Fig. 12

Load–displacement curves

Table 3 Cross-sectional and material properties

8.4 Rotating C-shaped section cantilever beam loaded by an eccentric transverse force at the free end

In order to demonstrate the applicability of the present beam formulation for multibody systems, we consider the post-buckling analysis of the C-shaped section cantilever (in Example 8.3) attached to a rigid hub of radius $R=0$, rotating at constant angular velocity $\varOmega$, see Fig. 13 (left). The load–displacement curves of point $D$, where load $F$ is applied at point $D$, are shown in Fig. 13 (right) for four angular velocities: $\varOmega =0$, $\varOmega =1.3\,\omega _{0}$, $\varOmega =2\,\omega _{0}$ and $\varOmega =3\,\omega _{0}$. The value of $\omega _{0}$ is presented in Table 3. The beam is divided into 10 BEAMW elements of equal length. The analyses illustrate the capability of the present beam formulation to account for the effect of the centrifugal stiffening on the post-buckling behaviour of a rotating C-shaped section cantilever beam.

Fig. 13

(Left) Rotating cantilever beam; (Right) Load–displacement curves

8.5 Buckling analysis of a flexible cross-hinge mechanism

The final example concerns the buckling analysis of a flexible cross-hinge mechanism shown in Fig. 14. The mechanism consists of a shuttle and two attached leaf springs of equal dimensions (80 mm length) with narrow rectangular cross-sections (30 mm width, 0.35 mm thickness), made of steel ($E=210$ GPa, $\nu =0.29$). The shuttle guides motion with respect to the base about the indicated rotation axis; in the other directions, it constrains motion. A translational guidance in the vertical direction, represented by the rollers for the bottom leaf spring, makes the system is free to rotate when displacement $v_{0}=0$. Since the vertical motion of the guidance is associated with high support stiffness of the leaf springs, a small misalignment displacement $v_{0}$ leads to large internal stresses. The leaf springs are essentially loaded by a transverse force $F_{0}$ under fixed-guided boundary conditions. Figure 14 depicts the distribution of the von Mises stress in the leaf springs due to misalignment displacement $v_{0}$. The leaf springs are essentially loaded by a transverse force $F_{0}$ under fixed-guided boundary conditions. These stresses can cause a decrease of stiffness even up to a level at which the leaf springs buckle. A flexible multibody model with varying number of BEAMW elements per leaf spring is created in which the cross-hinge shuttle is assumed to be rigid. The only elastic effects that are modelled are due to the leaf spring deformations. The torsion and bending deformations in the plane of lowest rigidity are allowed; elongation and bending in the plane of highest rigidity are excluded. Figure 15 shows that the critical buckling load of the mechanism increases by 55$\%$ from a converged value of 58.6 N (free-warping) to 90.63 N when warping is modelled (constrained warping). An explanation for this result is that warping is constrained at the clamped and fixed-guided ends of the leaf springs, effectively reducing the length over which the leaf spring can twist; since the buckling mode exhibits torsion deformation, the buckling load increases. An analytical buckling analysis of the mechanism at hand [38] predicts corresponding buckling loads of 58.25 N (free-warping) and 90.63 N (constrained warping). Thus the critical load matches with those of the analytical model.

Fig. 14

Flexible cross-hinge mechanism

Fig. 15

Convergence of critical buckling load

9 Summary and conclusions

The generalized strain beam formulation is used to derive a three-dimensional beam element for pre- and post-buckling analysis of thin-walled open-section beams within a multibody based frame work. The elastic and geometric stiffness matrices are derived from a nonlinear theory of non-uniform torsion of thin-walled open section beams in which axial shortening and nonlinear Wagner’s stiffening torques are accounted for. Coupling among bending and torsion due to non-coinciding centroid and shear centre is incorporated using a transformation matrix which accounts for the eccentricity of the shear centre. Second order approximations for the axial elongation and bending curvatures are included using a set of modified deformations. A lumped mass formulation describing twist, rotary and warping inertia of the beam’s cross-section is presented. Numerical results of the proposed element have shown to be in excellent agreement with analytical and finite element solutions from literature. The following conclusions can be drawn:

• 1. Due to the definition of physically meaningful deformations and a sound inclusion of nonlinear geometrical effects, only a relatively small number of beam elements are needed to accurately model thin-walled beams undergoing large deflections and angles of twist.

• 2. By employing the concept of modified deformations, the existing relationships between discrete stress resultants and deformations remain valid and as a consequence a separation of stiffness due to elongation, torsion, warping and bending in two directions is retained, so deformations associated with a large stiffness can be eliminated by constraining them to be zero.

• 3. The effect of nonlinear bending curvatures represented by the quadratic terms in the modified deformations is investigated in the post-buckling analyses in Examples 8.2 and 8.3. It appears that including these terms reduces the number of finite elements required to accurately model thin-walled beams and thus allows a significant reduction in computation time.

• 4. The low computational cost of the present beam model makes it particularly attractive for optimization studies, since it opens the possibility for using cross-sectional parameters as design parameters without re-meshing the beam cross-section.

Appendix A: Stiffness matrices of beam element

The 8×8 structural stiffness matrix $\mathbf{S}^{\,s}$ in Eq. (55a) and geometric stiffness matrix $\mathbf{G}^{s}_{\mathrm{T}}$ in Eq. (58) of the BEAMW element based on the Hermitian interpolation is written in block matrix form as

where

$$ [S^{s}_{11}] = \frac{EA}{l_{0}}\,, \qquad \begin{bmatrix} S^{\,s}_{33} & S^{\,s}_{34} \\ S^{\,s}_{43} & S^{\,s}_{44} \end{bmatrix} = \frac{EI_{y}}{l_{0}^{3}} \left [ \textstyle\begin{array}{r@{\quad }r} 4 & -2 \\ -2 & 4 \end{array}\displaystyle \right ]\,,\quad \begin{bmatrix} S^{\,s}_{55} & S^{\,s}_{56} \\ S^{\,s}_{65} & S^{\,s}_{66} \end{bmatrix} = \frac{EI_{z}}{l_{0}^{3}} \left [ \textstyle\begin{array}{r@{\quad }r} 4 & -2 \\ -2 & 4 \end{array}\displaystyle \right ]\,, $$

$$\begin{aligned} \begin{bmatrix} S^{\,s}_{22}& S^{\,s}_{27}& S^{\,s}_{28} \\ & S^{\,s}_{77} & S^{\,s}_{78} \\ \mathrm{symm} & & S^{\,s}_{88} \end{bmatrix} &= \frac{GI_{\mathrm{t}}}{30l_{0}^{3}} \begin{bmatrix} 36 & -3 & -3 \\ & 4& -1 \\ \mathrm{symm} & & 4 \end{bmatrix} + \frac{EI_{\omega }}{l_{0}^{5}} \begin{bmatrix} 12 & -6 & -6 \\ &4 & 2 \\ \mathrm{symm}& & 4 \end{bmatrix}\,, \end{aligned}$$

(89)

where

$$\begin{aligned} \begin{bmatrix} G^{\,s}_{12} & G^{\,s}_{17} & G^{\,s}_{18} \end{bmatrix} &= -\frac{EI_{p}}{30l_{0}^{4}} \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} 36\varepsilon _{2}-3\varepsilon _{7}-3\varepsilon _{8} & -3 \varepsilon _{2}+4\varepsilon _{7}-\varepsilon _{8} & -3 \varepsilon _{2}-\varepsilon _{7}+4\varepsilon _{8} \end{array}\displaystyle \right ], \end{aligned}$$

$$\begin{aligned} & \begin{bmatrix} G^{\,s}_{22} & G^{\,s}_{27} & G^{\,s}_{28} \\ & G^{\,s}_{77} & G^{\,s}_{78} \\ \mathrm{symm}& & G^{\,s}_{88} \end{bmatrix} \\ & \quad =+\frac{\beta _{z}EI_{y}}{15l_{0}^{5}} \begin{bmatrix} -36(\varepsilon _{7}-\varepsilon _{8}) &~ -6(3\varepsilon _{2}+2 \varepsilon _{7}) &~ 6(3\varepsilon _{2}+2\varepsilon _{8}) \\ &~ \phantom{-} 0 &~ 0 \\ \mathrm{symm} &~ &~ 0 \end{bmatrix} \\ &\qquad +\frac{\beta _{y}EI_{z}}{15l_{0}^{5}} \begin{bmatrix} \phantom{-} 36(\varepsilon _{7}-\varepsilon _{8}) &~ 6(3\varepsilon _{2}+2 \varepsilon _{7}) &~ -6(3\varepsilon _{2}+\varepsilon _{8}) \\ &~ 0 &~ \phantom{-} 0 \\ \mathrm{symm} &~ &~ \phantom{-} 0 \end{bmatrix} \\ &\qquad +\,\frac{\beta _{\omega }EI_{\omega }}{15l_{0}^{6}} \begin{bmatrix} -36(\varepsilon _{7}-\varepsilon _{8}) &~ -6(\varepsilon _{7}-3 \varepsilon _{2}+2\varepsilon _{8}) &~ \phantom{-} 6(2\varepsilon _{7}-3\varepsilon _{2}+\varepsilon _{8}) \\ &~ \phantom{-} 2(6\varepsilon _{2}+10\varepsilon _{7}-\varepsilon _{8}) &~ \varepsilon _{7}-\varepsilon _{8} \\ \mathrm{symm} &~ &~ -2(6\varepsilon _{2}-\varepsilon _{7}+10 \varepsilon _{8}) \end{bmatrix} \\ &\qquad + \frac{E\big(I_{\mathrm{n}}-\textstyle \frac{I_{0}^{2}}{A}+\textstyle \frac{I_{p}^{2}}{A}\big)}{420l_{0}^{7}} \Bigg\{ 3\varepsilon _{2}\\ &\qquad \times \begin{bmatrix} 72(4\varepsilon _{2}-\varepsilon _{7}-\varepsilon _{8}) ~& -12(3 \varepsilon _{2}-2\varepsilon _{7}) ~& -12(3\varepsilon _{2}-2 \varepsilon _{8}) \\ ~& \phantom{-}2(6\varepsilon _{2}+\varepsilon _{7}-\varepsilon _{8}) ~& -2( \varepsilon _{7}+\varepsilon _{8}) \\ \mathrm{symm} ~& ~& \phantom{-}2(6\varepsilon _{2}-\varepsilon _{7}+\varepsilon _{8}) \end{bmatrix} \\ &\qquad + \begin{bmatrix} 36(\varepsilon _{7}^{2}+\varepsilon _{8}^{2})~&3\varepsilon _{7}^{2}-6 \varepsilon _{7}\varepsilon _{8}-3\varepsilon _{8}^{2}) ~&-3 \varepsilon _{7}^{2}-6\varepsilon _{7}\varepsilon _{8}+3\varepsilon _{8}^{2}) \\ ~& 2(12\varepsilon _{7}^{2}-3\varepsilon _{7}\varepsilon _{8}+ \varepsilon _{8}^{2})~&-3\varepsilon _{7}^{2}+4\varepsilon _{7} \varepsilon _{8}-3\varepsilon _{8}^{2} \\ \mathrm{symm} ~& ~&\phantom{-}2(\varepsilon _{7}^{2}-3\varepsilon _{7}\varepsilon _{8}+12 \varepsilon _{8}^{2}) \end{bmatrix}\Bigg\} \,, \end{aligned}$$

$$\begin{aligned} \begin{bmatrix} G^{\,s}_{23} & G^{\,s}_{24} \end{bmatrix} &= \phantom{-} \frac{\beta _{z}EI_{y}}{15l_{0}^{5}} \left [ \textstyle\begin{array}{c@{\quad }c} \!\! \phantom{-} 6(6\varepsilon _{2}+\varepsilon _{7}-2\varepsilon _{8}) &~ 6(6 \varepsilon _{2}-2\varepsilon _{7}+\varepsilon _{8})\!\! \\ \end{array}\displaystyle \right ]\,, \\ \begin{bmatrix} G^{\,s}_{25} & G^{\,s}_{26} \end{bmatrix} &= -\frac{\beta _{y}EI_{z}}{15l_{0}^{5}} \left [ \textstyle\begin{array}{c@{\quad }c} \!\! \phantom{-} 6(6\varepsilon _{2}+\varepsilon _{7}-2\varepsilon _{8}) &~ 6(6 \varepsilon _{2}-2\varepsilon _{7}+\varepsilon _{8})\!\! \\ \end{array}\displaystyle \right ]\,, \\ \begin{bmatrix} G^{\,s}_{37} & G^{\,s}_{38} \\ G^{\,s}_{47} & G^{\,s}_{48} \end{bmatrix} &= \phantom{-} \frac{\beta _{z}EI_{y}}{15l_{0}^{5}} \left [ \textstyle\begin{array}{c@{\quad }c} \!\! \phantom{-} 6\varepsilon _{2}+10\varepsilon _{7}-\varepsilon _{8} &~ -(12 \varepsilon _{2}+\varepsilon _{7}+2\varepsilon _{8})\!\! \\ \!\!-(12\varepsilon _{2}+\varepsilon _{7}+\varepsilon _{8}) &~ \phantom{-} 6\varepsilon _{2}-\varepsilon _{7}+10\varepsilon _{8}\!\! \end{array}\displaystyle \right ]\,, \\ \begin{bmatrix} G^{\,s}_{57} & G^{\,s}_{58} \\ G^{\,s}_{67} & G^{\,s}_{68} \end{bmatrix} &= -\frac{\beta _{y}EI_{z}}{15l_{0}^{5}} \left [ \textstyle\begin{array}{c@{\quad }c} 6\varepsilon _{2}+10\varepsilon _{7}-\varepsilon _{8} &~ -(12 \varepsilon _{2}+\varepsilon _{7}+2\varepsilon _{8}) \\ -(12\varepsilon _{2}+\varepsilon _{7}+\varepsilon _{8}) & 6\varepsilon _{2}-\varepsilon _{7}+10\varepsilon _{8} \end{array}\displaystyle \right ]\,. \end{aligned}$$