# A finite element/quaternion/asymptotic numerical method for the 3D simulation of flexible cables

Emmanuel Cottanceau, Olivier Thomas, Philippe Véron, Marc Alochet, Renaud Deligny

► To cite this version:

Emmanuel Cottanceau, Olivier Thomas, Philippe Véron, Marc Alochet, Renaud Deligny. A finite element/quaternion/asymptotic numerical method for the 3D simulation of flexible cables. Finite Elements in Analysis and Design, 2018, 139, pp.14-34. 10.1016/j.finel.2017.10.002. hal-01989778

HAL Id: hal-01989778

https://hal.science/hal-01989778

Submitted on 22 Jan 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A finite element/quaternion/asymptotic numerical method for the 3D simulation of flexible cables

Emmanuel Cottanceau<sup>a, b</sup>, Olivier Thomas<sup>a, *</sup>, Philippe Véron<sup>a</sup>, Marc Alochet<sup>b</sup>, Renaud Deligny<sup>b</sup>

<sup>a</sup>Arts et Métiers ParisTech, LSIS UMR CNRS 7296, 8 bd. Louis XIV, 59046 Lille, France

<sup>b</sup>Technocentre Renault, 1 av. Du Golf, 78084 Guyancourt, France

## ARTICLE INFO

**Keywords:**

Cables

Finite element method

Asymptotic numerical method

Quaternions

## ABSTRACT

In this paper, a method for the quasi-static simulation of flexible cables assembly in the context of automotive industry is presented. The cables geometry and behavior encourage to employ a geometrically exact beam model. The 3D kinematics is then based on the position of the centerline and on the orientation of the cross-sections, which is here represented by rotational quaternions. Their algebraic nature leads to a polynomial form of equilibrium equations. The continuous equations obtained are then discretized by the finite element method and easily recast under quadratic form by introducing additional slave variables. The asymptotic numerical method, a powerful solver for systems of quadratic equations, is then employed for the continuation of the branches of solution. The originality of this paper stands in the combination of all these methods which leads to a fast and accurate tool for the assembly process of cables. This is proved by running several classical validation tests and an industry-like example.

## 1. Introduction

During the last decades, the room available in car vehicles (e.g. in engine compartment) has plummeted because of the rapid development of on-board electronics. As a result, a need for very accurate numerical tools for design has appeared in automotive industry. In the meantime, a fast computation is necessary so that design duration remains suitable for industry. In this context, flexible pieces represent an outstanding challenge since, unlike most of car pieces, they cannot be modeled as rigid body solids in CAD software. This paper focuses on a specific type of flexible piece, namely electrical cables. Cables have a complex structure. A wire is made up of copper filaments wrapped in an elastomer duct. These wires are most of the time gathered in bundles which are themselves surrounded by various protections such as tape, PVC tube ... Moreover, the full cable is often constituted of several drifted cable pieces forming a system with a complex geometry.

Due to its slender shape and its flexibility, one can consider a simple cable as a beam undergoing large displacements and large rotations. There exist several theories accounting for nonlinear beams, see Refs. [1-3]. The most widely used theory is the geometrically exact beam model whose founding principles were established by Reissner [4,5] and further generalized by Simo [6]. Various finite element formula-

tions (FEM) have been presented to solve these equations numerically. The notable works of [7-11] and more recently [12,13] can, among others, be listed.

At the heart of all these formulations, the rotation parameterization is of paramount importance in the numerical models. The 3D rotations modeling indeed is not an easy task especially when computational efficiency is sought. The rotational vector-like parameterization used by many authors features only 3 parameters (minimum set in 3D). However, as no parameterization of less than 4 parameters can be singularity-free [14], this choice poses several numerical limitations and lacks robustness. A powerful alternative consists in using quaternions, a set of 4 singularity-free parameters. Firstly used only for storage in the numerical models, Zupan et al. [12] have recently shown their utility when used as primary variables. They also have developed a model without rotation matrices exploiting the high potential of quaternion algebra, and very efficient for numerical purposes [15]. In addition, quaternions offer an original description of rotations since they substitute the usual trigonometric functions by algebraic variables and lead to polynomial equilibrium equations.

The finite element method, very adapted to the assembly of complex geometries such as electrical cables ones, is applied to the continuous equilibrium equations. The algebraic system obtained is then

<footnote>* Corresponding author.
E-mail address: olivier.thomas@ensam.eu (O. Thomas).</footnote>
generally solved by a classical predictor-corrector method (PCM) [16]. Even if the appearance of arc-length methods [17] has considerably enhanced the robustness of this type of method, they often require to choose a step size which may be very tricky for the user. A sufficiently small step size allows to compute even highly nonlinear part of equilibrium branches but in return may impractically increase the computation time, while a larger step size may spoil the convergence. Taking advantage of the polynomial form of the system of equations obtained when using quaternion parameters an alternative consists in replacing the PCM by the asymptotic numerical method (ANM) firstly presented by Damil and Potier-Ferry [18] and by Cochelin [19]. The ANM is a very powerful solver for quadratic problems and it overcomes all the drawbacks of the PCM. This technique indeed is very robust, does not require any tuning parameters and is thus well suited for an industrial use. In addition, Cochelin and Medale [20] have equipped the method with a bifurcation detector and improved its efficiency in the vicinity of bifurcation points.

Combining quaternions with the ANM has already been set up on a rod model discretized with a finite difference scheme by Lazarus et al. [21], which have got very promising results. We propose here to set up the technique on the finite-element based geometrically exact beam model, in what constitutes the main originality of this paper. Validations and illustrations of the method on very intricate problems are provided and discussed. A critical evaluation and future researches are presented by way of conclusion.

## 2. Governing equations

In this section, the classical quasi-static formulation of the geometrically exact beam model, based on the rotation vector, is firstly presented. It enables to explain all the main ingredients of the model and to discuss their physical meaning. Secondly, the equations are modified by using quaternions instead of the rotation vector. This leads to the formulation which serves our numerical model, presented in part 3.

### 2.1. The geometrically exact beam model

In the geometrically exact beam model, a beam is described by defining a family of cross-sections whose centroids form a curve called

the centerline of the beam. The kinematic variables essential to this description are introduced following the notations used in Ref. [22]. Let us consider a beam of initial length L, with an arbitrary cross-section Ω, as depicted Fig. 1. Let $(u_1, u_2, u_3)$ be a fixed Cartesian (global) frame. The current position of the centerline in this frame is described by the vector field $x_0(s)$ which is a function of the curvilinear abscissa $s \in [0, L]$ along the beam axis. A material (local) frame $(e_1(s), e_2(s), e_3(s))$ is introduced to define the cross-section at abscissa $s$. The vectors $e_2$ and $e_3$ span the cross-section while the vector $e_1$ remains normal to the cross-section for every deformed configurations. It is essential to note that $e_1$ is not necessarily tangent to the centerline of the beam such that shear deformation is taken into account (Timoshenko model). The reference configuration is chosen such that the beam is unstressed. However, in this state, the reference position of the centerline is an arbitrary curve and not necessarily a straight line, thus accounting for the initial curvature of the beam. Its initial position is then a function of $s$ denoted $X_0(s)$, such that the displacement of any point of the centerline is $x_0(s) - X_0(s)$. Similarly, the material frame in the reference configuration depends on $s$ and is denoted $(E_1(s), E_2(s), E_3(s))$. To end up with the definition of the kinematic variables, let us introduce the two rotation operators $R_0(s)$ and $R(s)$ which depict the orientation of the material frames in the global frame in the initial and in the current configuration respectively, that being $E_i(s) = R_0(s)u_i$ and $e_i(s) = R(s)u_i$, $i \in \{1, 2, 3\}$. As $R_0(s)$ defines the initial curvature of the beam it is constant through the deformation. With these notations, the position vector in the spatial frame of any point $M'$ of the undeformed beam located in the section at $s$ writes

$$X(s, X_2, X_3) = X_0(s) + R_0(s)Y(X_2, X_3), \quad (1)$$

where $Y(X_2, X_3) = [0 \ X_2 \ X_3]^T$ is the material position of $M'$ in the cross-section. Under the assumption that cross-sections remain plane and do not undergo any deformations along the transformation, $Y(X_2, X_3)$ is constant and $M'$ after deformation becomes $M''$ whose position is given by

$$x(s, X_2, X_3) = x_0(s) + R(s)Y(X_2, X_3). \quad (2)$$

The current configuration of the beam is then completely characterized by the position of the centerline $x_0(s)$ and the orientation of the cross-sections $R(s)$. One recovers here that the kinematic variables depend solely on the curvilinear abscissa $s$ as for any beam model. As it is

<figure><img src="image_2.png" /><figcaption>Fig. 1. Kinematics of the geometrically exact beam model.</figcaption></figure>
clearly implied, the dependency in s is dropped in the following. Fur-
thermore, the derivative with respect to s is from now on denoted with
the prime symbol ′.

For any rotation operator **R**, it is shown that the matrix **R**^T**R**' is skew-symmetric [22]. Any 3 × 3 skew-symmetric matrix **A** is made of only 3 independent components which can be gathered in a vector **a** = **vec**(**A**) = [a₁ a₂ a₃]ᵀ. The skew-symmetric matrix associated to **a** is denoted **ā** and is equal to:

$$
\tilde{\mathbf{a}} = \begin{bmatrix} 0 & -a_3 & a_2 \\ a_3 & 0 & -a_1 \\ -a_2 & a_1 & 0 \end{bmatrix}. \tag{3}
$$

With these notations and hypothesis, the material translational and
rotational strain measures are equal to [23].

$$
\Gamma = [\Gamma_1 \ \Gamma_2 \ \Gamma_3]^T = \mathbf{R}^T \mathbf{x}'_0 - \Gamma_0 \quad (4a)
$$

$$
\mathbf{K} = [\kappa_1 \ \kappa_2 \ \kappa_3]^T = \mathrm{vec}(\mathbf{R}^T \mathbf{R}') - \mathbf{K}_0, \quad (4\mathrm{b})
$$

where the intrinsic translational strain Γ₀ = **R**'₀**X**'₀ = **u**₁ and the intrinsic
rotational strain **K**₀ = **vec**(R̂'₀**R**'₀) solely depend on the initial configu-
ration of the beam. Let us point out that in the reference configuration,
the strain measures Γ and **K** are both equal to **0**: the assumption of an
initially unstressed configuration is thus well recovered by these mea-
sures definition.

To interpret these measures, let us consider the straightforward case
in which the beam centerline is a straight line in the reference configu-
ration. The intrinsic strains then obey Γ₀ = **u**₁ and **K**₀ = **0**. **x**'₀(s) is, by
definition, a vector tangent to the centerline of the beam at abscissa s
in the current configuration. In expression (4a), the operator **R**^T maps
back this vector projection from the material frame to the fixed frame.
The measure Γ then consists in comparing this quantity to its reference
value **u**₁. A fully equivalent interpretation is made by recalling that
the scalar product of 2 vectors may be geometrically interpreted as the
projection of one vector in the direction of the other. It follows that
Γ₁ = **e**'₀**x**'₀ - 1, Γ₂ = **e**'₂**x**'₀ and Γ₃ = **e**'₃**x**'₀ are the comparison of the pro-
jection of **x**'₀(s) on the material frame to **e**₁ in the 3 material directions.
As the beam is parameterized by the curvilinear abscissa in the initial
configuration, the norm of **x**'₀(s) is equal to 1 if the length of the center-
line remains unchanged, is greater than 1 in case of elongation of the
centerline and smaller than 1 in case of shortening. With regard to all
these observations, we eventually infer that Γ₁ is the axial strain of the
neutral axis, while Γ₂ and Γ₃ represent the shear strains in directions
**e**₂ and **e**₃.

In the case of a pure rotational transformation **R**, **R**^T **R**' represents
the gradient of the transformation along the beam axis. As this quan-
tity is a skew-symmetric matrix, it is equivalent to write **R**^T **R**' **Y** and
**vec**(**R**^T **R**') × **Y** = **K** × **Y**. It follows directly that in (4b) κ₁ depicts the
torsional strain while κ₂ and κ₃ represent the curvature strains around
**e**₂ and **e**₃ respectively.

The expression of the aforementioned rotation operator **R** is gener-
ally obtained from a reduced representation. The choice of this rep-
resentation is actually a crucial issue since various criteria compete
against each other to get the best set of parameters: numerical efficiency
(number of parameters), mathematical form, existence of singularities,
geometric interpretation. Let us recall that any 3D rotation of a rigid
body solid is made up of 3 independent parameters. The 3-parameters
rotational vector **ϑ** presented in this paragraph thus constitutes the min-
imal set of parameters. We choose it here for its easy geometrical inter-
pretation since it is defined by

$$
\mathfrak{g} = \mathfrak{g} \mathbf{n} = [\vartheta_1 \quad \vartheta_2 \quad \vartheta_3]^T, \quad (5)
$$

for a rotation of angle $\vartheta$ around the unit axis $\mathbf{n}$. The rotation operator then writes in terms of this parameterization as [22]

$$
\mathbf{R}(\vartheta) = I_3 + \frac{\sin(\vartheta)}{\vartheta} \tilde{\vartheta} + \frac{1 - \cos(\vartheta)}{\vartheta^2} \tilde{\vartheta} \tilde{\vartheta}. \quad (6)
$$

2.2. Virtual work principle

Using the quantities introduced in the previous section and the
Cartesian rotational vector **ϑ** defined in (5) as rotational parameters,
the virtual work principle for a 3D shear elastic beam is stated here.
Following Reissner's work [4,5] and the generalization of Simo [6], it
writes

$$
\int_{0}^{L} (\mathbf{N} \cdot \delta\mathbf{\Gamma} + \mathbf{M} \cdot \delta\mathbf{K}) \, ds = \int_{0}^{L} (\mathbf{n}_{e} \cdot \delta\mathbf{x}_{0} + \mathbf{m}_{e} \cdot \delta\mathbf{\vartheta}) \, ds \\
+ [\mathbf{N}_{e} \cdot \delta\mathbf{x}_{0}]_{0}^{L} + [\mathbf{M}_{e} \cdot \delta\mathbf{\vartheta}]_{0}^{L}, \qquad (7)
$$

with the generalized stress resultant in force and moment **N** and **M**,
the external applied force and moment per unit length **n**<sub>e</sub> and **m**<sub>e</sub>, the
external end force and moment **N**<sub>e</sub> and **M**<sub>e</sub> and the virtual positional
and rotational vectors δ**x**<sub>0</sub> and δ**θ**. The left-hand of the equation is
the opposite of the virtual work of internal forces −δ**ω**<sub>i</sub> and the right-
hand is the virtual work of external forces δ**ω**<sub>e</sub>. **N** = [N T<sub>2</sub> T<sub>3</sub>]<sup>T</sup> and
**M** = [M<sub>t</sub> M<sub>2</sub> M<sub>3</sub>]<sup>T</sup> are work conjugates to the virtual translational strain
δ**Γ** and the virtual rotational strain δ**K** respectively. Therefore, **N** stands
for the axial force resultant, T<sub>2</sub> and T<sub>3</sub> the transverse force resultants
in respective directions **e**<sub>2</sub> and **e**<sub>3</sub>, **M**<sub>t</sub> the torsional moment and M<sub>2</sub> and
M<sub>3</sub> the bending moments around axis **e**<sub>2</sub> and **e**<sub>3</sub>.

In the scope of this paper, only linear elastic materials are consid-
ered. In this particular case, the stress-strain relationship states pro-
portionality between stress resultant vector **Σ** = [**M**^T **N**^T]^T and strain
resultant vector **ε** = [**Γ**^T **K**^T]^T that be

$$
\Sigma = C\mathcal{E} = \operatorname{diag}(C_N, C_M)\mathcal{E}. \tag{8}
$$

In (8), $C_N = \operatorname{diag}(EA, GA_2, GA_3)$ and $C_M = \operatorname{diag}(GJ, EI_2, EI_3)$ are the two
$3 \times 3$ stress-strain submatrices corresponding solely to the force compo-
nents and the moment components respectively. In these expressions,
$E$ and $G$ are respectively the Young's modulus and the shear modulus
of the material; $A$ is the cross-sectional area; $A_2$ and $A_3$ are the shear
areas in directions $e_2$ and $e_3$; $I_2$, $I_3$ are the second moments of area with
respect to axes $e_2$, $e_3$; $J$ is the polar moment of area.

Let us notice that this model can take into account the case of the
origin of the material frame not being on the centerline effortlessly. In
this case, out-of-diagonal terms simply appear in the stress-strain matrix
**C**.

Finally, the virtual strains in (7) are given in terms of the kinematic
variables by the equations [22]:

$$
\begin{align}
\delta \Gamma &= \mathbf{R}^T \delta \mathbf{x}'_0 + \mathbf{R}^T \mathbf{x}'_0 T \delta \mathbf{\vartheta} = \mathbf{R}^T \delta \mathbf{x}'_0 + (\tilde{\Gamma} + \tilde{\Gamma}_0) T \delta \mathbf{\vartheta}, \tag{9} \\
\delta K &= (\mathbf{R}^T \mathbf{R}' T + T') \delta \mathbf{\vartheta} + T \delta \mathbf{\vartheta}' = ((\tilde{\mathbf{K}} + \tilde{\mathbf{K}}_0)T + T') \delta \mathbf{\vartheta} + T \delta \mathbf{\vartheta}'. \nonumber
\end{align}
$$

In (9), **T** is a tangent operator which can be expressed as a function of
**ϑ** only under the form:

$$
T = I_3 + \frac{\cos(\vartheta) - 1}{\vartheta^2} \tilde{\vartheta} + \frac{1}{\vartheta^2} \left(1 - \frac{\sin(\vartheta)}{\vartheta}\right) \tilde{\vartheta} \tilde{\vartheta}. \quad (10)
$$

2.3. Strong form of equilibrium equations

Introducing expressions (9) into the virtual work principle (7), inte-
grating by parts with some rearranging and using the fundamental
lemma of calculus of variations lead to the following equilibrium equa-
tions:

$$
\left\{
\begin{array}{l}
(\mathbf{R}\mathbf{N})' + \mathbf{n}_e = \mathbf{0}, \\
(\mathbf{R}\mathbf{M})' + \mathbf{x}'_0 \times (\mathbf{R}\mathbf{N}) + \mathbf{m}_e = \mathbf{0},
\end{array}
\right.
\quad \forall s \in [0, L] \qquad (11)
$$

along with the natural boundary conditions

$$
\left\{
\begin{array}{ll}
\mathbf{R}\mathbf{N} - \mathbf{N}_e = \mathbf{0} & \text{at } s = 0, L, \\
\mathbf{R}\mathbf{M} - \mathbf{M}_e = \mathbf{0} & \text{at } s = 0, L,
\end{array}
\right.
\qquad (12)
$$
or the essential boundary conditions

$$ \begin{cases} \mathbf{x}_0 = \mathbf{x}_d & \text{at } s = 0, L, \\ \vartheta = \vartheta_d & \text{at } s = 0, L, \end{cases} \qquad (13) $$

where $\mathbf{x}_d$ and $\vartheta_d$ are respectively imposed positions and rotations.

## 2.4. Quaternion parameterization of rotations

As explained at the end of paragraph 2.1, the choice of rotation parameters is of paramount importance in numerical models. In a total Lagrangian framework (i.e. using total rotations), the rotational vector $\vartheta$ prevents from having angles greater than $2\pi$. It is clear from equation (10) that the tangent operator $T$ indeed becomes rank-deficient around $2\pi$-angles and leads to a singularity. Procedures exist to circumvent this flaw [9], however, we find it more expedient to use an other parameterization not requiring any special procedure. As demonstrated in Ref. [14], the minimal set of parameters avoiding all singularities is made up of 4 parameters. With regard to this principle, the 4-parameters quaternions appear as an optimal choice for the representation of rotations. In addition, quaternions feature a special algebra totally equivalent to the algebra of rotation operators but which proves computationally more efficient [24]. Finally, they offer a polynomial representation of rotations in comparison to classical trigonometric ones. We will show in the following part how this property is used as a leverage in our numerical solution.

The next developments of this paragraph aim at giving the fundamentals on quaternions necessary to the comprehension of the method. The interested reader may consult for instance [25] for more details. The notations of [12] on quaternion algebra are used here. A quaternion $\hat{a}$ is defined by

$$ \hat{a} = a_0 + ia_1 + ja_2 + ka_3, \qquad (14) $$

where $a_0, a_1, a_2, a_3 \in \mathbb{R}$ and the imaginary numbers $i, j$ and $k$ are linked by the identities

$$ i^2 = j^2 = k^2 = ijk = -1. \qquad (15) $$

1, $i, j$ and $k$ form a basis of the set of quaternions, denoted $\mathbb{H}$. As a result, a quaternion may also be seen as an element of the four-dimensional Euclidean linear space $\mathbb{R}^4$, whose components are $a_0, a_1, a_2$ and $a_3$ in the Euclidean basis made up of the 4 vectors $[1 \ 0 \ 0 \ 0]^T$, $[0 \ 1 \ 0 \ 0]^T$, $[0 \ 0 \ 1 \ 0]^T$, $[0 \ 0 \ 0 \ 1]^T$. In this paper, no distinctions are made between the quaternions and their projection into a basis of the Euclidean space so that we can smoothly write a quaternion under the vector form $\hat{a} = [a_0 \ a_1 \ a_2 \ a_3]^T$. To distinguish quaternions from the three-dimensional vectors of $\mathbb{R}^3$, the hat symbol is used: $\mathbf{a} \in \mathbb{R}^3$, $\hat{a} \in \mathbb{R}^4$.

Let us introduce the handy notation

$$ \hat{a} = a_0 + \mathbf{a}. \qquad (16) $$

In this expression, $a_0$ and $\mathbf{a} = ia_1 + ja_2 + ka_3$ are called respectively the scalar and the vector part of the quaternion. This notation helps drawing a parallel between quaternions and complex numbers, the scalar part amounting to the real part and the vector part amounting to the imaginary part. It is also helpful when one wants to express the extension of a vector of $\mathbb{R}^3$ in the quaternion space $\mathbb{H}$. Indeed, any basis $(i_1, i_2, i_3)$ of the three-dimensional Euclidean space can be extended into a basis of the four-dimensional Euclidean space $(\hat{i}_0, \hat{i}_1, \hat{i}_2, \hat{i}_3)$, with $\hat{i}_0 = 1 + 0$ and $\hat{i}_k = 0 + i_k$, $k \in \{1, 2, 3\}$. As a result, any vector $\mathbf{v}$ of $\mathbb{R}^3$ has a quaternion counterpart in $\mathbb{H}$, $\hat{\mathbf{v}} = 0 + \mathbf{v}$. This quaternion with a zero scalar part is called a pure quaternion and depicts the extension of a vector into a quaternion. Conversely, a restriction operation may be defined, which transforms a pure quaternion into a vector. This operation is denoted $\text{Vec}(\hat{\mathbf{p}}) = \mathbf{p} \in \mathbb{R}^3$. In the following, no distinctions

will be made between pure quaternions and other quaternions in the notation. For a better comprehension, the reader just has to keep in mind that any vector extension in the quaternion space is a pure quaternion.

Using the analogy with the complex numbers, a conjugate quaternion is defined as $\hat{a}^* = a_0 - \mathbf{a}$. The set of quaternions is equipped with the inner product $\hat{a} \cdot \hat{b} = a_0 b_0 + \mathbf{a} \cdot \mathbf{b}$ and the quaternion norm $\|\hat{a}\| = \sqrt{\hat{a} \cdot \hat{a}}$. Furthermore, the three following elementary operations are also defined: the addition $\hat{a} + \hat{b} = (a_0 + b_0) + (\mathbf{a} + \mathbf{b})$, the scalar multiplication $\lambda \hat{a} = \lambda a_0 + \lambda \mathbf{a}$ for $\lambda \in \mathbb{R}$ and the quaternion multiplication, denoted $\circ$, which follows directly from the definition of a quaternion (14) and from the relations (15):

$$ \hat{a} \circ \hat{b} = a_0 b_0 - \mathbf{a} \cdot \mathbf{b} + (b_0 \mathbf{a} + a_0 \mathbf{b} + \mathbf{a} \times \mathbf{b}), \qquad (17) $$

where $\circ$ and $\times$ are the scalar product and the cross product of vectors. This operation is associative but not commutative, so left and right multiplication must be distinguished.

Finally, it is convenient for the further presented formulation to extend the cross product of vectors of $\mathbb{R}^3$ to pure quaternions. For two pure quaternions $\hat{p}$ and $\hat{q}$, it is defined as

$$ \hat{p} \times \hat{q} = 0 + \mathbf{p} \times \mathbf{q} = \hat{p} \times \hat{q}. \qquad (18) $$

To perform the rotation of a vector $\mathbf{x}$ into a vector $\mathbf{y}$, that be $\mathbf{y} = \mathbb{R}\mathbf{x}$, a special set of quaternions, called Euler parameters, is employed. For a rotation of angle $\vartheta$ around an axis oriented by the unit vector $\mathbf{n}$, Euler parameters are defined by

$$ \hat{q} = \cos\left(\frac{\vartheta}{2}\right) + n \sin\left(\frac{\vartheta}{2}\right). \qquad (19) $$

The rotation is then applied using the following formula:

$$ \hat{\mathbf{y}} = \hat{\mathbf{q}} \circ \hat{\mathbf{x}} \circ \hat{\mathbf{q}}^*, \qquad (20) $$

which is quadratic with respect to $\hat{\mathbf{q}}$ and is totally equivalent to using a rotation matrix (see Appendix F). Furthermore, if one wants to apply a combination of two successive rotations of respective quaternions $\hat{\mathbf{q}}_1$ and $\hat{\mathbf{q}}_2$, the rotational quaternion $\hat{\mathbf{q}}_2 \circ \hat{\mathbf{q}}_1$, evidencing the multiplicative nature of quaternions, needs to be used in (20), that be:

$$ \hat{\mathbf{y}} = (\hat{\mathbf{q}}_2 \circ \hat{\mathbf{q}}_1) \circ \hat{\mathbf{x}} \circ (\hat{\mathbf{q}}_2 \circ \hat{\mathbf{q}}_1)^*. \qquad (21) $$

Let us notice for practical purposes that the norm of the rotational quaternion defined in equation (19) is 1 (the so-called unit quaternions) obeying

$$ \hat{\mathbf{q}} \circ \hat{\mathbf{q}}^* = \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}} = \hat{\mathbf{1}} = 1 + \mathbf{0}. \qquad (22) $$

Besides, deriving this relation with respect to $s$ leads to the following useful relations

$$ \hat{\mathbf{q}}' \circ \hat{\mathbf{q}}^* = -(\hat{\mathbf{q}} \circ \hat{\mathbf{q}}')^*, \qquad (23) $$

$$ \hat{\mathbf{q}}'^* \circ \hat{\mathbf{q}} = -(\hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}'). $$

In order to manipulate equations in quaternion space, let us also introduce the following helpful relations, stemming from (17), (18) and the commutativity of inner product:

* for arbitrary quaternions $\hat{p}, \hat{q}$ and $\hat{r}$:

$$ \hat{p} \cdot (\hat{q} \circ \hat{r}) = \hat{q} \cdot (\hat{p} \circ \hat{r}^*) = (\hat{p} \circ \hat{r}^*) \cdot \hat{q}, \qquad (24) $$

$$ \hat{p} \cdot (\hat{r} \circ \hat{q}) = \hat{q} \cdot (\hat{r}^* \circ \hat{p}) = (\hat{r}^* \circ \hat{p}) \cdot \hat{q}. $$
<table><caption>**Table 1**
Free-end position and rotation of the cantilever beam under the end-moment $M_L = 3\pi \frac{Eh}{L}$ for several values of n and Ne on the left; deformed shape during loading on the right.</caption><thead><tr><th>n</th><th>N<sub>e</sub></th><th>x<sub>1</sub>(L)</th><th>x<sub>2</sub>(L)</th><th>θ<sub>3</sub>(L)</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>2</td><td>0.03376</td><td>0.29412</td><td>9.65337</td></tr><tr><td>5</td><td>0.00104</td><td>0.21346</td><td>9.43450</td></tr><tr><td>10</td><td>0.00007</td><td>0.21228</td><td>9.42542</td></tr><tr><td>20</td><td><strong>0.00000</strong></td><td><strong>0.21221</strong></td><td>9.42482</td></tr><tr><td rowspan="4">2</td><td>2</td><td>-0.00001</td><td>0.21100</td><td>9.42464</td></tr><tr><td>5</td><td>0.00000</td><td>0.21224</td><td>9.42484</td></tr><tr><td>10</td><td><strong>0.00000</strong></td><td><strong>0.21221</strong></td><td><strong>9.42478</strong></td></tr><tr><td>5</td><td>2</td><td>0.00000</td><td>0.21221</td><td>9.42478</td></tr><tr><td colspan="2">Exact</td><td>0</td><td>0.21221</td><td>9.42478</td></tr></tbody></table>

<figure><img src="image_3.png" /></figure>

• for pure quaternions $\hat{a}$ and $\hat{b}$:

$$ \hat{a} \circ \hat{b} + \hat{b}^* \circ \hat{a} = 2\hat{a} \times \hat{b}. \qquad (25) $$

## 2.5. Rewriting of governing equations

The weak form of equilibrium equations (7) along with (4) and (9) are now expressed in terms of the quaternion parameters. The rewriting operation consists in replacing vectors by their quaternion equivalent $x \to \hat{x}$ and the rotation operation by its expression in the quaternion space $\mathbf{y} = \mathbf{R}x \to \hat{\mathbf{y}} = \hat{\mathbf{q}} \circ \hat{x} \circ \hat{\mathbf{q}}^*$. Most of the relations introduced in the paragraph 2.4 serve this change of rotation parameters.

First, the strain measures (4) are rewritten by combining equations (22), (23) and (25) to get the pure quaternion strains (see Appendix A for the full demonstration)

$$ \hat{\Gamma} = \hat{q}^* \circ \hat{x}'_0 \circ \hat{q} - \hat{\Gamma}_0 \qquad (26a) $$

$$ \hat{\mathbf{K}} = 2\hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}' - \hat{\mathbf{K}}_0, \qquad (26b) $$

with the intrinsic translational strain $\hat{\Gamma}_0 = \hat{u}_1$ and rotational strain $\hat{\mathbf{K}}_0 = 2\hat{\mathbf{q}}'_0 \circ \hat{\mathbf{q}}'_0$ expressed in quaternion space. The virtual strain expressions (9) in the quaternion space are obtained by differentiating equations (26a) and (26b) and rearranging them with (25) (see Appendix B or [12] for the full demonstration) leading to

$$ \delta \hat{\Gamma} = \hat{q}^* \circ \delta \hat{x}'_0 \circ \hat{q} + 2 \hat{q}^* \circ (\hat{x}'_0 \times (\delta \hat{q} \circ \hat{q}^*)) \circ \hat{q}, \qquad (27a) $$

$$ \delta \hat{\mathbf{K}} = 2 \hat{\mathbf{q}}^* \circ (\delta \hat{\mathbf{q}} \circ \hat{\mathbf{q}}^*)' \circ \hat{\mathbf{q}}, \qquad (27b) $$

with the virtual positional (pure) quaternion $\delta \hat{x}_0$ and the virtual rotational quaternion $\delta \hat{q}$. The resultant of the external moment is conjugated to the virtual rotational vector in the virtual work principle (7). To express the resultant of the external moment in the quaternion space, a relationship between the virtual rotational vector and the virtual rotational quaternion must be established. One demonstrates (see Appendix D) that we have

$$ \delta \hat{\vartheta} = 2 \hat{q}^* \circ \delta \hat{q}. \qquad (28) $$

Substituting the quaternion expressions (26)–(28) to their vectorial counterpart in the virtual work principle (7) and using the relations (24) to obtain the resultant of the external moment, the new virtual work principle extended in quaternion space writes

$$ \int_0^L (\hat{\mathbf{N}} \cdot \delta \hat{\Gamma} + \hat{\mathbf{M}} \cdot \delta \hat{\mathbf{K}}) \, ds = \int_0^L (\hat{\mathbf{n}}_e \cdot \delta \hat{\mathbf{x}}_0 + 2(\hat{\mathbf{q}} \circ \hat{\mathbf{m}}_e) \cdot \delta \hat{\mathbf{q}}) \, ds + [\hat{\mathbf{N}}_e \cdot \delta \hat{\mathbf{x}}_0]^L_0 + [2(\hat{\mathbf{q}} \circ \hat{\mathbf{M}}_e) \cdot \delta \hat{\mathbf{q}}]^L_0. \qquad (29) $$

As noticed in paragraph 2.4, the four components of a rotational quaternion $\hat{\mathbf{q}}$ are not independent but constrained by the condition of unity $\hat{\mathbf{q}} \cdot \hat{\mathbf{q}} - 1 = 0$. This constraint is taken into account by the method of Lagrangian multipliers. Following the approach of Zupan et al. [15], the Lagrangian multiplier appended to this constraint and denoted $\mu$ is considered as a variable of the problem. The expression intro-

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Influence of element order on the algorithm accuracy for the example of paragraph 5.1 (number of elements circled on second curve).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Cantilever 45-degree bend subject to a terminal out-of-plane force.</figcaption></figure>

duced in the new virtual work principle (29) hence is the term $\mu(\hat{q} \cdot \hat{q} - 1)$ differentiated and integrated along the length of the beam, that be:

$$ \delta C_{\mu} = \delta \int_{0}^{L} \mu(\hat{q} \cdot \hat{q} - 1) ds = \int_{0}^{L} ((\hat{q} \cdot \hat{q} - 1)\delta\mu + 2\mu\hat{q} \cdot \delta\hat{q}) ds. \quad (30) $$

First, introducing the expressions of $\delta\hat{\Gamma}$ and $\delta\hat{K}$ (27a) and (27b), then rearranging (29) with equations (24)–(25), integrating by parts and applying the fundamental lemma of calculus of variations yield the strong form of equilibrium equations (see Appendix C for the full demonstration):

$$ \left\{ \begin{array}{ll} \text{Vec}((\hat{q} \circ \hat{N} \circ \hat{q}^*)') + n_e & = \mathbf{0}, \\ \text{Vec}((\hat{q} \circ \hat{M} \circ \hat{q}^*)' + \hat{x}'_0 \times (\hat{q} \circ \hat{N} \circ \hat{q}^*)) + \hat{m}_e & = \mathbf{0}, \\ \mu & = \mathbf{0}, \\ \hat{q} \cdot \hat{q} - 1 & = \mathbf{0}, \end{array} \right. \quad \forall s \in [0, L] \quad (31) $$

along with the natural boundary conditions

$$ \left\{ \begin{array}{ll} \text{Vec}(\hat{q} \circ \hat{N} \circ \hat{q}^*) - \mathbf{N}_e = \mathbf{0} & \text{at } s = 0, L, \\ \text{Vec}(\hat{q} \circ \hat{M} \circ \hat{q}^*) - \hat{M}_e = \mathbf{0} & \text{at } s = 0, L, \end{array} \right. \quad (32) $$

or the essential boundary conditions

$$ \left\{ \begin{array}{ll} x_0 = x_d, & \text{at } s = 0, L, \\ \text{Vec}(\hat{q}) = \text{Vec}(\hat{q}_d) & \text{at } s = 0, L, \end{array} \right. \quad (33) $$

where $\hat{q}_d$ is the prescribed rotation under quaternion form (19). It ensues from (31) that $\mu(s) = 0$ is solution of the continuous problem while the remaining equations are strictly equivalent to the original form (11). Introducing the Lagrangian multiplier in the virtual work principle thus do not alter the equilibrium equations. Besides, it should be noticed that only the vector part of the quaternion is prescribed in the boundary conditions (33). The scalar component is indeed automatically known from the fulfillment of the unity constraint (31d). Another equivalent possibility is to replace (33b) by $\hat{q} = \hat{q}_d$ with $|\hat{q}_d| = 1$ and $\mu = 0$ at $s = 0, L$.

At this point, the continuous form of equilibrium equations of the problem is established. The equations are expressed in both weak form through equations (29)–(30) and strong form through equations

<table><caption>Table 2
Free-end position of the bent cantilever beam for F = 600.</caption><thead><tr><th>Formulation</th><th>Order</th><th>x<sub>1</sub>(L)</th><th>x<sub>2</sub>(L)</th><th>x<sub>3</sub>(L)</th></tr></thead><tbody><tr><td rowspan="4">Present</td><td>n = 0</td><td>47.3061</td><td>15.7009</td><td>53.3364</td></tr><tr><td>n = 1</td><td>47.1507</td><td>15.6845</td><td>53.4751</td></tr><tr><td>n = 2</td><td>47.1501</td><td>15.6847</td><td>53.4755</td></tr><tr><td>n = 3</td><td>47.1501</td><td>15.6847</td><td>53.4755</td></tr><tr><td>Bathe and Blorchi [41]</td><td>n = 0</td><td>47.2</td><td>15.9</td><td>53.4</td></tr><tr><td>Simo and Vu-Quoc [10]</td><td>n = 0</td><td>47.23</td><td>15.79</td><td>53.37</td></tr><tr><td>Zupan et al. [15]</td><td>n = 15</td><td>47.4159</td><td>15.2861</td><td>53.4725</td></tr></tbody></table>

(31) along with (32), knowing relations (26a)–(27b). A vector equation is made up of 3 scalar equations. As a result, the set of governing equations in strong form (31) is composed of 8 (3 + 3 + 1 + 1) scalar equations. Besides, the unknowns of the problem are the 8 following fields of s: the 4 components of the rotational quaternion $\hat{q} = [q_0(s) \quad q_1(s) \quad q_2(s) \quad q_3(s)]^T$ depicting the rotation of the sections, the 3 components of the current position vector of the cross-sections centroids $x_0(s) = [x_1(s) \quad x_2(s) \quad x_3(s)]^T$ and the Lagrangian multiplier $\mu(s)$. As a result, our problem is a well-posed 8 equations – 8 unknowns problem and can henceforth be solved numerically.

For convenience, let us gather the 8 scalar unknown fields under a single vector field

$$ z(s) = [x_0(s)^T \quad \hat{q}(s)^T \quad \mu(s)]^T. \quad (34) $$

This allows to write the weak form of equilibrium equations (29)–(30) under the compact form

$$ \int_0^L \boldsymbol{\Sigma} \cdot \delta\boldsymbol{\varepsilon} \, ds = \int_0^L (f_e + f_\mu) \cdot \delta z \, ds + [F_e \cdot \delta z]_0^L, \quad (35) $$

with

$$ \begin{align} f_e &= [n_e^T (2\hat{q} \circ \hat{m}_e)^T 0]^T, \quad F_e = [N_e^T (2\hat{q} \circ \hat{M}_e)^T 0]^T, \tag{36} \\ f_\mu &= [0 (2\mu\hat{q})^T \hat{q} \cdot \hat{q} - 1]^T. \end{align} $$

In equation (35), the strain and the stress resultant vector are still equal to $\boldsymbol{\varepsilon} = [\Gamma^T \quad \mathbf{K}^T]^T$ and $\boldsymbol{\Sigma} = [\mathbf{N}^T \quad \mathbf{M}^T]^T$ but from this point forward are expressed in quaternion space, that be $\Gamma = \text{Vec}(\hat{\Gamma})$, $\mathbf{K} = \text{Vec}(\hat{\mathbf{K}})$, $\mathbf{N} = \text{Vec}(\hat{\mathbf{N}}) = \mathbf{C}_N \text{ Vec}(\hat{\Gamma})$ and $\mathbf{M} = \text{Vec}(\hat{\mathbf{M}}) = \mathbf{C}_M \text{ Vec}(\hat{\mathbf{K}})$, with $\hat{\Gamma}$ and $\hat{\mathbf{K}}$ expressed as in (26a) and (26b).

Equations (26a) and (26b) show clearly that $\hat{\Gamma}, \hat{\mathbf{K}}$ then $\hat{\mathbf{N}}, \hat{\mathbf{M}}$ through (8) and as a consequence $\boldsymbol{\varepsilon}, \boldsymbol{\Sigma}$ are polynomial functions of $z$ of degree 3. Similarly, equations (27a) and (27b) expose that $\delta\boldsymbol{\varepsilon}$ is a polynomial function of $z$ of degree 4. This in addition to the quadratic form of expressions (36) evidences that the governing equation (35) is a polynomial equation of $z$ of degree 7. Thus, using quaternion parameters instead of the rotational vector has allowed to transform the non-polynomial system of 6 equations/unknowns (11) (it contains trigonometric functions due to the expression of the rotation and the tangent operators (6) and (10)) in a polynomial system of 8 equations – 8 unknowns. How our method takes advantage of this specificity is shown in the following parts.

### 3. Finite element method

#### 3.1. Discrete form of equations

The finite element method is now applied to the weak form of governing equations (35) to obtain the discrete equilibrium equations. The beam is cut in $N_e$ elements, which for sake of simplicity are supposed of equal length $L_e = L/N_e$. Each element is composed of $n + 2$ equally-spaced nodes: $n$ internal nodes and one node at each extremity of
the element. In each element $e$, the unknown field $z^{(e)}$ value at node $i \in \{1, \dots, n+2\}$ is denoted $z^i$.

The virtual work principle (35) is discretized using a standard Galerkin approach in which the trial functions $x_0, \hat{q}$ and $\mu$ and the test functions $\delta x_0, \delta \hat{q}$ and $\delta \mu$ are interpolated similarly. This choice is not obvious when dealing with rotational quantities and for this reason it is discussed briefly in section 3.4. The same interpolation order is used for all variables, since numerical studies have shown that a different order for each variable may lead to numerical divergence [26]. The interpolation then writes

$$z^{(e)}(s) = \sum_{i=1}^{n+2} N_i(s) z^i. \qquad (37)$$

We use isoparametric elements so that the interpolation functions $N_i$ in equation (37) are the same as the shape functions and are chosen as standard Lagrange-type polynomials. Equation (37) may be written under the compact matrix form

$$z^{(e)}(s) = P(s)u^{(e)}, \qquad (38)$$

where, if $I_8$ is the $8 \times 8$ identity matrix,

$$P(s) = [N_1(s)I_8 \quad N_2(s)I_8 \quad \dots \quad N_{n+2}(s)I_8] \qquad (39)$$

is an interpolation matrix of size $8 \times 8(n+2)$ and $u^{(e)}$ contains the unknowns nodal values of element $e$:

$$u^{(e)} = [(z^1)^T \quad (z^2)^T \quad \dots \quad (z^{n+2})^T]^T. \qquad (40)$$

As the elementary fields all depend on $s$ and as all the equations are written inside the element, except when needed for a better comprehension, the notations $e$ and $s$ are dropped in the following. As said earlier, the test functions are interpolated similarly to the trial functions so that we also have $\delta z^{(e)}(s) = P(s)\delta u^{(e)}$. However, it is convenient for the writing to introduce a 2nd interpolation matrix $Q$ of size $12 \times 8(n+2)$ for the test functions such that

$$[(\delta x_0')^T \quad (\delta \hat{q}')^T \quad (\delta \hat{q})^T \quad \delta \mu]^T = Q \delta u^{(e)}, \qquad (41)$$

with

$$Q = \begin{bmatrix} N'_1 I_3 & 0 & 0 & \dots & N'_{(n+2)} I_3 & 0 & 0 \\ 0 & N'_1 I_4 & 0 & \dots & 0 & N'_{(n+2)} I_4 & 0 \\ 0 & N_1 I_4 & 0 & \dots & 0 & N_{(n+2)} I_4 & 0 \\ 0 & 0 & N_1 & \dots & 0 & 0 & N_{(n+2)} \end{bmatrix}. \qquad (42)$$

We can henceforth write under discrete form the continuous equation (35). To get the discrete version of the left-hand side, we search for an expression of the virtual strains $\delta \mathcal{E}$ as a function of $\delta u^{(e)}$. From the

expression of the virtual strains (27a) and (27b) and from (41), one shows (see Appendix E) that

$$\begin{bmatrix} N \\ M \end{bmatrix} \cdot \delta \mathcal{E} = Q^T g(N, M) \cdot \delta u^{(e)}, \qquad (43)$$

with the vector of 12 components

$$g(N, M) = \begin{bmatrix} \mathrm{Vec} \left[ \hat{q}^* \circ \hat{N} \circ \hat{q} \right] \\ 2\hat{q} \circ \hat{M} \\ 2\hat{q} \circ \left( \hat{N} \times (\hat{\Gamma} + \hat{\Gamma}_0) + \hat{M} \times (\hat{K} + \hat{K}_0) \right) + 2\hat{q}' \circ \hat{M} \end{bmatrix}. \qquad (44)$$

$Q^T g(\cdot, \cdot)$ may be interpreted as an elementary discrete gradient of the internal work, since it links $\Sigma \cdot \delta \mathcal{E}$ to the virtual nodal values.

The discrete version of the right-hand side of equation (35) in a given element stems directly from equations (38) and (41) and the virtual work principle may now be written in each element in a discrete manner. The elementary virtual work of internal forces, the elementary virtual constraint and the elementary virtual work of external forces respectively write

$$\begin{align} \delta W_i^{(e)} &= - \int_0^{L_e} \Sigma \cdot \delta \mathcal{E} \ ds = \delta u^{(e)} \cdot \left( - \int_0^{L_e} Q^T g(N, M) ds \right), \nonumber \\ \delta C_\mu^{(e)} &= - \int_0^{L_e} \mathbf{f}_\mu \cdot \delta \mathbf{z} \ ds = \delta u^{(e)} \cdot \left( - \int_0^{L_e} Q^T c_\mu \ ds \right), \tag{45} \\ \delta W_e^{(e)} &= \int_0^{L_e} \mathbf{f}_e \cdot \delta \mathbf{z} \ ds = \delta u^{(e)} \cdot \int_0^{L_e} \mathbf{P}^T \mathbf{f}_e \ ds. \nonumber \end{align}$$

In (45), we have introduced the vector $c_\mu = [0 \ 0 \ 2\mu \hat{q}^T \hat{q} \ \hat{q} - 1]^T$. Besides, for the special cases of end elements, the loads on end-sections must be added to external work. They write respectively $-F_e(0) \cdot \delta z(0) = \delta u^{(1)} \cdot (-P^T(0)F_e(0))$ and $F_e(L) \cdot \delta z(L) = \delta u^{(N_e)} \cdot (P^T(L_e)F_e(L))$ for the first and the last element. Summing these virtual quantities, the elementary virtual work principle takes the form

$$\delta W_i^{(e)} + \delta W_e^{(e)} + \delta C_\mu^{(e)} = \delta u^{(e)} \cdot (\mathbf{F}_e^{(e)} - \mathbf{F}_i^{(e)}), \qquad (46)$$

where $\mathbf{F}_i^{(e)}$ and $\mathbf{F}_e^{(e)}$ are respectively the elementary internal force vector (in which is included the unity constraint) and the elementary external force vector. They write

$$\mathbf{F}_i^{(e)} = \int_0^{L_e} Q^T \mathbf{f}_i^{(e)} ds, \quad \mathbf{F}_e^{(e)} = \int_0^{L_e} \mathbf{P}^T \mathbf{f}_e ds. \qquad (47)$$

As in (47), the vector $\mathbf{F}_e^{(e)}$ requires a special treatment at beam ends: the end forces $-P^T(0)F_e(0)$ and $P^T(L_e)F_e(L)$ must be included in the first and

<figure><img src="image_6.png" /><figcaption>Fig. 4. Bifurcation diagrams for the deep circular arch. Displacement of the force application point in $u_2$ direction (left) and $u_3$ direction (right) vs. the force applied. The 4 branches are numbered and their corresponding deformed shapes are found on Fig. 5.</figcaption></figure>
the last element respectively. In this equation, the vector $\mathbf{f}_i^{(e)}$ is equal to

$$ \mathbf{f}_i^{(e)} = \begin{bmatrix} \mathrm{Vec}[\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*] \\ 2\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \\ 2\hat{\mathbf{q}} \circ (\hat{\mathbf{N}} \times (\hat{\mathbf{\Gamma}} + \hat{\mathbf{\Gamma}}_0) + \hat{\mathbf{M}} \times (\hat{\mathbf{K}} + \hat{\mathbf{K}}_0)) + 2\hat{\mathbf{q}}' \circ \hat{\mathbf{M}} + 2\mu \hat{\mathbf{q}} \\ \hat{\mathbf{q}} \cdot \hat{\mathbf{q}} - 1 \end{bmatrix}. \quad (48) $$

Few comments need to be made on this expression because it is at the heart of the computer code, since the internal force vector is directly linked to $\mathbf{f}_i^{(e)}$. First, let us notice that it does not make use of the rotation matrix $\mathbf{R}$ which is replaced by computationally more efficient quaternion products (see for instance [24]). This product is thus directly implemented in the code. Afterward, recalling that $\mathbf{N} = \mathbf{C}_N\mathbf{\Gamma}, \mathbf{M} = \mathbf{C}_M\mathbf{K}$ and that the strains expressions are given by (26a) and (26b), one evidences that this expression is still polynomial of degree 7 but this time function of elementary fields. Finally, as quaternion expressions have 4 components and vectorial expression have 3, one observes that $\mathbf{f}_i^{(e)}$ is a vector of 12 (3 + 4 + 4 + 1) components. After being multiplied by $\mathbf{Q}^T$, it gives the vector $\mathbf{F}_i^{(e)}$ whose size is equal to the number of nodal values of the element, namely $8(n + 2)$.

In equations (47), the integrals over the element are evaluated numerically by a classical Gauss-Legendre quadrature of order $N_g$. Integration is performed in the parent element characterized by the coordinate $\xi$ varying between -1 and 1. The Jacobian of the transformation from physical to parent space is the same for all elements and is equal to $J = ds/d\xi = L_e/2$. One of the major problems of Timoshenko-like elements is the presence of shear locking [27]. In order to avoid this numerical deficiency, reduced integration is used. It consists in choosing a quadrature order lower than the one required for exact numerical integration. In our case, for an element with $n + 2$ nodes, reduced integration is employed by using $N_g = n + 1$ points in the quadrature formula [28].

At this point, all the necessary information has been given to compute elementary quantities. The last remaining step of the finite element method is to make an assembly operation on the elements to get the global discrete problem equations. For that purpose, let us define the global unknowns vector $\mathbf{U}$ containing all the degrees of freedom of the problem:

$$ \mathbf{U} = [(\mathbf{x}_0^1)^T \quad (\hat{\mathbf{q}}^1)^T \quad \mu^1 \quad \dots \quad (\mathbf{x}_0^{N_n})^T \quad (\hat{\mathbf{q}}^{N_n})^T \quad \mu^{N_n}]^T, \quad (49) $$

where $N_n$ is the total number of nodes. Assembling expression (44) for all elements and stating that the virtual work principle holds for any function $\delta\mathbf{U}$, the global discrete problem simply writes

$$ \mathbf{R}(\mathbf{U}) = \mathbf{F}_i(\mathbf{U}) - \mathbf{F}_e(\mathbf{U}) = \mathbf{0}, \quad (50) $$

with the global internal force vector $\mathbf{F}_i$ and the global external force vector $\mathbf{F}_e$.

With eight degrees of freedom per node and $N_n = (n + 1)N_e + 1$ nodes, the finite element method has transformed the 8 continuous equilibrium equations (31) in a set of $N_{eq} = 8(n + 1)N_e + 8$ equations forming the nonlinear global problem to solve.

<table><caption>**Table 3**

Critical forces obtained numerically with our code and the code of several authors for in-plane buckling of the deep circular arch.</caption><thead><tr><th></th><th>N<sub>e</sub></th><th>10</th><th>12</th><th>20</th><th>40</th></tr></thead><tbody><tr><td>Present</td><td>n = 1</td><td>905.3</td><td>901.3</td><td>897.8</td><td>897.3</td></tr><tr><td>Zupan and Saje [34]</td><td>n = 1</td><td></td><td>897.3</td><td></td><td></td></tr><tr><td>Cesarek et al. [42]</td><td>n = 1</td><td></td><td></td><td>906.57</td><td>899.69</td></tr><tr><td>Ibrahimbegovic et al. [9]</td><td>n = 1</td><td></td><td></td><td>897.3</td><td></td></tr><tr><td>Simo and Vu-Quoc [10]</td><td>n = 0</td><td></td><td></td><td></td><td>905.28</td></tr></tbody></table>

## 3.2. User control of a quasi-static problem

In the process of controlling equation (50) with a parameter $\lambda$, we can make the distinction between 2 main types of control: load control (force moment) and displacement/rotation control. In the former case, the control parameter is inserted into the known external force vector $\mathbf{F}_e$ such that the problem (50) becomes

$$ \mathbf{R}(\mathbf{U}, \lambda) = \mathbf{F}_i(\mathbf{U}) - \mathbf{F}_e(\mathbf{U}, \lambda) = \mathbf{0}. \quad (51) $$

In the latter case, $\lambda$ is included in the expression of the prescribed displacement/rotation degrees of freedom of vector $\mathbf{U}$ (a weight might be associated to $\lambda$ when there are several prescribed displacements/rotations). As a degree of freedom is removed of the system for each imposed displacement/rotation, the corresponding equations are also removed of the system: the system would otherwise be overdetermined. These equations are not necessary to get the equilibrium configuration but are not hollow though since they give the reaction force at locked degrees of freedom. Mathematically, we can express control of the j-th degree of freedom with a prescribed displacement $u_d$ through

$$ \begin{gathered} \mathbf{R}(\mathbf{U}^*, \lambda) = \mathbf{F}_i^*(\mathbf{U}^*) - \mathbf{F}_e^*(\mathbf{U}^*) = \mathbf{0}, \\ \mathbf{U}^* = [u_1 \quad u_2 \quad \dots \quad u_{j-1} \quad u_d(\lambda) \quad u_{j+1} \quad \dots \quad u_{N_n}]^T, \end{gathered} \quad (52) $$

where $\mathbf{F}_i^*$ and $\mathbf{F}_e^*$ are the force vectors from which the j-th component has been removed.

Finally, let us comment the implementation of boundary conditions (clamped or pinned ends) in the presence of quaternions. Boundary conditions are applied classically in a Boolean manner: the equations corresponding to locked degrees of freedom are removed from the system. However, because the 4 components of a quaternion are constrained, the rotational degrees of freedom need to be dealt with care. In particular, it stems from the expression of essential boundary conditions (33) that only the vector part of the quaternion needs to be prescribed: the fourth component is then determined by the unity constraint. It follows that to block all the rotations at one node, the three equations corresponding to $q_1, q_2$ and $q_3$ at this node are to be removed. To partially block the rotations (two or less rotation directions), the equations corresponding to these quaternion components are classically removed. For instance, to block the rotations around $e_2$ and $e_3$ at one node (and supposing for sake of simplicity that the axes of the corresponding section are aligned with global axes), it suffices to remove the equations corresponding to $q_2$ and $q_3$ at this node.

## 3.3. Assembly of beams/rigid joints

With the purpose of simulating cables, it is necessary to be able to assemble several beams. Indeed, in automotive industry, cables often have a harness-type geometry. There exist several great references which explain how to model the different types of joints for flexible bodies, for instance [22] or [29]. However, it is less often expressed in terms of quaternions. We thus explain here how we implemented a rigid joint in our simulation as a basis for modeling other types of joints.

Let us assume that we want to bind the extremities of two beams that we will denote with roman numbers I and II. Each beam, discretized with the finite element method, has one node at its extremity. To create a rigid joint between the ends of the two beams, one must constrain the position of the 2 nodes to remain the same and the rotation between the cross-sections attached at each node to not vary along the deformation. If we denote respectively $\mathbf{z}^I = [(\mathbf{x}_0^I)^T \quad (\hat{\mathbf{q}}^I)^T \quad \mu^I]^T$ and $\mathbf{z}^II = [(\mathbf{x}_0^II)^T \quad (\hat{\mathbf{q}}^II)^T \quad \mu^II]^T$ the degrees of freedom of nodes I and II in the current state, the constraint writes for the position:

$$ \mathbf{x}_0^I - \mathbf{x}_0^II = \mathbf{0}. \quad (53) $$

The rotation constraint is a bit more tricky to set up. Let us introduce the quaternions at each node in the reference configuration $\hat{\mathbf{q}}_0^I$ and
$\hat{q}_0^{II}$. The variation of the rotation between the reference and the actual configuration for the nodes I and II may be represented by the two quaternions $\Delta\hat{q}^I$ and $\Delta\hat{q}^{II}$ which, following the rule of combination of rotations (21), then meet

$$ \hat{q}^I = \Delta\hat{q}^I \circ \hat{q}_0^I, \quad \hat{q}^{II} = \Delta\hat{q}^{II} \circ \hat{q}_0^{II}. \qquad (54) $$

As the variation of rotation should be the same for the bound nodes, that being $\Delta\hat{q}^I = \Delta\hat{q}^{II}$, and using relations (54) and (22) the constraint finally writes as a function of the variables of the problem

$$ \hat{q}^I \circ \hat{q}_0^{I*} - \hat{q}^{II} \circ \hat{q}_0^{II*} = \hat{0}. \qquad (55) $$

As recalled in the previous section, only three components of the quaternions need to be prescribed. Consequently, only three constraints have to be added to the equations and we choose here to use only the three components corresponding to the vector part of (55). The Lagrangian method is chosen to insert these constraints in the equations. Introducing the vector of constraints $C_{jt}$ for the joint which is then equal to

$$ C_{jt} = \left[ \mathrm{Vec} \left( \hat{q}^I \circ \hat{q}_0^{I*} - \hat{q}^{II} \circ \hat{q}_0^{II*} \right) \right], \qquad (56) $$

<figure><img src="image_7.png" /><figcaption>(1) In-plane solution.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(3) Second out-of-plane solution.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(2) First out-of-plane solution.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(4) Third out-of-plane solution.</figcaption></figure>

and the associated Lagrangian multiplier $\mu_{jt}$, the constraint is included in the model by adding the quantity $\mu_{jt} \cdot \delta C_{jt} + C_{jt} \cdot \delta \mu_{jt}$ in the discrete form of the virtual work principle (29).

### 3.4. Comments on trial functions interpolation

As the virtual work principle (29) holds for any test functions $\delta x_0$, $\delta q$ and $\delta \mu$, their interpolation can be carried out without precautions using standard schemes. On the contrary, more careful attention needs to be paid to the interpolation of trial functions [23]. The position $x_0$ and the Lagrangian multiplier $\mu$ both belong to a linear vector space, so the classical interpolation (37) is adequate. However, it is not the case for the rotations $\hat{q}$ which belong to a non-Euclidean space and thus are non-additive quantities.

For this reason, many articles focus on the subject [30–32,13]. Among the existing interpolation, the most promising is the Spherical Linear Interpolation (SLERP), widely used in computer graphics [31]. However, this interpolation is thus far limited to linear elements ($n = 0$) and the attempt made by Ref. [13] to extend it to superior orders has unfortunately been inconclusive. Recently, the work of [33] on this matter looks interesting but it is not fully developed yet. Other non-additive methods have been presented such as interpolating the incre-

<caption>Fig. 5. Deep circular arch example. Deformed shape at several solution points for each branch. Corresponding branch number in Fig. 4 in parenthesis.</caption>
<figure><img src="image_11.png" /><figcaption>Fig. 6. Hinged right-angle frame under conservative vertical load: bifurcation diagram and deformed shapes for several intensities of the force. On the bifurcation diagram, the displacement of the load application point is displayed.</figcaption></figure>

mental rotation or using the curvature as primary variable [34,35] but are not adapted to the quadratic formalism of the ANM.

Departing from these more elaborate interpolations, our use of a classical high order interpolation has been motivated by several arguments. The first reason is that this interpolation does not spoil the polynomial shape of the system (50) and allows to use the ANM. Then, the conclusions of the work of Romero [31] show that no matter what interpolation is chosen, the inaccuracies are always cured by increasing the order of interpolation and the mesh refinement. Thus, only the computational cost is affected by a choice of interpolation less adapted to the rotation space. Finally, the classical additive interpolation benefits from an easy implementation while the interpolation devised on the above quoted article may reveal very tricky.

## 4. Asymptotic numerical method

To the best knowledge of the authors, up to now, solving the nonlinear systems (51) or (52) has always been carried out using classical predictor-corrector methods. In this part, the Asymptotic Numerical Method (ANM) introduced by Cochelin [19] and Damil and Potier-Ferry

[18] is presented. It is a high order perturbation technique that allows to compute a large part of a nonlinear branch with only one matrix inversion. The method principle is, starting from a solution point, to seek the branch of solution as an asymptotic expansion of a path parameter. This leads to a semi-analytical solution, whose range of validity is automatically calculated. Applied step-by-step, it allows to describe the full branch of solution. Application of the method to the herein described mechanical problem is presented in the following.

### 4.1. *Quadratic recast*

As a first step, the nonlinear system of equilibrium equations (50) is recast into a quadratic form. The ANM framework allows a systematic solving of this particular form of problems thus providing a powerful solver method for a large class of physical problems [19]. This recasting operation is set up by introducing additional variables, referred as the *slave variables* $V$, to the initial *primary variables* contained in the vector $U$. This operation is facilitated by the use of quaternions. As explained earlier, given the expressions of strains (26a) and (26b), the linear material law (8) and the expression of $f_i^{(e)}$ (48), the system (50) is

<figure><img src="image_12.png" /><figcaption>Fig. 7. Hinged right-angle frame under following vertical load: bifurcation diagram and deformed shapes for several intensities of the force. On the bifurcation diagram, the displacement of the load application point is displayed.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 8. Equilibrium curve of the clamped-clamped beam for in-plane and out-of-plane (3D) buckling. Displacement of the center of the beam (s = L/2) in u₂ direction on the left and maximum displacement in u₃ direction on the right.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

If $R_Q^{m+1}$ designates the coefficient of the power series of order $m + 1$, using the approximation $||R_Q(W(a))|| \approx a^{m+1} ||R_Q^{m+1}||$ (see Ref. [36] for the entire proof), the maximum step size is defined by

$$a_{\max} = \left( \frac{\epsilon}{||R_Q^{m+1}||} \right)^{\frac{1}{m+1}}. \qquad (71)$$

A section of the branch of solution is thus defined as $W(a)$, $a \in [0, a_{\max}]$. The next section is calculated by using $W(a_{\max})$ as the new starting point. The full diagram can thereby be calculated iteratively.

### 4.3. Further comments

* The ANM presents the following advantages:

- The branch of solution is known analytically, section by section.

- The tangent matrix $K_t$ needs to be calculated only for order $p = 1$ and can then be reused for superior orders. It considerably outperforms the classical predictor-corrector methods which require an actualization of the Jacobian matrix at each iteration.

- The calculation of an analytical expression of the tangent matrix $K_t$ is not required, since the quadratic form of the problem gives a direct access to its expression through the definition (67).

- The range of validity is calculated automatically and without tuning parameter. Unlike predictor-corrector method, both parameters $\epsilon$ and $m$ do not alter the convergence but solely the accuracy. From empirical results, setting $\epsilon = 1 \times 10^{-7}$ and $15 \le m \le 20$ is found optimal to get the best trade-off between accuracy and speed and do not need any further tuning.

- The ANM benefits from the same robustness as the arc-length method even for sharp turns, due to the use of a pseudo arc-length measure as path parameter. Moreover, for such smooth nonlinearities as in our equilibrium equations, a correction step is not necessary: in practice, it is rare to see a residue greater than $10\epsilon$ [19]. If ever necessary, for instance for non-polynomial nonlinearities, a Newton-Raphson type correction step can easily be coupled to the ANM [37].

* The maximum step size defined in (71) in practice hardly differs from the radius of convergence $a_{rc}$.

* One of the main advantages of the ANM is its ability to follow a branch even when a bifurcation occurs. This remarkable property allows to describe a full bifurcation diagram. To switch from a branch to another, a special procedure is however required. One possible method is to perturb the initial problem with a constant vector of small norm $cP$ such that the problem (63) becomes

$$R_p(W) = C + L(W) + Q(W, W) + cP = 0. \qquad (72)$$

$R_p(W)$ is the new perturbed problem, $P$ a normalized vector of constant random numbers and $c$ is the intensity of the perturbation, assumed small.

### 4.4. The special case of imposed rotations by quaternions in the ANM

When using quaternions, imposing a rotation at one node proves to be a delicate task. Because of the redundancy of quaternions ($\hat{q}$ and $-\hat{q}$ represent the same rotation), applying a rotation in a classical way, as presented in paragraph 3.2, does not permit to impose rotations of angle greater than $\pi$. It is thus necessary to come back to the rotational quaternion definition (19) to control the rotation. One way to perform a rotation of angle $\lambda$ around a unit axis $n$ at node $i$ is under the form of a Lagrangian constraint:

$$C_{\text{rot}} = \hat{q}^i - \begin{bmatrix} \cos(\lambda/2) \\ \sin(\lambda/2) \, n \end{bmatrix} = 0. \qquad (73)$$

However, an additional difficulty arises from using the ANM. This constraint is not in quadratic form and does not fall into the quadratic framework of the ANM. Fortunately, there exists a procedure [38,21] to overcome this difficulty and to quadratically recast equation (73). Let us introduce the two new variables

$$\begin{aligned} s &= \sin(\lambda/2), \\ c &= \cos(\lambda/2). \end{aligned} \qquad (74)$$

added to the global vector of unknowns $W$. Differentiating these expressions leads to

$$\begin{aligned} ds &= \frac{c}{2} d\lambda, \\ dc &= -\frac{s}{2} d\lambda, \end{aligned} \qquad (75)$$

which are quadratic equations. Let us now define the linear and bilinear vector-valued operators $L_h$ and $Q_h$ such that (71) writes

$$L_h(dW) = Q_h(W, dW), \qquad (76)$$

with $L_h(dW) = [ds \, dc]^T$ and $Q_h(W, dW) = [\frac{c}{2}d\lambda - \frac{s}{2}d\lambda]^T$. Karkar [38] noticed that if the vector $W$ is written as the power series expansion (64), differentiating it reads
<figure><img src="image_15.png" /><figcaption>Fig. 9. Deformed shape of the clamped-clamped beam for an out-of-plane buckling.</figcaption></figure>

$$d\mathbf{W} = (\mathbf{W}_1 + 2a\mathbf{W}_2 + \dots + ma^{m-1}\mathbf{W}_m) da. \quad (77)$$

Inserting the power series expansions (64) and (77) in (76) and equating the terms of each order to 0 in a similar manner to paragraph 4.2, leads to the m linear equations

$$\mathcal{J}_h \cdot \mathbf{W}_1 = \mathcal{L}_h(\mathbf{W}_1) - Q_h(\mathbf{W}_0, \mathbf{W}_1) = \mathbf{0} \quad \text{(order 0)},$$

$$\mathcal{J}_h \cdot \mathbf{W}_{p+1} - \sum_{i=1}^{p} \frac{p+1-i}{p+1} Q_h(\mathbf{W}_i, \mathbf{W}_{p+1-i}) = \mathbf{0} \quad \text{(order } 1<p \le m-1), \quad (78)$$

with the tangent operator $\mathcal{J}_h = \mathcal{L}_h(\bullet) - Q_h(\mathbf{W}_0, \bullet)$. By adjoining the system of equations (78) to the main system of equations (66) and solving it, the values of s and c are known on the branch of solution. The prescribed rotations are then finally applied by adding the constraint (73) in the virtual work principle under its virtual form

$$\mu_{\text{rot}} \cdot \delta \mathcal{C}_{\text{rot}} + \mathcal{C}_{\text{rot}} \cdot \delta \mu_{\text{rot}} = \mu_{\text{rot}} \cdot \delta \hat{\mathbf{q}}^i + \left( \hat{\mathbf{q}}^i - \begin{bmatrix} c \\ sn \end{bmatrix} \right) \cdot \delta \mu_{\text{rot}}, \quad (79)$$

with the vector $\mu_{\text{rot}}$ containing the 4 Lagrangian multipliers associated to the constraints.
<figure><img src="image_16.png" /><figcaption>Fig. 10. Initial configuration of the experiment of Lazarus et al. [21]. At the top, pictures of the bench from Ref. [21]. Below, results of our simulation for the 2 initialization steps of the curved rod: (1) unrolling (2) positioning of the right end.</figcaption></figure>

## 5. Validation

The method presented above has been implemented in ManLab [39], an open-source Matlab package. In this section, we show the results of our code on several reference numerical experiments for which an analytical solution is known or which have been tested by other authors. These examples serve to demonstrate the validity and the accuracy of our model as well as its abilities. In these tests, the centerline is chosen such that it passes through the centroid of all the sections (no additional out-of-diagonal terms in equation (8)). For the comparison, both element order n and element number N_e can be tuned.

In figures of deformed shapes, the color map represents the distribution of Von Mises stresses on the surface of the beam.

### 5.1. Cantilever beam under end-moment

We consider an initially straight beam, clamped at s = 0 and sub-

<figure><img src="image_17.png" /><figcaption>Fig. 11. Deformed shape of the initially straight rod for the Lazarus et al. experiment. Corresponding point on the equilibrium curve are found on Fig. 13.</figcaption></figure>

jected to an end-moment $M_L = 3\pi EI_3/L$ around $e_3$ at s = L. This example, studied in numerous papers [40,10], is known analytically to give an arc of a circle. In particular, the analytical position and rotation of the right-end section

$$ 
\begin{aligned}
x_1(L) &= \frac{EI_3}{M_L} \sin \left( \frac{M_L L}{EI_3} \right) = 0, & \theta_3(L) &= \frac{M_L L}{EI_3} = 3\pi, \\
x_2(L) &= \frac{EI_3}{M_L} \left( 1 - \cos \left( \frac{M_L L}{EI_3} \right) \right) = \frac{2L}{3\pi}.
\end{aligned} 
\quad (80) $$

are compared to our numerical solution. The geometrical and material properties chosen for the simulation are: length L = 1 m, circular cross-section of radius R = 4 × 10⁻² m, Young's modulus E = 105 × 10⁹ Pa and shear modulus G = 40.3 × 10⁹ Pa. The numerical results for various element orders n and number of elements N_e are gathered in Table 1.

This example demonstrates a good agreement of our method with the analytical solution since it converges toward the exact solution when the number of elements is increased. Additionally, it shows that increasing element order allows to converge more rapidly. These results are confirmed by Fig. 2 which displays the absolute error on vertical position $x_2(L)$ with respect to both the number of elements and the computation time for n = 0 to n = 3. One also sees on this figure that, for our computer code and on this example, it is worthier to increase the element order than the number of elements to get a small error. This remark is not a general conclusion for our method and should be regarded with precautions since the computation time result is both code dependent and case dependent. However, it is a good indication on our code behavior.

### 5.2. Cantilever 45° bend

This example first considered in Ref. [41] has become a classical numerical experiment since it tests all together bending, shear, extensional and torsional deformations. An initially bent cantilever beam of radius R = 100 m and of angle φ = π/4 rad (thus of length L = φR), contained in the $(u_1, u_2)$ plane as pictured Fig. 3, is subjected to an end force F = 600 N normal to this plane. The cross-section is a unit square while the Young's and shear moduli are respectively E = 10⁷ Pa and

<figure><img src="image_18.png" /><figcaption>Fig. 12. Deformed shape of the initially curved rod for the Lazarus et al. [21] experiment. Corresponding point on the equilibrium curve are found on Fig. 13.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 13. Equilibrium curve of straight (left) and curved (right) Lazarus et al. [21] rods.</figcaption></figure>

$G = E/2$. The mesh is constituted of 8 elements.

The experiment shows a good accordance with the numerical results of other authors, see Table 2. In addition, our algorithm proves to have a good convergence rate for increasing order since beyond $n = 2$ there is no more change on the 4th decimal.

### 5.3. Deep circular arch

In this well-known example, we consider an initially curved beam whose centerline forms an arc of a circle of angle $\phi = 215^\circ$ and of radius $R = 100$ m. The arch right-end is clamped while the left-end is hinged leading to a non-symmetrical structure. The remaining properties are $EA = GA_2 = GA_3 = 10^6$ N and $GJ = EI_2 = EI_3 = 10^8$ Nm$^2$. A vertical load $F_2$ is applied at the middle of the beam.

The in-plane buckling and post-buckling of the beam are then studied. This experiment allows to test the ability of our code to detect accurately a buckling load. The results are compared to several results of the literature [9,34,10,42] and to the reference solution of DaDeppo and Schmidt [43], $F_2^1 = 897$ N accurate up to 3 significant digits. The results shown on Table 3 totally concurs with other authors numerical findings. One notices that the numerical solution converges faster toward the analytical solution for Zupan and Saje [34] and Ibrahimovic [9] though. This difference is mainly due to the different type of interpolation chosen (see paragraph 3.4).

Besides, as reported by Cesarek et al. [42], this example has been traditionally studied only as an in-plane problem. Nevertheless, when considered in 3D, this example is much more complex and the arch experiments out-of-plane buckling. Consequently, it appears as a perfect test for the branch-following feature of the ANM. As illustrated on Fig. 4, the ANM allows to describe the full diagram. To do so the perturbation coefficient $c$ of equation (72) is set to $10^{-6}$. This perturbation is applied just before the bifurcation and then removed when the equilibrium point is “far enough” from the unstable branch. The corresponding deformed shapes are illustrated on Fig. 5.

On Fig. 4, the branch (1) is the in-plane branch of solution. Then, one can notice three “pairs” of pitchfork-type bifurcated branches. Each pair contains two symmetric branches, a $u_3$ positive and a $u_3$ negative, due to the symmetry of the problem. The branch (2) is very similar in shape to the one found in Ref. [42] but the critical load differs since we find a value of $F_2^3 = 78$ N. In the absence of other comparisons, we will consider the exact value remains outstanding. The two other bifurcation points (branches (3) and (4)), to our best knowledge never mentioned in any other papers, are found at $F_2^2 = 643$ N before the in-plane turning point and at $F_2^4 = 790$ N right after this point.

### 5.4. In-plane hinged right-angle frame

In this example, for which a numerical and an analytical solution

can be found in Simo and Vu-Quoc [10] and Argyris and Symeonidis [44], the buckling and post-buckling behavior of a frame hinged at both ends is under study. The frame is made of two perpendicular legs linked by a rigid joint as presented in paragraph 3.3. Both legs are identical of length $L = 120$ m. Their cross-section area is $A = 2$ m$^2$, their quadratic moment $I = 6$ m$^4$, their Young's modulus $E = 7.2 \cdot 10^6$ Pa and their Poisson's ratio $\nu = 0.3$. The loading consists in a vertical force applied at one fifth of the length of the horizontal leg. The example is looked at for both a conservative and a following load (see Figs. 6 and 7).

When running the example for a mesh of 10 quadratic elements per leg, a critical load $F_{crit}$ of 18552 N is obtained in the conservative case which matches with 3 significant digits the load of 18532 N obtained in Ref. [10] and the load of 18550 N obtained in Ref. [44]. For the following load case, we get a critical load of 35334 N against 35447 N in Ref. [10] and 36000 N in Ref. [44]. As our results concur with the results of other authors, it shows the validity of our method for a complex structure having joints. Once again, let us notice that the ANM follows the equilibrium effortlessly, even for the nonconservative load and despite the presence of a joint.

Besides, although it is not displayed in this paper for brevity, this example exhibits out-of-plane equilibrium branches which could be followed easily with the ANM.

### 5.5. Out-of-plane buckling of a clamped-clamped beam

The classical example illustrated on Fig. 9, is a buckling problem for

<figure><img src="image_20.png" /><figcaption>Fig. 14. Flexible cable designed shape in a CAD environment.</figcaption></figure>
which the beam is clamped at both ends and a force F applied in the axis of the beam on the right end. The beam measures L = 1 m, has a circular cross-section of radius R = 0.02 m and material properties E = 105 GPa and ν = 0.3. As the out-of-plane buckling is under study, a very small displacement ε = 10⁻⁶ in **u**₃ direction is applied before loading. In addition, for practical purposes the buckling is favored in **u**₂ direction by applying a vertical force F₂ = 10⁻⁶F at the center of the beam during loading.

The critical force for this example is known analytically to be $F_{crit} = \frac{4EI\pi^2}{L^2} = 5.2 \ 10^5$ N (Euler load). The numerical experiment matches well with this expected value: on the left equilibrium curve of Fig. 8, it corresponds to the first change of slope. The post-buckling is usually studied in-plane and the blue equilibrium curve obtained. However, there exists an out-of-plane solution more stable than the in-plane one. The ability of the ANM to switch from one branch of solution to an other allows to get this solution effortlessly. The bifurcation occurs at $F = 7.3 \ 10^5$ N which is visible on both equilibrium curves. The corresponding deformed shape right after bifurcation is exposed on Fig. 9. The three following deformed shapes show the evolution of this solution when pursuing the loading. An interesting phenomenon appears: a loop is formed and then rotates on itself, which physically translates as a coupling between torsion and bending. This illustrates the ability of the beam model to take into account complex phenomena.

### 5.6. Extremely twisted elastic rod

This example, more unusual, stems from an article of Lazarus et al. [21] for which an experimental procedure has been devised. For this experiment, two types of elastomeric rods have been fabricated under thoroughly controlled conditions. Thanks to this procedure, all the geometrical and material properties of the rods are accurately known. In particular, the two rods of same length 300 mm have been manufactured such that one is naturally straight and the other has a constant initial (natural) curvature $κ₃ = 44.84 \ m^{-1}$, forming the rolls at the top of Fig. 10. The other parameters of the experiment are the following: a circular cross-section of radius R = 1.55 mm, a Young's modulus E = 1300 kPa and a shear modulus G = 433 kPa. The setup, illustrated

on Fig. 10, consists in clamping both ends of the rod in two horizontally aligned chucks separated by a distance of 220 mm. Then, the rotation angle of the right-end is quasi-statically increased and the deformed shape observed by a camera.

Simulation-wise, the preparation of the setup requires preliminary loading steps (see Fig. 10). The curved rod is in a first step unrolled by applying a pure moment at one extremity while the other is clamped: as explained in paragraph 5.1, the analytical solution is known for this case and we can thus apply the predetermined moment such that the rod is straight at the end of the step. As a second step, a displacement in the axis of the rod $u₁ = -80$ mm is applied to the right-end of both rods. During this step, the gravity load is taken into account such that the beam buckles in the direction $-u₂$ and takes the same V-shape as in the real experiment. At this point, both rods are in the same configuration and solely the prestress differs between them because of the different initial curvature. The main step then consists in applying a rotation angle $θ₁$ using the special procedure described in paragraph 4.4.

The deformed shapes obtained for the initially straight rod (Fig. 11) and for the initially curved rod (Fig. 12) match perfectly the experimental observations of [21]. One observes that both deformed shapes are very different. In the first case, a plectoneme (loop) forms at the center of the rod after $θ = 11.3π$. In the second case, the rod takes a wavy shape with one twist per wave from $2.2π$ to $16π$. At the latter angle, a plectoneme also forms but this time on the right end. These qualitative observations illustrate the paramount importance of the initial curvature on the behavior of the beam.

This example also validates our method for rotation control and shows its ability to catch a torsion-bending coupling with ease. Besides, both equilibrium curves exhibited Fig. 13 are very intricate with a lot of turning points and unstable branches. As experimentally exposed by Lazarus et al. [21], the branch indeed becomes unstable after the turning point at $θ = 11.3π$ in the first case and is unstable from 0 to $2.2π$ in the second case. In spite of that, thanks to the ANM, the continuation of the branch is realized automatically and without efforts for the user.

Finally, let us highlight the fact that this experiment as a whole confronts our numerical solution to a real test case (with nevertheless thoroughly controlled conditions) and shows a very good agreement between the simulation and the reality.

<figure><img src="image_21.png" /><figcaption>Fig. 15. Numerical simulation results on an industrial example. Evolution of the deformed shape during loading.</figcaption></figure>
## 6. Application on a practical example

The examples of the previous section give a good overview of the robustness and accuracy of the herein presented method. To show the potential of the method in an industrial environment, we henceforth present an application to a real example. In this example pictured on Fig. 14, a cable harness made of a main fixed cable which first divides itself in two parts themselves fixed and then in 3 movable pieces is studied. The mounting operation consists in plugging these 3 cables to 3 connectors of known positions and orientations. The needs for the designers are to know accurately the final shape of the cable, the trajectory during the assembly, the stress inside the cable (to avoid damaging) and the feasibility of the operation.

The first step for the simulation is to build the geometry of the initial configuration. The 3 fixed parts of the cable harness are included in the geometry here only for illustration purposes but could be removed without affecting the results. Thanks to the choice of the FEM (see paragraph 3.3), the building of the complex geometry of the cable (6 pieces) is very easy and does not need any further efforts. For this simulation, the 3 movable cables are supposed to be initially straight. Then, for the end of the 3 cables, the final positions and orientations are set from CAD data. The loading thus consists in displacing and rotating linearly these ends from their initial ($\lambda = 0$) to their final configuration ($\lambda = 1$).

The shape of the cable during loading as well as the stress inside are displayed on Fig. 15. Of course, these data are known during the whole loading process and the trajectory of the cable can be retrieved from the code. Here, an arbitrary material is chosen but when the right material is defined, the simulation gives access to the right reaction force at the ends of the 3 cables. Finally, let us notice that the initial shape guessed in the CAD environment is not far from the “physical” shape obtained with our code. However, on the left cable of Fig. 14, the abrupt change of curvature at the joint seems unrealistic on the contrary

of our results, and could lead to unpredicted problems. In addition, this type of mockup of course does not give access to the physical data previously mentioned, i.e. inside stress, reaction force ...

## 7. Conclusion

In this paper, a finite element formulation of the geometrically exact beam model using rotational quaternions is presented. Their singularity-free property, their numerical efficiency and their algebraic nature are the main arguments for this choice. In particular, the latter trait leads to a polynomial formulation of the discrete equilibrium equations that we exploit through the use of the asymptotic numerical method. This method is a powerful solver for systems of quadratic non-linear equations which, among a lot of benefits, is fully automatic and is very robust since it gives a semi-analytical representation of the equilibrium branches. It replaces the classical predictor-corrector methods and the combination of the finite element method, the quaternions and the asymptotic numerical method then forms a remarkable and innovative tool for the simulation of nonlinear beams. Clues for using this method as an industrial tool for the simulation of cables in automotive industry are given, namely the assembly of several pieces of cables and the displacement/rotation control of the algorithm. Application of the method on several numerical and experimental examples exhibits great robustness and accuracy, which makes it very adapted for industry.

Further developments will consist in testing the method on a great variety of cables. In particular, the determination of material parameters which are inputs of the model appears to be a great challenge, especially for complex bundle-type structures. Refined constitutive laws may need to be devised. Finally, adding a stability studying tool to the numerical code will be a necessary step, in order to better comprehend the real phenomena.

# Appendix A. Strain measures in quaternion algebra

Let us prove here the expressions of the measures of strain in quaternion space. For the rotational strain measure, we need to come back to its definition and transpose it to quaternion space. Let us first define the local curvature $\kappa_L$ which stems from Frenet-Serret formulas through

$$ \forall i \in \{1, 2, 3\} \qquad \mathbf{e}'_i = \kappa_L \times \mathbf{e}_i. \qquad (A.1) $$

The expression of this curvature in the global frame, denoted $\kappa_G$, is obtained by mapping back this expression from the material to the global frame and then expressing it with the global coordinates which writes

$$ \mathbf{R}^T \mathbf{e}'_i = \kappa_G \times \mathbf{u}_i. \qquad (A.2) $$

Extending this definition in the quaternion space leads to

$$ \hat{\mathbf{q}}^* \circ \hat{\mathbf{e}}'_i \circ \hat{\mathbf{q}} = \hat{\kappa}_G \times \hat{\mathbf{u}}_i. \qquad (A.3) $$

Knowing the relation between the global and the local frame in quaternion space $\hat{\mathbf{e}}_i = \hat{\mathbf{q}} \circ \hat{\mathbf{u}}_i \circ \hat{\mathbf{q}}^*$, deriving it and finally using the relation (25) we get

$$ \begin{aligned} \hat{\mathbf{q}}^* \circ \hat{\mathbf{e}}'_i \circ \hat{\mathbf{q}} &= \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}' \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{e}}_i \circ \hat{\mathbf{q}} + \hat{\mathbf{q}}^* \circ \hat{\mathbf{e}}_i \circ \hat{\mathbf{q}} \circ \hat{\mathbf{q}}'^* \circ \hat{\mathbf{q}} \\ &= \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}' \circ \hat{\mathbf{u}}_i + \hat{\mathbf{u}}_i \circ \hat{\mathbf{q}}'^* \circ \hat{\mathbf{q}}. \end{aligned} \qquad (A.4) $$ $$ \begin{aligned} &= \hat{\mathbf{u}}_i \times (2\hat{\mathbf{q}}'^* \circ \hat{\mathbf{q}}), \\ &= (2\hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}') \times \hat{\mathbf{u}}_i. \end{aligned} $$

Consequently, the curvature in global coordinates writes

$$ \hat{\kappa}_G = 2\hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}'. \qquad (A.5) $$

As the rotational strain measure must equate 0 in the reference configuration, it is then obvious to define the strain measure as

$$ \hat{\mathbf{K}} = 2\hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}' - \hat{\mathbf{k}}_G^0. \qquad (A.6) $$
with $\hat{\kappa}_G^0 = 2\hat{q}_0^* \circ \hat{q}'_0$ the intrinsic curvature of the beam.

The expression of the translational strain measure is much more easier to obtain and derives simply from the classic expression (4). It suffices to recall that the rotation in terms of quaternions is defined by (20) and thus replacing $R^T$ by it gives directly

$$ \hat{\Gamma} = \hat{q}^* \circ \hat{x}'_0 \circ \hat{q} - \hat{\Gamma}_0. \qquad (A.7) $$

## Appendix B. Virtual strains in quaternion space

Let us prove here the expressions in quaternion space of the virtual strains (27a) and (27b). For the translational strain, let us on one hand differentiate simply the expression of the strain (26a) from which we get

$$ \delta \hat{\Gamma} = \delta \hat{q}^* \circ \hat{x}'_0 \circ \hat{q} + \hat{q}^* \circ \delta \hat{x}'_0 \circ \hat{q} + \hat{q}^* \circ \hat{x}'_0 \circ \delta \hat{q}. \qquad (B.1) $$

On the other hand, using the relation (25) allows to develop the final expression (27a) in

$$ \delta \hat{\Gamma} = \hat{q}^* \circ \delta \hat{x}'_0 \circ \hat{q} + \hat{q}^* \circ (\hat{x}'_0 \circ \delta \hat{q} \circ \hat{q}^* + \hat{q} \circ \delta \hat{q}^* \circ \hat{x}'_0) \circ \hat{q}. \qquad (B.2) $$

Finally using the relations (22) in (B.2), we recover the expression (B.1) and in that way prove the final expression (27a).

We proceed similarly to get the virtual rotational strain: we first differentiate the expression of the strain (26b) and get

$$ \delta \hat{K} = 2 \delta \hat{q}^* \circ \hat{q}' + 2 \hat{q}^* \circ \delta \hat{q}'. \qquad (B.3) $$

Then, developing the final expression (27b), we have

$$ \begin{aligned} \delta \hat{K} &= 2\hat{q}^* \circ (\delta \hat{q}' \circ \hat{q}^* + \delta \hat{q} \circ \hat{q}'^*) \circ \hat{q}, \\ &= 2\hat{q}^* \circ \delta \hat{q}' + 2\hat{q}^* \circ \delta \hat{q} \circ \hat{q}'^* \circ \hat{q}. \end{aligned} \qquad (B.4) $$

Finally, using relations (23) twice for the second term, we get

$$ \begin{aligned} \delta \hat{K} &= 2\hat{q}^* \circ \delta \hat{q}' - 2\hat{q}^* \circ \delta \hat{q} \circ \hat{q}^* \circ \hat{q}', \\ &= 2\hat{q}^* \circ \delta \hat{q}' + 2\delta \hat{q}^* \circ \hat{q}'. \end{aligned} \qquad (B.5) $$

and thus recover expression (B.3).

## Appendix C. Strong form of equilibrium equations in quaternion space

Let us demonstrate here how to deduce the strong form of equilibrium equations (31) from the weak form (29)–(30). Firstly, let us notice that the differential form of equations (23) writes

$$ \begin{aligned} \delta \hat{q} \circ \hat{q}^* &= -(\hat{q} \circ \delta \hat{q}^*), \\ \delta \hat{q}^* \circ \hat{q} &= -(\hat{q}^* \circ \delta \hat{q}). \end{aligned} \qquad (C.1) $$

Then, from the direct differentiation of the virtual translational strain $\delta \hat{q}^* \circ \hat{x}'_0 \circ \hat{q} + \hat{q}^* \circ \delta \hat{x}'_0 \circ \hat{q} + \hat{q}^* \circ \hat{x}'_0 \circ \delta \hat{q}$, and by using several times the relations (24) and (C.1) we show that the left part of the internal work is equal to

$$ \begin{aligned} \hat{N} \cdot \delta \hat{\Gamma} &= (\hat{N} \circ \hat{q}^*) \cdot (\hat{q}^* \circ \delta \hat{x}'_0) + (\hat{x}'_0^* \circ \hat{q} \circ \hat{N}) \cdot \delta \hat{q} + (\hat{N} \circ \hat{q}^* \circ \hat{x}'_0) \cdot \delta \hat{q}^*, \\ &= (\hat{q} \circ \hat{N} \circ \hat{q}^*) \cdot \hat{x}'_0 + (\hat{x}'_0^* \circ \hat{q} \circ \hat{N}) \cdot \delta \hat{q} + (\hat{N} \circ \hat{q}^* \circ \hat{x}'_0) \cdot (-\hat{q}'^* \circ \delta \hat{q} \circ \hat{q}^*), \\ &= (\hat{q} \circ \hat{N} \circ \hat{q}^*) \cdot \hat{x}'_0 + (\hat{x}'_0^* \circ \hat{q} \circ \hat{N}) \cdot \delta \hat{q} - (\hat{q} \circ \hat{N} \circ \hat{q}^* \circ \hat{x}'_0^*) \cdot (\delta \hat{q} \circ \hat{q}^*), \\ &= (\hat{q} \circ \hat{N} \circ \hat{q}^*) \cdot \hat{x}'_0 + (\hat{x}'_0^* \circ \hat{q} \circ \hat{N} - \hat{q} \circ \hat{N} \circ \hat{q}^* \circ \hat{x}'_0^* \circ \hat{q}^*) \cdot \delta \hat{q}. \end{aligned} \qquad (C.2) $$

Then factorizing on the right by $\hat{q}$ and using that, as a pure quaternion, $\hat{x}'_0$ meets $\hat{x}'_0^* = -\hat{x}'_0$, we have

$$ \hat{N} \cdot \delta \hat{\Gamma} = (\hat{q} \circ \hat{N} \circ \hat{q}^*) \cdot \delta \hat{x}'_0 + ((\hat{x}'_0^* \circ \hat{q} \circ \hat{N} \circ \hat{q}^* + \hat{q} \circ \hat{N} \circ \hat{q}^* \circ \hat{x}'_0) \circ \hat{q}) \cdot \delta \hat{q}. \qquad (C.3) $$

Finally, with relation (25), we get

$$ \hat{N} \cdot \delta \hat{\Gamma} = (\hat{q} \circ \hat{N} \circ \hat{q}^*) \cdot \delta \hat{x}'_0 + ((2(\hat{q} \circ \hat{N} \circ \hat{q}^*) \times \hat{x}'_0) \circ \hat{q}) \cdot \delta \hat{q}. \qquad (C.4) $$

Besides, still using relations (24) several times and the expression (23) in the direct differentiation of (27b) $\delta \hat{K} = 2\delta \hat{q}^* \circ \hat{q}' + 2\hat{q}^* \circ \delta \hat{q}'$, the right part of the internal work develops as

$$ \begin{aligned} \hat{M} \cdot \delta \hat{K} &= 2(\hat{q} \circ \hat{M}) \cdot \delta \hat{q}' - 2\hat{M} \cdot (\hat{q}^* \circ \delta \hat{q} \circ \hat{q}^* \circ \hat{q}'), \\ &= 2(\hat{q} \circ \hat{M}) \cdot \delta \hat{q}' - 2(\hat{M} \circ \hat{q}'^* \circ \hat{q}) \cdot (\hat{q}^* \circ \delta \hat{q}), \\ &= 2(\hat{q} \circ \hat{M}) \cdot \delta \hat{q}' - 2(\hat{q} \circ \hat{M} \circ \hat{q}'^* \circ \hat{q}) \cdot \delta \hat{q}, \\ &= 2(\hat{q} \circ \hat{M}) \cdot \delta \hat{q}' + 2(\hat{q} \circ \hat{M} \circ \hat{q}^* \circ \hat{q}') \cdot \delta \hat{q}. \end{aligned} \qquad (C.5) $$
Then, integrating by parts the virtual work of internal forces in which expressions (C.4) and (C.5) have been inserted we get

$$ \int_{0}^{L} (\hat{\mathbf{N}} \cdot \delta\hat{\mathbf{\Gamma}} + \hat{\mathbf{M}} \cdot \delta\hat{\mathbf{K}}) \, ds = [ (\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*)' \cdot \delta\hat{\mathbf{x}}_0 + 2(\hat{\mathbf{q}} \circ \hat{\mathbf{M}}) \cdot \delta\hat{\mathbf{q}} ]_{0}^{L} - \int_{0}^{L} \left( (\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*)' \cdot \delta\hat{\mathbf{x}}_0 + 2 (\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^*)' + \hat{\mathbf{x}}'_0 \times (\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*) ) \circ \hat{\mathbf{q}} \right) \cdot \delta\hat{\mathbf{q}} \, ds. \quad (C.6) $$

Finally, inserting this expression in the weak form (29)–(30) and using the fundamental lemma of calculus of variations leads to

$$ \left\{ \begin{array}{l} \mathrm{Vec}[(\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*)'] + \mathbf{n}_e = \mathbf{0}, \\[1.5ex] \left( (\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^*)' + \hat{\mathbf{x}}'_0 \times (\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*) + \hat{\mathbf{m}}_e - \mu \hat{\mathbf{1}} \right) \circ \hat{\mathbf{q}} = \hat{\mathbf{0}}, \quad \forall s \in [0, L] \\[1.5ex] \hat{\mathbf{q}} \cdot \hat{\mathbf{q}} - 1 = 0, \end{array} \right. \quad (C.7) $$

along with

$$ \left\{ \begin{array}{l} \mathrm{Vec}[(\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*)'] - \mathbf{N}_e = \mathbf{0} \quad \text{at } s=0, L, \\[1.5ex] \left( (\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^*) - \hat{\mathbf{M}}_e \right) \circ \hat{\mathbf{q}} = \hat{\mathbf{0}} \quad \text{at } s=0, L. \end{array} \right. \quad (C.8) $$

As $|\hat{\mathbf{q}}| \neq 0$, multiplying the second equation of (C.7) on the right by $\hat{\mathbf{q}}^*$ leads to the fully equivalent equation:

$$ (\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^*)' + \hat{\mathbf{x}}'_0 \times (\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*) + \hat{\mathbf{m}}_e - \mu \hat{\mathbf{1}} = \hat{\mathbf{0}}. \quad (C.9) $$

It is noteworthy that the term $\mu \hat{\mathbf{1}}$ in (C.9) is a scalar quaternion while the other terms of this equation are pure quaternions. It results that this equation is fully equivalent to $\mu(s) = 0$ and $\mathrm{Vec}[(\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^*)' + \hat{\mathbf{x}}'_0 \times (\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*) + \hat{\mathbf{m}}_e] = \mathbf{0}$, leading to the strong form of equilibrium equations (31). Similarly, the second equation of (C.8) may be multiplied on the right by $\hat{\mathbf{q}}^*$, leading to the final form of natural boundary conditions (32).

## Appendix D. Expression of the virtual rotational vector in terms of quaternions

The virtual rotational vector $\delta\theta$ is defined as the rate of change of the material basis expressed in the global frame [12] which means

$$ \mathbf{R}^T \delta \mathbf{e}_i = \delta \theta \times \hat{\mathbf{u}}_i. \quad (D.1) $$

This expression is fully analogous to the definition of the global curvature. Expressed in terms of quaternions, it then writes

$$ \hat{\mathbf{q}}^* \circ \delta \hat{\mathbf{e}}_i \circ \hat{\mathbf{q}} = \delta \theta \times \hat{\mathbf{u}}_i. \quad (D.2) $$

Using the relation between the material and the global frame as well as the equation (25), we demonstrate

$$ \begin{aligned} \hat{\mathbf{q}}^* \circ \delta \hat{\mathbf{e}}_i \circ \hat{\mathbf{q}} &= \hat{\mathbf{q}}^* \circ \delta \hat{\mathbf{q}} \circ \hat{\mathbf{u}}_i + \hat{\mathbf{u}}_i \circ \delta \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}} \\ &= \hat{\mathbf{u}}_i \times (2\delta \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}) \\ &= (2\hat{\mathbf{q}}^* \circ \delta \hat{\mathbf{q}}) \times \hat{\mathbf{u}}_i. \end{aligned} \quad (D.3) $$

We thereby deduce the relation (28): $\delta\hat{\theta} = 2\hat{\mathbf{q}}^* \circ \delta\hat{\mathbf{q}}$.

## Appendix E. Discrete expression of the virtual strains

Using the results obtained in Appendix C, we have

$$ \hat{\mathbf{N}} \cdot \delta\hat{\mathbf{\Gamma}} + \hat{\mathbf{M}} \cdot \delta\hat{\mathbf{K}} = (\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*) \cdot \delta\hat{\mathbf{x}}'_0 + 2(\hat{\mathbf{q}} \circ \hat{\mathbf{M}}) \cdot \delta\hat{\mathbf{q}}' + (\hat{\mathbf{x}}'_0^* \circ \hat{\mathbf{q}} \circ \hat{\mathbf{N}} + \hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{x}}'_0 \circ \hat{\mathbf{q}} + 2(\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^*) \circ \hat{\mathbf{q}}') \cdot \delta\hat{\mathbf{q}}. \quad (E.1) $$

Now developing the third component of (44) with equations (26a), (26b) and (25), we show that

$$ \begin{aligned} & 2\hat{\mathbf{q}} \circ (\hat{\mathbf{N}} \times (\hat{\mathbf{\Gamma}} + \hat{\mathbf{\Gamma}}_0) + \hat{\mathbf{M}} \times (\hat{\mathbf{K}} + \hat{\mathbf{K}}_0)) + 2\hat{\mathbf{q}}' \circ \hat{\mathbf{M}} \\ &= 2\hat{\mathbf{q}} \circ (\hat{\mathbf{N}} \times (\hat{\mathbf{q}}^* \circ \hat{\mathbf{x}}'_0 \circ \hat{\mathbf{q}})) + 2\hat{\mathbf{q}} \circ (\hat{\mathbf{M}} \times (2\hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}')) + 2\hat{\mathbf{q}}' \circ \hat{\mathbf{M}} \\ &= \hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{x}}'_0 \circ \hat{\mathbf{q}} + \hat{\mathbf{q}} \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{x}}'_0^* \circ \hat{\mathbf{q}} \circ \hat{\mathbf{N}} + 2\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}' + 2\hat{\mathbf{q}} \circ \hat{\mathbf{q}}'^* \circ \hat{\mathbf{q}} \circ \hat{\mathbf{M}} + 2\hat{\mathbf{q}}' \circ \hat{\mathbf{M}} \\ &= \hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{x}}'_0 \circ \hat{\mathbf{q}} + \hat{\mathbf{x}}'_0^* \circ \hat{\mathbf{q}} \circ \hat{\mathbf{N}} + 2\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}' - 2\hat{\mathbf{q}}' \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}} \circ \hat{\mathbf{M}} + 2\hat{\mathbf{q}}' \circ \hat{\mathbf{M}} \\ &= \hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{x}}'_0 \circ \hat{\mathbf{q}} + \hat{\mathbf{x}}'_0^* \circ \hat{\mathbf{q}} \circ \hat{\mathbf{N}} + 2\hat{\mathbf{q}} \circ \hat{\mathbf{M}} \circ \hat{\mathbf{q}}^* \circ \hat{\mathbf{q}}'. \end{aligned} \quad (E.2) $$

We thus recover the term factor of $\delta\hat{\mathbf{q}}$ in equation (E.1). Moreover, stating that $\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*$ and $\delta\hat{\mathbf{x}}'_0$ are pure quaternions, we can write

$$ (\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*) \cdot \delta\hat{\mathbf{x}}'_0 = \mathrm{Vec}(\hat{\mathbf{q}} \circ \hat{\mathbf{N}} \circ \hat{\mathbf{q}}^*) \cdot \delta\mathbf{x}'_0. \quad (E.3) $$
Finally, with the vector $g(N, M)$ defined in (44), we can rewrite (E.1) as

$$ \hat{N} \cdot \delta\hat{\Gamma} + \hat{M} \cdot \delta\hat{K} = g(N, M) \cdot \begin{bmatrix} \delta x'_{0} \\ \delta \hat{q}' \\ \delta \hat{q} \\ \delta \mu \end{bmatrix}, \qquad (E.4) $$

and with equation (41) it results the final relation (44).

## Appendix F. Expression of the rotation matrix with quaternions

Using the definition of the quaternion multiplication (17), one can express left and right quaternion multiplication under matrix form. They respectively write

$$ \hat{a} \circ \hat{b} = \Phi_L(\hat{a})\hat{b} = \begin{bmatrix} a_0 & -a^T \\ a & a_0I_3 + \tilde{a} \end{bmatrix} \begin{bmatrix} b_0 \\ b \end{bmatrix} \qquad (F.1) $$

and

$$ \hat{a} \circ \hat{b} = \Phi_R(\hat{b})\hat{a} = \begin{bmatrix} b_0 & -b^T \\ b & b_0I_3 - \tilde{b} \end{bmatrix} \begin{bmatrix} a_0 \\ a \end{bmatrix}. \qquad (F.2) $$

It is then trivial to express the rotation operation (20) under matrix form, that being

$$ \hat{y} = \hat{q} \circ \hat{x} \circ \hat{q}^* = \Phi_R(\hat{q}^*) \Phi_L(\hat{q}) \hat{x} = \begin{bmatrix} 1 & 0^T \\ 0 & R \end{bmatrix} \hat{x}. \qquad (F.3) $$

The matrix $R$ in (F.3) is the classical rotation matrix expressed with quaternion parameters through expression

$$ R = q_0^2 I_3 + 2q_0 \tilde{q} + \tilde{q} \tilde{q} + qq^T. \qquad (F.4) $$

The equivalence between quaternion operation and the classical rotation operator is then exhibited clearly under this form (see also Appendix B of [10] for an other proof).

## References

[1] A. Shabana, H. Hussien, J. Escalona, Application of the absolute nodal coordinate formulation to large rotation and large deformation problems, J. Mech. Des. 120 (2) (1998) 188–195.

[2] G. Hegemier, S. Nair, A nonlinear dynamical theory for heterogeneous, anisotropic, elastic rods, AIAA J. 15 (1) (1977) 8–15.

[3] D.H. Hodges, A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams, Int. J. Solids Struct. 26 (11) (1990) 1253–1273.

[4] E. Reissner, On one-dimensional finite-strain beam theory: the plane problem, Z. für Angew. Math. Phys. ZAMP 23 (5) (1972) 795–804.

[5] E. Reissner, On finite deformations of space-curved beams, Z. für Angew. Math. Phys. ZAMP 32 (6) (1981) 734–744.

[6] J.C. Simo, A finite strain beam formulation. the three-dimensional dynamic problem. part i, Comput. Methods Appl. Mech. Eng. 49 (1) (1985) 55–70.

[7] G. Jelenić, M. Crisfield, Geometrically exact 3d beam theory: implementation of a strain-invariant finite element for statics and dynamics, Comput. Methods Appl. Mech. Eng. 171 (1) (1999) 141–171.

[8] A. Cardona, M. Geradin, A beam finite element non-linear theory with finite rotations, Int. J. Numer. Methods Eng. 26 (11) (1988) 2403–2438.

[9] A. Ibrahimbegović, F. Frey, I. Kožar, Computational aspects of vector-like parametrization of three-dimensional finite rotations, Int. J. Numer. Methods Eng. 38 (21) (1995) 3653–3673.

[10] J.C. Simo, L. Vu-Quoc, A three-dimensional finite-strain rod model. part ii: computational aspects, Comput. Methods Appl. Mech. Eng. 58 (1) (1986) 79–116.

[11] M. Ritto-Corréa, D. Camotim, On the differentiation of the Rodrigues formula and its significance for the vector-like parameterization of Reissner–Simo beam theory, Int. J. Numer. Methods Eng. 55 (9) (2002) 1005–1032.

[12] E. Zupan, M. Saje, D. Zupan, The quaternion-based three-dimensional beam theory, Comput. Methods Appl. Mech. Eng. 198 (49) (2009) 3944–3956.

[13] S. Ghosh, D. Roy, Consistent quaternion interpolation for objective finite element approximation of geometrically exact beam, Comput. Methods Appl. Mech. Eng. 198 (3) (2008) 555–571.

[14] J. Stuelpnagel, On the parametrization of the three-dimensional rotation group, SIAM Rev. 6 (4) (1964) 422–430.

[15] E. Zupan, M. Saje, D. Zupan, On a virtual work consistent three-dimensional Reissner–Simo beam formulation using the quaternion algebra, Acta Mech. 224 (8) (2013) 1709–1729.

[16] M. Criesfield, Non Linear Finite Element Analysis of Solids and Structures, 1, Wiley, New York, 1991.

[17] E. Riks, An incremental approach to the solution of snapping and buckling problems, Int. J. Solids Struct. 15 (7) (1979) 529–551.

[18] N. Damil, M. Potier-Ferry, A new method to compute perturbed bifurcations: application to the buckling of imperfect elastic structures, Int. J. Eng. Sci. 28 (9) (1990) 943–957.

[19] B. Cochelin, A path-following technique via an asymptotic-numerical method, Comput. Struct. 53 (5) (1994) 1181–1192.

[20] B. Cochelin, M. Medale, Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations, J. Comput. Phys. 236 (2013) 594–607.

[21] A. Lazarus, J. Miller, P. Reis, Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method, J. Mech. Phys. Solids 61 (8) (2013) 1712–1736.

[22] M. Géradin, A. Cardona, Flexible Multibody Dynamics: a Finite Element Approach, Wiley, 2001.

[23] M.A. Crisfield, G. Jelenić, Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation, in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 455, The Royal Society, 1999, pp. 1125–1147.

[24] E. Salamin, Application of Quaternions to Computation with Rotations, Tech. rep., Working Paper, 1979.

[25] J. Vince, Quaternions for Computer Graphics, Springer Science & Business Media, 2011.

[26] C. Bottasso, A non-linear beam space-time finite element formulation using quaternion algebra: interpolation of the Lagrange multipliers and the appearance of spurious modes, Comput. Mech. 10 (5) (1992) 359–368.

[27] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method: Solid Mechanics, 2, Butterworth-heinemann, 2000.

[28] M. Oak-Key, K. Yong-Woo, Reduced minimization theory in beam elements, Int. J. Numer. Methods Eng. 37 (12) (1994) 2125–2145.

[29] A.A. Shabana, Dynamics of Multibody Systems, Cambridge university press, 2013.

[30] A. Ibrahimbegovic, R.L Taylor, On the role of frame-invariance in structural mechanics models at finite rotations, Comput. Methods Appl. Mech. Eng. 191 (45) (2002) 5159–5176.

[31] I. Romero, The interpolation of rotations and its application to finite element models of geometrically exact rods, Comput. Mech. 34 (2) (2004) 121–133.

[32] O.A. Bauchau, A. Epple, S. Heo, Interpolation of finite rotations in flexible multi-body dynamics simulations, Proc. Inst. Mech. Eng. Part K: J. Multi-body Dyn. 222 (4) (2008) 353–366.

[33] E. Zupan, D. Zupan, On higher order integration of angular velocities using quaternions, Mech. Res. Commun. 55 (2014) 77–85.

[34] D. Zupan, M. Saje, Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures, Comput. Methods Appl. Mech. Eng. 192 (49) (2003) 5209–5248.
[35] M. Schulz, F.C. Filippou, Non-linear spatial Timoshenko beam element with
curvature interpolation, Int. J. Numer. Methods Eng. 50 (4) (2001) 761-785.

[36] B. Cochelin, N. Damil, M. Potier-Ferry, Méthode Asymptotique Numérique,
Hermes Lavoissier, 2007.

[37] S. Karkar, B. Cochelin, C. Vergez, A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems, J. Sound Vib. 333 (12) (2014) 2554-2567.

[38] S. Karkar, Méthodes numériques pour les systèmes dynamiques non linéaires: application aux instruments de musique auto-oscillants (Ph.D. thesis), Aix-Marseille, 2012.

[39] S. Karkar, R. Arquier, A. Lazarus, O. Thomas, C. Vergez, B. Cochelin, Manlab - an Interactive Path-following and Bifurcation Analysis Software, 2011. http://manlab.lma.cnrs-mrs.fr/.

[40] M. Geradin, A. Cardona, Kinematics and dynamics of rigid and flexible
mechanisms using finite elements and quaternion algebra, Comput. Mech. 4 (2)
(1988) 115-135.

[41] K.-J. Bathe, S. Bolourchi, Large displacement analysis of three-dimensional beam structures, Int. J. Numer. Methods Eng. 14 (7) (1979) 961-986.

[42] P. Češarek, M. Saje, D. Zupan, Kinematically exact curved and twisted strain-based beam, Int. J. Solids Struct. 49 (13) (2012) 1802-1817.

[43] D. DaDeppo, R. Schmidt, Instability of clamped-hinged circular arches subjected to a point load, J. Appl. Mech. 42 (4) (1975) 894-896.

[44] J. Argyris, S. Symeonidis, Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. part i. quasistatic problems, Comput. Methods Appl. Mech. Eng. 26 (1) (1981) 75-123.