The beam finite element is one the main elements proposed by structural finite
element analysis software. This element deserves a full article to explain what
it is about. I will focus on the
Timoshenko beam element
because it is the most
common beam element implemented in FEA software. First, I will expose the
**assumptions underlying this element, as well as how to derive the beam element
stiffness matrix**. Finally, I will conclude this article with an example of
a cantilever beam model simulated using the
SesamX finite element analysis software.

## Quick remark on notations

Before starting, let’s define some notations that are used through this article:

vectors are denoted with an underline $ \underline{u} $ ,

matrices are represented with brackets $ \lbrack \ \rbrack $ ,

Einstein summation convention is used on repeated indices,

partial derivatives are denoted with the comma notation $ u_{,x} $ ,

infinitesimal strain and stress tensors are represented in column matrix notation $ \lbrack \varepsilon \rbrack = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \varepsilon_{12} \\ \varepsilon_{23} \\ \varepsilon_{13} \end{bmatrix} $ and $ \lbrack \sigma \rbrack = \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{23} \\ \sigma_{13} \end{bmatrix} $ , where $ \lbrack \varepsilon \rbrack $ represents the

**engineering strains**.

## Introducing the beam element

The beam element
is relevant to use when we aim at analyzing a **slender structure**
undergoing forces and moments in any direction. For instance, it makes it the
perfect element to analyze the support of a slab or a plate stiffener.

Slab support |
Stiffened plate |

A beam can be more simplistically represented as follows.

We are interested here in a beam whose
neutral axis
is the same as the
shear axis.
In other words, **bending is supposed uncoupled from shearing**. From beam theory
(or as it will appear during the following derivation), the relevant geometric
parameters of the beam are:

The

**length**of the beam: $ L $The

**area**of the beam cross-section: $ A $The

**area moments of inertia**of the beam cross-section:$$ I_y = \iint_A z^2 dydz $$

$$ I_z = \iint_A y^2 dydz $$

The

**product moment of inertia**of the beam cross-section:$$ I_{yz} = \iint_A yz dydz $$

The

**polar moment of inertia**:$$ I_{x} = \iint_A (y^2 + z^2) dydz $$

## Timoshenko beam element assumptions

We can view a beam element as a simplification of a more complex 3D structure. When designing such an element, the aim is to decrease the computational cost, while making relevant assumptions on how the element behaves. The main assumptions behind the Timoshenko beam element are given here with some explanations.

### 1. Stress assumption

As the beam must be able to sustain transverse loading, the stress state must allow for shear stress, as well as axial stress in order to sustain bending and axial loading. It is then sufficient to enforce:

$$ \sigma_{22} = \sigma_{33} = \sigma_{23} = 0 \tag{1} $$

Remark that unlike the truss element, the beam element makes no assumption about the shear stresses $ \sigma_{12} $ and $ \sigma_{13} $ .

### 2. Kinematic assumption

The main assumption associated to the beam element, is that **beam cross section
remains straight during deformation**. However, the
Timoshenko beam theory
does not enforce that a surface normal to the beam neutral axis remains normal
during the deformation (as in classical
Bernoulli beam theory).

Timoshenko beam is preferred because it makes looser assumptions on the beam
kinematics. In fact, Bernoulli beam is considered accurate for cross-section
typical dimension less than ^{1}⁄_{15} of the beam length. Whereas **Timoshenko beam
is considered accurate for cross-section typical dimension less than ^{1}⁄_{8} of
the beam length**.

The displacement inside the beam is defined from the displacements on the neutral axis $ u_1^0(x), u_2^0(x), u_3^0(x) $ , as well as the cross-section rotations $ \theta_1(x), \theta_2(x), \theta_3(x) $ :

$$ \underline{u} = \begin{cases} u_1(x, y, z) = u_1^0(x) + z\theta_2(x) - y\theta_3(x) \\ u_2(x, y, z) = u_2^0(x) - z\theta_1(x) \\ u_3(x, y, z) = u_3^0(x) + y\theta_1(x) \tag{2} \end{cases} $$

Once we know the displacements and rotations on the beam axis, we can compute
the displacement over the whole beam. Therefore, the beam element is a
1-dimensional element. The element presented here is the linear beam element.
Hence, this element consist of **2 nodes connected together through a segment**.

We can see from the assumed displacement given in (2), that it **does not take into
account the
warping effect**
of beams with non circular cross-section. To
approximately account for this effect, we can simply replace the beam polar
inertia with the
**beam torsion constant**
in our equations, leaving our assumptions
unchanged. The link provided previously explains how to compute the torsion
constant, given the cross-section dimensions.

### 3. Inertia assumption

We make the assumption that the beam cross-section principal axes are along the $ y $ and $ z $ axis. Which translate, to

$$ I_{yz} = 0 \tag{3} $$

## Beam element derivation

As mentioned previously, we can represent the beam element as shown in the following figure, with its local axis system.

The element local axis system is defined by the axis $ (x, y, z) $ whose basis vectors are respectively denoted as $ \underline{e_{1}^{elem}}, \underline{e_{2}^{elem}}, \underline{e_{3}^{elem}} $ . (Here $ x $ is directed from node 1 to node 2 however in SesamX the user is free to chose this orientation).

The translational and rotational degrees of freedom are required on each node of the element. Hence, the displacement vector on node $ I $ will be denoted $ \underline{u^I} = {u_j}^I \underline{e_{j}^{elem}} $ , and the rotation vector will be denoted $ \underline{\theta^I} = {\theta_j}^I \underline{e_{j}^{elem}}$ where $ I $ takes the values 1 or 2.

Applying the classical finite element method we can derive the beam element stiffness matrix.

### Strain tensor

The displacements and rotations are known at the nodes. To get the displacements
and rotation inside the element (but still on the neutral axis) **we
interpolate linearly** the quantities at the nodes:

$$ \begin{cases} \begin{aligned} {u^0_j}(x) &= N^I(x) {u_j}^I \\ {\theta_j}(x_1) &= N^I(x) {\theta_j}^I \end{aligned} \end{cases} $$

Where the $ N^I $ represent the shape functions of the element:

$$ \begin{cases} N^1(x) = 1 - \cfrac{x}{L} \\ N^2(x) = \cfrac{x}{L} \end{cases} $$

Using assumption (2) we can then compute the strains:

$$ \lbrack \varepsilon \rbrack = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \varepsilon_{12} \\ \varepsilon_{23} \\ \varepsilon_{13} \end{bmatrix} = \begin{bmatrix} u^0_{1,x} + z\theta_{2,x} - y\theta_{3,x} \\ 0 \\ 0 \\ u^0_{2,x} - z\theta_{1,x} - \theta_{3} \\ 0 \\ u^0_{3,x} + y\theta_{1,x} + \theta_{2} \end{bmatrix} $$

Since the actual shear strain in the beam is not constant over the
cross-section, **Timoshenko beam theory usually introduces
shear correction factors**
(depending on the cross-section shape)
$ k_2 $
and $ k_3 $
on $ \varepsilon_{12} $
and $ \varepsilon_{13} $
.

To make things easier, we define the following quantities:

the curvatures $ \kappa_1 = \theta_{2,x} $ and $ \kappa_2 = \theta_{3,x} $

the shear strains $ \gamma_1 = u_{2,x} - \theta_3 $ and $ \gamma_2 = u_{3,x} + \theta_2 $

the axial strain $ \varepsilon_{ax} = u_{1,x} $

the twist $ \psi = \theta_{1,x} $

And rewrite the strains as:

$$ \lbrack \varepsilon \rbrack = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \varepsilon_{12} \\ \varepsilon_{23} \\ \varepsilon_{13} \end{bmatrix} = \begin{bmatrix} \varepsilon_{ax} + z\kappa_1 - y\kappa_2 \\ 0 \\ 0 \\ k_2\gamma_1 - z\psi \\ 0 \\ k_3\gamma_2 + y\psi \end{bmatrix} \tag{4} $$

Where the curvatures, shear strains, axial strain and twist are provided with the help of the shape functions:

$$ \begin{aligned} \kappa_1 &= N^I_{,x} {\theta_2}^I \\[5pt] \kappa_2 &= N^I_{,x} {\theta_3}^I \\[5pt] \gamma_1 &= N^I_{,x} {u_2}^I - N^I {\theta_3}^I \\[5pt] \gamma_2 &= N^I_{,x} {u_3}^I + N^I {\theta_2}^I \\[5pt] \varepsilon_{ax} &= N^I_{,x} {u_1}^I \\[5pt] \psi &= N^I_{,x} {\theta_1}^I \end {aligned} \tag{5} $$

### Stress tensor

To relate the stresses to the strains we need to apply Hooke’s law for linear elastic material.

$$ \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \varepsilon_{12} \\ \varepsilon_{23} \\ \varepsilon_{13} \end{bmatrix} = \frac{1}{E} \begin{bmatrix} 1 && -\nu && -\nu && 0 && 0 && 0 \\ -\nu && 1 && -\nu && 0 && 0 && 0 \\ -\nu && -\nu && 1 && 0 && 0 && 0 \\ 0 && 0 && 0 && 2+2\nu && 0 && 0 \\ 0 && 0 && 0 && 0 && 2+2\nu && 0 \\ 0 && 0 && 0 && 0 && 0 && 2+2\nu \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{23} \\ \sigma_{13} \end{bmatrix} \\ $$

Using assumption (1) on the right hand side, as well as the result we got in (4) on the left hand side, leads to a peculiar relation:

$$ \underbrace{ \begin{bmatrix} \varepsilon_{ax} + z\kappa_1 - y\kappa_2 \\ 0 \\ 0 \\ k_2\gamma_1 - z\psi \\ 0 \\ k_3\gamma_2 + y\psi \end{bmatrix} }_{\text{kinematic assumption}} = \underbrace{ \frac{1}{E} \begin{bmatrix} \sigma_{11} \\ -\nu\sigma_{11} \\ -\nu\sigma_{11} \\ (2+2\nu)\sigma_{12} \\ 0 \\ (2+2\nu)\sigma_{13} \end{bmatrix} }_{\text{stress assumption}} $$

As for the the truss element,
this relation cannot hold. However using a similar
explanation, we can find a way out.
When making the kinematic assumption we were interested in the macroscopic
behavior of the beam. Whereas the stress assumption relates more to a microscopic
behavior. It turns out that even if $ \varepsilon_{22} $
and $ \varepsilon_{33} $
are not 0 (**microscopic scale**) their
effect on the displacement (**macroscopic scale**) is negligible compare to what
comes from the $ \varepsilon_{11}, \varepsilon_{12}, \varepsilon_{13} $
terms. Of course, this explanation becomes questionable as the slenderness of
the beam degrades.

We can then simplify this relation and write:

$$ \begin{bmatrix} \varepsilon_{ax} + z\kappa_1 - y\kappa_2 \\ 0 \\ 0 \\ k_2\gamma_1 - z\psi \\ 0 \\ k_3\gamma_2 + y\psi \end{bmatrix} = \frac{1}{E} \begin{bmatrix} \sigma_{11} \\ 0 \\ 0 \\ (2+2\nu)\sigma_{12} \\ 0 \\ (2+2\nu)\sigma_{13} \end{bmatrix} $$

And we can invert it to get the stresses from the strains:

$$ \begin{bmatrix} \sigma_{11} \\ \sigma_{12} \\ \sigma_{13} \end{bmatrix} = \begin{bmatrix} E && 0 && 0 \\ 0 && \cfrac{E}{2(1+\nu)} && 0 \\ 0 && 0 && \cfrac{E}{2(1+\nu)} \end{bmatrix} \begin{bmatrix} \varepsilon_{ax} + z\kappa_1 - y\kappa_2 \\ k_2\gamma_1 - z\psi \\ k_3\gamma_2 + y\psi \end{bmatrix} $$

## Beam element stiffness matrix

The element stiffness matrix is obtained through the expression of the virtual work. Denoting the virtual strains as $ \lbrack \overline{\varepsilon} \rbrack $ we have for the virtual work over the element volume:

$$ \overline{W} = \iiint_V \lbrack \overline{\varepsilon} \rbrack^T \lbrack \sigma \rbrack dV $$

After expanding we get:

$$ \begin{alignedat}{3} \overline{W} &= &E\iiint_V &\lbrack \overline{\varepsilon_{ax}}\varepsilon_{ax} + z^2 \overline{\kappa_1}\kappa_1 + y^2 \overline{\kappa_2}\kappa_2 + \overline{\varepsilon_{ax}}(z\kappa_1 - y\kappa_2) + \\ &&& (z\overline{\kappa_1} - y\overline{\kappa_2})\varepsilon_{ax} - zy \overline{\kappa_1} \kappa_2 - yz \overline{\kappa_2} \kappa_1 \rbrack dV \\ &+&\cfrac{E}{2(1+\nu)} \iiint_V &\lbrack k_2\overline{\gamma_1}\gamma_1 + z^2 \overline{\psi}\psi - z \overline{\gamma_1}\psi - k_2z \overline{\psi}\gamma_1 \rbrack dV \\ &+&\cfrac{E}{2(1+\nu)} \iiint_V &\lbrack k_3\overline{\gamma_2}\gamma_2 + y^2 \overline{\psi}\psi + y \overline{\gamma_2}\psi + k_3y \overline{\psi}\gamma_2 \rbrack dV \end{alignedat} $$

Using the fact that the beam is meshed on its neutral axis ($ \iiint_V zdV = \iiint_V ydV =0 $ ) and the inertia assumption (3), we can simplify this expression:

$$ \begin{alignedat}{4} \overline{W} &=& EA\int_L &\overline{\varepsilon_{ax}}\varepsilon_{ax} dx &\to \textbf{axial contribution} \\ &+& \cfrac{EI_x}{2(1+\nu)} \int_L &\overline{\psi}\psi dx &\to \textbf{twist contribution} \\ &+& E\int_L & \lbrack I_y \overline{\kappa_1}\kappa_1 dx + I_z \overline{\kappa_2}\kappa_2 \rbrack dx &\to \textbf{bending contribution} \\ &+& \cfrac{EA}{2(1+\nu)} \int_L & \lbrack k_2 \overline{\gamma_1}\gamma_1 + k_3 \overline{\gamma_2}\gamma_2 \rbrack dx &\to \textbf{shear contribution} \end{alignedat} $$

Before we compute the integral, it is interesting to analyse the polynomial order of each integrand:

the axial, twist and bending integrands are of

**order 0**, because they introduce only the $ N^I_{,x} N^J_{,x} $ products,the shear integrand is of

**order 2**, because it introduces the $ N^I N^J $ products (as well as the $ N^I_{,x} N^J_{,x} $ products).

For reasons not explained here, we can show that if the shear terms are integrated
exactly, it will lead to
element shear locking:
the element will display a very stiff behavior as the length
decreases. To avoid this phenomenon we degrade the shear terms integration
(called **reduced integration**).
The 1 point Gaussian integration rule
integrates exactly polynomials up
to the order 1. Thus using this rule, we will integrate exactly the
axial, twist and bending integrands and we will reduce the integration on the
shear terms.

Eventually, using (5), we can relate each of the terms obtained to the nodes displacements and rotations.

## Example of a simple structure using beam elements

As an example, the following cantilever beam model is set up and simulated with the SesamX finite element software.

The beam is $ 1 m $ long. It is made of 101 nodes and 100 elements. The beam orientation is such that it is exactly aligned with the global coordinate system. The beam properties come from a square section of $ 10 cm $ side:

- area of $ 10000 mm^2 $
- $ J = 1.406e7 mm^4 $
- $ I22 = I33 = 8.33e6 mm^4 $
- $ K2 = K3 = 0.44 $

A linear elastic material is applied with $ E = 200 GPa $ and $ \nu = 0.33 $ . The far left node is clamped while, arbitrary force and moment vectors are applied at the other end.

This case is then solved with a **linear static resolution**.

And the following figure gives an overview of the displacement of the model, as well as the node numbers.

To reproduce this simulation on your own, you need:

- at least the free version of the SesamX software installed,
- the beam model mesh (in Abaqus format),
- and the SesamX input cards (defining the element properties, the loading and so on).

If you are looking to extend this example, you can get help from the SesamX documentation.

## Conclusion

I hope this article has provided you a better understanding of the beam finite element. If you are looking for more information, do not hesitate to send me and email.

#### Did you like this content?

Register to our newsletter and get notified of new articles