频域疲劳（一）：从时域到频域 10.21

1. 米勒法则 (Miner's Rule)

米勒法则的物理核心是“线性累积损伤”。假设一个零件承受了多种不同大小的交变应力(\(S_{a1}, S_{a2}, S_{a3}...\)) ，每种应力作用的次数分别是 \(n_1, n_2, n_3...\)，对应这种应力下材料能撑住的极限寿命是 \(N_1, N_2, N_3...\)。

那么，总疲劳损伤 \(D\) 就是每一次微小损伤的总和：

\[D = \sum_{i=1}^{k} \frac{n_i}{N_i} \]

(当 \(D \ge 1\) 时，零件发生疲劳破坏)

2.S-N 曲线

根据材料的疲劳特性（Basquin 方程），应力幅 \(S_a\) 和极限寿命 \(N\) 的关系为：

\[N = \frac{C}{S_a^b} \]

把它代入米勒公式中，单个应力循环造成的损伤就变成了 \(\frac{1}{N} = \frac{S_a^b}{C}\)。
那么，\(n_i\) 次循环造成的损伤就是：

\[D = \sum_{i=1}^{k} n_i \frac{S_{a,i}^b}{C} \]

3.“离散”到“连续”

应力信号是连续且随机的。我们很难数出“精确的几次”，可以通过概率密度函数 \(p(S_a)\) 来描述应力幅值的分布。

1. 总循环次数：假设应力循环的平均频率是 \(\nu\) (次/秒)，总时间是 \(T\) (秒)，那么总循环次数 \(N_{total} = \nu T\)。
2. 某个应力区间的次数：在总次数中，大小处于 \([S_a, S_a + dS_a]\) 区间的循环次数，就变成了 \(N_{total} \times p(S_a)dS_a\)。
3. 连续化求和（积分）：将上述区间的循环次数，乘以该应力下的单次损伤 \(\frac{S_a^b}{C}\)，然后因为应力是连续的，把求和号 \(\Sigma\) 变成积分号 \(\int\)：

\[D = \int_{0}^{\infty} (\nu T \cdot p(S_a)dS_a) \cdot \left( \frac{S_a^b}{C} \right) \]

把常数 \(\nu T\) 提取到最外面：

\[D = \nu T \int_{0}^{\infty} \frac{S_a^b}{C} p(S_a)dS_a \]

4.从时域到频域

这两个总公式如下（注意：两者的数学结构完全一样，区别仅在于 \(p(S_a)\)的获取方式和 \(T\)的物理意义）：

1. 时域总疲劳损伤：

\[D = \nu T \int_0^\infty \frac{S_a^b}{C} p(S_a) dS_a \]

(其中 \(\nu T = N_{total}\)为总循环次数，\(p(S_a)\)为雨流计数得到的概率密度)

2. 频域总疲劳损伤：

\[D = \nu T \int_0^\infty \frac{S_a^b}{C} p(S_a) dS_a \]

(其中 \(\nu = \sqrt{m_4/m_2}\)为峰值率，\(p(S_a)\)为基于谱矩假设的解析概率密度，如 Rayleigh 或 Dirlik)

5.米勒法则+SN 画图

5.1 画图

5.2 代码

import numpy as np
import matplotlib.pyplot as plt

# 设置中文字体（解决图表中文乱码）
plt.rcParams['font.sans-serif'] = ['SimHei']  
plt.rcParams['axes.unicode_minus'] = False

# ================= 1. 模拟数据与 S-N 曲线参数 =================
stress_levels = np.array([100, 200, 300, 400, 500]) # 离散的应力水平
N_i = np.array([100000, 30000, 8000, 2000, 500])    # 对应的极限寿命
n_i = np.array([20000, 15000, 6000, 1500, 400])     # 实际循环次数

# 假设 Basquin 方程参数：N = C / (S_a)^b
b = 2.5 
C = 100000 * (100**b) # 反推常数 C

# 生成连续的 S-N 曲线数据
S_curve = np.linspace(80, 550, 200)
N_curve = C / (S_curve ** b)

# ================= 2. 计算权重与总损伤 =================
damage_weights = n_i / N_i
total_damage_D = np.sum(damage_weights)

# ================= 3. 开始绘图 =================
fig, ax1 = plt.subplots(figsize=(11, 6))

# ---------- Y轴1：直方图 (左) ----------
width = 40
rects = ax1.bar(stress_levels, n_i, width, color='skyblue', edgecolor='black', 
                label=r'实际循环次数 $n_i$', zorder=3)
ax1.set_xlabel(r'应力幅值 $S_{a}$ (MPa)', fontsize=12)

# 将左Y轴设为对数坐标
ax1.set_yscale('log')
ax1.set_ylabel(r'循环次数 $n_i$ / 寿命 $N_i$ (对数坐标)', fontsize=12, color='tab:blue')
ax1.tick_params(axis='y', labelcolor='tab:blue')
ax1.set_xticks(stress_levels)
ax1.set_ylim(100, 1000000) 

# ---------- Y轴1 副图：S-N 曲线 (左，共享坐标) ----------
line_sn = ax1.plot(S_curve, N_curve, color='darkgreen', linestyle='--', linewidth=2, 
                   label=r'S-N曲线 (材料极限寿命 $N$)', zorder=1)

# 标记材料在对应离散应力下的极限寿命点
ax1.scatter(stress_levels, N_i, color='darkgreen', marker='^', s=60, zorder=4, 
            label=r'对应应力下的极限寿命 $N_i$')

# ---------- Y轴2：损伤权重 (右) ----------
ax2 = ax1.twinx()
line = ax2.plot(stress_levels, damage_weights, color='crimson', marker='o', 
                linestyle='-', linewidth=2, markersize=8, label=r'区间损伤权重 $n_i/N_i$', zorder=5)
ax2.set_ylabel(r'局部损伤权重 ($n_i / N_i$)', fontsize=12, color='crimson')
ax2.tick_params(axis='y', labelcolor='crimson')
ax2.set_ylim(0, max(damage_weights) * 1.3)

# 标注局部权重数值
for i in range(len(stress_levels)):
    ax2.annotate(f'{damage_weights[i]:.3f}', 
                 xy=(stress_levels[i], damage_weights[i]), 
                 xytext=(0, 10), textcoords='offset points',
                 ha='center', color='crimson', fontsize=10)

# ---------- 文本框与图例 ----------
# 【已修复】去掉易导致乱码的 Emoji 和 \mathbf
status_text = '已发生疲劳破坏 (破坏)' if total_damage_D >= 1.0 else '安全 (未破坏)'

# 将纯文本与 LaTeX 公式合理拼接
textstr = '\n'.join((
    r'总疲劳损伤 $D = \sum (n_i / N_i)$',
    f'总损伤 $D = {total_damage_D:.3f}$',
    f'状态: {status_text}'
))

props = dict(boxstyle='round', facecolor='wheat', alpha=0.9)
ax1.text(0.02, 0.95, textstr, transform=ax1.transAxes, fontsize=12,
        verticalalignment='top', bbox=props)

# 合并图例 (ax1 和 ax2 的图例合并)
lines_1, labels_1 = ax1.get_legend_handles_labels()
lines_2, labels_2 = ax2.get_legend_handles_labels()
ax1.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper right')

plt.title('米勒法则综合可视化：应力循环、S-N曲线与寿命损伤权重', fontsize=14, pad=15)
plt.grid(axis='y', linestyle='--', alpha=0.5, which='both', zorder=0)

plt.tight_layout()
plt.show()