第九十七篇：热疲劳寿命预测

摘要

热疲劳是工程结构在循环热载荷作用下产生的渐进性损伤现象，是导致航空发动机涡轮叶片、核反应堆压力容器、汽车排气系统等高温部件失效的主要机制之一。本文系统阐述热疲劳的物理机理、损伤累积理论和寿命预测方法，重点介绍基于应变-寿命法和能量法的预测模型，包括Coffin-Manson方程、Smith-Watson-Topper(SWT)参数和Ostergren能量模型。通过Python仿真实现热机械耦合分析、循环应力-应变响应计算、损伤累积追踪和剩余寿命预测，为高温部件的可靠性设计和安全评估提供理论基础和数值工具。

关键词

热疲劳，寿命预测，损伤累积，Coffin-Manson方程，应变-寿命法，热机械耦合，蠕变-疲劳交互作用，Miner法则

---

1. 热疲劳基础理论

1.1 热疲劳的定义与特征

热疲劳(Thermal Fatigue)是指材料或结构在循环变化的热载荷作用下，由于温度梯度引起的热应力反复作用而产生的疲劳损伤现象。与机械疲劳不同，热疲劳具有以下显著特征：

温度-应力耦合特性：热疲劳过程中温度场和应力场相互耦合，温度变化产生热应力，而材料力学性能又随温度变化。这种耦合效应使得热疲劳分析比等温疲劳更为复杂。

非等温循环特征：热疲劳通常涉及温度幅值和平均温度的同时变化，材料在不同温度水平下的力学响应存在显著差异，包括弹性模量、屈服强度、延伸率等参数的温度依赖性。

蠕变-疲劳交互作用：在高温条件下，蠕变变形与疲劳损伤相互促进，形成复杂的损伤机制。蠕变损伤在保载阶段累积，而疲劳损伤主要在加载-卸载过渡阶段产生。

氧化与环境影响：高温环境下的氧化作用会加速表面裂纹的萌生与扩展，环境因素对热疲劳寿命有显著影响。

1.2 热疲劳损伤机理

热疲劳损伤的物理过程可分为三个主要阶段：

第一阶段：裂纹萌生

在循环热应力作用下，材料表面或近表面的应力集中区域（如夹杂物、晶界、相界）首先发生塑性变形累积。随着循环次数增加，位错密度增加，形成 Persistent Slip Bands (PSB)，最终在材料表面形成挤入-挤出特征，萌生微裂纹。

裂纹萌生寿命 $N_i$ 通常占总寿命的较大比例，对于光滑试件可达总寿命的 $50\%$ 以上。裂纹萌生寿命与应变幅 $\Delta \epsilon$ 的关系可用Coffin-Manson方程描述：

$$\frac{\Delta {\epsilon }_p}{2}={\epsilon }_f^{\prime }(2N_i{)}^c$$

其中，$\Delta {\epsilon }_p$ 为塑性应变幅，${\epsilon }_f^{\prime }$ 为疲劳延性系数，$c$ 为疲劳延性指数（通常为 $-0.5$ 至 $-0.7$）。

第二阶段：裂纹扩展

微裂纹萌生后，在循环热应力作用下沿晶界或穿晶扩展。裂纹扩展速率 $da/dN$ 与应力强度因子幅 $\Delta K$ 的关系可用Paris公式描述：

$$\frac{da}{dN}=C(\Delta K{)}^m$$

其中，$C$ 和 $m$ 为材料常数。在热疲劳条件下，由于温度变化导致材料断裂韧性变化，裂纹扩展行为呈现复杂的温度依赖性。

第三阶段：快速断裂

当裂纹扩展至临界尺寸时，剩余截面无法承受外载荷，发生快速断裂。临界裂纹尺寸 $a_c$ 可由断裂韧性 $K_{IC}$ 和最大应力 ${\sigma }_{max}$ 确定：

$$a_c=\frac{1}{\pi }{\left(\frac{K_{IC}}{Y{\sigma }_{max}}\right)}^2$$

其中，$Y$ 为几何修正因子。

1.3 热应力循环特性

热疲劳中的应力-应变循环呈现复杂的非线性特征：

应力-应变迟滞回线：在循环热载荷作用下，材料的应力-应变响应形成闭合的迟滞回线。回线面积代表每个循环中单位体积材料耗散的塑性应变能，与疲劳损伤直接相关。

循环软化与硬化：材料在循环载荷作用下可能表现出循环软化（应力幅随循环次数增加而降低）或循环硬化（应力幅增加）行为。这种循环稳定性对热疲劳寿命预测有重要影响。

平均应力效应：热疲劳循环中可能存在非零平均应力，拉伸平均应力会加速疲劳损伤，而压缩平均应力可能延长寿命。平均应力效应可通过多种模型考虑，如Goodman图、Gerber抛物线或Smith-Watson-Topper(SWT)参数。

---

2. 热疲劳寿命预测理论

2.1 应变-寿命法（ε-N法）

应变-寿命法是目前工程上应用最广泛的热疲劳寿命预测方法，基于局部应变概念，认为疲劳裂纹萌生寿命由缺口根部的局部应变幅控制。

总应变-寿命方程：

Manson和Coffin分别独立提出了总应变幅与疲劳寿命的关系，后由Basquin补充弹性应变项，形成完整的总应变-寿命方程：

$$\frac{\Delta \epsilon }{2}=\frac{\Delta {\epsilon }_e}{2}+\frac{\Delta {\epsilon }_p}{2}=\frac{{\sigma }_f^{\prime }}{E}(2N_f{)}^b+{\epsilon }_f^{\prime }(2N_f{)}^c$$

其中：

• $\Delta \epsilon /2$ 为总应变幅
• ${\sigma }_f^{\prime }$ 为疲劳强度系数
• $b$ 为疲劳强度指数（通常为 $-0.05$ 至 $-0.12$）
• $E$ 为弹性模量
• ${\epsilon }_f^{\prime }$ 为疲劳延性系数
• $c$ 为疲劳延性指数（通常为 $-0.5$ 至 $-0.7$）
• $N_f$ 为疲劳寿命（循环次数）

过渡寿命：

弹性应变幅与塑性应变幅相等时的寿命称为过渡寿命 $N_t$：

$$N_t=\frac{1}{2}{\left(\frac{{\sigma }_f^{\prime }}{E{\epsilon }_f^{\prime }}\right)}^{\frac{1}{c-b}}$$

当 $N_f<N_t$ 时，塑性应变主导，为低周疲劳区；当 $N_f>N_t$ 时，弹性应变主导，为高周疲劳区。

温度修正：

对于热疲劳问题，材料参数随温度变化，需进行温度修正。常用方法包括：

1. 等效温度法：采用循环中的平均温度或加权平均温度确定材料参数
2. 分段积分法：将温度循环分为若干等温段，分别计算各段损伤后累积
3. 连续变化法：考虑材料参数的连续温度依赖性

2.2 能量法

能量法基于疲劳损伤与耗散塑性应变能相关的概念，认为每个循环中耗散的应变能是驱动疲劳裂纹萌生与扩展的物理量。

Ostergren能量模型：

Ostergren提出了基于拉伸迟滞能的寿命预测模型：

$$\Delta W={\int }_{{\epsilon }_{min}}^{{\epsilon }_{max}}\sigma d\epsilon$$

寿命预测方程为：

$$\Delta W\cdot N_f^{\alpha }=C$$

或表示为：

$$N_f={\left(\frac{C}{\Delta W}\right)}^{1/\alpha }$$

其中，$\Delta W$ 为每个循环的拉伸迟滞能，$\alpha$ 和 $C$ 为材料常数。

Smith-Watson-Topper(SWT)参数：

SWT参数综合考虑了应变幅和最大应力的影响：

$$P_{SWT}={\sigma }_{max}\cdot \frac{\Delta \epsilon }{2}=\frac{{\sigma }_f^{\prime 2}}{E}(2N_f{)}^{2b}+{\sigma }_f^{\prime }{\epsilon }_f^{\prime }(2N_f{)}^{b+c}$$

SWT参数特别适用于存在平均应力的工况，能够很好地关联不同平均应力水平下的疲劳数据。

2.3 损伤累积理论

线性损伤累积法则（Miner法则）：

Miner提出了线性损伤累积假设，认为各级载荷下的损伤可以线性叠加：

$$D=\sum _{i=1}^k\frac{n_i}{N_{fi}}$$

当 $D\ge 1$ 时，认为发生疲劳失效。其中，$n_i$ 为第 $i$ 级载荷的实际循环次数，$N_{fi}$ 为第 $i$ 级载荷下的疲劳寿命。

Miner法则形式简单，但存在明显局限性：

• 未考虑载荷顺序效应
• 未考虑平均应力影响
• 未考虑载荷交互作用

非线性损伤累积模型：

为克服Miner法则的不足，研究者提出了多种非线性损伤累积模型：

1. Marco-Starkey模型：
   $$D=\sum _{i=1}^k{\left(\frac{n_i}{N_{fi}}\right)}^{{\alpha }_i}$$
   其中，${\alpha }_i$ 为与应力水平相关的指数。

2. Chaboche连续损伤力学模型：
   $$\frac{dD}{dN}={\left[1-(1-D{)}^{1-\alpha }\right]}^{\beta }{\left(\frac{\Delta \sigma }{M(1-D)}\right)}^{\gamma }$$
   其中，$D$ 为损伤变量（$0\le D\le 1$），$\alpha$、$\beta$、$\gamma$、$M$ 为材料常数。

2.4 蠕变-疲劳交互作用模型

高温热疲劳中，蠕变损伤与疲劳损伤相互耦合。常用的交互作用模型包括：

线性累积损伤法则：

$$D_{total}=D_{fatigue}+D_{creep}=\sum _i\frac{n_i}{N_{fi}}+\sum _j\frac{t_j}{t_{rj}}$$

其中，$t_j$ 为第 $j$ 个保载时间，$t_{rj}$ 为对应应力温度条件下的蠕变断裂时间。

频率修正法：

考虑加载频率对蠕变-疲劳交互作用的影响：

$$N_f=A\cdot f^{\beta }\cdot \Delta {\epsilon }_p^{\alpha }$$

其中，$f$ 为加载频率，$\beta$ 为频率敏感指数。

应变范围划分法（SRP）：

将非弹性应变范围划分为时间相关（蠕变）和时间无关（塑性）两部分，分别计算损伤：

$$\Delta {\epsilon }_{in}=\Delta {\epsilon }_p+\Delta {\epsilon }_c$$

$$\frac{1}{N_f}=\frac{1}{N_p}+\frac{1}{N_c}$$

其中，$N_p$ 为纯疲劳寿命，$N_c$ 为纯蠕变寿命。

---

3. 热机械耦合分析方法

3.1 热-结构耦合方程

热疲劳分析需要求解热-结构耦合问题，控制方程包括：

热传导方程：

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+\dot{q}$$

其中，$\rho$ 为密度，$c_p$ 为比热容，$k$ 为导热系数，$\dot{q}$ 为内热源。

热弹性力学方程：

$$\nabla \cdot \mathbf{\sigma }+\mathbf{f}=\rho \frac{{\partial }^2\mathbf{u}}{\partial t^2}$$

热弹性本构关系：

$$\mathbf{\sigma }=\mathbf{D}:(\mathbf{\epsilon }-{\mathbf{\epsilon }}^{th})$$

$${\mathbf{\epsilon }}^{th}=\alpha (T-T_0)\mathbf{I}$$

其中，$\mathbf{D}$ 为弹性张量，$\alpha$ 为热膨胀系数，$T_0$ 为参考温度。

3.2 循环塑性本构模型

为准确描述循环载荷下的材料响应，需采用循环塑性本构模型：

Chaboche粘塑性模型：

Chaboche模型通过引入多个背应力分量描述材料的循环硬化/软化行为：

$${\dot{\mathbf{\epsilon }}}^p=\frac{3}{2}\dot{p}\frac{\mathbf{s}-\mathbf{X}}{J_2(\mathbf{s}-\mathbf{X})}$$

$${\dot{\mathbf{X}}}_i=\frac{2}{3}{C}_i{\dot{\mathbf{\epsilon }}}^p-{\gamma }_i{\mathbf{X}}_i\dot{p}$$

$$\dot{p}={⟨\frac{J_2(\mathbf{s}-\mathbf{X})-{\sigma }_y-R}{K}⟩}^n$$

其中，$\mathbf{X}$ 为背应力，$R$ 为各向同性硬化变量，$C_i$、${\gamma }_i$、$K$、$n$ 为材料参数。

---

4. Python仿真实现

4.1 热疲劳寿命预测仿真框架

"""
主题097：热疲劳寿命预测仿真
包含：热应力分析、循环塑性、损伤累积、寿命预测
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch, Circle
import matplotlib
matplotlib.use('Agg')
import warnings
warnings.filterwarnings('ignore')
import os
from scipy.optimize import fsolve
from scipy.integrate import odeint

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\传热学仿真\主题097'
os.makedirs(output_dir, exist_ok=True)

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("热疲劳寿命预测仿真系统")
print("=" * 60)

# ============================================
# 材料参数定义（以镍基高温合金为例）
# ============================================
class SuperalloyProperties:
    """镍基高温合金材料属性"""
    def __init__(self):
        # 弹性模量 (MPa) - 温度相关
        self.E_room = 210000  # 室温弹性模量
        self.E_slope = -50    # 温度系数 (MPa/°C)
        
        # 热膨胀系数 (1/°C)
        self.alpha = 1.3e-5
        
        # 导热系数 (W/m·K)
        self.k = 15
        
        # 密度 (kg/m³)
        self.rho = 8200
        
        # 比热容 (J/kg·K)
        self.cp = 420
        
        # 疲劳性能参数
        self.sigma_f_prime = 1500  # 疲劳强度系数 (MPa)
        self.b = -0.08             # 疲劳强度指数
        self.epsilon_f_prime = 0.6  # 疲劳延性系数
        self.c = -0.58              # 疲劳延性指数
        
        # 循环应力-应变曲线参数 (Ramberg-Osgood)
        self.K_prime = 1800  # 循环强度系数 (MPa)
        self.n_prime = 0.15  # 循环应变硬化指数
        
        # 屈服强度 (MPa) - 温度相关
        self.sigma_y_room = 800
        self.sigma_y_slope = -0.4  # MPa/°C
        
        # 蠕变参数 (Norton定律)
        self.A_creep = 1e-15
        self.n_creep = 5
        self.Q_creep = 250000  # 激活能 (J/mol)
        
    def E(self, T):
        """温度相关的弹性模量"""
        return max(self.E_room + self.E_slope * T, 50000)
    
    def sigma_y(self, T):
        """温度相关的屈服强度"""
        return max(self.sigma_y_room + self.sigma_y_slope * T, 200)
    
    def thermal_diffusivity(self):
        """热扩散系数"""
        return self.k / (self.rho * self.cp)

# 实例化材料
material = SuperalloyProperties()
print("\n✓ 材料参数已初始化（镍基高温合金）")

# ============================================
# 1. 温度循环与热应力分析
# ============================================
print("\n" + "=" * 60)
print("1. 温度循环与热应力分析")
print("=" * 60)

def thermal_stress_analysis():
    """
    热应力分析：计算温度循环引起的热应力
    假设：受约束的杆件，温度变化产生热应力
    """
    # 时间参数
    t_max = 1000  # 总时间 (s)
    dt = 1.0      # 时间步长
    t = np.arange(0, t_max, dt)
    n_cycles = 5  # 循环次数
    
    # 温度循环参数
    T_mean = 600      # 平均温度 (°C)
    T_amp = 400       # 温度幅值 (°C)
    freq = n_cycles / t_max  # 频率 (Hz)
    
    # 生成温度循环
    T = T_mean + T_amp * np.sin(2 * np.pi * freq * t)
    T = np.maximum(T, 100)  # 最低温度限制
    
    # 约束条件：完全约束（热应变完全转化为应力）
    # 热应力 = E * alpha * (T - T_ref)
    # 假设参考温度为平均温度
    T_ref = T_mean
    
    # 计算热应力（弹性假设）
    sigma_thermal = np.zeros_like(t)
    for i, Ti in enumerate(T):
        Ei = material.E(Ti)
        sigma_thermal[i] = Ei * material.alpha * (Ti - T_ref)
    
    # 考虑屈服限制（简单弹塑性模型）
    sigma_effective = np.zeros_like(t)
    epsilon_plastic = np.zeros_like(t)
    epsilon_plastic_cumulative = 0
    
    for i in range(len(t)):
        sigma_y_i = material.sigma_y(T[i])
        
        if abs(sigma_thermal[i]) <= sigma_y_i:
            # 弹性范围
            sigma_effective[i] = sigma_thermal[i]
        else:
            # 塑性范围
            sigma_effective[i] = np.sign(sigma_thermal[i]) * sigma_y_i
            # 计算塑性应变增量
            epsilon_elastic = sigma_y_i / material.E(T[i])
            epsilon_thermal = material.alpha * (T[i] - T_ref)
            epsilon_plastic[i] = abs(epsilon_thermal) - epsilon_elastic
            epsilon_plastic_cumulative += epsilon_plastic[i]
    
    # 绘制结果
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：温度循环
    ax1 = axes[0, 0]
    ax1.plot(t, T, 'r-', linewidth=1.5, label='Temperature')
    ax1.axhline(y=T_mean, color='gray', linestyle='--', alpha=0.5, label='Mean Temp')
    ax1.fill_between(t, T_mean - T_amp, T_mean + T_amp, alpha=0.1, color='red')
    ax1.set_xlabel('Time (s)', fontsize=11)
    ax1.set_ylabel('Temperature (°C)', fontsize=11)
    ax1.set_title('Thermal Cycle', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(0, t_max)
    
    # 图2：热应力演变
    ax2 = axes[0, 1]
    ax2.plot(t, sigma_thermal/1e3, 'b--', linewidth=1, alpha=0.7, label='Elastic Stress')
    ax2.plot(t, sigma_effective/1e3, 'b-', linewidth=1.5, label='Actual Stress')
    sigma_y_curve = [material.sigma_y(Ti)/1e3 for Ti in T]
    ax2.plot(t, sigma_y_curve, 'g:', linewidth=1.5, label='Yield Strength')
    ax2.plot(t, -np.array(sigma_y_curve), 'g:', linewidth=1.5)
    ax2.set_xlabel('Time (s)', fontsize=11)
    ax2.set_ylabel('Stress (MPa)', fontsize=11)
    ax2.set_title('Thermal Stress Evolution', fontsize=12, fontweight='bold')
    ax2.legend(loc='upper right')
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, t_max)
    
    # 图3：应力-应变迟滞回线
    ax3 = axes[1, 0]
    epsilon_total = material.alpha * (T - T_ref)
    epsilon_elastic = sigma_effective / np.array([material.E(Ti) for Ti in T])
    
    # 绘制多个循环的迟滞回线
    cycle_period = int(1 / freq / dt)
    for cycle in range(min(5, n_cycles)):
        start_idx = cycle * cycle_period
        end_idx = min((cycle + 1) * cycle_period, len(t))
        color_intensity = 0.3 + 0.7 * cycle / min(5, n_cycles)
        ax3.plot(epsilon_total[start_idx:end_idx]*100, 
                sigma_effective[start_idx:end_idx]/1e3, 
                color=(0, 0, color_intensity), linewidth=1.5, 
                label=f'Cycle {cycle+1}' if cycle < 3 else '')
    
    ax3.set_xlabel('Strain (%)', fontsize=11)
    ax3.set_ylabel('Stress (MPa)', fontsize=11)
    ax3.set_title('Stress-Strain Hysteresis Loops', fontsize=12, fontweight='bold')
    ax3.legend(loc='upper right')
    ax3.grid(True, alpha=0.3)
    
    # 图4：累积塑性应变
    ax4 = axes[1, 1]
    cum_plastic = np.cumsum(epsilon_plastic)
    ax4.plot(t, cum_plastic * 100, 'purple', linewidth=2)
    ax4.set_xlabel('Time (s)', fontsize=11)
    ax4.set_ylabel('Cumulative Plastic Strain (%)', fontsize=11)
    ax4.set_title('Cumulative Plastic Strain', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    ax4.set_xlim(0, t_max)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/thermal_stress_analysis.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 热应力分析完成")
    print(f"  - 温度范围: {T.min():.1f}°C ~ {T.max():.1f}°C")
    print(f"  - 最大热应力: {np.max(np.abs(sigma_thermal))/1e3:.1f} MPa")
    print(f"  - 累积塑性应变: {cum_plastic[-1]*100:.4f}%")
    
    return t, T, sigma_effective, epsilon_total, epsilon_plastic

t, T, sigma, epsilon_total, epsilon_plastic = thermal_stress_analysis()

# ============================================
# 2. 应变-寿命曲线（Coffin-Manson）
# ============================================
print("\n" + "=" * 60)
print("2. 应变-寿命曲线分析")
print("=" * 60)

def strain_life_analysis():
    """
    基于Coffin-Manson方程的应变-寿命分析
    """
    # 寿命范围
    Nf = np.logspace(1, 7, 500)  # 10 ~ 10^7 循环
    
    # 温度（假设为平均温度）
    T_avg = 600  # °C
    E_avg = material.E(T_avg)
    
    # 材料参数
    sigma_f_prime = material.sigma_f_prime
    b = material.b
    epsilon_f_prime = material.epsilon_f_prime
    c = material.c
    
    # 计算应变分量
    elastic_strain = (sigma_f_prime / E_avg) * (2 * Nf) ** b
    plastic_strain = epsilon_f_prime * (2 * Nf) ** c
    total_strain = elastic_strain + plastic_strain
    
    # 过渡寿命
    N_t = 0.5 * (sigma_f_prime / (E_avg * epsilon_f_prime)) ** (1 / (c - b))
    epsilon_t = 2 * (sigma_f_prime / E_avg) * (2 * N_t) ** b
    
    # 创建图形
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 图1：应变-寿命曲线（双对数坐标）
    ax1 = axes[0]
    ax1.loglog(Nf, elastic_strain * 100, 'b-', linewidth=2, label='Elastic Strain')
    ax1.loglog(Nf, plastic_strain * 100, 'r-', linewidth=2, label='Plastic Strain')
    ax1.loglog(Nf, total_strain * 100, 'k-', linewidth=2.5, label='Total Strain')
    
    # 标记过渡寿命点
    ax1.scatter([N_t], [epsilon_t * 100], color='green', s=100, zorder=5, 
               label=f'Transition Life (Nt={N_t:.0f})')
    ax1.axvline(x=N_t, color='green', linestyle='--', alpha=0.5)
    
    # 标注区域
    ax1.axvspan(10, N_t, alpha=0.1, color='red', label='Low Cycle Fatigue')
    ax1.axvspan(N_t, 1e7, alpha=0.1, color='blue', label='High Cycle Fatigue')
    
    ax1.set_xlabel('Fatigue Life Nf (cycles)', fontsize=11)
    ax1.set_ylabel('Strain Amplitude (%)', fontsize=11)
    ax1.set_title('Strain-Life Curve (Coffin-Manson)', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, which="both", ls="-", alpha=0.3)
    ax1.set_xlim(10, 1e7)
    ax1.set_ylim(0.01, 10)
    
    # 图2：不同温度下的应变-寿命曲线
    ax2 = axes[1]
    temperatures = [20, 400, 600, 800]  # °C
    colors = ['blue', 'green', 'orange', 'red']
    
    for T, color in zip(temperatures, colors):
        E_T = material.E(T)
        # 假设疲劳参数随温度变化
        sigma_f_T = sigma_f_prime * (1 - 0.0005 * T)
        epsilon_f_T = epsilon_f_prime * (1 - 0.0003 * T)
        
        elastic_T = (sigma_f_T / E_T) * (2 * Nf) ** b
        plastic_T = epsilon_f_T * (2 * Nf) ** c
        total_T = elastic_T + plastic_T
        
        ax2.loglog(Nf, total_T * 100, color=color, linewidth=2, 
                  label=f'T = {T}°C')
    
    ax2.set_xlabel('Fatigue Life Nf (cycles)', fontsize=11)
    ax2.set_ylabel('Strain Amplitude (%)', fontsize=11)
    ax2.set_title('Temperature Effect on Strain-Life', fontsize=12, fontweight='bold')
    ax2.legend(loc='upper right')
    ax2.grid(True, which="both", ls="-", alpha=0.3)
    ax2.set_xlim(10, 1e7)
    ax2.set_ylim(0.01, 10)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/strain_life_curve.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 应变-寿命曲线分析完成")
    print(f"  - 过渡寿命 Nt = {N_t:.0f} 循环")
    print(f"  - 低周疲劳区 (N < Nt): 塑性应变主导")
    print(f"  - 高周疲劳区 (N > Nt): 弹性应变主导")
    
    return Nf, elastic_strain, plastic_strain, total_strain, N_t

Nf, elastic_strain, plastic_strain, total_strain, N_t = strain_life_analysis()

# ============================================
# 3. 循环应力-应变响应
# ============================================
print("\n" + "=" * 60)
print("3. 循环应力-应变响应分析")
print("=" * 60)

def cyclic_stress_strain_response():
    """
    分析循环载荷下的应力-应变响应
    包括循环硬化/软化行为
    """
    # 应变幅值范围
    strain_amplitudes = [0.002, 0.005, 0.01, 0.02]  # 应变幅
    n_cycles = 100
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    for idx, delta_epsilon in enumerate(strain_amplitudes):
        ax = axes[idx // 2, idx % 2]
        
        # 简化的循环应力-应变响应模型
        # 假设材料表现出循环硬化
        cycles = np.arange(1, n_cycles + 1)
        
        # 循环硬化模型
        sigma_sat = material.K_prime * (delta_epsilon) ** material.n_prime
        sigma_initial = 0.7 * sigma_sat
        
        # 指数硬化规律
        beta = 0.1  # 硬化速率
        stress_amplitude = sigma_sat - (sigma_sat - sigma_initial) * np.exp(-beta * cycles)
        
        # 绘制应力幅随循环次数变化
        ax.plot(cycles, stress_amplitude / 1e3, 'b-', linewidth=2)
        ax.axhline(y=sigma_sat / 1e3, color='r', linestyle='--', 
                  label=f'Saturation: {sigma_sat/1e3:.0f} MPa')
        
        ax.set_xlabel('Number of Cycles', fontsize=10)
        ax.set_ylabel('Stress Amplitude (MPa)', fontsize=10)
        ax.set_title(f'Strain Amplitude = {delta_epsilon*100:.2f}%', fontsize=11, fontweight='bold')
        ax.legend(loc='lower right')
        ax.grid(True, alpha=0.3)
        ax.set_xlim(1, n_cycles)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/cyclic_response.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    # 循环应力-应变曲线
    fig, ax = plt.subplots(figsize=(10, 7))
    
    strain_amps = np.logspace(-4, -1, 100)
    
    # 单调加载曲线
    sigma_monotonic = material.K_prime * strain_amps ** material.n_prime
    
    # 循环曲线（通常比单调曲线低或高，取决于硬化/软化）
    # 假设循环硬化
    sigma_cyclic = 0.9 * material.K_prime * strain_amps ** material.n_prime
    
    ax.loglog(strain_amps * 100, sigma_monotonic / 1e3, 'b-', linewidth=2, 
             label='Monotonic Curve')
    ax.loglog(strain_amps * 100, sigma_cyclic / 1e3, 'r-', linewidth=2, 
             label='Cyclic Curve')
    
    ax.set_xlabel('Strain Amplitude (%)', fontsize=11)
    ax.set_ylabel('Stress Amplitude (MPa)', fontsize=11)
    ax.set_title('Cyclic Stress-Strain Curve', fontsize=12, fontweight='bold')
    ax.legend(loc='lower right')
    ax.grid(True, which="both", ls="-", alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/cyclic_stress_strain_curve.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 循环应力-应变响应分析完成")
    print(f"  - 分析了 {len(strain_amplitudes)} 个应变幅值水平")
    print(f"  - 材料表现出循环硬化特性")

cyclic_stress_strain_response()

# ============================================
# 4. 损伤累积与寿命预测
# ============================================
print("\n" + "=" * 60)
print("4. 损伤累积与寿命预测")
print("=" * 60)

def damage_accumulation_analysis():
    """
    损伤累积分析与寿命预测
    使用Miner法则和连续损伤力学方法
    """
    # 多级载荷谱
    load_levels = [
        {'strain_amp': 0.01, 'cycles': 500},
        {'strain_amp': 0.008, 'cycles': 1000},
        {'strain_amp': 0.006, 'cycles': 2000},
        {'strain_amp': 0.004, 'cycles': 5000},
    ]
    
    T_avg = 600
    E_avg = material.E(T_avg)
    
    # 计算各级载荷下的寿命
    for level in load_levels:
        delta_epsilon = level['strain_amp']
        
        # 使用总应变-寿命方程求解寿命
        def strain_life_eq(N):
            elastic = (material.sigma_f_prime / E_avg) * (2 * N) ** material.b
            plastic = material.epsilon_f_prime * (2 * N) ** material.c
            return elastic + plastic - delta_epsilon
        
        N_f = fsolve(strain_life_eq, 1000)[0]
        level['N_f'] = max(N_f, 1)
        level['damage'] = level['cycles'] / level['N_f']
    
    # 累积损伤
    cumulative_damage = np.cumsum([level['damage'] for level in load_levels])
    
    # 创建图形
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：载荷谱
    ax1 = axes[0, 0]
    level_names = [f'Level {i+1}' for i in range(len(load_levels))]
    strain_amps = [level['strain_amp'] * 100 for level in load_levels]
    cycles = [level['cycles'] for level in load_levels]
    
    x_pos = np.arange(len(load_levels))
    bars = ax1.bar(x_pos, strain_amps, color='steelblue', alpha=0.7)
    ax1.set_xlabel('Load Level', fontsize=11)
    ax1.set_ylabel('Strain Amplitude (%)', fontsize=11)
    ax1.set_title('Multi-Level Load Spectrum', fontsize=12, fontweight='bold')
    ax1.set_xticks(x_pos)
    ax1.set_xticklabels(level_names)
    ax1.grid(True, alpha=0.3, axis='y')
    
    # 添加循环次数标注
    for i, (bar, cyc) in enumerate(zip(bars, cycles)):
        ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.05,
                f'{cyc} cycles', ha='center', va='bottom', fontsize=9)
    
    # 图2：各级载荷损伤
    ax2 = axes[0, 1]
    damages = [level['damage'] for level in load_levels]
    colors = ['green' if d < 0.1 else 'orange' if d < 0.3 else 'red' for d in damages]
    bars2 = ax2.bar(x_pos, damages, color=colors, alpha=0.7)
    ax2.set_xlabel('Load Level', fontsize=11)
    ax2.set_ylabel('Damage per Level', fontsize=11)
    ax2.set_title('Damage at Each Load Level', fontsize=12, fontweight='bold')
    ax2.set_xticks(x_pos)
    ax2.set_xticklabels(level_names)
    ax2.axhline(y=1.0, color='red', linestyle='--', linewidth=2, label='Failure')
    ax2.legend()
    ax2.grid(True, alpha=0.3, axis='y')
    
    # 图3：累积损伤
    ax3 = axes[1, 0]
    ax3.plot(range(1, len(load_levels) + 1), cumulative_damage, 'bo-', 
            linewidth=2, markersize=8)
    ax3.axhline(y=1.0, color='red', linestyle='--', linewidth=2, label='Failure Threshold')
    ax3.fill_between(range(1, len(load_levels) + 1), 0, cumulative_damage, alpha=0.3)
    ax3.set_xlabel('Load Level', fontsize=11)
    ax3.set_ylabel('Cumulative Damage', fontsize=11)
    ax3.set_title('Cumulative Damage (Miner Rule)', fontsize=12, fontweight='bold')
    ax3.set_xticks(range(1, len(load_levels) + 1))
    ax3.set_xticklabels(level_names)
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_ylim(0, max(1.2, max(cumulative_damage) * 1.1))
    
    # 图4：寿命预测对比
    ax4 = axes[1, 1]
    
    # 不同应变幅下的预测寿命
    strain_range = np.logspace(-3, -1, 50)
    lives_predicted = []
    
    for eps in strain_range:
        def eq(N):
            elastic = (material.sigma_f_prime / E_avg) * (2 * N) ** material.b
            plastic = material.epsilon_f_prime * (2 * N) ** material.c
            return elastic + plastic - eps
        
        N_pred = fsolve(eq, 1000)[0]
        lives_predicted.append(max(N_pred, 1))
    
    ax4.loglog(strain_range * 100, lives_predicted, 'b-', linewidth=2, label='Predicted Life')
    
    # 标记实际载荷点
    for level in load_levels:
        ax4.scatter([level['strain_amp'] * 100], [level['N_f']], 
                  s=100, color='red', zorder=5)
        ax4.annotate(f"N={level['N_f']:.0f}", 
                    xy=(level['strain_amp'] * 100, level['N_f']),
                    xytext=(10, 10), textcoords='offset points', fontsize=9)
    
    ax4.set_xlabel('Strain Amplitude (%)', fontsize=11)
    ax4.set_ylabel('Predicted Life (cycles)', fontsize=11)
    ax4.set_title('Life Prediction vs Strain Amplitude', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, which="both", ls="-", alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/damage_accumulation.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 损伤累积分析完成")
    print(f"  - 总累积损伤: {cumulative_damage[-1]:.4f}")
    print(f"  - 剩余寿命比例: {max(0, 1 - cumulative_damage[-1]):.2%}")
    
    return load_levels, cumulative_damage

load_levels, cum_damage = damage_accumulation_analysis()

# ============================================
# 5. 平均应力效应与SWT参数
# ============================================
print("\n" + "=" * 60)
print("5. 平均应力效应与SWT参数")
print("=" * 60)

def mean_stress_effect_analysis():
    """
    分析平均应力对疲劳寿命的影响
    使用Smith-Watson-Topper(SWT)参数
    """
    # 应变幅值
    strain_amps = np.array([0.002, 0.004, 0.006, 0.008, 0.01])
    
    # 平均应力范围
    mean_stresses = np.array([0, 100, 200, 300])  # MPa
    
    T_avg = 600
    E_avg = material.E(T_avg)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 图1：不同平均应力下的寿命
    ax1 = axes[0]
    colors = ['blue', 'green', 'orange', 'red']
    
    for sigma_m, color in zip(mean_stresses, colors):
        lives = []
        for delta_eps in strain_amps:
            # SWT参数
            # 近似计算最大应力
            sigma_max = sigma_m + E_avg * delta_eps / 2
            sigma_max = min(sigma_max, material.sigma_f_prime)
            
            P_SWT = sigma_max * delta_eps / 2
            
            # 由SWT参数求解寿命
            def swt_eq(N):
                term1 = (material.sigma_f_prime ** 2 / E_avg) * (2 * N) ** (2 * material.b)
                term2 = material.sigma_f_prime * material.epsilon_f_prime * (2 * N) ** (material.b + material.c)
                return (term1 + term2) - P_SWT
            
            N_pred = fsolve(swt_eq, 1000)[0]
            lives.append(max(N_pred, 10))
        
        ax1.loglog(strain_amps * 100, lives, 'o-', color=color, linewidth=2, 
                  markersize=6, label=f'σm = {sigma_m} MPa')
    
    ax1.set_xlabel('Strain Amplitude (%)', fontsize=11)
    ax1.set_ylabel('Fatigue Life (cycles)', fontsize=11)
    ax1.set_title('Mean Stress Effect on Fatigue Life', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, which="both", ls="-", alpha=0.3)
    
    # 图2：SWT参数与寿命关系
    ax2 = axes[1]
    
    # 生成SWT参数范围
    P_SWT_range = np.logspace(1, 5, 100)  # MPa * strain
    
    # 由SWT参数计算寿命
    def swt_to_life(P_SWT):
        def eq(N):
            term1 = (material.sigma_f_prime ** 2 / E_avg) * (2 * N) ** (2 * material.b)
            term2 = material.sigma_f_prime * material.epsilon_f_prime * (2 * N) ** (material.b + material.c)
            return (term1 + term2) - P_SWT
        
        return fsolve(eq, 1000)[0]
    
    lives_swt = [swt_to_life(P) for P in P_SWT_range]
    
    ax2.loglog(P_SWT_range, lives_swt, 'b-', linewidth=2, label='SWT Life Prediction')
    
    # 标记一些参考点
    ref_points = [100, 500, 1000, 5000]
    for P in ref_points:
        N_ref = swt_to_life(P)
        ax2.scatter([P], [N_ref], s=80, color='red', zorder=5)
        ax2.annotate(f'N={N_ref:.0f}', xy=(P, N_ref), 
                    xytext=(10, 5), textcoords='offset points', fontsize=9)
    
    ax2.set_xlabel('SWT Parameter (MPa)', fontsize=11)
    ax2.set_ylabel('Fatigue Life (cycles)', fontsize=11)
    ax2.set_title('Smith-Watson-Topper Parameter', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, which="both", ls="-", alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/mean_stress_effect.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 平均应力效应分析完成")
    print(f"  - 分析了 {len(mean_stresses)} 个平均应力水平")
    print(f"  - 拉伸平均应力显著降低疲劳寿命")

mean_stress_effect_analysis()

# ============================================
# 6. 蠕变-疲劳交互作用
# ============================================
print("\n" + "=" * 60)
print("6. 蠕变-疲劳交互作用分析")
print("=" * 60)

def creep_fatigue_interaction():
    """
    蠕变-疲劳交互作用分析
    """
    # 温度
    T_celsius = 700
    T_kelvin = T_celsius + 273.15
    
    # 应力水平
    stresses = np.linspace(200, 600, 100)  # MPa
    
    # 计算纯蠕变寿命（Larson-Miller参数）
    # LMP = T * (log(tr) + C)
    C_larson = 20
    LMP = 22000  # 典型值
    
    creep_lives = []
    for sigma in stresses:
        # 简化的Larson-Miller关系
        LMP_sigma = LMP - 5 * sigma  # 简化的线性关系
        log_tr = LMP_sigma / T_kelvin - C_larson
        tr = 10 ** log_tr
        creep_lives.append(max(tr, 0.1))
    
    # 计算纯疲劳寿命（基于应变）
    E_T = material.E(T_celsius)
    fatigue_lives = []
    for sigma in stresses:
        delta_epsilon = 2 * sigma / E_T  # 假设完全反向加载
        
        def eq(N):
            elastic = (material.sigma_f_prime / E_T) * (2 * N) ** material.b
            plastic = material.epsilon_f_prime * (2 * N) ** material.c
            return elastic + plastic - delta_epsilon
        
        N_f = fsolve(eq, 1000)[0]
        fatigue_lives.append(max(N_f, 1))
    
    # 交互作用图（蠕变-疲劳交互图）
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 图1：寿命-应力曲线
    ax1 = axes[0]
    ax1.semilogy(stresses, creep_lives, 'r-', linewidth=2, label='Creep Life')
    ax1.semilogy(stresses, fatigue_lives, 'b-', linewidth=2, label='Fatigue Life')
    
    # 找到交点（临界应力）
    diff = np.abs(np.log10(creep_lives) - np.log10(fatigue_lives))
    critical_idx = np.argmin(diff)
    critical_stress = stresses[critical_idx]
    
    ax1.axvline(x=critical_stress, color='green', linestyle='--', 
               label=f'Critical Stress: {critical_stress:.0f} MPa')
    ax1.scatter([critical_stress], [creep_lives[critical_idx]], 
               s=150, color='green', zorder=5)
    
    ax1.set_xlabel('Stress (MPa)', fontsize=11)
    ax1.set_ylabel('Life (hours or cycles)', fontsize=11)
    ax1.set_title('Creep vs Fatigue Life at T=700°C', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, alpha=0.3)
    
    # 图2：蠕变-疲劳交互图
    ax2 = axes[1]
    
    # 归一化损伤
    # 假设不同保载时间和循环次数的组合
    fatigue_ratios = np.linspace(0, 1, 100)
    
    # 线性累积损伤法则
    creep_linear = 1 - fatigue_ratios
    
    # 非线性交互（考虑交互作用增强）
    # 使用椭圆准则
    creep_nonlinear = np.sqrt(1 - fatigue_ratios ** 2)
    
    # 考虑交互作用增强
    interaction_factor = 0.8  # < 1 表示交互作用加速失效
    creep_interaction = interaction_factor * np.sqrt(1 - (fatigue_ratios/interaction_factor) ** 2)
    creep_interaction = np.maximum(creep_interaction, 0)
    
    ax2.plot(fatigue_ratios, creep_linear, 'b--', linewidth=2, label='Linear (Miner)')
    ax2.plot(fatigue_ratios, creep_nonlinear, 'g-', linewidth=2, label='Elliptical')
    ax2.fill_between(fatigue_ratios, 0, creep_nonlinear, alpha=0.2, color='green')
    
    ax2.set_xlabel('Fatigue Damage Ratio (Df)', fontsize=11)
    ax2.set_ylabel('Creep Damage Ratio (Dc)', fontsize=11)
    ax2.set_title('Creep-Fatigue Interaction Diagram', fontsize=12, fontweight='bold')
    ax2.set_xlim(0, 1)
    ax2.set_ylim(0, 1)
    ax2.legend(loc='upper right')
    ax2.grid(True, alpha=0.3)
    ax2.set_aspect('equal')
    
    # 添加安全区域标注
    ax2.text(0.3, 0.3, 'Safe Region', fontsize=12, ha='center', 
            bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.5))
    ax2.text(0.7, 0.7, 'Failure Region', fontsize=12, ha='center',
            bbox=dict(boxstyle='round', facecolor='lightcoral', alpha=0.5))
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/creep_fatigue_interaction.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 蠕变-疲劳交互作用分析完成")
    print(f"  - 临界应力: {critical_stress:.0f} MPa")
    print(f"  - 低于临界应力: 蠕变主导")
    print(f"  - 高于临界应力: 疲劳主导")

creep_fatigue_interaction()

# ============================================
# 7. 热疲劳寿命预测综合案例
# ============================================
print("\n" + "=" * 60)
print("7. 热疲劳寿命预测综合案例")
print("=" * 60)

def comprehensive_life_prediction():
    """
    综合热疲劳寿命预测案例
    模拟涡轮叶片的热疲劳问题
    """
    # 涡轮叶片简化模型
    # 假设：叶片在启动-停机循环中经历温度变化
    
    # 循环参数
    n_startup = 1000  # 启动-停机循环次数
    t_cycle = 3600    # 每个循环周期 (s)
    
    # 温度循环
    T_min = 300   # 最低温度 (°C) - 停机
    T_max = 900   # 最高温度 (°C) - 满负荷
    T_mean = (T_min + T_max) / 2
    T_amp = (T_max - T_min) / 2
    
    # 时间数组
    t = np.linspace(0, t_cycle, 1000)
    
    # 温度曲线（简化的启动-停机循环）
    # 启动阶段
    T = T_mean + T_amp * np.sin(np.pi * t / t_cycle)
    
    # 材料参数随温度变化
    E_values = np.array([material.E(Ti) for Ti in T])
    sigma_y_values = np.array([material.sigma_y(Ti) for Ti in T])
    
    # 热应力计算（假设约束系数为0.8）
    constraint_factor = 0.8
    T_ref = T_min
    thermal_strain = material.alpha * (T - T_ref)
    elastic_stress = constraint_factor * E_values * thermal_strain
    
    # 弹塑性应力计算
    actual_stress = np.clip(elastic_stress, -sigma_y_values, sigma_y_values)
    
    # 计算应变幅（用于寿命预测）
    delta_epsilon = np.max(thermal_strain) - np.min(thermal_strain)
    delta_epsilon_mechanical = delta_epsilon * constraint_factor
    
    # 预测寿命
    E_avg = material.E(T_mean)
    
    def life_eq(N):
        elastic = (material.sigma_f_prime / E_avg) * (2 * N) ** material.b
        plastic = material.epsilon_f_prime * (2 * N) ** material.c
        return elastic + plastic - delta_epsilon_mechanical
    
    predicted_life = fsolve(life_eq, 1000)[0]
    predicted_life = max(predicted_life, 1)
    
    # 损伤累积
    damage_per_cycle = 1 / predicted_life
    total_damage = n_startup * damage_per_cycle
    remaining_life = max(0, (1 - total_damage) * predicted_life)
    
    # 创建综合结果图
    fig = plt.figure(figsize=(16, 12))
    gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)
    
    # 图1：温度循环
    ax1 = fig.add_subplot(gs[0, 0])
    ax1.plot(t/60, T, 'r-', linewidth=2)
    ax1.set_xlabel('Time (min)', fontsize=10)
    ax1.set_ylabel('Temperature (°C)', fontsize=10)
    ax1.set_title('Temperature Cycle', fontsize=11, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    ax1.fill_between(t/60, T, alpha=0.2, color='red')
    
    # 图2：热应力演变
    ax2 = fig.add_subplot(gs[0, 1])
    ax2.plot(t/60, elastic_stress/1e3, 'b--', alpha=0.5, label='Elastic')
    ax2.plot(t/60, actual_stress/1e3, 'b-', linewidth=2, label='Actual')
    ax2.plot(t/60, sigma_y_values/1e3, 'g:', label='Yield')
    ax2.plot(t/60, -sigma_y_values/1e3, 'g:')
    ax2.set_xlabel('Time (min)', fontsize=10)
    ax2.set_ylabel('Stress (MPa)', fontsize=10)
    ax2.set_title('Thermal Stress', fontsize=11, fontweight='bold')
    ax2.legend(fontsize=8)
    ax2.grid(True, alpha=0.3)
    
    # 图3：应力-应变迟滞回线
    ax3 = fig.add_subplot(gs[0, 2])
    strain_pct = thermal_strain * constraint_factor * 100
    ax3.plot(strain_pct, actual_stress/1e3, 'purple', linewidth=2)
    ax3.fill(strain_pct, actual_stress/1e3, alpha=0.2, color='purple')
    ax3.set_xlabel('Strain (%)', fontsize=10)
    ax3.set_ylabel('Stress (MPa)', fontsize=10)
    ax3.set_title('Hysteresis Loop', fontsize=11, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 图4：材料性能温度依赖性
    ax4 = fig.add_subplot(gs[1, 0])
    T_range = np.linspace(20, 1000, 100)
    E_range = [material.E(Ti) for Ti in T_range]
    sigma_y_range = [material.sigma_y(Ti) for Ti in T_range]
    
    ax4.plot(T_range, np.array(E_range)/1e3, 'b-', linewidth=2, label='Elastic Modulus')
    ax4_twin = ax4.twinx()
    ax4_twin.plot(T_range, np.array(sigma_y_range)/1e3, 'r-', linewidth=2, label='Yield Strength')
    ax4.set_xlabel('Temperature (°C)', fontsize=10)
    ax4.set_ylabel('E (GPa)', color='blue', fontsize=10)
    ax4_twin.set_ylabel('σy (MPa)', color='red', fontsize=10)
    ax4.set_title('Temperature-Dependent Properties', fontsize=11, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    
    # 图5：寿命预测结果
    ax5 = fig.add_subplot(gs[1, 1])
    categories = ['Predicted\nLife', 'Cycles\nExperienced', 'Remaining\nLife']
    values = [predicted_life, n_startup, remaining_life]
    colors = ['green', 'orange', 'blue']
    bars = ax5.bar(categories, values, color=colors, alpha=0.7)
    ax5.set_ylabel('Cycles', fontsize=10)
    ax5.set_title('Life Prediction Results', fontsize=11, fontweight='bold')
    ax5.set_yscale('log')
    ax5.grid(True, alpha=0.3, axis='y')
    
    for bar, val in zip(bars, values):
        ax5.text(bar.get_x() + bar.get_width()/2, bar.get_height() * 1.1,
                f'{val:.0f}', ha='center', va='bottom', fontsize=9, fontweight='bold')
    
    # 图6：损伤累积
    ax6 = fig.add_subplot(gs[1, 2])
    cycles_array = np.arange(0, n_startup + 1, 10)
    damage_array = cycles_array * damage_per_cycle
    ax6.plot(cycles_array, damage_array, 'b-', linewidth=2)
    ax6.axhline(y=1.0, color='red', linestyle='--', linewidth=2, label='Failure')
    ax6.fill_between(cycles_array, 0, damage_array, alpha=0.3)
    ax6.set_xlabel('Number of Cycles', fontsize=10)
    ax6.set_ylabel('Cumulative Damage', fontsize=10)
    ax6.set_title('Damage Accumulation', fontsize=11, fontweight='bold')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    ax6.set_xlim(0, n_startup)
    
    # 图7：安全系数
    ax7 = fig.add_subplot(gs[2, 0])
    safety_factor = 1 / total_damage if total_damage > 0 else float('inf')
    sf_display = min(safety_factor, 10)
    
    colors_sf = ['red' if safety_factor < 1 else 'orange' if safety_factor < 2 else 'green']
    ax7.barh(['Safety Factor'], [sf_display], color=colors_sf, alpha=0.7, height=0.5)
    ax7.axvline(x=1, color='red', linestyle='--', linewidth=2)
    ax7.axvline(x=2, color='orange', linestyle='--', linewidth=2)
    ax7.set_xlim(0, 10)
    ax7.set_xlabel('Factor', fontsize=10)
    ax7.set_title('Safety Assessment', fontsize=11, fontweight='bold')
    ax7.text(sf_display/2, 0, f'SF = {safety_factor:.2f}', 
            ha='center', va='center', fontsize=12, fontweight='bold')
    
    # 图8：应变-寿命位置
    ax8 = fig.add_subplot(gs[2, 1])
    Nf_range = np.logspace(1, 6, 100)
    E_avg = material.E(T_mean)
    elastic_line = (material.sigma_f_prime / E_avg) * (2 * Nf_range) ** material.b
    plastic_line = material.epsilon_f_prime * (2 * Nf_range) ** material.c
    total_line = elastic_line + plastic_line
    
    ax8.loglog(Nf_range, total_line * 100, 'b-', linewidth=2, label='Strain-Life Curve')
    ax8.scatter([predicted_life], [delta_epsilon_mechanical * 100], 
               s=200, color='red', zorder=5, marker='*', 
               label=f'Operating Point (N={predicted_life:.0f})')
    ax8.set_xlabel('Fatigue Life (cycles)', fontsize=10)
    ax8.set_ylabel('Strain Amplitude (%)', fontsize=10)
    ax8.set_title('Operating Point on S-N Curve', fontsize=11, fontweight='bold')
    ax8.legend()
    ax8.grid(True, which="both", ls="-", alpha=0.3)
    
    # 图9：风险评估矩阵
    ax9 = fig.add_subplot(gs[2, 2])
    
    # 创建风险评估矩阵
    risk_levels = np.array([
        [1, 2, 3],
        [2, 3, 4],
        [3, 4, 5]
    ])
    
    im = ax9.imshow(risk_levels, cmap='RdYlGn_r', aspect='auto', vmin=1, vmax=5)
    
    # 标注
    labels = [['Low', 'Low-Med', 'Medium'],
              ['Low-Med', 'Medium', 'Med-High'],
              ['Medium', 'Med-High', 'High']]
    
    for i in range(3):
        for j in range(3):
            ax9.text(j, i, labels[i][j], ha='center', va='center', 
                    fontsize=9, fontweight='bold')
    
    # 确定当前风险等级
    if safety_factor > 3:
        risk_pos = (0, 0)
    elif safety_factor > 2:
        risk_pos = (1, 0)
    elif safety_factor > 1.5:
        risk_pos = (1, 1)
    elif safety_factor > 1:
        risk_pos = (2, 1)
    else:
        risk_pos = (2, 2)
    
    ax9.scatter([risk_pos[1]], [risk_pos[0]], s=300, color='blue', 
               marker='*', edgecolors='white', linewidths=2, zorder=5)
    
    ax9.set_xticks([0, 1, 2])
    ax9.set_yticks([0, 1, 2])
    ax9.set_xticklabels(['Low', 'Medium', 'High'])
    ax9.set_yticklabels(['Low', 'Medium', 'High'])
    ax9.set_xlabel('Damage Severity', fontsize=10)
    ax9.set_ylabel('Occurrence Frequency', fontsize=10)
    ax9.set_title('Risk Assessment Matrix', fontsize=11, fontweight='bold')
    
    plt.savefig(f'{output_dir}/comprehensive_life_prediction.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 综合寿命预测完成")
    print(f"  - 预测寿命: {predicted_life:.0f} 循环")
    print(f"  - 已服役: {n_startup} 循环")
    print(f"  - 剩余寿命: {remaining_life:.0f} 循环")
    print(f"  - 累积损伤: {total_damage:.4f}")
    print(f"  - 安全系数: {safety_factor:.2f}")
    
    return predicted_life, total_damage, safety_factor

predicted_life, total_damage, safety_factor = comprehensive_life_prediction()

# ============================================
# 8. 不同预测方法对比
# ============================================
print("\n" + "=" * 60)
print("8. 不同寿命预测方法对比")
print("=" * 60)

def prediction_methods_comparison():
    """
    对比不同寿命预测方法的结果
    """
    # 应变幅值范围
    strain_amps = np.logspace(-3, -1, 50)
    T_avg = 600
    E_avg = material.E(T_avg)
    
    lives_coffin_manson = []
    lives_swt = []
    lives_ostergren = []
    
    for delta_eps in strain_amps:
        # 1. Coffin-Manson方法
        def cm_eq(N):
            elastic = (material.sigma_f_prime / E_avg) * (2 * N) ** material.b
            plastic = material.epsilon_f_prime * (2 * N) ** material.c
            return elastic + plastic - delta_eps
        N_cm = fsolve(cm_eq, 1000)[0]
        lives_coffin_manson.append(max(N_cm, 1))
        
        # 2. SWT方法（假设R=-1，即sigma_max = E*delta_eps/2）
        sigma_max = E_avg * delta_eps / 2
        P_SWT = sigma_max * delta_eps / 2
        
        def swt_eq(N):
            term1 = (material.sigma_f_prime ** 2 / E_avg) * (2 * N) ** (2 * material.b)
            term2 = material.sigma_f_prime * material.epsilon_f_prime * (2 * N) ** (material.b + material.c)
            return (term1 + term2) - P_SWT
        N_swt = fsolve(swt_eq, 1000)[0]
        lives_swt.append(max(N_swt, 1))
        
        # 3. Ostergren能量法（简化）
        # 假设迟滞能正比于应变幅的平方
        delta_W = 0.5 * E_avg * delta_eps**2
        alpha_ost = 0.8
        C_ost = 500
        N_ost = (C_ost / delta_W) ** (1 / alpha_ost)
        lives_ostergren.append(max(N_ost, 1))
    
    # 创建对比图
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 图1：三种方法对比
    ax1 = axes[0]
    ax1.loglog(strain_amps * 100, lives_coffin_manson, 'b-', linewidth=2, 
              label='Coffin-Manson', marker='o', markersize=4, markevery=5)
    ax1.loglog(strain_amps * 100, lives_swt, 'r-', linewidth=2, 
              label='SWT Parameter', marker='s', markersize=4, markevery=5)
    ax1.loglog(strain_amps * 100, lives_ostergren, 'g-', linewidth=2, 
              label='Ostergren Energy', marker='^', markersize=4, markevery=5)
    
    ax1.set_xlabel('Strain Amplitude (%)', fontsize=11)
    ax1.set_ylabel('Predicted Life (cycles)', fontsize=11)
    ax1.set_title('Comparison of Life Prediction Methods', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, which="both", ls="-", alpha=0.3)
    ax1.set_xlim(0.1, 10)
    
    # 图2：方法间差异分析
    ax2 = axes[1]
    ratio_swt_cm = np.array(lives_swt) / np.array(lives_coffin_manson)
    ratio_ost_cm = np.array(lives_ostergren) / np.array(lives_coffin_manson)
    
    ax2.semilogx(strain_amps * 100, ratio_swt_cm, 'r-', linewidth=2, 
                label='SWT / C-M Ratio')
    ax2.semilogx(strain_amps * 100, ratio_ost_cm, 'g-', linewidth=2, 
                label='Ostergren / C-M Ratio')
    ax2.axhline(y=1.0, color='black', linestyle='--', linewidth=1.5, label='Reference')
    ax2.fill_between(strain_amps * 100, 0.5, 2.0, alpha=0.1, color='gray', 
                    label='±50% Band')
    
    ax2.set_xlabel('Strain Amplitude (%)', fontsize=11)
    ax2.set_ylabel('Life Ratio', fontsize=11)
    ax2.set_title('Method-to-Method Variation', fontsize=12, fontweight='bold')
    ax2.legend(loc='upper right')
    ax2.grid(True, alpha=0.3)
    ax2.set_ylim(0, 3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/methods_comparison.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 预测方法对比完成")
    print(f"  - 对比了3种主要预测方法")
    print(f"  - 各方法在低/高周区预测结果存在差异")
    print(f"  - 建议根据实际工况选择合适方法")

prediction_methods_comparison()

print("\n" + "=" * 60)
print("热疲劳寿命预测仿真全部完成")
print("=" * 60)
print(f"\n生成的图表文件：")
print(f"  1. thermal_stress_analysis.png - 热应力分析")
print(f"  2. strain_life_curve.png - 应变-寿命曲线")
print(f"  3. cyclic_response.png - 循环应力响应")
print(f"  4. cyclic_stress_strain_curve.png - 循环应力-应变曲线")
print(f"  5. damage_accumulation.png - 损伤累积分析")
print(f"  6. mean_stress_effect.png - 平均应力效应")
print(f"  7. creep_fatigue_interaction.png - 蠕变-疲劳交互")
print(f"  8. comprehensive_life_prediction.png - 综合寿命预测")
print(f"  9. methods_comparison.png - 方法对比")
print(f"\n输出目录: {output_dir}")

4.2 仿真2：多轴热疲劳分析

"""
仿真2：多轴热疲劳分析
考虑复杂应力状态下的热疲劳寿命预测
"""

fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 多轴应力状态
sigma_x = np.linspace(0, 400, 50)  # MPa
sigma_y = np.linspace(0, 400, 50)
SIGMA_X, SIGMA_Y = np.meshgrid(sigma_x, sigma_y)

# 等效应力（von Mises）
sigma_vm = np.sqrt(SIGMA_X**2 + SIGMA_Y**2 - SIGMA_X*SIGMA_Y)

# 基于von Mises应力的寿命预测
def life_from_stress(sigma_eq, sigma_f=800, b=-0.08):
    """从等效应力预测寿命"""
    Nf = 0.5 * (sigma_eq / sigma_f) ** (1/b)
    return np.clip(Nf, 1, 1e7)

Nf_biaxial = life_from_stress(sigma_vm)

# 绘制寿命等高线图
im1 = axes[0, 0].contourf(SIGMA_X, SIGMA_Y, np.log10(Nf_biaxial), levels=20, cmap='jet')
axes[0, 0].set_xlabel('σx (MPa)', fontsize=11)
axes[0, 0].set_ylabel('σy (MPa)', fontsize=11)
axes[0, 0].set_title('双轴应力状态下的疲劳寿命', fontsize=12)
plt.colorbar(im1, ax=axes[0, 0], label='log10(Nf)')

# 不同应力比的影响
stress_ratios = [0, 0.5, 1.0, -0.5, -1.0]
sigma_max = 300  # MPa
colors_sr = ['blue', 'green', 'red', 'purple', 'orange']

for sr, color in zip(stress_ratios, colors_sr):
    sigma_a = sigma_max * (1 - sr) / 2  # 应力幅
    sigma_m = sigma_max * (1 + sr) / 2  # 平均应力
    
    # Goodman修正
    sigma_ar = sigma_a / (1 - sigma_m/1200)  # 修正后的应力幅
    Nf_sr = life_from_stress(sigma_ar)
    
    axes[0, 1].semilogy([sr], [Nf_sr], 'o', color=color, markersize=10, 
                       label=f'R={sr}')

axes[0, 1].set_xlabel('应力比 R', fontsize=11)
axes[0, 1].set_ylabel('疲劳寿命 (循环次数)', fontsize=11)
axes[0, 1].set_title('应力比对疲劳寿命的影响', fontsize=12)
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 非比例加载效应
phase_angles = np.linspace(0, 90, 19)  # 相位角
Nf_proportional = life_from_stress(250)  # 比例加载寿命

# 非比例硬化因子
f_np = 1 + 0.3 * np.sin(np.radians(phase_angles)) ** 2
Nf_nonprop = Nf_proportional / f_np

axes[1, 0].plot(phase_angles, Nf_nonprop/Nf_proportional, 'b-', linewidth=2)
axes[1, 0].set_xlabel('相位角 (°)', fontsize=11)
axes[1, 0].set_ylabel('寿命比 (Nf/Nf_比例)', fontsize=11)
axes[1, 0].set_title('非比例加载效应', fontsize=12)
axes[1, 0].grid(True, alpha=0.3)

# 临界平面法
angles_plane = np.linspace(0, 180, 37)  # 平面角度
normal_stress = 200 * np.cos(np.radians(angles_plane)) ** 2
shear_stress = 100 * np.sin(np.radians(2*angles_plane))

# Fatemi-Socie损伤参数
gamma_max = shear_stress / 80000  # 最大剪应变
sigma_n_max = normal_stress  # 最大正应力
k_fs = 0.3
DP_fs = gamma_max * (1 + k_fs * sigma_n_max / 300)

axes[1, 1].plot(angles_plane, DP_fs, 'r-', linewidth=2)
axes[1, 1].set_xlabel('平面角度 (°)', fontsize=11)
axes[1, 1].set_ylabel('Fatemi-Socie损伤参数', fontsize=11)
axes[1, 1].set_title('临界平面分析', fontsize=12)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/multiaxial_thermal_fatigue.png', dpi=150, bbox_inches='tight')
plt.close()
print("多轴热疲劳分析图已保存")

print("\n多轴热疲劳分析结果:")
print(f"比例加载寿命: {Nf_proportional:.0f} 次")
print(f"90°相位角寿命: {Nf_nonprop[-1]:.0f} 次")
print(f"非比例加载损伤因子: {f_np[-1]:.2f}")

4.3 仿真结果分析

上述Python仿真实现了完整的热疲劳寿命预测流程，主要包含以下分析模块：

热应力分析模块：模拟温度循环引起的热应力演变，考虑材料弹塑性行为和温度相关的力学性能。应力-应变迟滞回线直观展示了每个循环的能量耗散。

应变-寿命分析模块：基于Coffin-Manson方程，绘制总应变-寿命曲线，区分弹性应变主导的高周疲劳区和塑性应变主导的低周疲劳区。同时分析温度对疲劳性能的影响。

循环响应分析模块：模拟材料在循环载荷下的硬化/软化行为，展示应力幅随循环次数的演变规律。

损伤累积模块：基于Miner线性损伤累积法则，计算多级载荷下的累积损伤，评估结构剩余寿命。

平均应力效应模块：使用SWT参数分析平均应力对疲劳寿命的影响，展示不同平均应力水平下的寿命预测结果。

蠕变-疲劳交互模块：在高温条件下分析蠕变损伤与疲劳损伤的交互作用，绘制交互作用图确定临界应力水平。

综合预测模块：以涡轮叶片为例，展示完整的热疲劳寿命预测流程，包括温度循环、应力分析、寿命预测、损伤累积和安全评估。

---

5. 工程应用与案例分析

5.1 航空发动机涡轮叶片热疲劳

航空发动机涡轮叶片是典型的高温热疲劳部件，工作条件极为苛刻：

工作条件：

• 温度范围：600°C ~ 1050°C
• 离心应力：200 ~ 400 MPa
• 气动弯曲应力：±100 MPa
• 启动-停机循环：数千次

热疲劳特征：

• 叶片前缘和尾缘存在严重温度梯度
• 气膜冷却孔边缘应力集中
• 涂层/基体界面热应力不匹配

寿命预测要点：

1. 准确获取叶片温度场分布（CFD+热传导分析）
2. 考虑材料性能的温度依赖性
3. 计入涂层对热应力的影响
4. 采用应变-寿命法进行寿命预测
5. 考虑蠕变-疲劳交互作用

5.2 核反应堆压力容器热疲劳

核反应堆压力容器在运行过程中承受温度瞬态引起的热疲劳：

瞬态类型：

• 冷启动/热启动
• 紧急停堆
• 负荷跟踪
• 安注系统启动

热疲劳敏感区域：

• 进水管与筒体连接处
• 安全端焊缝
• 堆芯筒体过渡区

分析方法：

1. 建立详细的瞬态温度场分析模型
2. 考虑流体混合引起的温度波动（热分层）
3. 采用弹塑性有限元分析计算应力-应变响应
4. 使用疲劳监测系统进行在线损伤评估

5.3 汽车排气系统热疲劳

汽车排气系统承受发动机循环产生的热疲劳载荷：

载荷特征：

• 温度循环：200°C ~ 800°C
• 循环频率：与驾驶工况相关
• 振动载荷叠加

关键部件：

• 排气歧管
• 涡轮增压器壳体
• 催化转化器载体

设计考虑：

1. 材料选择（铸铁、不锈钢、高温合金）
2. 结构优化（柔性设计、膨胀节）
3. 表面强化处理
4. 疲劳寿命验证试验

---

6. 热疲劳寿命预测的不确定性分析

6.1 不确定性来源

热疲劳寿命预测存在多种不确定性来源：

材料性能分散性：

• 疲劳性能参数的批次差异
• 温度相关性能数据的测试误差
• 材料微观结构的不均匀性

载荷不确定性：

• 实际运行温度波动
• 载荷谱的统计变异
• 环境因素的影响

模型不确定性：

• 寿命预测模型的固有误差
• 数值模拟的离散化误差
• 边界条件假设的偏差

6.2 可靠性分析方法

概率寿命预测：

将材料参数和载荷视为随机变量，采用Monte Carlo方法进行概率寿命预测：

$$P_f=P(N_{actual}<N_{design})$$

其中，$P_f$ 为失效概率，$N_{actual}$ 为实际寿命，$N_{design}$ 为设计寿命。

安全系数确定：

考虑不确定性，安全系数应满足：

$$SF=\frac{N_{predicted}}{N_{required}}\ge S{F}_{min}$$

对于核级设备，$S{F}_{min}$ 通常取 2.0 ~ 3.0；对于航空发动机，取 1.5 ~ 2.0。

---

7. 总结与展望

7.1 本文总结

本文系统阐述了热疲劳寿命预测的理论基础和数值方法，主要贡献包括：

1. 理论体系：建立了完整的热疲劳损伤理论框架，涵盖裂纹萌生、扩展和断裂三个阶段，阐明了温度-应力耦合、蠕变-疲劳交互等关键机制。

2. 预测方法：详细介绍了应变-寿命法（Coffin-Manson方程）、能量法（Ostergren模型）、SWT参数法和损伤累积理论（Miner法则），并讨论了各方法的适用范围和局限性。

3. 数值实现：开发了完整的热疲劳寿命预测Python仿真程序，实现了热应力分析、循环塑性响应计算、损伤累积追踪和剩余寿命预测等功能。

4. 工程应用：通过航空发动机涡轮叶片、核反应堆压力容器和汽车排气系统等典型案例，展示了热疲劳寿命预测在工程实践中的应用。

7.2 发展趋势

热疲劳寿命预测技术的发展趋势包括：

多尺度建模：从宏观连续介质向微观晶体塑性、介观晶粒尺度延伸，建立跨尺度的热疲劳损伤模型。

机器学习应用：利用深度学习方法建立数据驱动的寿命预测模型，提高预测精度和计算效率。

数字孪生技术：构建热疲劳数字孪生系统，实现结构热疲劳状态的实时监测和寿命预测。

先进试验技术：发展原位测试技术，在微观尺度上观察热疲劳损伤的演化过程，为模型验证提供数据支撑。

不确定性量化：建立完善的概率寿命预测方法，考虑材料分散性、载荷变异性和模型不确定性，为可靠性设计提供科学依据。

热疲劳寿命预测是保障高温结构安全可靠运行的关键技术，随着计算能力的提升和试验技术的进步，预测精度将不断提高，为重大装备的延寿和安全评估提供更加可靠的技术支撑。