频域疲劳计算（Dirlik 法）完整技术文档

频域疲劳计算（Dirlik 法）完整技术文档

目录

  1. 引言
  2. 频域疲劳计算总体流程
  3. 随机过程与谱矩基础
  4. 穿越理论详细推导（Rice 公式）
  5. Dirlik 法理论基础
  6. S-N 曲线与 Miner 损伤
  7. Python 代码实现
  8. 代码逐行数学解释
  9. 数值保护与稳定性
  10. 常见问题
  11. 关键注意事项
  12. 总结

1. 引言

频域疲劳分析是随机振动载荷下结构疲劳寿命评估的主流方法。相比时域雨流计数，频域法计算量小、效率高，特别适合长寿命随机振动问题。

Dirlik 法是目前工程上应用最广泛的频域疲劳方法，它用一个由指数分布与瑞利分布组合而成的经验公式，逼近雨流计数得到的应力幅值分布，对宽带和窄带随机过程均有良好精度。

2. 频域疲劳计算总体流程

输入载荷 PSD ──► 频响传递函数 H(f) ──► 应力 PSD S_σ(f)
                                          │
                                          ▼
                                    谱矩 m_n = ∫ f^n S_σ(f) df
                                          │
                                          ▼
                                    带宽参数 γ, X_m, ν_p
                                          │
                                          ▼
                                    Dirlik 参数 D1,D2,D3,R,Q
                                          │
                                          ▼
                                    幅值 PDF p(S)
                                          │
                                          ▼
                          Miner 损伤积分 D = ν_p T ∫ p(S)/N(S) dS
                                          │
                                          ▼
                                    疲劳寿命 T_f = T / D

3. 随机过程与谱矩基础

3.1 应力响应 PSD

由有限元频响分析得到传递函数 $$ H(f) $$，结合输入载荷 PSD $$ S_x(f) $$：

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

3.2 自相关函数与谱矩

设 $$ x(t) $$ 是零均值平稳随机过程，自相关函数：

\[R_x(\tau) = E[x(t) x(t+\tau)] \]

由 Wiener-Khinchin 定理：

\[S_x(f) = \int_{-\infty}^{\infty} R_x(\tau) e^{-j2\pi f\tau} d\tau \]

谱矩定义为：

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

自相关函数在 $$ \tau = 0 $$ 处的导数与谱矩的关系：

\[R_x(0) = m_0 \]

\[R_x'(0) = 0 \]

\[R_x''(0) = -(2\pi)^2 m_2 \]

\[R_x''''(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 各阶谱矩的物理意义

谱矩 公式 物理意义
$$ m_0 $$ $$ \int S_\sigma(f) df $$ 应力方差 $$ \sigma^2 $$
$$ m_1 $$ $$ \int f S_\sigma(f) df $$ 频率重心
$$ m_2 $$ $$ \int f^2 S_\sigma(f) df $$ 与过零率相关
$$ m_3 $$ $$ \int f^3 S_\sigma(f) df $$ 通常不参与 Dirlik
$$ m_4 $$ $$ \int f^4 S_\sigma(f) df $$ 与峰值率相关

4. 穿越理论详细推导（Rice 公式）

4.1 什么是穿越

在随机振动信号 $$ x(t) $$ 中：

  • 过零穿越：信号从负到正（或从正到负）穿过零点
  • 峰值穿越：信号出现局部最大值

单位时间内：

  • 过零次数 → 过零率 $$ \nu_0 $$
  • 峰值次数 → 峰值率 $$ \nu_p $$

4.2 过零率推导

4.2.1 穿越条件

在 $$ [t, t+dt] $$ 内发生正穿越：

\[x(t) < 0 \quad \text{且} \quad x(t+dt) > 0 \]

泰勒展开 $$ x(t+dt) \approx x(t) + \dot{x}(t) dt $$，条件变为：

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

4.2.2 期望穿越次数

单位时间正穿越次数：

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

4.2.3 高斯过程联合概率密度

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

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

\[p(x=0, \dot{x}) = \frac{1}{\sqrt{2\pi m_0}} \cdot \frac{1}{\sqrt{2\pi} \sigma_{\dot{x}}} e^{-\dot{x}^2/(2\sigma_{\dot{x}}^2)} \]

4.2.4 积分求解

\[\nu_0^+ = \int_0^\infty \dot{x} \cdot \frac{1}{\sqrt{2\pi m_0}} \cdot \frac{1}{\sqrt{2\pi} \sigma_{\dot{x}}} e^{-\dot{x}^2/(2\sigma_{\dot{x}}^2)} d\dot{x} \]

令 $$ u = \dot{x}2/(2\sigma_{\dot{x}}2) $$：

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

最终过零率（Hz）：

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

4.3 峰值率推导

峰值出现时 $$ \dot{x} = 0 $$ 且 $$ \ddot{x} < 0 $$。峰值率即速度过程 $$ \dot{x}(t) $$ 的过零率。

对 $$ \dot{x}(t) $$：

  • 方差：$$ E[\dot{x}^2] = (2\pi)^2 m_2 $$
  • 导数 $$ \ddot{x} $$ 的方差：$$ E[\ddot{x}^2] = (2\pi)^4 m_4 $$

套用过零率公式：

\[\nu_p = \sqrt{\frac{E[\ddot{x}^2]}{E[\dot{x}^2]}} = \sqrt{\frac{(2\pi)^4 m_4}{(2\pi)^2 m_2}} = 2\pi\sqrt{\frac{m_4}{m_2}} \]

用 Hz 表示：

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

4.4 不规则因子

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

$$ \gamma $$ 谱形状 信号特征
$$ \to 1 $$ 窄带 像正弦波，一峰一循环
$$ \to 0 $$ 宽带 像噪声，一循环多峰

4.5 数值验证

正弦波 $$ x(t) = A\sin(2\pi f_0 t) $$：

\[\nu_0 = f_0, \quad \nu_p = f_0, \quad \gamma = 1 \quad ✓ \]

白噪声 $$ S_x(f) = S_0, f \in [0,B] $$：

\[\nu_0 = \frac{B}{\sqrt{3}}, \quad \nu_p = B\sqrt{\frac{3}{5}}, \quad \gamma = \frac{\sqrt{5}}{3} \approx 0.745 \]

5. Dirlik 法理论基础

5.1 带宽参数

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

5.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 $$。

5.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} $$ 标准瑞利分布 大幅值循环

6. S-N 曲线与 Miner 损伤

6.1 S-N 曲线

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

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

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

6.2 Miner 线性累积损伤

每个幅值 $$ S $$ 作用 $$ n(S) $$ 次，消耗损伤 $$ n(S)/N(S) $$。总损伤：

\[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 \]

6.3 疲劳寿命

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

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

6.4 平均应力修正（Goodman）

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

7. Python 代码实现

7.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 法则）"""
    # 1. 谱矩
    moments = spectral_moments(f, S_sigma, 4)
    m0, m2, m4 = moments[0], moments[2], moments[4]

    # 2. 峰值率 nu_p = sqrt(m4/m2)，单位 Hz
    nu_p = np.sqrt(m4 / m2)

    # 3. Dirlik 参数
    D1, D2, D3, R, Q, gamma, Xm = dirlik_parameters(moments)

    # 4. 损伤积分
    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)

    # Goodman 平均应力修正
    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 * ∫ p(S)/N(S) dS
    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)

7.2 调用示例

if __name__ == "__main__":
    f = np.linspace(0, 200, 2000)
    S_sigma = 1e-4 * np.exp(-((f - 50)**2) / (2 * 10**2)) + 1e-7
    C, k = 1e12, 3.0
    T, sigma_m, sigma_u = 3600.0, 20.0, 500.0

    T_f, D, moments, params = fatigue_life_dirlik(
        f, S_sigma, T, C, k, sigma_m, sigma_u
    )

    print("谱矩 m0~m4:")
    for i, val in enumerate(moments):
        print(f"  m{i} = {val:.6e}")

    D1, D2, D3, R, Q, gamma, Xm = params
    print(f"\nDirlik 参数:")
    print(f"  gamma = {gamma:.6f}, Xm = {Xm:.6f}")
    print(f"  D1 = {D1:.6f}, D2 = {D2:.6f}, D3 = {D3:.6f}")
    print(f"  R = {R:.6f}, Q = {Q:.6f}")

    print(f"\n总损伤 D = {D:.6e}")
    print(f"疲劳寿命 T_f = {T_f:.2f} s = {T_f/3600:.2f} h")

7.3 最简版本

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

8. 代码逐行数学解释

8.1 谱矩计算

moments = spectral_moments(f, S_sigma, 4)

\[m_n = \int_0^\infty f^n S_\sigma(f) df \]

8.2 峰值率

nu_p = np.sqrt(m4 / m2)

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

8.3 Dirlik 参数

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

对应 5.1、5.2 节公式。

8.4 损伤积分

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)

\[S_{\max} = S_{\text{max\_factor}} \cdot \sqrt{m_0} \]

\[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] \]

8.5 Goodman 修正

S_eq_vals = S_vals / (1.0 - sigma_m / sigma_u)

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

8.6 损伤积分与寿命

integral = trapz(pdf_vals * S_eq_vals**k / C, S_vals)
D = nu_p * T * integral
T_f = T / D

\[D = \frac{\nu_p T}{C} \int p(S) S_{\text{eq}}^k dS, \quad T_f = \frac{T}{D} \]

8.7 对照表

代码 数学公式
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 $$

9. 数值保护与稳定性

判断 保护对象 触发条件 处理
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 \]

10. 常见问题

10.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

10.2 单位一致性

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

10.3 S-N 曲线形式转换

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

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

10.4 积分上限

默认 $$ S_{\max} = 15\sqrt{m_0} $$，覆盖 PDF 主要尾部；若尾部较长可适当增大 S_max_factor。

10.5 多轴疲劳

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

11. 关键注意事项

  1. 方法选择：Dirlik 法是宽带随机振动的主流方法；窄带过程可用更简单的窄带法或 Steinberg 三区间法。
  2. 平均应力修正：S-N 曲线通常基于对称循环（R=-1），非零平均应力需用 Goodman 公式修正。
  3. 积分上限：默认 $$ S_{\max} = 15\sqrt{m_0} $$，覆盖 PDF 主要尾部；若尾部较长可适当增大。
  4. 多轴疲劳：需先合成等效应力 PSD（如 Von Mises），再调用本函数。
  5. 计算效率：频域法仅处理 PSD 数据，比时域瞬态分析快几个数量级。
  6. 数值稳定性：若出现异常，检查 PSD 合理性或增大 n_points。

12. 总结

12.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

12.2 核心公式汇总

量 公式
谱矩 $$ m_n = \int 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} $$
Dirlik PDF $$ p(S) = \frac{1}{2\sqrt{m_0}}[\cdots] $$
Miner 损伤 $$ D = \nu_p T \int p(S)/N(S) dS $$
疲劳寿命 $$ T_f = T/D $$

12.3 一句话总结

谱矩 $$ m_2, m_4 $$ 通过 Rice 公式给出峰值率 $$ \nu_p $$，峰值率决定单位时间内的疲劳循环次数，再结合 Dirlik 分布 $$ p(S) $$ 和 S-N 曲线 $$ N(S) $$，就能用 Miner 法则算出总损伤和疲劳寿命 $$ T_f = T/D $$。

文档结束

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