Timoshenko beam element explained

Learn how to derive the beam element stiffness matrix from kinematic assumptions. And discover how the SesamX element compares to Abaqus B31 element.

The beam element is one the main elements used in a structural finite element model. It makes it a must have for SesamX. In this article, I will discuss the assumptions underlying this element, as well as the derivation of the stiffness matrix implemented in SesamX. Finally, I will present the SesamX data cards useful to define a beam element, and a comparison of results with Abaqus B31 element.

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 \epsilon \rbrack = \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \epsilon_{12} \\ \epsilon_{23} \\ \epsilon_{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 \epsilon \rbrack $ represents the engineering strains.

Introduction

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

Slab support Stiffened plate
Slab support
Stiffened plate

A beam can be more simplistically represented as follows.

Beam slender structure

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 $$

Beam cross-section

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 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} $$

2. Kinematic assumption

The main assumption related to beam, is that beam cross section remains straight during deformation. However, we do not enforce that a surface normal to the beam neutral axis remains normal during the deformation (as in classical Bernoulli beam). In other words, the beam detailed in this article is a Timoshenko beam.

Timoshenko and Bernoulli beam

Timoshenko beam is chosen in SesamX because it makes looser assumptions on the beam kinematics. In fact, Bernoulli beam is considered accurate for cross-section typical dimension less than 115 of the beam length. Whereas Timoshenko beam is considered accurate for cross-section typical dimension less than 18 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} $$

Beam cross-section rotation

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. More information can be found here as well as 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.

Beam element local coordinate

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.

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 \epsilon \rbrack = \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \epsilon_{12} \\ \epsilon_{23} \\ \epsilon_{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 $ \epsilon_{12} $ and $ \epsilon_{13} $ (more information here).

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 $ \epsilon_{ax} = u_{1,x} $

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

And rewrite the strains as:

$$ \lbrack \epsilon \rbrack = \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \epsilon_{12} \\ \epsilon_{23} \\ \epsilon_{13} \end{bmatrix} = \begin{bmatrix} \epsilon_{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] \epsilon_{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} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \epsilon_{12} \\ \epsilon_{23} \\ \epsilon_{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} \epsilon_{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 $ \epsilon_{22} $ and $ \epsilon_{33} $ are not 0 (microscopic scale) their effect on the displacement (macroscopic scale) is negligible compare to what comes from the $ \epsilon_{11}, \epsilon_{12}, \epsilon_{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} \epsilon_{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} \epsilon_{ax} + z\kappa_1 - y\kappa_2 \\ k_2\gamma_1 - z\psi \\ k_3\gamma_2 + y\psi \end{bmatrix} $$

Element stiffness derivation

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

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

After expanding we get:

$$ \begin{alignedat}{3} \overline{W} &= &E\iiint_V &\lbrack \overline{\epsilon_{ax}}\epsilon_{ax} + z^2 \overline{\kappa_1}\kappa_1 + y^2 \overline{\kappa_2}\kappa_2 + \overline{\epsilon_{ax}}(z\kappa_1 - y\kappa_2) + \\ &&& (z\overline{\kappa_1} - y\overline{\kappa_2})\epsilon_{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{\epsilon_{ax}}\epsilon_{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} $

  • the shear integrand is of order 2, because it introduces the $ N^I N^J $ (as well as the $ N^I_{,x} N^J_{,x} $ )

For reasons not explained here, we can show that if the shear terms are integrated exactly, it will lead to element locking: the element will display a very stiff behavior as the thickness 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.


SesamX input cards

To define a beam element in SesamX the first step is to create a mesh. This is done with the CREATE-SUBMESH function:

CREATE-SUBMESH  
  MODEL: BEAM_STD_SIMPLE
  NODES
    - ID: 1  POINT: 0.,0.,0.
    - ID: 2  POINT: 10.,0.,0.
    - ID: 3  POINT: 20.,0.,0.
  ELEMENTS  
    TYPE: LN2
    - ID: 1  NODES: 1,2
    - ID: 1  NODES: 2,3

Here we define 3 nodes and we create 2 line elements to connect the nodes. To orientate the beam cross-section we need to define a vectors basis, with the DEFINE-BASIS function:

CREATE-BASIS  
  MODEL: BEAM_STD_SIMPLE  
  VECTORS  
    NAME: BEAM_BASIS
    U1: 1.,0.,0  
    U2: 1.,1.,0.

Here, we just define the directions for the 2 first vectors of the basis. SesamX will take care of the orthonormalization as well as computing $ \underline{u_3} $ to form an orthonormal basis. Next, we apply the beam property on the 2 elements defined.

CREATE-PROPERTIES  
  MODEL: BEAM_STD_SIMPLE
  BEAM-STANDARD  
    BEHAVIOR: LINEAR  
    NAME: PROP1
    AREA: 1.E4  
    J: 1.406E7  
    I22: 8.33E6  
    I33: 8.33E6  
    K2: 0.44  
    K3: 0.44
    BASIS: BEAM_BASIS  
    MATERIAL: STEEL
    ON-ELEMENTS-FROM: ALL_BEAM

Here we apply a BEAM-STANDARD property on the elements from ALL_BEAM providing a material name STEEL (that we defined previously, not shown here), the basis name BEAM_BASIS and the beam geometric properties (area, moments, torsion constant and shear correction factors). BEHAVIOR: LINEAR indicates that we apply a linear elastic material. This keyword may seem unnecessary at the moment, but it is a provision for future material implementations (such as hyper-elastic materials).

The basis that we provide here will define the element local basis:

  • $ \underline{u_1} $ gets projected on the element segment to give $ \underline{e_1^{elem}} $
  • $ \underline{u_2} $ is orthogonalized with $ \underline{e_1^{elem}} $ to get $ \underline{e_2^{elem}} $
  • $ \underline{e_3^{elem}} $ is computed such that $ (\underline{e_1^{elem}}, \underline{e_2^{elem}}, \underline{e_3^{elem}})$ forms an orthonormal basis.

Comparison with Abaqus B31 element

Eventually, the last part of this article focuses on the comparison of SesamX linear beam element against Abaqus equivalent B31 element. A beam assembly clamped on one end is subjected to an arbitrary load on the second end. The displacements at the nodes, obtained from a linear static resolution, are compared between SesamX and Abaqus.

Model

The model studied for this comparison is made of the following beam assembly.

Abaqus comparison model

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.

Results

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

Model displacement

And the table below gives the comparison of the nodal displacements between Abaqus and SesamX.

SesamX and Abaqus comparison

Although the error seems greater on the displacements compared to the rotations, it is still quite small. SesamX linear beam element implementation is very close to Abaqus implementation.

Following these links you have access to the SesamX input cards as well as the Abaqus input cards that I used for this comparison. I imported the Abaqus mesh and selections in SesamX (to make sure I have the same model description) and then I defined the properties, loads and boundary conditions.

Conclusion

I hope you had a pleasant reading. Starting from version V2020_01 of SesamX you can use the beam element I presented in this article. Feel free to share your thoughts or simply ask for more information!


Did you like this content?

Register to our newsletter and get notified of new articles