"""Clamped cantilever under transverse tip load — 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` and Poisson's ratio :math:`\\nu = 0.3`. Clamped at
:math:`x = 0`, loaded by total transverse force :math:`P` at
:math:`x = L` (distributed uniformly across the tip-face nodes).
Closed-form quantities
----------------------
For a slender beam (shear-deformation negligible) the
Euler-Bernoulli solution predicts:
* **Tip deflection**
.. math::
\\delta = \\frac{P L^{3}}{3 E I}, \\qquad
I = \\frac{w h^{3}}{12}.
* **Maximum normal stress** at the clamped face
.. math::
\\sigma_\\text{max} = \\frac{P L\\, c}{I}, \\qquad c = h / 2.
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, §9.2 — solid-element shear-locking caveats.
Implementation notes
--------------------
Uses 8-node hex elements in the ``enhanced_strain`` formulation
(Simo-Rifai 1990) to skirt the shear-locking that plain
isoparametric bilinear hexes suffer on slender bending. See
:doc:`/reference/elements/solid185` for the formulation
choice.
"""
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 CantileverEulerBernoulli(BenchmarkProblem):
"""Cantilever tip deflection / root stress vs Euler-Bernoulli."""
name: str = "cantilever_eb"
description: str = (
"Slender clamped cantilever under transverse tip load — "
"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 transverse tip load [N].
P: float = 1.0e3
@property
def published_values(self) -> tuple[PublishedValue, ...]:
I = self.width * self.height**3 / 12.0 # noqa: E741
delta = self.P * self.L**3 / (3.0 * self.E * I)
c = self.height / 2.0
sigma_max = self.P * self.L * c / I
return (
PublishedValue(
name="tip_deflection",
value=delta,
unit="m",
source="Timoshenko 1955 §5.4",
formula="delta = P L^3 / (3 E I)",
# HEX8 enhanced-strain (Simo-Rifai 1990) recovers the
# EB tip deflection to within ~0.5 % on the default
# 40 × 3 × 3 mesh.
tolerance=0.01,
),
PublishedValue(
name="root_stress_max",
value=sigma_max,
unit="Pa",
source="Timoshenko 1955 §5.4",
formula="sigma_max = P L c / I (c = h/2)",
# Root bending stress carries the standard hex-corner
# over-shoot from Gauss-point extrapolation; convergence
# tracks the displacement law.
tolerance=0.10,
),
)
[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]
f_per_node = self.P / len(tip)
for n in tip:
m.apply_force(int(n + 1), fy=f_per_node)
# Stash the tip / clamped node sets on the model so extract()
# can find them without recomputing.
m._bench_tip_nodes = tip # type: ignore[attr-defined]
m._bench_clamped_nodes = clamped # type: ignore[attr-defined]
return m