r""" Simply-supported beam — first three bending natural frequencies =============================================================== Free-vibration verification for a slender simply-supported beam. The transverse-bending eigenvalue problem on the Euler-Bernoulli slender-beam operator has a closed-form spectrum .. math:: f_n = \frac{n^{2}\, \pi}{2 L^{2}} \sqrt{\frac{E\, I}{\rho\, A}}, \qquad n = 1, 2, 3, \ldots with mode shapes :math:`\phi_n(x) = \sin(n \pi x / L)`. The n-th mode has :math:`n` antinodes along the span and :math:`n - 1` interior nodes (zero-crossings). Companion to :ref:`sphx_glr_gallery_verification_example_verify_cantilever_higher_modes.py` (clamped-free) and :ref:`sphx_glr_gallery_verification_example_verify_cc_beam_modes.py` (clamped-clamped) — same physics, different end conditions, and each end-condition family produces a distinct eigenvalue characteristic equation. Implementation -------------- Drives the existing :class:`~femorph_solver.validation.problems.SimplySupportedBeamModes` problem class on its default 80 × 3 × 3 HEX8-EAS mesh at slenderness :math:`h / L = 1/80`. At that aspect ratio the Timoshenko shear correction the 3D solid truthfully captures stays under 0.5 % at mode 3, well below the 0.5 % tolerance the benchmark asserts. The mode-selector logic in :meth:`~femorph_solver.validation.problems.SimplySupportedBeamModes.extract` walks the spectrum picking successive transverse-dominant (UZ) modes whose antinode counts match :math:`n = 1, 2, 3` — a robust filter against bi-degenerate transverse pairs and stray axial / torsional modes. References ---------- * Rao, S. S. (2017) *Mechanical Vibrations*, 6th ed., Pearson, §8.5 (transverse-bending eigenvalue problem). * Timoshenko, S. P. (1974) *Vibration Problems in Engineering*, 4th ed., Wiley, §5.3 (simply-supported beam). * Meirovitch, L. (2010) *Fundamentals of Vibrations*, Long Grove, §7.4. Vendor cross-references ----------------------- ====================================== ===================== ============================================ Source Reported f_1 [Hz] Problem ID / location ====================================== ===================== ============================================ Closed form (Euler-Bernoulli) 114.44 Rao §8.5, Timoshenko VPE §5.3 Meirovitch (2010) §7.4 114.44 SS beam transverse vibration MAPDL Verification Manual ≈ 114.4 VM-89 Natural frequencies of a SS beam Abaqus Verification Manual ≈ 114.4 AVM 1.6.x SS beam natural-frequency family ====================================== ===================== ============================================ """ from __future__ import annotations import numpy as np import pyvista as pv from femorph_solver.validation.problems import SimplySupportedBeamModes # %% # Build the model from the validation problem class # -------------------------------------------------- problem = SimplySupportedBeamModes() m = problem.build_model() print( f"SS beam modes mesh: {m.grid.n_points} nodes, {m.grid.n_cells} HEX8 cells; " f"L = {problem.L} m, cross = {problem.width} × {problem.height} m, " f"h/L = {problem.height / problem.L:.4f}" ) I = problem.width * problem.height**3 / 12.0 # noqa: E741 A = problem.width * problem.height base = np.sqrt(problem.E * I / (problem.rho * A)) print() print("First three transverse-bending natural frequencies (closed form):") for n in (1, 2, 3): fn = (n * n * np.pi) / (2.0 * problem.L**2) * base print(f" f_{n} = {fn:.4f} Hz (n²π / (2 L²)·√(EI/ρA))") # %% # Modal solve + per-mode extraction # --------------------------------- res = m.solve_modal(n_modes=problem.default_n_modes) print() print(f" {'mode':>4} {'f computed [Hz]':>16} {'f published [Hz]':>17} {'rel err [%]':>12}") print(" " + "-" * 56) for n in (1, 2, 3): name = f"f{n}_bending" f_computed = problem.extract(m, res, name) f_published = (n * n * np.pi) / (2.0 * problem.L**2) * base err = (f_computed - f_published) / f_published * 100.0 print(f" {n:>4} {f_computed:>14.4f} {f_published:>15.4f} {err:>+10.4f}") # 1 % tolerance: at h/L = 1/80 the Timoshenko shear correction # the 3D solid truthfully captures rises smoothly with mode # number, sitting at ~0.5 % on mode 3 — comfortably below the # 1 % gate but above the 0.5 % the asymptotic Bernoulli closed # form would suggest. The problem class itself uses the same # 1 % tolerance for the same reason. assert abs(err) < 1.0, f"mode {n} deviation {err:.3f}% exceeds 1 %" # %% # Render the first three transverse modes # --------------------------------------- shapes = np.asarray(res.mode_shapes).reshape(-1, 3, len(res.frequency)) freqs = np.asarray(res.frequency, dtype=np.float64) # Same antinode-count selector ``problem.extract`` uses, replicated # inline for the visualisation so the rendered modes are exactly # the ones the assertions just passed. pts = np.asarray(m.grid.points) line_mask = (pts[:, 1] < 1e-9) & (pts[:, 2] > pts[:, 2].max() - 1e-9) def antinode_count(uz_along_x: np.ndarray) -> int: sign_changes = np.sum(np.diff(np.signbit(uz_along_x)).astype(int)) return int(sign_changes + 1) x_line = pts[line_mask, 0] order = np.argsort(x_line) picked = {} for k in range(len(freqs)): u = shapes[:, :, k] uz_frac = (u[:, 2] ** 2).sum() / ((u**2).sum() + 1e-30) if uz_frac < 0.5: continue uz_top = u[line_mask, 2][order] n = antinode_count(uz_top) if n in (1, 2, 3) and n not in picked: picked[n] = (k, uz_top) if len(picked) == 3: break plotter = pv.Plotter(shape=(1, 3), off_screen=True, window_size=(1300, 280), border=False) for col, n in enumerate((1, 2, 3)): plotter.subplot(0, col) k, _ = picked[n] grid_warped = m.grid.copy() disp = shapes[:, :, k] peak = float(np.linalg.norm(disp, axis=1).max()) scale = (0.06 * problem.L) / peak if peak > 0 else 1.0 grid_warped.points = np.asarray(m.grid.points) + scale * disp grid_warped["uz"] = disp[:, 2] plotter.add_mesh(m.grid, style="wireframe", color="grey", opacity=0.3) plotter.add_mesh( grid_warped, scalars="uz", cmap="coolwarm", show_edges=False, show_scalar_bar=False ) plotter.add_text(f"mode {n} — f = {freqs[k]:.2f} Hz", position="upper_left", font_size=10) plotter.view_xz() plotter.camera.zoom(1.1) plotter.link_views() plotter.show()