"""Cantilever under tip torque — Saint-Venant torsion.
First torsion benchmark in the corpus. A prismatic cantilever
of length :math:`L`, cross-section :math:`h \\times h`, clamped
at one end, loaded by a pure tip torque :math:`T` about the
beam axis. Saint-Venant torsion gives the tip twist angle
.. math::
\\varphi_{\\text{tip}} = \\frac{T L}{G J},
with :math:`G = E / (2(1+\\nu))` the shear modulus and
:math:`J` the torsion constant of the cross-section.
For a **rectangular** cross-section :math:`b \\times t`
(``b ≥ t``), Saint-Venant's infinite series simplifies to
.. math::
J = \\beta(b/t) \\, b \\, t^{3},
where :math:`\\beta` approaches :math:`1/3` for long thin
sections and takes tabulated values for finite aspect ratios
(Timoshenko & Goodier §109). For a **square** section
``b = t``, :math:`\\beta \\approx 0.141`, giving
:math:`J \\approx 0.141 \\, h^{4}`.
Uses the ``BEAM2`` kernel — femorph-solver's 3D linear
beam with 6 DOFs per node. The torsional stiffness uses the
user-supplied ``J`` real constant directly, so the comparison
is a straightforward closed-form verification of the kernel's
torsion block.
References
----------
* Saint-Venant, A. J. C. B. *Mémoire sur la torsion des
prismes.* Académie des Sciences, 1853 — original torsion
theory.
* Timoshenko, S. P. and Goodier, J. N. *Theory of Elasticity*,
3rd ed., McGraw-Hill, 1970, §§104–109 — standard exposition,
including the :math:`\\beta` table for rectangles.
Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):
* **Abaqus AVM 1.5.x** cantilever-torsion family.
Implementation notes
--------------------
1D beam mesh — a straight line of ``BEAM2`` elements along
the global x-axis. Clamp all six DOFs at the left end, apply
``MX = T`` on the tip node (moment about the beam axis). The
twist angle is simply the tip node's ``ROTX`` component.
"""
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,
)
def _saint_venant_beta(b_over_t: float) -> float:
"""Saint-Venant torsion factor :math:`\\beta(b/t)` for a rectangle.
Tabulated values interpolated from Timoshenko & Goodier §109:
``β → 1/3`` as ``b/t → ∞``; ``β ≈ 0.141`` at ``b/t = 1``
(square).
"""
table = {
1.0: 0.141,
1.2: 0.166,
1.5: 0.196,
2.0: 0.229,
2.5: 0.249,
3.0: 0.263,
4.0: 0.281,
5.0: 0.291,
10.0: 0.312,
100.0: 0.333,
}
ratios = sorted(table.keys())
if b_over_t <= ratios[0]:
return table[ratios[0]]
if b_over_t >= ratios[-1]:
return table[ratios[-1]]
for i in range(len(ratios) - 1):
r0, r1 = ratios[i], ratios[i + 1]
if r0 <= b_over_t <= r1:
w = (b_over_t - r0) / (r1 - r0)
return table[r0] + w * (table[r1] - table[r0])
return 0.141 # defensive
[docs]
@dataclass
class CantileverTorsion(BenchmarkProblem):
"""Cantilever tip twist under a pure tip torque."""
name: str = "cantilever_torsion"
description: str = (
"Cantilever beam with a pure tip torque — Saint-Venant "
"closed-form tip-twist φ = T L / (G J)."
)
analysis: str = "static"
#: Cantilever length [m].
L: float = 1.0
#: Cross-section width [m] (``b`` in the Saint-Venant β table).
width: float = 0.05
#: Cross-section height [m] (``t`` in the Saint-Venant β table).
height: float = 0.05
#: Young's modulus [Pa].
E: float = 2.0e11
nu: float = 0.3
rho: float = 7850.0
#: Tip torque about the beam axis [N·m].
T: float = 10.0
@property
def published_values(self) -> tuple[PublishedValue, ...]:
b = max(self.width, self.height)
t = min(self.width, self.height)
beta = _saint_venant_beta(b / t)
J = beta * b * t**3
G = self.E / (2.0 * (1.0 + self.nu))
phi = self.T * self.L / (G * J)
return (
PublishedValue(
name="tip_twist",
value=phi,
unit="rad",
source="Timoshenko & Goodier 1970 §109",
formula="φ = T L / (G J), J = β(b/t) b t^3",
tolerance=0.02,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
n_elem = int(mesh_params.get("n_elem", 20))
xs = np.linspace(0.0, self.L, n_elem + 1)
pts = np.column_stack((xs, np.zeros_like(xs), np.zeros_like(xs)))
cells = np.empty((n_elem, 3), dtype=np.int64)
cells[:, 0] = 2
cells[:, 1] = np.arange(n_elem)
cells[:, 2] = np.arange(1, n_elem + 1)
grid = pv.UnstructuredGrid(
cells.ravel(),
np.full(n_elem, 3, dtype=np.int64), # VTK_LINE
pts,
)
m = femorph_solver.Model.from_grid(grid)
# BEAM2 real constants: (AREA, IZZ, IYY, J). Use the
# same β table as ``published_values`` for consistency.
b = max(self.width, self.height)
t = min(self.width, self.height)
beta = _saint_venant_beta(b / t)
J = beta * b * t**3
A = self.width * self.height
Iz = self.width * self.height**3 / 12.0
Iy = self.height * self.width**3 / 12.0
m.assign(
ELEMENTS.BEAM2,
material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
real=(A, Iz, Iy, J),
)
# Clamp all 6 DOFs at the left end.
m.fix(nodes=1, dof="ALL")
# Apply tip torque (moment about the beam x-axis).
m.apply_force(n_elem + 1, mx=self.T)
m._bench_tip_node_id = n_elem + 1 # type: ignore[attr-defined]
return m