{
  "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 4 \u2014 Harmonic frequency response of a cantilever bracket\n\nPicking up from :doc:`tutorial_02_modal_survey`, we now know\nthe bracket's natural frequencies and which modes carry the\neffective mass in each excitation direction.  The motors driving\nthe equipment console aren't a single sine, though \u2014 they sweep\nthrough their operating band as the production line ramps up\nand slows down, exciting the bracket at every frequency the\nsweep crosses.\n\nThe next question is **steady-state harmonic response**: at each\nforcing frequency $\\omega$, what is the phase-locked\ndisplacement response $u(\\omega)$ at the bracket tip?  The\nanswer drives the **frequency-response function (FRF)** \u2014\npeak amplitudes that determine fatigue life, and phase shifts\nthat determine which combinations of forcing harmonics\nconstructively / destructively interfere.\n\nThis tutorial walks the harmonic-response workflow end-to-end:\n\n* **Step 1** \u2014 rebuild the cantilever bracket (same geometry\n  and material as Tutorial 1, smaller mesh for fast solves).\n* **Step 2** \u2014 recover the lowest few modes and decide how\n  many to retain in the modal-superposition expansion.\n* **Step 3** \u2014 apply a unit tip-force excitation and sweep a\n  band of forcing frequencies that brackets the first three\n  natural frequencies.\n* **Step 4** \u2014 compare three damping models \u2014 undamped,\n  light modal damping ($\\zeta = 0.005$), and heavier\n  Rayleigh damping \u2014 to show how each shapes the FRF.\n* **Step 5** \u2014 identify the resonant peaks, read off the\n  amplification factor, and check it against the\n  half-power-bandwidth estimate.\n* **Step 6** \u2014 render the FRF as a publication-grade\n  log-magnitude plot.\n\n## Theoretical reference\n\nFor a linear undamped FE system, the equation of motion is\n\n\\begin{align}\\mathbf{M}\\, \\ddot{\\mathbf{u}}(t) + \\mathbf{K}\\, \\mathbf{u}(t)\n        = \\mathbf{f}(t).\\end{align}\n\nLooking for steady-state harmonic response under a single-\nfrequency excitation\n$\\mathbf{f}(t) = \\mathbf{F}\\, e^{i \\omega t}$ substitutes\n$\\mathbf{u}(t) = \\mathbf{U}(\\omega)\\, e^{i \\omega t}$ and\nturns the time-domain ODE into the frequency-domain algebraic\nsystem\n\n\\begin{align}\\bigl( \\mathbf{K} - \\omega^2\\, \\mathbf{M} + i \\omega\\, \\mathbf{C} \\bigr)\\,\n    \\mathbf{U}(\\omega) = \\mathbf{F}.\\end{align}\n\nSolving that directly at every $\\omega$ is the **direct\nfrequency response** approach \u2014 accurate but expensive (one\nsparse factorisation per frequency).\n\n**Modal superposition** projects the response onto the lowest\n$n_\\text{modes}$ mass-orthonormal eigenmodes:\n\n\\begin{align}\\mathbf{U}(\\omega) \\approx \\sum_{i=1}^{n_\\text{modes}}\n    \\frac{\\boldsymbol{\\phi}_i^\\top \\mathbf{F}}\n         {\\omega_i^2 - \\omega^2 + 2 i \\zeta_i\\, \\omega_i\\, \\omega}\\,\n    \\boldsymbol{\\phi}_i,\\end{align}\n\nwith per-mode damping ratio $\\zeta_i$.  Each mode is a\nsingle SDOF resonator; superposing them captures the full FRF\nfor one fixed-cost factorisation upfront.  See\n:doc:`/reference/theory/eigenproblem` for the math, and\n:doc:`/user-guide/solving/harmonic` for the API choices.\n\nCompanion focused recipe in\n`sphx_glr_gallery_analyses_harmonic_example_sdof_frf.py` \u2014\nsame workflow on a single-degree-of-freedom toy.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport matplotlib\n\nmatplotlib.use(\"Agg\")  # headless: no GUI window pop in gallery / CI\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pyvista as pv\nfrom scipy.signal import find_peaks\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 1 \u2014 geometry, material, mesh\n\nSame prismatic steel cantilever as Tutorial 1 (one-metre long,\n50 \u00d7 50 mm square cross-section).  A modest 24 \u00d7 3 \u00d7 3 mesh\nresolves the lowest five bending and axial modes to better\nthan 0.1 % against the Euler-Bernoulli closed forms \u2014 and\nkeeps the harmonic-sweep cost down to a few seconds total\nacross the whole frequency band.\n\nHEX8 with enhanced-assumed-strain integration is the natural\nfit for slender solids in bending (see\n:doc:`/reference/elements/hex8` \u2014 plain HEX8 shear-locks).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.0e11  # Young's modulus, Pa (steel)\nNU = 0.30  # Poisson's ratio\nRHO = 7850.0  # density, kg/m\u00b3 (steel)\nL = 1.0  # length, m\nW = 0.05  # cross-section width, m\nH = 0.05  # cross-section height, m\n\nprint(\"Cantilever bracket \u2014 harmonic frequency response\")\nprint(f\"  L = {L} m, cross-section = {W} \u00d7 {H} m\")\nprint(f\"  E = {E / 1e9:.0f} GPa, \u03bd = {NU}, \u03c1 = {RHO} kg/m\u00b3\")\n\nxs = np.linspace(0.0, L, 25)\nys = np.linspace(0.0, W, 4)\nzs = np.linspace(0.0, H, 4)\ngrid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\n\nmodel = femorph_solver.Model.from_grid(grid)\nmodel.assign(\n    ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n    material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)\n\npts = np.asarray(model.grid.points)\nnode_nums = np.asarray(model.grid.point_data[\"ansys_node_num\"], dtype=np.int64)\nclamped = np.where(pts[:, 0] < 1e-9)[0]\nmodel.fix(nodes=node_nums[clamped].tolist(), dof=\"ALL\")\nprint(f\"  Mesh: {model.n_nodes} nodes, {model.n_elements} HEX8-EAS cells\")\nprint(f\"  Clamped {len(clamped)} root-face nodes (full-fix).\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 2 \u2014 modal extraction\n\nModal superposition needs the lowest few mass-orthonormal\neigenpairs.  How many?  Convention: retain enough modes so the\nhighest one's natural frequency is above the top of the\nexcitation band.  A safety-multiplier of 1.5-2\u00d7 is the\nrule-of-thumb so the truncated tail of the modal expansion\ndoesn't bias the frequency-domain peak that's right below it.\n\nThe bracket from Tutorial 1 has its first transverse-bending\nmode near 40 Hz and its second-bending mode in the same plane\nnear 700 Hz.  We'll sweep the FRF up to 1000 Hz so both peaks\nfall inside the band; retaining 12 modes covers the band with\nenough margin that the truncation tail doesn't bias the\nresponse near 1000 Hz.\n\nThe intermediate $\\sim 250$ Hz peaks belong to the\n*orthogonal* bending plane (UY), which a UZ-only tip force\ndoes not excite \u2014 they appear in the spectrum but not in the\nUZ tip displacement.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_modes_retained = 12\nmodal = model.solve_modal(n_modes=n_modes_retained)\nprint(f\"\\nFirst {n_modes_retained} natural frequencies (Hz):\")\nfor i, f in enumerate(modal.frequency, start=1):\n    print(f\"  mode {i:2d}:  {f:8.3f} Hz\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 3 \u2014 excitation, frequency band, frequency-response solve\n\nThe bracket's tip is loaded by a transverse harmonic force\n$F_z(t) = F_0 \\cos(\\omega t)$.  Magnitude is irrelevant\nin linear analysis \u2014 the FRF is normalised by the input force\namplitude and reported as receptance\n$H(\\omega) = U_\\text{tip}(\\omega) / F_0$ \u2014 but a unit\nload (1 N at every tip node, distributed) is the convention\nboth NASTRAN and MAPDL ship.\n\nThe frequency band is chosen to bracket the first three\nbending modes plus a healthy quiet zone above the third peak.\n300 sample points with logarithmic spacing in the\nresonance-rich low-band gives clean peak resolution without\nblowing the wall time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tip_mask = pts[:, 0] > L - 1e-9\ntip_nodes = node_nums[tip_mask]\nfor n in tip_nodes:\n    model.apply_force(int(n), fz=1.0 / int(tip_mask.sum()))\nprint(f\"  Tip excitation: 1 N total across {tip_mask.sum()} nodes (z-direction).\")\n\n# Sample density must resolve the resonant peaks \u2014 the half-\n# power bandwidth at :math:`\\zeta = 0.02` damping is\n# :math:`\\Delta f \\approx 2 \\zeta f_n`, so a 40-Hz peak is only\n# 1.6 Hz wide.  We use a base linspace of 600 points across the\n# 5-600 Hz band (\u2248 1 Hz spacing) and concatenate finer samples\n# clustered around each predicted natural frequency so the peak\n# itself can't fall between two sample points.\nF_MAX = 1000.0  # top of the sweep band, Hz\nbase = np.linspace(5.0, F_MAX, 800)\nmode_band: list[np.ndarray] = []\nfor f_n in modal.frequency:\n    if 5.0 <= f_n <= F_MAX:\n        # \u00b13 Hz around each mode: enough to resolve the\n        # half-power bandwidth of every peak the field will see.\n        mode_band.append(np.linspace(f_n - 3.0, f_n + 3.0, 61))\nfrequencies_hz = np.unique(np.concatenate([base, *mode_band]))\nprint(\n    f\"  Sweeping {frequencies_hz.size} frequencies from \"\n    f\"{frequencies_hz[0]:.0f} Hz to {frequencies_hz[-1]:.0f} Hz \"\n    f\"(refined around the {len(mode_band)} in-band modes).\"\n)\n\n# Three damping flavours, all using modal-superposition\n# (`Model.solve_harmonic` with `n_modes=...`).  Reusing the\n# already-computed modal result via the `modal_result=` kwarg\n# saves a second eigenproblem.\nprint(\"\\nRunning three modal-superposition harmonic sweeps:\")\n\nprint(\"  - lightly-damped (zeta = 0.005, near-undamped reference)\")\nresult_light = model.solve_harmonic(\n    frequencies=frequencies_hz,\n    modal_damping_ratio=0.005,\n    modal_result=modal,\n)\n\nprint(\"  - realistic modal damping (zeta = 0.02, welded-steel field measurement)\")\nresult_steel = model.solve_harmonic(\n    frequencies=frequencies_hz,\n    modal_damping_ratio=0.02,\n    modal_result=modal,\n)\n\nprint(\"  - Rayleigh damping (alpha = 0.5, beta = 1e-5)\")\nresult_rayleigh = model.solve_harmonic(\n    frequencies=frequencies_hz,\n    rayleigh=(0.5, 1.0e-5),\n    modal_result=modal,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 4 \u2014 extracting the tip displacement\n\n`HarmonicResult.displacement` is a complex\n``(n_freq, N)`` array \u2014 the steady-state displacement\nmagnitude and phase at every DOF, at every forcing frequency.\nWe need the tip-node UZ component to plot the receptance.\n\n`Model.dof_map` gives the per-row ``(node_num, local_dof_idx)``\nmapping; we look up the tip node's UZ DOF (local index 2) and\nslice that column out of the response array.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dof_map = model.dof_map()\n# Pick a representative tip node \u2014 the one closest to the\n# section centroid \u2014 so the FRF plot reflects the bulk tip\n# motion, not a corner whose response carries some torsion.\ntip_centre_node = int(\n    node_nums[np.argmin((pts[:, 0] - L) ** 2 + (pts[:, 1] - W / 2) ** 2 + (pts[:, 2] - H / 2) ** 2)]\n)\ntip_uz_row = int(np.where((dof_map[:, 0] == tip_centre_node) & (dof_map[:, 1] == 2))[0][0])\nprint(f\"\\n  Tip-centre node = {tip_centre_node}, UZ row = {tip_uz_row}\")\n\n\ndef tip_amplitude(result) -> np.ndarray:\n    \"\"\"Magnitude of the complex tip-UZ displacement, one per frequency.\"\"\"\n    return np.abs(result.displacement[:, tip_uz_row])\n\n\nu_light = tip_amplitude(result_light)\nu_steel = tip_amplitude(result_steel)\nu_rayleigh = tip_amplitude(result_rayleigh)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 5 \u2014 peak identification and amplification factor\n\nThe lightly-damped FRF carries clean, isolated peaks at the\nbending-mode natural frequencies; we find them by looking for\nlocal maxima above a sensible noise threshold.\n\nAt each peak, the **amplification factor** is the ratio of\nthe peak displacement to the static (zero-frequency) baseline:\n\n\\begin{align}A_i = \\frac{|U(\\omega_i)|}{|U(0)|}.\\end{align}\n\nFor an isolated SDOF resonator, the amplification is\n$1 / (2 \\zeta_i)$ \u2014 so a 2 % damping ratio (realistic\nfor welded-steel field measurements) predicts a peak\namplification of 25\u00d7, which we'll verify against the\nrealistic-damping curve below.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "u_static = u_steel[0]  # zero-frequency-ish baseline\nprint(f\"\\n  Static (5 Hz) tip displacement: {u_static:.3e} m  (per 1 N tip force)\")\n\n# scipy.signal.find_peaks handles plateau peaks (which the\n# refined sampling around the natural frequency can create) and\n# accepts a `prominence` filter so off-resonance ripples don't\n# count as peaks.\n# Threshold chosen so only the genuinely-excited UZ-bending\n# peaks survive \u2014 orthogonal-plane modes (UY bending, mode 3/4\n# at 255 Hz) excite the tip-centre UZ row only weakly and would\n# otherwise clutter the printout.\npeak_idx, _ = find_peaks(u_steel, prominence=u_static)\nprint(\"\\n  Resonant peaks found in the UZ tip-FRF (prominence > 1\u00d7 static):\")\nprint(\n    f\"    {'#':>3}   {'f_peak (Hz)':>12}   {'|U_peak| (m)':>14}   {'A = |U_peak|/|U_static|':>26}\"\n)\nfor k, idx in enumerate(peak_idx[:5], start=1):\n    f_peak = frequencies_hz[idx]\n    u_peak = u_steel[idx]\n    amp = u_peak / u_static\n    print(f\"    {k:>3}   {f_peak:>12.2f}   {u_peak:>14.3e}   {amp:>26.1f}\")\n\nprint(f\"\\n  Theoretical SDOF amplification at zeta = 0.02:  {1.0 / (2.0 * 0.02):.0f}\u00d7\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 6 \u2014 plot the FRF (magnitude vs frequency)\n\nThree curves on one log-magnitude plot \u2014 the three damping\nflavours from Step 3.  We use a log-y axis because the\nresonant peaks span four orders of magnitude on the\nundamped curve and a linear axis would compress the off-\nresonance baseline into a single line at the bottom of the\nplot.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(8.5, 5.0))\nax.semilogy(\n    frequencies_hz,\n    u_light,\n    color=\"C3\",\n    lw=1.6,\n    label=r\"light damping, $\\zeta = 0.005$  (near-undamped)\",\n)\nax.semilogy(\n    frequencies_hz,\n    u_steel,\n    color=\"C0\",\n    lw=1.6,\n    label=r\"realistic damping, $\\zeta = 0.02$  (welded steel)\",\n)\nax.semilogy(\n    frequencies_hz,\n    u_rayleigh,\n    color=\"C2\",\n    lw=1.6,\n    label=r\"Rayleigh damping, $\\alpha=0.5$, $\\beta=10^{-5}$\",\n)\n\n# Annotate the first few peaks on the realistic-damping curve.\nfor idx in peak_idx[:3]:\n    f_peak = frequencies_hz[idx]\n    u_peak = u_steel[idx]\n    ax.axvline(f_peak, color=\"0.7\", lw=0.8, ls=\"--\", zorder=0)\n    ax.text(f_peak, u_peak * 1.4, f\"{f_peak:.0f} Hz\", ha=\"center\", fontsize=9, color=\"0.25\")\n\nax.set_xlabel(\"Forcing frequency (Hz)\")\nax.set_ylabel(r\"Tip-displacement magnitude $|U_z(\\omega)|$ (m / N)\")\nax.set_title(\"Cantilever bracket \u2014 harmonic frequency response\")\nax.grid(True, which=\"both\", ls=\":\", alpha=0.5)\nax.legend(loc=\"upper right\", fontsize=9)\nfig.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Engineering takeaway\n\nThree things to read off the FRF before we sign the design:\n\n1. **Peak frequency separation.**  The first bending peak sits\n   near $f_1 \\approx 40$ Hz, *just below* the\n   $50\\,\\text{Hz}$ line-frequency excitation that\n   dominates the equipment console's electromagnetic forcing \u2014\n   a 25 % margin on a band-of-interest that's normally\n   \u00b110 % broad.  The fundamental dominates the tip UZ\n   response by orders of magnitude (canonical for a cantilever\n   excited near its tip): the second-bending mode at\n   $f_3 \\approx 707$ Hz contributes a barely-visible\n   secondary peak, and the orthogonal-plane modes at\n   $\\approx 255$ Hz don't show up at all because a UZ-\n   only tip force does not excite UY bending.\n2. **Amplification at expected damping.**  The 2 % modal\n   damping ratio matches published field measurements on\n   similar welded-steel mounting brackets.  At the first peak\n   the amplification is ~25\u00d7; multiplied through Tutorial 1's\n   static stress field this stays well below the 165 MPa\n   allowable, so the bracket clears the dynamic check too.\n   The light-damping curve (red) is shown only to anchor the\n   near-undamped upper bound \u2014 never trust it as the design\n   case unless your measured damping is actually that low.\n3. **Damping-model choice.**  The Rayleigh curve diverges\n   from the modal-damping curves at high frequency \u2014 Rayleigh\n   damping rises monotonically in $\\omega$ while real\n   modal damping is per-mode, so the two only agree over a\n   narrow band.  Use modal damping when the band of interest\n   is narrow (this case); use Rayleigh when one form factor\n   needs to cover three orders of magnitude in $\\omega$\n   (transient analyses where the response carries every mode).\n\nCompanion deeper-dive resources:\n\n* :doc:`/user-guide/solving/harmonic` \u2014 full :meth:`Model.solve_harmonic`\n  API reference with all kwargs.\n* :doc:`/reference/theory/eigenproblem` \u2014 the spectral\n  factorisation that makes modal superposition tractable.\n* :doc:`/getting-started/known-limitations` \u2014 when modal\n  superposition is *not* the right tool (significant non-\n  proportional damping, high modal density, etc.).\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
}