Multi-DOF Vibration and Critical Speed: Method

Multi-DOF Vibration and Critical Speed: Method

1. The lumped-parameter chain

The chain is n masses m_i in a line. Row i carries the spring k_i and the dashpot c_i that join mass i to mass i−1; with the left end fixed, row 1 joins mass 1 to the wall, with it free, row 1's spring is ignored. An optional spring k_R joins the last mass to the right wall. The equations of motion are M ẍ + C ẋ + K x = F(t) with the diagonal mass matrix M = diag(m_i) and the tridiagonal K (k_i + k_i+1 on the diagonal, −k_i+1 beside it); C has the same pattern with the c_i.

Eigen-solve. The undamped free vibration x = φ sin ωt gives K φ = ω² M φ. With M diagonal the substitution v = M^½ φ turns it into the symmetric problem (M^−½ K M^−½) v = ω² v, which the tool solves by Householder reduction to tridiagonal form followed by the QL algorithm with implicit shifts (the tred2 / tqli pair of Press et al.), giving every eigenvalue and eigenvector at once to round-off. Each mode is normalised to its largest component. A system with no connection to ground has a rigid-body mode at ω = 0; the elastic modes are numbered from 1 above it.

Modal quantities. Modal mass m_r = φ_r^T M φ_r, modal stiffness k_r = ω_r² m_r. For a uniform motion of the base (every mass driven through its spring to the wall by the same ground acceleration) the modal force is −(φ_r^T M 1) ÿ_g, so the participation factor is Γ_r = φ_r^T M 1 / m_r and the effective modal mass Γ_r² m_r is the part of the total mass that mode r contributes to the base reaction; the effective masses of all modes sum to the total mass, which the table shows as shares.

Holzer's tabulation. At a trial frequency ω set x₁ = 1 and carry the displacement and the spring force along the chain: the force in the spring left of mass i is F_i, the force to its right is F_i+1 = F_i − m_i ω² x_i, and the next mass sits at x_i+1 = x_i + F_i+1/k_i+1. F₁ is k₁ x₁ for a fixed left end and 0 for a free one. At a natural frequency the condition at the right end is met exactly: F_n+1 = 0 for a free end, F_n+1 + k_R x_n = 0 with an end spring. The table is printed at the eigen-solve frequency of the chosen mode, so its residual (normalised by the largest force met) is a check of the two methods against each other, and the residual is also drawn over a sweep of trial frequencies, where its zeros are the natural frequencies; the x column is the mode shape in Holzer's normalisation x₁ = 1.

Rayleigh's quotient and Dunkerley's bound. Rayleigh's quotient ω² = φ^T K φ / φ^T M φ is an upper bound on ω₁² for any trial shape φ. The tool uses the static deflection y of the chain under the loads m_i g (the classical choice), for which φ^T K φ = Σ m_i g y_i and ω_R² = g Σ m_i y_i / Σ m_i y_i². Dunkerley's equation 1/ω_D² = Σ m_i δ_ii, with δ_ii the flexibility at mass i (the diagonal of K⁻¹, i.e. the deflection at i under a unit force at i), equals the trace of M K⁻¹, which is the sum of all 1/ω_r² and therefore exceeds 1/ω₁²: ω_D ≤ ω₁ ≤ ω_R. Each term m_i δ_ii is the inverse square frequency of mass i alone on the chain. Both need a supported chain (K invertible); a free-free chain has neither a static deflection nor a flexibility.

Harmonic response. For F(t) = F₀ sin Ωt on one mass the steady state is x = Re(X e^iΩt) with the complex amplitudes from (K − Ω² M + iΩ C) X = F, solved directly as a real system of twice the size. The same response by mode superposition is X = Σ_r φ_r (φ_r^T F) / [m_r (ω_r² − Ω² + 2iζ_rω_rΩ)], each term a single-degree-of-freedom response in the modal coordinate q_r. The damping enters in one of two ways. With every dashpot at zero the modal damping ratio ζ is given to every mode and the direct solve uses the matrix C = Σ_r 2ζω_r (Mφ_r)(Mφ_r)^T/m_r it implies, so the two solutions agree to round-off. With dashpots entered the direct solve uses them, and the modal solution uses the ratios ζ_r = φ_r^T C φ_r / (2ω_rm_r) they imply with the coupling terms φ_r^T C φ_s dropped; the two then differ where the damping is not proportional, and the table shows by how much. Phase lags are of the displacement behind the force, 0 to 360°. The receptance chart sweeps |X_i|/F₀ from 0 to 1.3 times the highest natural frequency on a logarithmic scale.

Transient by Newmark-β. The equations are integrated with the average-acceleration (trapezoidal) form, β = ¼ and γ = ½: unconditionally stable, no numerical damping, period error (ω Δt)²/12 at the highest mode, with Δt a fortieth of the highest-mode period. The free vibration from an initial displacement of the force mass runs over four periods of the first elastic mode. For the half-sine base shock the left wall moves with ÿ_g(t) = A sin(πt/τ) for 0 ≤ t ≤ τ; in the relative coordinates z = x − y_g the chain obeys M z̈ + C ż + K z = −M 1 ÿ_g(t), so the inertia loading acts on every mass and the record runs over the pulse plus four first-mode periods. The largest relative displacements are the spring strokes.

2. The torsional system

Discs with polar inertias J_i on shaft segments of torsional stiffness k_t,i = G π d_i⁴ / (32 L_i) (solid round segment; the entered stiffness may be used instead) form the same chain with J for m, k_t for k and the twist θ for x, so the eigen-solve, Holzer's table (torque T in place of force) and the Rayleigh and Dunkerley bounds are those of section 1. A rotating train is free at both ends and has a rigid-body mode at zero frequency; a disc held to a frame makes a fixed end.

Gears. A branch turning at n times the speed of the reference shaft (shaft 1) stores kinetic energy ½ J (n θ̇)² and strain energy ½ k_t (n Δθ)², so it is referred to the reference shaft by J′ = n² J and k_t′ = n² k_t; the tool applies the row's ratio n_i to the disc and to the segment on its side of the gear. Single-path trains only; a branch that splits into two paths is not covered.

Nodes and critical speeds. The twist is linear along a massless segment, so a node (a section that does not turn) lies where θ changes sign between two discs, at the fraction θ_i/(θ_i − θ_i+1) of the segment from disc i; mode r has r nodes. An excitation of order q (q cycles per revolution of the reference shaft, half orders for a four-stroke engine, the blade or tooth count for a propeller, pump or gear) at speed N rpm has the frequency q N/60 Hz and meets mode r at the critical speed N_cr = 60 f_r/q. The Campbell diagram draws the order lines q N/60 and the horizontal natural frequencies; their crossings are the critical speeds and the operating speed is a vertical marker. The separation of the operating speed from the nearest critical of each order is tabulated; what separation a design needs depends on the strength of that order's excitation and the damping, which the tool does not assess.

3. Shaft critical speeds

The lateral critical speeds of a shaft are its lateral natural frequencies, at which the synchronous unbalance excitation resonates. The tool takes a stepped round shaft (solid or hollow sections), discs given by their mass, and bearings that are rigid (deflection zero) or linear springs to ground, and ignores the gyroscopic moments of the discs and shear deflection: the results are the classical critical speeds of the non-rotating shaft, which coincide with the forward whirl of a thin-disc rotor only when the gyroscopic stiffening is small.

Lumped-station model. The shaft is divided into at least the entered number of elements with a station at every section boundary, disc and bearing. The mass of each element (ρ A l) is split equally to its two stations and the discs are added at theirs; between stations the shaft is a massless Euler beam. The stiffness matrix is assembled from the exact massless-beam (Hermite) element on the deflection and rotation of each station, the bearing springs are added to the deflections, the deflections at rigid bearings are removed, and the rotations, which carry no mass, are condensed statically, K_c = K_yy − K_yθ K_θθ⁻¹ K_θy; this is exact for the lumped model, so the only approximation is the lumping, which converges from below (60 elements give the first critical of a uniform shaft to 0.001 %). The natural frequencies are the eigenvalues of M^½ F M^½ with F = K_c⁻¹ the flexibility (ω = 1/√λ), solved by the same Householder–QL routine; the flexibility route condenses light stations exactly and gives the influence coefficients δ_ij between the discs directly.

Transfer-matrix pass (Myklestad–Prohl). The state (y, θ, M, V) with M = EI y″ and V = M′ is carried from the free left end across each massless element, y⁺ = y + θl + Ml²/(2EI) + Vl³/(6EI), θ⁺ = θ + Ml/EI + Vl²/(2EI), M⁺ = M + Vl, V⁺ = V, and across each station, V⁺ = V + m ω² y for the station mass and V⁺ = V − K y for a bearing spring (NASA TN D-7385, eqs. 29 and 30). The free left end leaves y₀ and θ₀ unknown, and each rigid bearing adds an unknown reaction R_j (a jump in V) together with the condition y = 0 there; with the free-end conditions M = V = 0 on the right this gives 2 + n_b homogeneous equations in 2 + n_b unknowns, and the critical speeds are the zeros of their determinant, swept over speed and refined by bisection. Carrying the reactions as extra basis columns is Prohl's treatment of intermediate rigid supports and keeps the determinant a smooth function of ω. The pass runs on the same stations as the eigen-solve, so the two must agree to the tolerance of the sweep; the table prints both, and the mode shape of the first critical by the transfer matrix is drawn beside the eigenvector.

Rayleigh and Dunkerley. Rayleigh's method takes the static deflection curve under the station weights W_i = m_i g (shaft and discs, bearing springs compressed) as the trial shape, ω_R² = g Σ W_i y_i / Σ W_i y_i², an upper bound; for a uniform shaft on rigid end bearings the static shape gives ω = √(3024/31) √(EI/μL⁴), 0.072 % above the exact π² √(EI/μL⁴). Dunkerley's equation 1/ω_D² = Σ m_i δ_ii is summed over every station (the trace of M F), a lower bound; the classical form by parts, each disc alone on the massless shaft plus the bare shaft's own first critical (found by a second solve with the discs removed), is printed too. The influence coefficients at the discs are the sub-matrix of F at the disc stations, symmetric by Maxwell's reciprocity.

Campbell line. With only synchronous excitation the Campbell diagram reduces to the line f = N/60 against the horizontal criticals; the crossings are the critical speeds and the operating speed is marked. Operating within 25 % of a critical is noted; the margin a rotor needs depends on its unbalance and damping.

4. Checks

The test script (test-scripts/test_mdof_critical_speed.php, 308 checks) compares the library with closed forms and with an independent numpy / scipy implementation (References/MDOF-Vibration/checks/mdof_check.py): the uniform fixed-free chain ω_r = 2√(k/m) sin((2r−1)π/(2(2n+1))) for n up to 20, the two-degree-of-freedom quadratic and its mode ratio, twenty seeded chains with every end condition against numpy.linalg.eigh, Holzer residuals of zero at every eigenvalue, the ordering ω_D ≤ ω₁ ≤ ω_R, the harmonic response against a direct complex solve and the mode superposition, Newmark against the modal closed form and scipy odeint, the two-disc torsional closed form and the node at J₂/(J₁ + J₂), a geared train, NACA TN 1522's lumped torsion example (ω² = 4985 for the ten-station lumped model against the report's 5050 by Stodola iteration on the distributed beam), the uniform shaft (nπ/L)² √(EI/ρA), a disc on a massless shaft (48EI/mL³ and a²b²/3EIL), two discs by the 2 × 2 flexibility problem, a rigid shaft on two springs √(2k/m), five seeded stepped shafts with flexible and rigid bearings and overhangs against a numpy Hermite model, and the transfer-matrix pass against the eigen-solve on every shaft.