{
  "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# Cyclic-symmetry travelling-wave pair \u2014 bladed rotor\n\nA bladed rotor's first non-rigid harmonic produces a\n**travelling-wave mode pair**: two modes at the same frequency\nwhose mode shapes differ by a 90\u00b0 phase shift around the\nrotor.  The pair is what energizes engine-order excitation\nunder operating speed \u2014 when the rotor's harmonic index\n$k$ matches the engine order, the travelling-wave\namplitude grows resonantly (Crandall & Mark 1963; Bathe\n\u00a710.3.4).\n\nThis example builds a 16-sector ring of stubby radial blades\non an annular hub, runs\n:meth:`CyclicModel.solve_modal` over every harmonic, and\nvisualises the lowest non-rigid mode pair at the harmonic\nwhere the rotor is softest \u2014 the real / imaginary parts of\nthe complex mode shape are the two travelling-wave snapshots\nseparated by a quarter-period.\n\n## References\n\n* Thomas, D. L. (1979) \"Dynamics of rotationally periodic\n  structures,\" *J. Sound Vib.* 66 (4), 585\u2013597.\n* Crandall, S. H. and Mark, W. D. (1963) *Random Vibration in\n  Mechanical Systems*, Academic Press, \u00a73.5\n  (travelling-wave / standing-wave decomposition).\n* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,\n  \u00a710.3.4 (cyclic-symmetry modal).\n* Wildheim, S. J. (1979) \"Excitation of rotationally periodic\n  structures,\" *J. Appl. Mech.* 46, 891\u2013893.\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": [
        "## Geometry \u2014 one base sector of a 16-blade rotor\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_SECTORS = 16\nALPHA = 2.0 * np.pi / N_SECTORS\nHUB_INNER, HUB_OUTER = 0.20, 0.30\nBLADE_TIP = 0.50\nBLADE_THICKNESS = 0.06  # circumferential, fraction of sector\nT_AXIAL = 0.02\n\nE = 2.0e11\nNU = 0.3\nRHO = 7850.0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the base sector mesh\n\nHub: annular sector (HUB_INNER \u2192 HUB_OUTER).\nBlade: radial extension from HUB_OUTER \u2192 BLADE_TIP, narrower\ncircumferentially than the hub (BLADE_THICKNESS \u00d7 ALPHA wide).\n\nBoth pieces are HEX8 grids; we concatenate their meshes via\npyvista.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_t_hub, n_r_hub, n_z = 4, 3, 1\ntheta_hub = np.linspace(0.0, ALPHA, n_t_hub + 1)\nr_hub = np.linspace(HUB_INNER, HUB_OUTER, n_r_hub + 1)\nz = np.linspace(0.0, T_AXIAL, n_z + 1)\nR, T, Z = np.meshgrid(r_hub, theta_hub, z, indexing=\"ij\")\nhub = pv.StructuredGrid(R * np.cos(T), R * np.sin(T), Z).cast_to_unstructured_grid()\n\nn_t_blade, n_r_blade = 2, 4\ntheta_lo = (0.5 - 0.5 * BLADE_THICKNESS) * ALPHA\ntheta_hi = (0.5 + 0.5 * BLADE_THICKNESS) * ALPHA\ntheta_blade = np.linspace(theta_lo, theta_hi, n_t_blade + 1)\nr_blade = np.linspace(HUB_OUTER, BLADE_TIP, n_r_blade + 1)\nR, T, Z = np.meshgrid(r_blade, theta_blade, z, indexing=\"ij\")\nblade = pv.StructuredGrid(R * np.cos(T), R * np.sin(T), Z).cast_to_unstructured_grid()\n\ngrid = (hub + blade).clean(tolerance=1e-8)\nprint(f\"sector mesh: {grid.n_points} nodes, {grid.n_cells} HEX8 cells\")\nprint(f\"full-rotor equivalent: {N_SECTORS * grid.n_points} nodes\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Wrap as Model with steel material\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m = femorph_solver.Model.from_grid(grid)\nm.assign(ELEMENTS.HEX8, material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Identify cyclic-face nodes by angular position\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pts = np.asarray(grid.points)\nnode_nums = np.asarray(grid.point_data[\"ansys_node_num\"], dtype=np.int64)\nangle = np.arctan2(pts[:, 1], pts[:, 0])\ntol = 1e-6\nlow_mask = np.abs(angle) < tol\nhigh_mask = np.abs(angle - ALPHA) < tol\nlow_nodes = node_nums[low_mask]\nhigh_nodes = node_nums[high_mask]\nlow_pts = pts[low_mask]\nhigh_pts = pts[high_mask]\norder_low = np.lexsort(\n    np.column_stack([np.linalg.norm(low_pts[:, :2], axis=1), low_pts[:, 2]]).T[::-1]\n)\norder_high = np.lexsort(\n    np.column_stack([np.linalg.norm(high_pts[:, :2], axis=1), high_pts[:, 2]]).T[::-1]\n)\nlow_paired = np.where(low_mask)[0][order_low]\nhigh_paired = np.where(high_mask)[0][order_high]\nprint(f\"cyclic faces: {len(low_paired)} paired nodes per face\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run the harmonic sweep\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_modes = 4\ncm = femorph_solver.CyclicModel.wrap(\n    m,\n    n_sectors=N_SECTORS,\n    low_face=low_paired,\n    high_face=high_paired,\n)\nresults = cm.solve_modal(n_modes=n_modes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Find the harmonic with the softest non-rigid mode\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(f\"{'k':>3}  {'f_min [Hz]':>12}\")\nsoftest_k, softest_f = None, float(\"inf\")\nfor r_k in results:\n    freqs = np.asarray(r_k.frequency, dtype=np.float64)\n    elastic = freqs[freqs > 1.0]\n    if elastic.size:\n        f_min = float(elastic[0])\n        print(f\"{r_k.harmonic_index:>3}  {f_min:12.2f}\")\n        if r_k.harmonic_index > 0 and f_min < softest_f:\n            softest_f = f_min\n            softest_k = r_k.harmonic_index\n\nprint()\nprint(f\"softest non-rigid harmonic: k = {softest_k}, f = {softest_f:.2f} Hz\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the travelling-wave pair (real + imaginary parts)\n\nFor a cyclic-symmetry mode at harmonic $k$ the\nreal and imaginary parts of $\\boldsymbol{\\phi}_{k}$\ncorrespond to the standing-wave snapshots\n$\\boldsymbol{\\phi}_{k}\\, \\cos(k\\, \\alpha\\, t)$ and\n$\\boldsymbol{\\phi}_{k}\\, \\sin(k\\, \\alpha\\, t)$ \u2014\nthe two together describe a wave travelling around the\nrotor's circumference.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "assert softest_k is not None\nchosen = next(r for r in results if r.harmonic_index == softest_k)\nphi = np.asarray(chosen.mode_shapes)[:, 0].reshape(-1, 3)\nphi_re = np.real(phi)\nphi_im = np.imag(phi)\namp = max(np.max(np.abs(phi_re)), np.max(np.abs(phi_im))) or 1.0\nwarp_factor = 0.05 / amp\n\nwarp_re = grid.copy()\nwarp_re.points = grid.points + warp_factor * phi_re\nwarp_re[\"uz\"] = phi_re[:, 2]\n\nwarp_im = grid.copy()\nwarp_im.points = grid.points + warp_factor * phi_im\nwarp_im[\"uz\"] = phi_im[:, 2]\n\nf_chosen = float(np.asarray(chosen.frequency)[0])\nplotter = pv.Plotter(shape=(1, 2), off_screen=True, window_size=(960, 480))\nplotter.subplot(0, 0)\nplotter.add_text(f\"Re(\u03c6) \u2014 k = {softest_k}, f = {f_chosen:.1f} Hz\", font_size=10)\nplotter.add_mesh(grid, style=\"wireframe\", color=\"grey\", opacity=0.4)\nplotter.add_mesh(warp_re, scalars=\"uz\", cmap=\"coolwarm\", show_edges=False, show_scalar_bar=False)\nplotter.subplot(0, 1)\nplotter.add_text(\"Im(\u03c6) \u2014 quarter-period later\", font_size=10)\nplotter.add_mesh(grid, style=\"wireframe\", color=\"grey\", opacity=0.4)\nplotter.add_mesh(warp_im, scalars=\"uz\", cmap=\"coolwarm\", show_edges=False, show_scalar_bar=False)\nplotter.link_views()\nplotter.view_isometric()\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
}