Random Vibration Response: Method

Random Vibration Response: PSD Integration, Miles' Equation, the Exact SDOF Response, Crossings and Peaks

The problem

Equipment on a launch vehicle, an aircraft, a truck or a shaker table sees a vibration that is random: no single frequency and no repeating waveform, only statistics. A stationary Gaussian random vibration is described by its one-sided power spectral density W(f), the mean-square acceleration per hertz in g²/Hz against frequency in Hz, whose integral over frequency is the mean-square acceleration; the overall level is g_rms = √(∫W df). Specifications such as the GEVS of GSFC-STD-7000, NASA-HDBK-7005 and MIL-STD-810H Method 514.8 give W(f) as a few breakpoints joined by straight lines on log-log paper, with the overall g rms printed beside the curve. A component mounted on that base amplifies the input near its natural frequency. The design questions are how large the response is, how often it crosses a level, how large its peaks get in the test or flight duration, what deflection and inertia load that means for the mounting, and whether the stress stays below the allowable. This tool answers them for a component idealized as a single-degree-of-freedom system, in the frequency domain, without generating a time history.

Units and conventions

PSD levels are in g²/Hz and frequencies in Hz throughout; they are plain numbers and the unit switch leaves them alone. The component mass (lbm or kg), the stress per g and the allowable stresses (ksi or MPa) go through the unit system and are converted to SI base units for the calculation; stresses are reported in ksi (USCS) or MPa (SI), displacements in inches or millimetres, the static load in pounds-force or newtons. Durations are in seconds. The spectral moments are taken with frequency in Hz, m_n = ∫ fⁿ G(f) df, which gives rates directly per second. g = 9.80665 m/s² = 386.09 in/s².

Input PSD and its rms

Between two breakpoints (f₁, W₁) and (f₂, W₂) the specification is a straight line on log-log axes, W = W₁(f/f₁)^m with m = log(W₂/W₁)/log(f₂/f₁). The slope in dB per octave is 10 log₁₀(2^m) = 3.01 m: +3 dB/octave is m = 1 (W proportional to f), +6 dB/octave is m = 2 (W proportional to f²), −6 dB/octave is m = −2. A specification often gives the slope instead of the second level (the GEVS table reads "20 Hz 0.026, 20–50 Hz +6 dB/oct, 50–800 Hz 0.16"); the tool accepts either, a line such as 50 +6 dB/oct placing that breakpoint on the slope from the previous one, and the typed table has a slope column for the same purpose. The area of a segment is exact (Irvine 2000, eq. 11): W₁f₁[(f₂/f₁)^m+1 − 1]/(m + 1), or W₁f₁ ln(f₂/f₁) when m = −1, and the overall level is the square root of the sum. The standard profiles on the dropdown reproduce the rms values printed in their sources: GSFC-STD-7000A/B Table 2.4-3 (qualification 14.1 g rms, acceptance 10.0) and Table 2.4-4 (the component minimum workmanship level, 6.8 g rms, the same spectrum as NASA-HDBK-7005 Table 8.9), NASA-HDBK-7005 Table 8.10 (USAF units, 6.1 g rms), MIL-STD-810H Figure 514.8E-1 (minimum integrity, 0.04 g²/Hz 20–1000 Hz then −6 dB/octave, 7.7 g rms, one hour per axis), Table 514.8C-I (common-carrier truck, vertical, 1.08 g rms for 1,000 miles in 60 minutes), Table 514.8C-II (the envelope for unknown orientation, 1.17 g rms), Table 514.8C-IX (jet aircraft cargo general exposure, 4.02 g rms) and NAVMAT P-9492 Figure 5 (6.0 g rms printed, 6.06 exact). A flat PSD of level W between f₁ and f₂ has rms √(W(f₂ − f₁)): the 0.01 g²/Hz floor of Table 2.4-4 from 20 to 2000 Hz gives the 4.4 g rms the table states.

Response of a single-degree-of-freedom component

A mass m on a spring k and damper c fixed to the vibrating base has the natural frequency f_n = (1/2π)√(k/m), the damping ratio ζ = c/(2√(km)) and the amplification Q = 1/(2ζ), the transmissibility at resonance. Its absolute acceleration x(t) is related to the base acceleration y(t) by the transmissibility |H(f)|² = [1 + (2ζr)²]/[(1 − r²)² + (2ζr)²], r = f/f_n, and because the system is linear the response PSD is G_out(f) = |H(f)|² W(f) (Crandall and Mark 1963; Bendat and Piersol 2010). The response is Gaussian like the input, so its rms is √(∫G_out df) and that one number, with the rates below, fixes every statistic.

Miles' equation. For a flat input W over all frequencies the integral is exact: ∫₀^∞|H|² df = (π/2) f_n Q (1 + 1/Q²). Miles (1954) dropped the small 1/Q² term and wrote g_rms = √[(π/2) f_n Q W(f_n)], the everyday estimate of a resonant response; NASA-HDBK-7005 (eq. 5.27) states it in the form σ² = (π/2) G₀ f_n Q and says it "is generally a good approximation if G₀ is taken to be the value at the natural frequency f_n when the random input is nonwhite". For Q = 10 the omitted term is 0.5 % of the rms. The tool prints Miles' equation with the substituted values and the exact white-noise value beside it.

The exact integral. The actual specification is neither flat nor unbounded, so the tool integrates |H|² W(f) over the profile's own band, splitting it at the breakpoints and at f_n(1 ± kζ) for k from ¼ to 32 and applying 10-point Gauss–Legendre quadrature on twelve logarithmically spaced pieces of each part; the result is converged to better than 1 part in 10⁸ (the test script compares it with an 800,000-point trapezoid integration in numpy and with the closed forms for a flat band of widening width). The ratio Miles / exact tells how far the hand estimate is off for the case: above one when the profile falls away from f_n or the band is cut close to it (the response of a component at 100 Hz to the GEVS level, which is flat from 50 to 800 Hz, is 15.74 g rms against Miles' 15.85), below one when the profile rises on either side of f_n. The cumulative-fraction chart shows where the mean-square response comes from; for a lightly damped component about half of it lies inside the half-power band f_n(1 ± 1/(2Q)) and the remainder is the quasi-static response to the input below the resonance and the ripple above it.

Spectral moments, crossings and peaks

Rice (1944, 1945) showed that a stationary Gaussian process with one-sided PSD G(f) crosses zero upward at the expected rate E[0] = √(m₂/m₀) and has peaks at the rate E[P] = √(m₄/m₂), with m_n = ∫ fⁿ G df, and that it crosses the level a upward at the rate ν_a = E[0] e^{−a²/(2σ²)}. The ratio γ = E[0]/E[P] is the irregularity factor, 1 for a narrow-band process and smaller when high-frequency ripple adds peaks that do not cross zero; for a resonant response under a broadband input the fourth moment is dominated by the ripple above the resonance, so E[P] and γ depend on the upper frequency limit of the specification, while E[0] stays close to f_n. The number of positive zero crossings in a duration T, E[0]T, is the number of response cycles. For the relative displacement z = x − y, whose transmissibility is 1/[(2πf_n)²(1 − r² + 2iζr)], the white-noise integrals ∫|H_z|² dr and ∫r²|H_z|² dr are both π/(4ζ), so its mean-square is (π/2) f_n Q W g²/(2πf_n)⁴ with no 1/Q² term and its zero-crossing rate is exactly f_n; the tool integrates it over the actual profile as well.

Design levels. The response is Gaussian: it lies within ±1σ 68.27 % of the time, within ±2σ 95.45 % and within ±3σ 99.73 %. The peaks of a narrow-band Gaussian process are Rayleigh distributed, P(peak > kσ) = e^−k²/2: 60.7 % of the peaks exceed 1σ, 13.5 % exceed 2σ, 1.11 % exceed 3σ; the mean peak is √(π/2) σ = 1.25σ and the level exceeded by 0.27 % of the peaks is 3.44σ, the "three-sigma peak value" of NASA-HDBK-7005 section 5.1.7. The 3σ level is the usual design level for random vibration (NASA-HDBK-7005 section 8 treats the 3σ response as a limit load); the tool takes the multiple k as an input and reports the design-level acceleration k G_rms, the relative displacement k Z_rms and its double amplitude, the inertia load P = m k G_rms g on the component mass for the mounting and fastener checks, and the stress k σ_1g G_rms.

Exceedance in a duration. The up-crossings of a high level are rare and nearly independent, so they are treated as a Poisson process (Rice; Davenport 1964): the probability of no exceedance of kσ in a duration T is e^−ν_kT, the probability of at least one exceedance 1 − e^−ν_kT, the mean first-passage time 1/ν_k, the time at which the exceedance probability reaches p is T_p = −ln(1 − p)/ν_k, and the level exceeded with probability p during T is k_p = √(2 ln[E[0]T/(−ln(1 − p))]). Davenport's expected largest peak in T is k_max = √(2 ln(E[0]T)) + 0.5772/√(2 ln(E[0]T)) in multiples of the rms. These statements are for positive peaks; for the magnitude |x| the rate doubles (2E[0]). For a two-minute GEVS qualification run at 100 Hz the response crosses 3σ about 135 times, so a 3σ design level is a statement about the stress at which the material spends 0.27 % of the time, not a level that is rarely reached; the level exceeded with 1 % probability in the run is about 5.3σ and the expected largest positive peak 4.5σ. The seeded time-domain simulation filed with the tool reproduces the rms, the zero-crossing rate and the 3σ crossing rate of the default case and gives a mean largest |x| of 4.56σ per 120-second window, between the one-sided and two-sided Davenport values.

Stress check

The stress at the critical point follows the response acceleration through a stress per g, σ_1g, obtained from a static 1 g analysis of the component on its mounting or a unit-acceleration load case in a finite-element model (for a mass on a member of stressed area A it is m g/A). The stress rms is σ_1g G_rms, the design-level stress k σ_1g G_rms, and the design-level stress is a limit stress: the kσ response of a random environment is the limit load of NASA practice (NASA-HDBK-7005 section 8), k being the peak factor (3 by default; some programmes use 4 or 4.5, or the level with a set exceedance probability). The tool therefore offers no load-level radio and writes the margins in the site's convention: strength margins against F_tu and F_cy at ultimate load, M.S. = F / (FS · k σ_1g G_rms) − 1 with FS the ultimate factor (default 1.5), labelled ULT; yield margins against F_ty and F_cy at limit load times the yield factor, M.S. = F / (FS_yld · k σ_1g G_rms) − 1 (default 1.15, the factor several agencies require on limit load for no detrimental permanent set), labelled LIM. The vibratory stress alternates and reaches both signs, so the compressive allowable is checked on the same magnitude, F_cy standing for the compression allowable of a stable element; a part that can buckle or cripple needs that allowable instead. Each box names the level of its allowable. Fatigue under the random stress usually governs before the peak stress does and is a separate question, answered by the Random Vibration Fatigue tool from the same PSD and component.

Limits

The model is linear, stationary and Gaussian with one degree of freedom. A component with several modes in the band needs a modal or finite-element random response (or the transmissibility of each mode treated in turn); non-Gaussian inputs (kurtosis above 3) have more large peaks than the Rayleigh distribution predicts; sine-on-random and non-stationary flight events are not represented by a single PSD; and the mass of the component loads the base, which force limiting (NASA-HDBK-7004) accounts for and a hard-mounted test does not. Miles' equation and the exact integral both assume the base motion is unaffected by the component. The exceedance statistics use the Poisson approximation, which is accurate for levels above about 2σ and conservative below.

Reading the results

The input section tabulates every segment with its slope in dB/octave and its exact area, and compares the overall g rms with the value printed in the source when a standard profile is used. The response section prints Miles' equation with the numbers, the exact white-noise value and the integrated value, the response PSD on log-log axes beside the input, the cumulative fraction of the mean-square response, the spectral moments, the crossing and peak rates and the number of cycles in the duration. The design-level table lists the 1σ, 2σ and 3σ acceleration, relative displacement, stress and static load. The exceedance section gives the 3σ crossing rate and count, the first-passage time, the level with the chosen exceedance probability and the expected maximum, with a chart of the three probability statements against the level. The stress check ends in a margin box, and the last table sets Miles' estimate against the exact integral for every reported quantity.