频域疲劳计算（Dirlik 法）——完整技术总结

频域疲劳计算（Dirlik 法）——完整技术总结

一、总体概述

频域疲劳分析是随机振动载荷下结构疲劳寿命评估的主流方法。其核心思想是：将随机载荷的功率谱密度（PSD）转化为应力 PSD，提取谱矩，利用统计分布模型估计应力幅值的概率密度，最后结合 S-N 曲线和 Miner 线性累积损伤法则计算疲劳寿命。

相比时域雨流计数，频域法计算量小、效率高，特别适合长寿命随机振动问题。Dirlik 法是目前工程上应用最广泛的频域疲劳方法。

二、完整逻辑链（7 步）

\[\underbrace{S_x(f) \xrightarrow{\times |H(f)|^2} S_\sigma(f)}_{\text{第 1 步}} \]

\[\underbrace{S_\sigma(f) \xrightarrow{\int f^n df} m_n \quad (n = 0, 1, 2, 4)}_{\text{第 2 步}} \]

\[\underbrace{m_n \xrightarrow{\text{Rice 公式}} \nu_p = \sqrt{\frac{m_4}{m_2}}}_{\text{第 3 步}} \]

\[\underbrace{m_n \xrightarrow{\text{Dirlik 公式}} D_1, D_2, D_3, R, Q}_{\text{第 4 步}} \]

\[\underbrace{D_1, D_2, D_3, R, Q \xrightarrow{\text{Dirlik PDF}} p(S)}_{\text{第 5 步}} \]

\[\underbrace{p(S), \nu_p, T, C, k \xrightarrow{\text{Miner 法则}} D = \frac{\nu_p T}{C} \int_0^\infty p(S) S^k \, dS}_{\text{第 6 步}} \]

\[\underbrace{D \xrightarrow{T_f = T/D} T_f}_{\text{第 7 步}} \]

步骤汇总表

步 输入 变换 输出
1 $$ S_x(f) $$ $$ \times |H(f)|^2 $$ $$ S_\sigma(f) $$
2 $$ S_\sigma(f) $$ $$ \int f^n df $$ $$ m_n $$
3 $$ m_n $$ Rice 公式 $$ \nu_p $$
4 $$ m_n $$ Dirlik 公式 $$ D_1, D_2, D_3, R, Q $$
5 Dirlik 参数 Dirlik PDF $$ p(S) $$
6 $$ p(S), \nu_p, T, C, k $$ Miner 法则 $$ D $$
7 $$ D $$ $$ T_f = T/D $$ $$ T_f $$

三、理论基础

3.1 应力响应 PSD

\[S_\sigma(f) = |H(f)|^2 \, S_x(f) \]

3.2 谱矩

\[m_n = \int_0^\infty f^n \, S_\sigma(f) \, df, \quad n = 0,1,2,3,4 \]

谱矩 物理意义
$$ m_0 $$ 应力方差 $$ \sigma^2 $$
$$ m_1 $$ 频率重心
$$ m_2 $$ 与过零率相关
$$ m_4 $$ 与峰值率相关

自相关函数与谱矩的关系：

\[R(0) = m_0, \quad R''(0) = -(2\pi)^2 m_2, \quad R''''(0) = (2\pi)^4 m_4 \]

因此：

\[E[x^2] = m_0, \quad E[\dot{x}^2] = (2\pi)^2 m_2, \quad E[\ddot{x}^2] = (2\pi)^4 m_4 \]

3.3 穿越理论（Rice 公式）

3.3.1 泰勒展开

在 $$ [t, t+dt] $$ 内对信号作一阶泰勒展开：

\[x(t+dt) \approx x(t) + \dot{x}(t) \, dt \]

正穿越条件：

\[x(t) < 0, \quad x(t+dt) > 0 \]

代入展开式：

\[-\dot{x}(t) \, dt < x(t) < 0, \quad \dot{x}(t) > 0 \]

3.3.2 期望穿越次数

\[\nu_0^+ = \int_0^\infty \dot{x} \, p(x=0, \dot{x}) \, d\dot{x} \]

对高斯过程，$$ x $$ 与 $$ \dot{x} $$ 独立：

\[\sigma_x^2 = m_0, \quad \sigma_{\dot{x}}^2 = (2\pi)^2 m_2 \]

积分得：

\[\nu_0^+ = \frac{\sigma_{\dot{x}}}{2\pi \sqrt{m_0}} = \sqrt{\frac{m_2}{m_0}} \]

3.3.3 过零率与峰值率

\[\boxed{\nu_0 = \sqrt{\frac{m_2}{m_0}} \quad (\text{Hz})} \]

\[\boxed{\nu_p = \sqrt{\frac{m_4}{m_2}} \quad (\text{Hz})} \]

3.3.4 不规则因子

\[\boxed{\gamma = \frac{\nu_0}{\nu_p} = \frac{m_2}{\sqrt{m_0 m_4}}} \]

$ \gamma \to 1 $ 为窄带，$ \gamma \to 0 $ 为宽带。

3.4 高斯过程正向穿越率

\[v_a^+ = \frac{1}{2\pi} \sqrt{\frac{m_2}{m_0}} e^{-a^2 / 2m_0} \]

当 $$ a = 0 $$ 时：

\[v_0^+ = \frac{1}{2\pi} \sqrt{\frac{m_2}{m_0}} \]

3.5 峰值概率密度与瑞利分布

峰值超过高度 $$ a $$ 的概率：

\[\int_a^\infty p_p(a) \, da = \frac{v_a^+}{v_0^+} \]

对 $$ a $$ 求导：

\[p_p(a) = -\frac{1}{v_0^+} \frac{d v_a^+}{da} \]

代入高斯穿越率，得瑞利分布：

\[\boxed{p_p(a) = \frac{a}{m_0} e^{-a^2 / 2m_0}, \quad a \ge 0} \]

峰值期望值：

\[E[a] = \sqrt{\frac{\pi m_0}{2}} \]

四、Dirlik 法理论基础

4.1 带宽参数

\[\gamma = \frac{m_2}{\sqrt{m_0 m_4}}, \qquad X_m = \frac{m_1}{m_0} \sqrt{\frac{m_2}{m_4}} \]

4.2 Dirlik 参数

\[D_1 = \frac{2(X_m - \gamma^2)}{1 + \gamma^2} \]

\[R = \frac{\gamma - X_m - D_1^2}{1 - \gamma - D_1 + D_1^2} \]

\[D_2 = \frac{1 - \gamma - D_1 + D_1^2}{1 - R} \]

\[D_3 = 1 - D_1 - D_2 \]

\[Q = \frac{1.25(\gamma - D_3 - D_2 R)}{D_1} \]

约束：$$ D_1 + D_2 + D_3 = 1 $$。

4.3 Dirlik 概率密度函数

\[p(S) = \frac{1}{2\sqrt{m_0}} \left[ \frac{D_1}{Q} e^{-Z/Q} + \frac{D_2 Z}{R^2} e^{-Z^2/(2R^2)} + D_3 Z e^{-Z^2/2} \right], \quad Z = \frac{S}{2\sqrt{m_0}} \]

项 分布类型 描述
$$ \frac{D_1}{Q} e^{-Z/Q} $$ 指数分布 小幅值循环
$$ \frac{D_2 Z}{R^2} e{-Z2/(2R^2)} $$ 瑞利分布 中幅值循环
$$ D_3 Z e{-Z2/2} $$ 标准瑞利分布 大幅值循环

五、S-N 曲线与 Miner 损伤

5.1 S-N 曲线

\[N(S) = C \, S^{-k} \]

若曲线形式为 $$ N = (S/S_f)^{1/b} $$，则转换：

\[C = S_f^{-1/b}, \quad k = -1/b \]

5.2 Miner 线性累积损伤

离散形式：

\[D = \sum_i \frac{n_i}{N_i} \]

连续形式：

\[D = \int_0^\infty \frac{n(S)}{N(S)} dS \]

在时间 $$ T $$ 内，幅值落在 $$ [S, S+dS] $$ 的循环次数：

\[n(S) dS = \nu_p \cdot T \cdot p(S) dS \]

代入：

\[D = \nu_p T \int_0^\infty \frac{p(S)}{N(S)} dS = \frac{\nu_p T}{C} \int_0^\infty p(S) S^k \, dS \]

5.3 疲劳寿命

当 $$ D = 1 $$ 时失效：

\[\boxed{T_f = \frac{T}{D}} \]

5.4 Goodman 平均应力修正

\[S_{\text{eq}} = \frac{S}{1 - \sigma_m / \sigma_u} \]

六、Python 代码实现

6.1 完整版本

import numpy as np
from scipy.integrate import trapezoid as trapz


def spectral_moments(f, S_sigma, n_max=4):
    """计算谱矩 m_n = ∫ f^n * S_sigma(f) df"""
    moments = []
    for n in range(n_max + 1):
        integrand = f**n * S_sigma
        moments.append(trapz(integrand, f))
    return np.array(moments)


def dirlik_parameters(moments):
    """根据谱矩计算 Dirlik 概率密度函数的参数"""
    m0, m1, m2, m4 = moments[0], moments[1], moments[2], moments[4]

    gamma = m2 / np.sqrt(m0 * m4)
    Xm = (m1 / m0) * np.sqrt(m2 / m4)

    D1 = 2.0 * (Xm - gamma**2) / (1.0 + gamma**2)

    denom_R = 1.0 - gamma - D1 + D1**2
    R = 0.0 if abs(denom_R) < 1e-12 else (gamma - Xm - D1**2) / denom_R
    D2 = 0.0 if abs(1.0 - R) < 1e-12 else (1.0 - gamma - D1 + D1**2) / (1.0 - R)
    D3 = 1.0 - D1 - D2

    if abs(D1) < 1e-12:
        Q = 0.0
    else:
        Q = 1.25 * (gamma - D3 - D2 * R) / D1
        if Q <= 0:
            Q = 1e-12
    return D1, D2, D3, R, Q, gamma, Xm


def dirlik_pdf(S, m0, D1, D2, D3, R, Q):
    """Dirlik 应力幅值概率密度函数 p(S)"""
    S = np.asarray(S, dtype=float)
    sqrt_m0 = np.sqrt(m0)
    Z = S / (2.0 * sqrt_m0)

    Q_safe = max(Q, 1e-12)
    R2_safe = max(R**2, 1e-12)

    term1 = (D1 / Q_safe) * np.exp(-Z / Q_safe)
    term2 = (D2 * Z / R2_safe) * np.exp(-Z**2 / (2.0 * R2_safe))
    term3 = D3 * Z * np.exp(-Z**2 / 2.0)

    return (1.0 / (2.0 * sqrt_m0)) * (term1 + term2 + term3)


def fatigue_life_dirlik(f, S_sigma, T, C, k, sigma_m=0.0, sigma_u=None,
                        S_max_factor=15.0, n_points=20000):
    """频域疲劳寿命计算（Dirlik 法 + Miner 法则）"""
    moments = spectral_moments(f, S_sigma, 4)
    m0, m2, m4 = moments[0], moments[2], moments[4]

    nu_p = np.sqrt(m4 / m2)

    D1, D2, D3, R, Q, gamma, Xm = dirlik_parameters(moments)

    S_vals = np.linspace(1e-6, S_max_factor * np.sqrt(m0), n_points)
    pdf_vals = dirlik_pdf(S_vals, m0, D1, D2, D3, R, Q)

    if sigma_u and sigma_u > 0:
        S_eq_vals = S_vals / (1.0 - sigma_m / sigma_u)
    else:
        S_eq_vals = S_vals

    integral = trapz(pdf_vals * S_eq_vals**k / C, S_vals)
    D = nu_p * T * integral
    T_f = T / D if D > 0 else np.inf
    return T_f, D, moments, (D1, D2, D3, R, Q, gamma, Xm)

6.2 最简版本

import numpy as np
from scipy.integrate import trapezoid as trapz

def freq_fatigue_dirlik(f, S_sigma, T, C, k, sigma_m=0, sigma_u=None):
    m = [trapz(f**n * S_sigma, f) for n in range(5)]
    m0, m1, m2, m4 = m[0], m[1], m[2], m[4]
    nu_p = np.sqrt(m4 / m2)
    gamma = m2 / np.sqrt(m0 * m4)
    Xm = (m1 / m0) * np.sqrt(m2 / m4)
    D1 = 2 * (Xm - gamma**2) / (1 + gamma**2)
    R = (gamma - Xm - D1**2) / (1 - gamma - D1 + D1**2)
    D2 = (1 - gamma - D1 + D1**2) / (1 - R)
    D3 = 1 - D1 - D2
    Q = 1.25 * (gamma - D3 - D2 * R) / D1
    S = np.linspace(1e-6, 15 * np.sqrt(m0), 20000)
    Z = S / (2 * np.sqrt(m0))
    p = (D1/Q * np.exp(-Z/Q) + D2*Z/R**2 * np.exp(-Z**2/(2*R**2)) + D3*Z * np.exp(-Z**2/2)) / (2 * np.sqrt(m0))
    if sigma_u:
        S = S / (1 - sigma_m / sigma_u)
    D = nu_p * T * trapz(p * S**k / C, S)
    return T / D, D

6.3 代码与公式对照

代码 数学公式
spectral_moments(f, S_sigma, 4) $$ m_n = \int f^n S_\sigma(f) df $$
nu_p = np.sqrt(m4 / m2) $$ \nu_p = \sqrt{m_4/m_2} $$
dirlik_parameters(moments) $$ \gamma, X_m, D_1, D_2, D_3, R, Q $$
S_vals = np.linspace(...) $$ S \in [0, 15\sqrt{m_0}] $$
dirlik_pdf(...) $$ p(S) $$
S_eq_vals = S_vals / (1 - sigma_m/sigma_u) $$ S_{\text{eq}} = \dfrac{S}{1-\sigma_m/\sigma_u} $$
trapz(pdf_vals * S_eq_vals**k / C, S_vals) $$ \int \dfrac{p(S) S_{\text{eq}}^k}{C} dS $$
D = nu_p * T * integral $$ D = \nu_p T \int \frac{p(S)}{N(S)} dS $$
T_f = T / D $$ T_f = T/D $$

七、Miner 法则结合详解

7.1 核心思想

Miner 法则：每个幅值消耗一点寿命，总损伤线性叠加：

\[D = \sum_i \frac{n_i}{N_i} \quad \longrightarrow \quad D = \int_0^\infty \frac{n(S)}{N(S)} dS \]

7.2 频域里各量的来源

量 含义 来源
$$ \nu_p T $$ 总循环次数 峰值率 × 时间
$$ p(S) $$ 幅值概率分布 Dirlik 公式
$$ n(S) = \nu_p T p(S) dS $$ 幅值 $$ S $$ 的循环次数 总循环数 × 概率
$$ N(S) = C S^{-k} $$ 幅值 $$ S $$ 的失效循环数 S-N 曲线

7.3 结合过程

\[D = \int_0^\infty \frac{n(S)}{N(S)} dS = \nu_p T \int_0^\infty \frac{p(S)}{N(S)} dS = \frac{\nu_p T}{C} \int_0^\infty p(S) S^k dS \]

\[T_f = \frac{T}{D} \]

7.4 直观理解（"寿命预算"记账）

概念 类比
S-N 曲线 $$ N(S) $$ 每个幅值的"寿命预算"
循环数 $$ n(S) $$ 每个幅值的"消费次数"
损伤 $$ n/N $$ 预算消耗比例
总损伤 $$ D $$ 预算消耗总和
$$ D = 1 $$ 预算花完，失效
$$ T_f = T/D $$ 预算能撑多久

八、数值保护

判断 保护对象 触发条件 处理
abs(denom_R) < 1e-12 R 的分母 $$ 1-\gamma-D_1+D_1^2 \approx 0 $$ $$ R=0 $$
abs(1.0 - R) < 1e-12 D2 的分母 $$ R \approx 1 $$ $$ D_2=0 $$
abs(D1) < 1e-12 Q 的分母 $$ D_1 \approx 0 $$ $$ Q=0 $$
Q <= 0 Q 本身 计算值为负或零 $$ Q=1e-12 $$

物理背景：

  • \[ \gamma \to 1 $$（窄带）：$$ D_1 \to 0 $$，Q 分母趋近 0 \]

九、常见问题

9.1 ImportError: cannot import name 'trapz'

SciPy 1.14.0 起移除了 trapz，更名为 trapezoid：

from scipy.integrate import trapezoid as trapz

或降级：pip install scipy==1.12.0

9.2 单位一致性

量 单位
频率 Hz
PSD 应力²/Hz
S-N 曲线 与 PSD 应力单位一致
时间 s

9.3 S-N 曲线形式转换

若为 $$ N = (S/S_f)^{1/b} $$：

\[C = S_f^{-1/b}, \quad k = -1/b \]

9.4 积分上限

默认 $$ S_{\max} = 15\sqrt{m_0} $$，覆盖 PDF 主要尾部。

9.5 多轴疲劳

需先合成等效应力 PSD（如 Von Mises），再调用本函数。

十、关键注意事项

  1. 方法选择：Dirlik 法适合宽带随机振动；窄带过程可用窄带法或 Steinberg 三区间法。
  2. 平均应力修正：非零平均应力需用 Goodman 公式修正。
  3. 积分上限：默认 $$ S_{\max} = 15\sqrt{m_0} $$，若尾部较长可适当增大。
  4. 多轴疲劳：需先合成等效应力 PSD。
  5. 计算效率：频域法比时域瞬态分析快几个数量级。
  6. 数值稳定性：异常时检查 PSD 合理性或增大 n_points。

十一、完整流程总结

11.1 完整计算链条

输入载荷 PSD
   │
   ▼
应力 PSD S_σ(f) = |H(f)|² S_x(f)
   │
   ▼
谱矩 m_n = ∫ f^n S_σ(f) df
   │
   ├─► m0 ──► σ² ──► 幅值尺度 √m0
   │
   ├─► m2, m4 ──► νp = √(m4/m2) ──► 每秒循环次数
   │
   └─► Dirlik ──► p(S) ──► 各幅值概率分布
                                │
                                ▼
             n(S) = νp · T · p(S) dS
                                │
                                ▼
             D = νp T ∫ p(S)/N(S) dS
                                │
                                ▼
             T_f = T / D

11.2 核心公式汇总表

步骤 公式
应力 PSD $$ S_\sigma(f) = |H(f)|^2 S_x(f) $$
谱矩 $$ m_n = \int_0^\infty f^n S_\sigma(f) df $$
过零率 $$ \nu_0 = \sqrt{m_2/m_0} $$
峰值率 $$ \nu_p = \sqrt{m_4/m_2} $$
不规则因子 $$ \gamma = m_2/\sqrt{m_0 m_4} $$
高斯穿越率 $$ v_a^+ = \frac{1}{2\pi}\sqrt{\frac{m_2}{m_0}} e{-a2/2m_0} $$
峰值分布 $$ p_p(a) = \frac{a}{m_0} e{-a2/2m_0} $$
Dirlik 参数 $$ D_1, D_2, D_3, R, Q $$
Dirlik PDF $$ p(S) = \frac{1}{2\sqrt{m_0}}[\cdots] $$
S-N 曲线 $$ N(S) = C S^{-k} $$
Miner 损伤 $$ D = \frac{\nu_p T}{C} \int_0^\infty p(S) S^k dS $$
疲劳寿命 $$ T_f = T/D $$
Goodman 修正 $$ S_{\text{eq}} = \frac{S}{1-\sigma_m/\sigma_u} $$

十二、一句话总结

从载荷 PSD 出发，通过传递函数得到应力 PSD，积分得到谱矩 $$ m_0, m_1, m_2, m_4 $$；由 $$ m_2, m_4 $$ 通过 Rice 公式得到峰值率 $$ \nu_p = \sqrt{m_4/m_2} $$；由谱矩通过 Dirlik 公式得到概率密度 $$ p(S) $$；再用 Miner 法则把总循环数 $$ \nu_p T $$ 按 $$ p(S) $$ 分摊到各幅值，除以 S-N 曲线 $$ N(S) = C S^{-k} $$，积分得总损伤

\[D = \frac{\nu_p T}{C} \int_0^\infty p(S) S^k dS \]

最终疲劳寿命

\[T_f = \frac{T}{D} \]

posted @ 2026-10-10 21:17  redufa  阅读(25)  评论(0)    收藏  举报