# A special finite element for static and dynamic study of mechanical systems under large motion, Part 2

Mircea Munteanu, André Barraco

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

## To cite this version:

Mircea Munteanu, André Barraco. A special finite element for static and dynamic study of mechanical systems under large motion, Part 2. Revue Européenne des Éléments Finis, 2002, 11 (6), pp.791-814. 10.3166/reef.11.791-814. hal-00021125

HAL Id: hal-00021125

https://hal.science/hal-00021125v1

Submitted on 15 Apr 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Special Finite Element for Static and Dynamic Study of Mechanical Systems under Large Motion, Part 2

Mircea Gh. Munteanu* — André Barraco**

* Université Transilvania de Brasov, Roumanie
29, B-dul Eroilor, 2200 Brasov, Roumanie
mmm@deltanet.ro

** LM2S, CNRS UPRES A 8007
ENSAM, 151 Bd de L'Hôpital
F-75013 Paris
barraco@paris.ensam.fr

**ABSTRACT.** In the first part of the paper the theory of the 3D dynamics of mechanical systems composed by elastic beams, structures and mechanisms, was studied. These systems are divided into so-called macro-elements and the movement equations of one macro-element were established. Only the Euler-Rodrigues parameters are used to describe the global motion of the system. In this second part of the paper a special finite element (SFET) having four degrees of freedom per node, the Euler-Rodrigues parameters, is described in details. The stiffness and mass matrices are expressed only in nodal Euler-Rodrigues parameters. The most important aspect of the proposed approach is that the exact equations, written for the deformed configuration, are solved. Therefore an extremely accurate and very fast convergent method results. To validate the SFET finite element finally several 2D and 3D, static and dynamic examples are presented and the accuracy of the results is discussed.

**RÉSUMÉ.** Dans la première partie de ce papier nous présentions l'approche théorique de l'étude dynamique de systèmes, structures ou mécanismes, composés de poutres élastiques. Ces systèmes étaient décomposés en macro-éléments et les équations du mouvement d'un macro-élément étaient établies. Dans cette seconde partie nous présentons, en détail, un élément fini spécial (SFET) qui a quatre degrés de liberté par nœud, les paramètres d'Euler-Rodrigues. Les matrices de raideur et de masse ne sont exprimées qu'en fonction de ces degrés de liberté. Le point le plus important à souligner est que les équations exactes, écrites sur la configuration déformée, sont résolues. Ainsi la méthode proposée est-elle très exacte et rapidement convergente. Des exemples en 2D et 3D, statiques et dynamiques sont présentés et démontrent la validité de l'élément SFET, la précision d'un tel élément est discutée et validée.

**KEYWORDS:** finite element method, large motion, mechanical systems.

**MOTS-CLÉS** : méthode des éléments finis, grand mouvement, systèmes mécaniques.
# 1. Introduction

In the first part of this paper, the equations describing the static and dynamic elastic behaviour of a mechanical system were developed. In this work we consider that “mechanical system” denotes a multi body system, in 2D or 3D, which may be structures or mechanisms, both involving large displacements.

The exact differential equations are written for the actual, deformed, configuration of the mechanical system and the equations are exact in the limits of the accepted hypotheses used to find out the constitutive equations. These equations are written in the global axes of the mechanical system. The non-linearity is provided mainly by the large displacements that could be so large that the initial geometry of the system is completely changed (geometrical non-linearity). This firstly means that the rotations are very large and they could not be considered as vectors anymore as in small displacement mechanical systems.

The basic idea of the method is that the actual configuration of any mechanical system might be uniquely described only by the rotations of some points called nodes, of course with the approximation of a rigid body motion. Because in 3D it is very complicated to work in rotations, the authors defined the rotation by using the quaternion or Euler-Rodrigues parameters. The constitutive equations and the variation of the total potential energy were expressed in these Euler-Rodrigues parameters. In the presented approach we try to find a quasi-exact solution of constitutive equations (in static or dynamic analysis) or, that is the same thing, to satisfy the variational principle, that is any geometrically admissible virtual variation of the displacements must lead to a zero total potential energy variation. The problem that we solve using a special finite element (called SFET in the paper) is the dynamic of a Bernoulli-Euler beam, the macro-element shown in the Figure 4 of the first part of this paper. In this second part the special finite element is described in details and finally several examples are presented.

# 2. Finite element approach

The proposed solution is based on the finite element method in space and the finite difference method in the time domain. An original special curvilinear finite element type (SFET) was elaborated with 4 degrees of freedom per node in 3D, the four component of Euler-Rodrigues quaternion, and only one degree of freedom per node in plane as the vector $n$ is known (perpendicular to the plane of the problem). The finite element is a curvilinear one and may have several nodes, but very good results were obtained using a two-nodes finite element and, especially, a 3-nodes element (Figure 1). The two-node finite element (SFET2) is very easy to use, but the three-node finite element (SFET3) is extremely accurate even if it is more difficult to program. The stiffness and mass matrices are found out starting from the variation of the total potential energy written for the actual configuration.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Two-node and three-node finite elements</figcaption></figure>

The shape function of the finite element is:

$$ \Box(s) \Box(1 - \frac{s}{h})\Box_i \Box \frac{s}{h}\Box_{i\Box 1} \quad [1] $$

where $\hat{l}_e$ represents the column matrix of nodal unknowns having a number of elements equal to four times the number of element nodes, that is 8 for SFET2 and 12 for SFET3.

SFET2, Figure 1, has the shape function for each Euler-Rodrigues parameter:

$$ \Box(s) \Box(1 - \frac{s}{h})\Box_i \Box \frac{s}{h}\Box_{i\Box 1} \quad [2] $$

that is:

$$ \Box(s) \Box(1 - \frac{s}{h})\Box_i \Box \frac{s}{h}\Box_{i\Box 1} \quad [3] $$

with the nodal unknowns matrix

$$ \hat{l}_e \Box \Box_{0,1} l_{1,1} l_{2,1} l_{3,1} l_{0,2} l_{1,2} l_{2,2} l_{3,2} \Box^T $$

and

$$ \Box(s) \Box(1 - \frac{s}{h})\Box_i \Box \frac{s}{h}\Box_{i\Box 1} $$

For SFET3, Figure 1, we can write:

$$ l_k(S_1) \Box \Box \frac{S_1^2}{2h^2} \Box \frac{S_1}{2h} \Box l_{k,i-1} \Box (\Box \frac{S_1^2}{h^2} \Box 1)l_{k,i} \Box (\Box \frac{S_1^2}{2h^2} \Box \frac{S_1}{2h})l_{k,i\Box 1}, \quad k \Box 0,1,2,3 \quad [4] $$

and in this case $N(S_1)$ matrix has the dimension 4x12:
$$ \Box(s) \Box (1 \Box \frac{s}{h}) \Box_i \Box \frac{s}{h} \Box_{i \Box 1} $$

The nodal unknowns of one finite element are grouped in the column matrix:

$$ \hat{l}_e \Box \Box_{0,1} l_{1,1} l_{2,1} l_{3,1} l_{0,2} l_{1,2} l_{2,2} l_{3,2} l_{0,3} l_{1,3} l_{2,3} l_{3,3} \Box^T $$
and

$$ N_1(S_1) \Box \Box \frac{S_1^2}{2h^2} \Box \frac{S_1}{2h} \Box; \quad N_2(S_1) \Box (\Box \frac{S_1^2}{h^2} \Box 1); \quad N_3(S_1) \Box (\Box \frac{S_1^2}{2h^2} \Box \frac{S_1}{2h}) $$

Now we have to establish a relation between the nodal unknowns $\hat{l}_e$ and the column matrix $\hat{l} \Box \frac{d}{ds} \hat{l}$ starting from the above shape functions:

$$ \hat{l} \Box N(s) \hat{l}_e \quad [6] $$

For SFET2 we get:

$$ \Box(s) \Box (1 \Box \frac{s}{h}) \Box_i \Box \frac{s}{h} \Box_{i \Box 1} $$
and then for SFET3:

$$ \Box(s) \Box (1 \Box \frac{s}{h}) \Box_i \Box \frac{s}{h} \Box_{i \Box 1} $$
where:

$$ N_1(s) \Box \Box \frac{s_1}{h^2} \Box \frac{1}{2h} \Box; \quad N_2(s) \Box \Box \frac{2s_1}{h^2}; \quad N_3(s) \Box (\frac{s_1}{h^2} \Box \frac{1}{2h}) $$

For a finite element, it is possible to write the variation of the curvature along the finite element:

$$ \Box_i \Box \hat{l}_e^T N(S_1) A_i N(S_1) \hat{l}_e \Box \frac{1}{2} \hat{l}_e^T N(S_1) A_i N(S_1) \Box N(S_1) A_i N(S_1) \Box \hat{l}_e \quad [7] $$

Or more simply:

$$ \Box_i \Box \frac{1}{2} \hat{l}_e^T T_i \Box s_1 \Box \hat{l}_e \quad [8] $$

where $\kappa_i$, $i=1...3$, represents the i-th element of the matrix $\kappa$. The matrices:

$$ T_i \Box s_1 \Box N(S_1) A_i N(S_1) \Box N(S_1) A_i N(S_1) $$
numerical integration method is exact for the polynomials of third degree or smaller, but very good results were obtained as will be shown later.

### 3. Static problem

The assembly process of the elemental stiffness matrices, $k_{el}$ into the macro-element stiffness matrix, $K$, follows the standard procedure. For the static problem the Equation [11] leads to:

$$ \hat{l}_s^T E \hat{l}_s \hat{l}_s^T K \hat{l}_s^T F \hat{l}_s = 0 \quad [13] $$

which has to be true for the any virtual values of nodal unknowns $\hat{l}_s$ geometrically admissible, that is the non-linear system of equations $E$ is got:

$$ E \hat{l}_s K \hat{l}_s F \hat{l}_s = 0 \quad [14] $$

The force vector $F$, which in fact contains concentrated moments reduced to the nodal Euler-Rodrigues parameters, is obtained from the last two terms in the Equation [11] and follows the standard procedure of assembling, too. In the above relation the vector of nodal unknowns for the whole macro-element was denoted $\hat{l}_s$. It is interesting to notice that the elements of the stiffness matrix $K$ are quadratic polynomials in the nodal unknowns $\hat{l}_s$.

The four nodal unknown for each node are not independent and therefore the secant stiffness matrix $K$ is singular. We have to impose the following condition for each node, that is:

$$ \hat{l}_j^T \hat{l}_j = 1, \quad j = 1, \dots, n_s \quad [15] $$

Obviously it results:

$$ \hat{l}_j^T \hat{l}_j = 0, \quad j = 1, \dots, n_s \quad [16] $$

From the Equation [13] we select only the terms referring to the node j

$$ \hat{l}_j^T E \hat{l}_j = \begin{array}{c} \hat{l}_{0,j} \\ \hat{l}_{1,j} \\ \hat{l}_{2,j} \\ \hat{l}_{3,j} \end{array} = \begin{array}{c} \hat{E}_{0,j} \\ \hat{E}_{1,j} \\ \hat{E}_{2,j} \\ \hat{E}_{3,j} \end{array} = 0 $$

we express $\hat{l}_{3,j}$ from the Equation [16] and we replace it in the above equation:
<figure><img src="image_4.png" /></figure>

One equation among each four equations was thus eliminated. Because to eliminate one nodal unknown, too, is quite impossible, we preferred to add the Equation [15]. Thus for each node we get the following four non-linear equations:

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

[18]

In order to avoid the singularities which could occur when $l_{3,j}=0$, in the Equations [18], at every iteration cycle and for every node, the maximum absolute value component of the nodal unknowns $\hat{l}_j$ is searched and the corresponding equation is eliminated as shown above and then replaced by the Equation [15].

The boundary conditions for the clamped end of the macro-element, point 0, are:

$$ \hat{l}_0 \square 1; \hat{l}_1 \square 0; \hat{l}_2 \square 0; \hat{l}_3 \square 0 \quad [19] $$

and are imposed in the usual conventional way.

A special mention has to be done for the 2D case. In this case $l_0 \square \cos \frac{\square}{2}$ and, because the normal **n** is perpendicular to the plane, it results that $l_1 \square 0, l_2 \square 0, l_3 \square \sin \frac{\square}{2}$. So, we could chose two ways to elaborate the finite element: (I) either to consider two degrees of freedom per nodes, $l_0$ and $l_3$, and all we presented above it is easy to particularise for the 2D case, or (II) we can chose an unique function that can describe completely the configuration of the deformed beam, the angle $\theta(S_1)$ which in this case is the slope of the tangent to the actual configuration in the point of curvilinear co-ordination $S_1$. Therefore the finite element will have only one degree of freedom. The first way has the advantage that avoids the calculation of trigonometric functions and use polynomials instead, but
leads to two degrees of freedom per node. In addition, it is easy to modify for 2D case the computer code written for the 3D case simply by eliminating the columns and rows corresponding to the $l_1$ and $l_2$ nodal unknowns that are known to be zero. The authors followed the two approaches and the results were identical.

To solve the non-linear system [14] the iterative Newton-Raphson method is applied:

$$ \hat{l}_s^{(i)} \triangleq \hat{l}_s^{(i)} \triangleq J^{\top} E \hat{l}_s^{(i)} $$

where $\hat{l}_s^{(i)}$ is the column matrix of the nodal unknowns after the i-th iteration and J is the Jacobian of the matrix $E, J \triangleq \frac{d}{d\hat{l}_s} E$, computed for the $\hat{l}_s^{(i)}$ values of the nodal unknowns. The method has proved to be extremely rapid convergent. It is rather easy to compute the Jacobian as the expressions E are polynomials in the nodal unknowns $\hat{l}_s^{(i)}$.

In the static domain the finite element type elaborated by the authors is a non-incremental one, usually only one load step is sufficient to reach the final value, but in any case the accuracy of the results does not depend on the number of load steps. This happens because the exact equations (obviously in the limits of the beam theory) were solved and they are written for the actual deformed configuration of the beam system, which is not possible for the classic beam finite element. The accuracy and the very rapid convergence might be explicated by the fact that using only Euler-Rodrigues parameters as nodal unknown the constitutive differential equations become simpler and their order smaller. In addition the obtained finite element equations are polynomials in the nodal unknowns $\hat{l}_s$.

## 4. Dynamic problem

In order to simplify the explanations we will consider that the masses are lumped at the nodes. Let's consider the mass at node j, $m_j$. On this mass, in the actual configuration, the inertia force $F_{i,j} \triangleq m_j \ddot{\mathbf{r}}_i$ acts and therefore we get:

$$ F_{i,j} = \frac{\ddot{\mathbf{r}}_i}{m_j} \triangleq \frac{\ddot{\mathbf{r}}_i}{\mathbf{r}_i^2} \triangleq \frac{\ddot{\mathbf{r}}_i(l_0^2 + l_1^2) - 1}{\mathbf{r}_i^2(l_0^2 + l_1^2)} \triangleq \frac{\ddot{\mathbf{r}}_i(l_0 l_2 + l_1 l_3) - 1}{\mathbf{r}_i^2(l_0 l_2 + l_1 l_3)} \triangleq \frac{\ddot{\mathbf{r}}_i(l_0 l_2 + l_1 l_3) - 1}{\mathbf{r}_i^2(l_0 l_2 + l_1 l_3)} \triangleq \mathbf{0} $$

We can write successively:
Some concentrated rotational inertia may exist in node, too. For instance in this way we could take into account the transversal dimension of the beam.

For instance let's consider that in the node j we have to add the rotational mass inertia given by:

$$ \begin{bmatrix} J_1 & 0 & 0 \\ 0 & J_2 & 0 \\ 0 & 0 & J_3 \end{bmatrix} $$

The above tensor is written in local reference frame $R_B(G, B_i)$. The concentrated rotational inertia is converted to the nodal unknowns, in the global axes system $R_e(O, e_1, e_2, e_3)$, according the relation:

$$ J_l \otimes G_j^T \otimes {}^T_{Bb_0,j} J \otimes {}^T_{Bb_0,j} G_j $$

and then assembled to the mass matrix in the conventional way.

Relations [25] and [26] seem to be very complicated, but from programming point of view it is very simply to write a code to compute the two mass matrices.

The variational formulation becomes for the finite element assembly:

$$ \begin{bmatrix} \hat{l}_s^T & M l_s & C & M_v & K \hat{l}_s & F(t, \hat{l}_s) \end{bmatrix} = 0 \quad [27] $$

for any geometrically admissible virtual Euler-Rodrigues nodal parameters. In order to impose the conditions [15] for every node, we will use a similar approach to that in the static case. Thus we will replace the virtual “displacement” column matrix $\hat{l}_s$ by those given by the following relation:

$$ \hat{l}_s \otimes H \otimes \hat{l}_s $$

$H$ is a $4n_s \times 4n_s$ quasi-diagonal matrix, $n_s$ being the total number of nodes:

$$ \begin{array}{c} \begin{bmatrix} H_1 \\ H \\ H_1 \\ H_1 \end{bmatrix} \otimes \begin{bmatrix} 0 & 0 & 0 & \frac{l_0}{l_3} \\ 1 & 0 & 0 & \frac{l_1}{l_3} \\ 0 & 1 & 0 & \frac{l_2}{l_3} \\ 0 & 0 & 0 & 0 \end{bmatrix} \otimes H_1 \otimes \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \end{bmatrix} \otimes H_1 \otimes \begin{bmatrix} 0 \\ 1 \\ 0 \\ 0 \end{bmatrix} \otimes H_1 \otimes \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \end{array} $$

where we could consider that $l_3$ is the absolute maximum of the four Euler-Rodrigues parameters for each node and found out at each iteration cycle. Thus all matrices, $M, M_v, C, K$ and $F$ become:
$M$ $\square$ $H^T M$; $M_v$ $\square$ $H^T M_v$; $C$ $\square$ $H^T C$; $K$ $\square$ $H^T K$; $F$ $\square$ $H^T F$

and will have one zero row at each node corresponding to absolute maximum of the four Euler-Rodrigues parameters of every node. Only in the matrix K instead of this zero row we add the condition [15] for each node. Finally the non-linear system of equations describing the movement of the assembly of SFET becomes:

$$ M \hat{l}_s^T \square C \square M_v \square \hat{l}_s^T \square K \hat{l}_s \square F \quad [28] $$

where C is the damping matrix which will not be discussed in this work. The column matrix F is found out in the same manner as in the case of static analysis.

The five matrices $M$, $Mv$, $C$, $K$ and $F$ depend on the configuration of the macro-element. In other words the mentioned matrices depend on the time $t$. Moreover, the square matrices $M$, $Mv$, $C$ and $K$ of the macro-element are full and non-symmetrical, but this is not necessarily a disadvantage because the number of the finite elements for each macro-element is rather small.

To solve the non-linear system [28] the well-known Newmark method was adopted. At every iteration cycle the following conditions are imposed for each node:

$$ \hat{l}_j^T \square \hat{l}_j \square 0; \quad \hat{l}_j^T \square \hat{l}_j \square \hat{l}_j^T \square 0 $$

If the studied macro-element is not clamped, that is the beam can perform translation movements, up to three extra degrees of freedom and their corresponding equations are to be added to the non-linear system [28]. These supplementary degrees of freedom represent the translation displacement of the point 0. Only the mass matrices $M$ and $Mv$ are modified. The matrix $F$ is modified only if external concentrated forces act on the beam, in the node 0.

## 5. Examples

Firstly it is interesting to point out that SFET3 is exact for small displacement of straight beams loaded statically with concentrated forces. In 2D case, for small displacements, we have $l_3$ $\square$ $\sin(\frac{s}{2})$ $\square$ $\frac{s}{2}$, where $\theta$ is the slope of the beam. According to the Equation [4] the slope $\theta$ is a polynomial of second degree in the curved coordinate $S_1$ that is exact following the Euler-Bernoulli beam theory. For instance a cantilever beam loaded with the concentrated force $F$ on the free end has:

$$ \square(s) \square (1 \square \frac{s}{h}) \square_i \square \frac{s}{h} \square_{i+1} $$

the free end displacement and slope, respectively. As we expected, only one SFET3 leads to the exact value for the slope. The accuracy for the displacement depends on
the numerical method used to compute the displacements starting from nodal unknowns, the slopes, not on SFET finite element, and it could be as small as we want.

To illustrate the performances of the special finite element type (SFET) we present several examples, in static and dynamic regime, in 2D and 3D. The first one is the 2D post-buckling static behaviour of the Euler beam problem. This example was chosen to check the proposed solution to the exact solution given by Timoshenko, [TIM 61]. The number of finite elements was chosen in order to get the displacements for the free end of the beam with a relative error of $10^{-4}$. The cantilever beam is represented on Figure 3. Figure 4 shows the various positions of the reference lines when the ratio of F over the Euler load $F_{cr}$ is increasing from 1.015 to 9.116. It can be observed that the displacement along the x-axis is almost twice the length of the beam and that the angle of rotation of the final cross section is about $\pi$. This example indicates that with space decomposition into 10 finite elements SFET3 (21 nodes) we obtain the exact solution for the free end rotation with a relative precision of $10^{-4}$. Table 1 gives some numerical results for displacements.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Post-buckling static behaviour of a cantilever</figcaption></figure>

The computing of the buckling force is an eigenvalue problem that could be solved by SFET, too. For instance the exact Euler buckling force and the corresponding values computed using 10 or 20 SFET3 are respectively:

<table><thead><tr><th></th><th></th><th colspan="7">□(s) □ (1 □ <sup>S</sup>/<sub>h</sub>) □<sub>i</sub> □ <sup>S</sup>/<sub>h</sub> □<sub>i</sub>□<sub>l</sub></th></tr><tr><th></th><th></th><th>F/F<sub>cr</sub></th><th>1.015</th><th>1.152</th><th>1.518</th><th>1.884</th><th>2.541</th><th>4.029</th><th>9.116</th></tr><tr><th></th><th></th><th>α</th><th>20</th><th>60</th><th>100</th><th>120</th><th>140</th><th>160</th><th>176</th></tr><tr><th>Timoshenk</th><th>u/L</th><th><strong>0.0300</strong></th><th><strong>0.259</strong></th><th><strong>0.651</strong></th><th><strong>0.877</strong></th><th><strong>1.107</strong></th><th><strong>1.1340</strong></th><th><strong>1.577</strong></th></tr><tr><th>o</th><th>v/L</th><th>0.220</th><th>0.593</th><th>0.792</th><th>0.803</th><th>0.750</th><th>0.625</th><th>0.421</th></tr><tr><th>10 SFET3</th><th>u/L</th><th><strong>0.0292</strong></th><th><strong>0.2585</strong></th><th><strong>0.6504</strong></th><th><strong>0.8763</strong></th><th><strong>1.1066</strong></th><th><strong>1.3403</strong></th><th><strong>1.5774</strong></th></tr><tr><th>21 nodes</th><th>v/L</th><th>0.2157</th><th>0.5928</th><th>0.7914</th><th>0.8032</th><th>0.7504</th><th>0.6244</th><th>0.4209</th></tr></thead></table>
<table><caption>**Table 1.** Comparison between analytical and numerical SFET results</caption><tr><td>20 SFET3</td><td>u/L</td><td><b>0.0302</b></td><td><b>0.2592</b></td><td><b>0.6506</b></td><td><b>0.8764</b></td><td><b>1.1066</b></td><td><b>1.3402</b></td><td><b>1.5772</b></td></tr><tr><td>41 nodes</td><td>v/L</td><td>0.2190</td><td>0.5934</td><td>0.7915</td><td>0.8032</td><td>0.7505</td><td>0.6246</td><td>0.4213</td></tr></table>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Several post-buckling equilibrium configuration</figcaption></figure>

The next example is shown in Figure 5, the cantilever loaded with a force and a moment on the free end. All the data are given in figure. The results were checked using the COSMOS 2.0 finite element code. The tolerance was 1e-5 for both, SFET3 program and COSMOS program. COSMOS could reach only 40% of the final load in 231 load steps, with several iteration cycles being done for each load step, Figure 6a. AUTOSTEP command was activated. For this load the results got by SFET3 and COSMOS are practically the same. On the contrary, SFET3 needed only one load step for 40% of the final load and 3 load steps and 25 iteration cycles for the final load. The final configuration of the beam is given in the Figure 6b. In order to verify the results, the sectional efforts (moments) are computed following two ways. Using the static equation applied on the deformed beam we can write:

$$ \Box(s) \Box(1 - \frac{s}{h}) \Box_i \Box \frac{s}{h} \Box_{i-1} $$

where $R_F$ is the position vector of the force applied on the free end of the beam and $F_1 \Box 0 \ 0 \ F_z \Box$; and $C_1 \Box 0 \ 0 \ C_z \Box$. Using the physical relation we obtain:

$$ \Box(s) \Box(1 - \frac{s}{h}) \Box_i \Box \frac{s}{h} \Box_{i-1} $$
where the curvature vector $\chi(s)$ is computed by means of the nodal Euler-Rodrigues parameters, $l_s$, found out after solving the FE problem. The average relative error is computed with the formulas:

$$ \Box(s) = \left(1 - \left(\frac{s}{h}\right)^{\frac{\epsilon_1}{\epsilon_2}}\right) \Box_i \left(\frac{s}{h}\right)^{\frac{\epsilon_2}{\epsilon_1}} \quad [29] $$

L=100
GJ=2000/2.6= 769.23
EI2=EIy=1500
EI3=EIz=1000
A=5(cross section area)
ρ=10⁻⁶ (mass density)
0.1.20 SFET3
Fz=0.25
Cz=120

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Cantilever beam loaded with a force and a moment</figcaption></figure>

The moments were computed in local axis system $R_B(G, B_i)$. For the above example and for the final load these errors are $\epsilon_1=0.06\%$ and $\epsilon_2=0.5\%$. The error $\epsilon_2$ might seem too big, but we have to consider that the beam is highly distorted, the torsion angle included.

The third statical example is shown in the Figure 7, a helicoidal spring, an initial 3D curved beam subjected to a concentrated force. All the data, the initial configuration and the final one are shown in the Figure 7. Remarkable is that to reach the final value of the load only one load step and 9 iterations were needed, for an absolute tolerance of $10^{-5}$ for the Euler-Rodrigues parameters. In this case the errors computed accordingly the Relations [29] are $\epsilon_1=0.6\%$ and $\epsilon_2=0.8\%$ respectively, that is an extremely good accuracy for stresses (not for displacements).
<figure><img src="image_9.png" /><figcaption>a) 40% of final load</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>b) Final load F<sub>z</sub>=0.25 and M<sub>z</sub>=120</figcaption></figure>

<caption>**Figure 6.** *Deformed configuration of the cantilever*</caption>

R=10 (spring radius)
Thread=10
n=2 (number of threads)
F=10 (vertical direction)
GJ=2000/2.6=769.23
EI2=1500
EI3=1000
40 SFET3

<figure><img src="image_11.png" /><figcaption>Undeformed configuration</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Deformed configuration</figcaption></figure>

<caption>**Figure 7.** *3D curved beam*</caption>

Figures 8 and 9 present the accuracy of the method. In Figure 8, the strong non-linear dependence of force versus displacement at the point of application of the force is represented and shows that very good results are got with only 10 SFET3. Figure 9 shows how fast the accuracy increases with the number of iterations for 40 SFET3 model: at each iteration step the error diminishes over ten times.
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Force-displacement dependence</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Error versus number of iterations plot</figcaption></figure>

The first dynamic example is a very simple one, the eigenfrequencies and eigenmodes of the straight cantilever represented in Figure 5. This problem was solved in the small displacements hypothesis, the stiffness and mass matrices were computed as shown in this work. The problem was solved as 3D one, the numerical results using SFET3 being compared to the analytical ones. The results are given in the Table 2 and a extremely good concordance might be noticed. Figure 10 contains the first 6 eigenmodes in the xy plane.

<table><caption>**Table 2.** Comparison between eigenfrequencies computed numerically (SFET3) and analytically</caption><thead><tr><th>Mode</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th><th>#7</th><th>#8</th><th>#9</th><th>#10</th></tr></thead><tbody><tr><td>Anal.</td><td>4.972</td><td>6.090</td><td>31.16</td><td>38.17</td><td>87.24</td><td>106.8</td><td>171.0</td><td>209.4</td><td>282.6</td><td>346.1</td></tr><tr><td>SFET3</td><td>4.971</td><td>6.089</td><td>31.15</td><td>38.15</td><td>87.24</td><td>106.9</td><td>171.0</td><td>209.5</td><td>283.1</td><td>346.6</td></tr></tbody></table>

Figure 11 shows one example extracted from the paper [SIM 86], the “flying spaghetti”. All the data are presented in the Figure 11. Firstly the problem was solved exactly like in the work [SIM 86]. The problem is a 2D one as the load consists only in the force F and the moment $T_z$. The comparison between the results from the work [SIM 86] and those obtained using SFET2 is exposed in Figure 12. A very good concordance might be noticed. To get a 3D problem the moment $T_y$ was added to the load, Figure 11, and the results are shown in Figures 13 and 14.
<figure><img src="image_15.png" /><figcaption>**Figure 10.** First eigenmodes of the cantilever beam computed using SFET3</figcaption></figure>

The last example was studied in [SHA 98] and [TAK 99], the free falling of a flexible pendulum. The pendulum is assumed to fall under the effect of gravity, Figure 15. The beam has a length of 1.2 m, a cross sectional area of 0.0018 m², a mass density of 5540 Kg/m³, a modulus of elasticity of 0.700x10⁶ N/m² and the moment of inertia of 1.215x10⁻⁶ m⁴.

<figure><img src="image_16.png" /><figcaption>**Figure 11.** The “flying spaghetti”</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 12.** *The first 7.5 s of the movement of the “flying spaghetti” in 2D*</figcaption></figure>

The pendulum was studied firstly as a 2D problem and divided into 20 SFET2 finite elements (21 nodes) with only one degree of freedom per node, the slope. Then the pendulum was considered a 3D problem and was divided into 10 SFET3 finite elements (21 nodes) with four degree of freedom per node. The two results are practically identical and extremely close to those obtained in the above-mentioned papers. Figure15a shows 12 successive positions of the pendulum during 1.1 seconds with a time step of 0.1s.

<figure><img src="image_18.png" /><figcaption>**Figure 13.** *3D “flying spaghetti”: the first 15 s of the movement*</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 14.** 3D “flying spaghetti”: the projection of the successive positions</figcaption></figure>

In order to verify the accuracy of SFET two methods were applied. First the total energy of the pendulum was computed at each time step. It comprises three components: the elastic energy, the potential energy and the kinetic energy. Their sum, the total energy, has to be constant during the motion. Figure 15b shows the variation of all these energies with the time. The error was estimated using the relation:

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

where E is the total energy and $E_k$ is the kinetic energy. Its value depends on the time step and for this example we got 2.65% for 220 time steps, 1.05% for 550 time steps and 0.688% when 1100 time steps were used. The second method for computing the error uses the relations [29], but this time to compute the bending moment the gravity distributed force as well as the inertia forces were considered.

<figure><img src="image_21.png" /><figcaption>**Figure 15.** *Flexible pendulum*</figcaption></figure>

The errors $\epsilon_1$ and $\epsilon_2$ were evaluated at each time step. The maximum error do not exceed 1...2 % which is a very good accuracy taking into account that the error was computed for the bending moment.
# 6. Conclusions

A special finite element (SFET), which has four degrees of freedom per node, the nodal values of the Euler-Rodrigues parameters, was presented in this paper. The main advantage of SFET is that it allows without many complications to solve the exact movement equations of the mechanical systems composed by elastic beams. The proposed finite element may have several nodes, but SFET3 proved to be extremely accurate and simple enough from the programming point of view.

The method is very accurate and very rapidly convergent. This might be explained essentially by the following factors: (I) the exact equations are solved, (II) the unknowns are not the linear displacements, but their derivatives, and also (III) the stiffness and mass matrices are polynomials in Euler-Rodrigues parameters, thus avoiding to use trigonometric functions.

The accuracy for the displacement depends on the numerical method used to compute the displacements starting from nodal unknowns, the Euler-Rodrigues parameters. To compute the displacements the Relation [12] of the first part of the paper has to be used. This integration might be done introducing the shape functions [2] or [4] in the above mentioned relation and then using as many integration points as the accuracy requires. The accuracy for displacement was assessed directly by computing efforts considering both internal and external equilibrium of the beam (See the Relations [29]), and extremely good accuracy was again obtained.

The finite element formulation is written in the global coordinate system, the same for each macro-element in the mechanical system. This approach has several advantages one of the most important is that a simple expression for inertia forces is obtained. In addition due to the good accuracy and convergence of the method, the computer codes elaborated by the authors are very fast, a computer time of several minutes are enough for all 3D examples presented in the work. All the computer codes used in this paper and applying SFET were written by the authors in MATLAB 5.3.

Several examples are presented in this paper and they confirm the accuracy of the SFET. It is important to notice that the computed errors concern the efforts (bending and/or torsion moments) and not the nodal unknowns. The bending and torsion moments are expressed in derivatives of Euler-Rodrigues parameters. The errors for derivatives are much bigger than the errors for the functions, as it is known!

In the static domain the finite element type elaborated by the authors is a non-incremental one, usually only one load step is sufficient to reach the final value, but in any case the accuracy of the results does not depend on the number of load steps. This happens because the exact equations (obviously in the limits of the beam theory) were solved and they are written **the actual deformed configuration** for of the beam system. This is not possible for the classic beam finite element.
# 7. Bibliographie

[BAR 96] BARRACO A., "Dynamique des systèmes mécaniques complexes", *ENSAM, Centre d'enseignement et de recherche de Paris*, 1996.

[BAR 00] BARRACO A., "Mécanique des structures", *ENSAM, Centre d'enseignement et de recherche de Paris*, 2000.

[DVO 88] DVORKIN E. N., ONATE E., OLIVIER J., "On A Non-linear Formulation for Curved Timoshenko Beam Elements Considering Large Displacement/Rotation Increments", *Int. Journal for Numerical Methods in Engineering*, Vol. 26, pp. 1597-1613, 1988.

[MUN 96] MUNTEANU M. Gh., DE BONA F., ZELENIKA S., "An accurate non-linear analysis of very large displacements of beam systems", *Proceedings of the XXV AIAS National Conference, International Conference on Material Engineering, Gallipoli - Lecce*, pp. 59-66, 1996.

[MUN 00] MUNTEANU M. Gh., BARRACO A., CURTU I., DE BONA F., "Experimental Validation of a Special Finite Element Type", *Proceedings of the Danubia-Adria Symposium, Prague*, pp. 229-232, 2000.

[REI 73] REISSNER E., "On one dimensional large displacement finite strain beam theory", *Studies Appl. Math.* 52, pp. 87-95, 1973.

[SAJ 91] SAJE M., "Finite Element Formulation of Finite Planar Deformation of Curved Elastic Beams", *Computers & Structures*, Vol. 39, No 3/4, pp. 327-337, 1991.

[SHA 98] SHABANA A. A., HUSSIEN H. A., ESCALONA J. L., "Application of the Absolute Nodal Coordinate Formulation to Large Rotation and Large Deformation Problems", ASME, Journal of Mechanical Design, Vol. 120, pp. 188-195, 1998.

[SIM 86] SIMO J.C., VU-QUOC L., "On the dynamics of Flexible Beams Under Large Overall Motions – The Plane Case", *Journal of Applied Mechanics*, Vol 53, pp. 855-863, 1986.

[SIM 88] SIMO J.C., VU-QUOC L., "On the dynamics in space of rods undergoing large motions – a geometric exact approach", *Comput. Meth. Appl. Mech. Eng.* 66, pp. 125-161, 1988.

[TAK 99] TAKAHASHI Y., SHIMIZU N., "Study on Elastic Forces of the Absolute Nodal Coordinate Formulation for Deformable Beams", *Proceedings of the 1999 ASME Design Engineering Technical Conferences September 12-15, Las Vegas, Nevada*, 9 pp., 1999.

[TIM 61] TIMOSHENKO S. P., "Theory of Elastic Stability", *McGraw-Hill Book Company, New York*, 1961.