Cyclic-symmetry mode family across every harmonic

A 12-sector annular disk’s modal spectrum decomposes per harmonic index \(k \in \{0, 1, \ldots, N/2\}\). Each harmonic’s lowest non-rigid mode picks out a distinct nodal diameter pattern: \(k = 0\) is the breathing (axisymmetric) mode, \(k = 1\) is the rigid-tilt-like mode, \(k = 2\) is the diametrical-2 (two stationary radial nodes), \(k = 3\) adds a third nodal diameter, and so on up to \(k = N / 2 = 6\), the highest harmonic sustainable on a 12-sector rotor.

This example sweeps every harmonic, picks the lowest non- rigid mode at each \(k\), and renders all seven side-by- side in a single 4×2 viewport grid — the cyclic-symmetry mode family of the disk.

References

• Thomas, D. L. (1979) “Dynamics of rotationally periodic structures,” J. Sound Vib. 66 (4), 585–597.

• Bathe, K.-J. (2014) Finite Element Procedures, 2nd ed., §10.3.4 (cyclic-symmetry modal).

• Wildheim, S. J. (1979) “Excitation of rotationally periodic structures,” J. Appl. Mech. 46, 891–893.

• Singh, M. P. and Vakakis, A. F. (1993) “On the dynamics of cyclic-symmetric structures,” Mechanism and Machine Theory 28 (5), 627–637 (mode classification by harmonic index).

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver
from femorph_solver import ELEMENTS

Geometry — 12-sector annular disk

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

R_INNER = 0.30
R_OUTER = 0.55
T_DISK = 0.02
N_T, N_R, N_Z = 6, 4, 1

E = 2.0e11
NU = 0.3
RHO = 7850.0

theta = np.linspace(0.0, ALPHA, N_T + 1)
r = np.linspace(R_INNER, R_OUTER, N_R + 1)
z = np.linspace(0.0, T_DISK, N_Z + 1)
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()

m = femorph_solver.Model.from_grid(grid)
m.assign(ELEMENTS.HEX8, material={"EX": E, "PRXY": NU, "DENS": RHO})

Cyclic-face node pairing

pts = np.asarray(grid.points)
node_nums = np.asarray(grid.point_data["ansys_node_num"], dtype=np.int64)
angle = np.arctan2(pts[:, 1], pts[:, 0])
tol = 1e-6
low_mask = np.abs(angle) < tol
high_mask = np.abs(angle - ALPHA) < tol
low_pts, high_pts = pts[low_mask], pts[high_mask]
order_low = np.lexsort(
    np.column_stack([np.linalg.norm(low_pts[:, :2], axis=1), low_pts[:, 2]]).T[::-1]
)
order_high = np.lexsort(
    np.column_stack([np.linalg.norm(high_pts[:, :2], axis=1), high_pts[:, 2]]).T[::-1]
)
low_paired = node_nums[low_mask][order_low]
high_paired = node_nums[high_mask][order_high]

Run the harmonic sweep

results = m.cyclic_modal_solve(
    low_face=low_paired,
    high_face=high_paired,
    n_sectors=N_SECTORS,
    n_modes=4,
)

Pick the lowest non-rigid mode per harmonic

family: list[tuple[int, float, np.ndarray]] = []
for r_k in results:
    freqs = np.asarray(r_k.frequency, dtype=np.float64)
    shapes = np.asarray(r_k.mode_shapes)
    keep = np.where(freqs > 1.0)[0]
    if keep.size:
        idx = int(keep[0])
        phi = shapes[:, idx].reshape(-1, 3)
        family.append((r_k.harmonic_index, float(freqs[idx]), phi))

print(f"recovered {len(family)} non-rigid modes (one per harmonic):")
for k, f, _ in family:
    print(f"  k = {k:>2}   f = {f:8.2f} Hz")

recovered 7 non-rigid modes (one per harmonic):
  k =  0   f =   207.35 Hz
  k =  1   f =   307.51 Hz
  k =  2   f =    86.37 Hz
  k =  3   f =   240.57 Hz
  k =  4   f =   456.45 Hz
  k =  5   f =   731.79 Hz
  k =  6   f =  1065.60 Hz

Render all harmonics in a single 4 × 2 viewport grid

For visual consistency, normalise each mode shape to unit peak amplitude so the warp factor reads the same regardless of the eigsolver’s mass-orthonormalisation scaling.

n_panels = len(family)
rows = (n_panels + 1) // 2
plotter = pv.Plotter(shape=(rows, 2), off_screen=True, window_size=(960, 240 * rows))

for idx, (k, f, phi) in enumerate(family):
    row, col = divmod(idx, 2)
    plotter.subplot(row, col)
    phi_re = np.real(phi)
    peak = float(np.max(np.abs(phi_re))) or 1.0
    warp = grid.copy()
    warp.points = grid.points + (T_DISK * 4.0 / peak) * phi_re
    warp["UZ"] = phi_re[:, 2]
    plotter.add_mesh(grid, style="wireframe", color="grey", opacity=0.4)
    plotter.add_mesh(warp, scalars="UZ", cmap="coolwarm", show_edges=False, show_scalar_bar=False)
    plotter.add_text(f"k = {k}: f = {f:.0f} Hz", position="upper_edge", font_size=10)

# Hide any unused viewports.
for empty_idx in range(n_panels, rows * 2):
    row, col = divmod(empty_idx, 2)
    plotter.subplot(row, col)
    plotter.add_text("(unused)", position="upper_edge", font_size=8, color="white")

plotter.link_views()
plotter.view_isometric()
plotter.show()