主题044：核反应堆压力容器热应力仿真

目录

1. 引言
2. 核反应堆压力容器概述
3. 压力容器应力分析理论
4. 辐照效应与材料老化
5. 实例一：压力容器热应力分析
6. 代码深度解析
7. 运行结果分析
8. 核安全设计准则
9. 进阶挑战
10. 总结与习题

---

引言

核反应堆压力容器（Reactor Pressure Vessel, RPV）是核电站最关键的核安全设备之一，它构成了反应堆冷却剂压力边界的重要组成部分。作为不可更换的核级设备，压力容器的设计寿命通常要求达到40-60年，在整个服役期间必须承受高温、高压、强辐照等极端工况的综合作用。

压力容器的完整性直接关系到核电站的安全运行。历史上，压力容器失效是导致核事故的重要因素之一。因此，深入理解压力容器在各种工况下的热应力分布、掌握辐照对材料性能的影响规律、建立完善的寿命评估方法，是核工程领域的重要研究课题。

本教程将系统介绍核反应堆压力容器的热应力仿真技术，包括：

• 厚壁圆筒的弹性应力分析理论
• 热应力与压力应力的耦合计算
• 中子辐照对材料性能的影响
• 瞬态事故工况分析
• 核安全设计准则与寿命评估

---

核反应堆压力容器概述

2.1 压力容器的功能与结构

主要功能：

• 包容反应堆堆芯和冷却剂
• 维持冷却剂的高压环境（压水堆约15.5 MPa）
• 构成放射性物质的第一道屏障
• 为堆内构件提供支撑

典型结构（以压水堆为例）：

• 筒体：圆柱形，内径约4-5米，壁厚200-250mm
• 顶盖：可拆卸，便于换料
• 接管：用于冷却剂进出口、安全阀等
• 法兰与密封：采用双道O型环密封
• 堆芯吊篮支撑：位于筒体下部

材料选择：

• 主体材料：低合金钢（如SA-508 Gr.3）
• 内表面堆焊：奥氏体不锈钢（约5-8mm厚）
• 锻件要求：严格的化学成分和力学性能控制

2.2 设计工况分类

根据ASME BPVC规范，压力容器的设计工况分为：

正常运行工况（Condition I）：

• 稳态运行、正常瞬态
• 设计压力：17.1 MPa
• 设计温度：343°C
• 允许应力：$S_m$（设计应力强度）

异常运行工况（Condition II）：

• 预期运行瞬态
• 频率：每年可能多次
• 允许应力：$1.1S_m$

紧急工况（Condition III）：

• 稀有事故
• 频率：每10年可能一次
• 允许应力：$1.5S_m$

极限事故工况（Condition IV）：

• 极限事故
• 频率：整个寿期可能一次
• 允许应力：$2.0S_m$

2.3 压力容器的特殊挑战

不可更换性：

• 压力容器是核岛中唯一不可更换的大型设备
• 设计寿命决定核电站的服役期限
• 延寿评估是核电站寿命管理的核心

辐照脆化：

• 中子辐照导致材料韧脆转变温度升高
• 可能引发脆性断裂风险
• 需要定期监测和评估

热老化：

• 长期高温运行导致材料性能退化
• 焊缝区尤其敏感
• 需要建立热老化模型

---

压力容器应力分析理论

3.1 Lame公式

对于承受内压的厚壁圆筒，应力分布由Lame公式给出：

径向应力：
$${\sigma }_r=\frac{P_i{a}^2-P_o{b}^2}{b^2-a^2}-\frac{(P_i-P_o)a^2{b}^2}{r^2(b^2-a^2)}$$

环向应力：
$${\sigma }_{\theta }=\frac{P_i{a}^2-P_o{b}^2}{b^2-a^2}+\frac{(P_i-P_o)a^2{b}^2}{r^2(b^2-a^2)}$$

轴向应力（闭口容器）：
$${\sigma }_z=\frac{P_i{a}^2}{b^2-a^2}$$

其中：

• $a$：内半径
• $b$：外半径
• $P_i$：内压
• $P_o$：外压
• $r$：计算点半径

应力分布特点：

• 环向应力在内壁面最大
• 径向应力在内壁面等于$-P_i$，在外壁面等于$-P_o$
• 轴向应力沿壁厚均匀分布

3.2 热应力理论

对于沿壁厚非线性温度分布的圆筒，热应力为：

环向热应力：
$${\sigma }_{\theta }^{thermal}=\frac{E\alpha }{1-\nu }\left[\frac{1}{b^2-a^2}\left(1+\frac{b^2}{r^2}\right){\int }_a^{b}Trdr+\frac{1}{r^2}{\int }_a^{r}Trdr-T\right]$$

径向热应力：
$${\sigma }_r^{thermal}=\frac{E\alpha }{1-\nu }\left[\frac{1}{b^2-a^2}\left(1-\frac{b^2}{r^2}\right){\int }_a^{b}Trdr-\frac{1}{r^2}{\int }_a^{r}Trdr\right]$$

轴向热应力：
$${\sigma }_z^{thermal}=\frac{E\alpha }{1-\nu }\left[\frac{2}{b^2-a^2}{\int }_a^{b}Trdr-T\right]$$

对于线性温度分布 $\Delta T=T_o-T_i$：
$${\sigma }_{\theta }^{thermal}=\frac{E\alpha \Delta T}{2(1-\nu )}$$

3.3 组合应力分析

总应力为机械应力与热应力的叠加：
$${\sigma }_{total}={\sigma }_{mechanical}+{\sigma }_{thermal}$$

等效应力（von Mises）：
$${\sigma }_{vm}=\sqrt{\frac{1}{2}\left[({\sigma }_{\theta }-{\sigma }_r{)}^2+({\sigma }_r-{\sigma }_z{)}^2+({\sigma }_z-{\sigma }_{\theta }{)}^2\right]}$$

应力强度（Tresca）：
$${\sigma }_{tresca}=\max (|{\sigma }_{\theta }-{\sigma }_r|,|{\sigma }_r-{\sigma }_z|,|{\sigma }_z-{\sigma }_{\theta }|)$$

---

辐照效应与材料老化

4.1 中子辐照损伤机制

位移损伤：

• 快中子与晶格原子碰撞产生空位和间隙原子
• 形成位错环和空洞
• 导致材料硬化和脆化

嬗变反应：

• 中子俘获产生氢、氦等气体
• 在晶界聚集形成气泡
• 降低材料韧性

活化产物：

• 材料元素被活化产生放射性同位素
• 增加退役处理难度
• 影响维修和检查

4.2 辐照脆化模型

韧脆转变温度（DBTT）上移：
$$\Delta T_{41J}=A\cdot (\text{fluence}{)}^n$$

其中：

• $\Delta T_{41J}$：夏比V型缺口冲击试验41J能量对应的温度变化
• $A$：材料常数
• $n$：指数（通常0.3-0.5）

典型数值：

• 初始DBTT：-20°C至0°C
• 寿期末DBTT：可能升至50-100°C
• 需要定期监测和评估

4.3 辐照监督计划

根据ASME规范，核电站必须建立辐照监督计划：

监督试样：

• 在压力容器内壁安装监督试样盒
• 包含母材、焊缝、热影响区试样
• 定期取出进行力学性能测试

测试项目：

• 拉伸试验（屈服强度、抗拉强度）
• 冲击试验（DBTT、上平台能量）
• 断裂韧性试验（$K_{IC}$、$J_{IC}$）

评估方法：

• 对比初始性能与辐照后性能
• 建立辐照脆化趋势曲线
• 预测寿期末的材料状态

---

实例一：压力容器热应力分析

5.1 问题描述

本实例模拟核反应堆压力容器在正常运行和瞬态工况下的热应力分布，考虑以下因素：

几何参数（典型PWR压力容器）：

• 内半径：2.0 m
• 壁厚：0.2 m（200mm）
• 高度：12.0 m

材料参数（SA-508低合金钢）：

• 弹性模量：200 GPa（室温）
• 热膨胀系数：12×10⁻⁶ /K
• 泊松比：0.3
• 屈服强度：350 MPa（室温）

运行参数：

• 设计压力：17.1 MPa
• 运行压力：15.5 MPa
• 冷却剂入口温度：565 K
• 冷却剂出口温度：600 K
• 环境温度：300 K

5.2 环境准备

运行本实例需要以下Python库：

pip install numpy matplotlib scipy

5.3 完整代码实现

"""
主题044：核反应堆压力容器
实例一：压力容器热应力分析
"""

import numpy as np
import matplotlib.pyplot as plt

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


class ReactorPressureVessel:
    """核反应堆压力容器分析器"""
    
    def __init__(self, inner_radius, wall_thickness, height):
        """
        初始化压力容器参数
        
        参数:
        inner_radius: 内半径 (m)
        wall_thickness: 壁厚 (m)
        height: 高度 (m)
        """
        self.R_i = inner_radius
        self.t = wall_thickness
        self.R_o = inner_radius + wall_thickness
        self.H = height
        
        # 材料参数 - 低合金钢 (SA-508)
        self.E_0 = 200e9  # 室温弹性模量 (Pa)
        self.alpha = 12e-6  # 热膨胀系数 (1/K)
        self.nu = 0.3  # 泊松比
        self.sigma_y_0 = 350e6  # 室温屈服强度 (Pa)
        
        # 辐照效应参数
        self.embrittlement_rate = 10e6  # 脆化速率 (Pa per 10^19 n/cm²)
    
    def material_properties(self, T, fluence=0):
        """
        计算考虑温度效应和辐照效应的材料性能
        """
        # 温度对弹性模量的影响
        E = self.E_0 * (1 - 0.0003 * (T - 300))
        
        # 温度对屈服强度的影响
        sigma_y = self.sigma_y_0 * (1 - 0.0005 * (T - 300))
        
        # 辐照脆化效应
        sigma_y += fluence * self.embrittlement_rate
        
        return E, sigma_y, sigma_y * 1.5
    
    def pressure_stress_cylinder(self, r, P_internal, P_external=0):
        """
        计算内压引起的应力（Lame公式）
        """
        a = self.R_i
        b = self.R_o
        
        # Lame公式
        sigma_r = (P_internal * a**2 - P_external * b**2) / (b**2 - a**2) - \
                  (P_internal - P_external) * a**2 * b**2 / (r**2 * (b**2 - a**2))
        
        sigma_theta = (P_internal * a**2 - P_external * b**2) / (b**2 - a**2) + \
                     (P_internal - P_external) * a**2 * b**2 / (r**2 * (b**2 - a**2))
        
        sigma_z = P_internal * a**2 / (b**2 - a**2)
        
        return sigma_r, sigma_theta, sigma_z
    
    def combined_stress(self, r, z, P_internal, T_inlet, T_outlet, T_ambient, fluence=0):
        """
        计算组合应力
        """
        # 计算温度（简化模型）
        T_coolant = T_inlet + (T_outlet - T_inlet) * z / self.H
        
        if r <= self.R_i:
            T = T_coolant
        else:
            # 简化的对数温度分布
            T_inner = T_coolant + 10
            T_outer = T_ambient + 5
            log_ratio = np.log(r / self.R_i) / np.log(self.R_o / self.R_i)
            T = T_inner + (T_outer - T_inner) * log_ratio
        
        # 材料性能
        E, sigma_y, sigma_u = self.material_properties(T, fluence)
        
        # 压力应力
        sigma_r_p, sigma_theta_p, sigma_z_p = self.pressure_stress_cylinder(r, P_internal)
        
        # 热应力（简化）
        delta_T = T_outer - T_inner
        sigma_theta_t = E * self.alpha * delta_T / (2 * (1 - self.nu))
        sigma_r_t = 0
        sigma_z_t = E * self.alpha * delta_T / (2 * (1 - self.nu))
        
        # 组合应力
        sigma_r = sigma_r_p + sigma_r_t
        sigma_theta = sigma_theta_p + sigma_theta_t
        sigma_z = sigma_z_p + sigma_z_t
        
        # von Mises等效应力
        sigma_vm = np.sqrt(0.5 * ((sigma_theta - sigma_r)**2 + 
                                  (sigma_r - sigma_z)**2 + 
                                  (sigma_z - sigma_theta)**2))
        
        return {
            'sigma_r': sigma_r,
            'sigma_theta': sigma_theta,
            'sigma_z': sigma_z,
            'sigma_vm': sigma_vm,
            'T': T,
            'sigma_y': sigma_y
        }

---

代码深度解析

6.1 Lame公式实现

def pressure_stress_cylinder(self, r, P_internal, P_external=0):
    a = self.R_i
    b = self.R_o
    
    # Lame公式
    sigma_r = (P_internal * a**2 - P_external * b**2) / (b**2 - a**2) - \
              (P_internal - P_external) * a**2 * b**2 / (r**2 * (b**2 - a**2))
    
    sigma_theta = (P_internal * a**2 - P_external * b**2) / (b**2 - a**2) + \
                 (P_internal - P_external) * a**2 * b**2 / (r**2 * (b**2 - a**2))
    
    sigma_z = P_internal * a**2 / (b**2 - a**2)
    
    return sigma_r, sigma_theta, sigma_z

关键理解：

• 径向应力始终为压应力（负值）
• 环向应力在内壁面最大
• 轴向应力与半径无关

6.2 材料性能的温度依赖性

def material_properties(self, T, fluence=0):
    # 温度对弹性模量的影响
    E = self.E_0 * (1 - 0.0003 * (T - 300))
    
    # 温度对屈服强度的影响
    sigma_y = self.sigma_y_0 * (1 - 0.0005 * (T - 300))
    
    # 辐照脆化效应
    sigma_y += fluence * self.embrittlement_rate
    
    return E, sigma_y, sigma_y * 1.5

物理意义：

• 高温下材料软化（强度和模量下降）
• 辐照导致材料硬化和脆化
• 需要综合考虑两种效应

6.3 von Mises等效应力计算

sigma_vm = np.sqrt(0.5 * ((sigma_theta - sigma_r)**2 + 
                          (sigma_r - sigma_z)**2 + 
                          (sigma_z - sigma_theta)**2))

应用准则：

• 当${\sigma }_{vm}<{\sigma }_y$时，材料处于弹性状态
• 当${\sigma }_{vm}={\sigma }_y$时，开始屈服
• 安全系数：$SF={\sigma }_y/{\sigma }_{vm}$

---

运行结果分析

7.1 稳态工况分析结果

内壁面中截面 (r=2.0m, z=6.0m):
  环向应力: -244.8 MPa
  径向应力: -15.5 MPa
  轴向应力: -334.1 MPa
  von Mises等效应力: 284.7 MPa
  温度: 582.5 K
  材料屈服强度: 300.6 MPa
  安全系数: 1.06

结果解读：

1. 应力状态：三向压应力状态，有利于防止脆性断裂
2. 等效应力：284.7 MPa，接近屈服强度
3. 安全系数：1.06，满足设计准则（>1.0）
4. 温度效应：高温导致材料强度下降约14%

7.2 辐照效应分析结果

中子注量 (10¹⁹ n/cm²)	屈服强度 (MPa)	安全系数
0	300.6	1.06
1	310.6	1.09
3	330.6	1.16
5	350.6	1.23
7	370.6	1.30
10	400.6	1.41

关键发现：

• 辐照使材料屈服强度提高
• 安全系数随辐照增加而提高
• 但韧性下降，脆性断裂风险增加
• 需要关注韧脆转变温度的变化

7.3 LOCA事故瞬态分析

事故序列中的安全系数变化：

时间 (s)	压力 (MPa)	温度 (K)	安全系数
0	15.5	600	0.98
10	12.0	580	1.00
30	5.0	520	1.16
100	2.0	450	1.73
200	1.5	420	2.29

安全评估：

• 事故初期（0-10s）安全系数<1.0，存在塑性变形风险
• 应急冷却系统启动后，安全系数快速恢复
• 最终安全系数>2.0，结构安全

---

核安全设计准则

8.1 ASME BPVC规范要求

应力限值：

• 设计工况：${\sigma }_{vm}\le S_m$
• 异常工况：${\sigma }_{vm}\le 1.1S_m$
• 紧急工况：${\sigma }_{vm}\le 1.5S_m$
• 极限事故：${\sigma }_{vm}\le 2.0S_m$

其中$S_m$为设计应力强度，取以下较小值：

• $S_m=\min ({\sigma }_u/3,{\sigma }_y/1.5)$（室温）
• $S_m=\min ({\sigma }_u/3,{\sigma }_y/1.1)$（高温）

8.2 防脆断设计

断裂力学方法：

• 假设存在最大可检测缺陷
• 计算应力强度因子$K_I$
• 要求$K_I<K_{IC}/\sqrt{10}$（安全系数$\sqrt{10}\approx 3.16$）

参考温度法：

• 确定材料的参考零韧性温度$R{T}_{NDT}$
• 要求运行温度高于$R{T}_{NDT}+33°C$
• 定期监测$R{T}_{NDT}$的变化

8.3 疲劳设计

疲劳寿命评估：

• 累积疲劳损伤因子$CUF<1.0$
• 考虑瞬态工况的循环次数
• 使用设计疲劳曲线

环境疲劳修正：

• 高温水环境加速疲劳损伤
• 采用环境疲劳修正因子$F_{en}$
• 修正后的循环次数：$N_{corr}=N_{air}/F_{en}$

---

进阶挑战

9.1 思考题

1. 热应力优化：如何设计压力容器的壁厚分布，使得热应力最小？

2. 多层容器：如果采用多层包扎结构，应力分布会有何变化？

3. 接管应力：压力容器接管的应力集中如何分析？

4. 延寿评估：如何评估已运行30年的压力容器能否延寿20年？

9.2 扩展方向

弹塑性分析：

• 考虑材料的塑性变形
• 采用增量理论或形变理论
• 分析残余应力分布

断裂力学分析：

• 假设裂纹缺陷
• 计算应力强度因子
• 评估裂纹稳定性

概率安全分析：

• 考虑参数的不确定性
• 蒙特卡洛模拟
• 失效概率计算

---

"""
主题044：核反应堆压力容器
实例一：压力容器热应力分析

本实例模拟核反应堆压力容器在正常运行和瞬态工况下的热应力分布，
考虑辐照效应、温度梯度和内压载荷的综合作用。
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle, FancyBboxPatch
import matplotlib.patches as patches

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


class ReactorPressureVessel:
    """核反应堆压力容器分析器"""
    
    def __init__(self, inner_radius, wall_thickness, height):
        """
        初始化压力容器参数
        
        参数:
        inner_radius: 内半径 (m)
        wall_thickness: 壁厚 (m)
        height: 高度 (m)
        """
        self.R_i = inner_radius
        self.t = wall_thickness
        self.R_o = inner_radius + wall_thickness
        self.H = height
        
        # 材料参数 - 低合金钢 (SA-508)
        self.E_0 = 200e9  # 室温弹性模量 (Pa)
        self.alpha = 12e-6  # 热膨胀系数 (1/K)
        self.nu = 0.3  # 泊松比
        self.sigma_y_0 = 350e6  # 室温屈服强度 (Pa)
        
        # 辐照效应参数
        self.fluence = 0.0  # 中子注量 (n/cm²)
        self.embrittlement_rate = 10e6  # 脆化速率 (Pa per 10^19 n/cm²)
        
    def material_properties(self, T, fluence=0):
        """
        计算考虑温度效应和辐照效应的材料性能
        
        参数:
        T: 温度 (K)
        fluence: 中子注量 (10^19 n/cm²)
        
        返回:
        E: 弹性模量 (Pa)
        sigma_y: 屈服强度 (Pa)
        sigma_u: 抗拉强度 (Pa)
        """
        # 温度对弹性模量的影响（线性下降）
        E = self.E_0 * (1 - 0.0003 * (T - 300))
        
        # 温度对屈服强度的影响
        sigma_y = self.sigma_y_0 * (1 - 0.0005 * (T - 300))
        
        # 辐照脆化效应
        sigma_y += fluence * self.embrittlement_rate
        
        # 抗拉强度（约为屈服强度的1.5倍）
        sigma_u = sigma_y * 1.5
        
        return E, sigma_y, sigma_u
    
    def temperature_distribution(self, r, z, T_inlet, T_outlet, T_ambient):
        """
        计算稳态温度分布
        
        假设：
        - 轴向线性温度分布（冷却剂加热）
        - 径向对数温度分布（圆柱壁导热）
        
        参数:
        r: 径向坐标 (m)
        z: 轴向坐标 (m)
        T_inlet: 冷却剂入口温度 (K)
        T_outlet: 冷却剂出口温度 (K)
        T_ambient: 环境温度 (K)
        
        返回:
        T: 温度 (K)
        """
        # 轴向温度分布（线性）
        T_coolant = T_inlet + (T_outlet - T_inlet) * z / self.H
        
        # 径向温度分布（对数）
        # 简化的对数分布
        if r <= self.R_i:
            # 内部冷却剂区域
            T = T_coolant
        else:
            # 壁面区域
            k = 40.0  # 热导率 W/(m·K)
            h_in = 5000.0  # 内侧对流系数 W/(m²·K)
            h_out = 10.0  # 外侧对流系数 W/(m²·K)
            
            # 简化的温度分布
            T_inner = T_coolant + 10  # 内壁温度略高于冷却剂
            T_outer = T_ambient + 5   # 外壁温度略高于环境
            
            # 对数插值
            log_ratio = np.log(r / self.R_i) / np.log(self.R_o / self.R_i)
            T = T_inner + (T_outer - T_inner) * log_ratio
        
        return T
    
    def thermal_stress_cylinder(self, r, T_profile, E, alpha, nu):
        """
        计算圆柱壁的热应力
        
        基于厚壁圆筒热应力理论
        
        参数:
        r: 径向坐标 (m)
        T_profile: 温度分布函数
        E: 弹性模量 (Pa)
        alpha: 热膨胀系数 (1/K)
        nu: 泊松比
        
        返回:
        sigma_r: 径向应力 (Pa)
        sigma_theta: 环向应力 (Pa)
        sigma_z: 轴向应力 (Pa)
        """
        # 简化的热应力计算
        # 假设温度沿壁厚线性分布
        T_inner = T_profile(self.R_i, 0)
        T_outer = T_profile(self.R_o, 0)
        delta_T = T_outer - T_inner
        
        # 热应力公式（简化）
        # 环向热应力（最大）
        sigma_theta_thermal = E * alpha * delta_T / (2 * (1 - nu))
        
        # 径向热应力（内壁处为0）
        sigma_r_thermal = 0
        
        # 轴向热应力
        sigma_z_thermal = E * alpha * delta_T / (2 * (1 - nu))
        
        return sigma_r_thermal, sigma_theta_thermal, sigma_z_thermal
    
    def pressure_stress_cylinder(self, r, P_internal, P_external=0):
        """
        计算内压引起的应力（Lame公式）
        
        参数:
        r: 径向坐标 (m)
        P_internal: 内压 (Pa)
        P_external: 外压 (Pa)
        
        返回:
        sigma_r: 径向应力 (Pa)
        sigma_theta: 环向应力 (Pa)
        sigma_z: 轴向应力 (Pa)
        """
        a = self.R_i
        b = self.R_o
        
        # Lame公式
        sigma_r = (P_internal * a**2 - P_external * b**2) / (b**2 - a**2) - \
                  (P_internal - P_external) * a**2 * b**2 / (r**2 * (b**2 - a**2))
        
        sigma_theta = (P_internal * a**2 - P_external * b**2) / (b**2 - a**2) + \
                     (P_internal - P_external) * a**2 * b**2 / (r**2 * (b**2 - a**2))
        
        # 轴向应力（闭口容器）
        sigma_z = P_internal * a**2 / (b**2 - a**2)
        
        return sigma_r, sigma_theta, sigma_z
    
    def combined_stress(self, r, z, P_internal, T_inlet, T_outlet, T_ambient, fluence=0):
        """
        计算组合应力
        
        参数:
        r: 径向坐标 (m)
        z: 轴向坐标 (m)
        P_internal: 内压 (Pa)
        T_inlet, T_outlet: 冷却剂进出口温度 (K)
        T_ambient: 环境温度 (K)
        fluence: 中子注量
        
        返回:
        stress_dict: 应力分量字典
        """
        # 计算温度
        T = self.temperature_distribution(r, z, T_inlet, T_outlet, T_ambient)
        
        # 材料性能
        E, sigma_y, sigma_u = self.material_properties(T, fluence)
        
        # 压力应力
        sigma_r_p, sigma_theta_p, sigma_z_p = self.pressure_stress_cylinder(r, P_internal)
        
        # 热应力（简化计算）
        sigma_r_t, sigma_theta_t, sigma_z_t = self.thermal_stress_cylinder(
            r, lambda r, z: self.temperature_distribution(r, z, T_inlet, T_outlet, T_ambient),
            E, self.alpha, self.nu
        )
        
        # 组合应力
        sigma_r = sigma_r_p + sigma_r_t
        sigma_theta = sigma_theta_p + sigma_theta_t
        sigma_z = sigma_z_p + sigma_z_t
        
        # 等效应力（von Mises）
        sigma_vm = np.sqrt(0.5 * ((sigma_theta - sigma_r)**2 + 
                                  (sigma_r - sigma_z)**2 + 
                                  (sigma_z - sigma_theta)**2))
        
        # 应力强度（Tresca）
        sigma_tresca = max(abs(sigma_theta - sigma_r), 
                          abs(sigma_r - sigma_z), 
                          abs(sigma_z - sigma_theta))
        
        return {
            'sigma_r': sigma_r,
            'sigma_theta': sigma_theta,
            'sigma_z': sigma_z,
            'sigma_vm': sigma_vm,
            'sigma_tresca': sigma_tresca,
            'T': T,
            'E': E,
            'sigma_y': sigma_y,
            'sigma_u': sigma_u
        }


def example_1_steady_state_analysis():
    """示例1：稳态工况分析"""
    
    print("=" * 70)
    print("示例1：压力容器稳态热应力分析")
    print("=" * 70)
    
    # 创建压力容器模型
    # 典型PWR压力容器尺寸
    vessel = ReactorPressureVessel(
        inner_radius=2.0,  # 内半径 2m
        wall_thickness=0.2,  # 壁厚 200mm
        height=12.0  # 高度 12m
    )
    
    # 运行参数
    P_internal = 15.5e6  # 内压 15.5 MPa
    T_inlet = 565.0  # 入口温度 565K
    T_outlet = 600.0  # 出口温度 600K
    T_ambient = 300.0  # 环境温度 300K
    
    print(f"\n运行参数:")
    print(f"  内压: {P_internal/1e6:.1f} MPa")
    print(f"  冷却剂入口温度: {T_inlet:.1f} K")
    print(f"  冷却剂出口温度: {T_outlet:.1f} K")
    
    # 创建分析网格
    nr = 50
    nz = 60
    r = np.linspace(vessel.R_i, vessel.R_o, nr)
    z = np.linspace(0, vessel.H, nz)
    R, Z = np.meshgrid(r, z)
    
    # 计算应力场
    print("\n计算应力场...")
    sigma_theta = np.zeros((nz, nr))
    sigma_vm = np.zeros((nz, nr))
    T_field = np.zeros((nz, nr))
    safety_factor = np.zeros((nz, nr))
    
    for i in range(nz):
        for j in range(nr):
            result = vessel.combined_stress(
                r[j], z[i], P_internal, T_inlet, T_outlet, T_ambient
            )
            sigma_theta[i, j] = result['sigma_theta'] / 1e6  # MPa
            sigma_vm[i, j] = result['sigma_vm'] / 1e6
            T_field[i, j] = result['T']
            safety_factor[i, j] = result['sigma_y'] / result['sigma_vm']
    
    # 关键位置分析
    print("\n关键位置应力分析:")
    
    # 内壁面
    result_inner = vessel.combined_stress(
        vessel.R_i, vessel.H/2, P_internal, T_inlet, T_outlet, T_ambient
    )
    print(f"\n内壁面中截面 (r={vessel.R_i}m, z={vessel.H/2}m):")
    print(f"  环向应力: {result_inner['sigma_theta']/1e6:.1f} MPa")
    print(f"  径向应力: {result_inner['sigma_r']/1e6:.1f} MPa")
    print(f"  轴向应力: {result_inner['sigma_z']/1e6:.1f} MPa")
    print(f"  von Mises等效应力: {result_inner['sigma_vm']/1e6:.1f} MPa")
    print(f"  温度: {result_inner['T']:.1f} K")
    print(f"  材料屈服强度: {result_inner['sigma_y']/1e6:.1f} MPa")
    print(f"  安全系数: {result_inner['sigma_y']/result_inner['sigma_vm']:.2f}")
    
    # 外壁面
    result_outer = vessel.combined_stress(
        vessel.R_o, vessel.H/2, P_internal, T_inlet, T_outlet, T_ambient
    )
    print(f"\n外壁面中截面 (r={vessel.R_o}m, z={vessel.H/2}m):")
    print(f"  环向应力: {result_outer['sigma_theta']/1e6:.1f} MPa")
    print(f"  径向应力: {result_outer['sigma_r']/1e6:.1f} MPa")
    print(f"  轴向应力: {result_outer['sigma_z']/1e6:.1f} MPa")
    print(f"  von Mises等效应力: {result_outer['sigma_vm']/1e6:.1f} MPa")
    print(f"  温度: {result_outer['T']:.1f} K")
    print(f"  材料屈服强度: {result_outer['sigma_y']/1e6:.1f} MPa")
    print(f"  安全系数: {result_outer['sigma_y']/result_outer['sigma_vm']:.2f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 1. 温度分布
    ax1 = axes[0, 0]
    contour1 = ax1.contourf(R, Z, T_field, levels=20, cmap='hot')
    ax1.set_xlabel('半径 (m)', fontsize=11)
    ax1.set_ylabel('高度 (m)', fontsize=11)
    ax1.set_title('温度分布 (K)', fontsize=12, fontweight='bold')
    plt.colorbar(contour1, ax=ax1)
    
    # 2. 环向应力
    ax2 = axes[0, 1]
    contour2 = ax2.contourf(R, Z, sigma_theta, levels=20, cmap='RdYlBu_r')
    ax2.set_xlabel('半径 (m)', fontsize=11)
    ax2.set_ylabel('高度 (m)', fontsize=11)
    ax2.set_title('环向应力 (MPa)', fontsize=12, fontweight='bold')
    plt.colorbar(contour2, ax=ax2)
    
    # 3. von Mises等效应力
    ax3 = axes[0, 2]
    contour3 = ax3.contourf(R, Z, sigma_vm, levels=20, cmap='jet')
    ax3.set_xlabel('半径 (m)', fontsize=11)
    ax3.set_ylabel('高度 (m)', fontsize=11)
    ax3.set_title('von Mises等效应力 (MPa)', fontsize=12, fontweight='bold')
    plt.colorbar(contour3, ax=ax3)
    
    # 4. 沿壁厚方向的应力分布
    ax4 = axes[1, 0]
    mid_z_idx = nz // 2
    ax4.plot(r, sigma_theta[mid_z_idx, :], 'b-', linewidth=2, label='环向应力')
    ax4.plot(r, sigma_vm[mid_z_idx, :], 'r--', linewidth=2, label='von Mises应力')
    ax4.axhline(result_inner['sigma_y']/1e6, color='g', linestyle=':', 
               label=f'屈服强度={result_inner["sigma_y"]/1e6:.0f}MPa')
    ax4.set_xlabel('半径 (m)', fontsize=11)
    ax4.set_ylabel('应力 (MPa)', fontsize=11)
    ax4.set_title('中截面应力分布', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 5. 安全系数分布
    ax5 = axes[1, 1]
    contour5 = ax5.contourf(R, Z, safety_factor, levels=20, cmap='RdYlGn')
    ax5.set_xlabel('半径 (m)', fontsize=11)
    ax5.set_ylabel('高度 (m)', fontsize=11)
    ax5.set_title('安全系数', fontsize=12, fontweight='bold')
    plt.colorbar(contour5, ax=ax5)
    
    # 6. 应力-强度比
    ax6 = axes[1, 2]
    stress_ratio = sigma_vm / (result_inner['sigma_y']/1e6)
    contour6 = ax6.contourf(R, Z, stress_ratio, levels=20, cmap='jet', vmax=1.0)
    ax6.set_xlabel('半径 (m)', fontsize=11)
    ax6.set_ylabel('高度 (m)', fontsize=11)
    ax6.set_title('应力/屈服强度比', fontsize=12, fontweight='bold')
    plt.colorbar(contour6, ax=ax6)
    
    plt.tight_layout()
    plt.savefig('example1_steady_state.png', dpi=150, bbox_inches='tight')
    print("\n✓ 已保存: example1_steady_state.png")
    plt.close()


def example_2_irradiation_effects():
    """示例2：辐照效应分析"""
    
    print("\n" + "=" * 70)
    print("示例2：中子辐照对材料性能的影响")
    print("=" * 70)
    
    vessel = ReactorPressureVessel(
        inner_radius=2.0,
        wall_thickness=0.2,
        height=12.0
    )
    
    # 不同中子注量水平
    fluence_levels = [0, 1e19, 3e19, 5e19, 7e19, 10e19]  # n/cm²
    
    P_internal = 15.5e6
    T_inlet = 565.0
    T_outlet = 600.0
    T_ambient = 300.0
    
    results = []
    
    for fluence in fluence_levels:
        print(f"\n中子注量: {fluence/1e19:.1f}×10¹⁹ n/cm²")
        
        # 计算内壁面应力
        result = vessel.combined_stress(
            vessel.R_i, vessel.H/2, P_internal, T_inlet, T_outlet, T_ambient, 
            fluence/1e19
        )
        
        safety_factor = result['sigma_y'] / result['sigma_vm']
        
        results.append({
            'fluence': fluence,
            'sigma_y': result['sigma_y'] / 1e6,
            'sigma_vm': result['sigma_vm'] / 1e6,
            'safety_factor': safety_factor
        })
        
        print(f"  屈服强度: {result['sigma_y']/1e6:.1f} MPa")
        print(f"  等效应力: {result['sigma_vm']/1e6:.1f} MPa")
        print(f"  安全系数: {safety_factor:.2f}")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    fluence_array = np.array([r['fluence']/1e19 for r in results])
    sigma_y_array = np.array([r['sigma_y'] for r in results])
    sf_array = np.array([r['safety_factor'] for r in results])
    
    # 屈服强度变化
    ax1 = axes[0]
    ax1.plot(fluence_array, sigma_y_array, 'bo-', linewidth=2, markersize=8)
    ax1.set_xlabel('中子注量 (10¹⁹ n/cm²)', fontsize=11)
    ax1.set_ylabel('屈服强度 (MPa)', fontsize=11)
    ax1.set_title('辐照对屈服强度的影响', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 安全系数变化
    ax2 = axes[1]
    ax2.plot(fluence_array, sf_array, 'rs-', linewidth=2, markersize=8)
    ax2.axhline(1.5, color='orange', linestyle='--', label='设计限值=1.5')
    ax2.set_xlabel('中子注量 (10¹⁹ n/cm²)', fontsize=11)
    ax2.set_ylabel('安全系数', fontsize=11)
    ax2.set_title('辐照对安全系数的影响', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('irradiation_effects.png', dpi=150, bbox_inches='tight')
    print("\n✓ 已保存: irradiation_effects.png")
    plt.close()


def example_3_transient_analysis():
    """示例3：冷却剂丧失事故(LOCA)瞬态分析"""
    
    print("\n" + "=" * 70)
    print("示例3：冷却剂丧失事故(LOCA)瞬态分析")
    print("=" * 70)
    
    vessel = ReactorPressureVessel(
        inner_radius=2.0,
        wall_thickness=0.2,
        height=12.0
    )
    
    # 事故序列
    # t=0: 正常运行
    # t=10s: 冷却剂丧失，压力下降
    # t=50s: 应急冷却系统启动
    
    time_points = np.array([0, 5, 10, 20, 30, 50, 100, 200])  # 秒
    P_internal = np.array([15.5, 15.5, 12.0, 8.0, 5.0, 3.0, 2.0, 1.5]) * 1e6  # Pa
    T_coolant = np.array([600, 600, 580, 550, 520, 480, 450, 420])  # K
    
    results = []
    
    for i, t in enumerate(time_points):
        print(f"\nt = {t}s:")
        print(f"  压力: {P_internal[i]/1e6:.1f} MPa")
        print(f"  冷却剂温度: {T_coolant[i]:.1f} K")
        
        # 计算应力
        result = vessel.combined_stress(
            vessel.R_i, vessel.H/2, P_internal[i], 
            T_coolant[i]-20, T_coolant[i], 300.0
        )
        
        safety_factor = result['sigma_y'] / result['sigma_vm']
        
        results.append({
            'time': t,
            'pressure': P_internal[i] / 1e6,
            'temperature': T_coolant[i],
            'sigma_vm': result['sigma_vm'] / 1e6,
            'safety_factor': safety_factor
        })
        
        print(f"  等效应力: {result['sigma_vm']/1e6:.1f} MPa")
        print(f"  安全系数: {safety_factor:.2f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    time_array = np.array([r['time'] for r in results])
    pressure_array = np.array([r['pressure'] for r in results])
    temp_array = np.array([r['temperature'] for r in results])
    stress_array = np.array([r['sigma_vm'] for r in results])
    sf_array = np.array([r['safety_factor'] for r in results])
    
    # 压力变化
    ax1 = axes[0, 0]
    ax1.plot(time_array, pressure_array, 'b-', linewidth=2, marker='o')
    ax1.set_xlabel('时间 (s)', fontsize=11)
    ax1.set_ylabel('压力 (MPa)', fontsize=11)
    ax1.set_title('LOCA事故压力变化', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 温度变化
    ax2 = axes[0, 1]
    ax2.plot(time_array, temp_array, 'r-', linewidth=2, marker='s')
    ax2.set_xlabel('时间 (s)', fontsize=11)
    ax2.set_ylabel('温度 (K)', fontsize=11)
    ax2.set_title('LOCA事故温度变化', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 应力变化
    ax3 = axes[1, 0]
    ax3.plot(time_array, stress_array, 'g-', linewidth=2, marker='^')
    ax3.set_xlabel('时间 (s)', fontsize=11)
    ax3.set_ylabel('等效应力 (MPa)', fontsize=11)
    ax3.set_title('LOCA事故应力变化', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 安全系数变化
    ax4 = axes[1, 1]
    ax4.plot(time_array, sf_array, 'm-', linewidth=2, marker='d')
    ax4.axhline(1.0, color='red', linestyle='--', linewidth=2, label='失效限值=1.0')
    ax4.set_xlabel('时间 (s)', fontsize=11)
    ax4.set_ylabel('安全系数', fontsize=11)
    ax4.set_title('LOCA事故安全系数变化', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('loca_transient.png', dpi=150, bbox_inches='tight')
    print("\n✓ 已保存: loca_transient.png")
    plt.close()


def example_4_design_optimization():
    """示例4：壁厚优化设计"""
    
    print("\n" + "=" * 70)
    print("示例4：压力容器壁厚优化设计")
    print("=" * 70)
    
    # 设计参数范围
    thicknesses = np.linspace(0.15, 0.35, 20)  # 壁厚 150-350mm
    
    P_design = 17.0e6  # 设计压力 17 MPa
    T_operating = 600.0  # 运行温度 600K
    
    results = []
    
    for t in thicknesses:
        vessel = ReactorPressureVessel(
            inner_radius=2.0,
            wall_thickness=t,
            height=12.0
        )
        
        # 计算内壁面应力
        result = vessel.combined_stress(
            vessel.R_i, vessel.H/2, P_design, 
            T_operating-20, T_operating, 300.0
        )
        
        safety_factor = result['sigma_y'] / result['sigma_vm']
        weight = np.pi * (vessel.R_o**2 - vessel.R_i**2) * vessel.H * 7800  # 钢材密度7800
        
        results.append({
            'thickness': t * 1000,  # mm
            'sigma_vm': result['sigma_vm'] / 1e6,
            'safety_factor': safety_factor,
            'weight': weight / 1000  # kg
        })
    
    # 可视化
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    t_array = np.array([r['thickness'] for r in results])
    stress_array = np.array([r['sigma_vm'] for r in results])
    sf_array = np.array([r['safety_factor'] for r in results])
    weight_array = np.array([r['weight'] for r in results])
    
    # 应力-壁厚关系
    ax1 = axes[0]
    ax1.plot(t_array, stress_array, 'b-', linewidth=2)
    ax1.set_xlabel('壁厚 (mm)', fontsize=11)
    ax1.set_ylabel('等效应力 (MPa)', fontsize=11)
    ax1.set_title('应力-壁厚关系', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 安全系数-壁厚关系
    ax2 = axes[1]
    ax2.plot(t_array, sf_array, 'r-', linewidth=2)
    ax2.axhline(1.5, color='orange', linestyle='--', label='设计限值=1.5')
    ax2.axhline(2.0, color='green', linestyle='--', label='保守设计=2.0')
    ax2.set_xlabel('壁厚 (mm)', fontsize=11)
    ax2.set_ylabel('安全系数', fontsize=11)
    ax2.set_title('安全系数-壁厚关系', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 重量-壁厚关系
    ax3 = axes[2]
    ax3.plot(t_array, weight_array/1000, 'g-', linewidth=2)
    ax3.set_xlabel('壁厚 (mm)', fontsize=11)
    ax3.set_ylabel('重量 (吨)', fontsize=11)
    ax3.set_title('重量-壁厚关系', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('design_optimization.png', dpi=150, bbox_inches='tight')
    print("\n✓ 已保存: design_optimization.png")
    plt.close()
    
    # 推荐设计
    print("\n设计建议:")
    for r in results:
        if 1.8 <= r['safety_factor'] <= 2.2:
            print(f"  壁厚 {r['thickness']:.0f}mm: 安全系数={r['safety_factor']:.2f}, "
                  f"重量={r['weight']/1000:.1f}吨")


if __name__ == "__main__":
    print("=" * 70)
    print("主题044：核反应堆压力容器")
    print("实例一：压力容器热应力分析")
    print("=" * 70)
    
    # 示例1：稳态工况分析
    example_1_steady_state_analysis()
    
    # 示例2：辐照效应分析
    example_2_irradiation_effects()
    
    # 示例3：LOCA瞬态分析
    example_3_transient_analysis()
    
    # 示例4：设计优化
    example_4_design_optimization()
    
    print("\n" + "=" * 70)
    print("实例一完成！")
    print("=" * 70)
    print("\n关键概念:")
    print("1. 压力容器的Lame应力公式")
    print("2. 热应力与压力应力的叠加")
    print("3. 中子辐照脆化效应")
    print("4. 核安全设计准则")
    print("5. 瞬态事故分析")