Cantilever natural frequency — Rao 2017 Table 8.1

First transverse bending frequency of a clamped-free prismatic beam:

\[f_{1} = \frac{(\beta_{1} L)^{2}}{2 \pi L^{2}} \sqrt{\frac{E I}{\rho A}}, \qquad \beta_{1} L = 1.8751.\]

Reference: Rao, S. S. Mechanical Vibrations, 6th ed., Pearson, 2017, §8.5, Table 8.1 first entry.

from __future__ import annotations

import numpy as np
import pyvista as pv

from femorph_solver.validation import ConvergenceStudy
from femorph_solver.validation.problems import CantileverNaturalFrequency

Build + solve + convergence study

problem = CantileverNaturalFrequency()
study = ConvergenceStudy(
    problem=problem,
    refinements=[
        {"nx": 20, "ny": 3, "nz": 3},
        {"nx": 40, "ny": 3, "nz": 3},
        {"nx": 80, "ny": 3, "nz": 3},
    ],
)
records = study.run()
rec = records[0]
pub = rec.results[0].published

print(f"{problem.name}: {problem.description}")
print(f"\n  source : {pub.source}")
print(f"  formula: {pub.formula}")
print(f"  published: {pub.value:.6e} {pub.unit}\n")
print(f"  {'mesh':<22s} {'n_dof':>8s} {'f_1 (Hz)':>14s} {'rel err':>10s}  pass")
for r in rec.results:
    mesh_str = "  ".join(f"{k}={v}" for k, v in r.mesh_params.items())
    pass_str = "✓" if r.passed else "✗"
    print(
        f"  {mesh_str:<22s} {r.n_dof:>8d} {r.computed:+14.4f} "
        f"{r.rel_error * 100:+9.2f}%   {pass_str}"
    )
if rec.convergence_rate is not None:
    print(f"\n  fitted rate: |err| ∝ n_dof^(-{rec.convergence_rate:.2f})")

cantilever_nf: First transverse bending frequency of a clamped-free prismatic beam — Rao 2017 §8.5 Table 8.1.

  source : Rao 2017 §8.5 Table 8.1
  formula: f1 = (beta1 L)^2 / (2 pi L^2) sqrt(E I / (rho A))
  published: 4.076904e+01 Hz

  mesh                      n_dof       f_1 (Hz)    rel err  pass
  nx=20  ny=3  nz=3          1008       +41.0046     +0.58%   ✓
  nx=40  ny=3  nz=3          1968       +40.8698     +0.25%   ✓
  nx=80  ny=3  nz=3          3888       +40.8137     +0.11%   ✓

  fitted rate: |err| ∝ n_dof^(-1.20)

Render the first bending mode shape

m = problem.build_model(nx=40, ny=3, nz=3)
modal = m.modal_solve(n_modes=10)

shapes = np.asarray(modal.mode_shapes).reshape(-1, 3, modal.mode_shapes.shape[-1])
# Re-run the mode filter from the problem's extract to pick the
# actual first-bending mode (the square cross-section makes the
# UY/UZ bending pair degenerate).
selected = None
for k in range(shapes.shape[-1]):
    ux_rms = float(np.sqrt((shapes[:, 0, k] ** 2).mean()))
    uy_rms = float(np.sqrt((shapes[:, 1, k] ** 2).mean()))
    uz_rms = float(np.sqrt((shapes[:, 2, k] ** 2).mean()))
    if max(uy_rms, uz_rms) > 2.0 * ux_rms:
        selected = k
        break
if selected is None:
    selected = 0

phi = shapes[:, :, selected]
warped = m.grid.copy()
# Normalise the mode shape so the tip amplitude is ~10% of L,
# just for a visible figure.
scale = 0.1 * problem.L / max(np.linalg.norm(phi, axis=1).max(), 1.0e-30)
warped.points = m.grid.points + scale * phi
warped["|phi|"] = np.linalg.norm(phi, axis=1)

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(m.grid, style="wireframe", color="#888888", opacity=0.3)
plotter.add_mesh(warped, scalars="|phi|", cmap="plasma", show_edges=False)
plotter.view_isometric()
plotter.camera.zoom(1.1)
plotter.show()

Acceptance

finest = rec.results[-1]
assert finest.passed, (
    f"finest-mesh f1 failed: {finest.rel_error:.2%} "
    f"above tolerance {finest.published.tolerance:.0%}"
)