{ "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# Simply-supported plate \u2014 first transverse bending mode\n\nClassical thin-plate eigenvalue problem. A simply-supported\nsquare plate of side $a$ and thickness $h$ has\nKirchhoff natural frequencies (Timoshenko & Woinowsky-Krieger\n1959 \u00a763 eq. 151)\n\n\\begin{align}f_{mn} = \\frac{\\pi}{2}\n \\left(\\frac{m^{2}}{a^{2}} + \\frac{n^{2}}{b^{2}}\\right)\n \\sqrt{\\frac{D}{\\rho\\, h}},\n \\qquad\n D = \\frac{E\\, h^{3}}{12 (1 - \\nu^{2})},\\end{align}\n\nso the (1, 1) fundamental of a square plate $a = b$ is\n\n\\begin{align}f_{11} = \\pi \\frac{1}{a^{2}}\n \\sqrt{\\frac{D}{\\rho\\, h}}.\\end{align}\n\n## Implementation\n\nA 30 \u00d7 30 \u00d7 2 HEX8 enhanced-strain (Simo\u2013Rifai 1990) mesh on\nthe unit square plate at $L/h = 100$. Simply-supported\nboundary conditions pin $u_z = 0$ along all four edges\nthrough the full thickness; in-plane rigid-body translation is\nremoved at one corner. The mode-shape filter selects the first\nmode whose UZ root-mean-square dominates the in-plane components\n(i.e. the (1, 1) bending mode), in case any spurious in-plane\nmodes precede it on a coarse mesh.\n\nModal analysis runs through :meth:`Model.solve_modal` (ARPACK\nshift-invert via :mod:`femorph_solver.solvers.modal`).\n\n## References\n\n* Timoshenko, S. P. and Woinowsky-Krieger, S. (1959) *Theory\n of Plates and Shells*, 2nd ed., McGraw-Hill, \u00a763 eq. 151\n (rectangular plate eigenfrequencies).\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, \u00a712.5.\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_1 [Hz] Problem ID / location\n====================================== ===================== ============================================\nClosed form (Kirchhoff) ~278 Timoshenko & Woinowsky-Krieger \u00a763\nLeissa NASA SP-160 ~278 Table 4.3 (SSSS square plate)\nNAFEMS BENCHMARK-SSPlate modal ~278 NAFEMS FV52 (simply-supported thin square plate)\nMAPDL Verification Manual ~278 VM-62 (vibration of a thin plate)\nAbaqus Verification Manual ~278 AVM 1.4.6 (simply-supported plate \u2014 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\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "E = 2.0e11 # Pa (steel)\nNU = 0.3\nRHO = 7850.0 # kg/m^3\na = 1.0 # plate edge length [m]\nh = 0.01 # thickness [m]; L/h = 100\n\nNX = NY = 30\nNZ = 2\n\nD = E * h**3 / (12.0 * (1.0 - NU**2))\nf11_published = (math.pi / 2.0) * (1.0 / a**2 + 1.0 / a**2) * math.sqrt(D / (RHO * h))\nprint(f\"a = {a} m, h = {h} m, L/h = {a / h:.0f}\")\nprint(f\"E = {E / 1e9:.0f} GPa, \u03bd = {NU}, \u03c1 = {RHO} kg/m^3, D = {D:.3e} N m\")\nprint(f\"f_11 published (Timoshenko \u00a763) = {f11_published:.4f} Hz\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build a 30 \u00d7 30 \u00d7 2 HEX8 mesh\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "xs = np.linspace(0.0, a, NX + 1)\nys = np.linspace(0.0, a, NY + 1)\nzs = np.linspace(0.0, h, NZ + 1)\ngrid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\nprint(f\"mesh: {grid.n_points} nodes, {grid.n_cells} HEX8 cells\")\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": [ "## SS BCs + rigid-body suppression\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pts = np.asarray(m.grid.points)\ntol = 1e-9\nedge = (\n (np.abs(pts[:, 0]) < tol)\n | (np.abs(pts[:, 0] - a) < tol)\n | (np.abs(pts[:, 1]) < tol)\n | (np.abs(pts[:, 1] - a) < tol)\n)\nm.fix(nodes=(np.where(edge)[0] + 1).tolist(), dof=\"UZ\")\ncorner = int(np.where(np.linalg.norm(pts, axis=1) < tol)[0][0])\nm.fix(nodes=[corner + 1], dof=\"UX\")\nm.fix(nodes=[corner + 1], dof=\"UY\")\ncorner_x = int(\n np.where((np.abs(pts[:, 0] - a) < tol) & (np.abs(pts[:, 1]) < tol) & (np.abs(pts[:, 2]) < tol))[\n 0\n ][0]\n)\nm.fix(nodes=[corner_x + 1], dof=\"UY\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Modal solve + filter for (1, 1) bending mode\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_modes = 10\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\nf11_idx = None\nfor k in range(n_modes):\n ux_rms = float(np.sqrt((shapes[:, 0, k] ** 2).mean()))\n uy_rms = float(np.sqrt((shapes[:, 1, k] ** 2).mean()))\n uz_rms = float(np.sqrt((shapes[:, 2, k] ** 2).mean()))\n if uz_rms > 3.0 * max(ux_rms, uy_rms):\n f11_idx = k\n break\n\nassert f11_idx is not None, f\"no UZ-dominant mode found in first {n_modes}\"\n\nf11_computed = float(freqs[f11_idx])\nerr_pct = (f11_computed - f11_published) / f11_published * 100.0\nprint()\nprint(f\"first 5 frequencies (Hz) = {[round(float(x), 3) for x in freqs[:5]]}\")\nprint(f\"selected (1,1) bending mode idx = {f11_idx}\")\nprint(f\"f_11 computed = {f11_computed:.4f} Hz\")\nprint(f\"f_11 published = {f11_published:.4f} Hz\")\nprint(f\"relative error = {err_pct:+.2f} %\")\n\n# 0.5 % tolerance \u2014 HEX8 enhanced-strain on the default 30 \u00d7 30 \u00d7 2\n# mesh recovers the Kirchhoff fundamental to within ~0.06 % at the\n# default thin geometry (L/h = 100). Tighter convergence (sub-0.1 %)\n# would need a dedicated SHELL kernel (tracked separately).\nassert abs(err_pct) < 0.5, f\"f_11 deviation {err_pct:.2f}% exceeds 0.5%\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the (1, 1) bending mode shape\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "phi = shapes[:, :, f11_idx]\nwarped = m.grid.copy()\nwarped.points = m.grid.points + (h / np.max(np.abs(phi)) * 0.4) * phi\nwarped[\"UZ\"] = phi[:, 2]\n\nplotter = pv.Plotter(off_screen=True, window_size=(720, 480))\nplotter.add_mesh(\n m.grid,\n style=\"wireframe\",\n color=\"grey\",\n opacity=0.35,\n label=\"undeformed\",\n)\nplotter.add_mesh(\n warped,\n scalars=\"UZ\",\n cmap=\"coolwarm\",\n show_edges=False,\n scalar_bar_args={\"title\": f\"mode (1,1) \u2014 f = {f11_computed:.1f} Hz\"},\n)\nplotter.view_isometric()\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 }