{ "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 5 \u2014 Cyclic-symmetry rotor design check\n\nTutorials 1-3 walked the static + modal workflow on simply-\nshaped components. Real rotating-equipment design has an\nextra wrinkle: the structure is **cyclically symmetric**. An\n$N$-bladed rotor, an $N$-vaned diffuser, and an\n$N$-cylinder reciprocating engine all share the same\nalgebraic property \u2014 the spectrum factors by harmonic index\n$k \\in \\{0, 1, \\ldots, \\lfloor N / 2 \\rfloor\\}$, and a\nsingle base sector plus the cyclic-pairing constraint\nrecovers the full-rotor answer at a fraction of the DOF cost.\n\nYou're sizing a 16-blade impeller for a small turbomachine.\nThe rotor runs at 12 000 rpm continuous. The mechanical\nquestion is: do any of the impeller's natural frequencies sit\ninside an engine-order resonance band? ASME / API-style\ndesign guidance is to keep the lowest natural frequency at\nleast 20 % above the highest credible engine order at\noperating speed (the **Campbell diagram** clearance check).\n\nThis tutorial walks the design check end-to-end:\n\n* **Step 1** \u2014 build a single sector of the impeller.\n* **Step 2** \u2014 wrap it in\n :class:`~femorph_solver.CyclicModel`, which auto-detects\n the low / high cyclic faces from the rotation axis and\n sector angle.\n* **Step 3** \u2014 run a cyclic-symmetry modal solve sweeping\n all harmonic indices.\n* **Step 4** \u2014 assemble the per-(mode, k) frequency table \u2014\n the *cyclic spectrum*.\n* **Step 5** \u2014 compute the Campbell-diagram clearance\n against engine orders 1 \u2026 N at the operating speed.\n* **Step 6** \u2014 render the lowest mode pair as travelling-\n wave snapshots.\n* **Step 7** \u2014 compare DOF count + wall-time against the\n full-rotor solve a non-cyclic approach would have needed.\n\n## Theoretical reference\n\nFor a rotor with $N$ identical sectors, the cyclic-\nsymmetry constraint pairs each base-sector DOF $u_b$\nwith its rotated counterpart $u_c$ via\n\n\\begin{align}:label: cyclic-pairing\n\n u_c(\\alpha = 2\\pi / N) = e^{i\\, k\\, \\alpha}\\, R(\\alpha)\\, u_b,\\end{align}\n\nwhere $R(\\alpha)$ is the rotation matrix and $k$\nis the **harmonic index**. At each $k$ the per-sector\neigenproblem is real-symmetric (after the\nGrimes-Lewis-Simon real-2n augmentation \u2014\n:doc:`/reference/theory/eigenproblem`); sweeping $k$\nfrom 0 to $\\lfloor N/2 \\rfloor$ recovers the full-rotor\nspectrum.\n\nEach $(n, k)$ pair gives a **mode pair** at the same\nfrequency: real and imaginary parts of the complex\neigenvector are travelling-wave snapshots separated by a\nquarter-period. Engine-order excitation at speed\n$\\Omega$ and order $m$ resonates with mode shapes\nwhose harmonic index satisfies $m \\equiv k \\pmod N$\n(Wildheim 1979).\n\nReferences: Thomas 1979, Wildheim 1979, Crandall & Mark 1963\n\u00a73.5, Bathe 2014 \u00a710.3.4.\n\nCompanion gallery example:\n`sphx_glr_gallery_analyses_cyclic_example_cyclic_modes_bladed_rotor.py`\nexercises the same API on a similar rotor without the\ndesign-check angle.\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 one sector of the impeller\n\nSixteen identical wedges + radial blades. We build one\nsector spanning $\\alpha = 2\\pi / 16 = 22.5\u00b0$ and let\nCyclicModel handle the cyclic pairing. The blade is a\nthin radial fin attached to a thicker hub annulus.\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 = 0.05 # m, hub bore (shaft fit)\nHUB_OUTER = 0.10 # m, blade root\nBLADE_TIP = 0.20 # m, blade tip (outer rim)\nT_AXIAL = 0.02 # m, axial thickness\nE, NU, RHO = 2.0e11, 0.30, 7850.0 # steel impeller\n\n# Build a structured (\u03b8, r, z) grid for the hub + blade region.\nn_th = 4 # circumferential resolution within the sector\nn_r_hub = 4 # radial steps in the hub\nn_r_blade = 8 # radial steps in the blade\nn_z = 2 # axial through-thickness\n\n# Concatenated radial coordinate ladder: hub then blade.\nhub_radii = np.linspace(HUB_INNER, HUB_OUTER, n_r_hub + 1)\nblade_radii = np.linspace(HUB_OUTER, BLADE_TIP, n_r_blade + 1)[1:]\nr_ladder = np.concatenate([hub_radii, blade_radii])\n\ntheta = np.linspace(0.0, ALPHA, n_th + 1)\nzs = np.linspace(0.0, T_AXIAL, n_z + 1)\n\n\n# Assemble points + cells in (\u03b8, r, z) ordering. For this\n# tutorial we model the rotor as a uniform annular disk \u2014\n# every cell in the sector is solid material \u2014 to keep the\n# pedagogy focused on the cyclic-symmetry workflow rather\n# than the blade-shape modelling.\n\npts = []\nfor z in zs:\n for th in theta:\n for r in r_ladder:\n pts.append([r * np.cos(th), r * np.sin(th), z])\npts_arr = np.array(pts, dtype=np.float64)\n\nn_th_pts = len(theta)\nn_r_pts = len(r_ladder)\nplane_size = n_th_pts * n_r_pts\n\n\ndef node_idx(kz: int, ti: int, ri: int) -> int:\n return kz * plane_size + ti * n_r_pts + ri\n\n\ncells = []\nfor kz in range(n_z):\n for ti in range(n_th):\n for ri in range(n_r_hub + n_r_blade):\n n00 = node_idx(kz, ti, ri)\n n10 = node_idx(kz, ti, ri + 1)\n n11 = node_idx(kz, ti + 1, ri + 1)\n n01 = node_idx(kz, ti + 1, ri)\n cells.append(\n [\n 8,\n n00,\n n10,\n n11,\n n01,\n n00 + plane_size,\n n10 + plane_size,\n n11 + plane_size,\n n01 + plane_size,\n ]\n )\n\ncells_arr = np.array(cells, dtype=np.int64).ravel()\ncell_types = np.full(len(cells), 12, dtype=np.uint8)\ngrid = pv.UnstructuredGrid(cells_arr, cell_types, pts_arr)\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\n# Pin the hub-bore inner edge in all DOFs (shaft attachment).\nsector_pts = np.asarray(m.grid.points)\nnode_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"], dtype=np.int64)\ninner_radii = np.linalg.norm(sector_pts[:, :2], axis=1)\nhub_bore = np.where(inner_radii < HUB_INNER + 1e-9)[0]\nm.fix(nodes=node_nums[hub_bore].tolist(), dof=\"ALL\")\n\nprint(\"Bladed-impeller cyclic-symmetry modal survey\")\nprint(f\" N_sectors = {N_SECTORS}, sector angle \u03b1 = {np.degrees(ALPHA):.2f}\u00b0\")\nprint(f\" hub: r_inner = {HUB_INNER * 1e3:.0f} mm, r_outer = {HUB_OUTER * 1e3:.0f} mm\")\nprint(f\" blade: r_tip = {BLADE_TIP * 1e3:.0f} mm, axial thickness = {T_AXIAL * 1e3:.0f} mm\")\nprint(f\" Sector mesh: {m.grid.n_points} nodes, {m.grid.n_cells} HEX8-EAS cells\")\nprint(f\" Equivalent full-rotor DOF count: {m.grid.n_points * 3 * N_SECTORS:,}\")\nprint(f\" Cyclic-sector DOF count: {m.grid.n_points * 3:,} (\u00d7{N_SECTORS} reduction)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2 \u2014 wrap in CyclicModel\n\n:class:`~femorph_solver.CyclicModel` auto-detects the low /\nhigh cyclic faces from the rotation axis (default = z) and\nthe sector angle. The DOF-pair list it builds is what the\neigensolver applies as the cyclic constraint at each\nharmonic index.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cyc = femorph_solver.CyclicModel(m, n_sectors=N_SECTORS, axis=(0.0, 0.0, 1.0), angle=ALPHA)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3 \u2014 sweep harmonic indices\n\nFor an N-sector rotor, harmonic indices $k = 0, 1, \u2026,\n\\lfloor N/2 \\rfloor$ cover the full spectrum. Each\n``cyc.solve_modal(harmonic_index=k, n_modes=\u2026)`` call returns\nthe lowest few eigenpairs for that harmonic \u2014 together they\ntile out the (mode_n \u00d7 harmonic_k) frequency map.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N_MODES_PER_K = 4\nall_results = cyc.solve_modal(n_modes=N_MODES_PER_K)\n# Returns a list: one CyclicModalResult per harmonic 0 .. N/2.\nspectrum = np.zeros((N_MODES_PER_K, len(all_results)), dtype=np.float64)\nfor k, res_k in enumerate(all_results):\n spectrum[:, k] = np.sort(np.asarray(res_k.frequency, dtype=np.float64))[:N_MODES_PER_K]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 4 \u2014 assemble the cyclic spectrum table\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print()\nprint(\" Cyclic spectrum (frequencies in Hz):\")\nheader = \" mode\" + \"\".join(f\" {f'k={k}':>10}\" for k in range(N_SECTORS // 2 + 1))\nprint(header)\nprint(\" \" + \"-\" * (len(header) - 2))\nfor n in range(N_MODES_PER_K):\n row = f\" {n + 1:>4}\"\n for k in range(N_SECTORS // 2 + 1):\n row += f\" {spectrum[n, k]:>10.2f}\"\n print(row)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 5 \u2014 Campbell-diagram clearance check\n\nOperating speed \u03a9 = 12 000 rpm = 200 Hz. Engine-order\nexcitation produces forcing at $m \\cdot \\Omega$ for\n$m = 1, 2, 3, \u2026$. A mode at frequency $f_n^{(k)}$\nresonates with engine order $m$ at speed $\\Omega$\nwhen $f_n^{(k)} \\approx m \\cdot \\Omega$ *and*\n$m \\equiv k \\pmod{N_\\text{sectors}}$\n(Wildheim 1979 \u2014 only modes with matching harmonic index\ncouple to a given engine order).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "OMEGA_RPM = 12_000.0\nOMEGA_HZ = OMEGA_RPM / 60.0\nprint()\nprint(f\" Operating speed \u03a9 = {OMEGA_RPM:.0f} rpm = {OMEGA_HZ:.2f} Hz\")\n\n# 20 % design-clearance band around each engine-order resonance.\nclearance_pct = 0.20\n\nviolations = []\nfor k in range(N_SECTORS // 2 + 1):\n for n in range(N_MODES_PER_K):\n f_nk = spectrum[n, k]\n # Engine orders that match this k (m mod N == k or N-k for\n # k > 0; k == 0 only matches m mod N == 0).\n for m_eng in range(1, 2 * N_SECTORS + 1):\n mod = m_eng % N_SECTORS\n if not (mod == k or mod == (N_SECTORS - k) % N_SECTORS):\n continue\n f_eng = m_eng * OMEGA_HZ\n if abs(f_nk - f_eng) / f_nk < clearance_pct:\n violations.append((n + 1, k, m_eng, f_nk, f_eng))\n\nif violations:\n print()\n print(f\" ! Clearance violations (within \u00b1{clearance_pct * 100:.0f} % of an EO resonance):\")\n print(f\" {'mode':>4} {'k':>3} {'m_eng':>6} {'f_n^(k) [Hz]':>14} {'m\u00b7\u03a9 [Hz]':>10}\")\n print(\" \" + \"-\" * 50)\n for vrow in violations[:10]:\n print(f\" {vrow[0]:>4} {vrow[1]:>3} {vrow[2]:>6} {vrow[3]:>13.2f} {vrow[4]:>10.2f}\")\n if len(violations) > 10:\n print(f\" ... + {len(violations) - 10} more\")\nelse:\n print()\n print(\n f\" \u2713 No mode within \u00b1{clearance_pct * 100:.0f} % of any engine-order resonance \u2014 \"\n \"Campbell-diagram clearance OK.\"\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 6 \u2014 render the lowest mode pair (k = 1 travelling wave)\n\nAt harmonic $k = 1$ the lowest mode is a one-nodal-\ndiameter standing pattern that, in the rotating frame,\nbecomes a travelling wave with one nodal diameter. The two\ncomponents (real / imaginary parts of the complex\neigenvector) are quarter-period snapshots of that wave.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res_k1 = all_results[1]\nshapes_k1 = np.asarray(res_k1.mode_shapes, dtype=np.float64)\n\nplotter = pv.Plotter(shape=(1, 2), off_screen=True, window_size=(1100, 480), border=False)\nfor col, idx in enumerate((0, 1)):\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 dof_map = m.dof_map()\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_k1[row, idx]\n peak = float(np.linalg.norm(disp, axis=1).max())\n scale = (0.04 * BLADE_TIP) / 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, scalars=\"|disp|\", cmap=\"viridis\", show_edges=False, show_scalar_bar=False\n )\n plotter.add_text(\n f\"k = 1, mode {idx + 1}\\nf = {np.asarray(res_k1.frequency)[idx]:.1f} Hz\",\n position=\"upper_left\",\n font_size=11,\n )\n plotter.view_xy()\n plotter.camera.zoom(1.05)\nplotter.link_views()\nplotter.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 7 \u2014 DOF / wall-time savings vs full rotor\n\nThe cyclic-symmetry path solves $(N/2 + 1)$ sector-\nsized eigenproblems (one per harmonic) instead of one\nfull-rotor problem. For the impeller above:\n\n- Sector DOFs: ``m.grid.n_points \u00d7 3``\n- Full-rotor DOFs: same \u00d7 N_sectors\n- Direct factor cost scales as $O(N_\\text{dof}^{1.5})$\n for 3-D solid meshes; the cyclic path saves\n $\\sim N^{1.5} / (N/2 + 1)$ factor evaluations.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_sector_dofs = m.grid.n_points * 3\nn_rotor_dofs = n_sector_dofs * N_SECTORS\nfactor_savings = (N_SECTORS**1.5) / (N_SECTORS // 2 + 1)\nprint()\nprint(\" DOF savings:\")\nprint(\n f\" Sector eigen problems \u00d7 {N_SECTORS // 2 + 1} = \"\n f\"{(N_SECTORS // 2 + 1) * n_sector_dofs:,} cumulative DOFs\"\n)\nprint(f\" Full-rotor eigen problem (1 \u00d7) = {n_rotor_dofs:,} DOFs\")\nprint(f\" Factor-cost reduction \u2248 {factor_savings:.1f}\u00d7\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Take-aways\n\n- **Cyclic-symmetry modal** is the right discretisation for\n any rotor / vane row / cylinder array with N-fold\n rotational symmetry. CyclicModel auto-pairs the\n low / high faces; the user supplies axis + angle +\n N_sectors and the eigensolver runs the rest.\n- **The cyclic spectrum** is a (mode_n \u00d7 harmonic_k) table\n that recovers the full-rotor spectrum at a fraction of the\n DOF cost. Sweep $k = 0$ \u2026 $\\lfloor N/2 \\rfloor$.\n- **Campbell-diagram clearance** is the standard\n rotating-equipment design check: keep every\n $f_n^{(k)}$ at least 20 % away from any engine-order\n resonance $m \\cdot \\Omega$ whose order satisfies\n $m \\equiv k \\pmod N$. Wildheim's 1979 mod-rule\n restricts the relevant orders so the search is finite.\n- **Travelling-wave mode pairs** at $k > 0$ are two\n modes at the same frequency, mode shapes 90\u00b0 apart in the\n rotating frame. The complex eigenvector's real and\n imaginary parts are quarter-period snapshots.\n- **DOF / time savings** scale as $N^{1.5} / (N/2 +\n 1)$. For the 16-blade rotor here, the cyclic path is\n ~25 \u00d7 faster than a full-rotor eigen problem \u2014 and\n getting cheaper as N grows.\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 }