Euler-Bernoulli Beam Element

Last updated: 01 July 2026

1. Hypothesis

The Euler-Bernoulli beam theory states that cross-sections remain plane and perpendicular to the neutral axis after deformation. This implies that transverse shear deformation is zero. The theory directly describes axial deformation and bending; torsion is not part of the original kinematic hypothesis but can be superimposed independently via Saint-Venant’s torsion theory. The following derivations assume linear elasticity and small deformations. Extension to large deformations will be discussed later.

2. Kinematics

This section establishes the relationship between the nodal degrees of freedom and the strain field within the element. The fundamental strategy is to treat axial, bending, and torsion as decoupled deformation modes, apply their respective kinematic assumptions, and then superimpose the results using the principle of linear superposition.

2.1 Nodal Degrees of Freedom

The element possesses 12 degrees of freedom (DOFs), with 6 DOFs assigned to each of the two nodes:

The physical interpretations are:

• $u$: Axial displacement along the local $x$-axis.
• $v, w$: Transverse displacements along the local $y$- and $z$-axes.
• $\theta_x$: Rotational angle about the local $x$-axis (torsion).
• $\theta_y, \theta_z$: Rotational angles about the local $y$- and $z$-axes (bending slopes).

2.2 Kinematic Assumptions

To construct the displacement field, we define the governing assumptions for each deformation mode independently.

Axial Deformation

The axial displacement $u(x)$ is uniform over the cross-section and varies only along the beam axis.

Bending (xy-plane and xz-plane)

The Euler-Bernoulli beam theory states that plane cross-sections remain plane and perpendicular to the neutral axis after bending deformation. This implies that transverse shear deformation due to bending is zero.

For the $xy$-plane (bending about $z$), the slope of the neutral axis equals the cross-sectional rotation:

For the $xz$-plane (bending about $y$), following the right-hand rule, the rotation relates to the slope by:

The negative sign in $\theta_y$ is a consequence of the right-hand rule: a positive rotation about the $y$-axis corresponds to a negative slope of the transverse displacement $w$.

Torsion

The torsion deformation is superimposed independently. Adopting the Saint-Venant torsion theory, we assume:

• The cross-section rotates as a rigid body about the $x$-axis by an angle $\theta_x(x)$.
• The warping (out-of-plane axial displacement) is uniform along the beam axis, meaning it does not produce any axial normal strain $\varepsilon_{xx}$.
• The only strain components introduced by torsion are the shear strains $\gamma_{xy}$ and $\gamma_{xz}$ in the cross-section.

2.3 Displacement Field

The total displacement of a point $(x, y, z)$ is the superposition of contributions from each deformation mode. Cross-sectional rotations about the three axes produce additional displacements for points off the neutral axis ($y \neq 0$ or $z \neq 0$).

The signs follow from the right-hand rule.

2.4 Strain

With the displacement field established, we now derive the strain components. The infinitesimal strain tensor is defined as $\varepsilon_{ij} = \frac{1}{2}(\frac{\partial U_i}{\partial x_j} + \frac{\partial U_j}{\partial x_i})$.

Normal Strain

The axial normal strain is obtained by differentiating $U_x$ with respect to $x$:

Substituting the Euler-Bernoulli slope relations $\theta_z = v'$ and $\theta_y = -w'$:

Shear Strains (Critical Distinction)

We explicitly compute the full shear strain expressions before applying the bending assumptions.

For the $xy$-plane:

Notice the bracket $\left(-\theta_z + \frac{dv}{dx}\right)$. According to the Euler-Bernoulli hypothesis in Section 2.2, $\theta_z = dv/dx$. Therefore, the transverse shear deformation due to bending is identically zero:

Thus, the shear strain in the $xy$-plane reduces to:

Similarly, for the $xz$-plane:

Again, the bracket $\left(\theta_y + \frac{dw}{dx}\right)$ represents the transverse shear deformation due to bending in the $xz$-plane. Since $\theta_y = -dw/dx$, it vanishes:

Leaving only the torsional contribution:

Summary of Strains

The final strain field utilized for the element formulation is:

Remarks: The transverse shear strains generated by bending forces are completely omitted in this theory. The only shear strains present in the cross-section are entirely due to the Saint-Venant torsion. This confirms that the axial/bending and torsion modes are fully decoupled at the strain energy level.

4. Strain Energy and Weak Form

4.1 Strain Energy Functional

Total strain energy $U = \frac12 \int_V \boldsymbol{\sigma}^T \boldsymbol{\varepsilon}\,dV$. With $\sigma_{xx} = E\varepsilon_{xx}$ and $\tau = G\gamma$:

Integrating term by term over the cross-section:

Axial ($\int_A dA = A$): \(\frac12\int_0^L EA\,(u')^2 dx\)

xy-bending ($\int_A y^2 dA = I_z$): \(\frac12\int_0^L EI_z\,(v'')^2 dx\)

xz-bending ($\int_A z^2 dA = I_y$): \(\frac12\int_0^L EI_y\,(w'')^2 dx\)

Torsion ($\int_A r^2 dA = J$, the torsion constant): \(\frac12\int_0^L GJ\,(\theta_x')^2 dx\)

Important: For circular cross-sections, $J = I_p = I_y + I_z$. For non-circular sections, $J < I_p$ (determined by the Saint-Venant torsion solution). The stiffness matrix uses $J$, while the mass matrix uses $I_p$ — this distinction is critical and is detailed in Section 7.3.

The four modes are fully decoupled:

4.2 External Virtual Work

4.3 Weak Form

Applying the principle of virtual work $\delta U - \delta W = 0$ and integrating by parts yields the symmetric weak form:

5. Shape Functions

Introduce the natural coordinate $\xi = x/L \in [0, 1]$, so $dx = L d\xi$ and $d/dx = (1/L)\,d/d\xi$.

5.1 Axial and Torsion — $C^0$ Linear Lagrange Interpolation

Both $u(x)$ and $\theta_x(x)$ require only $C^0$ continuity (first derivative not needed in the weak form):

Properties: $N_1(0)=1, N_1(1)=0$; $N_2(0)=0, N_2(1)=1$.

5.2 Bending — $C^1$ Hermite Cubic Interpolation

The bending weak form contains $v''$, requiring $C^1$ continuity. A cubic polynomial $v(\xi) = a_0 + a_1\xi + a_2\xi^2 + a_3\xi^3$ is determined by four boundary conditions:

Solving for $a_0, a_1, a_2, a_3$ and grouping by DOFs yields:

Verification (only $v_1=1$, all other DOFs zero at $\xi=0$):

• $H_1(0) = 1$, $H_1'(0) = 0$ $\Rightarrow$ $v=1, \theta_z=0$ at node 1 $\checkmark$
• $H_2(0) = 0$, $H_2'(0)/L = 1$ $\Rightarrow$ $v=0, \theta_z=1$ at node 1 $\checkmark$

For xz-plane bending, since $\theta_y = -dw/dx$, the shape functions for $\theta_y$ are $-H_2$ and $-H_4$.

6. Element Stiffness Matrix

6.1 Axial Stiffness

Strain-displacement matrix $\mathbf{B}_a$:

Element stiffness matrix:

6.2 Torsional Stiffness

Identical form, replacing $EA$ with $GJ$:

6.3 xy-Plane Bending Stiffness (about z-axis)

Curvature-displacement matrix $\mathbf{B}_{bz}$:

Second derivatives of $H_i$ with respect to $\xi$:

Therefore:

The stiffness matrix is obtained by integration:

Example entry $k_{11}$:

All entries:

6.4 xz-Plane Bending Stiffness (about y-axis)

Since $\theta_y = -dw/dx$, the shape functions for $\theta_y$ are $-H_2$ and $-H_4$:

This sign change propagates to the stiffness matrix:

Note the sign difference from $\mathbf{k}_{bz}$: $v$-$\theta_z$ coupling is positive ($+6EI_z/L^2$), $w$-$\theta_y$ coupling is negative ($-6EI_y/L^2$).

6.5 Complete 12×12 Stiffness Matrix

With DOF ordering $[u_1, v_1, w_1, \theta_{x1}, \theta_{y1}, \theta_{z1}, u_2, v_2, w_2, \theta_{x2}, \theta_{y2}, \theta_{z2}]$:

$\mathbf{K}_{11}$ (Node 1 × Node 1):

$\mathbf{K}_{12}$ (Node 1 × Node 2):

$\mathbf{K}_{22}$ (Node 2 × Node 2):

Rigid-Body Verification

For $y$-direction rigid translation $[0,1,0,0,0,0,0,1,0,0,0,0]^T$:

• $v$-force at node 1: $\frac{12EI_z}{L^3}\cdot1 + (-\frac{12EI_z}{L^3})\cdot1 = 0$ $\checkmark$
• Moment about $z$ at node 1: $\frac{6EI_z}{L^2}\cdot1 + (-\frac{6EI_z}{L^2})\cdot1 = 0$ $\checkmark$

All rigid-body modes produce zero internal forces. The matrix is symmetric ($\mathbf{K}_{11}^T = \mathbf{K}_{11}$, $\mathbf{K}_{22}^T = \mathbf{K}_{22}$).

7. Element Consistent Mass Matrix

7.1 Kinetic Energy Formulation

The consistent mass matrix is derived from the kinetic energy using the same shape functions as the stiffness matrix. The kinetic energy is:

With $\mathbf{U} = \mathbf{N}\dot{\mathbf{u}}^{(e)}$:

7.2 Axial Mass

$u(\xi) = N_1 u_1 + N_2 u_2$, kinetic energy:

Mass matrix:

Using $\int_0^1 (1-\xi)^2 d\xi = \frac13$, $\int_0^1 \xi(1-\xi)d\xi = \frac16$, $\int_0^1 \xi^2 d\xi = \frac13$:

7.3 Torsional Mass — Fundamental Derivation

The kinetic energy of torsion comes from the rotational motion of the cross-section about the $x$-axis. The velocity field is $\dot{U}_y = -z\,\dot{\theta}_x$, $\dot{U}_z = y\,\dot{\theta}_x$. The squared velocity magnitude is:

Integrating over the cross-section:

Critical distinction:

• For stiffness ($GJ$): $J$ is the torsion constant, determined by the Saint-Venant torsion problem (stress-based). For non-circular sections, $J < I_p$.
• For mass ($\rho I_p$): $I_p = I_y + I_z$ is the polar moment of inertia (geometry-only, independent of stress). It is always simply $\int(y^2+z^2)dA$.
• For circular sections: $J = I_p$ (they coincide).
• For non-circular sections: $J \neq I_p$. Example: a rectangular section $b \times h$ has $I_p = b^3h/12 + bh^3/12$, while $J \approx bh^3\left[\frac13 - 0.21\frac{h}{b}(1 - \frac{h^4}{12b^4})\right]$ — always less than $I_p$.

Using the linear shape functions $N_1, N_2$ for $\theta_x(\xi) = (1-\xi)\theta_{x1} + \xi\theta_{x2}$:

Therefore:

The torsional off-diagonal is $+\rho I_p L/6$, identical in sign to the axial off-diagonal since both use the same linear shape functions.

7.4 xy-Plane Bending Mass

$v(\xi) = H_1 v_1 + H_2 \theta_{z1} + H_3 v_2 + H_4 \theta_{z2}$:

Mass matrix entries $m_{ij} = \rho AL\int_0^1 H_i(\xi)H_j(\xi)\,d\xi$. The Hermite integrals are:

Assembled matrix:

7.5 xz-Plane Bending Mass

Since $\theta_y = -dw/dx$, the shape functions for $\theta_y$ are $-H_2$ and $-H_4$, flipping the signs of all $w$-$\theta_y$ coupling terms relative to the xy-plane:

7.6 Complete 12×12 Mass Matrix

Assembled in block form with $m_{total} = \rho A L$ and $r_x^2 = I_p/A = (I_y + I_z)/A$:

$\mathbf{M}_a$ (Node 1 × Node 1):

$\mathbf{M}_b$ (Node 1 × Node 2 coupling):

$\mathbf{M}_c$ (Node 2 × Node 2):

$\mathbf{M}_b[3,3] = +r_x^2/6$, not $-r_x^2/6$. Both axial and torsional off-diagonal coupling masses are positive, since both use the same linear Lagrange shape functions with identical sign conventions.

Rigid-Body Verification

Torsion ($\dot{\theta}_{x1} = \dot{\theta}_{x2} = \omega$):

Continuum: $T = \frac12\int_0^L \rho I_p \omega^2 dx = \frac12 \rho I_p L \omega^2$ $\checkmark$

Axial ($\dot{u}_1 = \dot{u}_2 = \omega$):

Continuum: $T = \frac12 \rho A L \omega^2$ $\checkmark$

y-direction translation ($\dot{v}_1 = \dot{v}_2 = \omega$, $\dot{\theta}_{z1} = \dot{\theta}_{z2} = 0$):

where $\frac{26}{35} + \frac{18}{70} = \frac{52}{70} + \frac{18}{70} = 1$ $\checkmark$

8. 3D Coordinate Transformation

8.1 Local Coordinate System

The local $x$-axis is along the element:

The local $z$-axis is perpendicular to $\mathbf{e}_1$ and lies in the $\mathbf{e}_1$-$\mathbf{v}$ plane, where $\mathbf{v}$ is a reference vector (typically the global $Z$-axis $[0,0,1]^T$):

8.2 Vertical Element Singularity

When $\mathbf{e}_1 \parallel \mathbf{v}$ (vertical element), $\mathbf{v} \times \mathbf{e}_1 = \mathbf{0}$. Use the global $X$-axis $[1,0,0]^T$ as the reference instead:

8.3 Transformation Matrix

The 3×3 rotation matrix from global to local coordinates:

Properties: $\mathbf{T}_i^T = \mathbf{T}_i^{-1}$, $\det(\mathbf{T}_i) = 1$.

The 12×12 block-diagonal transformation matrix:

Matrix transformation:

9. Global Assembly

Using the direct stiffness method, element matrices are summed into the global system via DOF index mapping:

Numerical implementation using precomputed index arrays:

global_mat.at[I, J].add(local_mat.flatten())

Dirichlet boundary conditions are applied by zeroing rows/columns and setting diagonal entries to 1:

10. Common Implementation Errors

The stiffness matrix is correct in both cases — it correctly uses $GJ$.