Timoshenko Beam Element for Dynamic Analysis

Authors: Guisela Eloisa Ortega Chulde, Anastazja Broniatowska

Supervisors: Felix Schneider, Sebastian Schopper

Introduction

Modelling and simulation of structures are fundamental steps in engineering. Systems of equations are turned into ordinary differential equations (ODE) analytically or numerically and solved via various algorithms. Usually, to solve high derivative ODEs, these are turned into first derivative order matrices. As the system becomes more complex, the required order of the model gets higher. One of the most common structures to model are beams. There exist many techniques for the modelling of the mechanical behavior of a beam. Among these techniques, the Timoshenko Beam theory and the Euler Bernoulli theory are the most prominent.

This project implements a finite element model for the dynamic behavior of the Timoshenko beam. The model is generated in MATLAB and adapted into the FEM toolbox.  The finite element model for the dynamic behavior of the Timoshenko Beam can be used as a benchmark of parametric Model Order Reduction. This becomes useful and practical when reducing the order of a model to maintain the accuracy of the results. 

Timoshenko Beam Theory Theoretical Formulation

Unlike the Euler-Bernoulli theory, the Timoshenko beam theory considers the shear deformation from more pronounce bending. In this case, the shear deformation  $//<![CDATA[ \begin{array}{l}\tau (x,y)\end{array} //]]>$ is not longer assume as 0. Thus, the beam is no perpendicular to its normal axis after deformation (see Figure 1). The Timoshenko Beam theory includes cases related with small deflections based on shear deformations for short beams. In contrast, the Euler-Bernoulli beam theory neglects shear deformations and focus on the behavior of flexure dominated beams. It is expected that in a mid-length range both theories are equivalent. 

Generally, the Euler Bernoulli theory is utilized for its simplicity as it models a stiffer beam. Bernoulli tends to overestimate the deformation of the beam. However, shorter beams are better approximated with the Timoshenko Beam Theorem. 

General FEM Model

The displacement field for the Timoshenko beam for pure bending is:

where $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ is the transverse deflection and $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$ is the rotation of a transverse normal line about the y axis. For the dynamic case, the transverse deflection is the displacement of the mid-surface in the z axis and the rotation is the angle between the normal to the mid-surface of the beam. The angular deflection in this theory is a variable and not a slope approximation of the deflection like in the Euler-Bernoulli formulation. From the previous assumption the Timoshenko beam theory can be described using the principal of minimum total potential energy. The corresponding strains and stresses are derivate as follow:

For linear elastic, isotropic, homogeneous beam of constant cross-section, the equilibrium is reached at the following equations:

In matrix form: 

For normal modes (3)(4) can be solved, but its numerical approximation can present some difficulties. When modelling thin structures, a locking phenomenon presents during the application of finite element methods such as the Timoshenko Beam theory. 

Approaches to solve Locking Problem

Reduced Integration Element

The most well-known approach to alleviate locking in the Timoshenko beam element is the reduced integration of the element. Here we interpolate displacement and rotation fields by the linear shape functions. This initially yields the following equilibrium equation:

Where:

For the thin beam limit, we should receive no shear force, which yields additional constraint, namely:

Enforcing this constraint leads to the evaluation of zero bending energy to satisfy all the equations from (6), which of course is a trivial solution.

A simple way to overcome this problem is the reduced integration of $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$ fields. In numerical solutions, integration is done with the use of Gaussian quadrature, where for a linear function two points are needed for the exact evaluation. In the reduced integration element (RIE), w is integrated with two quadrature points and $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$ with only one for evaluation of the shear strain. This leads to a changed stiffness matrix, namely:

This formulation removes the problem of locking but also has its flaws. The element has additional zero-energy modes and the evaluated displacements are still not exact (although, with enough elements, this drawback can be minimized). For the Reduced Integration Element, the full equilibrium equation for the dynamic case looks as follows:

where $//<![CDATA[ \begin{array}{l}I_0 = \rho A, I_2 = \frac{\rho Ah^2}{12}\end{array} //]]>$ for a rectangular cross-section, constant along the x axis. 

Consistent Interpolation Element

Another idea to resolve the shear locking effect is a consistent interpolation of $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ and  $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$ fields. As already mentioned, at least linear interpolation of these fields is needed to suffice the continuity requirement.  Nevertheless, the Kirchhoff constraint requires w to have one order higher continuity, namely quadratic interpolation. This way 5x5 matrix is received with an additional degree of freedom $//<![CDATA[ \begin{array}{l}w_c\end{array} //]]>$. The system can be reduced to 4x4 by eliminating this DOF, yielding the same stiffness matrix as for RIE, but a modified load vector, so that the equilibrium equation looks as follows:

Due to algebraic manipulations done to reduce the system, this element cannot be directly applied to a dynamic case and will not be further discussed.

Interdependent Interpolation Element

Another way to use consistent interpolation is to approximate $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ with a cubic function and $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$ with a quadratic function. Here the classical interpolation would lead to a 7x7 stiffness matrix, which then again would require some manipulations to reduce the system to 4x4. Instead, Hermite cubic splines are used for interpolation of $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$, and a related quadratic interpolation for $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$. This method yields a super convergent element. As for any distribution of the transverse shear load $//<![CDATA[ \begin{array}{l}q(x)\end{array} //]]>$ with constant bending and shear stiffness, It results in the exact nodal displacements. In the thin beam limit, the equations reduce to the Euler - Bernoulli formulations without any additional constraints. For the static case the equilibrium equation looks as follows:

And for a dynamic case, the mass term is added using following matrix:

Compared to the RIE,  the Interdependent Interpolation Element shows slightly superior behavior for pure flexural frequencies, yielding more accurate results for the same number of elements, nevertheless, RIE also converges to the correct result. As for pure shear natural frequencies, IIE performs worse and does not converge to the correct solution. 

Assumed Strain - Displacement models

Assumed Strain - Displacement models are based on a different approach. Here displacement fields $//<![CDATA[ \begin{array}{l}w, \phi\end{array} //]]>$ are treated independently from strain fields, namely $//<![CDATA[ \begin{array}{l}\kappa, \gamma\end{array} //]]>$. 

ASD - LLCC

For this model, both displacement fields are interpolated with linear functions, while strain fields are constant across the element. This yields the following static equilibrium equation:

Where $//<![CDATA[ \begin{array}{l}\Omega\end{array} //]]>$ is defined same as in previous cases. After some reformulations, it is apparent that the equation is exactly the same as for RIE. This model can also be directly applicable to dynamic analysis.

ASD HQLC

Second possibility is to use Hermite cubic interpolation of displacement $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$, quadratic interpolation of rotation $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$, linear interpolation of the curvature $//<![CDATA[ \begin{array}{l}\kappa\end{array} //]]>$, and constant interpolation of shear strain $//<![CDATA[ \begin{array}{l}\gamma\end{array} //]]>$. This formulation again requires system size reduction, as it must contain a middle node for quadratic interpolation. Therefore, the formulation is not appropriate to implement in dynamic analysis and will not be discussed further. 

MATLAB implementation of the 3D Reduced Integration Element 

For the implementation a Reduced Integration Element formulation has been chosen, as it yields good results both for static and dynamic analysis, while the evaluation of mass and stiffness matrices remains quite simple. The beam element has been modeled to carry the load in all three dimensions. Therefore, each node has 6 degrees of freedom and a nodal displacement vector is defined as follows $//<![CDATA[ \begin{array}{l}\underline{u_i}=\{u_i, v_i, w_i, \phi_x_i, \phi_y_i, \phi_z_i\}^T\end{array} //]]>$.

As the implementation uses a linear interpolation of all displacement fields, each element has two nodes, which means that element matrices are of size 12x12. For the axial and torsional DOFs, the classical formulation has been implemented. These types of loading influence DOFs no. 1, 7, and 4, 10 respectively, and will not be discussed further in the scope of this project. The RIE formulation has been implemented for bending in two axes, so that DOFs no. 2, 6, 8, and 12 are influenced by bending around the z-axis, and DOFs no. 3, 5, 9, and 11 by bending around the y-axis. 

To validate the implementation's performance, tests for static and dynamic analysis have been created in MATLAB.

Testing and Validation

Static Test - Tip Load

The first test, for validation of the static behavior, is a cantilever beam with a point load applied at the tip (see Figure 2). Since matrix entries for bending along y- or z-axis are the same, for the sake of simplicity, only the scenario of a tip load pointing in the negative z-direction is considered. 

For the given load case, the bending moment and the shear force are given by:

Euler - Bernoulli solution

For the Euler - Bernoulli formulation  where the angle of rotation $//<![CDATA[ \begin{array}{l}\phi_y = -\frac{dw}{dx}\end{array} //]]>$, we have $//<![CDATA[ \begin{array}{l}M(x)=-EI\frac{d^2w}{dx^2}\end{array} //]]>$, while the shear force is simply$//<![CDATA[ \begin{array}{l}Q(x) = \frac{dM(x)}{dx}\end{array} //]]>$  (where $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ is the displacement in the z-direction). Finally, after the integration with respect to x, we get the function of z-displacement along the beam given by:

Timoshenko Solution

Timoshenko beam formulation additionally takes into account the shear deformation, so that: $//<![CDATA[ \begin{array}{l}\phi_y = -\frac{dw}{dx}+\gamma\end{array} //]]>$. Now we have the bending moment and the shear force defined as follows: $//<![CDATA[ \begin{array}{l}M(x) = -EI \frac{d\phi}{dx}, Q(x)=\kappa GA(\phi+ \frac{dw}{dx})\end{array} //]]>$. Integrating first equation, to get $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$ and substituting it in the second equation and, finally, integrating again, we receive the formula for the displacement:

Thin Beam Limit

As visible from equations (16) and (17), the solutions of the tip displacement for the two formulations differ by a term $//<![CDATA[ \begin{array}{l}\frac {PL}{\kappa AG}\end{array} //]]>$, therefore, as a condition of a thin beam limit, where both solutions should yield approximately the same value, we define:

which then leads to a factor:

Test results

The tests have been performed for the following sets of parameters:

As visible, the bending stiffness $//<![CDATA[ \begin{array}{l}EI_y\end{array} //]]>$ and cross-section area $//<![CDATA[ \begin{array}{l}A\end{array} //]]>$have been kept constant, and the shear correction factor has been chosen as for a rectangular cross-section $//<![CDATA[ \begin{array}{l}\kappa=\frac{5}{6}\end{array} //]]>$. A convergence study has been performed to verify, whether the formulation yields the results of the Euler - Bernoulli formulation in the thin beam limit. 

From Figure 4, it can be observed that with an increasing number of elements the solution converges properly and that for thicker beams the displacement is getting greater than the Euler - Bernoulli result. This can be easily explained, by the additional shear term in the Timoshenko solution. Based on these results, we define $//<![CDATA[ \begin{array}{l}c=1e-3\end{array} //]]>$ as a limiting value for thin beams, where the Euler - Bernoulli solution is valid with an error below 0.1% (see Figure 5). It has been also verified that the element implementation converges to Timoshenko analytical solution for thicker beams.

Eigenfrequencies Test

To verify if the implemented element can be used for dynamic analyses additional tests have been performed, the first of which was the eigenfrequency test. The cantilever beam has been chosen for verification and the first ten eigenfrequencies checked. To simplify the solution and check only the bending formulation of the implemented element, all degrees of freedom apart from z-displacement and rotation around the y-axis, have been fixed. This way we can receive only bending eigenmodes. Since the Timoshenko beam shows softer behavior, than Euler - Bernoulli, it is expected that the eigenfrequencies will be lower.

Euler - Bernoulli Formulation

The equation of motion for free vibrations is given as:

The boundary conditions for a cantilever beam are:

The equation can be solved by the separation of variables, which together with the boundary condition yields an equation

which can be solved numerically for $//<![CDATA[ \begin{array}{l}\lambda\end{array} //]]>$.  The wave number is defined:  $//<![CDATA[ \begin{array}{l}\lambda^4 = \frac{\rho A}{EI_y}\omega^2\end{array} //]]>$, from which eigenfrequencies can be computed, as: $//<![CDATA[ \begin{array}{l}f_i = \frac{\omega_i}{2\pi}\end{array} //]]>$. This way we get:

Test results

Tests in the thin beam limit have been performed with the following parameters:

To visualize the results, the absolute error with respect to Euler - Bernoulli solution has been plotted against the eigenfrequency number and implementations with a different number of elements have been shown. 

It is clearly visible, that with an increasing number of elements the solution converges properly, and the error increases with the eigenfrequency number. Nevertheless, for 90 elements the relative error is below 0.5%. Already for this thin beam, we see a small deviation to smaller eigenfrequencies, because of the softer behavior of the Timoshenko beam. 

Thick Beam Case

The thick beam case was compared to the results of the same structure discretized with Reissner - Mindlin plate elements. The analyzed structure was of size 10x1x0.3 m. For the Timoshenko beam, the structure is discretized with 100 elements with a checked convergence. For Reissner - Mindlin elements, the discretization was with 20x2 elements, but the convergence could not be verified, as with more elements, their slenderness was becoming too small to consider elements as plates. Therefore, the comparison is not very reliable.

Dynamic Test - Tip Load

For the dynamic case, the weak form from (5) needs to be modified. The equilibrium equation correspond to the displacement and mixed FEM models of the Timoshenko beam theory. Due to the presence of the second-derivative terms of $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\phi\end{array} //]]>$, it is not possible to manipulate the equations in algebraic form as in the static case.

RIE Model for Dynamic Analysis

By incorporating the RIE method for the dynamic case, the FEM model results as follow: 

Where $//<![CDATA[ \begin{array}{l}I_0 = \rho A, I_2 = \frac {\rho A h^2}{12}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\rho\end{array} //]]>$ represents the mass density of the material. 

Test Results

Fort the dynamic case, the parameters from the static case (see Figure 3) are used as well for consistency.  Figure 9 shows the tip displacement for each case value of c from the thin beam limit theory. The displacement results from the dynamic test yields very accurate results once enough elements are incorporated. In general, the relative error of the tip displacement reaches below %0.1 for thinner beams at $//<![CDATA[ \begin{array}{l}c\leqslant 5.85e^{-3}\end{array} //]]>$ . The thinner the beam and the more elements are used in the dynamic tests, the model convergence performs better. 

Conclusion

The Timoshenko beam theory models a more accurate representation of true bending whereas the Euler-Bernoulli theory models a stiffer beam.  The strain-displacement formulation of the theory allows for the decomposition of the transverse deflection into bending and shear deflections. If the shear modulus in the beam approaches to infinity, the beam becomes rigid. Moreover, if rotational inertia is neglected too, the Timoshenko beam theory would converge into simple beam theory.

For the locking phenomenon, there are several methods to solve it once it appears at the modeling of thin beams. Some of the methods require algebraic manipulations to reduce the system of equations. However, the reduced Integration element method (RIE) shows a superior performance with a simpler formulation of the mass matrix. The eigenfrequency results indicates RIE can predict flexural and pure shear frequencies accurately once a certain number of elements are incorporated to the model.  This convergence is not monotonic, specially at higher frequencies. 

In the thin beam limit, the element equations of the Timoshenko beam theory with RIE reduces to only one constraint (i.e. Kirchhoff condition). However, for the dynamic case it does not yield into exact displacements as for the static one. Although, the beam element does not lock, a sufficient number of elements is needed in the dynamic case to obtain accurate deflections. Moreover, the thinner the beam, the more accurate the results yield. The implemented element shows a softer behavior than the Euler-Bernoulli element. 

References

[1] Fogang, V. (2020). Timoshenko beam theory exact solution for bending, second-order analysis, and stability. https://doi.org/10.20944/preprints202011.0457.v1

[2] Peuscher, Heiko & Hubele, Jörg & Eid, Rudy & Lohmann, Boris. (2009). Generating a Parametric Finite Element Model of a 3D Cantilever Timoshenko Beam Using MATLAB.

[3] Reedy J N. On the dynamic behavior of the Timoshenko beam finite elements.