"""Higher cantilever bending modes.
Extends the :class:`CantileverNaturalFrequency` fundamental-
mode coverage with the second and third bending natural
frequencies of a clamped-free prismatic beam.
Closed-form quantities
----------------------
Cantilever modes have frequencies
.. math::
f_n = \\frac{(\\beta_n L)^{2}}{2 \\pi L^{2}}
\\sqrt{\\frac{E I}{\\rho A}},
where :math:`\\beta_n L` are the roots of
:math:`1 + \\cos(\\beta L) \\cosh(\\beta L) = 0`:
* :math:`\\beta_1 L = 1.875104068712` — fundamental.
* :math:`\\beta_2 L = 4.694091132974`.
* :math:`\\beta_3 L = 7.854757438238`.
* :math:`\\beta_n L \\to (2n-1)\\pi/2` asymptotically for large
:math:`n`.
References
----------
* Rao, S. S. *Mechanical Vibrations*, 6th ed., Pearson, 2017,
§8.5 Table 8.1 — cantilever bending-mode coefficients.
* Timoshenko, S. P. *Vibration Problems in Engineering*,
4th ed., Wiley, 1974, §5.3.
Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):
* **NAFEMS FV-2** transverse modes of a slender cantilever
(tabulated β_n L values for n = 1..5) — the canonical
industry-neutral reference.
* Equivalent benchmarks ship in established commercial solvers
(cantilever-bending-modes family) but the analytical roots
above are the primary verification target.
Implementation notes
--------------------
Uses the same ``HEX8`` enhanced-strain mesh as
:class:`CantileverNaturalFrequency`. Identifies the n-th
bending mode by post-processing the eigen output: pick the
first mode that is (a) dominated by the non-axial
(transverse) displacement family and (b) has ``n``
antinodes along the top-face centerline. Higher modes get
successively larger tolerances to account for the
Timoshenko shear + rotary-inertia correction that grows with
frequency on a 3D solid mesh.
"""
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,
)
#: First four roots of :math:`1 + \\cos(\\beta L)\\cosh(\\beta L) = 0`.
BETA_L_ROOTS: tuple[float, float, float, float] = (
1.8751040687119611,
4.694091132974174,
7.854757438237613,
10.995540734875467,
)
[docs]
@dataclass
class CantileverHigherModes(BenchmarkProblem):
"""Cantilever second + third bending natural frequencies."""
name: str = "cantilever_higher_modes"
description: str = (
"Clamped-free prismatic beam — 2nd and 3rd transverse "
"bending modes. Euler-Bernoulli closed-form with "
"β_n L from the characteristic equation (Rao 2017 §8.5)."
)
analysis: str = "modal"
default_n_modes: int = 12
#: Cantilever length [m]. Set to 4.0 (h/L = 1/80) so the
#: Euler-Bernoulli closed form is asymptotically exact for the
#: first four bending modes — the FE result tracks the EB
#: prediction to within ~0.2 % across all four modes at the
#: default mesh, instead of drifting up to ~4 % from the
#: Timoshenko shear correction that a stockier beam picks up.
L: float = 4.0
width: float = 0.05
height: float = 0.05
E: float = 2.0e11
nu: float = 0.3
rho: float = 7850.0
@property
def published_values(self) -> tuple[PublishedValue, ...]:
A = self.width * self.height
I = self.width * self.height**3 / 12.0 # noqa: E741
base = np.sqrt(self.E * I / (self.rho * A))
return (
PublishedValue(
name="f2_bending",
value=BETA_L_ROOTS[1] ** 2 / (2.0 * np.pi * self.L**2) * base,
unit="Hz",
source="Rao 2017 §8.5 Table 8.1",
formula="f2 = (β2 L)^2 / (2 π L^2) sqrt(E I / (ρ A))",
tolerance=0.005,
),
PublishedValue(
name="f3_bending",
value=BETA_L_ROOTS[2] ** 2 / (2.0 * np.pi * self.L**2) * base,
unit="Hz",
source="Rao 2017 §8.5 Table 8.1",
formula="f3 = (β3 L)^2 / (2 π L^2) sqrt(E I / (ρ A))",
tolerance=0.005,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
# Default mesh sized for the L = 4 m, h = 0.05 m cantilever:
# 120 axial slices keeps the Galerkin-error h-bound below the
# ~0.5 % tolerance the published_values use for modes 2 and 3.
nx = int(mesh_params.get("nx", 120))
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")
return m