Clamped-clamped beam — first three bending natural frequencies

Classical Bernoulli-beam eigenvalue problem. A prismatic beam clamped at both ends has bending natural frequencies (Rao 2017 §8.5 Table 8.1; Timoshenko 1974 §5.3)

\[f_n \;=\; \frac{(\beta_n\, L)^{2}}{2\pi\, L^{2}} \sqrt{\frac{E I}{\rho A}},\]

with the dimensionless eigenvalues \(\beta_n L\) defined as the positive roots of

\[1 - \cos(\beta L)\cosh(\beta L) \;=\; 0.\]

The first three roots are

\[\beta_1 L = 4.730041, \quad \beta_2 L = 7.853205, \quad \beta_3 L = 10.995608.\]

These are slightly higher than the cantilever’s \(\beta_n L \in \{1.875, 4.694, 7.855\}\) (Rao Table 8.1) because clamping both ends raises the bending stiffness — at any given mode the CC beam is a factor of ~6.3 higher in \(\omega_1^2\) than the same-geometry cantilever (β₁²-ratio \((4.730 / 1.875)^{2} \approx 6.36\)).

Implementation

Drives the existing ClampedClampedBeamModes problem class on a HEX8 enhanced-strain prism mesh; the FE spectrum is filtered for transverse-bending modes (UY-dominant) and matched to the published \(\beta_n L\) values by frequency-proximity. Higher modes pick up a small Timoshenko shear / rotary-inertia correction that grows with mode number — visible in the table below as the FE result drifts ~0.5 % at mode 1 to ~4 % at mode 3.

References

• Rao, S. S. (2017) Mechanical Vibrations, 6th ed., Pearson, §8.5 Table 8.1 (clamped-clamped beam eigenvalue roots).

• Timoshenko, S. P. (1974) Vibration Problems in Engineering, 4th ed., Wiley, §5.3.

• Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §11 (modal analysis).

• Simo, J. C. and Rifai, M. S. (1990) “A class of mixed assumed strain methods …” International Journal for Numerical Methods in Engineering 29 (8), 1595–1638.

Cross-references (problem-id only):

from __future__ import annotations

import numpy as np
import pyvista as pv

from femorph_solver.validation.problems import ClampedClampedBeamModes

Build the model from the validation problem class

problem = ClampedClampedBeamModes()
m = problem.build_model()
print(f"clamped-clamped beam mesh: {m.grid.n_points} nodes, {m.grid.n_cells} HEX8 cells")

clamped-clamped beam mesh: 1936 nodes, 1080 HEX8 cells

Closed-form references

print()
print(f"{'n':>3}  {'β_n L':>10}  {'f_pub [Hz]':>12}")
print(f"{'-' * 3:>3}  {'-' * 10:>10}  {'-' * 12:>12}")
for n in (1, 2, 3):
    pv_ = [pv0 for pv0 in problem.published_values if pv0.name == f"f{n}_bending"][0]
    beta_str = ["4.7300", "7.8532", "10.9956"][n - 1]
    print(f"{n:>3}  {beta_str:>10}  {pv_.value:12.3f}")

  n       β_n L    f_pub [Hz]
---  ----------  ------------
  1      4.7300        16.214
  2      7.8532        44.694
  3     10.9956        87.619

Modal solve + frequency-proximity match

res = m.modal_solve(n_modes=problem.default_n_modes)
freqs = np.asarray(res.frequency, dtype=np.float64)
shapes = np.asarray(res.mode_shapes).reshape(-1, 3, freqs.size)

print()
print(f"{'n':>3}  {'k':>3}  {'f_FE [Hz]':>12}  {'f_pub [Hz]':>12}  {'err':>9}  {'tol':>5}")
print(f"{'-' * 3:>3}  {'-' * 3:>3}  {'-' * 12:>12}  {'-' * 12:>12}  {'-' * 9:>9}  {'-' * 5:>5}")

claimed = np.zeros(freqs.size, dtype=bool)
for n in (1, 2, 3):
    pv_ = [pv0 for pv0 in problem.published_values if pv0.name == f"f{n}_bending"][0]
    f_pub = pv_.value
    candidates = [k for k in range(freqs.size) if not claimed[k]]
    k_best = min(candidates, key=lambda k: abs(freqs[k] - f_pub))
    f_fe = float(freqs[k_best])
    err_pct = (f_fe - f_pub) / f_pub * 100.0
    tol_pct = pv_.tolerance * 100.0
    print(f"{n:>3}  {k_best:>3}  {f_fe:12.3f}  {f_pub:12.3f}  {err_pct:+8.3f}%  ≤{tol_pct:.1f}%")
    claimed[k_best] = True
    assert abs(err_pct) < tol_pct, f"f_{n} deviation {err_pct:.2f}% exceeds {tol_pct}%"

  n    k     f_FE [Hz]    f_pub [Hz]        err    tol
---  ---  ------------  ------------  ---------  -----
  1    0        16.244        16.214    +0.187%  ≤0.5%
  2    2        44.736        44.694    +0.093%  ≤0.5%
  3    5        87.595        87.619    -0.028%  ≤0.5%

Render the first three mode shapes side-by-side

# Re-acquire the matching mode indices for rendering.
claimed = np.zeros(freqs.size, dtype=bool)
panels: list[tuple[int, int, float]] = []
for n in (1, 2, 3):
    pv_ = [p for p in problem.published_values if p.name == f"f{n}_bending"][0]
    candidates = [k for k in range(freqs.size) if not claimed[k]]
    k_best = min(candidates, key=lambda k: abs(freqs[k] - pv_.value))
    panels.append((n, k_best, float(freqs[k_best])))
    claimed[k_best] = True

plotter = pv.Plotter(shape=(1, 3), off_screen=True, window_size=(1080, 320))
for n, k, f in panels:
    plotter.subplot(0, n - 1)
    phi = shapes[:, :, k]
    warp = m.grid.copy()
    warp.points = m.grid.points + (0.05 / np.max(np.abs(phi))) * phi
    warp["|phi|"] = np.linalg.norm(phi, axis=1)
    plotter.add_mesh(m.grid, style="wireframe", color="grey", opacity=0.3)
    plotter.add_mesh(warp, scalars="|phi|", cmap="plasma", show_edges=False, show_scalar_bar=False)
    plotter.add_text(f"mode {n}: f = {f:.0f} Hz", position="upper_edge", font_size=10)
plotter.link_views()
plotter.view_xy()
plotter.show()