Cantilever under uniformly distributed load — Euler–Bernoulli closed form

For a clamped-free prismatic cantilever of length \(L\), flexural rigidity \(EI\), and uniform transverse line load \(w\) (force per unit length), the tip deflection is

\[\delta_\text{tip} = \frac{w L^{4}}{8 E I}, \qquad I = \frac{w_b h^{3}}{12}\]

(Timoshenko 1955 §5.4, Cook 2002 §2.4–§2.5). The corresponding slope and curvature at the clamp are

\[\theta_\text{tip} = \frac{w L^{3}}{6 E I}, \qquad \kappa(0) = \frac{M(0)}{E I} = \frac{w L^{2}}{2 E I},\]

so the maximum bending stress at the root extreme fibre is \(\sigma_\text{max} = E\, \kappa(0)\, c = w L^{2} c / (2 I)\).

References

• Timoshenko, S. P. Strength of Materials, Part I: Elementary Theory and Problems, 3rd ed., Van Nostrand, 1955, §5.4.

• Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §2.4 (Hermite cubics) — the cantilever-UDL example is the canonical test of consistent-load integration on the beam element.

Implementation note

The cantilever is meshed as a 3D HEX8 plate at slenderness \(L/h = 20\) with three transverse layers; the enhanced_strain formulation (Simo-Rifai 1990) is selected to skirt the shear-locking that plain bilinear hexes suffer on slender bending. The UDL is distributed across the top-surface nodes as a per-node F_y proportional to the nodal influence area.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver
from femorph_solver import ELEMENTS

Geometry + material

L = 1.0
WIDTH = 0.05
HEIGHT = 0.05
E = 2.0e11
NU = 0.30
RHO = 7850.0

# Total load per unit length [N/m]; integrated over L gives the
# total force for the equivalent point-load comparison.
w = 1000.0
W_total = w * L

I = WIDTH * HEIGHT**3 / 12.0  # noqa: E741
# Closed-form magnitude is wL^4 / (8EI); sign convention here is
# negative because the applied UDL points in -y.
delta_published = -w * L**4 / (8.0 * E * I)

print(f"problem: L={L} m, w={w} N/m (downward), EI={E * I:.4e} N·m²")
print(f"published δ_tip = -w L^4 / (8 E I) = {delta_published:+.6e} m")

problem: L=1.0 m, w=1000.0 N/m (downward), EI=1.0417e+05 N·m²
published δ_tip = -w L^4 / (8 E I) = -1.200000e-03 m

Build a HEX8 plate cantilever

NX, NY, NZ = 40, 3, 3
xs = np.linspace(0.0, L, NX + 1)
ys = np.linspace(0.0, WIDTH, NY + 1)
zs = np.linspace(0.0, HEIGHT, NZ + 1)
grid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing="ij")).cast_to_unstructured_grid()

m = femorph_solver.Model.from_grid(grid)
m.assign(
    ELEMENTS.HEX8(integration="enhanced_strain"),
    material={"EX": E, "PRXY": NU, "DENS": RHO},
)

pts = np.asarray(m.grid.points)
clamped = np.where(pts[:, 0] < 1e-9)[0]
m.fix(nodes=(clamped + 1).tolist(), dof="ALL")

Distribute the UDL as nodal F_y on the top surface

A uniform line load w (force/length) on the beam axis maps to a per-node \(F_y\) proportional to the nodal influence area on the top surface. For a regular structured grid this is constant per interior node and half at the edges. Total distributed force = w * L is invariant.

top = np.where(pts[:, 2] > HEIGHT - 1e-9)[0]
# Equal weights per top-surface node — coarse but acceptable on
# a structured grid for this verification.
fy_per_node = -W_total / len(top)
for n in top:
    m.apply_force(int(n + 1), fy=fy_per_node)

Solve and compare

res = m.solve()
u = np.asarray(res.displacement).reshape(-1, 3)

tip = np.where(pts[:, 0] > L - 1e-9)[0]
delta_computed = float(u[tip, 1].mean())

rel_err = (delta_computed - delta_published) / delta_published
print(f"computed δ_tip               = {delta_computed:+.6e} m")
print(f"relative error vs textbook   = {rel_err * 100:+.2f} %")

# 10 % tolerance at this mesh — HEX8 EAS recovers the
# Euler-Bernoulli answer with ~5-8 % residual at L/h=20 on a
# 40x3x3 grid.
assert abs(rel_err) < 0.10, f"cantilever UDL deflection failed: {rel_err:.2%}"

computed δ_tip               = -1.202263e-03 m
relative error vs textbook   = +0.19 %

Render the deflected shape

warped = m.grid.copy()
warped.points = m.grid.points + u
warped["|u|"] = np.linalg.norm(u, axis=1)

plotter = pv.Plotter(off_screen=True, window_size=(640, 240))
plotter.add_mesh(m.grid, style="wireframe", color="grey", opacity=0.3)
plotter.add_mesh(warped, scalars="|u|", cmap="viridis", show_edges=False)
plotter.view_xy()
plotter.camera.zoom(1.05)
plotter.show()