频域疲劳（二）：频域疲劳的总体思路 10

2.1总体概述

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

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

---

2.2 完整逻辑链（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 步}} \]

3.3 步骤汇总表

步	输入	变换	输出
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.4 计算过程可视化

3.4.1 画图

3.4.2 代码

import numpy as np
from scipy.integrate import quad
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import warnings

# 屏蔽警告
warnings.filterwarnings("ignore")

# 设置中文字体（防止图表上的中文变方块）
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False 
plt.rcParams.update({'font.size': 11})

def run_optimized_visualization():
    # ==========================================
    # 0. 初始设定与计算
    # ==========================================
    f = np.linspace(0.1, 150, 2000) 
    Sx_f = 20 * np.exp(-((f - 20)**2) / 200) + 10 * np.exp(-((f - 80)**2) / 300)
    H_f_squared = 1.0 / (1.0 + (f/100.0)**4) * 2.0 
    T = 3600.0
    C = 2.0e12
    k = 3.0

    # Step 1
    So_f = Sx_f * H_f_squared
    So_func = interp1d(f, So_f, kind='linear', fill_value=0, bounds_error=False)

    # Step 2
    def integrand_mn(freq, n): return (freq**n) * So_func(freq)
    m0, _ = quad(integrand_mn, 0, 200, args=(0,), limit=500)
    m2, _ = quad(integrand_mn, 0, 200, args=(2,), limit=500)
    m4, _ = quad(integrand_mn, 0, 200, args=(4,), limit=500)
    m_list = [m0, m2, m4]

    # Step 3
    nu_p = np.sqrt(m4 / m2)

    # Step 4
    gamma = m2 / np.sqrt(m0 * m4)
    Xm = 0 
    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
    dirlik_params = [D1, D2, D3, R, Q]

    # Step 5
    def dirlik_pdf(S):
        if S <= 0: return 0.0
        Z = S / (2 * np.sqrt(m0))
        t1 = (D1 / Q) * np.exp(-Z / Q) if Q > 0 else 0
        t2 = (D2 * Z / R**2) * np.exp(-Z**2 / (2 * R**2)) if R > 0 else 0
        t3 = D3 * Z * np.exp(-Z**2 / 2) if D3 > 0 else 0
        pdf_val = (t1 + t2 + t3) / (2 * np.sqrt(m0))
        return max(0.0, pdf_val)

    # Step 6
    def damage_integrand(S): return dirlik_pdf(S) * (S**k)
    integral_val, _ = quad(damage_integrand, 0, 500, limit=500)
    D_total = (nu_p * T / C) * integral_val

    # Step 7
    Tf = T / D_total
    Tf_hours = Tf / 3600

    # ==========================================
    # 📊 优化后的图形排版 (4行2列)
    # ==========================================
    fig = plt.figure(figsize=(16, 14)) 
    fig.suptitle('频域疲劳寿命评估可视化报告 (频域 Dirlik 方法)', 
                 fontsize=20, fontweight='bold', y=0.96)
    
    gs = gridspec.GridSpec(4, 2, figure=fig, hspace=0.5, wspace=0.25, height_ratios=[1, 1, 1, 0.6])

    S_range = np.linspace(0.01, 120, 500)

    # --- 画图 1 (1行1列): PSD计算 ---
    ax1 = fig.add_subplot(gs[0, 0])
    # 修复：在字符串前面加上 r
    ax1.plot(f, Sx_f, label=r'输入 $S_x(f)$', color='#1f77b4', linestyle='--', alpha=0.8)
    ax1.plot(f, So_f, label=r'响应 $S_o(f)$', color='#d62728', linewidth=2.5)
    ax1.fill_between(f, So_f, color='#d62728', alpha=0.15)
    ax1.set_title('第1步: 系统频响 PSD 计算', fontsize=14, fontweight='bold', color='#333333')
    ax1.set_xlabel('频率 f (Hz)')
    ax1.set_ylabel('PSD (MPa²/Hz)')
    ax1.legend(loc='upper right')
    ax1.grid(True, linestyle=':', alpha=0.6)

    # --- 画图 2 (1行2列): 谱矩柱状图 ---
    ax2 = fig.add_subplot(gs[0, 1])
    bars = ax2.bar(['m0', 'm2', 'm4'], m_list, color=['#4C72B0', '#55A868', '#C44E52'], width=0.5)
    # 修复：在字符串前面加上 r
    ax2.set_title(r'第2步: 谱矩 $m_n$ 数值大小', fontsize=14, fontweight='bold', color='#333333')
    ax2.set_ylabel('数值 (对数坐标)')
    ax2.set_yscale('log') 
    for bar in bars:
        yval = bar.get_height()
        ax2.text(bar.get_x() + bar.get_width()/2.0, yval*1.2, f'{yval:.2e}', ha='center', va='bottom', fontsize=11, fontweight='bold')
    ax2.grid(True, axis='y', linestyle=':', alpha=0.6)

    # --- 画图 3 (2行1列): Rice 公式看板 ---
    ax3 = fig.add_subplot(gs[1, 0])
    ax3.axis('off')
    # 修复：在字符串前面加上 r (注意里面有 f-string 时，外部也要加 r)
    rice_text = (
        "【 第3步: Rice 公式计算 】\n\n"
        "期望峰值频率:\n"
        r"$\nu_p = \sqrt{m_4/m_2}$" "\n\n"
        f"$\\nu_p = $ {nu_p:.2f} Hz\n\n"
        "💡 物理意义:\n"
        "系统每秒发生的平均应力峰值次数"
    )
    ax3.text(0.5, 0.5, rice_text, fontsize=15, ha='center', va='center', linespacing=1.8,
             bbox=dict(boxstyle="round,pad=0.8", facecolor="#e8f4f8", edgecolor="#4C72B0", linewidth=2))

    # --- 画图 4 (2行2列): Dirlik 参数柱状图 ---
    ax4 = fig.add_subplot(gs[1, 1])
    colors = ['#2ca02c' if val > 0 else '#ff7f0e' for val in dirlik_params] 
    bars = ax4.bar(['D1', 'D2', 'D3', 'R', 'Q'], dirlik_params, color=colors, alpha=0.85)
    ax4.axhline(0, color='black', linewidth=1.5)
    ax4.set_title('第4步: Dirlik 形状参数 (橙色为负数修正项)', fontsize=14, fontweight='bold', color='#333333')
    ax4.set_ylabel('参数值')
    for i, v in enumerate(dirlik_params):
        ax4.text(i, v + (0.05 if v>0 else -0.05), f'{v:.2f}', ha='center', va='bottom' if v>0 else 'top', fontsize=11)
    ax4.grid(True, axis='y', linestyle=':', alpha=0.6)

    # --- 画图 5 (3行1列): 概率密度函数 PDF ---
    ax5 = fig.add_subplot(gs[2, 0])
    pdf_vals = [dirlik_pdf(S) for S in S_range]
    ax5.plot(S_range, pdf_vals, color='#9467bd', linewidth=3)
    ax5.fill_between(S_range, pdf_vals, color='#9467bd', alpha=0.25)
    # 修复：在字符串前面加上 r
    ax5.set_title(r'第5步: 应力幅值概率密度函数 $p(S)$', fontsize=14, fontweight='bold', color='#333333')
    ax5.set_xlabel('应力幅值 S (MPa)')
    ax5.set_ylabel(r'概率密度 $p(S)$')
    ax5.grid(True, linestyle=':', alpha=0.6)

    # --- 画图 6 (3行2列): Miner 积分面积 ---
    ax6 = fig.add_subplot(gs[2, 1])
    damage_vals = [damage_integrand(S) for S in S_range]
    ax6.plot(S_range, damage_vals, color='#ff7f0e', linewidth=3)
    ax6.fill_between(S_range, damage_vals, color='#ff7f0e', alpha=0.4, label='积分面积 = 总损伤基数')
    ax6.set_title('第6步: Miner 法则损伤积分域', fontsize=14, fontweight='bold', color='#333333')
    ax6.set_xlabel('应力幅值 S (MPa)')
    # 修复：在字符串前面加上 r
    ax6.set_ylabel(r'损伤贡献 $p(S) \cdot S^k$')
    ax6.legend(loc='upper right', fontsize=12)
    ax6.grid(True, linestyle=':', alpha=0.6)

    # --- 画图 7 (4行合并): 最终寿命结论 ---
    ax7 = fig.add_subplot(gs[3, :]) 
    ax7.axis('off')
    
    # 修复：在 f-string 前面加上 r
    life_text = (
        f"【 第7步: 疲劳寿命最终预测 】\n\n"
        rf"累计总损伤 $D = {D_total:.4e}$" "\n\n"
        rf"计算寿命 $T_f = T / D = {Tf:,.0f} \text{{ 秒}}$  ( 约 {Tf_hours:,.2f} 小时 )"
    )
    ax7.text(0.5, 0.5, life_text, fontsize=20, ha='center', va='center', fontweight='bold', color='#b8860b',
             bbox=dict(boxstyle="round,pad=1.0", facecolor="#fff8dc", edgecolor="#daa520", linewidth=3))

    plt.subplots_adjust(top=0.92, bottom=0.05, left=0.08, right=0.92)
    plt.show()

if __name__ == "__main__":
    run_optimized_visualization()