# A cell-centred finite volume formulation of geometrically-exact Simo-
Reissner beams with arbitrary initial curvatures

Seevani Bali  Željko Tuković  Philip<EMAIL_ADDRESS>Alojz
Ivanković  Vikram Pakrashi  School of Mechanical and Materials Engineering,
University College Dublin, Ireland  Faculty of Mechanical Engineering and
Naval Architecture, University of Zagreb, Croatia  Professor, School of
Mechanical and Materials Engineering, University College Dublin, Ireland  SFI
I-Form Centre, University College Dublin, Ireland  SFI MaREI Centre,
University College Dublin, Ireland  Bekaert University Technology Centre,
School of Mechanical and Materials Engineering, University College Dublin,
Ireland  Dynamical Systems and Risk Laboratory, School of Mechanical and
Materials Engineering, University College Dublin, Ireland  UCD Centre of
Adhesion and Adhesives, University College Dublin, Ireland  UCD Centre for
Mechanics, University College Dublin, Ireland

###### Abstract

This paper presents a novel total Lagrangian cell-centred finite volume
formulation of geometrically exact beams with arbitrary initial curvature
undergoing large displacements and finite rotations. The choice of rotation
parametrisation, the mathematical formulation of the beam kinematics,
conjugate strain measures and the linearisation of the strong form of
governing equations is described. The finite volume based discretisation of
the computational domain and the governing equations for each computational
volume are presented. The discretised integral form of the equilibrium
equations are solved using a block-coupled Newton-Raphson solution procedure.
The efficacy of the proposed methodology is presented by comparing the
simulated numerical results with classic benchmark test cases available in the
literature. The objectivity of strain measures for the current formulation and
mesh convergence studies for both initially straight and curved beam
configurations are also discussed.

###### keywords:

Finite volume method , total Lagrangian , geometrically exact beam , block-
coupled , Newton-Raphson

##  1 INTRODUCTION

The mathematical modelling of nonlinear beams has flourished in the past few
decades with their applicability spanning various fields of engineering. Meier
et al. [31] presented a concise review of the existing beam theories. These
formulations, in particular, and the numerical computations of solid mechanics
problems, in general, are mostly analysed using the finite element (FE)
approach; however, since the $1980$s, the applicability of the simple and
conservative finite volume (FV) methods to solid mechanics problems has become
increasingly popular and is evolving rapidly [3]. Particularly in the purview
of beams, the FV implementation of Euler-Bernoulli and shear-deformable
Timoshenko beam theories for both straight and curved beam configurations, can
be found in [32, 14, 22, 12, 13, 26]. Additionally, these works have been
extended to study the stability and buckling analysis of beam-columns as well
as plates [15, 11, 23, 9, 10, 16, 40]. In contrast to the FE implementation,
Fallah and Ghanbari [13] and Wheel [40] reported that the FV formulation did
not suffer from “shear-locking” for both thin Timoshenko beams and thin
Mindlin-Reissner plate analysis respectively. The geometrically exact Simo-
Reissner beam formulation is the most general nonlinear 3-D beam theory
capable of dealing with finite displacements and rotations. The FV formulation
of geometrically exact beams was first investigated by Tuković et al. [38].
The current work builds on the former developed formulations for quasi-static
shear-deformable geometrically exact beams and provides a lucid comparison
between the FE based classic benchmark cases and the proposed FV methodology.
This is the first article to present a FV formulation for beams subjected to
finite displacements and rotations.

The formulations of a beam in a geometrically exact sense considering 3-D
rotations and its FE approximation, was first proposed by Simo et al. [36,
37]. The involvement of finite rotations, elements of the nonlinear
differentiable manifold $SO(3)$ and the complexity of interpolating large
rotations [1, 20] led to three major techniques of addressing them, viz., (a)
incremental and total rotational vector based parametrisations [28, 21, 25,
29], (b) quaternion based rotation interpolation [18, 42, 17, 27] and, (c) co-
rotational beam formulations [7, 4], to name a few. Romero [34] presents a
concise review of different rotation interpolations used for geometrically
exact beams. In the current work, the finite rotations are parameterised using
rotational vectors and the rotational strain measures are updated using
incremental rotation vectors from the previously converged configuration of
the beam. The loss of objectivity in strain measures for FE approximation of
interpolated rotations was first pointed out be Jelenić et al. [6] which was
followed by strain invariant formulations of the geometrically exact beams
[24, 35, 41, 19]. For the current work, the objectivity of adopted strain
measures are numerically verified using a test case for the current work
(Section 4.1). As a recent contribution, Meier et al. [30] developed
geometrically exact Kirchhoff-Love formulation for slender rod geometries. For
a general review on other research advances in this field, like discrete
Cosserat rod kinematics, implementation of different time-stepping schemes in
beam multi-body dynamics, inclusion of nonlinear constitutive laws and
enhanced kinematics of beams, the interested readers are referred to articles
by Meier et al. [31] and Chadha et al. [5].

The article is structured as follows: Section 2 outlines the mathematical
model in total Lagrangian form. The FV discretisation of the mathematical
model is described in Section 3. In Section 4, the proposed methodology is
evaluated on five complementary benchmark test cases, where predictions are
compared to analytical solutions and existing FE benchmarks.

##  2 MATHEMATICAL MODEL

In this section, the mathematical formulation for a 3-D quasi-static, shear-
deformable geometrically-exact Simo-Reissner beam is summarised. Subsequently,
the governing equilibrium equations of spatial forces and moments, the
constitutive relations and their linearised formulation are described.

###  2.1 Kinematic description

A total Lagrangian formulation is adopted to describe the large deformations
of the beam model. Accordingly, a right-handed fixed, reference Cartesian
frame defined by the orthonormal basis vectors $\bm{e}_{1}=[1,0,0]^{\rm T}$,
$\bm{e}_{2}=[0,1,0]^{\rm T}$ and $\bm{e}_{3}=[0,0,1]^{\rm T}$ is specified
(Fig. 1). In the reference (i.e. material) configuration, the mean line of a
beam (i.e. line of centroids of the beam cross-sections) is straight and
parallel to the basis vector $\bm{e}_{1}$ and the cross-sections of the beam
are orthogonal to the basis vector $\bm{e}_{1}$. The principal axis of inertia
of the beam cross-sections are directed along the basis vectors $\bm{e}_{2}$
and $\bm{e}_{3}$.

An initial (stress-free) configuration of the beam mean line is defined by the
space curve $\bm{r}_{0}(s)$, where $s\in[0,L]$ is the beam-length and $L$ is
the initial length of the beam (Fig. 1). To describe the continuous kinematic
description of the moving mean line curve of the beam, a body-attached right
handed orthonormal base vectors $\bm{g}_{0,1}(s)$, $\bm{g}_{0,2}(s)$ and
$\bm{g}_{0,3}(s)$ are defined by the Frenet-Serret formulae, where the base
vector $\bm{g}_{0,1}(s)$ is directed along the initial mean line so that

$\bm{g}_{0,1}(s)=\bm{r}_{0}^{\prime}(s)$ (1)

and the base vectors $\bm{g}_{0,2}(s)$ and $\bm{g}_{0,3}(s)$ are directed
along the principal axis of inertia of the cross-section at $s$. The prime
$(\cdot)^{\prime}$ operator in Eq. 1 and hence forward wherever used, denotes
a derivative with respect to arc-length parameter $s$, i.e.
$(\cdot)^{\prime}\equiv\frac{\partial(\cdot)}{\partial s}$. For an initially
straight beam, the base vectors $\bm{g}_{0,i}$ and the reference bases
$\bm{e}_{i}$ coincide, and the $\bm{g}_{0,1}(s)\equiv\bm{e}_{1}(s)$.

Figure 1: Beam Kinematics: (a) reference, (b) initial and, (c) deformed
configurations of the beam respectively

The orthonormal basis $\bm{g}_{0,1}(s)$, $\bm{g}_{0,2}(s)$, $\bm{g}_{0,3}(s)$
of the local frame and the orthonormal bases $\bm{e}_{1}$, $\bm{e}_{2}$,
$\bm{e}_{3}$ of the global Cartesian frame are related through a linear
transformation $\bm{\Lambda}_{0}(s)$ as

$\bm{g}_{0,i}(s)=\bm{\Lambda}_{0}(s)\bm{e}_{i},\quad i=1,2,3$ (2)

where $\bm{\Lambda}_{0}(s)\in SO(3)$ is the initial two-point second-order
orthogonal rotational tensor field that defines the orientations of beam
cross-sections w.r.t the reference basis. Hence, the initial configuration of
the beam is fully defined by the position vector $\bm{r}_{0}(s)$ of the beam
mean line and orientation of the cross-section at $s$ via the orthogonal
rotation tensor $\bm{\Lambda}_{0}(s)$.

The deformed configuration of the beam mean line is defined by a space curve
$\bm{r}(s)$ and the orientation of the cross-sections via the moving spatial
bases $\bm{g}_{1}(s)$, $\bm{g}_{2}(s)$, $\bm{g}_{3}(s)$. Contrary to the
initial configuration, the basis vector $\bm{g_{1}(s)}$ need not be directed
along the deformed mean line since the beam model is capable of representing
shear deformations. The base vectors $\bm{g}_{2}(s)$ and $\bm{g}_{3}(s)$,
however, are still directed along the principal axes of the cross-section at
$s$ and the orthonormal bases $\bm{g}_{1}(s)$, $\bm{g}_{2}(s)$,
$\bm{g}_{3}(s)$ are related to $\bm{g}_{0,1}(s)$, $\bm{g}_{0,2}(s)$,
$\bm{g}_{0,3}(s)$ bases by the linear transformation

$\bm{g}_{i}(s)=\bm{\Lambda}(s)\bm{g}_{0,i}(s)=\bm{\Lambda}(s)\bm{\Lambda}_{0}(s)\bm{e}_{i},\quad
i=1,2,3$ (3)

where $\bm{\Lambda}(s)\in SO(3)$ is the relative rotational matrix which
rotates beam cross-sections from the initial to the deformed configuration.
Additionally, the deformed mean line of the beam can also be defined using the
mean line displacement vector $\bm{w}(s)$ and the mean position vector of the
initial configuration by,

$\bm{r}(s)=\bm{r}_{0}(s)+\bm{w}(s)$ (4)

The position vector of the deformed mean line $\bm{r}(s)$ and the orientation
of the orthonormal frame $\bm{\Lambda}(s)$ attached to the cross-section at
$s$, ($\bm{r}(s),\bm{\Lambda}(s)$) fully define the deformed configuration of
the beam. Orientation of the beam cross-sections in the deformed configuration
can also be defined w.r.t to the fixed reference bases by the total rotation
matrix,

$\bm{\Lambda}_{t}(s)=\bm{\Lambda}(s)\bm{\Lambda}_{0}(s)$ (5)

###  2.2 Balance equations and strain measures

Assuming that the deformation from the initial to the deformed configuration
is caused by distributed external forces and torques $\bm{f}$ and $\bm{t}$
respectively, the strong differential form of the balance equations are given
as [36],

$\bm{n}^{\prime}\ +\ \bm{f}=0$ (6) $\bm{m}^{\prime}\ +\
\bm{r}^{\prime}\times\bm{n}\ +\ \bm{t}=0$ (7)

where $\bm{n}$ and $\bm{m}$ are the vectors of spatial internal forces and
moments acting over the cross-section at $s$ and “$\times$” symbol denotes the
cross product between two vectors. The corresponding material counterparts
($\bm{N}$ and $\bm{M}$) are related to $\bm{n}$ and $\bm{m}$ via the pull-back
mapping $\bm{\Lambda}_{t}^{\rm T}$ as $\bm{N}=\bm{\Lambda}_{t}^{\rm T}\bm{n}$
and $\bm{M}=\bm{\Lambda}_{t}^{\rm T}\bm{m}$ [36].

The equivalent strong integral form of the above balance equations over a
length, $L$ can be expressed as,

$\int_{L}\bm{n}^{\prime}\textrm{d}L\ +\ \int_{L}\bm{f}\textrm{d}L=0$ (8)

$\int_{L}\bm{m}^{\prime}\textrm{d}L\ +\
\int_{L}(\bm{r}^{\prime}\times\bm{n})\textrm{d}L\ +\
\int_{L}\bm{t}\textrm{d}L=0$ (9)

The required strain measures follow from a geometrically exact beam theory,
where the relationships between the beam configuration and the strain measures
are consistent with the virtual work principle and the equilibrium equations
at a deformed state regardless of the magnitude of displacements, rotations
and strains [37]. To that end, the equivalent strain-configuration
relationships involve three translational strains and a skew-symmetric tensor
of the three rotational strains [36, 21] (from here on, the argument $(s)$ is
dropped from the terms for clarity):

$\bm{\Gamma}=\bm{\Lambda}_{t}^{\rm T}\bm{r}^{\prime}-\bm{\Lambda}_{t}^{\rm
T}\bm{g}_{1}$ (10) $\hat{\bm{K}}=\bm{\Lambda}_{t}^{\rm
T}\bm{\Lambda}_{t}^{\prime}-\bm{\Lambda}_{0}^{\rm
T}\bm{\Lambda}_{0}^{\prime}=\bm{\Lambda}_{t}^{\rm T}\bm{\Lambda}_{t}^{\prime}$
(11)

where $\bm{\Gamma}$ and $\bm{K}$ are (material) translational and rotational
strain measures respectively. Physically, the strain vector $\bm{\Gamma}$
represents axial tension (first entry) and shear deformation whereas the
(material) rotational strain vector $\bm{K}$ represents torsion (first entry)
and bending deformation of the beam body. These strain measures are the
energy-conjugate pairs to the stress resultants ($\bm{N}$ and $\bm{M}$).

The hat $\hat{(\cdot)}$ operator denotes a skew-symmetric matrix associated
with the corresponding (pseudo)-vector given by the relation,
$\bm{a}\times\bm{h}=\hat{\bm{a}}\bm{h}$, $\forall$ $\bm{h}\in\mathbb{R}^{3}$.
In this work, the relative rotation matrix $\bm{\Lambda}$ is parameterised in
terms of its rotational vector $\bm{\psi}$ as

$\bm{\Lambda}(\bm{\psi})=\exp(\hat{\bm{\psi}})=\mathbf{I}\ +\
\frac{\sin\psi}{\psi}\hat{\bm{\psi}}\ +\
\frac{1-\cos\psi}{\psi^{2}}\hat{\bm{\psi}}\hat{\bm{\psi}}$ (12)

where $\mathbf{I}$ is a $3\times 3$ identity matrix and $\psi$ is the
magnitude of rotation vector $\bm{\psi}$. The
$\hat{\bm{\psi}}\in\mathfrak{so(3)}$ is the skew-symmetric tensor living in
the tangent space of $SO(3)$ at $\bm{\Lambda}$ and its exponentiation yields
the finite rotation $\bm{\Lambda}\in SO(3)$. For a rotation vector based
parametrisation, an alternative expression for evaluating incremental
(material) rotational strain vector $\Delta\bm{K}$, as demonstrated in [20,
21, 24] is given by,

$\Delta\bm{K}=\bm{\Lambda}_{t}^{\rm T}\bm{T}^{\rm T}\bm{\psi}^{\prime}$ (13)

where the tangent operator $\bm{T}\in\mathfrak{so(3)}$ is defined as
(“$\otimes$” denotes the dyadic product between two vectors),

$\bm{T}=\frac{\sin\psi}{\psi}\mathbf{I}\ +\
\frac{1}{\psi^{2}}\left(1-\frac{\sin\psi}{\psi}\right)\bm{\psi}\otimes\bm{\psi}\
+\ \frac{1-\cos\psi}{\psi^{2}}\hat{\bm{\psi}}\ $ (14)

###  2.3 Constitutive relations

The present study is limited to linear hyperelastic materials whose length-
specific stored energy function is given by,

$\tilde{\Pi}_{\text{int}}=\frac{1}{2}\bm{\Gamma}^{\rm T}\bm{C}_{\rm
N}\bm{\Gamma}+\frac{1}{2}\bm{K}^{\rm T}\bm{C}_{\rm M}\bm{K}$ (15)

For this hyperelastic model, the material internal forces $\bm{N}$ and moments
$\bm{M}$ are linearly related to the material strain measures ($\bm{\Gamma}$
and $\bm{K}$) as, $\bm{N}=\bm{C}_{\rm N}\bm{\Gamma}$ and $\bm{M}=\bm{C}_{\rm
M}\bm{K}$ where $\bm{C}_{\rm N}=\rm diag{[EA,GA_{2},GA_{3}]}$ and $\bm{C}_{\rm
M}=\rm diag{[GJ,EI_{2},EI_{3}]}$ are constant diagonal constitutive matrices.
Here, $E$ and $G$ denote the elastic and shear moduli of the material
respectively, and $A_{\rm p}$ and $I_{\rm p}$ ($\rm p=1,2,3$), are the
effective (material) areas and area moments about the principal axes of
inertia respectively ($J\equiv I_{1}$).

###  2.4 Linearisation of balance equations

To apply the Newton-Raphson solution procedure, the spatial internal forces
and moments acting on the beam cross-section must be linearised. In this
paper, the procedure described by Marino [29] and the spatial internal forces
and moments,

$\bm{n}=\bm{\Lambda}_{t}\bm{C}_{\rm N}\bm{\Gamma}$ (16)
$\bm{m}=\bm{\Lambda}_{t}\bm{C}_{\rm M}\bm{K}$ (17)

are linearised as follows,

$\bm{n}=\underbrace{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*}}_{\bm{n}^{\rm*}}-\widehat{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*})}\Delta\bm{\psi}+\left(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\right)\Delta\bm{w}^{\prime}+\widehat{(\bm{r}^{*})^{\prime}}\Delta\bm{\psi}$
(18) $\bm{m}=\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}\bm{K}^{*}-\widehat{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}\bm{K}^{*})}\Delta\bm{\psi}+(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}(\bm{\Lambda}_{t}^{\rm*})^{\rm T})\Delta\bm{\psi}^{\prime}$ (19)

Here, the superscript $*$ represents the fields obtained from the previous
Newton-Raphson iteration and $\Delta\bm{w}$ and $\Delta\bm{\psi}$ are the
incremental correction displacement and rotational vectors respectively. The
linearised counterpart of the term $(\bm{r}^{\prime}\times\bm{n})$ in Eq. (7)
is given by,

$(\bm{r}^{\prime}\times\bm{n})=\widehat{(\bm{r}^{*})^{\prime}}\bm{n}^{*}-\widehat{\bm{n}^{*}}\Delta\bm{w}^{\prime}+\widehat{(\bm{r}^{*})^{\prime}}\big{(}\bm{n}-\bm{n}^{\rm*}\big{)}$
(20)

The details of the linearisation are provided in A. The correction vectors
$\Delta\bm{w}$ and $\Delta\bm{\psi}$ are used to calculate the new mean line
displacement vector $\bm{w}$ and the new rotation matrix $\bm{\Lambda}$ at the
end of each Newton-Raphson iteration according to the formulae,

$\bm{w}=\bm{w}^{*}+\Delta\bm{w}$ (21)
$\bm{\Lambda}=\exp(\widehat{\Delta\bm{\psi}})\bm{\Lambda}^{*}$ (22)

where the exponentiation of the skew-symmetric tensor
$\widehat{\Delta\bm{\psi}}\in\mathfrak{so(3)}$ is evaluated by the Rodrigues’
formula (Eq. 12) to compute the rotation matrix $\bm{\Lambda}\in SO(3)$.

Remark:

The type of rotation matrix interpolation adopted in this work corresponds to
the variant of Type-I as presented by Romero [34] and the objectivity of the
corresponding strain measures is verified using a numerical test case as
presented in Section 4.1.

##  3 NUMERICAL MODEL

In this section, cell-centred FV discretisation of the computational domain
and the governing balance equations is discussed. Unlike the standard FE
approach, where the governing equation in strong form is cast into its
equivalent weak form, the current FV method directly discretises the strong
form of the integral equation. The discretisation procedure is separated into
two distinct parts: discretisation of the solution domain and discretisation
of the governing equations.

###  3.1 Solution domain discretisation

For the quasi-static case, solution domain discretisation implies space
discretisation where loads are applied in pseudo-time increments. The beam
body in its reference configuration is divided into a finite number of uniform
segments or control volumes (CVs) as is shown in Fig. 2. A typical
computational stencil (Fig. 2) consists of the central CV (cell) of length
$L_{\rm C}$ with computational node $\rm C$ , located at the cell centroid,
bounded by two internal faces $\rm w$ and $\rm e$ shared with the
corresponding west and east neighbouring cells, with cell centroids at $\rm W$
and $\rm E$ and lengths, $L_{\rm w}$ and $L_{\rm e}$ from the node $\rm C$
respectively.

Figure 2: Beam body in reference configuration discretised by a finite set of
1-D CVs (cells).

###  3.2 Equation discretisation

For an isolated CV in the deformed configuration (Fig. 3), the integral form
of the balance equations (Eqs. 8 and 9) can be discretised over a CV as
following,

$\displaystyle\begin{split}\int_{L}\bm{n}^{\prime}\textrm{d}L&=\sum_{f}\bm{n}_{f}\bm{d}_{f}=\bm{n}_{\textrm{e}}-\bm{n}_{\textrm{w}}\quad;\quad\int_{L}\bm{f}\textrm{d}L\approx\bm{f}_{\rm
C}L_{\rm C}\\\\[5.69046pt] \implies&\bm{n}_{\rm e}-\bm{n}_{\rm w}+\bm{f}_{\rm
C}L_{C}=0\end{split}$ (23)

where the subscript $f$ denotes the values evaluated at cell faces and
subscript $\rm C$ represents the values at cell-centre $\rm C$; $\bm{d}_{f}$
is the unit normal of the face. The term $\bm{f}$ is assumed to have a linear
variation across the CV and hence, can be approximated by the mid-point rule.

Similarly, the discretised version of the moment balance equation about the
cell-centre $\rm C$ takes the form,

$\displaystyle\begin{split}\int_{L}\bm{m}^{\prime}\textrm{d}L=\sum_{f}\bm{m}_{f}\bm{d}_{f}&=\bm{m}_{\textrm{e}}-\bm{m}_{\textrm{w}}\\\\[5.69046pt]
\int_{L}(\bm{r}^{\prime}\times\bm{n})\textrm{d}L=\sum_{f}(\bm{r}^{\prime}_{\rm
f}\times\bm{n}_{\rm f})L_{f}\bm{d}_{f}&=\frac{1}{2}L_{\rm
e}(\bm{r}^{\prime}_{\rm e}\times\bm{n}_{\rm e})+\frac{1}{2}L_{\rm
w}(\bm{r}^{\prime}_{\rm w}\times\bm{n}_{\rm w})\\\\[5.69046pt]
\implies\bm{m}_{\rm e}-\bm{m}_{\rm w}+\frac{1}{2}L_{\rm
e}(\bm{r}^{\prime}_{\rm e}\times\bm{n}_{\rm e})&+\frac{1}{2}L_{\rm
w}(\bm{r}^{\prime}_{\rm w}\times\bm{n}_{\rm w})+\bm{t}_{\rm
C}L_{C}=0\end{split}$ (24)

where, $L_{f}$ is the distance from centre $\rm C$ to faces $f$; $\bm{t}$ is
also assumed to have a linear variation across the CV and hence is
approximated by the mid-point rule.

Figure 3: Balance of forces and moments on an isolated CV in the deformed
configuration

Unlike the traditional FE approach where the primary unknowns stored-at the
nodal points-are used along with element shape functions to interpolate values
to any other point, in the FV method, the unknown fields are computed at
computational nodal points (cell-centres) and interpolated, by appropriate
discretisation schemes, to the face centres and elsewhere. Subsequently, all
the terms in the discretised equilibrium equations (Eqs. 23 and 24) are
substituted by their linearised counterpart (Eqs. 18 \- 20).

To evaluate the internal forces and moments at cell-faces from the linearised
equations mentioned in Section 2, the mean line displacement correction vector
$\Delta\bm{w}$ and its derivative $\Delta\bm{w}^{\prime}$, and the value of
cross-section rotational correction vector $\Delta\bm{\psi}$ and its
derivative $\Delta\bm{\psi}^{\prime}$ have to be approximated at the face
centres in terms of the cell-centre values. All the cell-face derivatives are
approximated by the central finite difference scheme while cell-face values
are approximated by the linear interpolation. The cell-face values on the
internal faces $w$ and $e$ are linearly interpolated using cell-centre values
as,

$\Big{[}(\cdot)\Big{]}_{e}=\gamma_{\rm e}\Big{[}(\cdot)\Big{]}_{\rm
E}+(1-\gamma_{\rm e})\Big{[}(\cdot)\Big{]}_{\rm
C}\quad;\quad\Big{[}(\cdot)\Big{]}_{w}=\gamma_{\rm
w}\Big{[}(\cdot)\Big{]}_{\rm W}+(1-\gamma_{\rm w})\Big{[}(\cdot)\Big{]}_{\rm
C}$ (25)

where subscripts $\rm E$ and $\rm W$ represent values at the neighbouring
cell-centres; $\Big{[}(\cdot)\Big{]}_{e}$ and $\Big{[}(\cdot)\Big{]}_{w}$ are
interpolated values at internal faces $e$ and $w$ respectively and
$\gamma_{\rm e}$, $\gamma_{\rm w}$ are the weighing factors given by,

$\gamma_{\rm e}=\frac{1}{2}\frac{L_{\rm e}}{L_{\rm C}}\quad;\quad\gamma_{\rm
w}=\frac{1}{2}\frac{L_{\rm w}}{L_{\rm C}}$ (26)

The notations $\Big{[}(\cdot)\Big{]}_{\rm C}$, $\Big{[}(\cdot)\Big{]}_{\rm W}$
and $\Big{[}(\cdot)\Big{]}_{\rm E}$ denote the cell values at the location
$\rm C$ and the neighbouring cells $\rm W$ and $\rm E$ (Fig. 2). The cell face
derivatives at the internal faces $\rm w$ and $\rm e$ given by
$\Big{[}(\cdot)\Big{]}^{\prime}_{e}$ and $\Big{[}(\cdot)\Big{]}^{\prime}_{w}$
are approximated using the central finite difference scheme as,

$\Big{[}(\cdot)\Big{]}^{\prime}_{e}=\frac{\Big{[}(\cdot)\Big{]}_{\rm
E}-\Big{[}(\cdot)\Big{]}_{\rm
C}}{L_{C}}\quad;\quad\Big{[}(\cdot)\Big{]}^{\prime}_{w}=\frac{\Big{[}(\cdot)\Big{]}_{\rm
C}-\Big{[}(\cdot)\Big{]}_{\rm W}}{L_{C}}$ (27)

In the discretised equilibrium equations, all the terms are explicitly
computed except for the primary unknowns, the correction vectors of
incremental displacement and incremental rotation ($\Delta\bm{w}$ and
$\Delta\bm{\psi}$), which are treated implicitly and are evaluated at the
computational cell-centres. As mentioned in Section 2.4, at the end of every
Newton-Raphson iteration, $\Delta\bm{w}$ and $\Delta\bm{\psi}$ are used to
update the previously converged displacement and rotation fields (indicated by
superscript $*$) and evaluate the deformed mean line position vector
$\bm{r}(s)$ and the new rotation matrix $\bm{\Lambda}$ according to the Eqs.
21, 4 and 22. The interpolated cell-face values and their corresponding cell-
face gradients at internal faces $w$ and $e$ are used to compute the strain
measures, $\bm{\Gamma}$ and $\bm{K}$ (Eqs. 10 and 13), which in turn are used
to evaluate the spatial forces and moments at the cell-faces (Eqs. 16 and 17).
The employed discretisation provides a nominally $2$nd order accurate
approximation for displacements.

###  3.3 Initial and Boundary Conditions

External forces and moments are applied in pseudo-time increments. For the
beam body, there are two boundary faces, one at the left and the other at the
right of the beam. For a Dirichlet boundary condition, in the discretised
governing equations (Eqs. 23 and 24), the displacement/rotation component at
the face centre of the boundary are directly replaced by the user defined
value. For Neumann boundary conditions, the values of forces/moments at either
boundary locations, are directly specified in the Eqs. 23 and 24. Following
solution of the linear system, the corresponding displacements and rotations
at that boundary location are obtained by linear extrapolation from the
interiors of the solution domain using Eqs. 18 and 19. For instance, Fig. 4
shows the boundary face b and the neighbouring cell-centre C. For a specified
force $\bar{\bm{n}}$ and/or moment $\bar{\bm{m}}$ on the boundary, the
incremental displacements/rotations can be obtained from Eqs. 18 and 19 as,

Figure 4: CV at the right boundary with a specified force $\bar{\bm{n}}$ and
moment $\bar{\bm{m}}$

$\displaystyle\bar{\bm{n}}$ $\displaystyle=\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*}+\left(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\right)\Delta\bm{w}^{\prime}_{\textrm{b}}+\Big{[}\left(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\right)\widehat{(\bm{r}^{*})^{\prime}}-\widehat{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*})}\Big{]}\Delta\bm{\psi}_{\rm b}$ (28a)
$\displaystyle\bar{\bm{m}}$ $\displaystyle=\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}\bm{K}^{*}-\widehat{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}\bm{K}^{*})}\Delta\bm{\psi}_{\textrm{b}}+(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}(\bm{\Lambda}_{t}^{\rm*})^{\rm T})\Delta\bm{\psi}^{\prime}_{\textrm{b}}$
(28b)

where, the cell-face derivatives $\Delta\bm{w}^{\prime}_{\textrm{b}}$ and
$\Delta\bm{\psi}^{\prime}_{\textrm{b}}$ are given by,

$\displaystyle\Delta\bm{w}^{\prime}_{\textrm{b}}=\frac{\Delta\bm{w}_{\rm
b}-\Delta\bm{w}_{\rm C}}{\Delta x_{\rm
b}}\quad;\quad\Delta\bm{\psi}^{\prime}_{\textrm{b}}=\frac{\Delta\bm{\psi}_{\rm
b}-\Delta\bm{\psi}_{\rm C}}{\Delta x_{\rm b}}$

The value of $\Delta\bm{\psi}^{\prime}_{\rm b}$ is substituted into Eq. 28b
and the incremental rotation vector $\Delta\bm{\psi}_{\rm b}$ is evaluated as,

$\Delta\bm{\psi}_{\rm b}=\Bigg{(}\frac{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}(\bm{\Lambda}_{t}^{\rm*})^{\rm T}}{\Delta x_{\rm
b}}-(\widehat{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}\bm{K}^{\rm*}})\Bigg{)}^{-1}\Bigg{[}\bar{\bm{m}}-\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}\bm{K}^{*}+\frac{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}(\bm{\Lambda}_{t}^{\rm*})^{\rm T}}{\Delta x_{\rm b}}\Delta\bm{\psi}_{\rm
C}\Bigg{]}$

Using $\Delta\bm{\psi}_{\rm b}$ and the value of
$\Delta\bm{w}^{\prime}_{\textrm{b}}$ as shown above, the incremental
displacement $\Delta\bm{w}_{\rm b}$ is then calculated as,

$\Delta\bm{w}_{\rm b}=\Delta\bm{w}_{\rm
C}+\Bigg{(}\frac{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm T}}{\Delta x_{\rm
b}}\Bigg{)}^{-1}\Bigg{[}\bar{\bm{n}}-\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*}-\Big{[}\left(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\right)\widehat{(\bm{r}^{*})^{\prime}}-\widehat{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*})}\Big{]}\Delta\bm{\psi}_{\rm b}\Bigg{]}$

###  3.4 Solution procedure

The final form of the discretised equilibrium equations, with appropriate
discretisation schemes described in Section 3.2, for a typical computational
node $\rm C$ read as follows,

$\bm{A}_{\rm W}\begin{bmatrix}(\Delta\bm{w})_{\rm W}\\\ (\Delta\bm{\psi})_{\rm
W}\end{bmatrix}+\bm{A}_{\rm C}\begin{bmatrix}(\Delta\bm{w})_{\rm C}\\\
(\Delta\bm{\psi})_{\rm C}\end{bmatrix}+\bm{A}_{\rm
E}\begin{bmatrix}(\Delta\bm{w})_{\rm E}\\\ (\Delta\bm{\psi})_{\rm
E}\end{bmatrix}=\begin{bmatrix}(\bm{R}^{\bm{w}})_{\rm C}\\\
(\bm{R}^{\bm{\psi}})_{\rm C}\end{bmatrix}$ (29)

where $\bm{A}_{\rm C}$ is a coefficient matrix containing the contributions of
node $\rm C$ while the matrices $\bm{A}_{\rm W}$ and $\bm{A}_{\rm E}$
representing the interactions with the neighbouring cell centres $\rm W$ and
$\rm E$. The right hand side of the Eq. 29 is the source vector contribution.
All the coefficient matrices are $(6\times 6)$ dense coupled matrices with the
primary unknowns being $\Delta\bm{w}$ and $\Delta\bm{\psi}$. The three
components of the mean line displacement correction and cross-section rotation
vectors have to be solved in coupled manner. The detailed structure of the
diagonal and off-diagonal coefficient matrices is provided in the B.

The linearised Eqs. 29 are assembled for all CVs forming a system of equations
given by,

$\big{[}\bm{A}\big{]}\big{[}\bm{\phi}\big{]}=\big{[}\bm{R}\big{]}$ (30)

resulting in $6M\times 6M$ sparse matrix $[\bm{A}]$ with weak diagonal
dominance, with $M$ being the total number of CVs. The coefficients
$\bm{A}_{\rm C}$ constitute the diagonal whereas matrices $\bm{A}_{\rm W}$ and
$\bm{A}_{\rm E}$ contribute to the off-diagonal terms of $[\bm{A}]$. The
solution vector $\big{[}\bm{\phi}\big{]}$ contains the primary unknowns
$\Delta\bm{w}$ and $\Delta\bm{\psi}$, and $\big{[}\bm{R}\big{]}$ is the source
vector containing the explicit discretised terms and boundary condition
contributions. The final system of linearised algebraic equations obtained by
assembling Eqs. 29 for all control volumes in the mesh is successfully solved
using the block variant of the Thomas algorithm.

For every pseudo-time increment, the coupled equations are iteratively solved
by Newton-Raphson procedure, until a convergence tolerance of order $10^{-10}$
is achieved with the maximum number of allowable iterations set to $30$. For
calculating the normalised (with respect to their magnitudes) residuals, after
each iteration, for all CVs, a list of the magnitudes of the normalised
correction vectors, $\Delta\bm{w}$ and $\Delta\bm{\psi}$ is calculated and the
global maximum of the two magnitudes is selected
($\textrm{Residuals}=max\Big{[}|\Delta\bm{w}|,|\Delta\bm{\psi}|\Big{]}$) and
the equations are solved until the convergence tolerance is achieved. The
current method has been implemented in open-source software OpenFOAM ([39])
(version foam-extend-4.1), exploiting the developed object oriented FV
procedures. The overall solution procedure is summarised in Algorithm 1 and
the procedure to update the kinematic quantities and stress resultants is
given in Algorithm 2.

Algorithm 1 Solution procedure; $(\cdot)_{f}$ : fields at the face centres

1:Set $\bm{\Lambda}_{0}=\mathbf{I}$ $\triangleright$ Reference mean tangent
line and reference rotation matrix

2:if initially curved beam then

3: Set $\bm{r}_{0}$ and $\bm{\Lambda}_{0}$ at cell faces $\triangleright$
Initial undeformed mean line and initial rotation matrix; specific to test
case

4: $\implies(\bm{r}_{0})_{f},(\bm{\Lambda}_{0})_{f}$

5:end if

6:Set $\bm{\Lambda}=\mathbf{I}$, $\bm{\Lambda}_{f}=\mathbf{I}$ and
$\bm{\Lambda}_{t}=\bm{\Lambda}(\bm{\Lambda}_{0})_{f}$ $\triangleright$
Relative rotation and total rotation matrix

7:Calculate $\big{(}\bm{r}_{0}^{\prime}\big{)}_{f}$ from $(\bm{r}_{0})_{f}$

8:Set $\bm{g}_{0}=\big{(}\bm{r}_{0}^{\prime}\big{)}_{f}$ $\triangleright$
Initial tangent vector at face centres

9:for all pseudo-time steps do

10: while residuals are not converged do

11: Coupled linearised balance equations: assemble and solve for
$\Delta\bm{w}$ and $\Delta\bm{\psi}$ $\triangleright$ Eq. 30

12: Interpolate cell-centre displacements $\Delta\bm{w}$ and rotation vectors
$\Delta\bm{\psi}$ to face values $\triangleright$ Eq. 25

13: Update kinematic quantities and stress resultants $\triangleright$
Algorithm 2

14: end while

15:end for

Algorithm 2 Update of the kinematic quantities and stress resultants;
$(\cdot)_{f}$ : fields at the face centres
$(\cdot)^{*}$ : fields calculated in the previous Newton-Raphson iteration

1:$\bm{w}=\bm{w}^{*}+\Delta\bm{w}$ $\triangleright$ Update displacement vector
(Eq. 21)

2:$\bm{\psi}=\bm{\psi}^{*}+(\bm{\Lambda}^{*})^{T}\Delta\bm{\psi}$
$\triangleright$ Update rotation vector

3:Calculate $(\Delta\bm{\psi})_{f}$ from $\Delta\bm{\psi}$ $\triangleright$
Interpolate rotation vector to the face centre (Eq. 25)

4:$\Delta\bm{T}_{f}=\frac{\sin\|\Delta\bm{\psi}_{f}\|}{\|\Delta\bm{\psi}_{f}\|}\mathbf{I}+\Big{(}1-\frac{\sin\|\Delta\bm{\psi}_{f}\|}{\|\Delta\bm{\psi}_{f}\|}\Big{)}\frac{\Delta\bm{\psi}_{f}\cdot\Delta\bm{\psi}^{T}_{f}}{\|\Delta\bm{\psi}_{f}\|^{2}}+\Big{(}\frac{1-\cos\|\Delta\bm{\psi}_{f}\|}{\|\Delta\bm{\psi}_{f}\|^{2}}\Big{)}\widehat{\Delta\bm{\psi}_{f}}^{2}$
$\triangleright$ Calculate incremental tangent operator (Eq. 14)

5:

6:$\bm{K}_{f}=\bm{K}_{f}^{*}+\big{(}\bm{\Lambda}_{0}^{T}\big{)}_{f}(\bm{\Lambda}^{*})^{T}_{f}\Delta\bm{T}_{f}^{T}(\Delta\bm{\psi}^{\prime})$
$\triangleright$ Calculate (material) rotational strains (Eq. 13)

7:$(\Delta\bm{\Lambda})_{f}=\mathbf{I}+\Big{(}\frac{\sin\|\Delta\bm{\psi}_{f}\|}{\|\Delta\bm{\psi}_{f}\|}\Big{)}\Delta\bm{\psi}_{f}\cdot\Delta\bm{\psi}^{T}_{f}+\Big{(}\frac{1-\cos\|\Delta\bm{\psi}_{f}\|}{\|\Delta\bm{\psi}_{f}\|^{2}}\Big{)}\widehat{\Delta\bm{\psi}_{f}}^{2}$
$\triangleright$ Calculate face centre incremental rotation matrix (Eq. 12)

8:

9:$\bm{\Lambda}_{f}=(\Delta\bm{\Lambda})_{f}\bm{\Lambda}_{f}^{*}$
$\triangleright$ Update the face centre rotation matrix (Eq. 22)

10:$\Delta\bm{\Lambda}=\mathbf{I}+\Big{(}\frac{\sin\|\Delta\bm{\psi}\|}{\|\Delta\bm{\psi}\|}\Big{)}\Delta\bm{\psi}\cdot\Delta\bm{\psi}^{T}+\Big{(}\frac{1-\cos\|\Delta\bm{\psi}\|}{\|\Delta\bm{\psi}\|^{2}}\Big{)}\widehat{\Delta\bm{\psi}}^{2}$
$\triangleright$ Calculate the incremental cell centre rotation matrix (Eq.
12)

11:$\bm{\Lambda}=(\Delta\bm{\Lambda})\bm{\Lambda}^{*}$ $\triangleright$ Update
the cell centre rotation matrix

12:$\big{(}\bm{r}^{\prime}\big{)}_{f}=\big{(}\bm{r}_{0}\big{)}_{f}+\bm{w}^{\prime}$
$\triangleright$ Calculate derivative of deformed mean line at face centre

13:$\bm{\Gamma}_{f}=\big{(}\bm{\Lambda}_{0}^{T}\big{)}_{f}\bm{\Lambda}_{f}^{T}\big{(}\bm{r}^{\prime}\big{)}_{f}-\big{(}\bm{\Lambda}_{0}^{T}\big{)}_{f}\big{(}\bm{r}^{\prime}_{0}\big{)}_{f}$
$\triangleright$ Calculate (material) translational strains (Eq. 10)

14:$\bm{n}_{f}=\bm{\Lambda}_{f}\big{(}\bm{\Lambda}_{0}\big{)}_{f}\bm{C}_{\rm
N}\bm{\Gamma}$ $\triangleright$ Calculate spatial forces at face centres (Eq.
16)

15:$\bm{m}_{f}=\bm{\Lambda}_{f}\big{(}\bm{\Lambda}_{0}\big{)}_{f}\bm{C}_{\rm
M}\bm{K}$ $\triangleright$ Calculate spatial moments at face centres (Eq. 17)

##  4 VERIFICATION TEST CASES

In this section, the capabilities of the developed FV methodology are
investigated on the five complementary benchmark cases:

  1. 1.

Rigid rotation of an initially curved beam: this case verifies the objectivity
of the adopted strains and rotation interpolation measures.

  2. 2.

In-plane bending of a cantilever beam subjected to concentrated moment at one
end: this case provides an analytical solution to check the accuracy and order
of accuracy of the proposed methodology.

  3. 3.

Out-of-plane bending of the cantilever to form a helix due to a concentrated
moment and an out-of-plane force: this case with complex loading conditions,
tests the ability of the numerical solver to cater for large rotations and 3-D
deformation.

  4. 4.

A cantilever, initially bent into a $45^{\circ}$ arc in $xy$-plane, subjected
to a force along the $z$-direction: this is another benchmark curved beam case
used by many authors to establish the accuracy and order of accuracy of the
numerical model.

  5. 5.

Deep-circular arch with a concentrated in-plane force at the crown location:
this unsymmetrical circular arch case has exact solutions available and is
used to assess the ability of the solver to obtain the critical buckling load
value.

All the test cases have been executed using a quad core CPU with hyper-
threading (Intel Core(TM) i7-8565U CPU with base frequency $1.80$GHz and
maximum turbo frequency $4.6$GHz).

###  4.1 Rigid rotation of an initially curved beam

For the FV formulation presented in this paper and the type of the rotation
interpolation adopted, the objectivity of the conjugate strain measures is
verified using the numerical test case first presented in Meier et al. [30,
31] as shown in Fig. 5(a). For an initially curved beam with a centerline
configuration of a quarter circle with radius $R=100$ m, initial length,
$L=\pi R/2=157.079$ m, discretised into $10$ CVs, a rotation of $20\pi$ is
gradually applied about the global x-axis in $100$ load increments. The
mechanical properties for this test case are, $E=1$ GPa and $G=0.5$ GPa. The
nature of the boundary condition at the clamped left end of the beam is such
that, this prescribed rotation should only cause rigid body transformation in
the beam about the $x$-axis with no accumulation of internal strain energy and
no physical deformation, i.e. $w_{x}=w_{y}=w_{z}=0$ and $\psi_{y}=\psi_{z}=0$.

(a)

(b)

Figure 5: Rigid rotation of an initially curved beam: (a) test case setup, (b)
deformed configuration of the beam for a rotation angle, $\psi_{x}=19\pi$

Fig. 5(b) shows the deformed beam configuration for an angle $\psi_{x}=19\pi$
about the x-axis and it is evident that there is no physical deformation. For
the test case, execution time is approximately $1.5$ s and requires $9$ outer
iterations on average to converge per pseudo-time increment. For each CV, the
strains $\bm{\Gamma}$ and $\bm{K}$ are used to calculate the internal energy
as per Eq. 15 and summed over all the CVs to obtain the total accumulated
strain energy of the beam due to the applied rotation. The total energy is
found to be zero to machine precision (in the order of $10^{-31}$) and hence,
the objectivity of the adopted strain measures is confirmed.

###  4.2 In-plane bending of a cantilever beam subjected to concentrated
moment at one end

This in-plane pure flexural bending of a cantilever case has been investigated
by a number of authors [37, 21]. An initially straight cantilever of length
$L=10$ m is bent into a circle by applying a concentrated moment at one end.
The mechanical properties available in the literature for this test case are
$EA=10^{4}$ N , $GA_{2}=GA_{3}=5000$ N, $EI_{2}=EI_{3}=100$ $\textrm{Nm}^{2}$,
$GJ=100$ $\textrm{Nm}^{2}$ respectively. Hence, the cross-section radius of
the beam, $r=0.2$ m and the Young’s modulus $E=7.95\times 10^{4}$ Pa are
assumed in a way to achieve these desired numerical values. The Poisson’s
ratio $\nu$ is taken as zero. According to the classic Euler formula, the
analytical solution for pure beam-bending is given by,

$\psi_{z}=\frac{M_{z}L}{EI}\quad;\quad
w_{x}=L-\frac{L}{\psi_{z}/2}sin\frac{\psi_{z}}{2}cos\frac{\psi_{z}}{2}\quad;\quad
w_{y}=\frac{L}{\psi_{z}/2}\Big{(}sin\frac{\psi_{z}}{2}\Big{)}^{2}$ (31)

where $\psi_{z}$, $w_{x}$ and $w_{y}$ are the rotation and in-plane
displacements of the beam tip.

Figure 6: In-plane bending of a cantilever beam subjected to concentrated
moment at one end: deformed configuration for $\psi_{z}=2\pi$

The computational domain is discretised into $10$ CVs and a moment of
$M_{z}=20\pi$ Nm $(\psi_{z}=2\pi)$ is applied at the right end in $4$
increments and the Newton-Raphson solution procedure converges in an average
of $6$ outer iterations per pseudo-time increment. Fig. 6 shows the deformed
configuration of the cantilever for the applied moment. Although the solver
takes larger number of iterations to converge than is reported in the
literature (Simo et al. [37] reported only two iterations), the execution time
is less than $0.5$ s. For an applied moment, $M_{z}=2.5\pi$ Nm and a mesh of
$10$ CVs, the tip in-plane displacements are found to be $w_{x}=-1.00146$ m
and $w_{y}=3.72731$ m respectively, which vary from the analytical solution
(Eq. 31) by only $0.4\%$ and $0.05\%$ respectively. The percentage relative
error is calculated by the formula,

$\textrm{\% relative
error}=\bigg{|}\frac{\xi^{\textrm{num}}-\xi^{\textrm{ref}}}{\xi^{\textrm{ref}}}\bigg{|}\times
100\%$ (32)

where, $\xi^{\textrm{num}}$ denotes the numerical value and
$\xi^{\textrm{ref}}$ is the reference/analytical result. Fig. 7 shows the
percentage mesh error convergence of the in-plane displacements for successive
reduction of mesh sizes, i.e. $5,10,20$ and $40$ CVs; a quadratic order of
error convergence is observed.

Figure 7: In-plane bending of a cantilever beam subjected to concentrated
moment at one end: mesh convergence plot of in-plane displacements for an
applied moment $M=2.5\pi$ Nm

###  4.3 Out-of-plane bending of the cantilever to form a helix due to a
concentrated moment and an out-of-plane force

For this case, the previous problem is extended by applying a concentrated
out-of-plane force at the free end of the beam, $\bm{n}=[0,0,50]$ N along with
a rotation $\psi_{z}=20\pi$ ($\bm{m}=[0,0,200\pi]$ Nm) about the z-axis. Thus,
an initially straight cantilever beam is bent to a circular helix shape. The
constitutive matrices are taken similar to the previous case, i.e.
$\bm{C}_{\rm N}=\rm diag[10^{4},5\times 10^{3},5\times 10^{3}]$ N and
$\bm{C}_{\rm M}=\rm diag[10^{2},10^{2},10^{2}]$ Nm2. The cross-section radius
is reduced to $r=0.02$ m and the Young’s modulus is taken as, $E=7.95\times
10^{8}$ Pa, so that the beam does not overlap into its own body while forming
the helix. For comparing the achieved numerical results with those published
by Ibrahimbegovic [20], the beam of initial length $L=10$ m is discretised
into $100$ uniform CVs.

(a)

(b)

Figure 8: Out-of-plane bending of the cantilever to form a helix due to a
concentrated moment and an out-of-plane force: (a) initial beam configuration
(b) deformed beam for $20\%$ (left), $30\%$ (middle) and $40\%$ (right) of the
applied load Figure 9: Out-of-plane bending of the cantilever to form a helix
due to a concentrated moment and an out-of-plane force: free-end displacement
component in the direction of applied load ($w^{\textrm{ref}}_{z}$ adopted
from [20])

The external loads to the beam are applied in $2000$ pseudo-time increments
and the final load is reached at pseudo time-value of $10$ s. The Newton-
Raphson solution procedure required $6$ outer iterations on average to
converge per pseudo time increment and $60$ s of total execution time. Fig.
8(a) shows the initial beam and Fig. 8(b) shows the deformed shape of the beam
for $20\%$, $30\%$ and $40\%$ of the total loading respectively. Fig. 9
presents the free-end displacement component $w_{z}$ as a function of the
applied load. From Fig. 9, it is evident that with increasing rotations, the
out-of-plane displacement oscillates about the $z$-axis crossing the zero
value. This ability of the solver to capture the complex deformation and
oscillatory motion of the beam is possible because of the incremental rotation
vector based parametrisation and the results (Fig. 9) are in good agreement
with those published in [20]; note the same test case failed to converge when
a non-incremental formulation was initially trialed.

###  4.4 A cantilever, initially bent into a $45^{\circ}$ arc in $xy$-plane,
subjected to a force along the $z$-direction

This case consists of an initially curved cantilever beam with an applied
vertical concentrated force along $z$-axis. Several authors have reproduced
this benchmark case [2, 37, 4, 7, 21]. The cantilever, which has a unit square
cross-section, is initially bent into a $45^{\circ}$ arc of radius $100$ m,
length $L=\pi R/4=78.54$ m and then a vertical force, $n_{z}=600$ N is
applied. The mechanical properties adopted are $E=1\times 10^{7}$ $Pa$ and
$\nu=0$. For a spatial discretisation of $8$ CVs, the tip end displacements
values are $w_{x}=-23.5281$ m, $w_{y}=-13.5415$ m and $w_{z}=53.082$ m
respectively. Simo et al. [37] apply a load in three steps ($300$, $150$ and
$150$ N) and take 27 cumulative iterations to attain convergence. The current
model takes $34$ cumulative iterations for the same loading increment but less
than $0.5$ s of execution time. On the other hand, for a systematic loading in
$10$ increments and discretisation of $10$ CVs, it take an average of $6$
outer iterations to converge per pseudo time increment. Fig. 10 shows the
initial and the deformed configuration of the beam and the displacement
component $w_{z}$ for $10$ CVs.

Figure 10: A cantilever, initially bent into a $45^{\circ}$ arc in $xy$-plane, subjected to a force along the $z$-direction: initial and deformed configuration for an applied load of $600$ N Table 1: Displacement values observed for a curved cantilever in $xy-$plane subjected to a force along $z-$axis Numerical Results | $|w_{x}|$ | $|w_{y}|$ | $|w_{z}|$
---|---|---|---
Present | $23.540$ | $13.564$ | $53.225$
Bathe and Bolourchi [2] | $23.5$ | $13.4$ | $53.4$
Simo and Vu-Quoc [37] | $23.48$ | $13.50$ | $53.37$
Cardona and Geradin [4] | $23.67$ | $13.73$ | $53.50$
Crisfield [7] | $23.87$ | $13.63$ | $53.71$
Ibrahimbegović [21] | $23.697$ | $13.668$ | $53.498$
Figure 11: A cantilever, initially bent into a $45^{\circ}$ arc in $xy$-plane,
subjected to a force along the $z$-direction: mesh convergence of the
displacement components for the applied load, $n_{z}=600$ N Figure 12: A
cantilever, initially bent into a $45^{\circ}$ arc in $xy$-plane, subjected to
a force of $3000$ N along the $z$-direction: the end-tip displacements for
increasing load values Figure 13: A cantilever, initially bent into a
$45^{\circ}$ arc in $xy$-plane: mesh convergence of the force $n_{z}$ for an
applied displacement, $w_{x}=-23.5607$ m, $w_{y}=-13.6048$ m, $w_{z}=53.4756$
m

For this test case, the reported displacements in the literature are for a
discretisation of $8$ linear FE elements and the current numerical results
(for $10$ CVs) are found to be comparable with the ones presented in the
literature (Table 1). Fig. 11 presents the % mesh error convergence of the
displacement components for successive mesh size reductions, viz. $5,10,20$
and $40$ CVs; a quadratic order of error convergence is observed. The
analytical solutions for this test case are not available and the numerical
results reported in literature are similar. To calculate the % mesh error (Eq.
32), the displacement values corresponding to a finer mesh of $80$ CVs are
adopted as reference results here. For comparing the displacement versus load
curves presented in Simo et al [37], the beam is discretised into $8$ CVs and
an end force of $3000$ N is applied to its free end at $30$ N load increments.
Fig. 12 presents the displacement components of the free end of the beam and
the results are found to be in close agreement with those reported by Simo et
al. [37].

To investigate the order of accuracy of the solver in calculating forces, the
present test case is reproduced by applying displacements at the right end
instead of forces. To obtain mesh independent displacements, a fine mesh of
$1280$ CVs is selected and for an applied force of $600$ N along the $z$-axis,
the values observed are $w_{x}=-23.5607$ m, $w_{y}=-13.6048$ m,
$w_{z}=53.4756$ m. On applying these displacements to the same test case, the
exact force vector $(0\ 0\ 600)$ N should be retrieved. Consequently, the test
case is run by applying the displacements $(-23.5607,-13.6048,53.4756)$ m for
different mesh sizes, viz. $10$, $20$, $40$, $80$ and $160$ CVs. Fig. 13 shows
the % mesh error convergence of the vertical force ($n_{z}$) for successive
mesh size reductions; a quadratic order of error convergence is observed.
Since the true value of forces $n_{x}$ and $n_{y}$ expected for the given
displacements is zero, relative error is not defined; however, the force
components obtained for successive mesh reductions are seen to go to zero
(Table 2).

Table 2: Curved cantilever in $xy-$plane subjected to a force along $z-$axis: force values observed for different mesh sizes CVs |  $n_{x}$ (N) |  $n_{y}$ (N)
---|---|---
$10$ | $53.98$ | $24.76$
$20$ | $12.81$ | $5.91$
$40$ | $3.17$ | $1.47$
$80$ | $0.80$ | $0.37$
$160$ | $0.21$ | $0.09$

###  4.5 Deep-circular arch with a concentrated in-plane force at the crown
location

The instabilities of asymmetric clamped-hinged arches were first investigated
by DaDeppo and Schmidt [8] and then by many other authors [37, 30, 33]. In the
current work, a $215^{\circ}$ deep circular arch with a unit circular cross-
section, hinged at the left end ($w_{x}=w_{y}=w_{z}=0$) and clamped at the
right ($w_{x}=w_{y}=w_{z}=0$ and $\psi_{x}=\psi_{y}=\psi_{z}=0$), having a
radius $R=100$ m and length $375.245$ m is considered and a force-displacement
behaviour of the arch is studied until buckling. A concentrated point load
$P\equiv(0\quad n_{y}\quad 0)$ acts at the crown location of the arch. The
adopted mechanical properties are: $EI_{2}=EI_{3}=GJ=1\times 10^{4}$ Nm2.

Figure 14: Deep-circular arch with a concentrated in-plane force at the crown
location: initial and deformed configuration

In accordance with the literature [37], the entire spatial domain is
discretised into $40$ CVs and Fig. 14 shows the initial and the deformed
configuration of the arch. For an increasing load at $1$ N increments up to
$8$ N and $0.005$ N increments beyond $8$ N , the predictive critical buckling
load is found to be is $9.085$ N exact value where the value reported by
DaDeppo and Schmidt [8] is $8.97$ N. The solver, being quasi-static, does not
converge for the unstable post-buckling analysis of the circular arch; beyond
the critical buckling load, a dynamic analysis or special procedures would be
required. For the loading values less than $8$ N, the residuals converge at an
average of $8$ iterations per pseudo time increment; for loads higher than $8$
N, the residuals converge in an average of $17$ iterations. The total time of
computation is less than $1$ s. Fig. 15 shows the variation of the
dimensionless quantity $\frac{PR^{2}}{EI}$ versus the normalised displacements
($w_{x}$ and $w_{y}$ with respect to $R$) and rotation angle $\phi$ of the
tangent vector at the crown location. The results are found to be in good
agreement with the exact solutions presented by DaDeppo and Schmidt [8].

Figure 15: Deep-circular arch with a concentrated in-plane force at the crown
location: the normalised horizontal and vertical displacements of the crown
point ($w_{x}/R,w_{y}/R$) and the rotation angle ($10\phi$) of the crown point
for different loading values of the force

##  5 CONCLUDING REMARKS

This paper is the first to develop and verify a total Lagrangian cell-centred
finite volume methodology for geometrically exact beams with arbitrary initial
curvatures subjected to finite displacements and rotations. The mathematical
model and the corresponding FV-based spatial and equation discretisation are
described in detail. The potential of the developed methodology has been
tested using five complementary benchmark test cases, where the spatial
discretisation has been shown to be second-order accurate for displacements
and forces. For the cases examined, even though the computational time is
quite fast (less than $1$ minute), it is found that the results require a
higher number of outer iterations to converge as compared to those reported in
literature (average of 6 iterations instead of 2 for pure bending of
cantilever beam). Nevertheless, the predictions using the proposed FV
methodology are found to be in good agreement with the analytical solutions
and the FE based numerical results. Extension of the quasi-static FV
formulation to dynamic analysis will be considered in the future work.

## Acknowledgements

This publication has emanated from research [conducted with the financial
support of/supported in part by a grant from] Science Foundation Ireland under
Grant number RC2302_2. For the purpose of Open Access, the author has applied
a CC BY public copyright license to any Author Accepted Manuscript version
arising from this submission. Financial support is gratefully acknowledged
from the Irish Research Council through the Laureate program, grant number
IRCLA/2017/45. Additionally, the authors want to acknowledge project
affiliates, Bekaert, through the Bekaert University Technology Centre (UTC) at
University College Dublin (www.ucd.ie/bekaert), and I-Form, funded by Science
Foundation Ireland (SFI) Grant Number 16/RC/3872, co-funded under European
Regional Development Fund and by I-Form industry partners. Provision of
computational facilities and support from the DJEI/DES/SFI/HEA Irish Centre
for High-End Computing (ICHEC, www.ichec.ie) and ResearchIT Sonic cluster,
funded by UCD IT Services and the Research Office is gratefully acknowledged.

## Conflict of interest

The authors declare that they have no conflict of interest.

##  Appendix A Linearisation of stress-resultants

The linearisation of the stress-resultants involves a systematic use of the
directional derivatives of the kinematic quantities and the strain measures.
For a given beam configuration, the linearisation of perturbed kinematic
quantities
$\big{(}\bm{r}_{\epsilon}=\bm{r}_{0}+\epsilon\Delta\bm{w},(\bm{\Lambda}_{t})_{\epsilon}=\exp(\epsilon\widehat{\Delta\bm{\psi}})\bm{\Lambda}_{t}\big{)}$
yields,

$\frac{d}{d\epsilon}\Big{(}\bm{r}_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}=\Delta\bm{w}\qquad\qquad\frac{d}{d\epsilon}\Big{(}(\bm{\Lambda}_{t})_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}=\widehat{\Delta\bm{\psi}}\bm{\Lambda}_{t}^{\rm*}$
(A.1)

where $\epsilon$ is a small scalar perturbation. The linearised forms for the
perturbed (material) strain measures
$\bm{\Gamma}_{\epsilon}=(\bm{\Lambda}_{t}^{\rm
T})_{\epsilon}\bm{r}^{\prime}_{\epsilon}-\bm{e}_{1}$ and
$\bm{K}_{\epsilon}=(\bm{\Lambda}_{t}^{\rm
T})_{\epsilon}(\bm{\Lambda}_{t}^{\prime})_{\epsilon}$ [37, 29] are,

$\frac{d}{d\epsilon}\Big{(}\bm{\Gamma}_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}=(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\widehat{(\bm{r}^{*})^{\prime}}(\Delta\bm{\psi})+(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\Delta\bm{w}^{\prime}$ (A.2)

$\frac{d}{d\epsilon}\Big{(}\bm{K}_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}=(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\Delta\bm{\psi}^{\prime}$ (A.3)

Using the above values, the spatial force, $\bm{n}=\bm{\Lambda}_{t}\bm{C}_{\rm
N}\bm{\Gamma}$ is linearised as,

$\displaystyle\begin{split}\frac{d}{d\epsilon}\Big{(}\bm{n}_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}&=\frac{d}{d\epsilon}\Big{(}(\bm{\Lambda}_{t})_{\epsilon}\bm{C}_{\rm
N}\bm{\Gamma}_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}\\\\[5.69046pt]
&=\bm{\Lambda}_{t}^{\rm*}\bm{C}_{\rm
N}\bm{\Gamma}^{\rm*}+\frac{d}{d\epsilon}(\bm{\Lambda}_{t})_{\epsilon}\bigg{|}_{\epsilon=0}\bm{C}_{\rm
N}\bm{\Gamma}^{*}+(\bm{\Lambda}_{t}^{\rm*})\bm{C}_{\rm
N}\frac{d}{d\epsilon}(\bm{\Gamma}_{\epsilon})\bigg{|}_{\epsilon=0}\\\\[5.69046pt]
\implies\bm{n}&=\underbrace{\bm{\Lambda}_{t}^{\rm*}\bm{C}_{\rm
N}\bm{\Gamma}^{\rm*}}_{\bm{n}^{*}}-\big{(}\widehat{\bm{\Lambda}_{t}^{\rm*}\bm{C}_{\rm
N}\bm{\Gamma}^{\rm*}}\big{)}\Delta\bm{\psi}+\big{(}\bm{\Lambda}_{t}^{\rm*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\big{)}\Big{[}\Delta\bm{w}^{\prime}+\widehat{(\bm{r}^{*})^{\prime}}(\Delta\bm{\psi})\Big{]}\end{split}$
(A.4)

Similarly, the spatial moment, $\bm{m}=\bm{\Lambda}_{t}\bm{C}_{\rm M}\bm{K}$
can be linearised as,

$\displaystyle\begin{split}\frac{d}{d\epsilon}\Big{(}\bm{m}_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}&=\frac{d}{d\epsilon}\Big{(}(\bm{\Lambda}_{t})_{\epsilon}\bm{C}_{\rm
M}\bm{K}_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}\\\\[5.69046pt]
&=\bm{\Lambda}_{t}^{\rm*}\bm{C}_{\rm
M}\bm{K}^{\rm*}+\frac{d}{d\epsilon}(\bm{\Lambda}_{t})_{\epsilon}\bigg{|}_{\epsilon=0}\bm{C}_{\rm
M}\bm{K}^{*}+(\bm{\Lambda}_{t}^{\rm*})\bm{C}_{\rm
M}\frac{d}{d\epsilon}(\bm{K}_{\epsilon})\bigg{|}_{\epsilon=0}\\\\[5.69046pt]
\implies\bm{m}&=\bm{\Lambda}_{t}^{\rm*}\bm{C}_{\rm
M}\bm{K}^{\rm*}-\big{(}\widehat{\bm{\Lambda}_{t}^{\rm*}\bm{C}_{\rm
N}\bm{K}^{\rm*}}\big{)}\Delta\bm{\psi}+(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}(\bm{\Lambda}_{t}^{\rm*})^{\rm T})\Delta\bm{\psi}^{\prime}\end{split}$ (A.5)

The linearised counterpart of the term $(\bm{r}^{\prime}\times\bm{n})$ in Eq.
(7) is given by,

$\displaystyle\begin{split}\frac{d}{d\epsilon}\Big{(}\bm{r}^{\prime}\times\bm{n}\Big{)}\bigg{|}_{\epsilon=0}&=\frac{d}{d\epsilon}\Big{(}\bm{r}^{\prime}_{\epsilon}\times(\bm{\Lambda}_{t})_{\epsilon}\bm{C}_{\rm
N}\bm{\Gamma_{\epsilon}}\Big{)}\bigg{|}_{\epsilon=0}\\\\[5.69046pt]
&=(\bm{r}^{*})^{\prime}\times(\bm{\Lambda}_{t}^{*})\bm{C}_{\rm
N}\bm{\Gamma^{*}}+\frac{d}{d\epsilon}(\bm{r}^{\prime}_{\epsilon})\bigg{|}_{\epsilon=0}\times\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*}\\\\[5.69046pt]
&\qquad+(\bm{r}^{*})^{\prime}\times\frac{d}{d\epsilon}\Big{(}(\bm{\Lambda}_{t})_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}\bm{C}_{\rm
N}\bm{\Gamma}^{*}+(\bm{r}^{*})^{\prime}\times(\bm{\Lambda}_{t})^{*}\bm{C}_{\rm
N}\frac{d}{d\epsilon}\Big{(}\bm{\Gamma}_{\epsilon}\Big{)}\bigg{|}_{\epsilon=0}\\\\[5.69046pt]
\implies\big{(}\bm{r}^{\prime}\times\bm{n}\big{)}&=\widehat{(\bm{r}^{*})^{\prime}}(\bm{\Lambda}_{t}^{*})\bm{C}_{\rm
N}\bm{\Gamma^{*}}+\Delta\bm{w}^{\prime}\times\big{(}\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*}\big{)}+(\bm{r}^{*})^{\prime}\times\big{(}\widehat{\Delta\bm{\psi}}\big{)}\bm{\Lambda}_{t}^{\rm*}\bm{C}_{\rm
N}\bm{\Gamma^{*}}\\\\[5.69046pt]
&\qquad+(\bm{r}^{*})^{\prime}\times(\bm{\Lambda}_{t})^{*}\bm{C}_{\rm
N}\Big{[}\bm{\Lambda}_{t}^{\rm*,T}\widehat{(\bm{r}^{*})^{\prime}}(\Delta\bm{\psi})+(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\Delta\bm{w}^{\prime}\Big{]}\\\\[5.69046pt]
&=\widehat{(\bm{r}^{*})^{\prime}}\bm{n}^{*}-\big{(}\widehat{\bm{n}^{*}}\big{)}\Delta\bm{w}^{\prime}+\widehat{(\bm{r}^{*})^{\prime}}\big{(}\bm{n}-\bm{n}^{*}\big{)}\end{split}$
(A.6)

Here, $\big{(}\cdot\big{)}^{\rm*}$ denotes the known fields from the
previously converged Newton-Raphson iteration. The manipulation of skew-
symmetric tensors and the corresponding axial vector relation,
$\bm{\theta}\times\bm{h}=\hat{\bm{\theta}}\bm{h}=-(\bm{h}\times\bm{\theta})=-\hat{\bm{h}}\bm{\theta}$
is used in all the above expressions.

##  Appendix B Coefficients of the block linear system

###  B.1 Force Equilibrium

The linearised equation for spatial force $\bm{n}$ (Eq. 18) can be re-written
as,

$\bm{n}=\underbrace{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*}}_{\bm{C}_{\rm exp\rm
w}}+\underbrace{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm T})}_{\bm{C}_{\rm w\rm
w}}\Delta\bm{w}^{\prime}+\underbrace{\left[\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm
T}\widehat{(\bm{r}^{*})^{\prime}}-\widehat{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*})}\right]}_{\bm{C}_{\rm w\bm{\psi}}}\Delta\bm{\psi}$ (B.1)

Using the discretisation schemes for cell faces (Eq. 25) and cell-face
derivatives (Eq. 27), the discretised force equilibrium (Eq. 23) can be
expanded and written in matrix-form as,

$\big{[}\bm{A}_{\rm W}\big{]}_{\bm{n}}\begin{bmatrix}(\Delta\bm{w})_{\rm W}\\\
(\Delta\bm{\psi})_{\rm W}\end{bmatrix}+\big{[}\bm{A}_{\rm
C}\big{]}_{\bm{n}}\begin{bmatrix}(\Delta\bm{w})_{\rm C}\\\
(\Delta\bm{\psi})_{\rm C}\end{bmatrix}+\big{[}\bm{A}_{\rm
E}\big{]}_{\bm{n}}\begin{bmatrix}(\Delta\bm{w})_{\rm E}\\\
(\Delta\bm{\psi}){\rm E}\end{bmatrix}=(\bm{R^{\bm{w}}})_{\rm C}$ (B.2)

where the coefficient matrices $\big{[}\bm{A}_{\rm W}\big{]}_{\bm{n}}$,
$\big{[}\bm{A}_{\rm C}\big{]}_{\bm{n}}$ and $\big{[}\bm{A}_{\rm
E}\big{]}_{\bm{n}}$ are $3\times 6$ matrices containing the contributions of
$\rm W$, $\rm C$ and $\rm E$ cell-centres respectively are given by,

$\big{[}\bm{A}_{\rm W}\big{]}_{\bm{n}}=\begin{bmatrix}\frac{1}{L_{\rm
C}}\Big{(}\gamma_{\rm w}(\bm{C}_{\rm w\rm w})_{\rm W}+(1-\gamma_{\rm
w})(\bm{C}_{\rm w\rm w})_{\rm C}\Big{)}&\qquad-\gamma_{\rm
w}\Big{(}\gamma_{\rm w}(\bm{C}_{\rm w\bm{\psi}})_{\rm W}+(1-\gamma_{\rm
w})(\bm{C}_{\rm w\bm{\psi}})_{\rm C}\Big{)}\end{bmatrix}$ (B.3)

$\big{[}\bm{A}_{\rm C}\big{]}_{\bm{n}}=\begin{bmatrix}-\frac{1}{L_{\rm
C}}\Big{(}\gamma_{\rm e}(\bm{C}_{\rm w\rm w})_{\rm E}+(1-\gamma_{\rm
e})(\bm{C}_{\rm w\rm w})_{\rm C}+\gamma_{\rm w}(\bm{C}_{\rm w\rm w})_{\rm
W}+(1-\gamma_{\rm w})(\bm{C}_{\rm w\rm w})_{\rm C}\Big{)}\\\\[14.22636pt]
(1-\gamma_{\rm e})\Big{(}\gamma_{\rm e}(\bm{C}_{\rm w\bm{\psi}})_{\rm
E}+(1-\gamma_{\rm e})(\bm{C}_{\rm w\bm{\psi}})_{\rm C}\Big{)}-(1-\gamma_{\rm
w})\Big{(}\gamma_{\rm w}(\bm{C}_{\rm w\bm{\psi}})_{\rm W}+(1-\gamma_{\rm
w})(\bm{C}_{\rm w\bm{\psi}})_{\rm C}\Big{)}\end{bmatrix}^{T}$ (B.4)

$\big{[}\bm{A}_{\rm E}\big{]}_{\bm{n}}=\begin{bmatrix}\frac{1}{L_{\rm
C}}\Big{(}\gamma_{\rm e}(\bm{C}_{\rm w\rm w})_{\rm E}+(1-\gamma_{\rm
e})(\bm{C}_{\rm w\rm w})_{\rm C}\Big{)}&\qquad\gamma_{\rm e}\Big{(}\gamma_{\rm
e}(\bm{C}_{\rm w\bm{\psi}})_{\rm E}+(1-\gamma_{\rm e})(\bm{C}_{\rm
w\bm{\psi}})_{\rm C}\Big{)}\end{bmatrix}$ (B.5)

The residual forces evaluated at the cell centre $\rm C$ is given by,

$(\bm{R}^{\bm{w}})_{\rm C}=\gamma_{\rm w}(\bm{C}_{\rm exp\rm w})_{\rm
W}+(1-\gamma_{\rm w})(\bm{C}_{\rm exp\rm w})_{\rm C}-\gamma_{\rm
e}(\bm{C}_{\rm exp\rm w})_{\rm E}-(1-\gamma_{\rm e})(\bm{C}_{\rm exp\rm
w})_{\rm C}-\bm{f}_{\rm C}L_{\rm C}$ (B.6)

###  B.2 Moment Equilibrium

The discretised moment equilibrium Eq. 24 can be written as,

$\bm{m}=\underbrace{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm M}\bm{K}^{*}}_{\bm{C}_{\rm
exp\rm m}}+\underbrace{-\widehat{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}\bm{K}^{*})}}_{\bm{C}_{\rm
m\rm\bm{\psi}}}\Delta\bm{\psi}+\underbrace{(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
M}(\bm{\Lambda}_{t}^{\rm*})^{\rm T})}_{\bm{C}_{\rm m\rm\bm{\psi}\rm
2}}\Delta\bm{\psi}^{\prime}$ (B.7)

Using Eq. 18, the terms in the linearised counterpart of the term
$(\bm{r}^{\prime}\times\bm{n})$ in Eq. (20) can be rearranged and expressed
as,

$\begin{split}(\bm{r}^{\prime}\times\bm{n})=\underbrace{\widehat{(\bm{r}^{*})^{\prime}}(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*})}_{\bm{C}_{\rm exp\rm m\rm
w}}+\underbrace{\Big{[}\widehat{(\bm{r}^{*})^{\prime}}(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm T})-(\widehat{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*}})\Big{]}}_{\bm{C}_{\rm m\rm w}}\Delta\bm{w}^{\prime}\\\
+\underbrace{\Big{[}\widehat{(\bm{r}^{*})^{\prime}}(\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}(\bm{\Lambda}_{t}^{\rm*})^{\rm
T})\widehat{(\bm{r}^{*})^{\prime}}-\widehat{(\bm{r}^{*})^{\prime}}(\widehat{\bm{\Lambda}_{t}^{*}\bm{C}_{\rm
N}\bm{\Gamma}^{*}})\Big{]}}_{\bm{C}_{\rm m\rm
w\rm\bm{\psi}}}\Delta\bm{\psi}\end{split}$ (B.8)

The discretised moment equilibrium Eq. 24 thus can be expanded and written as,

$\big{[}\bm{A}_{\rm W}\big{]}_{\bm{m}}\begin{bmatrix}(\Delta\bm{w})_{\rm W}\\\
(\Delta\bm{\psi})_{\rm W}\end{bmatrix}+\big{[}\bm{A}_{\rm
C}\big{]}_{\bm{m}}\begin{bmatrix}(\Delta\bm{w})_{\rm C}\\\
(\Delta\bm{\psi})_{\rm C}\end{bmatrix}+\big{[}\bm{A}_{\rm
E}\big{]}_{\bm{m}}\begin{bmatrix}(\Delta\bm{w})_{\rm E}\\\
(\Delta\bm{\psi})_{\rm
E}\end{bmatrix}=\begin{bmatrix}(\bm{R}^{\bm{\bm{\psi}}})_{\rm C}\end{bmatrix}$
(B.9)

The coefficient matrices $\big{[}\bm{A}_{\rm W}\big{]}_{\bm{m}}$,
$\big{[}\bm{A}_{\rm C}\big{]}_{\bm{m}}$ and $\big{[}\bm{A}_{\rm
E}\big{]}_{\bm{m}}$ are $3\times 6$ matrices given as,

$\big{[}\bm{A}_{\rm W}\big{]}_{\bm{m}}=\begin{bmatrix}-\frac{1}{2}\frac{L_{\rm
w}}{L_{\rm C}}\Big{(}\gamma_{\rm w}(\bm{C}_{\rm m\rm w})_{\rm
W}+(1-\gamma_{\rm w})(\bm{C}_{\rm m\rm w})_{\rm C}\Big{)}\\\\[22.76228pt]
\begin{gathered}-\gamma_{\rm w}\Big{(}\gamma_{\rm w}(\bm{C}_{\rm
m\bm{\psi}})_{\rm W}+(1-\gamma_{\rm w})(\bm{C}_{\rm m\bm{\psi}})_{\rm
C}\Big{)}+\frac{1}{L_{\rm C}}\Big{(}\gamma_{\rm w}(\bm{C}_{\rm
m\bm{\psi}2})_{\rm W}+(1-\gamma_{\rm w})(\bm{C}_{\rm m\bm{\psi}2})_{\rm
C}\Big{)}\\\ +\frac{1}{2}L_{\rm w}\gamma_{\rm w}\Big{(}\gamma_{\rm
w}(\bm{C}_{\rm m\rm w\bm{\psi}})_{\rm W}+(1-\gamma_{\rm w})(\bm{C}_{\rm m\rm
w\bm{\psi}})_{\rm C}\Big{)}\end{gathered}\end{bmatrix}^{T}$ (B.10)

$\big{[}\bm{A}_{\rm C}\big{]}_{\bm{m}}=\begin{bmatrix}\frac{1}{2}\frac{L_{\rm
w}}{L_{\rm C}}\Big{(}\gamma_{\rm w}(\bm{C}_{\rm m\rm w})_{\rm
W}+(1-\gamma_{\rm w})(\bm{C}_{\rm m\rm w})_{\rm
C}\Big{)}-\frac{1}{2}\frac{L_{\rm e}}{L_{\rm C}}\Big{(}\gamma_{\rm
e}(\bm{C}_{\rm m\rm w})_{\rm E}-(1-\gamma_{\rm e})(\bm{C}_{\rm m\rm w})_{\rm
C}\Big{)}\\\\[22.76228pt] \begin{gathered}(1-\gamma_{\rm e})\Big{(}\gamma_{\rm
e}(\bm{C}_{\rm m\bm{\psi}})_{\rm E}+(1-\gamma_{\rm e})(\bm{C}_{\rm
m\bm{\psi}})_{\rm C}\Big{)}-\frac{1}{L_{\rm C}}\Big{(}\gamma_{\rm
e}(\bm{C}_{\rm m\bm{\psi}2})_{\rm E}+(1-\gamma_{\rm e})(\bm{C}_{\rm
m\bm{\psi}2})_{\rm C}\Big{)}\\\ -(1-\gamma_{\rm w})\Big{(}\gamma_{\rm
w}(\bm{C}_{\rm m\bm{\psi}})_{\rm W}+(1-\gamma_{\rm w})(\bm{C}_{\rm
m\bm{\psi}})_{\rm C}\Big{)}-\frac{1}{L_{\rm C}}\Big{(}\gamma_{\rm
w}(\bm{C}_{\rm m\bm{\psi}2})_{\rm W}+(1-\gamma_{\rm w})(\bm{C}_{\rm
m\bm{\psi}2})_{\rm C}\Big{)}\\\ +\frac{1}{2}L_{\rm e}(1-\gamma_{\rm
e})\Big{(}\gamma_{\rm e}(\bm{C}_{\rm m\rm w\bm{\psi}})_{\rm E}+(1-\gamma_{\rm
e})(\bm{C}_{\rm m\rm w\bm{\psi}})_{\rm C}\Big{)}\\\ +\frac{1}{2}L_{\rm
w}(1-\gamma_{\rm w})\Big{(}\gamma_{\rm w}(\bm{C}_{\rm m\rm w\bm{\psi}})_{\rm
W}+(1-\gamma_{\rm w})(\bm{C}_{\rm m\rm w\bm{\psi}})_{\rm
C}\Big{)}\end{gathered}\end{bmatrix}^{T}$ (B.11)

$\big{[}\bm{A}_{\rm E}\big{]}_{\bm{m}}=\begin{bmatrix}\frac{1}{2}\frac{L_{\rm
e}}{L_{\rm C}}\Big{(}\gamma_{\rm e}(\bm{C}_{\rm m\rm w})_{\rm
E}+(1-\gamma_{\rm e})(\bm{C}_{\rm m\rm w})_{\rm C}\Big{)}\\\\[22.76228pt]
\begin{gathered}\gamma_{\rm e}\Big{(}\gamma_{\rm e}(\bm{C}_{\rm
m\bm{\psi}})_{\rm E}+(1-\gamma_{\rm e})(\bm{C}_{\rm m\bm{\psi}})_{\rm
C}\Big{)}+\frac{1}{L_{\rm C}}\Big{(}\gamma_{\rm e}(\bm{C}_{\rm
m\bm{\psi}2})_{\rm E}+(1-\gamma_{\rm e})(\bm{C}_{\rm m\bm{\psi}2})_{\rm
C}\Big{)}\\\ +\frac{1}{2}L_{\rm e}\gamma_{\rm e}\Big{(}\gamma_{\rm
e}(\bm{C}_{\rm m\rm w\bm{\psi}})_{\rm E}+(1-\gamma_{\rm e})(\bm{C}_{\rm m\rm
w\bm{\psi}})_{\rm C}\Big{)}\end{gathered}\end{bmatrix}^{T}$ (B.12)

The residual moment evaluated at the cell centre $\rm C$ is given by,

$\begin{gathered}(\bm{R}^{\bm{\psi}})_{\rm C}\quad=\quad\gamma_{\rm
w}(\bm{C}_{\rm exp\rm m})_{\rm W}+(1-\gamma_{\rm w})(\bm{C}_{\rm exp\rm
m})_{\rm C}-\gamma_{\rm e}(\bm{C}_{\rm exp\rm m})_{\rm E}-(1-\gamma_{\rm
e})(\bm{C}_{\rm exp\rm m})_{\rm C}\\\ -\frac{1}{2}L_{\rm w}\Big{(}\gamma_{\rm
w}(\bm{C}_{\rm exp\rm m\rm w})_{\rm W}+(1-\gamma_{\rm w})(\bm{C}_{\rm exp\rm
m\rm w})_{\rm C}\Big{)}\\\ -\frac{1}{2}L_{\rm e}\Big{(}\gamma_{\rm
e}(\bm{C}_{\rm exp\rm m\rm w})_{\rm E}+(1-\gamma_{\rm e})(\bm{C}_{\rm exp\rm
m\rm w})_{\rm C}\Big{)}-\bm{t}_{\rm C}L_{\rm C}\end{gathered}$ (B.13)

The diagonal and off-diagonal coefficient matrix contributions from the force
and moment equilibrium equations can be assembled and the coupled block linear
system (Eq. 29) coefficient matrices take the form as follows,

$\bm{A}_{\rm W}=\begin{bmatrix}\big{[}\bm{A}_{\rm W}\big{]}_{\bm{n}}\\\ \\\
\big{[}\bm{A}_{\rm W}\big{]}_{\bm{m}}\end{bmatrix}\quad;\quad\bm{A}_{\rm
C}=\begin{bmatrix}\big{[}\bm{A}_{\rm C}\big{]}_{\bm{n}}\\\ \\\
\big{[}\bm{A}_{\rm C}\big{]}_{\bm{m}}\end{bmatrix}\quad;\quad\bm{A}_{\rm
E}=\begin{bmatrix}\big{[}\bm{A}_{\rm E}\big{]}_{\bm{n}}\\\ \\\
\big{[}\bm{A}_{\rm E}\big{]}_{\bm{m}}\end{bmatrix}$ (B.14)

## References

  * Argyris [1982] Argyris, J., 1982.  An excursion into large rotations.  Computer methods in applied mechanics and engineering 32, 85–155.
  * Bathe and Bolourchi [1979] Bathe, K.J., Bolourchi, S., 1979\.  Large displacement analysis of three-dimensional beam structures.  International journal for numerical methods in engineering 14, 961–986.
  * Cardiff and Demirdžić [2021] Cardiff, P., Demirdžić, I., 2021\.  Thirty years of the finite volume method for solid mechanics.  Archives of Computational Methods in Engineering , 1–60.
  * Cardona and Geradin [1988] Cardona, A., Geradin, M., 1988\.  A beam finite element non-linear theory with finite rotations.  International journal for numerical methods in engineering 26, 2403–2438.
  * Chadha and Todd [2020] Chadha, M., Todd, M.D., 2020\.  The mathematical theory of a higher-order geometrically-exact beam with a deforming cross-section.  International Journal of Solids and Structures 202, 854–880.
  * Crisfield and Jelenic [1999] Crisfield, M., Jelenic, G., 1999\.  Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 455, 1125–1147.
  * Crisfield [1990] Crisfield, M.A., 1990.  A consistent co-rotational formulation for non-linear, three-dimensional, beam-elements.  Computer methods in applied mechanics and engineering 81, 131–150.
  * Dadeppo and Schmidt [1975] Dadeppo, D., Schmidt, R., 1975\.  Instability of clamped-hinged circular arches subjected to a point load.  Journal of Applied Mechanics 42, 894\.
  * Fallah [2004] Fallah, N., 2004.  A cell vertex and cell centred finite volume method for plate bending analysis.  Computer methods in applied Mechanics and Engineering 193, 3457–3470.
  * Fallah [2005] Fallah, N., 2005.  A new approach in cell centred finite volume formulation for plate bending analysis, in: Int. Conf. Comput. Methods Sci. and Eng.(ICCMSE 2005), Lect. Ser. Comput. Comput. Sci., pp. 187–190.
  * Fallah [2013] Fallah, N., 2013.  Finite volume method for determining the natural characteristics of structures.  Journal of Engineering Science and Technology 8, 93–106.
  * Fallah and Ebrahimnejad [2014] Fallah, N., Ebrahimnejad, M., 2014\.  Finite volume analysis of adaptive beams with piezoelectric sensors and actuators.  Applied Mathematical Modelling 38, 722–737.
  * Fallah and Ghanbari [2018] Fallah, N., Ghanbari, A., 2018\.  A displacement finite volume formulation for the static and dynamic analysis of shear deformable circular curved beams.  Scientia Iranica 25, 999–1014.
  * Fallah and Hatami [2006a] Fallah, N., Hatami, F., 2006a.  A displacement formulation based on finite volume method for analysis of timoshenko beam, in: Proceedings of the 7th international conference on civil engineering, Tehran, Iran.
  * Fallah and Hatami [2006b] Fallah, N., Hatami, F., 2006b.  Extension of the finite volume method for instability analysis of columns with shear effects, in: Proceedings of the eighth international conference on computational structures technology, Stirlingshire, Scotland. Civil-Comp Press. Paper.
  * Fallah [2006] Fallah, N.A., 2006.  A finite volume method for plate buckling analysis, in: III European Conference on Computational Mechanics, Springer. pp. 690–690.
  * Geradin and Cardona [1988] Geradin, M., Cardona, A., 1988\.  Kinematics and dynamics of rigid and flexible mechanisms using finite elements and quaternion algebra.  Computational Mechanics 4, 115–135.
  * Ghosh and Roy [2008] Ghosh, S., Roy, D., 2008.  Consistent quaternion interpolation for objective finite element approximation of geometrically exact beam.  Computer Methods in Applied Mechanics and Engineering 198, 555–571.
  * Ghosh and Roy [2009] Ghosh, S., Roy, D., 2009.  A frame-invariant scheme for the geometrically exact beam using rotation vector parametrization.  Computational Mechanics 44, 103–118.
  * Ibrahimbegovic [1997] Ibrahimbegovic, A., 1997.  On the choice of finite rotation parameters.  Computer Methods in Applied Mechanics and Engineering 149, 49–71.
  * Ibrahimbegović et al. [1995] Ibrahimbegović, A., Frey, F., Kožar, I., 1995.  Computational aspects of vector-like parametrization of three-dimensional finite rotations.  International Journal for Numerical Methods in Engineering 38, 3653–3673.
  * Isić et al. [2007a] Isić, S., Doleček, V., Karabegović, I., 2007a.  A comparison between finite element and finite volume methods on the problem of stability of timoshenko beam, in: The 12th international conference on problems of material engineering, Jasna, Slovakia.
  * Isić et al. [2007b] Isić, S., Doleček, V., Karabegović, I., 2007b.  Numerical and experimental analysis of postbuckling behaviour of prismatic beam under displacement dependent loading, in: Proceedings of First Serbian Congress on Theoretical and Applied Mechanics.-Kopaonik, Serbia, pp. 331–338.
  * Jelenić and Crisfield [1999] Jelenić, G., Crisfield, M., 1999\.  Geometrically exact 3d beam theory: implementation of a strain-invariant finite element for statics and dynamics.  Computer methods in applied mechanics and engineering 171, 141–171.
  * Jelenić and Saje [1995] Jelenić, G., Saje, M., 1995\.  A kinematically exact space finite strain beam model—finite element formulation by generalized virtual work principle.  Computer Methods in Applied Mechanics and Engineering 120, 131–161.
  * Jing et al. [2016] Jing, L.L., Ming, P.J., Zhang, W.P., Fu, L.R., Cao, Y.P., 2016.  Static and free vibration analysis of functionally graded beams by combination timoshenko theory and finite volume method.  Composite structures 138, 192–213.
  * Kapania and Li [2003] Kapania, R., Li, J., 2003.  A formulation and implementation of geometrically exact curved beam elements incorporating finite strains and finite rotations.  Computational Mechanics 30, 444–459.
  * Mäkinen [2007] Mäkinen, J., 2007.  Total lagrangian reissner’s geometrically exact beam element without singularities.  International journal for numerical methods in engineering 70, 1009–1048.
  * Marino [2016] Marino, E., 2016.  Isogeometric collocation for three-dimensional geometrically exact shear-deformable beams.  Computer Methods in Applied Mechanics and Engineering 307, 383–410.
  * Meier et al. [2014] Meier, C., Popp, A., Wall, W.A., 2014.  An objective 3d large deformation finite element formulation for geometrically exact curved kirchhoff rods.  Computer Methods in Applied Mechanics and Engineering 278, 445–478.
  * Meier et al. [2019] Meier, C., Popp, A., Wall, W.A., 2019.  Geometrically exact finite element formulations for slender beams: Kirchhoff–love theory versus simo–reissner theory.  Archives of Computational Methods in Engineering 26, 163–243.
  * Moosavi et al. [2011] Moosavi, M., Delfanian, F., Khelil, A., 2011.  The orthogonal meshless finite volume method for solving euler–bernoulli beam and thin plate problems.  Thin-walled structures 49, 923–932.
  * Noor and Peters [1981] Noor, A.K., Peters, J.M., 1981\.  Mixed models and reduced/selective integration displacement models for nonlinear analysis of curved beams.  International Journal for Numerical Methods in Engineering 17, 615–631.
  * Romero [2004] Romero, I., 2004.  The interpolation of rotations and its application to finite element models of geometrically exact rods.  Computational mechanics 34, 121–133.
  * Romero and Armero [2002] Romero, I., Armero, F., 2002\.  An objective finite element approximation of the kinematics of geometrically exact rods and its use in the formulation of an energy–momentum conserving scheme in dynamics.  International Journal for Numerical Methods in Engineering 54, 1683–1716.
  * Simo [1985] Simo, J., 1985.  A finite strain beam formulation. the three-dimensional dynamic problem. part i.  Computer Methods in Applied Mechanics and Engineering 49, 55–70.
  * Simo and Vu-Quoc [1986] Simo, J.C., Vu-Quoc, L., 1986\.  A three-dimensional finite-strain rod model. part ii: Computational aspects.  Computer methods in applied mechanics and engineering 58, 79–116.
  * Tukovic et al. [2019] Tukovic, Ž., De Jaeger, P., Cardiff, P., Ivankovic, A., 2019\.  A finite volume solver for geometrically exact simo-reissner beams, in: ECCOMAS MSF Thematic conference. Sarajevo, Bosnia and Herzegovina.
  * Weller et al. [1998] Weller, H.G., Tabor, G., Jasak, H., Fureby, C., 1998\.  A tensorial approach to computational continuum mechanics using object-oriented techniques.  Computers in physics 12, 620–631.
  * Wheel [1997] Wheel, M., 1997.  A finite volume method for analysing the bending deformation of thick and thin plates.  Computer methods in applied mechanics and engineering 147, 199–208.
  * Zupan and Saje [2003] Zupan, D., Saje, M., 2003.  Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures.  Computer Methods in Applied Mechanics and Engineering 192, 5209–5248.
  * Zupan et al. [2009] Zupan, E., Saje, M., Zupan, D., 2009.  The quaternion-based three-dimensional beam theory.  Computer methods in applied mechanics and engineering 198, 3944–3956.