"""Clamped-clamped beam — first three bending natural frequencies.
Companion to :class:`SimplySupportedBeamModes` and
:class:`CantileverHigherModes`. Both ends clamped — the
characteristic equation has a different root set than the
SS beam (whose roots are simply :math:`n \\pi`).
The Euler-Bernoulli characteristic equation for fully-clamped
(CC) beam free vibration is
.. math::
1 - \\cos(\\beta L) \\cosh(\\beta L) = 0,
with first three roots (Rao 2017 Table 8.1):
* :math:`\\beta_1 L = 4.730040745`
* :math:`\\beta_2 L = 7.853204624`
* :math:`\\beta_3 L = 10.99560784`
The natural frequencies follow
.. math::
f_n = \\frac{(\\beta_n L)^{2}}{2 \\pi L^{2}}
\\sqrt{\\frac{E I}{\\rho A}}.
References
----------
* Rao, S. S. *Mechanical Vibrations*, 6th ed., Pearson, 2017,
§8.5 Table 8.1 — clamped-clamped beam coefficients.
* Timoshenko, S. P. *Vibration Problems in Engineering*,
4th ed., Wiley, 1974, §5.3.
Cross-references (public verification manuals — fair-use
problem-ID citations only; no vendor content vendored):
* **Abaqus AVM 1.6.x** clamped-clamped-beam NF family.
* **NAFEMS FV-3** clamped-clamped beam fundamental + higher
modes.
"""
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 three roots of :math:`1 - \\cos(\\beta L)\\cosh(\\beta L) = 0`.
CC_BETA_L_ROOTS: tuple[float, float, float] = (
4.730040744862704,
7.853204624095838,
10.995607838001671,
)
[docs]
@dataclass
class ClampedClampedBeamModes(BenchmarkProblem):
"""CC beam first three bending natural frequencies."""
name: str = "cc_beam_modes"
description: str = (
"Clamped-clamped prismatic beam, first three transverse "
"bending natural frequencies — Rao 2017 §8.5 Table 8.1."
)
analysis: str = "modal"
default_n_modes: int = 8
#: Beam span [m]. Set to 4.0 (h/L = 1/80) so the Euler-
#: Bernoulli closed form is asymptotically exact for the first
#: three bending modes — at the original L = 1 (h/L = 1/20) the
#: 3D solid mesh truthfully captures a Timoshenko shear /
#: rotary-inertia correction that grows up to ~4 % at mode 3.
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="f1_bending",
value=CC_BETA_L_ROOTS[0] ** 2 / (2.0 * np.pi * self.L**2) * base,
unit="Hz",
source="Rao 2017 §8.5 Table 8.1",
formula="f1 = (β1 L)² / (2 π L²) sqrt(E I / (ρ A))",
tolerance=0.005,
),
PublishedValue(
name="f2_bending",
value=CC_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 π L²) sqrt(E I / (ρ A))",
tolerance=0.005,
),
PublishedValue(
name="f3_bending",
value=CC_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 π L²) sqrt(E I / (ρ A))",
tolerance=0.005,
),
)
[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)
# Clamp both end faces fully.
ends = np.where((pts[:, 0] < 1e-9) | (pts[:, 0] > self.L - 1e-9))[0]
m.fix(nodes=(ends + 1).tolist(), dof="ALL")
return m