{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pyvista\npyvista.OFF_SCREEN = True\npyvista.set_jupyter_backend('static')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Fixed-free axial rod \u2014 natural frequencies\n\nThe longitudinal natural frequencies of a uniform rod fixed at\none end and free at the other are (Rao 2017 \u00a78.5; Timoshenko\n1974 \u00a75.1)\n\n\\begin{align}\\omega_n = \\frac{(2n - 1)\\pi}{2 L}\\sqrt{\\frac{E}{\\rho}},\n \\qquad n = 1, 2, 3, \\ldots\\end{align}\n\nso the fundamental is $f_1 = \\frac{1}{4 L}\\sqrt{E/\\rho}$,\nthe second $f_2 = 3 f_1$, the third $f_3 = 5 f_1$,\netc. The closed form is the eigenvalue problem of the 1D wave\nequation $\\rho \\ddot u = E u''$ with mixed Dirichlet\n($u(0) = 0$) / Neumann ($u'(L) = 0$) boundary\nconditions.\n\n## References\n* Rao, S. S. *Mechanical Vibrations*, 6th ed., Pearson, 2017,\n \u00a78.5 (longitudinal vibration of rods).\n* Timoshenko, S. P. *Vibration Problems in Engineering*, 4th\n ed., Wiley, 1974, \u00a75.1.\n\n## Implementation note\nWe mesh a slender 3D HEX8 prism with the **fixed end clamped\nin all three translations** and the **free end loaded in\ndisplacement only along the axial direction** (the cross-section\ncontraction from Poisson is left free). The first transverse\nbending mode appears below the axial fundamental at\n$L / h \\gtrsim 100$; for $L/h = 10$ (the geometry\nbelow) the bending and axial modes don't cross, so the lowest\naxial-dominant mode in the modal sweep is the analytical\n$f_1$.\n\n## Vendor cross-references\n\n================================ ===================== ===================================\nSource Reported f_1 [Hz] Problem ID / location\n================================ ===================== ===================================\nClosed form (wave equation) 1 261.886 Rao 2017 \u00a78.2\nMeirovitch (2010) \u00a76.6 1 261.886 fixed-free rod\nMAPDL Verification Manual 1 262 VM-47 *Torsional freq of a uniform shaft*\nAbaqus Verification Manual \u2248 1 261.9 AVM 1.6.x bar-NF family\n================================ ===================== ===================================\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport math\n\nimport numpy as np\nimport pyvista as pv\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Geometry + material\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L = 1.0\nA = 0.01 * 0.01 # cross-section area, m\u00b2\nWIDTH = HEIGHT = math.sqrt(A)\nE = 2.0e11\nNU = 0.30\nRHO = 7850.0\n\nc_axial = math.sqrt(E / RHO)\nf1_published = c_axial / (4.0 * L)\n\nprint(f\"problem: L={L} m, c_axial=\u221a(E/\u03c1)={c_axial:.1f} m/s\")\nprint(f\"f_1 = c / (4 L) = {f1_published:.4f} Hz\")\nfor n in (2, 3, 4):\n print(f\"f_{n} = (2n-1) f_1 = {(2 * n - 1) * f1_published:.4f} Hz\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Build a 3D prism with one axial direction much finer than the\ntransverse ones so axial waves are well-resolved.\n----------------------------------------------------------------\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "NX, NY, NZ = 40, 2, 2\nxs = np.linspace(0.0, L, NX + 1)\nys = np.linspace(0.0, WIDTH, NY + 1)\nzs = np.linspace(0.0, HEIGHT, NZ + 1)\ngrid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n ELEMENTS.HEX8,\n material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)\n\npts = np.asarray(m.grid.points)\nclamped = np.where(pts[:, 0] < 1e-9)[0]\nm.fix(nodes=(clamped + 1).tolist(), dof=\"ALL\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Modal solve and filter for the axial-dominant mode\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_modes = 12\nres = m.solve_modal(n_modes=n_modes)\nfreqs = np.asarray(res.frequency, dtype=np.float64)\nshapes = np.asarray(res.mode_shapes).reshape(-1, 3, n_modes)\n\n# Pick the mode whose UX RMS dominates UY + UZ \u2014 that's the\n# axial-stretching mode; UY/UZ-dominant modes are the bending\n# pair (which lie above the axial fundamental at L/h=10).\naxial_idx = None\nfor k in range(n_modes):\n ux = float(np.sqrt((shapes[:, 0, k] ** 2).mean()))\n uy = float(np.sqrt((shapes[:, 1, k] ** 2).mean()))\n uz = float(np.sqrt((shapes[:, 2, k] ** 2).mean()))\n if ux > 2.0 * max(uy, uz):\n axial_idx = k\n break\nassert axial_idx is not None, f\"no axial-dominant mode in first {n_modes}\"\n\nf_computed = float(freqs[axial_idx])\nrel_err = (f_computed - f1_published) / f1_published\nprint()\nprint(f\"axial mode index in spectrum : {axial_idx} of {n_modes}\")\nprint(f\"computed f_axial : {f_computed:.4f} Hz\")\nprint(f\"relative error vs Rao 2017 : {rel_err * 100:+.2f} %\")\n\n# 2 % tolerance \u2014 HEX8 axial mode at NX=40 elements per\n# wavelength is well within float64-eigenvalue precision of the\n# textbook answer.\nassert abs(rel_err) < 0.02, f\"axial f_1 failed: {rel_err:.2%}\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the deflected first axial mode\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "phi = shapes[:, :, axial_idx]\nwarped = m.grid.copy()\nscale = 0.05 * L / max(np.abs(phi).max(), 1e-30)\nwarped.points = m.grid.points + scale * phi\nwarped[\"|phi|\"] = np.linalg.norm(phi, axis=1)\n\nplotter = pv.Plotter(off_screen=True, window_size=(640, 240))\nplotter.add_mesh(m.grid, style=\"wireframe\", color=\"grey\", opacity=0.3)\nplotter.add_mesh(warped, scalars=\"|phi|\", cmap=\"plasma\", show_edges=False)\nplotter.view_xy()\nplotter.camera.zoom(1.05)\nplotter.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 0 }