Note
Go to the end to download the full example code.
Cyclic-symmetry travelling-wave pair — bladed rotor#
A bladed rotor’s first non-rigid harmonic produces a travelling-wave mode pair: two modes at the same frequency whose mode shapes differ by a 90° phase shift around the rotor. The pair is what energizes engine-order excitation under operating speed — when the rotor’s harmonic index \(k\) matches the engine order, the travelling-wave amplitude grows resonantly (Crandall & Mark 1963; Bathe §10.3.4).
This example builds a 16-sector ring of stubby radial blades
on an annular hub, runs
Model.cyclic_modal_solve() over every harmonic, and
visualises the lowest non-rigid mode pair at the harmonic
where the rotor is softest — the real / imaginary parts of
the complex mode shape are the two travelling-wave snapshots
separated by a quarter-period.
References#
Thomas, D. L. (1979) “Dynamics of rotationally periodic structures,” J. Sound Vib. 66 (4), 585–597.
Crandall, S. H. and Mark, W. D. (1963) Random Vibration in Mechanical Systems, Academic Press, §3.5 (travelling-wave / standing-wave decomposition).
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.
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 16-blade rotor#
Build the base sector mesh#
Hub: annular sector (HUB_INNER → HUB_OUTER). Blade: radial extension from HUB_OUTER → BLADE_TIP, narrower circumferentially than the hub (BLADE_THICKNESS × ALPHA wide).
Both pieces are HEX8 grids; we concatenate their meshes via pyvista.
n_t_hub, n_r_hub, n_z = 4, 3, 1
theta_hub = np.linspace(0.0, ALPHA, n_t_hub + 1)
r_hub = np.linspace(HUB_INNER, HUB_OUTER, n_r_hub + 1)
z = np.linspace(0.0, T_AXIAL, n_z + 1)
R, T, Z = np.meshgrid(r_hub, theta_hub, z, indexing="ij")
hub = pv.StructuredGrid(R * np.cos(T), R * np.sin(T), Z).cast_to_unstructured_grid()
n_t_blade, n_r_blade = 2, 4
theta_lo = (0.5 - 0.5 * BLADE_THICKNESS) * ALPHA
theta_hi = (0.5 + 0.5 * BLADE_THICKNESS) * ALPHA
theta_blade = np.linspace(theta_lo, theta_hi, n_t_blade + 1)
r_blade = np.linspace(HUB_OUTER, BLADE_TIP, n_r_blade + 1)
R, T, Z = np.meshgrid(r_blade, theta_blade, z, indexing="ij")
blade = pv.StructuredGrid(R * np.cos(T), R * np.sin(T), Z).cast_to_unstructured_grid()
grid = (hub + blade).clean(tolerance=1e-8)
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: 68 nodes, 20 HEX8 cells
full-rotor equivalent: 1088 nodes
Wrap as Model with steel material#
Identify cyclic-face nodes by angular position#
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]
low_pts = pts[low_mask]
high_pts = 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 = low_nodes[order_low]
high_paired = high_nodes[order_high]
print(f"cyclic faces: {len(low_paired)} paired nodes per face")
cyclic faces: 8 paired nodes per face
Run the harmonic sweep#
n_modes = 4
results = m.cyclic_modal_solve(
low_face=low_paired,
high_face=high_paired,
n_sectors=N_SECTORS,
n_modes=n_modes,
)
Find the harmonic with the softest non-rigid mode#
print()
print(f"{'k':>3} {'f_min [Hz]':>12}")
softest_k, softest_f = None, float("inf")
for r_k in results:
freqs = np.asarray(r_k.frequency, dtype=np.float64)
elastic = freqs[freqs > 1.0]
if elastic.size:
f_min = float(elastic[0])
print(f"{r_k.harmonic_index:>3} {f_min:12.2f}")
if r_k.harmonic_index > 0 and f_min < softest_f:
softest_f = f_min
softest_k = r_k.harmonic_index
print()
print(f"softest non-rigid harmonic: k = {softest_k}, f = {softest_f:.2f} Hz")
k f_min [Hz]
0 233.53
1 372.78
2 154.28
3 348.38
4 437.87
5 466.56
6 478.16
7 483.22
8 484.68
softest non-rigid harmonic: k = 2, f = 154.28 Hz
Render the travelling-wave pair (real + imaginary parts)#
For a cyclic-symmetry mode at harmonic \(k\) the real and imaginary parts of \(\boldsymbol{\phi}_{k}\) correspond to the standing-wave snapshots \(\boldsymbol{\phi}_{k}\, \cos(k\, \alpha\, t)\) and \(\boldsymbol{\phi}_{k}\, \sin(k\, \alpha\, t)\) — the two together describe a wave travelling around the rotor’s circumference.
assert softest_k is not None
chosen = next(r for r in results if r.harmonic_index == softest_k)
phi = np.asarray(chosen.mode_shapes)[:, 0].reshape(-1, 3)
phi_re = np.real(phi)
phi_im = np.imag(phi)
amp = max(np.max(np.abs(phi_re)), np.max(np.abs(phi_im))) or 1.0
warp_factor = 0.05 / amp
warp_re = grid.copy()
warp_re.points = grid.points + warp_factor * phi_re
warp_re["uz"] = phi_re[:, 2]
warp_im = grid.copy()
warp_im.points = grid.points + warp_factor * phi_im
warp_im["uz"] = phi_im[:, 2]
f_chosen = float(np.asarray(chosen.frequency)[0])
plotter = pv.Plotter(shape=(1, 2), off_screen=True, window_size=(960, 480))
plotter.subplot(0, 0)
plotter.add_text(f"Re(φ) — k = {softest_k}, f = {f_chosen:.1f} Hz", font_size=10)
plotter.add_mesh(grid, style="wireframe", color="grey", opacity=0.4)
plotter.add_mesh(warp_re, scalars="uz", cmap="coolwarm", show_edges=False, show_scalar_bar=False)
plotter.subplot(0, 1)
plotter.add_text("Im(φ) — quarter-period later", font_size=10)
plotter.add_mesh(grid, style="wireframe", color="grey", opacity=0.4)
plotter.add_mesh(warp_im, scalars="uz", cmap="coolwarm", show_edges=False, show_scalar_bar=False)
plotter.link_views()
plotter.view_isometric()
plotter.show()
Total running time of the script: (0 minutes 0.447 seconds)