DOI:10.14311/APP.2021.30.0087
Acta Polytechnica CTU Proceedings 30:87–92, 2021 © Czech Technical University in Prague, 2021

available online at http://ojs.cvut.cz/ojs/index.php/app

BEAM ELEMENT UNDER FINITE ROTATIONS

Emma La Malfa Ribolla∗, Milan Jirásek, Martin Horák

Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7,
166 29 Prague 6, Czech Republic

∗ corresponding author: emma.la.malfa.ribolla@fsv.cvut.cz

Abstract. The present work focuses on the 2-D formulation of a nonlinear beam model for slender
structures that can exhibit large rotations of the cross sections while remaining in the small-strain
regime. Bernoulli-Euler hypothesis that plane sections remain plane and perpendicular to the deformed
beam centerline is combined with a linear elastic stress-strain law.
The formulation is based on the integrated form of equilibrium equations and leads to a set of three
first-order differential equations for the displacements and rotation, which are numerically integrated
using a special version of the shooting method. The element has been implemented into an open-source
finite element code to ease computations involving more complex structures. Numerical examples show
a favorable comparison with standard beam elements formulated in the finite-strain framework and
with analytical solutions.

Keywords: Finite rotations, nonlinear beam, shooting method.

1. Introduction
Highly slender fiber- or rod-like components repre-
sent essential constituents of mechanical systems in
many fields of application such as civil, mechanical
and biomedical engineering. It is widely recognized
that slender bodies can be efficiently modeled apply-
ing a beam theory instead of a 3D continuum me-
chanics theory. Kirchoff proposed the first beam for-
mulation which includes large three-dimensional de-
formations [1, 2], and Reissner completed the theory
for both two-dimensional [3] and three-dimensional
cases [4] with two additional deformation measures
representing the shear distortion of beam segments.

Reissner’s finite-strain beam theory is one of the
most important geometrically nonlinear models, sub-
sequently extended and used by many other authors
for 2-D and 3-D analysis of static as well as dynamic
problems. Simo developed a dynamic formulation for
Reissner’s beam [5] and together with Vu-Quoc [6]
initiated the finite element implementation. He also
introduced the useful concept of a geometrically exact
beam, based on recasting Reissner’s theory in a form
which is valid for any magnitude of displacements and
rotations.

In this paper, a geometrically nonlinear beam
model is formulated for the 2-D case. This model is
applicable when the strains remain in a limited range
while the rotations of beam sections can become ar-
bitrarily large, and it properly accounts for the ef-
fect of curvature on the change of distance between
end sections measured along the chord. Therefore
the model considers large rotations and displacements
while strains are treated as small and simple Hooke’s
law is used (an extension to a nonlinear material law
would be relatively straightforward). Cross sections
are still assumed to remain planar and perpendicular

to the deformed beam axis. The formulation exploits
the equilibrium equations in their strong form. They
are combined with the kinematic equations and gen-
eralized material equations (linking internal forces to
the curvature and centerline extension) and then nu-
merically integrated.

2. Model description
2.1. Basic kinematic assumptions
In this section, the kinematics of the element is pre-
sented, leading to the strain-displacement equations.
Each beam is analyzed in its local coordinate sys-
tem, with the x-axis passing through the end joints
in the undeformed configuration and the z-axis ro-
tated by 90◦ clockwise, while the rotations are posi-
tive counterclockwise (see Figure 1). We use u and w
for horizontal and vertical displacements at an arbi-
trary point, while the displacements at the centerline
are denoted by us and ws.

Let us consider an infinitesimal segment of the
centerline of initial length dx, which is after defor-
mation transformed into a segment of length dx̄ =!

(dx + dus)2 + dw2s . The corresponding centerline
stretch can be expressed as

λs =
dx̄
dx

=
!

(1 + u′s)2 + w′2s (1)

where the prime denotes differentiation with re-
spect to x. The assumption of planarity of the de-
formed section and of its perpendicularity to the de-
formed centerline implies that the rotation angle is
linked to the centerline displacements by the relation

tan ϕ = −
dws

dx + dus
= −

w′s
1 + u′s

(2)

87

http://dx.doi.org/10.14311/APP.2021.30.0087
http://ojs.cvut.cz/ojs/index.php/app


E. La Malfa Ribolla, M. Jirásek, M. Horák Acta Polytechnica CTU Proceedings

Figure 1. Kinematics of deformation of the presented nonlinear beam model.

and that the displacements of a general point with
coordinate z can be expressed as

u = us + z sin ϕ (3)
w = ws − z(1 − cos ϕ) (4)

In contrast to the geometrically linear beam model
with u = us + zϕ and w = ws, here the vertical
displacement varies along the height of the section.
However, the change of distance from the centerline
due to deformation perpendicular to the beam axis is
still neglected.

Differentiating (3)–(4) with respect to x, we obtain

u′ = u′s + z cos ϕ · ϕ
′ (5)

w′ = w′s − z sin ϕ · ϕ
′ (6)

and we can express the longitudinal stretch at an ar-
bitrary location

λ =
!

(1 + u′)2 + w′2 =

=
!

λ2s + 2zλsϕ′ + z2ϕ′2 = λs + zϕ
′ (7)

It is worth noting that the stretch is a linear func-
tion of the spatial coordinate z. Consequently, if we
use a linear stress-strain law formulated in terms of
the Biot strain defined as ε = λ − 1, the resulting
work-conjugate stress will be a linear function of z,
too. The choice of working with linear stress-strain
relation between the Biot strain and stress can be
justified in a regime where the strains are assumed to
remain small. Making use of (7), we can express the
Biot strain as

ε = λ − 1 = λs − 1 + zϕ′ = εs + zκ (8)

where

εs = λs − 1 =
!

(1 + u′s)2 + w′2s − 1 (9)

is the strain at the centerline and

κ = ϕ′ (10)

is the curvature.

2.2. Stress resultants and equilibrium
Based on (8) combined with Hooke’s law σ = E ε,
quantities εs and κ that characterize the deformation
of an infinitesimal beam segment can be linked to
the normal force N and bending moment M by the
standard linear relations

N =
"

A

σ dA =
"

A

E(εs + zκ) dA = EA εs (11)

M =
"

A

zσ dA =
"

A

Ez(εs + zκ) dA = EI κ (12)

in which EA is the sectional axial stiffness and EI is
the sectional bending stiffness. As usual, E denotes
Young’s modulus, A is the sectional area and I is the
sectional moment of inertia.

The structure is assumed to be loaded only at its
joints, i.e., no external loads are applied on individ-
ual beams, and so the internal forces in each beam
can be expressed in terms of the end forces and mo-
ments, based on equilibrium conditions written for
the deformed configuration. From equilibrium of the
finite beam segment between the left end section and
a generic section x, one can, according to Fig. 2, de-
rive expressions

N (x) = −Xab cos ϕ(x) + Zab sin ϕ(x) (13)
M (x) = −Mab + Xabws(x) − Zab(x + us(x)) (14)

in which Xab, and Zab are the left-end forces (positive
to the right and downwards) and Mab is the left-end
moment (positive counterclockwise). For the sake of
brevity, we will consider the moment as one of the
(generalized) forces, and, in a similar spirit, the rota-
tion as one of the (generalized) displacements.

Since the internal forces are involved only alge-
braically, they can be easily eliminated and equations
(2) and (9)–(14) can be converted into the following
set of first-order differential equations for unknown
functions us, ws and ϕ:

ϕ′(x) =
−Mab + Xabws(x) − Zab(x + us(x))

EI
(15)

u′s(x) = −1 + cos ϕ(x)+

+
−Xab cos ϕ(x) + Zab sin ϕ(x)

EA
cos ϕ(x) (16)

88



vol. 30/2021 Beam element under finite rotations

M(x)
N(x)

x

z

Xab

Zab

Mab
x us

ws

!

x̃

z̃

ã

a

Figure 2. Schematic representation of beam internal
forces at a general cross section x and left-end forces
and moment.

w′s(x) = − sin ϕ(x)+

+
Xab cos ϕ(x) − Zab sin ϕ(x)

EA
sin ϕ(x) (17)

These equations are accurate even for very large rota-
tions. The only limitation is that the strain remains
small, which is true if εs ≪ 1 and hκ ≪ 1 where h is
the depth of the cross section.

In the context of a finite element analysis, the joint
displacements are global degrees of freedom that are
determined by equilibrium iterations on the struc-
tural level. In each iteration, one needs to determine,
for each beam, the end forces that correspond to the
given joint displacements and rotations. This means
that constants Xab, Zab and Mab in (15)–(17) are not
known in advance. On the other hand, the displace-
ments at both end sections are prescribed. At the left
end, conditions

ϕ(0) = ϕa (18)
us(0) = ua (19)
ws(0) = wa (20)

can be considered as initial conditions for first-order
differential equations (15)–(17). If the values of Xab,
Zab and Mab are somehow estimated, these equations
can be integrated, at least numerically, and the val-
ues of ϕ(L), us(L) and ws(L) can be determined. The
values of forces Xab, Zab and Mab then need to be ad-
justed such that the resulting kinematic quantities at
the right end satisfy the remaining boundary condi-
tions

ϕ(L) = ϕb (21)
us(L) = ub (22)
ws(L) = wb (23)

For a given set of end displacements, the initial esti-
mate of the left-end forces can be constructed based
on linear beam theory, or on the values at the end of
the previous step if the calculation is done in the con-
text of an incremental iterative structural analysis.

In fact, the suggested approach is a special ver-
sion of the shooting method. Normally, the shooting
method iterates on some of the initial values. Actu-
ally, the left-end forces would appear in static bound-
ary conditions at the left end, which make it possible
to integrate the differential equations of equilibrium
and proceed from external forces to internal forces.
Here, the equilibrium equations were directly pre-
sented in their integrated form (13)–(14), and the
unknown variables play the role of integration con-
stants instead of boundary values.

2.3. Numerical solution
The suggested numerical approach is based on the full
form of equations (15)–(17). The right-end displace-
ments ub, wb and ϕb can be evaluated if the left-end
displacements ua, wa and ϕa and the left-end forces
Xab, Zab and Mab are given. To this end, the deriva-
tives in (15)–(17) can be replaced by finite differences
and the interval [0, L] is divided into N equally sized
numerical segments.

Mapping of the right-end displacements, rotation,
forces and moment based on the left-end displace-
ments and rotation can formally be described by

ub = g(f ab, ua) (24)

where

ua =

#

$
ua
wa
ϕa

%

& , ub =

#

$
ub
wb
ϕb

%

& ,

(25)

f ab =

#

$
Xab
Zab
Mab

%

&

and g is a column matrix of three functions of six
variables. For given values of ua and ub, equation
(24) represents a set of three nonlinear algebraic equa-
tions for unknowns f ab. The solution is found by
the Newton-Raphson method, using the recursive for-
mula

δf
(k)
ab = G

−1
'

f
(k)
ab , ua

( '
ub − g

'
f

(k)
ab , ua

((
(26)

f
(k+1)
ab = f

(k)
ab + δf

(k)
ab , k = 0, 1, 2, . . . (27)

where

G =
∂g

∂f ab
(28)

is the Jacobi matrix of mapping g with respect to
variables f ab. The entries of the Jacobi matrix are
evaluated numerically using the linearized version of
the computational scheme. Once the left-end forces
have converged, the right-end forces are easily eval-
uated from equilibrium conditions written for the
whole beam.

The foregoing equations were formulated in a lo-
cal coordinate system aligned with the undeformed

89



E. La Malfa Ribolla, M. Jirásek, M. Horák Acta Polytechnica CTU Proceedings

M

Figure 3. Pure bending of a cantilever beam sub-
jected to end moment.

beam. Since the deformed shape of the beam and
the distribution of internal forces must be invariant
with respect to rigid-body translations and rotations,
one can also use a rotated coordinate system x̃ − z̃
with the origin located at the left end of the beam
in the deformed configuration (point ã in Fig. 2) and
with the x̃ axis in the direction of the tangent to
the deformed centerline at the left end. With respect
to this coordinate system, the left-end displacements
are always zero, and so ua can be omitted from the
list of arguments of function g. Of course, the right-
end displacements ub must be computed from the
global components of the joint displacements using an
appropriate transformation rule, and the end forces
obtained by the iterative process (26)–(27) must be
transformed back into the global coordinate system.
Consistent linearization of the relation between the
end forces and end displacements leads to the element
tangent stiffness matrix. Details will be presented in
a full journal paper.

3. Numerical examples
A nonlinear beam element based on the proposed ap-
proach has been implemented into OOFEM [7, 8], an
object-oriented finite element code. To verify the im-
plementation and demonstrate the potential of the
suggested approach, two examples involving simple
structures are shown next.

3.1. Pure bending of a cantilever beam
The first test is a cantilever of length L = 8 and bend-
ing stiffness EI = 1 loaded by a concentrated end mo-
ment M on its right-end. The exact solution to this
problem is a circular arc with radius R = EI/M . An
applied end moment, M = π/4, will force the rod to
deform into a full closed circle. In this example, the
moment is applied in six load steps, making the rod
wind around itself at the end of the sixth step. The
deformed shape of the beam at the end of each step
is depicted in Fig. 3. The solution of the presented

Figure 4. Close-up view of the solution at the end
of the sixth step using 6,10,20, and 50 segments.

nonlinear model is compared with the one obtained
by employing the geometrically exact finite beam el-
ement by Simo and Vu-Quoc [6] with a mesh of eight
elements. The exact solution is reported as well. The
overall agreement is good, and the simulation based
on the present model, which uses only one two-noded
element, is closer to the analytical solution.

The example demonstrates that the present model
allows for a dramatic reduction of the number of
global degrees of freedom, but of course the number
of segments N used for numerical integration of the
governing equations (15)–(17) must be chosen high
enough to provide a good approximation. The pre-
sented results have been obtained using 100 segments.
A close-up view of the sixth step circle is showed in
Fig. 5 for calculations in which 6, 10, 20 and 50 nu-
merical segments were employed. To ease the inter-
pretation of the results we plot straight segments in
between the integration points even though the cur-
vature is constant along the beam. It turns out that
20 segments provide an acceptable solution, which is
further improved and gets close to the analytical one
if 50 segments are used. It is also noted that even
with just 10 segments the method provides an accu-
rate estimation of displacements at the integration
points. To better quantify the previous comparisons,
we calculated the sum of the distances between the
the analytical and the computed solutions, normal-
ized with respect to the number of sample points.
The results for the sixth load step are reported in Ta-
ble 1. Already for 8 integration segments, the error
is smaller that 5%. The results indicate that compa-
rable accuracy can be obtained if the present method
replaces traditional finite elements by the same num-
ber of integration segments located within one single
finite element.

3.2. Clamped-hinged circular arch
subject to point load

The second example deals with an arch instability
after large asymmetrical deflections. This specific

90



vol. 30/2021 Beam element under finite rotations

Model Error [%]
Simo and Vu-Quoc [6] 2.98
Present approach, 101 segments 0.23
Present approach, 50 segments 0.29
Present approach, 20 segments 0.70
Present approach, 10 segments 2.13
Present approach, 8 segments 3.17
Present approach, 6 segments 5.39

Table 1. Sum of the absolute values of the differences between the exact and the computed displacement solutions
calculated for the sixth load step. The model of Simo and Vu-Quoc [6] uses an eight-element mesh while the presented
beam model uses one element with a variable number of segments of the integration scheme.

Figure 5. Circular arch geometry.

problem was investigated by many authors [6, 9] and
the exact solution based on the Kirchhoff-Love theory
was given by DaDeppo and Schmidt [10]. The arch
is circular with one boundary clamped and the other
boundary hinged, and is loaded by a point load, as
shown in Fig. 5. The cross section of the beam is a
rectangle of depth h = 2.289 and width t = 1.0. The
elastic modulus is set to E = 1.0 × 106 and the Pois-
son ratio is zero. Geometrical and constitutive pa-
rameters are dimensionless meaning that any length
or force scale can be used for these values.

Simo and Vu-Quoc [6] performed a FE analysis
with a forty-element mesh, leading to a buckling load
of 9.0528 EI/R2, while Wood and Zienkiewicz [9] cal-
culated a buckling load of 9.24 EI/R2. The exact
value of the buckling load reported by DaDeppo and
Schmidt [10] is 8.97 EI/R2. In our simulations, we
used forty elements as well, each with four integra-
tion segments, and we obtained a buckling load of
8.9874 EI/R2, which is very close to the analytical
solution and exhibits a relative error less than 0.02%.
The full comparison with the analytical solution is re-
ported in Fig. 6 where we plot the load-displacement
and the load-rotation curves relative to the top of the
arch (loaded point). Again, the analytical solution is
satisfactorily reproduced. The deformed configura-
tions for the four load levels numbered in the curves
are also shown together with the bending moments
indicated by colors. The calculation was completed
in 612 load steps using the arc-length method with a

1

2

3

4
1

2

3

4
1

2

3

1

2

3

4

Figure 6. Load-displacements and load-rotation of
the top of the arch (top). Grey-dashed line repre-
sents the analytical solution reported by DaDeppo
and Schmidt [10]. Deformed shapes and moment
trend along the arch relative to the four load levels
numbered in the top curves.

maximum number of 3 iterations per step.
It is noted that Simo and Vu-Quoc [6] and Wood

and Zienkiewicz [9] show a longer snap-through be-
havior together with the associated deformed shapes
even though the response is not reasonable because of
non-physical penetration of the left support into the
structure. This aspect was investigated by Simo et
al. [11], who showed that a contact constraint con-
dition on the left support needs to be introduced
to obtain a more realistic solution. For the aim of

91



E. La Malfa Ribolla, M. Jirásek, M. Horák Acta Polytechnica CTU Proceedings

this paper we did not consider contact activation and
the subsequent stiffening effect in the structure even
though this extension could be added.

4. Conclusions
We present an efficient formulation for a geometri-
cally nonlinear beam model that accounts for arbi-
trarily large rotations of the beam sections and the
effect of curvature on the change of distance be-
tween end sections measured along the chord. In the
first stage, the model is restricted to the small-strain
framework and cross sections are still assumed to re-
main planar and perpendicular to the deformed beam
axis. Equilibrium equations in their strong form are
combined with the kinematic equations and gener-
alized material equations and then numerically inte-
grated using a special version of the shooting method.

Two examples show the capabilities of the element
to replicate analytical results in a more accurate way
with respect to existing beam finite elements that
account for geometric nonlinearities. Moreover, the
h-refinement of finite elements, which increases the
number of global unknowns and is thus computation-
ally demanding, can be substituted by employing a
proper number of segments in the integration scheme,
while keeping the number of elements limited. The el-
ement will be extended into a formulation that takes
into account an initial curvature and also accounts
for a possible contact mechanism activation.

Acknowledgements
The authors acknowledge the support of the Czech Sci-
ence Foundation (project No. 19-26143X).

References
[1] G. Kirchhoff. Ueber das Gleichgewicht und die

Bewegung eines unendlich dünnen elastischen Stabes.
Journal für die reine und angewandte Mathematik
1859(56):285–313, 1859.

[2] E. H. Dill. Kirchhoff’s theory of rods. Archive for
History of Exact Sciences 44(1):1–23, 1992.

[3] E. Reissner. On one-dimensional finite-strain beam
theory: the plane problem. Zeitschrift für angewandte
Mathematik und Physik ZAMP 23(5):795–804, 1972.

[4] E. Reissner. On finite deformations of space-curved
beams. Zeitschrift für angewandte Mathematik und
Physik ZAMP 32(6):734–744, 1981.

[5] J. C. Simo. A finite strain beam formulation. The
three-dimensional dynamic problem. I. Computer
Methods in Applied Mechanics and Engineering
49(1):55–70, 1985.

[6] J. C. Simo, L. Vu-Quoc. A three-dimensional
finite-strain rod model. Part II: Computational
aspects. Computer Methods in Applied Mechanics and
Engineering 58(1):79–116, 1986.

[7] B. Patzák, Z. Bittnar. Design of object oriented finite
element code. Advances in Engineering Software
32(10-11):759–767, 2001.

[8] B. Patzák. OOFEM—an object-oriented simulation
tool for advanced modeling of materials and structures.
Acta Polytechnica 52(6), 2012.

[9] R. D. Wood, O. Zienkiewicz. Geometrically nonlinear
finite element analysis of beams, frames, arches and
axisymmetric shells. Computers & Structures
7(6):725–735, 1977.

[10] D. DaDeppo, R. Schmidt. Instability of
clamped-hinged circular arches subjected to a point
load. Journal of Applied Mechanics, Transactions
ASME 42(4):894–896, 1975.

[11] J. Simo, P. Wriggers, K. Schweizerhof, R. Taylor.
Finite deformation post-buckling analysis involving
inelasticity and contact constraints. International
journal for numerical methods in engineering
23(5):779–800, 1986.

92