{ "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# Cantilever beam \u2014 first four bending modes\n\nClassical Bernoulli-beam eigenvalue problem. A clamped-free\nprismatic cantilever has bending natural frequencies\n(Rao 2017 \u00a78.5 Table 8.1; Timoshenko 1974 \u00a75.3)\n\n\\begin{align}f_n = \\frac{(\\beta_n\\, L)^{2}}{2\\pi\\, L^{2}}\n \\sqrt{\\frac{E I}{\\rho A}},\\end{align}\n\nwith the dimensionless eigenvalues $\\beta_n L$ defined as\nthe positive roots of\n\n\\begin{align}1 + \\cos(\\beta L)\\cosh(\\beta L) = 0.\\end{align}\n\nThe first four roots are\n\n\\begin{align}\\beta_1 L = 1.875104, \\quad\n \\beta_2 L = 4.694091, \\quad\n \\beta_3 L = 7.854757, \\quad\n \\beta_4 L = 10.995541.\\end{align}\n\n## Implementation\n\nA long-aspect HEX8 enhanced-strain prism (40 axial \u00d7 3 \u00d7 3,\nsquare cross-section $b = h$) clamped at one end. The\nmode-shape filter selects the lowest mode whose UY (or UZ \u2014 the\nsquare section is bi-degenerate) RMS dominates the axial UX\nresponse, then walks up the spectrum picking successive\ntransverse-dominant modes.\n\nWe take the cantilever **slender enough** ($h / L = 1/80$)\nthat the Timoshenko shear / rotary-inertia correction the FE\ntruthfully captures stays under 0.5 % at mode 4 \u2014 bench geometries\nlike $h / L = 1/20$ would let the correction grow to ~4 %\nat the same mode, which is real physics that EB doesn't capture\nbut which would obscure the verification pass.\n\n## References\n\n* Rao, S. S. (2017) *Mechanical Vibrations*, 6th ed., Pearson,\n \u00a78.5 Table 8.1 (cantilever-beam eigenvalue roots).\n* Timoshenko, S. P. (1974) *Vibration Problems in Engineering*,\n 4th ed., Wiley, \u00a75.3.\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)\n *Concepts and Applications of Finite Element Analysis*, 4th\n ed., Wiley, \u00a711 (modal analysis).\n* Simo, J. C. and Rifai, M. S. (1990) \"A class of mixed\n assumed-strain methods \u2026\", *IJNME* 29 (8), 1595\u20131638.\n\n## Vendor cross-references\n\n================================ ===================== ===================================\nSource Reported f_2 [Hz] Problem ID / location\n================================ ===================== ===================================\nClosed form (Euler-Bernoulli) 255.54 Rao 2017 \u00a78.5 Table 8.1\nTimoshenko (1974) \u00a75.3 255.54 cantilever characteristic equation\nMAPDL Verification Manual \u2248 255 VM-57 family (cantilever-shaft modal)\nAbaqus Verification Manual \u2248 255 AVM 1.6.x cantilever-bending family\nNAFEMS FV-2 255.54 cantilever transverse modes\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": [ "## Problem data + closed-form eigenvalue roots\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "E = 2.0e11 # Pa\nNU = 0.30\nRHO = 7850.0\nb = 0.05 # square cross-section side [m]\nA = b * b\nI = b**4 / 12.0 # second moment about each principal axis # noqa: E741\nL = 4.0 # span [m]; h/L = 1/80 keeps EB asymptotic to within 0.5 %\n\nBETA_L = (1.8751040687119611, 4.694091132974174, 7.854757438237613, 10.995540734875467)\n\n\ndef f_published(n: int) -> float:\n return BETA_L[n - 1] ** 2 / (2.0 * math.pi * L**2) * math.sqrt(E * I / (RHO * A))\n\n\nprint(f\"L = {L} m, b = h = {b} m, EI = {E * I:.3e} N m^2, \u03c1A = {RHO * A:.3f} kg/m\")\nfor n in (1, 2, 3, 4):\n print(\n f\"f_{n} = (\u03b2_{n} L)^2 / (2\u03c0 L^2) \u221a(EI/\u03c1A) = {f_published(n):8.3f} Hz \"\n f\"(\u03b2_{n} L = {BETA_L[n - 1]:.6f})\"\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build a 120 \u00d7 3 \u00d7 3 HEX8 EAS cantilever prism\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "NX, NY, NZ = 120, 3, 3\nxs = np.linspace(0.0, L, NX + 1)\nys = np.linspace(0.0, b, NY + 1)\nzs = np.linspace(0.0, b, 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(integration=\"enhanced_strain\"),\n material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Clamp the x = 0 face\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pts = 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, filter for transverse-dominant modes\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_modes = 24\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# For a thin cantilever (L / b = 20) the four bending pairs sit\n# **below** any torsion or axial mode (~750 Hz for first torsion,\n# ~1265 Hz for first axial). Match each published \u03b2_n to the\n# computed mode whose frequency is closest, marking the matched\n# mode so it can't be claimed twice \u2014 this is order-independent\n# and survives the eigsolver's freedom to rotate the\n# (UY-bending, UZ-bending) degenerate pair into any orthonormal\n# basis at each frequency.\n#\n# After matching, we re-confirm each pick by checking that the\n# selected mode's RMS is dominated by the transverse direction\n# (ux \u226a uy + uz) \u2014 guards against mis-matching a torsion or axial\n# mode whose frequency happens to be nearby on a non-default mesh.\n\nclaimed = np.zeros(n_modes, dtype=bool)\nmatches: list[tuple[int, int, float]] = [] # (n, k, f)\nfor n in (1, 2, 3, 4):\n f_pub_n = f_published(n)\n candidates = [k for k in range(n_modes) if not claimed[k]]\n k_best = min(candidates, key=lambda k: abs(freqs[k] - f_pub_n))\n matches.append((n, k_best, float(freqs[k_best])))\n claimed[k_best] = True\n\nprint()\nprint(\"first four transverse-bending modes (computed, matched by closest frequency)\")\nprint(f\"{'n':>3} {'k':>3} {'f_FE [Hz]':>12} {'f_pub [Hz]':>12} {'err':>9} {'tol':>5}\")\nfor n, k, f in matches:\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 transverse_pure = (uy + uz) > 2.0 * ux\n f_pub = f_published(n)\n err = (f - f_pub) / f_pub * 100.0\n tol = (0.5, 0.5, 0.5, 0.5)[n - 1] # h/L = 1/80 keeps Timoshenko shear < 0.5 %\n print(\n f\" {n} {k:>3} {f:12.3f} {f_pub:12.3f} {err:+8.3f}% \u2264{tol:.1f}% \"\n f\"{'(transverse-dominant)' if transverse_pure else '(suspected non-bending)'}\"\n )\n assert abs(err) < tol, f\"mode {n} (k={k}) deviation {err:.2f}% exceeds {tol}%\"\n assert transverse_pure, (\n f\"mode {n} (k={k}) is not transverse-dominant: ux={ux:.2f} uy={uy:.2f} uz={uz:.2f}\"\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the first four mode shapes as a 2 \u00d7 2 plotter grid\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plotter = pv.Plotter(shape=(2, 2), off_screen=True, window_size=(960, 720))\nfor n, k, f in matches:\n row, col = divmod(n - 1, 2)\n plotter.subplot(row, col)\n phi = shapes[:, :, k]\n warp = m.grid.copy()\n warp.points = m.grid.points + (b * 5.0 / np.max(np.abs(phi))) * phi\n warp[\"|phi|\"] = np.linalg.norm(phi, axis=1)\n plotter.add_mesh(m.grid, style=\"wireframe\", color=\"grey\", opacity=0.3)\n plotter.add_mesh(warp, scalars=\"|phi|\", cmap=\"plasma\", show_edges=False, show_scalar_bar=False)\n plotter.add_text(\n f\"mode {n}: f = {f:.1f} Hz (\u03b2_{n} L = {BETA_L[n - 1]:.3f})\",\n position=\"upper_edge\",\n font_size=10,\n )\nplotter.link_views()\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 }