1. Background: Spectral Fatigue Analysis
Frequency-domain fatigue methods avoid the cost of time-domain rainflow counting by working directly with the one-sided stress PSD G(ƒ). The key quantities are the spectral moments:
m_n = ∫_0^∞ f^n G(f) df, n = 0, 1, 2, 4(1)From these moments several bandwidth descriptors are derived:
ν_0 = √(m_2/m_0) — zero-crossing rate (Hz) ν_P = √(m_4/m_2) — peak rate (Hz) α_1 = m_1 / √(m_0 m_2) — first bandwidth parameter α_2 = m_2 / √(m_0 m_4) — second bandwidth parameter (irregularity factor)(2)The RMS stress is σ = √m0. For a Gaussian process the peak amplitude PDF is fully characterised by σ and the bandwidth. The Dirlik method exploits these parameters to approximate the rainflow amplitude distribution without explicit time-domain simulation.
2. The Dirlik Rainflow Amplitude PDF
Dirlik (1985) fitted a three-term mixture density to Monte Carlo rainflow histograms generated from a wide range of PSDs:
p_DK(S) = [D_1/Q · exp(-Z/Q) + D_2 Z/R^2 · exp(-Z^2/2R^2) + D_3 Z · exp(-Z^2/2)] / (2√m_0)(3)where Z = S / (2√m0) is the normalised amplitude. The five Dirlik coefficients are:
x_m = m_1/m_0 · √(m_2/m_4) D_1 = 2(x_m – α_2^2) / (1 + α_2^2) R = (α_2 – x_m – D_1^2) / (1 – α_2 – D_1 + D_1^2) D_2 = (1 – α_2 – D_1 + D_1^2) / (1 – R) D_3 = 1 – D_1 – D_2 Q = 1.25(α_2 – D_3 – D_2 R) / D_1(4)2.1 Closed-Form Fatigue Damage Rate
For the Basquin S-N relation N Sb = C, the expected fatigue damage per unit time is:
Ḋ_DK = (ν_P / C) · (2√m_0)^b · [D_1 Q^b Γ(b+1) + (D_2 R^b + D_3) Γ((b+2)/2)](5)where Γ(·) is the complete Gamma function. Fatigue life is T = 1 / Ḋ.
3. Why Kurtosis Matters Even for Stationary Processes
A process is stationary if its statistical properties do not change with time—mean, variance, and higher moments are time-invariant. Stationarity does not imply Gaussianity. Many real-world stress histories are stationary yet leptokurtic (κ > 3):
- Aircraft structural loads during turbulence with intermittent gusts
- Road vehicle chassis loads on mixed surfaces
- Offshore structural responses in combined swell and wind seas
- High-cycle fatigue test rigs with occasional overloads
The PSD of a stationary non-Gaussian process is identical in form to that of a Gaussian process with the same variance. Therefore two signals can share the same PSD while having very different peak amplitude statistics. The Dirlik integral, which assumes Gaussian input, sees only the PSD and therefore underestimates damage when κ > 3.
3.1 Physical Mechanism
Excess kurtosis (κ − 3 > 0) indicates heavier-than-Gaussian tails in the amplitude PDF. In fatigue terms this means a larger fraction of cycles occur at high stress amplitudes. Because S-N damage scales as Sb with b = 3–8 for metals, even a modest shift of probability mass toward high amplitudes produces a disproportionate increase in cumulative damage.
For illustration, the probability that a stress amplitude exceeds 3σ in a Gaussian process is about 0.27%. For a process with κ = 5 this tail probability can be several times larger, depending on the specific distribution shape.
4. Correction Approaches
4.1 Method 1: Winterstein Hermite Polynomial Transformation
Winterstein (1988) proposed mapping the non-Gaussian process X(t) to an underlying Gaussian process U(t) via a monotone polynomial transformation:
X(t) = σ [U + h_3(U^2 – 1) + h_4(U^3 – 3U)](6)The coefficients h3 and h4 are estimated from the skewness γ1 and excess kurtosis (κ − 3) of the measured record:
h_3 ≈ γ_1 / 6 h_4 ≈ (κ – 3) / 24(7)The standard deviation of the underlying Gaussian is σU = (1 + 2h32 + 6h42)−1/2.
The Hermite model converts the stress PSD to an equivalent Gaussian PSD via a bandwidth-dependent scaling, then applies the standard Dirlik integral to the modified PSD. For zero-skewness processes (γ1 ≈ 0, so h3 ≈ 0) the correction is driven entirely by h4.
Applicability: Best for mild non-Gaussianity (3 ≤ κ ≤ 6) and broadband processes. Requires reliable estimates of γ1 and κ from measured data. Fails if the transformation is non-monotone, which occurs roughly when |h3| + 3|h4| > 1/3.
4.2 Method 2: Polynomial Kurtosis Correction Factor Γκ
Several authors (Benasciutti & Tovo 2006; Rizzi et al. 2011) have proposed a multiplicative correction to the Dirlik damage rate:
Ḋ_corrected = Γ_κ · Ḋ_DK(8)where Γκ is a polynomial in the excess kurtosis (κ − 3) calibrated against Monte Carlo rainflow results. A convenient two-parameter form is:
Γ_κ = 1 + c_1 (κ – 3) + c_2 (κ – 3)^2(9)Table 1 gives indicative values of c1 and c2 for common S-N slopes, based on Rizzi et al. (2011) for a moderately broadband process (α2 ≈ 0.6):
| S-N slope b | c1 | c2 | Γκ at κ = 5 |
|---|---|---|---|
| 3 | 0.045 | 0.003 | 1.102 |
| 4 | 0.062 | 0.005 | 1.144 |
| 5 | 0.078 | 0.007 | 1.184 |
| 6 | 0.095 | 0.010 | 1.230 |
| 7 | 0.112 | 0.013 | 1.276 |
| 8 | 0.130 | 0.017 | 1.328 |
Table 1. Polynomial correction factors Γκ for the Dirlik damage rate; α2 ≈ 0.6. Values adapted from Rizzi et al. (2011). Recalibrate for significantly different bandwidth.
Applicability: Straightforward to implement. Valid for 3 ≤ κ ≤ 8 and zero skewness. Coefficients should be recalibrated if the bandwidth parameter α2 differs substantially from the calibration dataset.
4.3 Method 3: Non-Gaussian Amplitude Distribution Substitution
A third approach replaces the Gaussian peak amplitude PDF implicit in Dirlik’s formulation with a heavier-tailed distribution whose kurtosis matches the measurement. Three options are commonly used:
- Generalised Gaussian distribution: Shape parameter p < 2 gives heavy tails; p = 2 recovers Gaussian. Kurtosis = Γ(5/p)Γ(1/p)/[Γ(3/p)]2.
- Student’s t-distribution: ν degrees of freedom; κ = 3 + 6/(ν − 4) for ν > 4. Analytically convenient.
- Johnson SU distribution: Four-parameter family; matches skewness and kurtosis simultaneously.
The damage integral is then evaluated numerically using the chosen distribution in place of the Rayleigh/Gaussian peak PDF. This is the most general method but requires numerical integration for each case.
5. Decision Guide: Which Method to Use?
| Kurtosis range | Recommended method | Notes |
|---|---|---|
| κ = 3 | Standard Dirlik | No correction needed |
| 3 < κ ≤ 4 | Method 2 (polynomial Γκ) | Correction < 5–10%; simple |
| 4 < κ ≤ 6 | Method 1 (Winterstein) or Method 2 | Check h4 monotonicity condition |
| 6 < κ ≤ 10 | Method 3 (distribution substitution) | Winterstein may be non-monotone |
| κ > 10 | Time-domain rainflow + Miner | Spectral methods unreliable; use simulation |
6. Worked Numerical Example
6.1 Problem Definition
A flat-topped stress PSD representative of an aerospace structural attachment is defined as:
G(ƒ) = 0.5 MPa^2/Hz, 20 Hz ≤ ƒ ≤ 80 Hz G(ƒ) = 0, otherwise(10)S-N parameters: b = 4, C = 2×1021 MPa4 cycle. The measured stress time history is stationary with kurtosis κ = 5.2 and negligible skewness.
6.2 Spectral Moments
m_0 = ∫_{20}^{80} 0.5 df = 0.5 × 60 = 30.0 MPa^2 ⇒ σ = 5.477 MPa m_1 = 0.5 × (80^2 – 20^2)/2 = 0.5 × 3000 = 1500.0 MPa^2·Hz m_2 = 0.5 × (80^3 – 20^3)/3 = 0.5 × 169333 = 84666.7 MPa^2·Hz^2 m_4 = 0.5 × (80^5 – 20^5)/5 = 0.5 × 655.4×10^6 = 327.7×10^6 MPa^2·Hz^4(11) x_m = (m_1/m_0) √(m_2/m_4) = 50.0 × 0.01607 = 0.804 α_2 = m_2 / √(m_0 m_4) = 84666.7 / √(30.0 × 327.7×10^6) = 0.854 ν_P = √(m_4/m_2) = √(3869) = 62.2 Hz(12)6.3 Dirlik Coefficients
D_1 = 2(0.804 – 0.854^2)/(1 + 0.854^2) = 2(0.074)/(1.729) = 0.086 R = (0.854 – 0.804 – 0.086^2)/(1 – 0.854 – 0.086 + 0.086^2) = 0.653 D_2 = (1 – 0.854 – 0.086 + 0.086^2)/(1 – 0.653) = 0.181 D_3 = 1 – 0.086 – 0.181 = 0.733 Q = 1.25(0.854 – 0.733 – 0.181×0.653)/0.086 = 0.700(13)6.4 Gaussian (Standard) Dirlik Damage Rate
s_0 = 2√m_0 = 2×5.477 = 10.954 MPa T_1 = D_1 Q^b Γ(b+1) = 0.086 × 0.700^4 × Γ(5) = 0.086 × 0.2401 × 24 = 0.496 T_2 = D_2 R^b Γ((b+2)/2) = 0.181 × 0.653^4 × Γ(3) = 0.181 × 0.1817 × 2.0 = 0.066 T_3 = D_3 Γ((b+2)/2) = 0.733 × 2.0 = 1.466 Ḋ_DK = (ν_P/C) × s_0^b × (T_1 + T_2 + T_3) = (62.2 / 2×10^21) × 10.954^4 × (0.496 + 0.066 + 1.466) = 3.11×10^-20 × 14404 × 2.028 = 9.09×10^-16 [damage/s](14)6.5 Non-Gaussian Correction (κ = 5.2, Method 2)
From Table 1 for b = 4: c1 = 0.062, c2 = 0.005. Excess kurtosis = 5.2 − 3 = 2.2:
Γ_κ = 1 + 0.062(2.2) + 0.005(2.2)^2 = 1 + 0.1364 + 0.0242 = 1.161(15) Ḋ_corrected = 1.161 × 9.09×10^-16 = 1.055×10^-15 [damage/s](16)6.6 Life Estimates
| Method | Damage rate (s−1) | Fatigue life |
|---|---|---|
| Gaussian Dirlik (uncorrected) | 9.09 × 10−16 | 34.8 years |
| Kurtosis-corrected (κ = 5.2) | 1.055 × 10−15 | 30.0 years |
| Life reduction | — | 13.8% |
7. Estimating Kurtosis from Measured Data
Reliable kurtosis estimation requires sufficient record length. The standard sample kurtosis estimator is:
κ_hat = [n(n+1)/((n-1)(n-2)(n-3))] × ∑[(x_i – x_bar)^4 / s^4] – 3(n-1)^2 / ((n-2)(n-3))(17)where n is the sample count and s is the sample standard deviation. As a rule of thumb, at least 5000 cycles (N ≥ 5000 peak-trough pairs) are needed to estimate kurtosis to within ±0.5 at the 95% confidence level for moderately leptokurtic processes.
Before applying any spectral method, verify stationarity by:
- Dividing the record into non-overlapping segments and testing for variance homogeneity (Levene test or Brown-Forsythe test).
- Checking that the running mean does not trend.
- Confirming the PSD shape is consistent across segments.
8. Limitations and Cautions
- Skewness. The correction methods above assume negligible skewness (γ1 ≈ 0). Compressive mean stress combined with asymmetric loading produces skewness that requires separate treatment.
- Narrowband processes. The Dirlik method itself is less accurate for very narrow bandwidth (α2 > 0.95). Apply the narrow-band correction (Rayleigh PDF) in that regime and superimpose the kurtosis factor.
- Plasticity. At very high kurtosis (κ > 6) peak amplitudes may approach yield. Linear damage accumulation (Miner’s rule) combined with spectral methods is then non-conservative; local plasticity redistributes damage.
- Non-stationarity. Intermittent high-amplitude bursts that cause κ > 3 may indicate non-stationarity (e.g., shock-on-random). Verify stationarity before applying these corrections; if the process is non-stationary, separate the burst and stationary components.
- Mean stress. These corrections do not account for mean stress effects (Goodman, Gerber). Apply mean stress correction to the S-N curve before computing damage.
9. Summary
The Dirlik spectral fatigue method is a well-validated tool for Gaussian stress processes. When measured stress records are stationary but leptokurtic (κ > 3), the standard method underestimates fatigue damage because it misses the excess probability mass in the high-amplitude tails. Three correction strategies are available: the Winterstein Hermite transformation (best for mild to moderate κ), the polynomial correction factor Γκ (simplest and robust to κ ≤ 8), and non-Gaussian amplitude distribution substitution (most general). For the worked aerospace example with κ = 5.2 and b = 4, ignoring kurtosis overestimates fatigue life by about 16%. For steeper S-N slopes the error is proportionally larger.
References
[1] Dirlik, T. (1985). Application of Computers in Fatigue Analysis. PhD Thesis, University of Warwick.
[2] Winterstein, S.R. (1988). Nonlinear vibration models for extremes and fatigue. Journal of Engineering Mechanics, 114(10), 1772–1790.
[3] Rizzi, S.A., Przekop, A., and Turner, T.L. (2011). On the response of a nonlinear structure to high kurtosis non-Gaussian random loadings. Eurodyn 2011 Proceedings.
[4] Benasciutti, D. and Tovo, R. (2006). Fatigue life assessment in non-Gaussian random loadings. International Journal of Fatigue, 28(7), 733–746.
[5] Johannesson, P. and Rychlik, I. (2013). Laplace processes for fatigue. Probabilistic Engineering Mechanics, 33, 1–10.
[6] Tovo, R. and Benasciutti, D. (2005). Spectral methods for lifetime prediction under wide-band stationary random processes. International Journal of Fatigue, 27(8), 867–877.
[7] Irvine, T. (2024). Fatigue Damage Spectrum. VibrationData ebook, vibrationdata.com.
[8] Newland, D.E. (1993). An Introduction to Random Vibrations, Spectral and Wavelet Analysis, 3rd ed. Longman.