"""Clamped cantilever under tip moment — Euler-Bernoulli beam.
Reference geometry
------------------
Slender cantilever of length :math:`L`, rectangular cross-section
:math:`w \\times h`, isotropic elastic material with Young's
modulus :math:`E`. Clamped at :math:`x = 0`, loaded by a pure
tip moment :math:`M` at :math:`x = L`.
The moment is distributed as a linearly-varying axial force on
the tip face — compression on the lower fibres, tension on the
upper fibres — producing a pure bending moment without a
resultant shear or axial force. This is the saint-Venant
equivalent of the textbook "pure tip moment" problem and the
integral recovers :math:`M = \\int_A \\sigma_{xx} \\, z \\, dA`
exactly for any cross-section that the tip traction linearly
samples.
Closed-form quantities
----------------------
For a slender beam (Euler-Bernoulli) the curvature is uniform
along the span (pure moment), so both the tip deflection and
the tip rotation follow directly from
:math:`\\kappa = M / (E I)`:
* **Tip deflection**
.. math::
\\delta = \\frac{M L^{2}}{2 E I}, \\qquad
I = \\frac{w h^{3}}{12}.
* **Tip rotation**
.. math::
\\theta = \\frac{M L}{E I}.
References
----------
* Timoshenko, S. P. *Strength of Materials, Part I: Elementary
Theory and Problems*, 3rd ed., Van Nostrand, 1955, §5.4 —
pure-moment cantilever (integrate :math:`EI w'' = M` twice).
* Gere, J. M. and Goodno, B. J. *Mechanics of Materials*, 9th
ed., Cengage, 2018, §9.3 — standard cantilever deflection /
rotation tables for distributed loadings and tip couples.
Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):
* **Abaqus AVM 1.5.x** cantilever-with-end-moment family.
Implementation notes
--------------------
The moment is applied as a linearly-varying axial traction on
the tip face: for a hex-grid mesh the tip-face nodes carry
``fx_i = k * (z_i - z_mid)`` with ``k`` chosen so the resultant
moment equals :math:`M`. On a structured grid the closed-form
:math:`k = M / \\sum_i z_i^2` gives the exact resultant. This
avoids the reader building a coupling equation or loading a
``CLOAD`` with a single rotational DOF that HEX8 (3-DOF
kernel, no rotational DOF) cannot carry.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Any
import numpy as np
import pyvista as pv
import femorph_solver
from femorph_solver import ELEMENTS
from femorph_solver.validation._benchmark import (
BenchmarkProblem,
PublishedValue,
)
[docs]
@dataclass
class CantileverTipMoment(BenchmarkProblem):
"""Cantilever tip deflection + rotation under a pure tip moment."""
name: str = "cantilever_tip_moment"
description: str = (
"Slender clamped cantilever under a pure tip moment — "
"Euler-Bernoulli closed-form (Timoshenko 1955 §5.4)."
)
analysis: str = "static"
#: Cantilever length [m].
L: float = 1.0
#: Cross-section width [m].
width: float = 0.05
#: Cross-section height [m].
height: float = 0.05
#: Young's modulus [Pa].
E: float = 2.0e11
#: Poisson's ratio.
nu: float = 0.3
#: Density [kg/m^3] — irrelevant for the static solution but
#: required by the material stamp.
rho: float = 7850.0
#: Total tip moment [N·m] (about the horizontal neutral axis).
M: float = 50.0
@property
def published_values(self) -> tuple[PublishedValue, ...]:
I = self.width * self.height**3 / 12.0 # noqa: E741
delta = self.M * self.L**2 / (2.0 * self.E * I)
theta = self.M * self.L / (self.E * I)
return (
PublishedValue(
name="tip_deflection",
value=delta,
unit="m",
source="Timoshenko 1955 §5.4",
formula="delta = M L^2 / (2 E I)",
# Solid hex EAS on a pure-moment problem converges
# faster than tip-force bending; 5 % at a coarse
# mesh is a conservative floor.
tolerance=0.05,
),
PublishedValue(
name="tip_rotation",
value=theta,
unit="rad",
source="Timoshenko 1955 §5.4",
formula="theta = M L / (E I)",
# Tip rotation is extracted from the through-
# thickness displacement gradient; accept 8 %
# on the reference mesh.
tolerance=0.08,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
nx = int(mesh_params.get("nx", 40))
ny = int(mesh_params.get("ny", 3))
nz = int(mesh_params.get("nz", 3))
xs = np.linspace(0.0, self.L, nx + 1)
ys = np.linspace(0.0, self.width, ny + 1)
zs = np.linspace(0.0, self.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": self.E, "PRXY": self.nu, "DENS": self.rho},
)
pts = np.asarray(m.grid.points)
clamped = np.where(pts[:, 0] < 1e-9)[0]
m.fix(nodes=(clamped + 1).tolist(), dof="ALL")
tip = np.where(pts[:, 0] > self.L - 1e-9)[0]
# Apply a linearly-varying axial traction ``fx_i = k * (z_i - z_mid)``
# that integrates to a pure moment M about the y-axis. For a
# discrete node set the exact ``k`` comes from setting
# ``M = sum_i (k * (z_i - z_mid)) * (z_i - z_mid) = k * Σ(z_i - z_mid)^2``.
z_tip = pts[tip, 2]
z_mid = 0.5 * (z_tip.min() + z_tip.max())
dz = z_tip - z_mid
denom = float(np.sum(dz * dz))
# Sign convention: we want ``+M`` to produce ``+z`` tip
# deflection. For a cantilever clamped at ``x = 0`` with
# free end at ``x = L``, this requires the upper fibres
# (``z > z_mid``) to be in compression and the lower
# fibres in tension — so ``fx_i = -k * (z_i - z_mid)``
# with ``k = M / Σ(z_i - z_mid)^2``.
k = self.M / denom
for n, z in zip(tip, z_tip, strict=True):
fx = -k * (z - z_mid)
m.apply_force(int(n + 1), fx=fx)
m._bench_tip_nodes = tip # type: ignore[attr-defined]
m._bench_clamped_nodes = clamped # type: ignore[attr-defined]
return m