Cyclic-symmetry modal on a rotor sector

A rotor with \(N\) identical sectors spanning 360° has a spectrum that factors by harmonic index \(k \in \{0, 1, \dots, N/2\}\). Analysing one base sector plus the constraint

\[u_\text{high} = e^{i k \alpha}\, R(\alpha)\, u_\text{low}, \qquad \alpha = 2\pi / N,\]

recovers the full-rotor spectrum at a fraction of the DOF count.

This example builds a simple 4-sector rotor from a pyvista annular sector, then runs solve_cyclic_modal() across all harmonic indices.

from __future__ import annotations

import numpy as np
import pyvista as pv
from scipy.spatial import cKDTree

import femorph_solver
from femorph_solver.solvers.cyclic import solve_cyclic_modal

Build a 90° annular sector

N_SECTORS = 4
ALPHA = 2.0 * np.pi / N_SECTORS

N_T, N_R, N_Z = 12, 6, 3
theta = np.linspace(0.0, ALPHA, N_T + 1)
r = np.linspace(0.8, 1.0, N_R + 1)
z = np.linspace(0.0, 0.1, N_Z + 1)

# (r, θ, z) order gives a CCW-from-above hex winding in the (r, θ)
# plane — swapping r and θ would flip the winding inside-out.
R, T, Z = np.meshgrid(r, theta, z, indexing="ij")
X = R * np.cos(T)
Y = R * np.sin(T)

grid = pv.StructuredGrid(X, Y, Z).cast_to_unstructured_grid()

Build the femorph-solver model and assemble

m = femorph_solver.Model.from_grid(grid)
m.et(1, "SOLID185")
m.mp("EX", 1, 2.1e11)
m.mp("PRXY", 1, 0.30)
m.mp("DENS", 1, 7850.0)

K = m.stiffness_matrix()
M = m.mass_matrix()

/home/runner/_work/solver/solver/examples/analyses/cyclic/example_rotor_sector.py:54: DeprecationWarning: Model.et(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).et(et_id, name)` for line-by-line APDL deck porting, or the native `Model.assign("HEX8", material)` for new code.
  m.et(1, "SOLID185")
/home/runner/_work/solver/solver/examples/analyses/cyclic/example_rotor_sector.py:55: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("EX", 1, 2.1e11)
/home/runner/_work/solver/solver/examples/analyses/cyclic/example_rotor_sector.py:56: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("PRXY", 1, 0.30)
/home/runner/_work/solver/solver/examples/analyses/cyclic/example_rotor_sector.py:57: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("DENS", 1, 7850.0)

Identify cyclic-face node pairs

pts = np.asarray(grid.points)
low = np.where(np.abs(np.arctan2(pts[:, 1], pts[:, 0])) < 1e-6)[0]
high = np.where(np.abs(np.arctan2(pts[:, 1], pts[:, 0]) - ALPHA) < 1e-6)[0]

c, s = np.cos(ALPHA), np.sin(ALPHA)
R_sector = np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]])
low_rot = pts[low] @ R_sector.T
_, pair_idx = cKDTree(pts[high]).query(low_rot)
high_paired = high[pair_idx]

dof_map = m.dof_map()
nn_col = dof_map[:, 0]
dof_col = dof_map[:, 1]
node_nums = np.asarray(grid.point_data["ansys_node_num"])
low_dofs = np.empty((len(low), 3), dtype=np.int64)
high_dofs = np.empty((len(low), 3), dtype=np.int64)
for i, (lo, hi) in enumerate(zip(low, high_paired)):
    nn_lo, nn_hi = int(node_nums[lo]), int(node_nums[hi])
    for d in range(3):
        low_dofs[i, d] = int(np.where((nn_col == nn_lo) & (dof_col == d))[0][0])
        high_dofs[i, d] = int(np.where((nn_col == nn_hi) & (dof_col == d))[0][0])

Solve for every harmonic index

results = solve_cyclic_modal(
    K,
    M,
    low_dofs,
    high_dofs,
    n_sectors=N_SECTORS,
    n_modes=4,
    pair_rotation=R_sector,
)

print(f"Cyclic modal on {N_SECTORS}-sector rotor:")
for res in results:
    fs = ", ".join(f"{f:.2f}" for f in res.frequency[:4])
    print(f"  k = {res.harmonic_index}: f [Hz] = [{fs}]")

Cyclic modal on 4-sector rotor:
  k = 0: f [Hz] = [0.00, 0.00, 503.32, 503.32]
  k = 1: f [Hz] = [0.00, 0.00, 262.09, 452.92]
  k = 2: f [Hz] = [91.75, 91.75, 164.68, 164.68]

Visualise mode 0 of k = 1

r1 = next(res for res in results if res.harmonic_index == 1)
phi = r1.mode_shapes[:, 0].reshape(-1, 3)
grid.point_data["mode_k1"] = np.abs(phi)

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(
    grid.warp_by_vector("mode_k1", factor=0.05 / (np.max(np.abs(phi)) + 1e-30)),
    scalars=np.linalg.norm(phi, axis=1),
    cmap="viridis",
    scalar_bar_args={"title": f"|mode| (k=1, f={r1.frequency[0]:.1f} Hz)"},
)
plotter.add_axes()
plotter.show()