Cyclic-symmetry modal analysis#

A rotor built from N identical sectors spanning 360° admits a cyclic-symmetry modal reduction: every full-rotor eigenmode satisfies

\[u_{\text{sector } j+1}(\text{global}) = e^{i k \alpha}\, R(\alpha)\, u_{\text{sector } j}(\text{global})\]

with \(\alpha = 2\pi / N\) and \(R(\alpha)\) the spatial rotation about the symmetry axis. The integer \(k \in \{0, 1, \dots, N/2\}\) is the harmonic index (or nodal diameter). Solving one sector per k reproduces the full-rotor spectrum at a cost of \(1 / N\) per harmonic — for a 15-sector bladed rotor the eight needed sector solves already beat a single full-rotor eigensolve.

API#

femorph_solver.solvers.cyclic.solve_cyclic_modal() takes the base-sector stiffness / mass plus a one-to-one DOF pairing between the low-angle and high-angle cyclic faces:

from femorph_solver.solvers.cyclic import solve_cyclic_modal

results = solve_cyclic_modal(
    K,                    # (N, N) sector stiffness, real SPD
    M,                    # (N, N) sector mass, real SPD
    low_node_dofs,        # (P, d) global DOFs at P paired low-face nodes
    high_node_dofs,       # (P, d) paired high-face DOFs (same labels, same order)
    n_sectors=N,
    n_modes=10,
    pair_rotation=R_alpha,  # (d, d) spatial rotation, or None for local frames
    harmonic_indices=None,  # defaults to [0, 1, …, N//2]
)

Each entry of results is a CyclicModalResult for one harmonic index, holding omega_sq, frequency, and the complex base-sector mode_shapes.

The constraint matrix is assembled directly in COO form: low-face and high-face DOFs are identified, the high-face DOFs are eliminated via the phase-rotation relation, and the reduced Hermitian generalised problem is solved with scipy.linalg.eigh() (dense, subset by index). For k = 0 and k = N/2 (even N) the phase is \(\pm 1\) and the reduced system is real; intermediate k give a complex Hermitian reduction.

Harmonic-index counting#

The cyclic sweep is equivalent to the full-rotor eigenproblem — every full-rotor frequency falls in exactly one harmonic index — but modes are degenerate in pairs for intermediate k:

  • k = 0: single real eigenmodes (the “in-phase” family)

  • k = N/2 (even N only): single real eigenmodes (anti-phase)

  • 0 < k < N/2: each eigenvalue represents a travelling-wave pair in the full rotor; count it twice when aggregating to the full-rotor multiset.

Rigid-body modes are partitioned across k = 0 (axial translation + axial rotation) and k = 1 (the in-plane translation pair).

Why the spatial rotation matters#

Stiffness / mass are usually assembled in a global XYZ frame, so the phase relation between neighbouring sectors picks up not just the scalar exp(i k α) but also the rotation \(R(\alpha)\) that takes sector j to sector j+1. Pass it as the (d, d) matrix acting on each paired node’s d translational DOFs (3 for a 3D solid sweeping about z). Omit it (None) only when face DOFs are already in a per-sector cylindrical frame, where \(R(\alpha)\) reduces to the identity.

Elements with both translational and rotational DOFs (shells, beams) use a block-diagonal R of two copies of the spatial rotation — one for the translations, one for the rotations. The solver itself is element-agnostic: it reads only the (P, d) DOF-index arrays and the (d, d) rotation.

Example — bladed-rotor sector#

The femorph-solver test suite ships the 15-sector bladed rotor sector from mapdl_archive.examples.sector_archive_file (101 SOLID185 hex elements per sector, α = 24°) along with a port of femorph’s rotation-matching pair detection. The end-to-end flow:

import mapdl_archive
import numpy as np
from scipy.spatial import cKDTree

import femorph_solver
from femorph_solver.solvers.cyclic import solve_cyclic_modal

a = mapdl_archive.Archive(mapdl_archive.examples.sector_archive_file)
grid = a.grid.extract_cells(np.where(a.grid.celltypes == 12)[0]).cast_to_unstructured_grid()

# Auto-detect N and low/high node pairing.
pts = np.asarray(grid.points)
tree = cKDTree(pts)
N = 15                                  # or scan N via rotation+match
alpha = 2 * np.pi / N
R = np.array([[np.cos(alpha), -np.sin(alpha), 0.0],
              [np.sin(alpha),  np.cos(alpha), 0.0],
              [0.0,            0.0,           1.0]])
d, idx = tree.query(pts @ R.T, k=1)
low_pt = np.where(d < 1e-4)[0]
high_pt = idx[low_pt]

# Build the sector model and look up UX/UY/UZ DOF indices per pair.
m = femorph_solver.Model.from_grid(grid)
m.et(185, "SOLID185")
m.mp("EX", 1, 16.9e6);  m.mp("PRXY", 1, 0.31);  m.mp("DENS", 1, 4.1408e-4)

dof_map = m.dof_map()
def node_dofs(nn):
    nn_col = dof_map[:, 0]
    return np.array([int(np.where((nn_col == n) & (dof_map[:, 1] == j))[0][0])
                     for n in nn for j in (0, 1, 2)]).reshape(-1, 3)

nns = np.asarray(grid.point_data["ansys_node_num"])
low_nd = node_dofs(nns[low_pt])
high_nd = node_dofs(nns[high_pt])

results = solve_cyclic_modal(
    m.stiffness_matrix(), m.mass_matrix(),
    low_nd, high_nd,
    n_sectors=N, n_modes=4,
    pair_rotation=R,
)
for r in results:
    print(f"k={r.harmonic_index}  f={r.frequency}")

That full example is exercised in tests/elements/solid185/test_solid185_cyclic_sector_archive.py, which verifies the sector sweep reproduces the full 15-sector rotor’s spectrum to within 0.5 % on the first 10 elastic modes.