Home > Theory > Lines Modelling > Beam Element

Beam element

The beam element has been developed at the Ecole Centrale de Paris. The mathematical model uses two different co-ordinate systems :

The global co-ordinate system $(O, \vec {e_1}, \vec {e_2}, \vec {e_3})$ common to all elements.

In the unstressed configuration $\Omega_{\xi}$ of each beam, the torsion axis is along the Z-axis of this global frame. In this frame the physical values are denoted with the exponent $\xi$.

The local co-ordinate system $(P, \vec{d_{x1}}, \vec{d_{x2}}, \vec{d_{x3}})$ at each node P. In the deformed configuration $\Omega_x$, this system of axes is connected to the beam section at this node. When asking for local outputs, direction $i$ coreesponds to vector $d_{xi}$.

Kinetics assumptions

A point of the structure in the global reference may be located by the vector $\overrightarrow {OP} = \vec {x_o} + \vec{x}^*$. The first vector gives its position on the beam mean line and the second one its position on the beam section. As the beam section is assumed to be rigid, $x^*$ is a simple rotation of the initial position on the beam $\vec{\xi}^*$ : $\vec{x}^* = R_x\vec{\xi}^*$.

With these kinetics assumptions, the independent unknowns are the position and the rotation vectors $(\vec{x_o}, \vec{\theta_x})$, both functions of the arc length $s$.

The Mindlin assumption is made i.e. the normal vector of a beam section is not necessary parallel to the beam mean line which leads to shear strains.

\[ \Omega_x:~ \begin{aligned} \overrightarrow {OP}&= \vec {x_o}(s) \\ \overrightarrow {OM}&= \vec {x_o}(s) + \vec x ^* = \vec {x_o}(s) + R(\vec{\theta_y}) \vec {\xi ^*} \end{aligned} \quad ; \quad \Omega_\xi:~ \begin{aligned} \overrightarrow {OP}&= s\vec{e_3} \\ \overrightarrow {OM}&= s\vec{e_3} + \vec {\xi}^* \end{aligned} \]

3D large rotations

A special thought is needed to deal with large rotations in space since they don't behave like true vectors and they are not actually commutative.

Rotations are defined by a rotation vector $\vec{\theta_x}$. Thanks to the Rodrigues formula, this vector gives the rotation matrix $R_x$ which defines the local reference $(d_{x1},d_{x2},d_{x3})$ at this node : $\vec{d_{xi}} =R_x(\vec {\theta_x})\vec{e_i}$. with $\vec{\theta_x} = \begin{pmatrix} \theta_1 \\ \theta_2\\ \theta_3 \end{pmatrix}$ and $\theta = \lVert \vec{\theta_x} \rVert$

\[ R_x(\vec \theta_x) = I_d + \frac{sin \theta}{\theta} \begin{pmatrix} 0 & -\theta_3 & \theta_2 \\ \theta_3 & 0 & -\theta_1 \\ -\theta_2 & \theta_1 & 0 \end{pmatrix} + \frac{1-cos\theta}{\theta^2}\begin{pmatrix} -(\theta_2^2+\theta_3^2) & \theta_1\theta_2 & \theta_1\theta_3 \\ \theta_2\theta_1 & -(\theta_1^2+\theta_3^2) & \theta_2\theta_3 \\ \theta_3\theta_1 & \theta_3\theta_2 & -(\theta_1^2+\theta_2^2) \end{pmatrix} \]

Therefore there is a non linear relation between the variation $\delta \vec{\theta_x}$ of $\vec{\theta_x}$ and the rotation increment $\delta\vec{\Omega}$ following the matrix $H(\vec{\theta_x})$:

\[ \delta \vec {\Omega} = H(\vec{\theta_x})\delta\vec{\theta_x} \]

\[ H(\vec{\theta_x}) = \frac{sin\theta}{\theta} I_d + \Big \{1 - \frac{sin\theta}{\theta} \Big \} \frac{\vec{\theta_x}\otimes \vec{\theta_x}}{\theta^2} + \frac{1-cos\theta}{\theta^2} \begin{pmatrix} 0 & -\theta_3 & \theta_2 \\ \theta_3 & 0 & -\theta_1 \\ -\theta_2 & \theta_1 & 0 \end{pmatrix} \]

Curvature

Based on the previous kinetics assumption, let derive the position vector with respect to $s$ :

\[ \frac {\partial \vec P}{\partial s} = \vec x_o^t + [D_{\theta}R_x\vec{\theta}_x^t]\vec{\xi^*} = \vec x_o^t + [D_{\theta}R_x\vec {\theta}_x^t]R_x^t \vec x^* \]

$K_x = [D_{\theta}R_x\vec {\theta}_x^t] R_x^t$ is an antisymetrical tensor and its eigenvector $\vec k_x$ may be defined :

\[ K_x \vec k_x = \vec 0 ~~\text{and} ~~K_x\vec x^{*} = \vec k_x \wedge\vec x^* \]

Hence : $\frac {\partial \vec P}{\partial s} = \vec x_o^t + \vec k_x \wedge \vec x^*$

It results from the last equation that $\vec k_x$ is the curvature vector in the global reference. This vector is composed of both torsion and flexion and it is impossible to distinguish one from another.

In the local reference, the curvature vector may be introduced in the same way as the eigenvector of $K_x^{\xi} = R_x^tK_x$. After calculations, this vector may be expressed as a non-linear function of the derivative of $\vec{\theta}_x$ with respect to $s$ :

\[\vec k_x^{\xi} = T(\vec {\theta_x}) \theta_x^t\]

\[ T(\vec {\theta_x}) = \frac {sin \theta}{\theta} I_d + \Big\{ 1- \frac{sin\theta}{\theta} \Big\} \frac{\vec{\theta_x}\otimes \vec{\theta_x}}{\theta^2} - \frac {1-cos\theta}{\theta^2}\begin{pmatrix} 0 & -\theta_3 & \theta_2 \\ \theta_3 & 0 & -\theta_1 \\ -\theta_2 & \theta_1 & 0 \end{pmatrix} \]

Here, the first two components of $\vec k_x^{\xi}$ are the curvature of the structure, the third one is the torsion.

Note that the curvature vector is a direct result of the calculations. It is a functiun of the nodal rotations which are the unknows of the system of equations.

Stresses and strains

Within the theoretical framework of large displacements and large rotations, the definition of a material law has to be done with great care, since problems of objectivity of the stress and strain variables may occur. For the elastic material behaviour, this problem has been studied and solved, by many authors (Simo). A correct formulation of the material law is obtained by replacing classical Cauchy stresses by their Lagrangian representation (see 1. and 2.). Since strains are assumed to remain small, the tensors then reduce to the representation of Eulerian in the rotated reference frame (see 2.).

\[ E = \frac 1 2 (F^TF-1) \]

$$ S = F^T s F \approx R_x^t s R_x $$ $E$, $F$ and $R_x$ designating respectively the Lagrange strain tensor, the displacement gradient, and rotation tensors. $s$ and $S$ represent the Cauchy and the Piola-Kirchhoff 2 stress tensors.

It comes :

\[ E_{\alpha3} = \frac 1 2 (\vec e_{\alpha}.R_x^t \vec x^t - \vec e_{3_o}) + \frac 1 2 (\vec e_{\alpha}. \vec k_x^{\xi}\wedge\vec {\xi}^*)_{\alpha = 1,2,3} \]

When a linear material behaviour is assumed : $S = \lambda * tr(E) + 2\mu E$, the stress vector is :

\[ \vec S = ({\lambda} ^{*}+ 2\mu )E_{33}(R_x^t \vec x_o^t + k_x^{\xi}\wedge\vec {\xi}^{*})+2\mu E_{\alpha 3}\vec e_{\alpha} \]

The internal efforts are calculated by integration on the beam section :

\[ \vec N_x^{\xi} = \int_{S_t} \vec S ds ~~~~~~and ~~~~~~~ \vec M_x^{\xi} = \int _ {S_t} \vec {\xi}^* \wedge \vec S ds \]

$\vec N_x^{\xi}. \vec e_3$ is the normal effort to the beam, $\vec N_x^{\xi}. \vec e_1$ and $\vec N_x^{\xi}. \vec e_2$ are the shear forces,

$\vec M_x^{\xi}. \vec e_3$ is the torsion momentum, $\vec M_x^{\xi}. \vec e_1$ and $\vec M_x^{\xi}. \vec e_2$ are the bending momentum.

\[ \vec M_x^{\xi} = A^{\xi}(R_x^t \vec - \vec e_3) \]

$$ M_x^{\xi} = J_I^{\xi} \vec k_x^{\xi} $$ here $A^{\xi} = \mu S[\vec e_{\alpha} \otimes \vec e_{\alpha}]_{\alpha=1,2} + ES \vec e_3 \otimes \vec e_3$

and $J^{\xi} = J_{E \alpha \beta} [\vec e_{\alpha} \otimes \vec e_{\beta}]_{\alpha,\beta=1,2}+ J_{\mu} \vec e_3 \otimes \vec e_3$.

Virtual works principle

From the kinetics assumptions, the different derivatives in time lead to express the inertia of the beam both in translation and rotation. Finally, the contribution of a beam element to the Virtual work principle writes :

Internal efforts contribution: $G_{internal}^{beam}(\vec x,\delta \vec x) = \Bigl<\vec N_x^{\xi}.\delta(R_x^t \vec x_o^t)\Bigr> + \Bigl< \vec M_x ^{\xi} . \delta \vec k_x^{\xi}\Bigr>$

Inertia contribution: $G_{internal}^{beam} = \Bigl<A_{\rho}^{\xi} \ddot {\vec x}_o , \delta \vec x_o\Bigr> + \Bigl<R_x \left(J_{\rho}^{\xi} \dot{ \vec {\omega}^{\xi}} + \vec {\omega}^{\xi} \wedge J_{\rho}^{\xi} \vec {\omega}^{\xi} \right), H \delta \vec {\theta}_x\Bigr>$

Several terms could be pointed out :

$A_{\rho}^{\xi} \ddot { \vec x_0}$ : beam inertia in translation
$J_{\rho}^{\xi} \ddot { \vec {\omega}_x^{\xi}}$ : beam inertia in rotation.
$\vec {\omega}^{\xi} \wedge J_{\rho}^{\xi} \vec {\omega}^{\xi}$: Gyroscopic effect