Engineering Structures 265 (2022) 114446

Contents lists available at ScienceDirect

Engineering Structures
journal homepage: www.elsevier.com/locate/engstruct

Nonlinear instability problem for geometrically exact beam under

conservative and non-conservative loads
Rosa Adela Mejia-Nava a, *, Ismar Imamovic c, 1, Emina Hajdo c, 2, Adnan Ibrahimbegovic a, b, 3
a
Université de Technologie Compiègne-Alliance Sorbonne Université, France
b
Institut Universitaire de France
c
Faculty of Civil Engineering, University Sarajevo, Bosnia and Herzegovina

A R T I C L E I N F O A B S T R A C T

Keywords: In this work we propose a numerical solution procedure for nonlinear instability problems for geometrically
Nonlinear instability problems exact beams (representing arbitrary large displacements and rotations) under either conservative or non-
Non-conservative load conservative load. Both static and dynamics frameworks are developed, with the latter needed for non-
geometrically exact Reissner and Kirchhoff
conservative case. The proposed model can also successfully deliver the analytic solution to linearized insta­
beam models
bility problems (with small pre-buckling displacements) for critical loads of Euler (no shear) and Timoshenko
(with shear) cantilever beam, which is used for validation for both conservative and non-conservative loads. We
also validated the proposed geometrically exact beam models against the analytic solutions to nonlinear problem
(with finite rotations) and pointed out their advantage in providing a more sound interpretation of instability for
both conservative and non-conservative loads. Finally, the proposed models performance is illustrated on more
complex problems with corresponding numerical solutions obtained by using both the geometrically exact finite
elements for Reissner beam (with shear) and for the Kirchhoff beam (no shear).

1. Introduction [9,10]).
The main novelty in this work is to have these geometrically exact
The application field of slender, lightweight structures is currently beam model formulations to bear upon the cases where the applied
expanding to many application areas, from structures typical of road and external forces are no longer conservative. In other words, these forces
air transportation capable of reducing the fuel consumption to biome­ do not derive from a potential, so that the total energy potential cannot
chanics and soft-robotics applications used for different soft materials be obtained. Thus, the fundamental postulate of instability (also known
reinforcements. The instability problem characterizing such slender as Lyapunov condition; e.g. [2]) in terms of disproportional increase of
structures, which can lead to very large pre-buckling displacements and total energy due to a small perturbation (e.g. see [14]) no longer applies.
rotations, is often referred to nonlinear. Furthermore, the same ability to Hence, we propose in this work a novel approach based upon dynamic
account for overall large displacements and rotations that come as the equilibrium equations that can be written by means of numerical models
consequence of geometric instabilities, where a small perturbation re­ for structures of arbitrary complexity, here illustrated upon the frame
sults with a disproportionally large response, is also required for the structures.
structural models of this kind. This paper deals with this class of prob­ The problem of instability under non-conservative forces can occur
lems, which can be successfully modelled and solved by using so-called in a number of different application domains, with one of the most
geometrically exact models of beams (e.g. [7,12,15,16,17]). Such a important being fluid–structure interaction (e.g. [6]) where the insta­
model has an ability to represent in the geometrically exact manner the bility problem is referred to as flutter (e.g. [22]). It was a fundamental
large overall motion, including the problems with instabilities (e.g. contribution of Bolotin [19,20], followed by a number of more recent

* Corresponding author. Tel.: +33 622138560.

E-mail addresses: rosa-adela.mejia-nava@utc.fr (R. Adela Mejia-Nava), imamovic.ismar@gmail.com (I. Imamovic), emina.hajdo@gmail.com (E. Hajdo), adnan.
ibrahimbegovic@utc.fr (A. Ibrahimbegovic).
1
Tel.: +387 61776806.
2
Tel.: +387 61578960.
3
Tel.: +33 607949780.

https://doi.org/10.1016/j.engstruct.2022.114446
Received 19 March 2021; Received in revised form 18 May 2022; Accepted 21 May 2022
Available online 8 June 2022
0141-0296/Published by Elsevier Ltd.
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

2. Geometrically exact 3d model for reissner shear deformable

beam

In this section, we first develop the theoretical formulation of a three-

dimensional geometrically exact Reissner beam model [12] that can
account for bending, axial and shear deformations, large displacements
and large rotations. The corresponding numerical model is obtained by
space discretization, carried out by the finite element method, and the
time discretization, implementing the Newmark time-stepping scheme.

2.1. Strong and weak formulations

The beam strain measures are denoted as: k regrouping bending and
Fig. 1. Three-dimensional curved beam initial and deformed configurations.
torsion, and ∊ regrouping axial and shear strains. For geometrically
exact beam, these can be written in vector notation as.
works (e.g. [23,25,26,27,28], among others), to show that the instability ′
∊ = φ − a ; k = ω’ (1)
problems under non-conservative loading have to be treated within the
dynamics framework. Here, one important result obtained with the where a is the local basis unit vector perpendicular to the cross-section
geometrically exact beam model concerns solving the Bolotin paradox in the deformed configuration and φ is the position vector of the beam
[19,20] that the presence of damping reduces the value of critical load. axis in the deformed configuration (see Fig. 1). In (1) above, superposed
The geometrically nonlinear response under non-conservative loads prime denotes the derivative with respect to the beam axis coordinate
cannot in general be computed analytically, with a few exceptions for ‘s’, i.e. (•)′ = d/ds (•). Finally, ω is the axial vector of the skew-
academic problems [18,19,20,30]. The nonlinear dynamics computa­ symmetric tensor.Ω
tions pertaining to instability problems of more complex structures can
Ωb = ω × b; ∀b ∈ R3 (2)
be carried out by using the time-integration schemes; our choice here
being the Newmark implicit scheme (e.g. [5,9,25]). which is computed as the derivative of the orthogonal tensor Λ, i.
We show that the geometrically nonlinear computations can also be e..Ω = ΛẤ ΛT
used to obtain nonlinear response under conservative loads in statics. The orthogonal tensor Λ rotates the cross-section in its new position
This serves to validate our geometrically exact approach against analytic to the deformed configuration, defined by a new position of unit vector
solution for Euler elastica [29], and further illustrate its capabilities to a, with a¼Λg. Hence, the most important advantage of the Reissner
deliver accurate solution under large overall motion. Here, we describe beam model is in separating the tensor representation of 3D finite
and compare two geometrically exact beam models capable of rotation and vector representation of displacement field, which both
describing finite rotations. The one is the Reissner beam, including shear contribute to geometrically exact strain measures in (1).
deformation, and the second is so called Kirchhoff beam, excluding Remark: We want to show that (1) represents objective spatial strain
shear deformation. The computed results with two different models can measures residing in the deformed configuration but parameterized with
serve not only to better understand the instability phenomena under respect to the initial configuration. To that end, we first consider a case
different modelling hypotheses, but also to quantify the relative where the current configuration has zero deformation, when obtained
importance of different model features and its ability to provide either by a (large) rigid body motion represented by a constant rotation
more precise results or more robust computations. The key de­ (orthogonal) tensor R (R’¼0) and a constant translation vector t (t’ ¼
velopments that allow both models to be used for efficient computations 0),
of instability problems under non-conservative loads pertain to corre­
sponding follower force tangent stiffness, which is here presented for φrbm = Rx + t→∊rbm := φ’ − a = R x’ − Rg = Rg − Rg = 0
both Reissner and Kirchhoff beam models.
The outline of the paper is as follows. In Section 2, we present the Λrbm = R→K rbm := R’RT = 0→krbm = 0
theoretical formulation of the Reissner geometrically exact beam model,
Thus, with ∊ and k remaining equal to zero in rigid body motion
including its weak form, the corresponding finite element discrete
confirms that (1) indeed represents strain measures for geometrically
approximation and the special follower force element used for
exact beam. Next, we consider a (large) rigid body motion superposed
geometrically exact treatment of non-conservative loads. In Section 3,
on a given deformed configuration (φ, Λ), with again R a constant
we specialize our beam formulation to so-called Kirchhoff beam where
rotation (orthogonal) tensor (R’¼0) and t a constant translation vector
the beam cross-section remains perpendicular to beam axis, and
(t’ ¼ 0); we thus get a new current configuration with the corresponding
consequently eliminates any contribution the shear deformation. We
stains given as.
note in passing that such is a more general model than the extension to
geometrically nonlinear regime of the classical engineering beam φsrbm = RΛx + t→∊srbm := Rφ’ − Ra = R∊
formulation in terms of co-rotational beam model (e.g. [3]). For further
clarity, the development is also specialized to two-dimensional case, Λsrbm = RΛ→Ksrbm := (RΛ)’ (RΛ)T = RΩ’ RT →ksrbm = Rk
which allows to simplify finite rotation representation and provide the
The last result thus confirms that the strain measures in (1) are
converged analytic solutions to model validation problems for geomet­
indeed objective in the sense that ensures the frame-indifference under
rically exact beam: pure bending with zero shear and axial deformation
superposed rigid body motion [38].
with deformed deformation as (part of) a circle, or instability of Euler
The equilibrium equations for the geometrically exact beam are set
elastica (e.g. see [29]). In Section 4, we present and discuss the results of
in terms of stress resultants regrouping axial and shear forces in n = (N,
several numerical simulations, including both static and dynamic non-
V2, V3)T and regrouping torsion and bending moments in m = (T, M2,
conservative loadings. Finally, in Section 5, we present conclusions of
M3)T, all acting in beam deformed configuration but parameterized with
our studies.
respect to initial configuration coordinate s; which we can write as.

n’ + f = 0; (3a)

2
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

the most efficient rotation update, by replacing the orthogonal tensors

(3b)
′
m’ + φ xn = 0
multiplication, by corresponding computation of resulting quaternion,
The final set of equations that define the connection between strain which can be defined as follows:
measures in (1) and stress resultant in (2) are the constitutive equations.
Λ1 := (q0 , q); Λ2 := (p0 , p), (7)
By making the simplest choice for beam linear elastic behavior, we can
write. Λ = Λ2 Λ1 := (q0 ←p0 q0 − p • q, q←p0 q − p × q) (8)
n = ΛCΛT ∊; C = diag(EA, GAk2 , GAk3 ) (4a) However, it is known [7] that by using the quaternion parameters for
representing orthogonal matrix, the consistent tangent will be non-
m = ΛDΛ k; T
D = diag(GJ, EI2 , EI3 ) (4b) symmetric, which can be repaired as shown here by using the rotation
vector parameterization for finite rotation. Defining the orthogonal
where EA is axial stiffness, with E as Younǵs modulus and A as beam
rotation matrix via rotation vector is possible due to the Euler theorem
cross section, GAki is shear stiffness, with G as the shear modulus and ki
for finite rotation (e.g. [11]) that states there exists a vector θ that re­
as stiffness reduction factor, J is a polar moment and I2 and I3 are mo­
mains invariant to multiplication by the corresponding orthogonal
ments of inertia of the beam cross section. In the simplest choice of linear
tensor Λ. Such vector is referred to as rotation vector and can be written
elastic constitutive equation for shear in (4a) above, following proposal
as.
in [31], we can easily obtain the corresponding value of shear reduction
factor for ki as the ratio of the maximum shear stress computed by θ := Λϑ = Iϑ (9)
continuum mechanics solution (e.g. see [33], p. 220) and average stress
in any particular cross-section. Yet better approach to find the shear where we took into account that both Λ and I are two-point tensors;
reduction factor value ki is by using the model reduction to match the hence, θ is the spatial vector corresponding to a material vector ϑ.
cantilever beam energy computed by continuum mechanics versus the Rotation vector can also be used for matrix representation of the
energy computed by beam theory, which allows accounting for elastic orthogonal rotation tensor Λ leading to the Rodriguez formula, which
beams with non-uniform cross-section (see [37]). Moreover, we can can be written as follows.
adapt such model reduction approach to proposed geometrically exact sinθ 1 − cosθ
model beam to define more complex inelastic constitutive equations in Λ = cosθI + Θ+ θ⨂θ (10)
θ θ2
stress resultants for geometrically exact beam; see our previous works
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
[24,35] for metallic beams and [21,36,39] for reinforced concrete
where θ = θ21 + θ22 + θ23 is the magnitude of the rotation vector θ and Θ
beams. In this work we will stay with the simplest set of constitutive
linear equations. Thus, the strong form of geometrically exact beam is a skew-symmetric tensor for which θ is the axial vector, i.e.
model is simply obtained by replacing kinematics in (1) and constitutive ⎡
0 − θ3 θ2
⎤ ⎛ ⎞
θ1
equation in (4) into equilibrium equation in (3), and imposing the 3
Θb = θ × b ∀ b ∈ R ; Θ = ⎣ θ3 0 − θ1 , θ = ⎝ θ2 ⎠
⎦
appropriate boundary conditions. − θ2 θ1 0 θ3
By considering, for simplicity, the homogeneous Dirichlet boundary (11)
conditions, we can obtain the displacement-type weak form of the [ ] 2
geometrically exact beam, or the virtual work equations, which can be By using the standard vector identity Θ(Θb) = θ⨂θ − θ I b, where b
written as. is an arbitrary constant vector, we recover from (10) above an alterna­
∫{ ∫ } ∫ tive form of the Rodrigues formula, which is proved to be a closed form
Gint − Gext = 0; Gint (⋅) := δ ⋅n + δk⋅m ds, Gext = δu⋅fds (5) solution of the exponential mapping (e.g. [11]).
L L sinθ 1 − cosθ
Λ := exp[Θ] = I + Θ+ ΘΘ (12)
Here, δ∊ and δk are virtual strain measures corresponding to the θ θ2
virtual displacement δu and the virtual rotation δΛ, which ought to be We note in passing the geometry aspects of this problem considering
multiplied respectively by the energy-conjugate stress resultants n and that the rotation of the beam cross-section can be very large, or finite;
couples m in order to define internal virtual work. Of special interest for thus, its intrinsic parameterization is given by an orthogonal tensor Λ,
developments presented herein are computations of the external virtual which is an element of so-called SO(3) group.
work for non-conservative loading, which will be described in Section { }
2.3. The virtual work equations in (5) are highly nonlinear, and have to SO(3) = Λ : R3 →R3 |ΛT Λ = I, detΛ = 1 (13)
be solved by incremental-iterative procedure, which is described next. The principal difficulty introduced by such a parameterization is due
to the fact that SO(3) is not a linear space (it is rather a manifold); hence,
the issues pertinent to the theoretical formulation, consistent lineari­
2.2. Finite rotations
zation and update procedure become more complex, as briefly discussed
here.
The main difficulty in dealing with geometrically exact beam con­
Angular velocities and accelerations are computed by exploiting the
cerns the proper parameterization of the finite rotation tensor Λ. In in­
notion of the rotation vector, as defined in the previous section. With Λ
cremental analysis the update of the finite rotations can be reduced to
the product of incremental orthogonal tensor by the corresponding one as a two-point tensor, in constructing the time derivative Λ̇, one should
computed in the previous time step. Such orthogonal tensor can be choose between its spatial and material representations [9]. For
represented by four quaternion parameters [7], which renders these example, the angular velocity field can be computed in a spatial repre­
updates computationally more efficient. Namely, a set of four parame­ sentation as.
ters includes a scalar part q0 and a vector part q, defining the quaternion Ẇ = Λ̇ΛT (14)
{q0 , q}. Such a four-parameter representation can be used to provide an
explicit form of an orthogonal matrix Λ written as: or in a material representation as.
( )
Λ = 2q20 − 1 I + 2q0 Q + 2q⨂q (6) Ψ̇ = ΛT Λ̇ (15)

with Q indicating skew-symmetric matrix for which vector part q is From (14) and (15) above, and from Λ being a two-point tensor, one
axial vector, i.e. Qb − (q × b) = 0. The quaternion representation offers can readily verify that Ẇ is a spatial, whereas Ψ̇ is a material object. The

3
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

last two expressions can be rewritten as the derivatives of.

⃒ ⃒ ⃒ Λ( ̃ ϑn+1 ) Λn T
̃ θn+1 ) = Λn Λ( (23)
d ⃒ d ⃒ [ ( )] d ⃒
Λ̇ : = ⃒⃒ [Λt ] = ⃒⃒ exp tẆ Λ = Λ ⃒⃒ [exp(tΨ̇) ]
dt t=0 dt t=0 dt t=0 ̃ ϑn+1 )= Λn T Λ(
Λ( ̃ θn+1 ) Λn (24)
(16)
Furthermore, a skew-symmetric tensor and the corresponding
which brings out more clearly the fact that Λ̇ is constructed by su­ orthogonal tensor obtained by its exponentiation, share the same ei­
perposing infinitesimal rotations, either Ẇ in spatial or Ψ̇ in material genvectors [11]. From (23) and (24) above, it follows that.
representation, onto the finite rotations represented by Λ. From the last
Θn+1 = Λn Θn+1 Λn T (25)
expression and orthogonality of Λ, we can establish a mutual relation­
ship between the material and spatial angular velocities skew-symmetric
Θn+1 = Λn T Θn+1 Λn (26)
tensors as.

Ẇ = ΛΨ̇ΛT (17) where Θn+1 b = θn+1 × b and Θn+1 b = ϑn+1 × b, ∀ b ∈ R3 . The last result
in return leads to the following relation between the spatial and material
Ψ̇ = ΛT ẆΛ (18) rotation vectors:

which further leads to the corresponding relationship between their θn+1 = Λn ϑn+1 , ϑn+1 = Λn T θn+1 (27)
axial vectors. It is important to note that θn+1 and ϑn+1 belong to the linear
ẇ = Λψ̇ (19) (tangent) spaces TΛ SO(3) and TI SO(3), respectively.
Within each increment, the final values of the state variables are
ψ̇ = ΛT ẇ (20) established by an iterative update procedure that ensures the satisfac­
tion of the weak form of the equilibrium equations. Details of iterative
where Ẇ b = ẇ × band Ψ̇b = ψ̇ × b , ∀ b ∈ R3 . updates of rotations and displacements are discussed next.
By taking the time derivatives of the expressions in (17) and (18) we At each iteration the incremental displacement update is performed
get, respectively, spatial and material form of the angular acceleration in standard, additive manner.
field.
u(i+1) (i) (i)
n+1 = un+1 + Δun+1 (28)
(18)
T T
Ẅ = Λ̈Λ + Λ̇Λ̇
where Δun+1 , is the (i)th iteration contribution to the incremental
(i)
and.
displacement field, with the superscript ‘(i)’ denoting a typical value of
T
Ψ̈ = ΛT Λ̈ + Λ̇ Λ̇ (19) iteration counter.
The iterative update of finite rotations is more involved in that not
only do we have to choose between spatial and material representations,
2.3. Incremental and iterative updates in discrete approximation but also between different iterative rotation parameters. To elaborate
further upon the latter choice, we first consider the material form of the
We can produce the discrete approximation of the Reissner beam iterative rotation parameters and the kind of the rotation update where
model by using the simplest choice of linear polynomials as shape we make use of the exponential mapping at each iteration; thus, we get.
functions (e.g. see [5]) for displacement field and incremental rotation ( )
(28)
(i+1) (i) ̃ (i)
field, thus reducing the problem to computing the time evolution of Λn+1 = Λn+1 Λ Δψ n+1
nodal values for displacements φ(t) and rotations Λ(t).
The solution procedure is based upon incremental-iterative analysis where Δψ n+1 , is material form of the iterative rotation vector. The same
(i)

where the evolution of configuration space variables is obtained by a

value of the total rotation, Λn+1 , can be obtained by making use of the
(i+1)
step-by-step integration scheme. To that end, the time interval of in­
material form of the incremental rotation vector ϑn+1 , and its iterative
(i)
terest [0, T] is partitioned into a number of time steps: 0 < t1 < t2 ⋅⋅⋅
< tn < tn+1 < Â⋅Â⋅Â⋅ < T . At a typical time tn the values of displace­ increment Δϑn+1 .Since both are the elements of a linear space, they can
(i)

ment and rotation components are denoted as. directly be superposed to get.
( )
dn = φ( tn ); Λn = Λ(tn ) (20) (i+1)
Λn+1 = Λn Λ ̃ ϑ(i) (i)
(29)
n+1 + Δϑn+1

As described further, one-step schemes are used to compute the

From the last two expressions and the incremental rotation vector
evolution of the state variables, so that their values at time tn+1 are
definition in (22), we can obtain.
computed solely on the basis of the corresponding values at time tn . For
( ) ( ) ( )
displacement vector we thus have. ̃ ϑ(i) ̃ (i) ̃ (i) (i)
(39)
Λ n+1 Λ εΔψ n+1 = Λ ϑn+1 + εΔϑn+1
dn+1 = dn + un+1 (21)
which, by the analogy with the procedure described in the previous
where un+1 are incremental displacements. section for angular velocity computation, further leads to.
For rotation update we need to take into account that Λ is a two-point ( )
(i) ̃ T ϑ(i) Δϑ(i)
Δψ n+1 = T (31)
tensor and thus choose between two possibilities corresponding to either n+1 n+1

spatial or material representations. If we denote the spatial incremental

rotation vector θn+1 and the material incremental rotation vector ϑn+1 , where T ̃ is the matrix connecting the rotation vector and rotation matrix
the rotation update can be written as. increments (e.g. see [7]).
By exchanging the order of existing and iterative rotations, the
̃ θn+1 ) Λn = Λn Λ(
Λn+1 = Λ( ̃ ϑn+1 ) (22) rotation update can be performed with the spatial rotation parameters
as.
where Λ(•)
̃ is given by the exponential mapping formula in (10). ( ) ( )
By considering that Λn is an orthogonal tensor from (22) we can (i+1)
Λn+1 = Λ ̃ θ(i) (i)
n+1 + Δθn+1
̃ Δw(i)
Λn = Λ
(i)
n+1 Λn+1 (32)
obtain that.

4
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 2. Nonlinear instability of geometrically exact Euler cantilever beam: initial and deformed configurations for the Reissner model; critical equilibrium state of the
Euler elastica [29]; successive deformed shapes obtained by analytic solution [29].

Contrary to the material rrepresentation of incremental rotation

updates, the mutual relationship between the two spatial representa­ un+1 = hvn + h2 [(0.5 − β)an + βan+1 ] (39)
tions of iterative rotation parameters is no longer analogous to the and.
angular velocity relation. In order to derive the relationship of this kind
we use. vn+1 = vn + h[(1 − γ)an + γan+1 ] (40)

Δw(i) (i) (i) (i) (i)

n+1 = Λn+1 Δψ n+1 ; Δθn+1 = Λn Δϑn+1 (33) 2.5. Follower force finite element for geometrically exact beams
which, along with the relationships in (22) and (31), further leads to.
Most of instability problem solutions are in general available for a
( ) ( ) ( )
(i) (i) ̃ T
Δwn+1 = Λn+1 T
(i) (i) ̃ θ(i)
ϑn+1 Λn T Δθn+1 = Λ ̃ (i) T (i) system under conservative load, such as the dead load, which remains of
n+1 Λn T ϑn+1 Λn Δθn+1
constant direction during deformation and motion of the body. How­
(34) ever, any force acting in the direction perpendicular to the cross section
From the material form in (31) we obtain the corresponding spatial or along beam axis can no longer be considered as conservative. It has
form in (34), along with an additional result. been noticed that the manner of application of the loading has a decisive
( ) ( ) influence upon instability of a system. Structures subjected to non-
̃ − 1 θ(i)
T ̃ − 1 ϑ(i)
n+1 = Λn T n+1 Λn
T
(35) conservative forces can be more sensitive to a loss of stability, which
happens in terms of flutter leading to large oscillations with increasing
( ) ( )
̃− 1 T ̃− 1 amplitudes even when triggered by a small perturbation force. In this
T ϑ(i)
n+1 = Λn T θ(i)
n+1 Λn
section, we study instability problem for non-conservative loads (tied to
which can be proved by taking into account that a tensor T and its the beam cross-section) that cannot be derived from a potential. The
inverse share the same eigenvectors [11], and hence they commute. dynamic framework offers the only appropriate physical model of the
system instability evolution under any applied non-conservative loading
2.4. Dynamics and Newmark implicit time-stepping scheme (e.g. [4,20]).
Here we present more refined development of the enhancement of
By using the rotation vector parameterization of large rotations, we such dynamics framework corresponding to the follower forces. More
can further simplify the construction of time-stepping schemes for dy­ precisely, let us consider the force evolution computed form initial load
namics, contrary to previous works [9]. The standard implementation of value p0 which varies with given a time function g(t), that can express as
the Newmark algorithm [12] can be used to compute the velocities and follows.
accelerations at time tn+1 with. Fn+1 = p0 g( tn+1 ) (41)

vn+1
γ
= un+1 +
β− γ
vn +
(β − 0.5γ)h
an (36) If we consider this load as a follower force tied to a beam cross-
βh β β section, it has to further be multiplied by the orthogonal matrix Λ,
which is increasing in the same way as the body motion under study. In
1 1 0.5 − β
an+1 = un+1 − vn − an (37) this manner, we can obtain the final format of the follower force at time
βh2 βh β
tn+1 :
⎡ ⎤
where h = tn+1 − tn is a typical time step, and β and γ are free Newmark p1
parameters. Replacing these approximations into the weak form of the pn+1 = Λn+1 Fn+1 ; pn+1 = ⎣ p2 ⎦ (42)
balance equations, we obtain a system of non-linear equations in in­ p3
cremental displacements as.
Here, Λn+1 is the corresponding rotation tensor that can be repre­
r(un+1 ) = f n+1 (38) sented with quaternion algebra (see Section 2.1.2), which allows for
truly finite rotations.
Typical choice for β = 14 and γ = 12 leads to the scheme of second-order
The corresponding contribution of such a follower force to the re­
accuracy. The corresponding algorithm is yet referred to as a trapezoidal sidual is defined in terms of the tangent stiffness. To that end, we need to
rule or average acceleration method, the reason for which can be made define the skew symmetric tensor Pn+1 , for which pn+1 is defined as axial
more transparent with an alternative, so-called acceleration form of the
vector, so that for any vector v we can write.
Newmark approximations is used, where we take.

5
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Table 1
Load-deflection data for a geometrically exact beam under compressive load.
0◦ 20 40 60 80 100 120 140 160 176◦
◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦
α
P 1 1.015 1.063 1.152 1.293 1.518 1.884 2.541 4.029 9.116
Pcr
xa 1 0.970 0.881 0.741 0.560 0.349 0.123 -0.107 -0.340 -0.577
L
ya 0 0.220 0.422 0.593 0.719 0.792 0.803 0.750 0.625 0.421
L

Fig. 3. Cantilever beam under concentrated moment.

Fig. 4. Successive deformed shapes under free-end bending moment: discrete approximations with 8 beam elements using either Reissner or Kirchhoff beam model,
starting from initial configuration of straight beam and bending into final deformed configuration of a full circle.

Fig. 5. Free-end displacement components vs. time (load increase parameter) computed by Reissner and Kirchhoff beam models in pure bending.

⎡ ⎤ ⎡ ⎤⎡ ⎤
p2 v3 − p3 v2 0 − p3 p2 v1
Gext (φn+1 , Λn+1 ; δφa ) : = δφa Λn+1 pn+1 (44)
⎣
P(tn+1 )v = pn+1 × v = p3 v1 − p1 v3 ⎦ = ⎣ p3 0 − p1 ⎦⎣ v2 ⎦
p1 v2 − p2 v1 − p2 p1 0 v3
where δφa is virtual displacement, Λn+1 is rotation matrix, pn+1 is the
(43)
current value of the follower force. The follower force will contribute to
The corresponding contribution of this kind of follower force applied the tangent operator according to.
at node ‘a’ to the virtual work principle can then be written as:
DG( • )(Δφa ) = δφa P(tn+1 )Δwa (45)

6
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 6. Successive deformed shapes in the first phase of load increase (top) and
in the second phase for imposed support (rigid body) rotations around hori­
zontal axis when using fixed force.

Fig. 7. Successive deformed shapes in the first phase of load increase (top) and
in the second phase for imposed support (rigid body) rotations around hori­
zontal axis when using follower force.

which should be taken into account writing a modified form of the

linearized form of the equilibrium equation:
[ ]
Lin Gint − Gext = Gint − Gext + Kn+1 un+1 + Kffn+1 un+1 (46)

In summary, we can see that the external follower loading also

contributes to the element tangent stiffness matrix Kffn+1 , for a node
where the follower force is applied. The contribution of the follower
external force to tangent stiffness is as follows:
⎡ ⎤
0 − p3 p2
ff
Kn+1 = Pn+1 := ⎣ p3 0 − p1 ⎦ (47)
− p2 p1 0 Fig. 8. Displacements components for fixed force vs. for follower force during
subsequent support rotations around x-axis.

3. Geometrically exact 2d beam reduced models of reissner and

kirchhoff
Table 2
Euler beam critical load for linearized instability for different FE mesh.
3.1. Strong form of 2D Reissner geometrically exact beam and validation
No. elem 2 10 20 100 200 1000
problems solution
Pcr [N] 1.64709 1.48654 1.48196 1.48050 1.48045 1.48044
For clarity, in this section we further study only 2D case with initially
straight geometrically exact beam aligned with x-axis (see Fig. 2), where
deformed configuration can be defined in a simple manner. More pre­
cisely, the place of each particle at beam axis in the deformed

7
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 9. Results of computation by arc-length method for different FE mesh.

Table 3
Critical load computed by Kirchhoff and Reissner beams for different FE mesh.
Beam type No. elem 2 10 20 100 200 2000

Kirchhoff Pcr [N] 1.5581 1.4837 1.4814 1.4805 1.4804 1.4804

Reissner Pcr [N] 1.6471 1.4866 1.4821 1.4807 1.4806 1.4801

Fig. 10. Successive deformed shapes computed with geometric exact beam models.

configuration can be described by its position vector with only two orthogonal to cross-section in initial configuration (aligned with unit
(large) displacement components u and v to write: vector ei) will be rotated to its new position denoted as (unit) vector ai,
which is performed by multiplication with 2D orthogonal matrix Λ.
φ = (x + u)e1 + ve2
⎛ ⎞
cosθ − sinθ 0 a1 = cosθe1 + sinθe2
Further crucial simplification for 2D case is that (large) rotations of
ai = Λei ; Λ = sinθ cosθ 0 ⎠ ⇔ a2 = − sinθe1 + cosθe2
⎝
beam cross-section no longer need (orthogonal) tensor representation, 0 0 1 a3 = e3
and can simply be defined in terms of rotation angle with respect to x-
axis, which is denoted as θ. With no need for tensor representation, we We can further obtain both axial strain Σ and shear strain Γ by
no longer insist on two-point nature of large rotation, and further use θ projecting the (vector) strain measures (1) to this rotated frame ai.
≅ θ. For 2D Reissner geometrically exact beam, the unit vector

8
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 11. Free-end force vs. angle α diagram computed with geometric exact beam.

( )
du dv
Γ= − 1+ sinθ + cosθ
dx dx
The axial and shear strains Σ and Γ are nonlinear functions, whereas
the curvature or bending strain K is a linear function of rotation angle.
We note that such rotated strains are corresponding material represen­
tation of strain measures in (1), which can also be obtained in more
general 3D case by pull-back (e.g. [38]) of spatial strain to initial
Fig. 12. Cantilever beam under follower load. configuration performed with the orthogonal tensor. Moreover, one can
define the rotated stress resultants corresponding to these strains in
Σ = aT1 φ’ − aT1 a1 ⇒aT1 φ’ = Σ + 1 terms of simple linear elastic constitutive equations, by projecting (4) to
rotated frame in order to obtain.
Γ = aT2 φ’ − aT2 a1 ⇒aT2 φ’ = Γ
N = EAΣ; V = GAΓ; M = EIK
The last result can also be recast as the rotated strain measure for 2D
geometrically exact beam model, as initially proposed by the Reissner Finally, the equilibrium equations in (3) can also be simplified in 2D
[15], which can be written as: case to corresponding scalar equations by projecting onto rotated frame
dθ basis ai to obtain:
K= ( )
dx ′ du dv
N = 0; V ’ = 0; M ’ + 1 + V− N=0
( ) dx dx
du dv
Σ= 1+ cosθ + sinθ − 1 (48)
dx dx We next consider two special cases of equilibrium equations for
geometrically exact beam, where it is possible to obtain corresponding

Fig. 13. Reissner beam free-end displacements in transverse direction.

9
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 14. Shear deformable beam under vertical follower load.

in (3) can be simplified to:

Table 4
Critical force value P** for damped case normalized with respect to critical value N = 0; V = 0; M ’ = 0⇒M = cst.
for non-damped case for different values of damping ratio ξ.
By combining the last two results and taking into account that axial
ξ 0 0.05 0.25 0.5 1 1.5 and shear strains are equal to zero Σ = 0 and Γ = 0, we can obtain well
P** 2.023 2.48 2.5 2.7 3 3.5 posed set of equations for computing beam deformed shape, with
boundary conditions chosen to impose zero rotation and displacement
components at built-in end of the beam.
analytic solution used for validation of our proposed models.
Pure bending: First we consider the pure bending of cantilever in dθ M
= ; θ(0) = 0
Fig. 2 under free-end moment (with zero axial and shear force compo­ dx EI
nents), for which the corresponding projection of equilibrium equations

10
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 15. Normalized critical load for geometrically exact beam for damped case
beam under vertical follower load for different values of damping ratio ξ.

Fig. 16. Frame under two vertical compressive loads, and two horizontal
disturbances.

( )
du
1+ = cosθ; u(0) = 0
dx

dv
= sinθ; v(0)
dx
It easy to see that curvature remain constant and that there is no
change of beam length in deformed configuration, and thus conclude
that the deformed shape is always a part of a circle with same circum­
ferential length as the initial length of the beam. We can integrate these
equations and easily compute the free-end rotation and displacement
components of geometrically exact beam in pure bending as:
ML
θ(L) =
EI

L θ(L) θ(L)
u(L) = L − θ(L)
sin cos (49)
2
2 2

L θ(L) θ(L)
v(L) = θ(L) sin sin
2
2 2

Euler instability: Second special case we consider corresponds to

Euler instability problem, where no restriction on small pre-buckling
displacements is imposed (see Fig. 2 for illustration). Here the corre­ Fig. 17. Critical load of frame under vertical conservative load for Reiss­
sponding projections of equilibrium equations in (4) can be written by ner beam.
consider that we again have statically determined problem. Hence, we
can also obtain the force equilibrium equation by decomposing the
applied force within the same local orthogonal frame in deformed

11
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Table 5
Values of conservative critical load semi-analytic vs. numerical solution.
Semi-Analytic Numerical Solution Numerical Solution
Solution Reissner beam Kirchhoff beam

Case 1.EIb = 0.75 0.74

Pcr =
1*EIc
2.712
*EI = 0.73
l2
Case 2.EIb = 0.3 0.26
Pcr =
0.01*EIc
1.572
*EI = 0.25
l2
Case 3.EIb = 1.0 1.0
Pcr =
100*EIc
3.142
*EI = 0.99
l2

configuration.
N = − λcosθ
n = Na1 + Va2 ⇒
V = λsinθ
thus defining normal force component N and shear force component
V. Furthermore, by replacing the last two results into the moment
equilibrium equation and by using the linear elastic constitutive equa­
tion for bending moment with M = EIθ’, we can write:

0 = m’ + φ’ x n|⇒0 = (EI)θ’’ + λ[(1 + Σ)sinθ + Γcosθ ]

For Euler’s beam, yet called Euler’s elastica, we assume that both
shear and axial strain components are equal to zero, which leads to
equilibrium (nonlinear) equation in deformed configuration defined in
terms of angle θ:

Γ = 0 and Σ = 0⇒EIθ’’ + λsinθ = 0

The corresponding linearized form or Euler instability, referred to as
Euler buckling, can be recovered from the last result by enforcing a small
pre-buckling displacement with.

sinθ ≈ θ⇒EIθ′′ + λθ = 0
By further using that small pre-buckling displacements hypothesis
with θ = dv/dx and λ = P, we recover the classical result for the strong
form of Euler linearized instability problem (e.g. [29]) leading to critical
force value Pcr = π 2 EI/(2 l)2. One can also obtain analytic solution to
nonlinear instability problem with a complete elliptic integral of the first
kind (see [29]), with values that are tabulated in Table 1 for different
values of free-end rotation (see Fig. 2), showing clearly that geometri­
cally exact beam can carry the load (much) superior to critical value Pcr.
for linearized instability.

3.2. Variational formulation for 2D geometrically exact Kirchhoff beam

model

In this section we specialize the previous developments for the case

of so-called Kirchhoff beam where no shear deformation is taken into
account (see [34]). Again, we study only 2D case, where the rotation
matrix can be defied by one parameter representing the angle between
cross section and x-axis, further denoted as θ. The rotated strain mea­
sures for such 2D geometrically exact beam model will further be
modified with respect to those proposed by Reissner [15] and defined in
(48). More precisely, by imposing the Kirchhoff constraint we enforce
that the beam section in the deformed configuration remains not only
plane but also perpendicular to the beam axis. This further results with
Fig. 18. Critical load of frame under vertical conservative load for Kirchh­ zero value of the shear strain (Γ = 0), and allows to recast the result in
off beam.
equation (48)3 as follows:
⎛ ⎞
dv dv
⎜ ⎟
tanθ = dx
⇒̃
θ = arctan⎝ dx du⎠
du
(50)
1+ dx
1 + dx

12
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 19. Critical load of frame under vertical non-conservative load for Reiss­
ner beam. Fig. 20. Critical load of frame under vertical non-conservative load for
Kirchhoff Beam.
Furthermore, the Kirchhoff constraint is also enforcing that the axial
deformation and the orthogonal rotation tensor Λ can be rewritten as a
function of the beam axis displacement components:

13
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Table 6
Values of non-conservative critical load semi-analytic vs. numerical solution.
Semi-Analytic Numerical Solution Numerical Solution
Solution Reissner beam Kirchhoff beam

Case 1.EIb = 1.54 1.55

1*EIc
Case 2.EIb = 2.02 2.05 2.1
0.01*EIc
Case 3.EIb = 1.11 1.25
100*EIc

Fig. 21. frame 3D case 1 of load.

Fig. 23. Critical load case 1 frame 3D For case 2, the displacements are
Fig. 22. deformed configuration case 1. different at both corners, because the orientations of small disturbances being
opposite to one another, as seen in Fig. 25.

( )
du dv
Σ= 1+ cos̃
θ + siñ
θ− 1 (51)
dx dx

14
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

dδu ̃ dδv ̃
δΣ = cosθ − sinθ (53)
dx dx
( ( 2 ))
d2 δu ̃ d2 δv ̃ d δu ̃ d2 δv ̃ 1
δK = − sinθ + cosθ − δ̃θ cosθ + sinθ (54)
dx2 dx2 dx2 dx2 Δl

We note that the virtual rotation δ̃

θ can be expressed in terms of the
first derivatives of virtual displacement field:
( )
dδu ̃ dδv ̃ 1
δ̃
θ= − sinθ + cosθ (55)
dx dx Δl
Finally, with this result in hand, we can express virtual curvature as:
( )
d2 δu ̃ d2 δv ̃ 1
δK = − 2
sinθ + 2 cosθ
dx dx Δl
( )( )
dδu ̃ dδv ̃ d2 δu ̃ d2 δv ̃ 1
− − sinθ + cosθ cosθ + sinθ (56)
dx dx dx2 dx2 Δl2

Fig. 24. frame 3D case 2 of load. For the Kirchhoff beam model, we choose the simplest set of linear
constitutive equations for finite rotated strain measures with the Biot
stress resultants N and M:
N = EA Σ; M = EI K (57)
The weak form of the equations of motion is provided by the
d’Alembert principle, which postulates that the snap-shot of motion
taken at time ‘t’ can be described formally with the equilibrium equa­
tions. Unlike the statics problem, these equilibrium equations should
also include the work of inertia forces, which are proportional to mass
and directed opposite to acceleration. In such an approach, the fixed
time ‘t’ corresponds to a particular deformed configuration, and hence
virtual displacement field is independent of time. Weak form of equi­
librium of equation in the material description, see [5], can be written
as.
∫⋅ ∫⋅ ∫⋅
G(a, ̂
a ) := δuT fi dx + δΣNdx + δKMdx − Gext = 0
L L L

where δuT = ( δu δv δθ ) is a vector of displacement components

and rotations, f i is the inertial force vector, which can be determined
from the kinetic energy by using the Lagrange’s equation of motion.
( ) ( )
d ∂Ek ∂Ek
fi = −
dt ∂u̇ ∂u
The kinetic energy Ek is expressed in terms of mass density ρ and
velocity components in deformed configuration, which can be derived
Fig. 25. deformed configuration case 2. from position vector φ related to reference coordinates.
(( ) ( )) ( ) ( )
x+u − sinθ u̇ − θ̇cosθ
⎡ ⎤ φ̇ = dtd φ = dtd +ζ =. +ζ where
v cosθ v̇ − θ̇sinθ
θ − siñ
̃ = ⎣ cos̃
Λ θ0siñ
θcos̃
θ0001⎦ ζ is the co-ordinate along the normal to the beam axis in the reference
configuration. With this result in hand, we can express kinetic energy.
∫ ∫
By exploiting the results in (50) and (48)1 we can obtain the corre­ 1 1
Ek = ρ φ̇T φ̇dAdx = u̇T Mu̇
sponding expression for the curvature of 2D geometrically exact 2 L A 2
Kirchhoff beam. ( )
where u̇T = u̇ v̇ θ̇ is velocity vector components along the co­
( )
d̃
θ d2 u ̃ d2 v ̃ 1 ordinate axes and rotational velocity, and M is constant mass matrix.
K= = − sinθ + cos θ (52)
dx dx2 dx2 Δl ∫
⎡ ⎤
mu,v 0 0 mu,v = ρdA
where. M=⎣ 0 mu,v 0 ⎦; ∫ A
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
( )2 ( )2̅ 0 0 mθ mu,v = ρζ2 dA
du dv A
Δl = 1+ +
dx dx By using previously defined kinetic energy, inertial forces cane be
expressed as.
The virtual strain measure, featuring in the weak form of equilibrium
( ) ( )
equation, can then be derived by taking the directional derivative of the d ∂Ek ∂Ek d
fk = − = (Mu̇) − 0 = Mü
strain measures in (51) and (52) in the direction of virtual displacements dt ∂u̇ ∂u dt
δu and δv, which can be written explicitly as:

15
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 26. Diagram time versus displacement direction 1, 2 and 3 in node 11 and 21.

16
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 27. Tower under a non-conservative load plus a small disturbance.

The final step needed for numerical implementation is the lineari­

zation of the weak form of equilibrium equations so that an iterative
strategy can be employed. It can be obtained by the consistent lineari­
zation of the expression (67) to get.
⃒
d ⃒
L(G)|a = G(a, ̂ ̂ + βΔa) ] ⃒⃒
a ) |a + [G( a,a
dβ β

⎛ ⎞
∫ ( ) d(Δa)
dd(̂
a) [ ]⎜ ⎟
= G(a, ̂
a) + d(̂
a)
Σ K
D + D ⎝ dd(Δa) ⎠dx (69)
L dx
dx
where DΣ and DK are tangent stiffness’s related to axial and bending
response. Further details on linearization of shear deformable 2D
geometrically exact beam are given in [13], and the modification for
Kirchhoff constraint in [34].

4. Numerical examples

In this section we present several numerical examples in order to

Fig. 28. Tower under follower load, diagram displacement vs. time computed
illustrate the performance of two proposed beam models of Reissner and
for the Kirchhoff beam. of Kirchhoff. The examples consider the instability problems under the
follower load and, when adequate, the comparisons against fixed load,
for both static and dynamic case.

4.1. Static analysis of pure bending under fixed moment

The first one is a simple validation example considering the pure

bending of a cantilever beam under concentrated moment applied at its
free end. Both proposed Reissner beam and Kirchhoff beam models can
deliver the exact solution, due to their ability to properly represent finite
rotation (this is not the case of many other beam models). The selected
mechanical and geometric properties chosen in both models are: length
100, axial stiffness EA = 106 , bending stiffness EI = 1000. For Reissner
beam we also choose the model shear stiffness GA = 106 and torsional
stiffness GJ = 2000, although this choice is not important given that the
corresponding strain components remain equal to zero. Below (see Fig. 3
and Fig. 4) we present a comparison of our numerical results obtained by
using these two models for beam deformed shapes and displacements at
the free end. Both of these results match perfectly with analytic solution
presented in Section 3, with deformed shapes represented as roll-up of
the cantilever beam under correctly chosen value of bending moment
equal to M = 2 π EI/L.
Fig. 29. Tower under follower load, diagram displacement vs. time computed The deformed configuration that is obtained in Fig. 4 remains exactly
for the Reissner beam.
the same for both beam models, and so are the free-end displacement

17
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

Fig. 30. Linear buckling analysis of the tower.

Fig. 31. Tower under conservative plus non-conservative loads, diagram

displacement vs. time computed for the Kirchhoff beam. Fig. 32. Tower under conservative plus non-conservative loads, diagram
displacement vs. time computed for the Reissner Beam.

components in Fig. 5. The results obtained for the Reissner and the
Kirchhoff beam are identical since there is no shear deformation in pure loaded by a vertical force F = 5 at the other, free-end. This force belongs
bending case. Moreover, the computed value of rotation is exactly the to two different load cases. The first one considers that the force is fixed,
same for both geometrically exact models, being directly proportional to increasing from zero to its final value in five equal time steps. The sec­
applied bending moment. With geometrically exact beam model, which ond load case, applied subsequently to the first one, seeks to represent
is using the finite strain measure in (1), it is easy to reproduce full-circle rigid body rotation (by the imposed support rotation around x-axis)
deformed shape. This requires one step further from the proof which we superposed upon the deformed configuration at the end of the first load
have already given that such strains remain zero under rigid body mo­ case. This kind of test is used to illustrate the discrete model perfor­
tion, to show that the strain remain constant (or zero) in pure bending. mance under superposed rigid body motion (here support rotation),
Thus, the result is exactly reproduced even by the lowest order of which should leave the deformed state – with product between strains
discrete approximation with 2-node geometrically exact beam elements. and stress resultants as the internal energy of the beam – the same
It is not as simple for many other beam element formulations. throughout entire superposed rigid body motion (just like the strain
measures in (1) remain invariant). However, in order to keep the
external energy the same, the support rotation will also affect the
4.2. Static analysis under fixed-follower force and imposed support
initially vertical position of the fixed force by forcing it to follow the
rotation
beam cross-section and turning it into a follower force. Thus, in repre­
senting the superposed rigid body motion, it is wrong to keep the di­
This example, adapted from [7] considers bending of an L-shape
rection of the force fixed. In such a case the results at the end of any full
frame, with length of 10 for each member. The mesh was constructed
turn are different from those obtained previously, in that there appears a
with four elements in each leg. The material and section properties are
constant shift from one turn to another; see Fig. 6. The main reason for
selected as: EA = GAky = GAkz = 106 and GJ = EIy = EIz = 103 . The
this difference pertains to a nontrivial external work contribution by the
frame is initially placed in a horizontal plane, built-in at one end and

18
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

fixed force during the support rotation phase. Namely, as it is clear that 106 , bending stiffness EI = 1000. For Reissner beam model we also
the rigid body motion of the support should not perturb the existing choose shear stiffness GA = 106 and torsional stiffness GJ = 2000. We
stress state established at the end of the first loading phase, this can be explore different values of applied force and small perturbation in order
achieved only if the force will subsequently follow the support rotation to find the critical load by disproportional vibration response of the
as well. By taking into account that in the second loading phase the cantilever beam, which mimics the Bolotin solution. For the Reissner
cantilever and applied load rotate through imposed rigid body motion, beam model with 8 finite elements, we can see in Fig. 13a that any load
with no change in internal stress state that remains constant under value smaller than 2 is not producing large displacement amplitudes in
support rotation. Our computations carried out for that kind of follower free vibrations triggered by a small perturbation. Hence, the cantilever
force confirm that the deformed shapes are indeed identical in each and beam remains stable. The value of the critical load is equal to 2.2, since
every turn, as shown in Fig. 7. Further illustration of the fundamental then such vibration amplitudes become very large (see Fig. 13). In the
difference in results between two load cases is given Fig. 8, where we case of refined mesh with 100 elements, we can see in Fig. 13b that the
present free-end displacement components at the end of each turn for critical load value is reduced to 2.1. We consider this result to be logical,
100 turns in total. since such a value is closer to the analytical solution, given that 100
elements mesh can better represent true deformed buckled shape of the
4.3. Cantilever beam instability under fixed force convergence studies beam obtained by analytical solution in terms of sin-function. However,
given that the difference between the numerical results obtained by
This example provides the numerical solution of the cantilever beam these two meshes is not extremely big, we will carry out parametric
instability under a fixed compressive load. This reduces to the classical studies by using only a coarse mesh.
Euler instability problem, if we consider that the beam deformation is We further explore the shear effects in the Reissner beam model by
constrained to have both axial and shear deformations equal to zero. In changing the value of the shear stiffness kGA, with a smaller value of k
such a case, one can obtain the analytic solution [29] by integrating the increasing the contribution of shear deformation. Here, we perform
governing equilibrium equations, to obtain the results presented in parametric studies with different values of parameter k. The marked
Table 1. difference from the case without shear is seen only for rather small
In computations to follow, we have chosen: EA = 106 [N], EI = 0.6 values of k, where the influence of shear deformation becomes more
10 [N cm2], L = 100 [cm]. The corresponding analytic solution for
4
significant (see Fig. 14).
critical force value for linearized instability is Pcr = 1.48044 [N] for We next repeat the same analysis, but in the presence of damping
Euler beam and Pcr = 1.48007 [N] for Timoshenko beam if we take GAk that can influence the free vibration phase. Apparently, according to
= EA. Bolotin [20], the presence of damping produces a paradox-type result
Here, we first carried out the computations for this problem by using that by adding the damping one can reduce the corresponding value of
the geometrically exact Kirchhoff beam element, and then recomputed critical force. A detailed development the analytic solution is presented
the results corresponding to using the Reissner beam element, where we in [20]. Here, we repeated the computations above by using the same
can illustrate the influence of shear and eventual difference from the model with an addition of viscous dampers placed at each and every
corresponding critical load reduction due to shear for linearized insta­ degree of freedom to produce the mass matrix proportional part of
bility case (see [4]). In order to better illustrate the influence of discrete Rayleigh damping (e.g. [41]), with the corresponding value of damping
approximation, we present in Table 2 the computed value of critical load ratio obtained for the first vibration mode of the cantilever beam. We
for linearized instability obtained by approach presented in [4]. have performed the numerical computations to carry out parametric
The computations can be performed by using the arc-length method studies by varying the damping ratio starting from zero, which can
in presence of instabilities (e.g. see [1,8]), which allows computing the recover the reference solution for non-damped case. The results of our
value of critical load with higher precision (see Fig. 9) and Table 3), computations with geometrically exact beam mode are presented in
which can be compared against the linearized instability computational Table 4 for critical force value for damped case, starting with the critical
results in Table 2 above. An alternative computation that can done by value for non-damped case, which corresponds to Bolotin solution with
imposing the vertical displacement of the free-end of the cantilever P*=2.023. In Fig. 15, we illustrated the computation results for critical
beam (to get the force value as corresponding reaction) is more robust, load P** obtained with geometrically exact beam for damped case
but the can lead to insufficiently precise identification of the critical normalized by the corresponding result for non-damped case presented
force. in [20] for linearized instability case. We can see that contrary linearized
The computed deformed shapes are plotted in Fig. 10a and 10b for theory, which predicts the reduction of critical load for increasing values
the Kirchhoff and the Reissner beam models, respectively, and the cor­ of damping (see [20], p. 99), the geometrically exact beam model leads
responding force values vs. free-end cross-section rotation angle values to increasing value of critical load for higher damping values. The same
are plotted in Fig. 11a and 11b. The results are normalized with respect conclusion is also made in [23], based upon simplified model with two-
to theoretical values of critical loads for linearized instability with shear bar system studied in [30].
for the Reissner beam and without shear for the Kirchhoff beam model. Fig. 15 Normalized critical load for geometrically exact beam for
damped case beam under vertical follower load for different values of
4.4. Cantilever beam instability under follower force convergence studies damping ratio ξ.

This example provides the numerical solution of the cantilever beam

instability under a follower compressive load (see Fig. 12). Such insta­ 4.5. Frame structure instability-convergence studies
bility problem has to be placed within dynamics framework where
instability is triggered by a small disturbance. For the simplest case of In this example we propose a frame under two compressive loads
linearized instability we can use the analytical solution proposed by applied in vertical direction, along with two small disturbances in hor­
Bolotin [19,20], and its generalization to shear deformable beam given izontal direction, as shown in Fig. 16. We are searching for the value of
in our recent work on linearized instability with small pre-buckling the critical load by following the same procedure as in the previous
displacements (Hajdo et al. [4]). Here, we test this problem with example. We study three different cases, first for conservative loads and
Reissner and Kirchhoff beam. In order to find the critical load for the first then for non-conservative loads. We vary the values of bending stiffness
model, we use two different meshes. The first mesh is constructed with as follows: i) all elements have the same bending stiffness EI = 1000, ii)
eight elements and the second with 100 elements. The chosen values for columns have the same bending stiffness EI = 1000 and the beam is
material and geometric properties are: length 100, axial stiffness EA = weaker with only 1% of the bending stiffness of columns, iii) columns

19
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

have the same bending stiffness EI = 1000 and beam is stronger with 100 corresponds to a real structure of this size (see Fig. 27a). It was first
times bigger bending stiffness than columns. In each case we obtain the designed to withstand standard mechanical loads, including the self-
value of the critical load. The length of each column is equal to 10, the weight, and the fixed load. Here, we check the instability of such a
same as for the beam. The axial stiffness for all the cases remains the structure under an external non-conservative load, applied as a follower
same and equal to EA = 106 . force acting at deformed configuration (see Fig. 27b), which causes
First, we study the frame under conservative load. In order to obtain failure of the top of the tower structure, while in the rest of the structure
the critical value, we try with different values of load in each case, as large displacement (caused by instability) does not occur. We performed
shown in Fig. 17 for the Reissner beam and in Fig. 18 for the Kirchhoff two different analyses by using two different finite elements beam
beam. models presented herein, the one by Reissner and the other by Kirchhoff.
We present Table 5 with summary of the values computed for critical We also performed analysis for two different load histories. The first
load of the frame under conservative force. In the same table we also loading scenario is similar to cantilever beam in example 4.4 to establish
show the values of semi-analytic solution for each case, which was the geometric instability under a vertical non-conservative load applied
presented in our recent work on linearized theory with small pre- at the top of the tower, followed by a small disturbance in horizontal
buckling displacements (Hajdo et al. [4]). One can compare the ana­ direction. The second load scenario is closer to real life load application,
lytic solution versus the corresponding results obtained by numerical where the first load applied is the self-weight of the structure, followed
solution, (see also Figs. 17 and 18) as summarized in this table. In such a by a follower compressive force applied on top, and finally a small
comparison we can see that values are very close to analytic results, perturbation load.
which confirms the validity of our computational results. For the first load case under a vertical non-conservative load plus
In Fig. 19 we can see the numerical solution for the critical force of small disturbance (see Fig. 27) we performed several computations for
the frame shown in Fig. 17 under non-conservative load for the Reissner the critical load of the tower. The computed value of the critical load for
beam and in Fig. 20 for Kirchhoff beam. Kirchhoff beam is equal to 42 KN, whereas for the Reissner beam equals
We also present in Table 6 a summary of results obtained for all 44KN (see Figs. 28 and 29).
numerical solution and simulations compared with the semi-analytical In order to illustrate an important difference between present
solution obtained in Hajdo et al. [4]. nonlinear instability approach and linearized buckling that we studied
previously Hajdo [4], we recomputed the critical load for this structure
by using linearized buckling approach. This is carried out by setting up
4.6. 3D frame structure instability
the linear eigenvalue problem, with material part of the stiffness matrix
and geometric part of the stiffness matrix, and solving for critical load
In order to also illustrate the behaviour of the Reissner beam in 3D
and corresponding critical mode (see [4] for details). The obtained value
case, we next present in this example an out-of-plane instability for the
of the critical load is 122 kN, which is much higher than the values of 42
same planar frame used in previous example, but now loaded by out-of-
kN or 44 kN obtained with present approach of nonlinear instability. The
plane forces. The frame is built with two columns and one beam, each
difference of this kind can be explained by comparing the corresponding
with length equal to 10. The columns are fixed at the free end. The
critical modes for linearized buckling (Fig. 30b) and the one for
chosen material and geometric properties are: length 100, axial stiffness
nonlinear instability (Fig. 27b), with the former spreading globally
EA = 106 , bending stiffness EI = 1000. For Reissner model we also
throughout the structures whereas the latter much more local. The main
choose shear stiffness GA = 106 and torsional stiffness GJ = 2000. We
concern is not only the important difference between two values, but
study two load cases, the first is compressive loads at both corners
also the fact that (easier) computation of linearized buckling is not on
applied as a follower load in direction 2, while in direction 3 we apply a
the safe side, with much overestimating the critical load value. Hence,
small follower disturbance equal to 0.10 of the compressive loads (see
although more laborious, the nonlinear instability computations clearly
Fig. 21). The second case is the compressive loads at both corners with a
find its use in this case to provide better estimate of true critical load.
follower load in direction 2, plus in direction 3 a small follower
Going one more step further in computing the most reliable estimate
disturbance equal to 0.10 of the main load, the second corner is loaded
of the critical load, we study the second load scenario that is closer to
with a follower load in direction 2, plus in direction 3 a small follower
different extreme loads brought by man-made accidents (airplane crash,
disturbance equal to 0.10 of the main load in the opposite direction to
studied in [40]) or stormy wind effects. Here, the dynamic analysis of
that of the first corner Fig. 22.
the tower gathers the action of three different loads, including the dead-
For the first case of frame in Fig. 22, the displacements are the same
load, the follower force and a small perturbation. For its own weight
in both corners. Thus, we present displacement graph versus time in
(conservative load), this tower can be loaded to carry 10KN at the top of
direction one, two and tree only at one corner. With the purpose to
tower (the value according to Eurocode), which is subsequently kept
obtain the critical load we test with different load values until relatively
constant in time. Second, we apply a follower load (non-conservative
large displacements with respect to the previous ones are observed. In
load) followed by a small disturbance. The computed value for the
Fig. 23 we can see that the critical load for case 1 is equal to 1.8.
critical load for the Kirchhoff beam is equal to 32 KN (see Fig. 31) and for
In order to determine the value of the critical load we try different
the Reissner beam equals 34KN (see Fig. 32).
values of compressive force, and we observe the behaviour at both
corner nodes in direction 3 (see Fig. 24). The value of the critical load is
5. Conclusions
now equal to 2.5, as we can see in Fig. 26.
In this work, we developed novel numerical models capable of
4.7. Illustrative (real-life) example with non-conservative loading solving nonlinear instability problems for frame structures. This implies
that the proposed beam models can handle arbitrary large displacements
In order to illustrate an instability analysis for more complex (closer and rotations that may occur before reaching critical equilibrium state,
to real-life) structures under non-conservative loads and show the po­ and are thus not limited to linearized instability problems studied pre­
tential of the proposed theory, we further present an instability analysis viously [4]. We confirmed this finding by computing the corresponding
of a tower (somewhat simplified Eiffel tower) that can be struck by a numerical solution to Euler instability problem with geometrically exact
follower force. This example is motivated by a real-life inspired attack beam, which delivers high accuracy for not only the bifurcation load
on structure by an airplane, which was studied previously with linear­ under fixed (conservative) force, but for the computed results repre­
ized buckling (small pre-buckling displacements) under combined me­ senting complete nonlinear response with very large displacements and
chanical and thermal loads (e.g. [32]). The geometry of the tower rotation in geometrically exact beam final deformed configuration.

20
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

The most important finding concerns the risk for a significant [4] Hajdo E., Mejia-Nava R.A., Imamovic I., A. Ibrahimbegovic, ‘Linearized instability
analysis of frame structures under nonconservative loads: Static and dynamic
reduction in critical instability load when taking into account the
approach’, Coupled Systems Mechanics, 10, 79-102, (2021).
nonlinear instability phenomena, with large displacements and rota­ [5] Ibrahimbegovic A. ‘Nonlinear solid mechanics. In: theoretical formulations and
tions with respect to linearized instability assuming only small pre- finite element solution methods’. Berlin: Springer; 2009. pp. 1-571,.
buckling displacements. As illustrated in our final example, such dif­ [6] Ibrahimbegovic A. (ed.), ‘Computational Methods for Solids and Fluids: Multiscale
Analysis, Probability Aspects and Model Reduction’, Springer, (ISBN 978-3-319-
ference in linearized and nonlinear instability critical loads can be an 27994-7), pp. 1-493, (2016).
order-of-magnitude large, and furthermore not placed on the safe side in [7] Ibrahimbegovic A, Taylor RL. On the role of frame-invariance of structural
the sense that nonlinear instability predicts lower value of critical load. mechanics models at finite rotations. Comput. Methods Appl. Mech. Eng. 2002;
191:5159–76.
The main reason most likely pertains to significantly different deformed [8] Ibrahimbegović A, Al Mikdad M. Quadratically convergent direct calculation of
shapes of the critical equilibrium state between nonlinear and linearized critical points for 3d structures undergoing finite rotations. Comput. Methods Appl.
instability case. Thus, for any such case one needs the full power of Mech. Eng. 2000;189(1):107–20.
[9] Ibrahimbegović A, Mikdad MA. Finite rotations in dynamics of beams and implicit
geometrically exact beam models in order to capture the true deformed time-stepping schemes. Int. J. Numer. Meth. Eng. 1998;41(5):781–814.
shape for nonlinear instability with significant contribution of pre- [10] Ibrahimbegović A, Shakourzadeh H, Batoz J-L, AI Mikdad M, Guo Y-Q. ‘On the role
buckling deformations. of geometrically exact and second order theories in buckling and post-buckling
analysis of three-dimensional beam structures’. J. Comput. Struct. 1996;61(6):
We have developed two finite element models, the Kirchhoff beam 1101–14.
and the Reissner beam, that can both successfully solve such nonlinear [11] Ibrahimbegovic A, Frey F, Kozar I. ‘Computational Aspects of Vector-like
instability problems. The main difference between these two models is Parameterization of Three-Dimensional Finite Rotations. Int. J. Numer. Meth. Eng.
1995;38:3653–73.
in either including the shear deformation for the Reissner beam or
[12] Ibrahimbegovic A. On FE implementation of geometrically nonlinear reissner’s
excluding it for the Kirchhoff beam. Two beam models only give the beam theory: three-dimensional curved beam elements. Comput. Methods Appl.
same result when no shear contribution is present, such as in the first Mech. Eng. 1995;122:11–26.
validation example of pure bending. Also, with no surprise, the Reissner [13] Ibrahimbegović A, Frey F. Finite element analysis of linear and nonlinear planar
deformations of elastic initially curved beams. Int. J. Numerical Methods Eng.
beam model typically corrects the critical load result provided by the 1993;36(19):3239–58.
Kirchhoff beam model for any case where a significant shear strain [14] Mejia Nava AR, Ibrahimbegovic A, Lozano R. Instability phenomena and their
contribution occurs (although in typical slender beams, such contribu­ control in statics and dynamics: application to deep and shallow truss and frame
structures. Coupled Systems Mechan. 2020;9:47–62.
tion in general remains small.). [15] Reissner E. On finite deformations of space curved beams. J. Appl. Math. Phys.
Another important finding in this paper concerns the computation of 1981;32:734–44.
nonlinear instability problems for non-conservative load, such as fol­ [16] Simo JC, Vu-Quoc L. A three-dimensional finite-strain rod model, Part 11:
computational aspects. Comput. Methods Appl. Mech. Engrg. 1986;58:79–116.
lower force. First, in the validation example, we have shown that the [17] Saje M. Finite element formulation of finite planar deformation of curved elastic
geometrically exact beam models developed herein can deal with large beams. J. Comput. Struct. 1991;39(3-4):327–37.
rotations under follower load that is needed to truly represent the rigid [18] Beck M. Die knicklast des einseitig eingespannten, tangential gedrückten stabes”.
Z Angew Math. Phys. 1952;3(3):225–8.
body motion under support rotation. For full validation of our approach, [19] Bolotin V.V., The Dynamic Stability of Elastic Systems, Holden-Day Inc., (1964).
we have successfully applied the proposed approach to the instability [20] Bolotin VV. Nonconservative problems of theory of elastic stability. Pergamon
problems under non-conservative follower load, which were solved Press; 1963.
[21] Bui NN, Ngo M, Nikolic M, Brancherie D, Ibrahimbegovic A. Enriched timoshenko
previously only for linearized buckling problem in Bolotin [19] by using
beam finite element for modeling bending and shear failure of reinforced concrete
the Euler beam (linearized version of the Kirchhoff beam) and in Hajdo frames. Comput. Struct. 2014;143:9–18.
[4] by using the Timoshenko beam (linearized version of the Reissner [22] Farhat C, Chiu E-y, Amsallem D, Schotté J-S, Ohayon R. Modeling of fuel sloshing
beam). However, we have also successfully applied our models to the and its physical effects on flutter. AIAA Journal 2013;51(9):2252–65.
[23] Jeronen J., R. Kouhia, On the effect of damping on stability of nonconservative
case where the linearized instability solution is not sufficient, with either systems, Proceedings XII Finish Mechanics Days, (eds. R. Kouhia et al.), pp. 77-82,
largely overestimating the critical load, such in the final example, or (2015).
predicting reduction of critical load due to addition of damping. [24] Dujc J, Brank B, Ibrahimbegovic A. Multi-scale computational model for failure
analysis of metal frames that includes softening and local buckling. Comput.
Namely, we have shown that the geometrically exact beam model Methods Appl. Mech. Eng. 2010;199(21-22):1371–85.
removes the Bolotin paradox, which is apparently the consequence of [25] Lacarbonara W. Nonlinear structural mechanics: theory, dynamical phenomena,
using the linearized instability approximation for nonlinear dynamics and modeling. Springer; 2013.
[26] Langthjem MA, Sugiyama Y. Dynamics stability of column subjected to follower
vibrations under critical value of follower (non-conservative) load. load. J. Sound Vib 2000;238(5):809–51.
[27] McHugh KA, Dowell EH. Nonlinear response of an inextensible, free-free beam
subjected to a nonconservative follower force. J. Comput. Nonlinear Dynam. 2020;
Declaration of Competing Interest 15(2):021003.
[28] Sugiyama Y, Langthjem MA, Katayama K. dynamic stability of columns under
nonconservative forces : theory and experiment. Springer; 2019.
The authors declare that they have no known competing financial [29] Timoshenko S, Gere JM. Theory of elastic stability. McGraw Hill; 1961.
interests or personal relationships that could have appeared to influence [30] Ziegler H. The stability criterion in elastomechanics (in German : Die
the work reported in this paper. stabilitatskriterien der elastomechanik). Ing. Arch. 1952;20:49–56.
[31] Timoshenko S.P., G.H. MacCullouoh, Elements of Strength of Materials, van
Nostrand Co., (1949).
Acknowledgment [32] Ibrahimbegovic A, Hajdo E, Dolarevic S. Linear instability or buckling problems for
mechanical and coupled thermomechanical extreme conditions. Coupled Systems
Mechan. 2013;2(4):349–74.
This work was supported by funding from ANR (project MS3C), [33] Hughes TJR. The Finite element method. Linear Static and Dynamic Finite Element
MEAE (project CESPA) and IUF (project MS1479). All this support is Analysis: Prentice-Hall; 1987.
gratefully acknowledged. [34] Imamovic I, Ibrahimbegovic A, Hajdo E. Geometrically exact initially curved
Kirchhoff’s planar elasto-plastic beam. Coupled Systems Mechanics 2019;8:
537–53.
References [35] Imamovic I, Ibrahimbegovic A, Knopf-Lenoir C, Mesic E. Plasticity-damage model
parameters identification for structural connections. Coupled Systems Mechanics
2015;4(4):337–64.
[1] Brank B, Stanic A, Korelc J. On path-following method for structural failure
[36] Jukić M, Brank B, Ibrahimbegović A. Embedded discontinuity finite element
problems. Comput. Mech. 2016;58:281–306.
formulation for failure analysis of planar reinforced concrete beams and frames.
[2] Brogliato B, Lozano R, Landau ID. New relationships between Lyapunov functions
Int. J. Eng. Struct. 2013;50:115–25.
and the passivity theorem. Int. J. of Adaptive Control and Signal Processing 1993;7
[37] Medic S, Dolarevic S, Ibrahimbegovic A. Beam Model Refinement and Model
(5):353–65.
Reduction. Int. J. Eng. Struct. 2013;50:158–69.
[3] Crisfield MA. A consistent co-rotational formulation for non-linear, three-
[38] Marsden JE, Hughes TJR. Mathematical foundations of elasticity. Prentice Hall;
dimensional beam elements. Comput. Methods Appl. Mech. Engrg. 1990;81:
1983.
131–50.

21
R. Adela Mejia-Nava et al. Engineering Structures 265 (2022) 114446

[39] Pham BH, Brancherie D, Davenne L, Ibrahimbegovic A. Stress-resultant models for [40] Ibrahimbegovic A, Herve G, Villon P. Nonlinear impact dynamics and field transfer
ultimate load analysis of reinforced concrete frames and its multi-scale parameter suitable for parametric design studies. Int. J. Eng. Comput. 2009;26:185–204.
estimates. Comput. Mech. 2013;51:347–60. [41] Ibrahimbegovic A, Ademovic N. Nonlinear dynamics of structures under extreme
transient loads. CRC Press; 2019.

22