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.

High-level API#

The cyclic-aware front door is femorph_solver.CyclicModel: it wraps a base-sector femorph_solver.Model plus the cyclic parameters (n_sectors and the rotation axis) and detects the low / high cyclic-face node pairing internally. Downstream code does not need to roll its own KD-tree or hand-build DOF index arrays.

import femorph_solver as fs

cm = fs.CyclicModel.from_pv(
    fs.examples.cyclic_bladed_rotor_sector_path(),
    n_sectors=15,
    axis="z",       # default; pass "x"/"y"/"z" or a 3-vector
)
results = cm.modal_solve(n_modes=4)
for r in results:
    print(f"k={r.harmonic_index}  f={r.frequency}")

Each entry of results is a CyclicModalResult for one harmonic index, holding omega_sq, frequency, and the complex base-sector mode_shapes. When you want the full-rotor frequency multiset (each intermediate-k mode counted twice), call cm.aggregated_frequencies(n_modes=4).

Manual face override#

When auto-detection misclassifies the cyclic faces — for example, on a mesh whose interface nodes don’t sit within pair_tol of their rotated partners because of upstream meshing tolerance — pass the low / high faces explicitly as grid-point-index arrays:

cm = fs.CyclicModel.from_pv(
    fs.examples.cyclic_bladed_rotor_sector_path(),
    n_sectors=15,
    low_face=low_indices,
    high_face=high_indices,
)

The arrays must have the same length and pass an R(α)-image-within- pair_tol consistency check; misuses raise ValueError with a clear message naming the offending pair.

Loading a Model from a .pv file#

The bundled rotor sector ships as a single native .pv file: geometry plus element-type assignments, material properties, and unit-system stamp baked in. This is the same format save() writes, so any Model you serialise can come back through from_pv() ready-to-solve, with no APDL setup:

import femorph_solver as fs

m = fs.Model.from_pv("rotor_sector.pv")
r = m.modal_solve(n_modes=4)         # full-sector modal (no cyclic BC)

cm = fs.CyclicModel.from_pv("rotor_sector.pv", n_sectors=15)
r = cm.modal_solve(n_modes=4)        # cyclic sweep across harmonics

Underneath: the building blocks#

femorph_solver.solvers.cyclic.solve_cyclic_modal() is the lower- level entry point if you have a model assembled outside this package or want to drive the cyclic reduction directly:

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,
    harmonic_indices=None,
)

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.

CyclicModel does this assembly + DOF lookup for you: the inner Model exposes stiffness_matrix() and mass_matrix(), the wrapper translates its cached point-pair arrays into global DOF indices, and the spatial rotation is built from the user-supplied axis / angle. Call cm.modal_solve(...) and you get the same CyclicModalResult list back.

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. CyclicModel builds it from the axis / angle kwargs (default: α = 2π / n_sectors about +z) and threads it through the lower-level solver automatically; only when DOFs are already in a per-sector cylindrical frame would you bypass CyclicModel and pass pair_rotation=None to solve_cyclic_modal() directly.

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 lower-level solver itself is element-agnostic: it reads only the (P, d) DOF-index arrays and the (d, d) rotation.

Bundled example — bladed-rotor sector#

The femorph-solver wheel ships a 15-sector × 24° bladed-rotor base sector (101 HEX8 cells) at femorph_solver.examples.cyclic_bladed_rotor_sector_path(). The file is a complete save() snapshot: the HEX8 element-type and isotropic material assignments are baked in, so the demo above runs without any external download or APDL setup.

Cyclic faces on this mesh are not θ=const planes — they’re curved surfaces that map to each other under a 24° rotation about z. CyclicModel resolves the pairing uniformly on the rotated point-cloud match.

The deeper-history version of this example (driven directly off mapdl_archive.examples.sector_archive_file with hand-rolled cKDTree pair detection) is exercised in tests/elements/solid185/test_solid185_cyclic_sector_archive.py as a regression against the legacy code path.