Timoshenko Beam Elements
This document summarizes the theoretical background of the Timoshenko beam element formulation implemented in numgeo, with distinct bending inertias about local axes \(y\) and \(z\), and distinct shear correction factors \(\kappa_y\) and \(\kappa_z\). The element supports:
- 2-node (linear) and 3-node (quadratic) formulations,
- 2D and 3D configurations (2D beams about one major axis; 3D beams can have bending about two principal axes).
In 3D, each element can have:
- \(I_{yy}\) and \(I_{zz}\) as the area moments of inertia,
- \(\kappa_{y}\) and \(\kappa_{z}\) as distinct shear correction factors in the local \(y\)- and \(z\)-directions,
- \(J\) as the torsional constant,
- \(A\) as the cross‐sectional area,
- \(E\) for axial Young's modulus,
- \(G\) for the shear modulus.
In 2D, typically you only have one bending inertia \(I_{zz}\) (or \(I\)) about the out‐of‐plane axis and one shear factor \(\kappa_z\).
Comparison with Euler-Bernoulli beams
Unlike Euler-Bernoulli beams which neglect shear deformation, Timoshenko beams account for both bending and shear deformations. This makes them particularly suitable for:
- Short, stubby beams (length-to-depth ratio < 10)
- Composite beams with significant shear flexibility
- Applications where shear deformation is significant
1. Kinematics and Displacement Field
2-node Linear Element
A 2-node Timoshenko beam has nodes \(1\) and \(2\). In:
- 2D: each node has \(\left(u_x,\,u_y,\,r_z\right)\). That's 3 DOFs per node, so 6 total.
- 3D: each node has \(\left(u_x,\,u_y,\,u_z,\,r_x,\,r_y,\,r_z\right)\). That's 6 DOFs per node, so 12 total.
The displacements in the local axis \(\ell\) are usually interpolated by linear shape functions:
where \(L\) is the beam length. The transverse displacement and rotation (and in 3D, the out‐of‐plane displacements and rotations) follow a similar linear interpolation.
Shear Deformation
In Timoshenko theory, the cross section is not constrained to remain perpendicular to the neutral axis. The shear strain in each principal direction (say, local \(y\) or local \(z\)) depends on the difference:
in 2D (just one shear). In 3D, we analogously have:
(depending on the chosen local coordinate definitions). Each shear mode has its own correction factor \(\kappa_y\) or \(\kappa_z\).
Total deformation
The total deformation in a Timoshenko beam is a combination of both bending and shear contributions. The total deflection \(w\) at any point can be expressed as:
where \(w_b(x)\) is the deflection due to bending and \(w_s(x)\) is the deflection due to shear. This additive nature distinguishes Timoshenko beams from Euler-Bernoulli beams, which only account for bending deformation.
3-node Quadratic Element
For a 3-node beam we use quadratic shape functions in a natural coordinate \(\xi \in [-1, 1]\). Each node still has the same DOFs as in the linear element.
Node Ordering Convention
In the 3-node quadratic element, the physical node ordering is: - Node 1: Left node (ξ = -1) - Node 2: Right node (ξ = 1) - Node 3: Middle node (ξ = 0)
The shape functions mapped to the physical nodes are:
Derivatives are:
These derivatives are used to form the strain‐displacement matrix \(\mathbf{B}\). Integrating along the beam axis with multiple Gauss points captures the bending and shear distribution more accurately than the simpler 2‐node approach.
Strain-Displacement Matrices
For the 3-node quadratic element, the B-matrices (strain-displacement matrices) relate nodal displacements to strains:
- Axial strain (\(\varepsilon\)):
- Bending curvature (\(\omega\)):
- Shear strain (\(\gamma\)):
This operates on both transverse displacement and rotation DOFs.
- Torsional strain (in 3D):
2. Constitutive Model
The current code in numgeo uses a linear elastic Timoshenko beam model with:
- \(I_{yy}\) and \(I_{zz}\) for bending about the local \(y\) and \(z\) axes,
- \(\kappa_{y}\) and \(\kappa_{z}\) as the shear correction factors in each direction.
Hence, the beam can have different bending rigidities \(E\,I_{yy}\) vs. \(E\,I_{zz}\) – for instance, if the cross section is rectangular with different side lengths in local \(y\) vs. local \(z\). Likewise, shear stiffness is \(G\,\kappa_y\,A\) about \(y\) and \(G\,\kappa_z\,A\) about \(z\).
Additionally, you specify:
- \(E\): Young's modulus (axial/bending).
- \(G\): Shear modulus \(G = \dfrac{E}{2\,(1 + \nu)}\), where \(\nu\) is Poisson's ratio.
- \(A\): Cross-sectional area (for axial and shear).
- \(J\): Torsional constant (for 3D torsion).
2D simulations
In 2D simulations, the model typically just needs \(I_{zz}\) (or "\(I\)") and \(\kappa_{z}\) (the single out‐of‐plane direction). Simply set \(I_{yy}=I_{zz}\) and \(\kappa_{y}=\kappa_{z}\).
Shear & Bending Formula
In standard 2‐node Timoshenko derivations, the local bending–shear block is often expressed with the parameter
But now we do it twice in 3D:
Each bending axis has a corresponding shear axis and separate factors.
3. Internal Force Vector and Stiffness Matrix
2-node Linear Element
- Axial stiffness: \(\frac{E\,A}{\,L\,}\).
- Bending about local \(y\) (in 3D) or local "out‐of‐plane" axis (in 2D) includes the shear factor \(\phi,\psi\). If in 3D, the same logic is repeated for bending about local \(z\).
- Torsion: \(\frac{G\,J}{\,L\,}\) in 3D.
- Shear: \(\left(G\,\kappa_y\,A\right)\) in local‐\(y\), \(\left(G\,\kappa_z\,A\right)\) in local‐\(z\).
Putting these into a local 6×6 matrix for 2D or 12×12 for 3D gives sub‐blocks:
- \(K_{\mathrm{axial}}\),
- \(K_{\mathrm{bend}+y}\),
- \(K_{\mathrm{bend}+z}\),
- \(K_{\mathrm{torsion}}\),
- Coupling terms that appear in Timoshenko (e.g., \(\pm6\,E\,I\,\psi/L^{2}\) blocks, etc.).
Example (3D, 2-node)
You often see sub‐matrices like:
(and so forth) with off‐diagonal terms from Timoshenko coupling. For a 2D 2-node Timoshenko beam element in its local coordinate system, the complete 6×6 stiffness matrix is:
where the DOF ordering is \([u_{x}^1, u_{y}^1, r_{z}^1, u_{x}^2, u_{y}^2, r_{z}^2]\).
3-node Quadratic Element
For the 3-node beam, we do a Gauss integration along the element axis in local coordinate \(\xi\). At each Gauss point:
- The strain–displacement matrix \(\mathbf{B}\) is formed for axial, shear in y, shear in z, bending about y, bending about z, torsion.
- The material (constitutive) matrix \(\mathbf{D}\) uses \((E,\,I_{yy},\,I_{zz},\,\kappa_y,\,\kappa_z,\,J)\).
- Summing over 3 Gauss points yields the final \(\mathbf{K}_{\mathrm{local}}\) and \(\mathbf{f}_{\mathrm{local}}\).
The local stiffness matrix is typically formed by integrating:
Where \(|J|\) is the determinant of the Jacobian matrix (for beams, simply \(|J| = L/2\)).
When using selective reduced integration, the shear terms are integrated using a reduced scheme (typically 1-point for quadratic elements) to prevent shear locking.
4. Coordinate Transformation to Global
Because the above stiffness computations are in the local beam system, the final step is to:
- Build a transformation matrix \(\mathbf{T}_{\mathrm{trans}}\) from local → global.
- In 3D, define a local frame \(\{\ell,\,y,\,z\}\) from the unit tangent \(\mathbf{n}\), plus two perpendicular vectors \(\mathbf{n}_y,\mathbf{n}_z\).
-
In 2D, define \(\mathbf{n}\) and a single perpendicular \(\mathbf{n}_\perp\).
-
Compute
These \(\mathbf{K}_{\mathrm{global}}\) and \(\mathbf{f}_{\mathrm{global}}\) then assemble into the overall global system.
5. Numerical Integration Scheme
For the quadratic (3-node) element, a 3-point Gauss integration scheme is used over the natural coordinate \(\xi \in [-1,1]\). The integration points and weights are:
| Integration Point | \(\xi\) | Weight \(w\) |
|---|---|---|
| 1 | \(-\sqrt{\frac{3}{5}}\) | \(\frac{5}{9}\) |
| 2 | \(0\) | \(\frac{8}{9}\) |
| 3 | \(\sqrt{\frac{3}{5}}\) | \(\frac{5}{9}\) |
The mapping from the natural coordinate \(\xi\) to the physical coordinate \(x\) is given by
with the Jacobian
This Jacobian is used in the numerical integration to transform the integrals from the physical coordinate system to the natural coordinate system. The element stiffness matrix components are calculated as:
Where \(n_{gp}\) is the number of Gauss points, \(w_k\) is the weight at Gauss point \(k\), and \(\frac{L}{2}\) is the Jacobian determinant.
For selective reduced integration of shear terms, a single Gauss point at \(\xi = 0\) with weight \(w = 2\) is commonly used.
6. Output of forces and moments
The naming convention and orientation of structural forces and bending moments are summarised in Figure 2.
The relationships between the different forces / bending moments and the corresponding strains are defined as follows:
Forces and bending moments are calculated at the integration points of the beam elements. If the output is requested at the nodes, extrapolation of the integration point values to the nodes is performed (see Output for more information).
Output forces and bending moments
See the Reference Manual to learn how to request forces and bending moments as output in your simulations.
7. Notes on shear correction factors κ
The Timoshenko beam theory accounts for shear deformation effects that are neglected in the classical Euler-Bernoulli beam theory. To correctly incorporate these shear effects, a shear correction factor (κ) is introduced to adjust the shear stiffness. This factor accounts for the non-uniform distribution of shear stresses across the cross-section and varies significantly depending on the cross-sectional geometry.
Rectangular Cross-Sections
For homogeneous rectangular cross-sections, the shear correction factor is typically taken as κ = 5/6 ≈ 0.8333
This value is derived from the exact solution of the three-dimensional elasticity theory and represents the ratio between the average shear strain energy and the strain energy calculated from the assumed constant shear stress distribution. The factor 5/6 has been widely accepted in engineering practice and incorporated into most standard references and textbooks.
I-Sections and IPE Profiles
For I-shaped sections such as IPE profiles, the shear correction factor is considerably lower than for rectangular sections due to the concentration of material in the flanges and the relatively thin web that carries most of the shear force:
- For strong-axis bending (load perpendicular to the web): κ ≈ 0.3-0.4
- For weak-axis bending (load parallel to the web): κ ≈ 0.5-0.6
These values reflect the reduced effective shear area of I-sections. The precise value depends on:
- The ratio of flange width to web thickness
- The ratio of flange thickness to web thickness
- The depth-to-width ratio of the section
For standard IPE profiles, a value of κ = 0.4 is often appropriate for strong-axis bending calculations.
Other Common Cross-Sections
Circular Sections
- Solid circular section: κ = 0.9
- Thin-walled tube: κ = 0.5-0.6 (depending on the ratio of outer to inner diameter)
Box Sections
- Square hollow section with thin walls: κ = 0.4-0.5
- Rectangular hollow section: κ varies from 0.35-0.5 depending on aspect ratio
Channel Sections
- C-shaped sections: κ ≈ 0.35-0.45 for strong-axis bending
- U-shaped sections: κ ≈ 0.45-0.55 for weak-axis bending
T-Sections
- T-shaped sections: κ ≈ 0.3-0.4 (similar to I-sections but asymmetric)
Angle Sections
- L-shaped sections: κ ≈ 0.45-0.55 (varies with leg length ratio)
Calculation Methods
For more precise applications, the shear correction factor can be calculated using several methods:
-
Energy-based methods (Cowper, 1966): κ = coefficient derived from equating strain energies in the simplified beam theory to the exact three-dimensional theory
-
Stress-based methods (Bach and Baumann): κ = ratio of maximum shear stress to average shear stress
-
Displacement-based methods (Timoshenko): κ = coefficient derived to match deflections from simplified theory with exact solutions
Practical Considerations
-
For preliminary analysis, the values provided above are generally sufficient.
-
For critical applications requiring high precision, finite element analysis or specialized analytical methods should be employed to determine section-specific correction factors.
-
For composite sections or sections with complex geometry, numerical integration over the cross-section may be necessary.
-
The importance of accurate shear correction factors increases as the span-to-depth ratio decreases (L/h < 10), especially for deep beams or short spans where shear deformation becomes significant.
References
- S. P. Timoshenko, Strength of Materials, Part I & II, 3rd ed.
- O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method.
- K.-J. Bathe, Finite Element Procedures.
- R. D. Cook, Concepts and Applications of Finite Element Analysis.