.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/verification/example_verify_axial_rod_nf.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_verification_example_verify_axial_rod_nf.py: Fixed-free axial rod — natural frequencies ========================================== The longitudinal natural frequencies of a uniform rod fixed at one end and free at the other are (Rao 2017 §8.5; Timoshenko 1974 §5.1) .. math:: \omega_n = \frac{(2n - 1)\pi}{2 L}\sqrt{\frac{E}{\rho}}, \qquad n = 1, 2, 3, \ldots so the fundamental is :math:`f_1 = \frac{1}{4 L}\sqrt{E/\rho}`, the second :math:`f_2 = 3 f_1`, the third :math:`f_3 = 5 f_1`, etc. The closed form is the eigenvalue problem of the 1D wave equation :math:`\rho \ddot u = E u''` with mixed Dirichlet (:math:`u(0) = 0`) / Neumann (:math:`u'(L) = 0`) boundary conditions. References ---------- * Rao, S. S. *Mechanical Vibrations*, 6th ed., Pearson, 2017, §8.5 (longitudinal vibration of rods). * Timoshenko, S. P. *Vibration Problems in Engineering*, 4th ed., Wiley, 1974, §5.1. Implementation note ------------------- We mesh a slender 3D HEX8 prism with the **fixed end clamped in all three translations** and the **free end loaded in displacement only along the axial direction** (the cross-section contraction from Poisson is left free). The first transverse bending mode appears below the axial fundamental at :math:`L / h \gtrsim 100`; for :math:`L/h = 10` (the geometry below) the bending and axial modes don't cross, so the lowest axial-dominant mode in the modal sweep is the analytical :math:`f_1`. .. GENERATED FROM PYTHON SOURCE LINES 40-51 .. code-block:: Python from __future__ import annotations import math import numpy as np import pyvista as pv import femorph_solver from femorph_solver import ELEMENTS .. GENERATED FROM PYTHON SOURCE LINES 52-54 Geometry + material ------------------- .. GENERATED FROM PYTHON SOURCE LINES 54-70 .. code-block:: Python L = 1.0 A = 0.01 * 0.01 # cross-section area, m² WIDTH = HEIGHT = math.sqrt(A) E = 2.0e11 NU = 0.30 RHO = 7850.0 c_axial = math.sqrt(E / RHO) f1_published = c_axial / (4.0 * L) print(f"problem: L={L} m, c_axial=√(E/ρ)={c_axial:.1f} m/s") print(f"f_1 = c / (4 L) = {f1_published:.4f} Hz") for n in (2, 3, 4): print(f"f_{n} = (2n-1) f_1 = {(2 * n - 1) * f1_published:.4f} Hz") .. rst-class:: sphx-glr-script-out .. code-block:: none problem: L=1.0 m, c_axial=√(E/ρ)=5047.5 m/s f_1 = c / (4 L) = 1261.8862 Hz f_2 = (2n-1) f_1 = 3785.6585 Hz f_3 = (2n-1) f_1 = 6309.4308 Hz f_4 = (2n-1) f_1 = 8833.2031 Hz .. GENERATED FROM PYTHON SOURCE LINES 71-74 Build a 3D prism with one axial direction much finer than the transverse ones so axial waves are well-resolved. ---------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 74-91 .. code-block:: Python NX, NY, NZ = 40, 2, 2 xs = np.linspace(0.0, L, NX + 1) ys = np.linspace(0.0, WIDTH, NY + 1) zs = np.linspace(0.0, 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, material={"EX": E, "PRXY": NU, "DENS": RHO}, ) pts = np.asarray(m.grid.points) clamped = np.where(pts[:, 0] < 1e-9)[0] m.fix(nodes=(clamped + 1).tolist(), dof="ALL") .. GENERATED FROM PYTHON SOURCE LINES 92-94 Modal solve and filter for the axial-dominant mode -------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 94-125 .. code-block:: Python n_modes = 12 res = m.modal_solve(n_modes=n_modes) freqs = np.asarray(res.frequency, dtype=np.float64) shapes = np.asarray(res.mode_shapes).reshape(-1, 3, n_modes) # Pick the mode whose UX RMS dominates UY + UZ — that's the # axial-stretching mode; UY/UZ-dominant modes are the bending # pair (which lie above the axial fundamental at L/h=10). axial_idx = None for k in range(n_modes): ux = float(np.sqrt((shapes[:, 0, k] ** 2).mean())) uy = float(np.sqrt((shapes[:, 1, k] ** 2).mean())) uz = float(np.sqrt((shapes[:, 2, k] ** 2).mean())) if ux > 2.0 * max(uy, uz): axial_idx = k break assert axial_idx is not None, f"no axial-dominant mode in first {n_modes}" f_computed = float(freqs[axial_idx]) rel_err = (f_computed - f1_published) / f1_published print() print(f"axial mode index in spectrum : {axial_idx} of {n_modes}") print(f"computed f_axial : {f_computed:.4f} Hz") print(f"relative error vs Rao 2017 : {rel_err * 100:+.2f} %") # 2 % tolerance — HEX8 axial mode at NX=40 elements per # wavelength is well within float64-eigenvalue precision of the # textbook answer. assert abs(rel_err) < 0.02, f"axial f_1 failed: {rel_err:.2%}" .. rst-class:: sphx-glr-script-out .. code-block:: none axial mode index in spectrum : 11 of 12 computed f_axial : 1263.0670 Hz relative error vs Rao 2017 : +0.09 % .. GENERATED FROM PYTHON SOURCE LINES 126-128 Render the deflected first axial mode ------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 128-141 .. code-block:: Python phi = shapes[:, :, axial_idx] warped = m.grid.copy() scale = 0.05 * L / max(np.abs(phi).max(), 1e-30) warped.points = m.grid.points + scale * phi warped["|phi|"] = np.linalg.norm(phi, axis=1) plotter = pv.Plotter(off_screen=True, window_size=(640, 240)) plotter.add_mesh(m.grid, style="wireframe", color="grey", opacity=0.3) plotter.add_mesh(warped, scalars="|phi|", cmap="plasma", show_edges=False) plotter.view_xy() plotter.camera.zoom(1.05) plotter.show() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_axial_rod_nf_001.png :alt: example verify axial rod nf :srcset: /gallery/verification/images/sphx_glr_example_verify_axial_rod_nf_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.215 seconds) .. _sphx_glr_download_gallery_verification_example_verify_axial_rod_nf.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_verify_axial_rod_nf.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_verify_axial_rod_nf.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_verify_axial_rod_nf.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_