"""Simply-supported beam — first three bending natural frequencies.
Reference geometry
------------------
Prismatic slender beam of length :math:`L`, cross-section
:math:`w \\times h`, simply supported at both ends (knife-edge
vertical support on the bottom face). We solve for the first
three transverse bending modes.
Closed-form quantities
----------------------
The Euler-Bernoulli equation :math:`EI\\,w'''' + \\rho A\\,\\ddot w
= 0` with SS BCs :math:`w(0)=w(L)=0` and
:math:`w''(0)=w''(L)=0` separates into mode shapes
:math:`w_n(x) = \\sin(n \\pi x / L)` and natural frequencies
.. math::
f_n = \\frac{n^{2}\\pi}{2 L^{2}} \\sqrt{\\frac{E I}{\\rho A}},
\\qquad n = 1, 2, 3, \\ldots
For a square section, :math:`I / A = h^{2} / 12`, so
:math:`f_n` scales with :math:`h / L^{2}` and is independent of
``w`` (width).
References
----------
* Rao, S. S. *Mechanical Vibrations*, 6th ed., Pearson, 2017,
§8.5 — simply-supported Euler-Bernoulli beam.
* Meirovitch, L. *Fundamentals of Vibrations*, Waveland, 2010,
§7.4 — transverse vibration of a uniform beam.
Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):
* **Abaqus AVM 1.6.x** simply-supported-beam NF family.
* **NAFEMS FV-5** transverse natural frequencies of a thin square
plate — 1D limit matches the SS beam formula here.
"""
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 SimplySupportedBeamModes(BenchmarkProblem):
"""Simply-supported beam first three bending natural frequencies."""
name: str = "ss_beam_modes"
description: str = (
"Prismatic simply-supported beam, first three transverse "
"bending natural frequencies — Euler-Bernoulli closed-form "
"(Rao 2017 §8.5)."
)
analysis: str = "modal"
default_n_modes: int = 6
#: Beam length [m]. Set to 4.0 (h/L = 1/80) so the Euler-
#: Bernoulli closed form is asymptotically exact across the first
#: three bending modes — at L = 1 (h/L = 1/20) the 3D solid mesh
#: truthfully captured a Timoshenko shear correction up to ~2.6 %
#: at mode 3.
L: float = 4.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
nu: float = 0.3
rho: float = 7850.0
@property
def published_values(self) -> tuple[PublishedValue, ...]:
I = self.width * self.height**3 / 12.0 # noqa: E741
A = self.width * self.height
base = np.sqrt(self.E * I / (self.rho * A))
coeffs = [1, 2, 3]
vals = []
for n in coeffs:
fn = (n * n * np.pi) / (2.0 * self.L**2) * base
vals.append(
PublishedValue(
name=f"f{n}_bending",
value=fn,
unit="Hz",
source="Rao 2017 §8.5",
formula=f"f{n} = (n^2 pi / 2 L^2) sqrt(E I / (rho A)), n = {n}",
# Slender geometry (h/L = 1/80) keeps Timoshenko
# shear drift below 0.5 % across the first three
# modes; uniform tolerance applies.
tolerance=0.005,
)
)
return tuple(vals)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
# 120 axial slices match the 4 m span; ny / nz unchanged.
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)
# SS supports: pin UZ along the full-width bottom line at
# each end, plus the minimum extra constraints needed to
# kill the axial and lateral rigid-body modes. Pinning UY
# globally would impose plane-strain (stiffening f₁ by ~4 %
# for ν = 0.3); instead suppress UY on a single span-wise
# line and rely on mode-shape filtering in ``extract`` to
# identify the x-z bending modes among the returned
# frequencies.
left_line = np.where((pts[:, 0] < 1e-9) & (pts[:, 2] < 1e-9))[0]
for n in left_line:
m.fix(nodes=int(n + 1), dof="UZ", value=0.0)
left_pin = np.where((pts[:, 0] < 1e-9) & (pts[:, 1] < 1e-9) & (pts[:, 2] < 1e-9))[0]
for n in left_pin:
m.fix(nodes=int(n + 1), dof="UX", value=0.0)
right_line = np.where((pts[:, 0] > self.L - 1e-9) & (pts[:, 2] < 1e-9))[0]
for n in right_line:
m.fix(nodes=int(n + 1), dof="UZ", value=0.0)
# Suppress y-sway: pin UY along the x-axis centreline (y = 0
# edge of the bottom face). This kills the rigid-body y-
# translation + the rigid-body y-rotation about the x-axis
# without imposing plane-strain anywhere off the line.
# The y-bending modes still exist but now have a unique
# mode shape we can filter in ``extract``.
yline = np.where((pts[:, 1] < 1e-9) & (pts[:, 2] < 1e-9))[0]
for n in yline:
m.fix(nodes=int(n + 1), dof="UY", value=0.0)
m._bench_L = self.L # type: ignore[attr-defined]
return m