Cyclic-symmetry modal sweep — annular disk

A 12-sector annular disk has its global eigenspectrum factored exactly by harmonic index \(k \in \{0, 1, \ldots, N/2\}\) (Thomas 1979, Wildheim 1979). Each \(k\) produces a reduced complex-Hermitian eigenproblem on a single base sector; femorph-solver solves the real-symmetric 2n-augmentation form (Grimes, Lewis & Simon 1994) so scipy.sparse.linalg.eigsh() can return mode pairs at every harmonic.

This example builds a 12-sector disk, runs Model.cyclic_modal_solve() over every harmonic, and prints the lowest non-rigid frequency for each — illustrating how the harmonic sweep recovers the full free-vibration spectrum at \(1 / N\) the DOF count of an explicit full-rotor model.

References

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

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

• Grimes, R. G., Lewis, J. G. and Simon, H. D. (1994) “A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems,” SIAM J. Matrix Analysis and Applications 15, 228–272 (real 2n augmentation for the complex-Hermitian sub-problem).

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

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver
from femorph_solver import ELEMENTS

Geometry — one base sector of a 12-sector annular disk

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

R_INNER = 0.40
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

Build the base sector mesh

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()
print(f"sector mesh: {grid.n_points} nodes, {grid.n_cells} HEX8 cells")
print(f"full-rotor equivalent: {N_SECTORS * grid.n_points} nodes")

sector mesh: 70 nodes, 24 HEX8 cells
full-rotor equivalent: 840 nodes

Wrap as a Model with steel material

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

Identify low / high cyclic-face nodes by angular position

The base sector spans \(\theta \in [0, \alpha]\). Nodes at \(\theta = 0\) form the low face; nodes at \(\theta = \alpha\) form the high face.

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_nodes = node_nums[low_mask]
high_nodes = node_nums[high_mask]

# Pair the two faces by (radius, z) coordinates so the cyclic
# constraint links matching points across the sector boundary.
low_pts = pts[low_mask]
high_pts = pts[high_mask]
low_radii_z = np.column_stack([np.linalg.norm(low_pts[:, :2], axis=1), low_pts[:, 2]])
high_radii_z = np.column_stack([np.linalg.norm(high_pts[:, :2], axis=1), high_pts[:, 2]])
order_low = np.lexsort(low_radii_z.T[::-1])
order_high = np.lexsort(high_radii_z.T[::-1])
low_nodes_paired = low_nodes[order_low]
high_nodes_paired = high_nodes[order_high]
print(f"cyclic faces: {len(low_nodes_paired)} paired nodes per face")

cyclic faces: 10 paired nodes per face

Run the harmonic sweep

Solve the cyclic-symmetry modal problem on the base sector at every harmonic \(k = 0, 1, \ldots, N/2\). Each CyclicModalResult carries the lowest n_modes eigenvalues for that harmonic.

n_modes = 8
results = m.cyclic_modal_solve(
    low_face=low_nodes_paired,
    high_face=high_nodes_paired,
    n_sectors=N_SECTORS,
    n_modes=n_modes,
)

Tabulate the lowest non-rigid frequency at each harmonic

print()
print(f"{'k':>3}  {'f_min [Hz]':>12}  {'mode shape':>16}")
print(f"{'-':>3}  {'-' * 12:>12}  {'-' * 16:>16}")
for r_k in results:
    freqs = np.asarray(r_k.frequency, dtype=np.float64)
    # Discard rigid-body modes (f < 1 Hz) and report the lowest
    # bending-like frequency instead.
    elastic = freqs[freqs > 1.0]
    f_min = float(elastic[0]) if elastic.size else float("nan")
    mode_kind = "rigid (k=0)" if r_k.harmonic_index == 0 else f"k={r_k.harmonic_index} pair"
    print(f"{r_k.harmonic_index:>3}  {f_min:12.2f}  {mode_kind:>16}")

k    f_min [Hz]        mode shape
-  ------------  ----------------
0        328.21       rigid (k=0)
1        438.02          k=1 pair
2         76.12          k=2 pair
3        218.20          k=3 pair
4        421.23          k=4 pair
5        684.49          k=5 pair
6       1008.36          k=6 pair

Render the lowest non-rigid mode at the most interesting harmonic

k = 2 typically yields the lowest non-rigid mode on an annular disk — the classical “diametrical-2” pattern with two stationary radial nodes on the disc’s rim.

target_k = 2
chosen = next((r for r in results if r.harmonic_index == target_k), results[1])
phi = np.asarray(chosen.mode_shapes)[:, 0].reshape(-1, 3)
# k > 0 mode shapes are complex; take the real part for the
# stationary snapshot.  The travelling-wave companion lives in
# the imaginary part — see example_cyclic_modes_bladed_rotor.py.
phi_real = np.real(phi)
warp = grid.copy()
warp.points = grid.points + (T_DISK * 4.0 / np.max(np.abs(phi_real))) * phi_real
warp["UZ"] = phi_real[:, 2]

f_chosen = float(np.asarray(chosen.frequency)[0])
plotter = pv.Plotter(off_screen=True, window_size=(720, 480))
plotter.add_mesh(grid, style="wireframe", color="grey", opacity=0.4)
plotter.add_mesh(
    warp,
    scalars="UZ",
    cmap="coolwarm",
    show_edges=False,
    scalar_bar_args={"title": f"k = {chosen.harmonic_index} mode 1 — f = {f_chosen:.1f} Hz"},
)
plotter.view_isometric()
plotter.camera.zoom(1.05)
plotter.show()