{ "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# Tutorial 2 \u2014 Modal survey of a cantilever bracket\n\nPicking up from :doc:`tutorial_01_cantilever_combined`, the\nstatic analysis confirmed the bracket carries its design load\nwith a comfortable safety margin. But the equipment console it\nsupports has variable-speed motors that excite the bracket at\nmultiple harmonics across the operating band. The next\nquestion is **dynamic**: do any of the bracket's natural\nfrequencies fall inside the excitation band, and if so which\nmodes need to be included in a response-spectrum analysis to\ncapture 90 % of the modal mass?\n\nThis tutorial walks the modal-survey workflow end-to-end:\n\n* **Step 1** \u2014 build the same cantilever geometry from\n Tutorial 1, this time as a 3-D solid that carries its own\n mass (no separate tip mass, just the bracket's structural\n mass distributed by ``\u03c1 \u00d7 V``).\n* **Step 2** \u2014 run a modal solve and inspect the spectrum.\n* **Step 3** \u2014 identify each mode's character (transverse\n bending, axial, torsional) by inspecting the recovered\n eigenvectors.\n* **Step 4** \u2014 compute modal participation factors and\n effective modal mass per excitation direction (X, Y, Z).\n* **Step 5** \u2014 read the cumulative-effective-mass curve to\n pick $n_\\text{modes}$ for a response-spectrum analysis\n meeting the ASCE 7 / Eurocode 8 90 % threshold.\n* **Step 6** \u2014 render the lowest few mode shapes for the\n design-review packet.\n\nThe post-processing recipes for participation factors and\nmode-shape animation live in\n`sphx_glr_gallery_post-processing_example_modal_participation.py`\nand\n`sphx_glr_gallery_post-processing_example_mode_shape_animation.py`\n\u2014 they exercise the same APIs in narrower scope.\n\n## Theoretical reference\n\nA free-vibration modal analysis solves the generalised\nsymmetric-definite eigenvalue problem\n(:doc:`/reference/theory/eigenproblem`)\n\n\\begin{align}\\mathbf{K}\\, \\boldsymbol{\\phi}_i = \\omega_i^{2}\\,\n \\mathbf{M}\\, \\boldsymbol{\\phi}_i,\\end{align}\n\nwith $\\mathbf{K}$ and $\\mathbf{M}$ the assembled\nstiffness and mass matrices. Each eigenpair returns one mode\nshape $\\boldsymbol{\\phi}_i$ and one natural frequency\n$f_i = \\omega_i / (2\\pi)$.\n\nFor base-excitation analysis (the bracket bolted to a vibrating\nfloor), the **modal participation factor** for mode $i$\nin direction $d$ is\n\n\\begin{align}\\Gamma_i^{(d)} = \\frac{\\boldsymbol{\\phi}_i^\\top \\mathbf{M}\\, \\mathbf{r}^{(d)}}\n {\\boldsymbol{\\phi}_i^\\top \\mathbf{M}\\, \\boldsymbol{\\phi}_i},\\end{align}\n\nwhere $\\mathbf{r}^{(d)}$ is the **influence vector** \u2014\nones on every translational DOF along direction $d$,\nzeros elsewhere. The corresponding **effective modal mass**\nis\n\n\\begin{align}M_\\text{eff,i}^{(d)} = (\\Gamma_i^{(d)})^{2} \\cdot\n \\boldsymbol{\\phi}_i^\\top \\mathbf{M}\\, \\boldsymbol{\\phi}_i,\\end{align}\n\nand the partition-of-unity property\n$\\sum_i M_\\text{eff,i}^{(d)} = M_\\text{total}$ lets the\n*cumulative* effective-mass curve answer the design-code\nquestion \"how many modes do I need?\".\n\nReferences: Bathe 2014 \u00a79.3.4, Chopra 2017 \u00a713.1, Cook et al.\n2002 \u00a711.6.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport numpy as np\nimport pyvista as pv\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1 \u2014 build the bracket\n\nSame prismatic-cantilever geometry and material as Tutorial 1,\nthis time at slightly higher slenderness to push the first\ntransverse-bending frequency lower. We keep HEX8-EAS so the\nbending modes don't shear-lock.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "E = 2.0e11 # Pa, steel\nNU = 0.30\nRHO = 7850.0 # kg/m\u00b3\nL = 1.0 # m\nWIDTH = 0.05\nHEIGHT = 0.05\n\nNX, NY, NZ = 60, 3, 3\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(integration=\"enhanced_strain\"),\n material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)\n\npts = np.asarray(m.grid.points)\nnode_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"], dtype=np.int64)\nclamped = np.where(pts[:, 0] < 1e-9)[0]\nm.fix(nodes=node_nums[clamped].tolist(), dof=\"ALL\")\n\ntotal_mass = RHO * L * WIDTH * HEIGHT\nprint(\"Cantilever bracket \u2014 modal survey\")\nprint(f\" L = {L} m, cross = {WIDTH} \u00d7 {HEIGHT} m, \u03c1 = {RHO} kg/m\u00b3\")\nprint(f\" Mesh: {NX} \u00d7 {NY} \u00d7 {NZ} HEX8-EAS = {m.grid.n_cells} cells, {m.grid.n_points} nodes\")\nprint(f\" Total structural mass M_total = \u03c1\u00b7L\u00b7A = {total_mass:.3f} kg\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2 \u2014 modal solve and spectrum inspection\n\nTwelve modes is enough for a slender cantilever to clear the\n90 % cumulative-mass threshold in the bending directions \u2014\nmore than we need but it gives a clean spectrum to study.\n\n:meth:`Model.solve_modal `\ndispatches to :func:`scipy.sparse.linalg.eigsh` via shift-\ninvert Lanczos at $\\sigma = 0$\n(:doc:`/reference/theory/eigenproblem`). ``n_modes`` controls\nhow many lowest eigenpairs come back.\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, dtype=np.float64)\n# ``mode_shapes`` is shaped ``(n_dof, n_modes)``. Each column is a\n# mass-orthonormalised eigenvector.\n\nprint()\nprint(f\" {'mode':>4} {'f [Hz]':>10}\")\nprint(\" \" + \"-\" * 18)\nfor i in range(N_MODES):\n print(f\" {i + 1:>4} {freqs[i]:>10.3f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3 \u2014 mode-shape character via translation-energy split\n\nEach mode mixes translational components. A *transverse-\nbending* mode has nearly all its kinetic energy in $u_y$\n(or $u_z$); an *axial* mode has it in $u_x$; a\n*torsional* mode has the warping pattern characteristic of\nSaint-Venant torsion.\n\nWe classify each mode by the per-direction translation-energy\nfraction. The square cross-section makes $y$ / $z$\nbending pairs degenerate, so transverse modes appear in\npairs at the same frequency.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dof_map = m.dof_map()\nclassifications = []\nfor i in range(N_MODES):\n phi = shapes[:, i]\n energy_total = (phi**2).sum()\n energy_per_dof = np.zeros(3)\n for row, (_node, dof_idx) in enumerate(dof_map.tolist()):\n if dof_idx < 3:\n energy_per_dof[int(dof_idx)] += phi[row] ** 2\n fractions = energy_per_dof / (energy_total + 1e-30)\n # Pick the dominant direction.\n label = (\"axial (x)\", \"transverse-y\", \"transverse-z\")[int(np.argmax(fractions))]\n classifications.append((label, fractions))\n\nprint()\nprint(f\" {'mode':>4} {'f [Hz]':>9} {'character':>14} {'%U_x':>6} {'%U_y':>6} {'%U_z':>6}\")\nprint(\" \" + \"-\" * 53)\nfor i, (label, frac) in enumerate(classifications):\n print(\n f\" {i + 1:>4} {freqs[i]:>9.3f} {label:>14} \"\n f\"{100 * frac[0]:>5.1f}% {100 * frac[1]:>5.1f}% {100 * frac[2]:>5.1f}%\"\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 4 \u2014 modal participation factors and effective mass\n\nBuild the X / Y / Z influence vectors once, then compute\n\u0393 and M_eff for every mode \u00d7 every direction in one\nvectorised pass. The recipe matches\n`sphx_glr_gallery_post-processing_example_modal_participation.py`.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M = m.mass_matrix()\nn_dof = M.shape[0]\n\nr_x = np.zeros(n_dof)\nr_y = np.zeros(n_dof)\nr_z = np.zeros(n_dof)\nfor row, (_node, dof_idx) in enumerate(dof_map.tolist()):\n if dof_idx == 0:\n r_x[row] = 1.0\n elif dof_idx == 1:\n r_y[row] = 1.0\n elif dof_idx == 2:\n r_z[row] = 1.0\n\ndirections = ((\"X\", r_x), (\"Y\", r_y), (\"Z\", r_z))\ngamma = np.zeros((N_MODES, 3))\nm_eff = np.zeros((N_MODES, 3))\nfor i in range(N_MODES):\n phi = shapes[:, i]\n M_phi = M @ phi\n modal_mass_i = float(phi @ M_phi) # \u2248 1 for mass-orthonormal\n for d_idx, (_, r) in enumerate(directions):\n L_i = float(phi @ M @ r)\n gamma_i = L_i / modal_mass_i\n gamma[i, d_idx] = gamma_i\n m_eff[i, d_idx] = gamma_i**2 * modal_mass_i\n\ncum_eff = np.cumsum(m_eff, axis=0)\ncum_pct = cum_eff / total_mass * 100.0\n\nprint()\nprint(\" Cumulative effective mass [%] per direction:\")\nprint(f\" {'mode':>4} {'f [Hz]':>9} {'cum_X':>7} {'cum_Y':>7} {'cum_Z':>7}\")\nprint(\" \" + \"-\" * 48)\nfor i in range(N_MODES):\n print(\n f\" {i + 1:>4} {freqs[i]:>9.3f} \"\n f\"{cum_pct[i, 0]:>6.2f}% {cum_pct[i, 1]:>6.2f}% {cum_pct[i, 2]:>6.2f}%\"\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 5 \u2014 design-code mode-count threshold\n\nASCE 7 \u00a712.9 and Eurocode 8 \u00a74.3.3.3 both require **at least\n90 % cumulative effective mass** in every direction the\nresponse-spectrum input acts. Read the thresholds off the\ncumulative table.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "target_pct = 90.0\nfor d_idx, dname in enumerate((\"X\", \"Y\", \"Z\")):\n crossings = np.where(cum_pct[:, d_idx] >= target_pct)[0]\n if len(crossings) > 0:\n n_required = int(crossings[0]) + 1\n print(\n f\" {dname}: needs \u2265 {n_required} modes for {target_pct:.0f} % cumulative effective mass \"\n f\"({cum_pct[n_required - 1, d_idx]:.2f} % at mode {n_required})\"\n )\n else:\n print(\n f\" {dname}: {N_MODES} modes only reach {cum_pct[-1, d_idx]:.2f} % \u2014 \"\n f\"needs more modes to clear {target_pct:.0f} % threshold\"\n )\n\n# X direction has no bending modes in the first 12 (all bending\n# is transverse y / z), so X stays at 0 % cumulative mass \u2014 that's\n# the correct physics for a slender cantilever where axial modes\n# only appear far above the bending fundamentals. Asserts only\n# Y / Z which are the design-relevant directions.\ny_n = (\n int(np.where(cum_pct[:, 1] >= target_pct)[0][0]) + 1\n if (cum_pct[:, 1] >= target_pct).any()\n else N_MODES\n)\nz_n = (\n int(np.where(cum_pct[:, 2] >= target_pct)[0][0]) + 1\n if (cum_pct[:, 2] >= target_pct).any()\n else N_MODES\n)\nn_for_response_spectrum = max(y_n, z_n)\nprint()\nprint(\n f\" \u2192 Use n_modes = {n_for_response_spectrum} for a 3-direction \"\n f\"response-spectrum analysis (max over Y / Z; X is axial-only and \"\n f\"out of the bending band the equipment excites).\"\n)\nassert n_for_response_spectrum <= N_MODES, \"default 12 modes was too few \u2014 increase the solve\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 6 \u2014 render the lowest four mode shapes\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mode_pairs = []\nseen_freqs: list[float] = []\nfor i in range(N_MODES):\n f = freqs[i]\n if any(abs(f - sf) < 1e-3 for sf in seen_freqs):\n continue\n seen_freqs.append(float(f))\n mode_pairs.append(i)\n if len(mode_pairs) >= 4:\n break\n\nplotter = pv.Plotter(shape=(1, 4), off_screen=True, window_size=(1400, 320), border=False)\nfor col, idx in enumerate(mode_pairs):\n plotter.subplot(0, col)\n grid_warped = m.grid.copy()\n n_pts = grid_warped.n_points\n disp = np.zeros((n_pts, 3))\n for row, (node, dof_idx) in enumerate(dof_map.tolist()):\n if dof_idx < 3:\n pt_idx = np.where(node_nums == int(node))[0]\n if len(pt_idx):\n disp[pt_idx[0], int(dof_idx)] += shapes[row, idx]\n peak = float(np.linalg.norm(disp, axis=1).max())\n scale = 0.10 * L / peak if peak > 0 else 1.0\n grid_warped.points = np.asarray(m.grid.points) + scale * disp\n grid_warped[\"|disp|\"] = np.linalg.norm(disp, axis=1)\n plotter.add_mesh(m.grid, style=\"wireframe\", color=\"grey\", opacity=0.3)\n plotter.add_mesh(\n grid_warped,\n scalars=\"|disp|\",\n cmap=\"viridis\",\n show_edges=False,\n show_scalar_bar=False,\n )\n plotter.add_text(\n f\"mode {idx + 1}\\nf = {freqs[idx]:.2f} Hz\\n{classifications[idx][0]}\",\n position=\"upper_left\",\n font_size=9,\n )\n plotter.view_isometric()\n plotter.camera.zoom(1.05)\nplotter.link_views()\nplotter.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Take-aways\n\n- **The modal solve** is a single ``Model.solve_modal(n_modes=N)``\n call; the eigensolver picks shift-invert Lanczos automatically\n (:doc:`/reference/solvers/modal`).\n- **Mode classification** by translation-energy fraction is a\n robust filter against degenerate transverse pairs and stray\n torsional modes \u2014 works on any geometry.\n- **Modal participation + effective mass** is computed from\n first principles via ``Model.mass_matrix()`` plus the\n per-direction influence vector built from ``Model.dof_map()``.\n- **Mode-count selection** for response-spectrum analysis\n reads directly off the cumulative-effective-mass table:\n ASCE 7 / Eurocode 8 require \u2265 90 % per direction.\n- **Mode-shape rendering** uses the same scatter pattern as\n `sphx_glr_gallery_post-processing_example_mode_shape_plot.py`\n \u2014 display amplitude is scaled to a fraction of the bounding\n box (mode shapes carry no absolute magnitude).\n\nNext steps:\n\n- Tutorial 3 (planned) \u2014 pressure-vessel design-by-analysis\n on a Lam\u00e9 thick cylinder.\n- For a worked **response-spectrum** computation that builds\n on this modal survey, see the planned Tutorial 4 \u2014 it\n convolves a synthetic ASCE 7 design spectrum with the\n participation-factor table generated above.\n\n" ] } ], "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 }