{ "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 modal sweep \u2014 annular disk\n\nA 12-sector annular disk has its global eigenspectrum factored\nexactly by **harmonic index** $k \\in \\{0, 1, \\ldots, N/2\\}$\n(Thomas 1979, Wildheim 1979). Each $k$ produces a\nreduced complex-Hermitian eigenproblem on a single base sector;\nfemorph-solver solves the real-symmetric 2n-augmentation form\n(Grimes, Lewis & Simon 1994) so :func:`scipy.sparse.linalg.eigsh`\ncan return mode pairs at every harmonic.\n\nThis example builds a 12-sector disk, runs\n:meth:`CyclicModel.solve_modal` over every harmonic, and prints\nthe lowest non-rigid frequency for each \u2014 illustrating how the\nharmonic sweep recovers the full free-vibration spectrum at\n$1 / N$ the DOF count of an explicit full-rotor model.\n\n## References\n\n* Thomas, D. L. (1979) \"Dynamics of rotationally periodic\n structures,\" *Journal of Sound and Vibration* 66 (4),\n 585\u2013597.\n* Wildheim, S. J. (1979) \"Excitation of rotationally periodic\n structures,\" *Journal of Applied Mechanics* 46, 891\u2013893.\n* Grimes, R. G., Lewis, J. G. and Simon, H. D. (1994)\n \"A shifted block Lanczos algorithm for solving sparse\n symmetric generalized eigenproblems,\" *SIAM J. Matrix\n Analysis and Applications* 15, 228\u2013272 (real 2n\n augmentation for the complex-Hermitian sub-problem).\n* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,\n \u00a710.3.4 (cyclic-symmetry modal analysis).\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 12-sector annular disk\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N_SECTORS = 12\nALPHA = 2.0 * np.pi / N_SECTORS\n\nR_INNER = 0.40\nR_OUTER = 0.55\nT_DISK = 0.02\nN_T, N_R, N_Z = 6, 4, 1\n\nE = 2.0e11\nNU = 0.3\nRHO = 7850.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build the base sector mesh\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "theta = np.linspace(0.0, ALPHA, N_T + 1)\nr = np.linspace(R_INNER, R_OUTER, N_R + 1)\nz = np.linspace(0.0, T_DISK, N_Z + 1)\nR, T, Z = np.meshgrid(r, theta, z, indexing=\"ij\")\nX = R * np.cos(T)\nY = R * np.sin(T)\n\ngrid = pv.StructuredGrid(X, Y, Z).cast_to_unstructured_grid()\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 a 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(\n ELEMENTS.HEX8,\n material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Identify low / high cyclic-face nodes by angular position\n\nThe base sector spans $\\theta \\in [0, \\alpha]$. Nodes\nat $\\theta = 0$ form the **low** face; nodes at\n$\\theta = \\alpha$ form the **high** face.\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]\n\n# Pair the two faces by (radius, z) coordinates so the cyclic\n# constraint links matching points across the sector boundary.\nlow_pts = pts[low_mask]\nhigh_pts = pts[high_mask]\nlow_radii_z = np.column_stack([np.linalg.norm(low_pts[:, :2], axis=1), low_pts[:, 2]])\nhigh_radii_z = np.column_stack([np.linalg.norm(high_pts[:, :2], axis=1), high_pts[:, 2]])\norder_low = np.lexsort(low_radii_z.T[::-1])\norder_high = np.lexsort(high_radii_z.T[::-1])\nlow_pt_idx_paired = np.where(low_mask)[0][order_low]\nhigh_pt_idx_paired = np.where(high_mask)[0][order_high]\nprint(f\"cyclic faces: {len(low_pt_idx_paired)} paired nodes per face\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run the harmonic sweep\n\nSolve the cyclic-symmetry modal problem on the base sector\nat every harmonic $k = 0, 1, \\ldots, N/2$. Each\n:class:`CyclicModalResult` carries the lowest n_modes\neigenvalues for that harmonic.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_modes = 8\ncm = femorph_solver.CyclicModel.wrap(\n m,\n n_sectors=N_SECTORS,\n low_face=low_pt_idx_paired,\n high_face=high_pt_idx_paired,\n)\nresults = cm.solve_modal(n_modes=n_modes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tabulate the lowest non-rigid frequency at each harmonic\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print()\nprint(f\"{'k':>3} {'f_min [Hz]':>12} {'mode shape':>16}\")\nprint(f\"{'-':>3} {'-' * 12:>12} {'-' * 16:>16}\")\nfor r_k in results:\n freqs = np.asarray(r_k.frequency, dtype=np.float64)\n # Discard rigid-body modes (f < 1 Hz) and report the lowest\n # bending-like frequency instead.\n elastic = freqs[freqs > 1.0]\n f_min = float(elastic[0]) if elastic.size else float(\"nan\")\n mode_kind = \"rigid (k=0)\" if r_k.harmonic_index == 0 else f\"k={r_k.harmonic_index} pair\"\n print(f\"{r_k.harmonic_index:>3} {f_min:12.2f} {mode_kind:>16}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the lowest non-rigid mode at the most interesting harmonic\n\nk = 2 typically yields the lowest non-rigid mode on an annular\ndisk \u2014 the classical \"diametrical-2\" pattern with two\nstationary radial nodes on the disc's rim.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "target_k = 2\nchosen = next((r for r in results if r.harmonic_index == target_k), results[1])\nphi = np.asarray(chosen.mode_shapes)[:, 0].reshape(-1, 3)\n# k > 0 mode shapes are complex; take the real part for the\n# stationary snapshot. The travelling-wave companion lives in\n# the imaginary part \u2014 see example_cyclic_modes_bladed_rotor.py.\nphi_real = np.real(phi)\nwarp = grid.copy()\nwarp.points = grid.points + (T_DISK * 4.0 / np.max(np.abs(phi_real))) * phi_real\nwarp[\"UZ\"] = phi_real[:, 2]\n\nf_chosen = float(np.asarray(chosen.frequency)[0])\nplotter = pv.Plotter(off_screen=True, window_size=(720, 480))\nplotter.add_mesh(grid, style=\"wireframe\", color=\"grey\", opacity=0.4)\nplotter.add_mesh(\n warp,\n scalars=\"UZ\",\n cmap=\"coolwarm\",\n show_edges=False,\n scalar_bar_args={\"title\": f\"k = {chosen.harmonic_index} mode 1 \u2014 f = {f_chosen:.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 }