How the LabyrinthSeal Code Works

How the LabyrinthSeal Code Works#

This tutorial is a developer / code walkthrough of ross.seals.labyrinth_seal.LabyrinthSeal. Unlike tutorial_seal, which shows how to use the seal models, this notebook explains how the code is built — the architecture, the solver pipeline, and what each method does — with runnable examples that open up the internal state at every stage.

By the end you should be able to read labyrinth_seal.py from top to bottom and know exactly what every method contributes to the final rotordynamic coefficients.

Physics in one paragraph. A labyrinth seal limits leakage by forcing the gas through a series of narrow throttlings (teeth). Between teeth, the flow expands into cavities and dissipates kinetic energy. ROSS models this as a 1-D compressible bulk-flow problem: it first solves the base flow (leakage rate, pressure / temperature / density at every cavity, and the tangential swirl velocity), then perturbs the rotor position by a small whirl orbit and linearizes the governing equations to extract stiffness and damping coefficients (kxx, kxy, cxx, cxy, …).

1. Architecture#

1.1 Two classes — a solver and an element#

Element                  (ross/element.py, abstract base)
  └─ BearingElement      (ross/bearing_seal_element.py) — stores k/c/m coefficients vs speed (and whirl frequency),
       │                   provides K(), C(), M(), G(), plot(), interpolation
       └─ SealElement     (ross/bearing_seal_element.py) — a bearing handled specially by Rotor
            │              (no static reaction force; removed/replaced in Level-1 stability)
            └─ LabyrinthSeal  (ross/seals/labyrinth_seal.py) — thin element wrapper

LabyrinthSolver           (ross/seals/labyrinth_seal.py) — the physics, held as seal.solver

The responsibilities are split:

  1. LabyrinthSolver owns the physics: given geometry, gas and operating conditions it computes the leakage and dynamic coefficients for one rotor speed at a time (solve(speed, frequency=None)), keeping all mutable flow state (pressures, swirl velocities, gradient tables) on itself.

  2. LabyrinthSeal is the ROSS element: it resolves gas properties, builds one solver, maps it over the requested speeds, and hands the coefficient lists to SealElement.__init__. The resulting object is a regular ROSS element — drop it into rs.Rotor(...) and it inherits K(), C(), plot(), coefficient interpolation, etc.

The solver holds only plain data, so multi-speed runs can pickle it to worker processes.

1.2 The constructor pipeline#

LabyrinthSeal.__init__ (decorated with @check_units so pint Q_ inputs are converted to SI) does four things:

  1. Resolve gas properties. If gas_composition is given, it builds two ccp.State objects (inlet, and outlet at constant enthalpy) and derives the molar mass, gamma, and the two-state references reference_temperatures / reference_viscosities. Otherwise you pass molar_mass, gamma, reference_temperatures, reference_viscosities directly. The specific gas constant is R = 8314 / molar_mass. The thermodynamic backend (IdealGas or RealGas, selected by gas_model) comes from ross/seals/gas_model.py.

  2. Validate and store the geometry (seal_type, n_teeth >= 2, tooth dimensions) and build a LabyrinthSolver with it. The solver arrays are sized from the geometry (n_stations = n_teeth + 1), so there is no fixed tooth-count limit.

  3. Solve once per speed. Unless coefficients were passed explicitly via kwargs, it maps solver.solve over speed (solve_frequencies uses a multiprocessing.Pool for more than 4 speeds and stays sequential below that). Each solve returns a dict; these are stacked into per-coefficient lists aligned with speed. With frequency= given, solver.solve_grid is used instead and the lists become 2-D (see section 6).

  4. Hand the results to SealElement. super().__init__(n, speed=..., frequency=..., **coefficients_dict) stores everything as speed-aligned tables and sets up interpolation.

1.3 The solve() pipeline#

LabyrinthSolver.solve(speed, frequency=None) is the heart of the model. For a single rotor speed (and whirl frequency, by default equal to the speed) it executes a fixed sequence of physics stages:

solve(speed, frequency=None)
  ├─ _solve_base_flow(speed)
  │    ├─ inlet_swirl_velocity = preswirl · Ω · shaft_radius
  │    ├─ _reset_state()                 → fresh state arrays, overall pressure ratio
  │    ├─ _vermes_leakage()              → Vermes leakage estimate → initial mass-flux guess (mdot)
  │    ├─ _solve_pressure_distribution() → pressure / temperature / density per cavity; refines mdot
  │    └─ _solve_swirl_velocities()      → cavity swirl + wall shear; builds the cg / cx gradient tables
  └─ _solve_whirl(frequency)           (whirl frequency ω; defaults to Ω)
       └─ _solve_perturbation_system()   → perturb + linearize → kxx, kxy, cxx, cxy
  → returns {kxx, kyy, kxy, kyx, cxx, cyy, cxy, cyx, seal_leakage, pressure, ...}

Notice the two-phase structure: the first three stages solve the steady base flow, and the perturbation stage linearizes it. The map below ties each physical concept to the method that owns it.

Stage

Method

Physical role

Key state it produces

Reset

_reset_state()

fresh state arrays for a new speed

overall_pressure_ratio

Leakage seed

_vermes_leakage() (+ _solve_choked_flow_function())

Vermes discharge / carry-over model, choked ratio

discharge_coefficient, carryover_factor, mdot (initial)

Base pressure

_solve_pressure_distribution()

compressible flow through each throttle, mass-flux matching

p[], t[], rho[], w[], mdot (final)

Swirl

_solve_swirl_velocities()

tangential momentum balance with wall shear

v[], taur[], taus[], cg, cx

Coefficients

_solve_perturbation_system()

perturb continuity + momentum, LU solve, integrate forces

kxx, kxy, kyx, cxx, cxy, cyx

2. Building a seal#

We use the same configuration as the test suite (ross/tests/test_labyrinth.py). With gas_composition the thermodynamic properties are derived automatically through ccp.

import plotly.io as pio

pio.renderers.default = "notebook"

import numpy as np
from ross.seals.labyrinth_seal import LabyrinthSeal
from ross.units import Q_

seal = LabyrinthSeal(
    n=0,
    shaft_diameter=Q_(145, "mm"),
    radial_clearance=Q_(0.3, "mm"),
    n_teeth=16,
    pitch=Q_(3.175, "mm"),
    tooth_height=Q_(3.175, "mm"),
    tooth_width=Q_(0.1524, "mm"),
    seal_type="inter",
    inlet_pressure=308_000.0,
    outlet_pressure=94_300.0,
    inlet_temperature=283.15,
    speed=Q_([8000], "RPM"),
    preswirl=0.98,
    gas_composition={"Nitrogen": 0.7812, "Oxygen": 0.2096, "Argon": 0.0092},
)
seal

/home/raphaelts/ross/.venv/lib/python3.14/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

LabyrinthSeal(n=0, n_link=None,
 kxx=[-50410.238024114755], kxy=[35509.21176636041],
 kyx=[-35509.21176636041], kyy=[-50410.238024114755],
 kzz=[0.0], cxx=[23.798628094904444],
 cxy=[56.20254034252011], cyx=[-56.20254034252011],
 cyy=[23.798628094904444], czz=[0.0],
 mxx=[0.0], mxy=[0.0],
 myx=[0.0], myy=[0.0],
 mzz=[0.0],
 speed=[837.75804096], frequency=None, tag=None)

2.1 Gas properties resolved in __init__#

When gas_composition is supplied, the constructor calls ccp.State.define(...) to get real-gas properties and stores the values the solver actually needs. Everything downstream uses only R, gamma, reference_temperatures, reference_viscosities — so the gas_composition and the “manual” paths are interchangeable.

print(f"molar mass : {seal.molar_mass:.5f} g/mol")
print(f"R          : {seal.R:.3f} J/(kg·K)   (= 8314 / molar mass)")
print(f"gamma      : {seal.gamma:.5f}")
print(f"reference temperatures [K]    : {seal.reference_temperatures}")
print(f"reference viscosities  [Pa·s] : {seal.reference_viscosities}")

molar mass : 28.95860 g/mol
R          : 287.100 J/(kg·K)   (= 8314 / molar mass)
gamma      : 1.40652
reference temperatures [K]    : [283.15, 282.6082973077877]
reference viscosities  [Pa·s] : [1.7788536533753193e-05, 1.7727919920975653e-05]

2.2 Geometry and the solver#

The element stores the geometry as plain scalars and precomputes the axial coordinate z of each station. The physics lives in seal.solver, a LabyrinthSolver sized from the geometry: n_stations = n_teeth + 1 stations, n_cavities = n_teeth − 1 interior cavities carrying dynamic unknowns, and ndof = 8 · n_cavities perturbation degrees of freedom.

solver = seal.solver

print(f"radial_clearance     : {seal.radial_clearance}  m")
print(f"pitch                : {seal.pitch}  m")
print(f"axial position z[:5] : {np.round(seal.z[:5], 5)}  m")
print(f"n_stations  (teeth + 1)  : {solver.n_stations}")
print(f"n_cavities  (teeth - 1)  : {solver.n_cavities}")
print(f"ndof (= 8 · n_cavities)  : {solver.ndof}")

radial_clearance     : 0.0003  m
pitch                : 0.003175  m
axial position z[:5] : [0.      0.00318 0.00635 0.00952 0.0127 ]  m
n_stations  (teeth + 1)  : 17
n_cavities  (teeth - 1)  : 15
ndof (= 8 · n_cavities)  : 120

3. Walking through solve() step by step#

The constructor already executed the whole pipeline. To inspect the intermediate state, we will re-drive the solver one stage at a time for a single rotor speed (with synchronous whirl) — exactly mirroring the body of LabyrinthSolver.solve(). The element itself is untouched by this: all mutable flow state lives on seal.solver.

omega = Q_(8000, "RPM").to("rad/s").m
solver.speed = omega  # rotor speed: sets the base flow
solver.whirl_frequency = (
    omega  # whirl frequency: sets the perturbation (synchronous here)
)
solver.inlet_swirl_velocity = solver.preswirl * omega * solver.shaft_radius

print(f"Ω (operating speed)   : {omega:.2f} rad/s")
print(
    f"inlet swirl velocity  : {solver.inlet_swirl_velocity:.3f} m/s "
    f"(= preswirl · Ω · shaft radius)"
)

Ω (operating speed)   : 837.76 rad/s
inlet swirl velocity  : 59.523 m/s (= preswirl · Ω · shaft radius)

3.1 _reset_state() — a clean slate#

_reset_state() re-creates the per-run state arrays (sized n_stations = n_teeth + 1) and the run constants:

  • pr, p, w, v, rho, t — pressure ratios and the base-flow fields per station; taur / taus the wall shear stresses; cg / cx the gradient tables filled by the swirl stage.

  • overall_pressure_ratio = outlet_pressure / inlet_pressure.

  • The perturbation amplitudes are constants of the solver: perturbation_eccentricity = 0.6 and the whirl displacement amplitudes pert_amplitude_direct / pert_amplitude_cross derived from it.

solver._reset_state()

print(f"perturbation eccentricity : {solver.perturbation_eccentricity}")
print(f"direct whirl amplitude    : {solver.pert_amplitude_direct:.3e} m")
print(f"n_cavities (interior)     : {solver.n_cavities}")
print(f"n_stations                : {solver.n_stations}")
print(f"ndof (= 8·n_cavities)     : {solver.ndof}")
print(f"overall pressure ratio    : {solver.overall_pressure_ratio:.4f}")

perturbation eccentricity : 0.6
direct whirl amplitude    : 1.800e-04 m
n_cavities (interior)     : 15
n_stations                : 17
ndof (= 8·n_cavities)     : 120
overall pressure ratio    : 0.3062

3.2 _vermes_leakage() — leakage seed and discharge model#

_vermes_leakage() produces a first estimate of the mass flux using the Vermes leakage correlation. It computes:

  • discharge_coefficient — from the ratio of tooth (tip) width to clearance.

  • carryover_factor — the kinetic-energy carry-over factor between cavities (forced to 0 for interlocking "inter" seals, where carry-over is suppressed).

  • _solve_choked_flow_function() — finds, with scipy.optimize.brentq, the choked pressure ratio that maximizes the Vermes flow function, and evaluates the flow function at the operating pressure ratio (or at the choked ratio when the seal is choked).

  • mdot — the resulting initial mass-flux guess (per unit circumference).

This mdot is only a starting point; the pressure stage refines it next.

solver._vermes_leakage()

print(f"discharge coefficient : {solver.discharge_coefficient:.5f}")
print(f"carry-over factor     : {solver.carryover_factor}      (0 for 'inter' seals)")
print(f"choked pressure ratio : {solver.choked_pressure_ratio:.5f}")
print(f"mdot (initial guess)  : {solver.mdot:.6f} kg/s per unit circumference")

discharge coefficient : 0.65382
carry-over factor     : 0      (0 for 'inter' seals)
choked pressure ratio : 0.03101
mdot (initial guess)  : 0.048660 kg/s per unit circumference

3.3 _solve_pressure_distribution() — base pressure and mass-flux matching#

_solve_pressure_distribution() solves the steady compressible flow. There are two nested solves:

  1. Per throttle. For each tooth, _throttle_pressure_ratio() finds the pressure ratio that passes the current mdot: the root of the throttle mass-flux balance, bracketed between the critical (choked) ratio and 1 and solved with scipy.optimize.brentq. If the residual at the critical ratio shows the tooth cannot pass mdot, the tooth is choked and the outer loop reacts. From the solved ratio, the code propagates pressure p, throat velocity w, density rho and temperature t (isentropic relations) to the next station.

  2. Globally. The tooth chain is repeated while bisecting mdot until the computed outlet pressure matches the prescribed outlet_pressure (or the last throttle chokes). This enforces overall mass conservation.

All thermodynamic relations go through the gas backend (IdealGas or RealGas), so the same solver handles both models.

solver._solve_pressure_distribution()

print(f"converged mdot (leakage) : {solver.mdot:.6f} kg/s per unit circumference\n")
print("station |   p [kPa] |   T [K] | rho [kg/m3]")
for i in range(solver.n_stations):
    print(
        f"  {i:2d}    | {solver.p[i] / 1e3:8.2f} | {solver.t[i]:7.3f} | {solver.rho[i]:8.4f}"
    )

converged mdot (leakage) : 0.051931 kg/s per unit circumference

station |   p [kPa] |   T [K] | rho [kg/m3]
   0    |   308.00 | 283.150 |   3.7888
   1    |   298.43 | 280.579 |   3.7047
   2    |   288.63 | 277.884 |   3.6178
   3    |   278.57 | 275.048 |   3.5277
   4    |   268.22 | 272.056 |   3.4340
   5    |   257.56 | 268.887 |   3.3364
   6    |   246.56 | 265.514 |   3.2344
   7    |   235.16 | 261.907 |   3.1274
   8    |   223.32 | 258.024 |   3.0146
   9    |   210.96 | 253.814 |   2.8950
  10    |   198.01 | 249.207 |   2.7675
  11    |   184.34 | 244.107 |   2.6302
  12    |   169.79 | 238.375 |   2.4809
  13    |   154.12 | 231.796 |   2.3158
  14    |   136.94 | 224.013 |   2.1292
  15    |   117.55 | 214.344 |   1.9102
  16    |    94.30 | 201.116 |   1.6332

The pressure decays in steps from inlet to outlet — one drop per tooth. The element exposes plot_pressure_distribution() (shared by all flow-model seals via SealElement), which reads the per-frequency profiles stored at construction:

fig_p = seal.plot_pressure_distribution(pressure_units="kPa", length_units="mm")
fig_p.show()

3.4 _solve_swirl_velocities() — cavity swirl and the gradient tables#

_solve_swirl_velocities() solves the tangential momentum balance in each cavity to get the swirl velocity v[i]. The driving physics:

  • Wall shear stresses on the rotor (taur) and stator (taus) using a Blasius-type friction law, tau = ½ ρ V² · 0.079 · Re^(−0.25), evaluated with the cavity hydraulic diameter.

  • Wetted-area ratios that depend on seal_type ("rotor", "stator", "inter").

  • The inlet pre-swirl v[0] = inlet_swirl_velocity sets the boundary condition; each cavity’s swirl is the root of the momentum residual, bracketed by the local sonic speed and solved with scipy.optimize.brentq.

  • With use_jenny_kanki=True, the Jenny–Kanki momentum factor (looked up per seal type) reduces the fraction of the through-flow momentum exchanged in each cavity; the classic model corresponds to a factor of 1.

Crucially, this stage also fills the gradient tables cg (9 rows) and cx (8 rows) — the partial derivatives of the continuity and momentum equations with respect to the field variables. They are precisely what the perturbation stage assembles into its linear system.

solver._solve_swirl_velocities()

print("cavity |  swirl v [m/s] |  taur [Pa] |  taus [Pa]")
for i in range(1, 7):
    print(
        f"  {i:2d}   | {solver.v[i]:13.3f} | {solver.taur[i]:10.3f} | {solver.taus[i]:10.3f}"
    )

print(
    f"\ncg column for cavity 1 (9 base-flow gradients):\n{np.round(solver.cg[:, 1], 4)}"
)
print(
    f"\ncx column for cavity 1 (8 momentum gradients):\n{np.round(solver.cx[:, 1], 4)}"
)

cavity |  swirl v [m/s] |  taur [Pa] |  taus [Pa]
   1   |        55.643 |      0.491 |     32.223
   2   |        52.323 |      1.156 |     28.303
   3   |        49.486 |      1.877 |     25.077
   4   |        47.062 |      2.576 |     22.399
   5   |        44.990 |      3.210 |     20.155
   6   |        43.218 |      3.759 |     18.252

cg column for cavity 1 (9 base-flow gradients):
[ 0.0000e+00  0.0000e+00  6.0000e-04  0.0000e+00 -0.0000e+00 -1.1800e-02
 -9.0276e+00 -0.0000e+00 -0.0000e+00]

cx column for cavity 1 (8 momentum gradients):
[ 2.000000e-04  0.000000e+00  3.140000e-02  5.940000e-02 -5.190000e-02
  0.000000e+00  0.000000e+00  6.785897e+02]

3.5 _solve_perturbation_system() — linear solve and the coefficients#

This is where the rotordynamic coefficients are born. The rotor is given a small whirl orbit, the continuity and momentum equations are linearized about the base flow, and the resulting linear system is solved. Step by step:

  1. _assemble_perturbation_system() places 8 unknowns per interior cavity — the cosine and sine components of the pressure (DOFs 0–3) and swirl velocity (DOFs 4–7) perturbations for the two whirl directions — into a dense (ndof, ndof) matrix. Rows 0–3 are the linearized continuity equations and rows 4–7 the tangential momentum equations; the off-diagonal blocks couple neighbouring cavities through the cg/cx tables. The right-hand sides for the direct and cross whirl perturbations are built from the imposed clearance perturbation amplitudes.

  2. LU-factor and solve. scipy.linalg.lu_factor / lu_solve solve both right-hand sides in a single call. The condition number is checked first: a near-singular system raises a ValueError; a high condition number warns.

  3. Integrate the perturbed pressures around the circumference to get the net forces, summed over all cavities and scaled by π · R · pitch. This yields kxx, kxy (with kyx = −kxy) and, when the whirl frequency ω ≠ 0, cxx, cxy (with cyx = −cxy, and a 1/ω factor). At ω = 0 the damping terms are set to zero. The whirl frequency enters only here: the unsteady terms of the perturbation matrix (cf1, cf4, cf7, cf9) and of the forcing scale with ω, while the base flow behind the cg/cx tables is set by the rotor speed Ω (see section 6).

The model is isotropic, so the symmetry kyy = kxx, cyy = cxx, kyx = −kxy, cyx = −cxy holds.

solver._solve_perturbation_system()

print(f"condition number   : {solver.pert_condition_number:.3e}")
print(f"rcond              : {solver.pert_rcond:.3e}\n")
print(f"kxx = kyy : {solver.kxx:12.3f} N/m")
print(f"kxy       : {solver.kxy:12.3f} N/m   (kyx = -kxy = {solver.kyx:.3f})")
print(f"cxx = cyy : {solver.cxx:12.4f} N·s/m")
print(f"cxy       : {solver.cxy:12.4f} N·s/m   (cyx = -cxy = {solver.cyx:.4f})")

condition number   : 1.326e+05
rcond              : 7.541e-06

kxx = kyy :   -50410.238 N/m
kxy       :    35509.212 N/m   (kyx = -kxy = -35509.212)
cxx = cyy :      23.7986 N·s/m
cxy       :      56.2025 N·s/m   (cyx = -cxy = -56.2025)

A few physical remarks:

  • kxx < 0 (negative direct stiffness) is typical for labyrinth seals — gas seals are not load-carrying like bearings.

  • kxy > 0 is the cross-coupled stiffness: driven by the swirl, it produces a follower force that can destabilize the rotor. Reducing inlet pre-swirl (e.g. with swirl brakes) lowers kxy — the whole motivation for modelling it.

  • The cross-coupling and the symmetry kyx = −kxy come straight out of the perturbation integration.

4. From solver output to a ROSS element#

solve() packs the results into a dictionary, mapping the isotropic solver values to the coefficient names BearingElement expects (note kyy = kxx, cyy = cxx):

return {
    "kxx": self.kxx, "kyy": self.kxx, "kxy": self.kxy, "kyx": self.kyx,
    "cxx": self.cxx, "cyy": self.cxx, "cxy": self.cxy, "cyx": self.cyx,
    "pressure": self.p,
    "seal_leakage": self._circumferential_leakage(self.mdot),
    ...
}

Back in LabyrinthSeal.__init__, these dicts (one per speed) are stacked into lists and passed to SealElement.__init__, which stores them as speed-aligned tables and builds interpolators. From then on the object behaves like any ROSS element.

# K() and C() use the coefficient interpolators (BearingCoefficient) built at construction by
# SealElement.__init__,
# so they return the stored coefficients regardless of the manual re-drive above.
spd = Q_(8000, "RPM").to("rad/s").m
print(f"K({spd:.1f} rad/s) =\n{seal.K(spd)}")
print(f"\nC({spd:.1f} rad/s) =\n{seal.C(spd)}")

K(837.8 rad/s) =
[[-50410.23802411  35509.21176636      0.        ]
 [-35509.21176636 -50410.23802411      0.        ]
 [     0.              0.              0.        ]]

C(837.8 rad/s) =
[[ 23.79862809  56.20254034   0.        ]
 [-56.20254034  23.79862809   0.        ]
 [  0.           0.           0.        ]]

Note: re-driving the solver by hand, as we did here, only touches seal.solver — the element’s stored coefficient lists are unaffected. In a normally constructed seal these are lists aligned with speed, e.g. seal.kxx[0], and seal.plot(coefficients=[...]) draws them against rotor speed. See tutorial_seal for the usage-oriented view, including attaching the seal to a Rotor.

5. Leakage#

The converged leakage is the base-flow mdot, exposed on the element as seal_leakage (a list aligned with speed):

print(f"leakage (mdot per unit circumference): {seal.seal_leakage[0]:.6f} kg/s")

leakage (mdot per unit circumference): 0.023705 kg/s

Coefficients are always computed#

There is no leakage-only shortcut: solve() always runs the perturbation stage, so stiffness, damping and seal_leakage are all available after construction. This keeps the model to a single, well-defined behavior.

Be aware that the perturbation stage is the most expensive one — for each whirl frequency it assembles and LU-solves the 8 · n_cavities perturbation system. The base-flow stages are comparatively cheap, which is why the element parallelizes the per-speed solves across processes when there are more than four frequencies.

6. Rotor speed versus whirl frequency#

Two different frequencies enter the model, and the pipeline above keeps them apart:

  • The rotor speed Ω only appears in the base flow: the inlet swirl preswirl · Ω · R (_solve_base_flow) and the rotor surface velocity R · Ω that drives the wall shear in _solve_swirl_velocities. Everything the base flow produces — leakage, cavity pressures, swirl velocities and the cg/cx gradient tables — depends on Ω alone.

  • The whirl frequency ω only appears in the perturbation: the unsteady coefficients of the linearized continuity and momentum equations (cf1 = ω·cg0 + cg1, cf4 = −ω·cg0 + cg1, cf7 = ω·cx1 + cx2, cf9 = −ω·cx1 + cx2), the forcing terms of _assemble_perturbation_system, and the damping extraction c = f / ω in _solve_perturbation_system.

solve(speed, frequency=None) therefore runs _solve_base_flow(speed) and then _solve_whirl(frequency), with frequency defaulting to speed — the synchronous whirl used for the 1-D tables. Because the expensive part (the bracketed root-finding of the base flow) does not depend on ω, one base flow can serve several whirl frequencies:

  • solve_row(speed, frequencies) solves the base flow once and the perturbation once per whirl frequency;

  • solve_grid(speeds, frequencies) maps solve_row over the speeds (in parallel when there are more than parallel_threshold), which is what LabyrinthSeal(speed=..., frequency=...) uses to build a 2-D (speed, frequency) coefficient table.

The diagonal of that grid — whirl frequency equal to the rotor speed — is exactly the synchronous solve:

speeds = Q_([5000, 8000], "RPM").to("rad/s").m
whirls = np.sort(np.concatenate([0.5 * speeds, speeds]))

rows = solver.solve_grid(speeds, whirls)  # rows[i][j]: speeds[i], whirls[j]
synchronous = [solver.solve(w) for w in speeds]

for i, w in enumerate(speeds):
    j = int(np.where(whirls == w)[0][0])
    print(
        f"Ω = {w:7.2f} rad/s | grid diagonal kxy = {rows[i][j]['kxy']:12.3f}"
        f" | synchronous kxy = {synchronous[i]['kxy']:12.3f}"
        f" | equal: {rows[i][j]['kxy'] == synchronous[i]['kxy']}"
    )

print("\nkxy at Ω = 8000 RPM for whirl / Ω =", np.round(whirls / speeds[1], 3))
print([round(float(r["kxy"]), 1) for r in rows[1]])
print(
    "leakage is a base-flow quantity, identical along the row:",
    len({r["seal_leakage"] for r in rows[1]}) == 1,
)

Ω =  523.60 rad/s | grid diagonal kxy =    16600.326 | synchronous kxy =    16600.326 | equal: True
Ω =  837.76 rad/s | grid diagonal kxy =    35509.212 | synchronous kxy =    35509.212 | equal: True

kxy at Ω = 8000 RPM for whirl / Ω = [0.312 0.5   0.625 1.   ]
[29121.7, 30108.7, 31050.6, 35509.2]
leakage is a base-flow quantity, identical along the row: True

On the element, the 2-D table is exposed through the same interface as the 1-D one: seal_2d.kxy is a nested list of shape (len(speed), len(frequency)), seal_2d.kxy_interpolated has kind == "grid" and evaluates (frequency, speed) pairs, and a single value still returns the synchronous diagonal. Rotor.run_modal(speed, frequency=...) and Rotor.run_modal(speed, matched_whirl=True) are the analyses that put the whirl axis to use; see tutorial_seal section 4.

7. Summary#

  • The physics and the element are separate objects. LabyrinthSolver computes; LabyrinthSeal wraps it and turns the results into a standard speed- (and whirl-frequency-) dependent coefficient element.

  • solve() is a fixed physics pipeline: _reset_state → _vermes_leakage → _solve_pressure_distribution → _solve_swirl_velocities → _solve_perturbation_system, executed once per speed (in parallel when there are many); the perturbation stage can be repeated at several whirl frequencies on the same base flow (solve_row, solve_grid).

  • Base flow first, perturbation second. The leakage, pressure and swirl fields are solved with bracketed root-finding (scipy.optimize.brentq); the perturbation stage linearizes them to get stiffness and damping.

  • The cg/cx tables are the bridge between the base-flow solution and the perturbation matrix.

  • Two frequencies, two roles: the rotor speed Ω drives the base flow (swirl, surface velocity), the whirl frequency ω drives the unsteady perturbation terms and the damping.

  • Symmetry and signs are physical: isotropy gives kyy = kxx, cyy = cxx; the cross-coupled kxy (swirl-driven) is the destabilizing term that this whole model exists to quantify.

Contents