第027篇：热-流-固耦合分析

摘要

热-流-固耦合分析是工程仿真中最具挑战性的多物理场问题之一，涉及热传导、流体流动和结构变形三个物理场的强耦合作用。本主题系统阐述热-流-固耦合的基本理论、数学模型和数值求解方法。首先介绍热-流-固耦合的物理机制，包括流体热对流、热应力效应和流固界面传热等核心概念；然后建立完整的耦合控制方程组，涵盖能量守恒、动量守恒和结构动力学方程；接着讨论耦合界面的连续性条件和数值处理方法；最后通过六个典型工程案例展示热-流-固耦合分析的实际应用，包括热交换器管板应力分析、航空发动机叶片热疲劳、核反应堆燃料棒热流固耦合、电子器件散热变形、汽车排气系统热应力和深海管道热屈曲等。每个案例均配有详细的Python仿真代码和可视化结果，帮助读者深入理解热-流-固耦合问题的建模与求解技术。

关键词

热-流-固耦合，热弹性力学，流固耦合，热对流，界面传热，耦合算法，有限元方法，热应力，热疲劳，工程应用

---

1. 热-流-固耦合的物理机制

1.1 耦合现象概述

热-流-固耦合（Thermo-Fluid-Structure Interaction, TFSI）描述的是热场、流场和结构场三者之间的相互作用过程。这种耦合在航空航天、能源动力、核工程、电子封装等领域广泛存在，是工程实践中最为复杂的多物理场问题之一。

热-流-固耦合的核心特征在于三个物理场之间存在双向或多向的相互作用：

热场对流场的影响：流体温度分布改变流体的密度、粘度和导热系数等热物性参数，进而影响流体的流动特性。温度梯度还会引起自然对流，这是由浮力驱动的流动现象。

流场对热场的影响：流体流动通过对流换热影响固体表面的温度分布。对流换热系数取决于流体的速度、物性和流动状态（层流或湍流）。高速流动区域的换热通常更为强烈。

热场对结构场的影响：温度变化引起材料的热膨胀或收缩，产生热应变和热应力。不均匀的温度分布会导致结构变形，严重时可能引发热屈曲或热疲劳失效。

结构场对流场的影响：结构变形改变流道的几何形状，影响流体的流动路径和流动特性。大变形情况下，这种影响可能非常显著，如柔性管道的流致振动。

流场对结构场的影响：流体压力作用于结构表面，产生流体动力载荷。非定常流动还会引起结构的振动响应，如涡激振动和颤振现象。

1.2 热对流机理

热对流是热-流-固耦合中的关键传热机制，分为自然对流和强制对流两种类型。

自然对流是由浮力驱动的流动，当流体中存在温度梯度导致密度不均匀时产生。自然对流的强度用格拉晓夫数（Grashof number）表征：

$$Gr=\frac{g\beta \Delta T{L}^3}{{\nu }^2}$$

其中，$g$为重力加速度，$\beta$为热膨胀系数，$\Delta T$为特征温差，$L$为特征长度，$\nu$为运动粘度。

自然对流与强制对流的相对重要性由理查森数（Richardson number）决定：

$$Ri=\frac{Gr}{R{e}^2}$$

当$Ri\gg 1$时，自然对流占主导；当$Ri\ll 1$时，强制对流占主导；当$Ri\approx 1$时，两种对流机制同等重要。

强制对流是由外部驱动力（如泵、风机）引起的流动。强制对流的换热强度用努塞尔数（Nusselt number）表征：

$$Nu=\frac{hL}{k_f}$$

其中，$h$为对流换热系数，$k_f$为流体导热系数。努塞尔数与雷诺数和普朗特数相关：

$$Nu=C\cdot R{e}^m\cdot P{r}^n$$

对于层流管内流动，常用Dittus-Boelter关联式：

$$Nu=0.023\cdot R{e}^{0.8}\cdot P{r}^{0.4}$$

1.3 热弹性力学基础

热弹性力学研究温度变化引起的结构变形和应力问题。热应变与温度变化的关系由热膨胀定律描述：

$${\epsilon }^{th}=\alpha \Delta T$$

其中，$\alpha$为热膨胀系数，$\Delta T=T-T_{ref}$为相对于参考温度的温差。

对于各向同性材料，三维热应变张量为：

$${\mathbf{\epsilon }}^{th}=\alpha \Delta T\mathbf{I}$$

总应变由机械应变和热应变组成：

$$\mathbf{\epsilon }={\mathbf{\epsilon }}^{mech}+{\mathbf{\epsilon }}^{th}$$

热应力本构关系为：

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

展开形式为：

$${\sigma }_{ij}=2G{\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}-\frac{E\alpha \Delta T}{1-2\nu }{\delta }_{ij}$$

其中，$G$为剪切模量，$\lambda$为拉梅常数，$E$为弹性模量，$\nu$为泊松比。

热应力平衡方程为：

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

展开后包含热应力项：

$$G{\nabla }^2\mathbf{u}+(\lambda +G)\nabla (\nabla \cdot \mathbf{u})-\frac{E\alpha }{1-2\nu }\nabla T+\mathbf{f}=\rho \frac{{\partial }^2\mathbf{u}}{\partial t^2}$$

1.4 界面传热与耦合条件

热-流-固耦合问题的关键在于正确处理界面处的物理量传递。在流固界面上，需要满足以下连续性条件：

温度连续性：

$$T_f=T_s\text{on}{\Gamma }_{fs}$$

热流连续性：

$$-k_f\frac{\partial T_f}{\partial n}=-k_s\frac{\partial T_s}{\partial n}\text{on}{\Gamma }_{fs}$$

位移连续性（无滑移条件）：

$${\mathbf{u}}_f={\mathbf{u}}_s\text{on}{\Gamma }_{fs}$$

力平衡：

$${\mathbf{\sigma }}_f\cdot \mathbf{n}={\mathbf{\sigma }}_s\cdot \mathbf{n}\text{on}{\Gamma }_{fs}$$

对于对流换热边界，常用牛顿冷却定律：

$$q=h(T_s-T_f)$$

其中，$h$为对流换热系数，取决于流动状态和流体物性。

---

2. 热-流-固耦合的数学模型

2.1 控制方程组

热-流-固耦合问题由以下控制方程组描述：

流体域控制方程：

连续性方程（质量守恒）：

$$\frac{\partial {\rho }_f}{\partial t}+\nabla \cdot ({\rho }_f{\mathbf{u}}_f)=0$$

对于不可压缩流动：

$$\nabla \cdot {\mathbf{u}}_f=0$$

动量方程（Navier-Stokes方程）：

$${\rho }_f\left(\frac{\partial {\mathbf{u}}_f}{\partial t}+{\mathbf{u}}_f\cdot \nabla {\mathbf{u}}_f\right)=-\nabla p+\mu {\nabla }^2{\mathbf{u}}_f+{\rho }_f\mathbf{g}$$

能量方程：

$${\rho }_f{c}_p\left(\frac{\partial T_f}{\partial t}+{\mathbf{u}}_f\cdot \nabla T_f\right)=k_f{\nabla }^2{T}_f+\Phi$$

其中，$\Phi$为粘性耗散函数：

$$\Phi =2\mu \mathbf{S}:\mathbf{S}$$

$\mathbf{S}$为应变率张量：

$$S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$

固体域控制方程：

热传导方程：

$${\rho }_s{c}_s\frac{\partial T_s}{\partial t}=k_s{\nabla }^2{T}_s+Q$$

其中，$Q$为内热源。

结构动力学方程：

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

热弹性本构关系：

$${\sigma }_{ij}=C_{ijkl}\left({\epsilon }_{kl}-{\alpha }_{kl}\Delta T\right)$$

2.2 无量纲化分析

为便于分析不同尺度的问题，引入无量纲变量：

$${\mathbf{u}}^{\ast }=\frac{\mathbf{u}}{U},p^{\ast }=\frac{p}{\rho U^2},T^{\ast }=\frac{T-T_0}{\Delta T}$$

$$t^{\ast }=\frac{tU}{L},{\mathbf{x}}^{\ast }=\frac{\mathbf{x}}{L}$$

无量纲控制方程揭示了几个关键的无量纲数：

雷诺数（Reynolds number）：惯性力与粘性力之比

$$Re=\frac{\rho UL}{\mu }$$

普朗特数（Prandtl number）：动量扩散与热扩散之比

$$Pr=\frac{\mu c_p}{k}=\frac{\nu }{\alpha }$$

毕奥数（Biot number）：内部热阻与外部热阻之比

$$Bi=\frac{hL}{k_s}$$

傅里叶数（Fourier number）：热传导与热储存之比

$$Fo=\frac{\alpha t}{L^2}$$

埃克特数（Eckert number）：动能与焓变之比

$$Ec=\frac{U^2}{c_p\Delta T}$$

弗劳德数（Froude number）：惯性力与重力之比

$$Fr=\frac{U}{\sqrt{gL}}$$

2.3 耦合强度分析

热-流-固耦合的强度可以通过耦合系数来量化：

热-流耦合系数：

$$C_{tf}=\frac{\text{对流换热}}{\text{导热}}=\frac{hL}{k_f}=Nu$$

流-固耦合系数：

$$C_{fs}=\frac{\text{流体动压}}{\text{结构刚度}}=\frac{\rho U^2}{E}$$

热-固耦合系数：

$$C_{ts}=\frac{\text{热应力}}{\text{机械应力}}=\frac{E\alpha \Delta T}{{\sigma }_{mech}}$$

根据耦合系数的相对大小，可以将热-流-固耦合问题分为：

• 弱耦合：$C\ll 1$，各物理场可以顺序求解
• 中等耦合：$C\approx 1$，需要迭代求解
• 强耦合：$C\gg 1$，需要全耦合求解

2.4 简化模型与假设

根据具体问题的特点，可以对热-流-固耦合模型进行合理简化：

稳态假设：当瞬态效应不重要时，可以忽略时间导数项：

$$\frac{\partial }{\partial t}=0$$

不可压缩假设：对于低速流动（$Ma<0.3$），可以假设流体不可压缩：

$$\rho =\text{const}$$

** Boussinesq近似**：对于自然对流问题，仅在浮力项中考虑密度变化：

$$\rho ={\rho }_0[1-\beta (T-T_0)]$$

线性热弹性假设：对于小变形和小温差，采用线性热弹性理论：

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

准静态假设：当惯性效应不重要时，忽略结构动力学方程中的加速度项：

$${\rho }_s\frac{{\partial }^2{\mathbf{u}}_s}{\partial t^2}\approx 0$$

---

3. 数值求解方法

3.1 分区耦合法

分区耦合法（Partitioned Coupling）将热-流-固问题分解为三个子问题分别求解，通过界面信息传递实现耦合。这种方法的优点是可以利用现有的单物理场求解器，灵活性高。

求解流程：

1. 求解流体域，获得流场和温度场分布
2. 将流体域的壁面热流传递给固体域作为边界条件
3. 求解固体域热传导，获得固体温度场
4. 将固体壁面温度传递给流体域
5. 计算固体热应力和变形
6. 更新流体域网格（如有必要）
7. 检查收敛性，如未收敛则返回步骤1

迭代格式：

对于第$n$次迭代：

$$T_s^{n+1}=\mathcal{S}(q_f^n)$$

$$T_f^{n+1}=\mathcal{F}(T_s^{n+1})$$

$$q_f^{n+1}=-k_f\frac{\partial T_f^{n+1}}{\partial n}$$

其中，$\mathcal{S}$和$\mathcal{F}$分别表示固体和流体求解算子。

松弛技术：为提高收敛性，引入松弛因子：

$$T_s^{n+1}=\omega \cdot T_s^{new}+(1-\omega )\cdot T_s^n$$

通常取$\omega =0.5\sim 0.8$。

3.2 整体耦合法

整体耦合法（Monolithic Coupling）将所有控制方程联立求解，形成一个大型代数方程组。这种方法收敛性好，但计算成本高。

离散化后的方程组：

$$\left[\begin{array}{ccc}{\mathbf{K}}_{ff}& {\mathbf{K}}_{fs}& 0\\ {\mathbf{K}}_{sf}& {\mathbf{K}}_{ss}& {\mathbf{K}}_{st}\\ 0& {\mathbf{K}}_{ts}& {\mathbf{K}}_{tt}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_f\\ {\mathbf{u}}_s\\ \mathbf{T}\end{array}\right]=\left[\begin{array}{c}{\mathbf{f}}_f\\ {\mathbf{f}}_s\\ {\mathbf{f}}_t\end{array}\right]$$

其中，下标$f$、$s$、$t$分别表示流体、结构和温度自由度。

求解策略：

• 直接求解：适用于中小规模问题
• 迭代求解：GMRES、BiCGSTAB等Krylov子空间方法
• 多重网格：加速收敛
• 预处理器：块Jacobi、块Gauss-Seidel等

3.3 界面插值与数据传递

在分区耦合中，流体和固体网格通常不匹配，需要进行界面数据插值。

常用插值方法：

1. 最近邻插值：将最近的节点值赋予目标点

   $${\varphi }_{target}={\varphi }_{nearest}$$

2. 线性插值：基于三角形或四面体单元的形函数插值

   $${\varphi }_{target}=\sum _{i=1}^n{N}_i{\varphi }_i$$

3. 径向基函数（RBF）插值：

   $${\varphi }_{target}=\sum _{j=1}^m{\alpha }_j\phi (||\mathbf{x}-{\mathbf{x}}_j||)$$

   其中，$\phi (r)=e^{-(\epsilon r{)}^2}$为高斯RBF。

守恒性要求：

界面数据传递需要满足守恒性：

$${\int }_{\Gamma }{q}_{f}d\Gamma ={\int }_{\Gamma }{q}_{s}d\Gamma$$

3.4 动网格技术

当结构变形较大时，需要更新流体域网格。

弹簧近似法：将网格边视为弹簧，节点位移通过求解弹簧系统获得：

$$\sum _{j\in N_i}{k}_{ij}({\mathbf{u}}_i-{\mathbf{u}}_j)=0$$

其中，$k_{ij}=1/||{\mathbf{x}}_i-{\mathbf{x}}_j||$为弹簧刚度。

弹性体法：将网格视为弹性体，求解线弹性方程：

$$\nabla \cdot {\mathbf{\sigma }}_{mesh}=0$$

径向基函数法：使用RBF插值网格位移：

$${\mathbf{u}}_{mesh}(\mathbf{x})=\sum _{j=1}^n{\mathbf{\alpha }}_j\phi (||\mathbf{x}-{\mathbf{x}}_j||)+\mathbf{p}(\mathbf{x})$$

其中，$\mathbf{p}(\mathbf{x})$为多项式项。

---

4. 工程应用案例

4.1 案例1：热交换器管板应力分析

问题描述：管壳式热交换器是化工和能源行业常用的换热设备。高温流体在管内流动，通过管壁将热量传递给管外流体。管板作为支撑管束的关键部件，承受着热应力和流体压力的联合作用。

几何参数：

• 管板直径：$D=0.6$ m
• 管板厚度：$t=0.05$ m
• 管孔数量：$N=100$
• 管孔直径：$d=0.025$ m
• 管间距：$P=0.04$ m（三角形排列）

材料参数：

• 管板材料：碳钢
• 弹性模量：$E=200$ GPa
• 泊松比：$\nu =0.3$
• 热膨胀系数：$\alpha =12\times 10^{-6}$ /K
• 导热系数：$k=45$ W/(m·K)

工况条件：

• 管程流体温度：$T_{tube}=300$°C
• 壳程流体温度：$T_{shell}=150$°C
• 管程压力：$p_{tube}=5$ MPa
• 壳程压力：$p_{shell}=2$ MPa

仿真目标：

1. 计算管板温度分布
2. 分析热应力和机械应力的叠加效应
3. 确定最大应力位置和安全裕度

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例1：热交换器管板应力分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题027'
os.makedirs(output_dir, exist_ok=True)

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

print("\n" + "="*60)
print("案例5：汽车排气系统热应力分析")
print("="*60)

# 排气歧管几何参数
L_manifold = 0.4  # 歧管长度 (m)
D_pipe = 0.045  # 管径 (m)
t_wall = 0.003  # 壁厚 (m)
N_branches = 4  # 分支数
t_flange = 0.015  # 法兰厚度 (m)

# 材料参数（不锈钢）
E_steel = 200e9  # 弹性模量 (Pa)
nu_steel = 0.3  # 泊松比
alpha_steel = 16e-6  # 热膨胀系数 (/K)
k_steel = 25  # 导热系数 (W/(m·K))
rho_steel = 7850  # 密度 (kg/m³)

# 工况条件
T_exhaust_peak = 850  # 废气峰值温度 (°C)
T_exhaust_steady = 600  # 废气稳态温度 (°C)
p_exhaust = 0.15e6  # 废气压力 (Pa)
T_amb = 30  # 环境温度 (°C)
h_inner = 200  # 内侧换热系数 (W/(m²·K))
h_outer = 30  # 外侧换热系数 (W/(m²·K))

# 创建歧管网格
n_x = 50
n_theta = 36
x = np.linspace(0, L_manifold, n_x)
theta = np.linspace(0, 2*np.pi, n_theta)
X, Theta = np.meshgrid(x, theta)

# 计算温度场（考虑热瞬态）
# 简化模型：峰值温度工况
time_steps = 100
t_max = 300  # 最大时间 (s)
time = np.linspace(0, t_max, time_steps)

# 径向温度分布（圆柱坐标）
r_inner = D_pipe / 2
r_outer = r_inner + t_wall
n_r = 20
r = np.linspace(r_inner, r_outer, n_r)

# 计算稳态温度分布
T_inner = T_exhaust_peak - (T_exhaust_peak - T_amb) / (1 + h_inner/h_outer + h_inner*t_wall/k_steel)
T_outer = T_inner - h_inner * (T_exhaust_peak - T_inner) * t_wall / k_steel

# 温度分布
T_wall = np.zeros((n_theta, n_x))
for i in range(n_theta):
    for j in range(n_x):
        # 沿长度方向的温度衰减
        T_local = T_inner * np.exp(-x[j] / 0.5) + T_amb
        T_wall[i, j] = T_local

# 计算热应力
# 热应力 = E * alpha * delta_T / (1 - nu)
delta_T_wall = T_wall - T_amb
sigma_th_manifold = E_steel * alpha_steel * delta_T_wall / (1 - nu_steel)

# 内压引起的应力（薄壁圆筒）
sigma_p_hoop = p_exhaust * r_inner / t_wall  # 环向应力
sigma_p_long = p_exhaust * r_inner / (2 * t_wall)  # 纵向应力

# 组合应力（von Mises）
sigma_vm_manifold = np.sqrt(sigma_th_manifold**2 + sigma_p_hoop**2 - sigma_th_manifold*sigma_p_hoop)

print(f"歧管内壁温度: {T_inner:.1f}°C")
print(f"歧管外壁温度: {T_outer:.1f}°C")
print(f"最大热应力: {np.max(sigma_th_manifold)/1e6:.2f} MPa")
print(f"内压环向应力: {sigma_p_hoop/1e6:.2f} MPa")
print(f"最大von Mises应力: {np.max(sigma_vm_manifold)/1e6:.2f} MPa")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 排气歧管几何示意图
ax1 = axes[0, 0]
# 绘制主管
y_offset = 0
for i in range(N_branches):
    angle = i * 2 * np.pi / N_branches
    x_branch = np.linspace(0, L_manifold*1000, 20)
    y_branch = np.sin(angle) * 20
    ax1.plot(x_branch, [y_branch]*20, 'b-', linewidth=8, alpha=0.7)
    # 分支入口
    ax1.plot([0, 0], [0, y_branch], 'b-', linewidth=6, alpha=0.7)

# 法兰
flange = Rectangle((-10, -30), 20, 60, fill=True, color='gray', alpha=0.5)
ax1.add_patch(flange)

ax1.set_xlim(-20, L_manifold*1000 + 20)
ax1.set_ylim(-40, 40)
ax1.set_xlabel('Length (mm)', fontsize=11)
ax1.set_ylabel('Radial Position (mm)', fontsize=11)
ax1.set_title('Exhaust Manifold Geometry', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.set_aspect('equal')

# 2. 温度分布
ax2 = axes[0, 1]
contour2 = ax2.contourf(X*1000, Theta*180/np.pi, T_wall, levels=20, cmap='hot')
ax2.set_xlabel('Length (mm)', fontsize=11)
ax2.set_ylabel('Circumferential Angle (deg)', fontsize=11)
ax2.set_title('Wall Temperature Distribution (°C)', fontsize=12, fontweight='bold')
cbar2 = plt.colorbar(contour2, ax=ax2)
cbar2.set_label('Temperature (°C)', fontsize=10)

# 3. 热应力分布
ax3 = axes[1, 0]
contour3 = ax3.contourf(X*1000, Theta*180/np.pi, sigma_th_manifold/1e6, levels=20, cmap='jet')
ax3.set_xlabel('Length (mm)', fontsize=11)
ax3.set_ylabel('Circumferential Angle (deg)', fontsize=11)
ax3.set_title('Thermal Stress (MPa)', fontsize=12, fontweight='bold')
cbar3 = plt.colorbar(contour3, ax=ax3)
cbar3.set_label('Stress (MPa)', fontsize=10)

# 4. von Mises等效应力
ax4 = axes[1, 1]
contour4 = ax4.contourf(X*1000, Theta*180/np.pi, sigma_vm_manifold/1e6, levels=20, cmap='jet')
ax4.set_xlabel('Length (mm)', fontsize=11)
ax4.set_ylabel('Circumferential Angle (deg)', fontsize=11)
ax4.set_title('von Mises Stress (MPa)', fontsize=12, fontweight='bold')
cbar4 = plt.colorbar(contour4, ax=ax4)
cbar4.set_label('Stress (MPa)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case5_exhaust_system.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例5完成，图片已保存")

# 热疲劳寿命估算
delta_T_cycle = T_exhaust_peak - T_exhaust_steady
C_fatigue = 5e12  # 材料常数
m_fatigue = 3  # 指数
N_fatigue = C_fatigue / (delta_T_cycle**m_fatigue)
print(f"\n热疲劳分析:")
print(f"温度循环幅值: {delta_T_cycle:.1f}°C")
print(f"估算热疲劳寿命: {N_fatigue:.0f} 循环")

4.6 案例6：深海管道热屈曲分析

问题描述：深海油气输送管道在海底低温环境和输送高温油气的联合作用下，可能产生热屈曲失稳。当管道轴向热应力超过临界屈曲载荷时，管道会发生侧向屈曲或上浮屈曲，导致结构失效。

管道几何：

• 管道外径：$D_o=0.3239$ m（12英寸）
• 管道壁厚：$t=0.0127$ m
• 管道长度：$L=1000$ m
• 埋深：$H=1.5$ m（部分埋设）

材料参数（X65钢）：

• 弹性模量：$E=207$ GPa
• 屈服强度：${\sigma }_y=450$ MPa
• 热膨胀系数：$\alpha =12\times 10^{-6}$ /K
• 导热系数：$k=45$ W/(m·K)
• 泊松比：$\nu =0.3$

工况条件：

• 输送介质温度：$T_{in}=120$°C
• 海水温度：$T_{sw}=4$°C
• 输送压力：$p=15$ MPa
• 土壤摩擦系数：$\mu =0.6$
• 土壤覆盖压力：$q=50$ kN/m

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例6：深海管道热屈曲分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题027'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "="*60)
print("案例6：深海管道热屈曲分析")
print("="*60)

# 管道几何参数
D_o = 0.3239  # 管道外径 (m)
t_pipe = 0.0127  # 壁厚 (m)
L_pipe = 1000  # 管道长度 (m)
H_buried = 1.5  # 埋深 (m)

D_i = D_o - 2 * t_pipe  # 内径
r_o = D_o / 2
r_i = D_i / 2

# 材料参数（X65钢）
E_pipe = 207e9  # 弹性模量 (Pa)
sigma_y = 450e6  # 屈服强度 (Pa)
alpha_pipe = 12e-6  # 热膨胀系数 (/K)
k_pipe = 45  # 导热系数 (W/(m·K))
nu_pipe = 0.3  # 泊松比

# 工况条件
T_inlet = 120  # 输送介质温度 (°C)
T_seawater = 4  # 海水温度 (°C)
p_oper = 15e6  # 操作压力 (Pa)
mu_soil = 0.6  # 土壤摩擦系数
q_soil = 50e3  # 土壤覆盖压力 (N/m)

# 计算温度场
# 简化模型：沿管道长度方向温度衰减
x_pipe = np.linspace(0, L_pipe, 200)
h_conv = 100  # 总换热系数 (W/(m²·K))
U_overall = 1 / (1/h_conv + t_pipe/k_pipe)  # 总传热系数

# 温度沿管道衰减
T_pipe = T_seawater + (T_inlet - T_seawater) * np.exp(-x_pipe / 5000)

# 计算热应力
# 约束热应力 = E * alpha * delta_T
delta_T_pipe = T_pipe - T_seawater
sigma_thermal = E_pipe * alpha_pipe * delta_T_pipe

# 截面特性
A_cross = np.pi * (r_o**2 - r_i**2)  # 截面积
I_moment = np.pi * (r_o**4 - r_i**4) / 4  # 惯性矩

# 轴向刚度
EA = E_pipe * A_cross

# 屈曲分析
# 侧向屈曲临界载荷（基于弹性地基梁理论）
k_soil = 100e6  # 地基刚度 (N/m²)
P_cr_lateral = 2 * np.sqrt(k_soil * E_pipe * I_moment)

# 上浮屈曲临界载荷
P_cr_upheaval = q_soil * L_pipe / 10  # 简化模型

# 实际轴向力（热应力 + 内压引起的泊松效应）
nu_effect = nu_pipe * p_oper * np.pi * r_i**2 / A_cross  # 泊松效应
P_actual = sigma_thermal * A_cross + nu_effect * A_cross

# 安全系数
SF_lateral = P_cr_lateral / np.max(P_actual)
SF_upheaval = P_cr_upheaval / np.max(P_actual)

print(f"管道入口温度: {T_inlet:.1f}°C")
print(f"管道出口温度: {T_pipe[-1]:.1f}°C")
print(f"最大热应力: {np.max(sigma_thermal)/1e6:.2f} MPa")
print(f"侧向屈曲临界载荷: {P_cr_lateral/1e6:.2f} MN")
print(f"上浮屈曲临界载荷: {P_cr_upheaval/1e6:.2f} MN")
print(f"最大实际轴向力: {np.max(P_actual)/1e6:.2f} MN")
print(f"侧向屈曲安全系数: {SF_lateral:.2f}")
print(f"上浮屈曲安全系数: {SF_upheaval:.2f}")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 管道系统示意图
ax1 = axes[0, 0]
# 海底
seabed = Rectangle((-50, -2), L_pipe + 100, 2, fill=True, color='sandybrown', alpha=0.7)
ax1.add_patch(seabed)
# 管道
pipe = Rectangle((0, 0), L_pipe, D_o*10, fill=True, color='steelblue', alpha=0.8)
ax1.add_patch(pipe)
# 海水
water = Rectangle((-50, D_o*10), L_pipe + 100, 50, fill=True, color='lightblue', alpha=0.3)
ax1.add_patch(water)

ax1.set_xlim(-50, L_pipe + 50)
ax1.set_ylim(-5, 60)
ax1.set_xlabel('Pipeline Length (m)', fontsize=11)
ax1.set_ylabel('Elevation (m)', fontsize=11)
ax1.set_title('Subsea Pipeline System', fontsize=12, fontweight='bold')
ax1.text(L_pipe/2, 30, 'Seawater', ha='center', fontsize=10)
ax1.text(L_pipe/2, -1, 'Seabed', ha='center', fontsize=10)

# 2. 沿管道温度分布
ax2 = axes[0, 1]
ax2.plot(x_pipe, T_pipe, 'r-', linewidth=2, label='Pipe Temperature')
ax2.axhline(y=T_seawater, color='b', linestyle='--', label='Seawater Temperature')
ax2.fill_between(x_pipe, T_seawater, T_pipe, alpha=0.3, color='red')
ax2.set_xlabel('Pipeline Length (m)', fontsize=11)
ax2.set_ylabel('Temperature (°C)', fontsize=11)
ax2.set_title('Temperature Distribution Along Pipeline', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 3. 热应力分布
ax3 = axes[1, 0]
ax3.plot(x_pipe, sigma_thermal/1e6, 'b-', linewidth=2)
ax3.axhline(y=sigma_y/1e6, color='r', linestyle='--', label='Yield Strength')
ax3.fill_between(x_pipe, 0, sigma_thermal/1e6, alpha=0.3, color='blue')
ax3.set_xlabel('Pipeline Length (m)', fontsize=11)
ax3.set_ylabel('Thermal Stress (MPa)', fontsize=11)
ax3.set_title('Thermal Stress Distribution', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 4. 轴向力与屈曲临界值对比
ax4 = axes[1, 1]
ax4.plot(x_pipe, P_actual/1e6, 'g-', linewidth=2, label='Actual Axial Force')
ax4.axhline(y=P_cr_lateral/1e6, color='r', linestyle='--', label='Lateral Buckling Limit')
ax4.axhline(y=P_cr_upheaval/1e6, color='m', linestyle='--', label='Upheaval Buckling Limit')
ax4.fill_between(x_pipe, 0, P_actual/1e6, alpha=0.3, color='green')
ax4.set_xlabel('Pipeline Length (m)', fontsize=11)
ax4.set_ylabel('Axial Force (MN)', fontsize=11)
ax4.set_title('Buckling Analysis', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case6_subsea_pipeline.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例6完成，图片已保存")

# 屈曲风险评估
print(f"\n屈曲风险评估:")
if SF_lateral > 2.0 and SF_upheaval > 2.0:
    print("✓ 管道设计安全，无屈曲风险")
elif SF_lateral > 1.5 or SF_upheaval > 1.5:
    print("⚠ 管道存在潜在屈曲风险，建议增加约束")
else:
    print("✗ 管道屈曲风险高，需要重新设计")

print("\n" + "="*60)
print("所有案例仿真完成！")
print("="*60)

---

5. 总结与展望

5.1 本章小结

本主题系统阐述了热-流-固耦合分析的理论基础、数学模型和工程应用。主要内容包括：

1. 物理机制：深入分析了热场、流场和结构场之间的双向耦合关系，包括热对流机理、热弹性力学基础和界面传热条件。

2. 数学模型：建立了完整的热-流-固耦合控制方程组，包括流体域的Navier-Stokes方程和能量方程，以及固体域的热传导方程和结构动力学方程。通过无量纲化分析揭示了雷诺数、普朗特数、毕奥数等关键无量纲数的物理意义。

3. 数值方法：介绍了分区耦合法和整体耦合法两种求解策略，讨论了界面插值、数据传递和动网格技术等关键数值处理方法。

4. 工程应用：通过六个典型案例展示了热-流-固耦合分析在热交换器、航空发动机、核反应堆、电子器件、汽车排气和深海管道等领域的实际应用。

5.2 技术挑战与发展趋势

热-流-固耦合分析仍面临诸多挑战：

计算效率：全耦合求解计算成本高昂，需要发展更高效的耦合算法和预处理技术。机器学习方法在加速耦合求解方面展现出潜力。

不确定性量化：材料参数、边界条件和几何尺寸的 uncertainties 对耦合分析结果影响显著，需要发展鲁棒的 uncertainty quantification 方法。

多尺度问题：从微观材料行为到宏观结构响应，跨尺度耦合建模仍是研究热点。

实验验证：高温、高压环境下的耦合场测量技术有待发展，以提供更可靠的验证数据。

工业应用：将先进的耦合分析方法集成到工程设计流程中，发展面向特定行业的专用仿真工具和标准规范。

5.3 学习建议

对于希望深入学习热-流-固耦合分析的读者，建议：

1. 夯实基础：掌握流体力学、传热学和固体力学的基本理论，理解各物理场的控制方程和本构关系。

2. 熟悉数值方法：学习有限差分法、有限元法和有限体积法等数值离散技术，理解耦合界面的处理方法。

3. 实践练习：通过修改本主题提供的Python代码，探索不同参数对耦合行为的影响，加深对物理机制的理解。

4. 阅读文献：关注国际期刊如Journal of Computational Physics、International Journal of Heat and Mass Transfer、Computer Methods in Applied Mechanics and Engineering等发表的最新研究成果。

5. 使用商业软件：学习ANSYS、COMSOL、ABAQUS等商业多物理场仿真软件，了解工业界的实际应用需求。

热-流-固耦合分析是连接基础理论与工程实践的桥梁，掌握这一技术对于解决复杂的工程问题具有重要意义。随着计算能力的提升和数值方法的进步，热-流-固耦合分析将在更多领域发挥关键作用。
axes.unicode_minus’] = False

print(“=”*60)
print(“案例1：热交换器管板应力分析”)
print(“=”*60)

几何参数

D = 0.6 # 管板直径 (m)
t = 0.05 # 管板厚度 (m)
N_tubes = 100 # 管孔数量
d_tube = 0.025 # 管孔直径 (m)
P = 0.04 # 管间距 (m)

材料参数

E = 200e9 # 弹性模量 (Pa)
nu = 0.3 # 泊松比
alpha = 12e-6 # 热膨胀系数 (/K)
k = 45 # 导热系数 (W/(m·K))

工况条件

T_tube = 300 # 管程温度 (°C)
T_shell = 150 # 壳程温度 (°C)
p_tube = 5e6 # 管程压力 (Pa)
p_shell = 2e6 # 壳程压力 (Pa)
T_ref = 20 # 参考温度 (°C)

创建管板网格

n_r = 50
n_theta = 60
r = np.linspace(0, D/2, n_r)
theta = np.linspace(0, 2*np.pi, n_theta)
R, Theta = np.meshgrid(r, theta)
X = R * np.cos(Theta)
Y = R * np.sin(Theta)

生成管孔位置（三角形排列）

n_rings = int((D/2) / P) + 1
tube_positions = []
for i in range(n_rings):
r_ring = i * P
if r_ring > D/2 - d_tube:
break
n_tubes_ring = max(1, int(2 * np.pi * r_ring / P))
for j in range(n_tubes_ring):
angle = 2 * np.pi * j / n_tubes_ring
if i % 2 == 1:
angle += np.pi / n_tubes_ring
x_tube = r_ring * np.cos(angle)
y_tube = r_ring * np.sin(angle)
if np.sqrt(x_tube2 + y_tube2) < D/2 - d_tube:
tube_positions.append((x_tube, y_tube))

tube_positions = tube_positions[:N_tubes]
print(f"实际布置管孔数量: {len(tube_positions)}")

计算温度场（简化模型）

假设管孔区域温度接近管程温度，其他区域温度渐变

T_field = np.zeros_like(X)
h_eff = 1000 # 有效换热系数 (W/(m²·K))

for i in range(X.shape[0]):
for j in range(X.shape[1]):
x, y = X[i,j], Y[i,j]
# 计算到最近管孔的距离
min_dist = float(‘inf’)
for tx, ty in tube_positions:
dist = np.sqrt((x-tx)**2 + (y-ty)**2)
if dist < min_dist:
min_dist = dist

    # 温度插值
    if min_dist < d_tube/2:
        T_field[i,j] = T_tube
    else:
        # 指数衰减的温度分布
        T_field[i,j] = T_shell + (T_tube - T_shell) * np.exp(-min_dist / 0.05)

计算热应力（简化模型）

热应力与温度梯度成正比

sigma_th = E * alpha * (T_field - T_ref) / (1 - nu)

压力引起的应力（简化模型）

delta_p = p_tube - p_shell
sigma_p = delta_p * D / (4 * t) # 薄膜应力

组合应力（考虑应力集中系数）

K_t = 3.0 # 管孔应力集中系数
sigma_total = K_t * (sigma_th + sigma_p)

限制最大应力

sigma_max = 300e6 # 材料许用应力 (Pa)
sigma_total = np.clip(sigma_total, 0, sigma_max * 1.5)

print(f"管板最大热应力: {np.max(sigma_th)/1e6:.2f} MPa")
print(f"压力引起的应力: {sigma_p/1e6:.2f} MPa")
print(f"组合最大应力: {np.max(sigma_total)/1e6:.2f} MPa")

绘制结果

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

1. 管板几何和管孔分布

ax1 = axes[0, 0]

绘制管板外圆

circle = plt.Circle((0, 0), D/2, fill=False, color=‘black’, linewidth=2)
ax1.add_patch(circle)

绘制管孔

for tx, ty in tube_positions:
hole = plt.Circle((tx, ty), d_tube/2, fill=True, color=‘lightblue’, alpha=0.7)
ax1.add_patch(hole)
ax1.set_xlim(-D/21.1, D/21.1)
ax1.set_ylim(-D/21.1, D/21.1)
ax1.set_aspect(‘equal’)
ax1.set_xlabel(‘x (m)’, fontsize=11)
ax1.set_ylabel(‘y (m)’, fontsize=11)
ax1.set_title(‘Tube Sheet Geometry’, fontsize=12, fontweight=‘bold’)
ax1.grid(True, alpha=0.3)

2. 温度分布

ax2 = axes[0, 1]
contour2 = ax2.contourf(X, Y, T_field, levels=20, cmap=‘hot’)
ax2.set_aspect(‘equal’)
ax2.set_xlabel(‘x (m)’, fontsize=11)
ax2.set_ylabel(‘y (m)’, fontsize=11)
ax2.set_title(‘Temperature Distribution (°C)’, fontsize=12, fontweight=‘bold’)
cbar2 = plt.colorbar(contour2, ax=ax2)
cbar2.set_label(‘Temperature (°C)’, fontsize=10)

3. 热应力分布

ax3 = axes[1, 0]
contour3 = ax3.contourf(X, Y, sigma_th/1e6, levels=20, cmap=‘jet’)
ax3.set_aspect(‘equal’)
ax3.set_xlabel(‘x (m)’, fontsize=11)
ax3.set_ylabel(‘y (m)’, fontsize=11)
ax3.set_title(‘Thermal Stress (MPa)’, fontsize=12, fontweight=‘bold’)
cbar3 = plt.colorbar(contour3, ax=ax3)
cbar3.set_label(‘Stress (MPa)’, fontsize=10)

4. 组合应力分布

ax4 = axes[1, 1]
contour4 = ax4.contourf(X, Y, sigma_total/1e6, levels=20, cmap=‘jet’)
ax4.set_aspect(‘equal’)
ax4.set_xlabel(‘x (m)’, fontsize=11)
ax4.set_ylabel(‘y (m)’, fontsize=11)
ax4.set_title(‘Total Stress (MPa)’, fontsize=12, fontweight=‘bold’)
cbar4 = plt.colorbar(contour4, ax=ax4)
cbar4.set_label(‘Stress (MPa)’, fontsize=10)

plt.tight_layout()
plt.savefig(f’{output_dir}/case1_heat_exchanger.png’, dpi=150, bbox_inches=‘tight’)
plt.close()
print(“✓ 案例1完成，图片已保存”)

安全评估

safety_factor = sigma_max / np.max(sigma_total)
print(f"\n安全评估:“)
print(f"材料许用应力: {sigma_max/1e6:.2f} MPa”)
print(f"实际最大应力: {np.max(sigma_total)/1e6:.2f} MPa")
print(f"安全系数: {safety_factor:.2f}")
if safety_factor > 1.5:
print(“✓ 管板设计安全”)
else:
print(“⚠ 管板应力超标，需要优化设计”)

### 4.2 案例2：航空发动机叶片热疲劳分析

**问题描述**：航空发动机涡轮叶片在高温燃气环境中工作，承受着复杂的热应力和机械应力循环。叶片前缘温度可达1400°C以上，而冷却通道内的冷却空气温度仅约600°C，巨大的温差导致严重的热应力。同时，发动机转速变化引起离心力和气动力的周期性变化，导致热疲劳失效。

**叶片几何**：
- 叶片长度：$L = 0.12$ m
- 叶根弦长：$C_{root} = 0.04$ m
- 叶尖弦长：$C_{tip} = 0.025$ m
- 最大厚度：$t_{max} = 0.006$ m
- 冷却通道数：$N_c = 5$

**材料参数**（镍基高温合金）：
- 弹性模量：$E = 120$ GPa（随温度变化）
- 泊松比：$\nu = 0.3$
- 热膨胀系数：$\alpha = 14 \times 10^{-6}$ /K
- 导热系数：$k = 15$ W/(m·K)
- 密度：$\rho = 8200$ kg/m³

**工况条件**：
- 燃气温度：$T_g = 1400$°C
- 冷却空气温度：$T_c = 600$°C
- 转速：$\omega = 15000$ rpm
- 换热系数（燃气侧）：$h_g = 2000$ W/(m²·K)
- 换热系数（冷却侧）：$h_c = 800$ W/(m²·K)

**仿真目标**：
1. 计算稳态温度场分布
2. 分析热应力和离心应力的叠加
3. 评估热疲劳寿命

**Python仿真代码**：

```python
# -*- coding: utf-8 -*-
"""
案例2：航空发动机叶片热疲劳分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon, FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题027'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "="*60)
print("案例2：航空发动机叶片热疲劳分析")
print("="*60)

# 叶片几何参数
L_blade = 0.12  # 叶片长度 (m)
C_root = 0.04  # 叶根弦长 (m)
C_tip = 0.025  # 叶尖弦长 (m)
t_max = 0.006  # 最大厚度 (m)
N_c = 5  # 冷却通道数

# 材料参数
E_0 = 120e9  # 室温弹性模量 (Pa)
nu = 0.3  # 泊松比
alpha = 14e-6  # 热膨胀系数 (/K)
k = 15  # 导热系数 (W/(m·K))
rho = 8200  # 密度 (kg/m³)

# 工况条件
T_g = 1400  # 燃气温度 (°C)
T_c = 600  # 冷却空气温度 (°C)
omega = 15000 * 2 * np.pi / 60  # 角速度 (rad/s)
h_g = 2000  # 燃气侧换热系数 (W/(m²·K))
h_c = 800  # 冷却侧换热系数 (W/(m²·K))

# 创建叶片网格
n_span = 50
n_chord = 60
span = np.linspace(0, L_blade, n_span)
chord = np.linspace(0, 1, n_chord)

# 叶片截面形状参数化（NACA四位数翼型近似）
def airfoil_shape(x_c, t_c, chord_len):
    """生成翼型形状"""
    # NACA 0012近似
    y_t = 5 * t_c * (0.2969*np.sqrt(x_c) - 0.1260*x_c - 0.3516*x_c**2 + 
                     0.2843*x_c**3 - 0.1015*x_c**4)
    return y_t * chord_len

# 计算温度场
T_blade = np.zeros((n_span, n_chord))
for i, s in enumerate(span):
    # 弦长随展向变化
    C_local = C_root + (C_tip - C_root) * s / L_blade
    
    # 当地燃气温度（考虑径向温度梯度）
    T_g_local = T_g - 50 * s / L_blade
    
    for j, x_c in enumerate(chord):
        # 到表面的距离（简化模型）
        y_t = airfoil_shape(x_c, 0.12, 1.0)
        
        # 温度分布：表面接近燃气温度，内部接近冷却温度
        # 使用一维热传导近似
        Bi_g = h_g * t_max / k  # 燃气侧毕奥数
        Bi_c = h_c * t_max / k  # 冷却侧毕奥数
        
        # 简化温度分布
        T_surface = (Bi_g * T_g_local + Bi_c * T_c) / (Bi_g + Bi_c)
        T_blade[i, j] = T_surface - (T_surface - T_c) * y_t * 2

# 计算热应力
# 热应力与温度梯度和约束有关
T_avg = np.mean(T_blade)
delta_T = T_blade - T_avg

# 弹性模量随温度变化
E_T = E_0 * (1 - 0.0003 * (T_blade - 20))  # 简化模型
E_T = np.clip(E_T, E_0 * 0.5, E_0)

# 热应力（简化模型）
sigma_th_blade = E_T * alpha * delta_T / (1 - nu)

# 离心应力（随半径增加）
radius = 0.2 + span  # 从叶根开始的半径 (m)
sigma_centrifugal = rho * omega**2 * radius[:, np.newaxis] * C_root / 2
sigma_centrifugal = np.tile(sigma_centrifugal, (1, n_chord))

# 组合应力
sigma_total_blade = sigma_th_blade + sigma_centrifugal

# 计算热疲劳参数
# 使用Manson-Coffin公式估算疲劳寿命
delta_epsilon = alpha * (np.max(T_blade) - np.min(T_blade))
epsilon_f = 0.5  # 疲劳延性系数
c = -0.6  # 疲劳延性指数
N_f = (epsilon_f / delta_epsilon)**(1/abs(c))

print(f"叶片最高温度: {np.max(T_blade):.1f}°C")
print(f"叶片最低温度: {np.min(T_blade):.1f}°C")
print(f"最大热应力: {np.max(np.abs(sigma_th_blade))/1e6:.2f} MPa")
print(f"最大离心应力: {np.max(sigma_centrifugal)/1e6:.2f} MPa")
print(f"最大组合应力: {np.max(sigma_total_blade)/1e6:.2f} MPa")
print(f"估算热疲劳寿命: {N_f:.0f} 循环")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 叶片几何形状
ax1 = axes[0, 0]
# 绘制叶片截面
s_plot = [0, L_blade/2, L_blade]
for idx, s in enumerate(s_plot):
    C_local = C_root + (C_tip - C_root) * s / L_blade
    x_af = np.linspace(0, 1, 50) * C_local
    y_upper = airfoil_shape(np.linspace(0, 1, 50), 0.12, C_local)
    y_lower = -y_upper
    
    offset = idx * C_root * 1.5
    ax1.plot(x_af + offset, y_upper + s*1000, 'b-', linewidth=2)
    ax1.plot(x_af + offset, y_lower + s*1000, 'b-', linewidth=2)
    ax1.fill_between(x_af + offset, y_lower + s*1000, y_upper + s*1000, alpha=0.3)
    
    # 绘制冷却通道
    for k in range(N_c):
        x_cool = C_local * (0.2 + 0.6 * k / (N_c-1)) + offset
        y_cool = s * 1000
        circle = plt.Circle((x_cool, y_cool), 0.002, fill=True, color='cyan', alpha=0.8)
        ax1.add_patch(circle)

ax1.set_xlabel('Chord Position (mm)', fontsize=11)
ax1.set_ylabel('Span Position (mm)', fontsize=11)
ax1.set_title('Turbine Blade Geometry', fontsize=12, fontweight='bold')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# 2. 温度分布
ax2 = axes[0, 1]
Span, Chord = np.meshgrid(span*1000, chord)
im2 = ax2.contourf(Span, Chord, T_blade.T, levels=20, cmap='hot')
ax2.set_xlabel('Span (mm)', fontsize=11)
ax2.set_ylabel('Normalized Chord', fontsize=11)
ax2.set_title('Temperature Distribution (°C)', fontsize=12, fontweight='bold')
cbar2 = plt.colorbar(im2, ax=ax2)
cbar2.set_label('Temperature (°C)', fontsize=10)

# 3. 热应力分布
ax3 = axes[1, 0]
im3 = ax3.contourf(Span, Chord, sigma_th_blade.T/1e6, levels=20, cmap='RdBu_r')
ax3.set_xlabel('Span (mm)', fontsize=11)
ax3.set_ylabel('Normalized Chord', fontsize=11)
ax3.set_title('Thermal Stress (MPa)', fontsize=12, fontweight='bold')
cbar3 = plt.colorbar(im3, ax=ax3)
cbar3.set_label('Stress (MPa)', fontsize=10)

# 4. 组合应力分布
ax4 = axes[1, 1]
im4 = ax4.contourf(Span, Chord, sigma_total_blade.T/1e6, levels=20, cmap='jet')
ax4.set_xlabel('Span (mm)', fontsize=11)
ax4.set_ylabel('Normalized Chord', fontsize=11)
ax4.set_title('Total Stress (MPa)', fontsize=12, fontweight='bold')
cbar4 = plt.colorbar(im4, ax=ax4)
cbar4.set_label('Stress (MPa)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case2_turbine_blade.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例2完成，图片已保存")

4.3 案例3：核反应堆燃料棒热流固耦合

问题描述：核反应堆燃料棒是核裂变能量的来源，燃料芯块产生的热量通过包壳传递给冷却剂。燃料棒承受着高热梯度和冷却剂压力的联合作用，是反应堆安全分析的关键部件。

燃料棒几何：

• 燃料芯块直径：$D_f=8.2$ mm
• 包壳外径：$D_c=9.5$ mm
• 包壳厚度：$t_c=0.57$ mm
• 燃料棒长度：$L=3.8$ m
• 芯块-包壳间隙：$g=0.08$ mm

材料参数：

• 燃料（UO₂）：
  • 导热系数：$k_f=2.5$ W/(m·K)
  • 热膨胀系数：${\alpha }_f=10\times 10^{-6}$ /K
  • 弹性模量：$E_f=200$ GPa
• 包壳（Zircaloy-4）：
  • 导热系数：$k_c=13$ W/(m·K)
  • 热膨胀系数：${\alpha }_c=6.5\times 10^{-6}$ /K
  • 弹性模量：$E_c=75$ GPa

工况条件：

• 线功率密度：$q^{\prime }=40$ kW/m
• 冷却剂温度：$T_{cool}=300$°C
• 冷却剂压力：$p=15.5$ MPa
• 冷却剂流速：$v=4.5$ m/s
• 换热系数：$h=12000$ W/(m²·K)

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例3：核反应堆燃料棒热流固耦合
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Wedge
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题027'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "="*60)
print("案例3：核反应堆燃料棒热流固耦合")
print("="*60)

# 燃料棒几何参数
D_f = 8.2e-3  # 燃料芯块直径 (m)
D_c = 9.5e-3  # 包壳外径 (m)
t_c = 0.57e-3  # 包壳厚度 (m)
L = 3.8  # 燃料棒长度 (m)
g_gap = 0.08e-3  # 芯块-包壳间隙 (m)

R_f = D_f / 2  # 燃料半径
R_ci = D_c / 2 - t_c  # 包壳内半径
R_co = D_c / 2  # 包壳外半径

# 材料参数
# 燃料 (UO2)
k_f = 2.5  # 导热系数 (W/(m·K))
alpha_f = 10e-6  # 热膨胀系数 (/K)
E_f = 200e9  # 弹性模量 (Pa)
nu_f = 0.25  # 泊松比

# 包壳 (Zircaloy-4)
k_c = 13  # 导热系数 (W/(m·K))
alpha_c = 6.5e-6  # 热膨胀系数 (/K)
E_c = 75e9  # 弹性模量 (Pa)
nu_c = 0.33  # 泊松比

# 工况条件
q_prime = 40e3  # 线功率密度 (W/m)
T_cool = 300  # 冷却剂温度 (°C)
p_cool = 15.5e6  # 冷却剂压力 (Pa)
v_cool = 4.5  # 冷却剂流速 (m/s)
h_cool = 12000  # 换热系数 (W/(m²·K))

# 创建径向网格
n_r = 100
r = np.linspace(0, R_co, n_r)

# 计算温度分布（圆柱坐标稳态热传导）
# 燃料芯块：有内热源的圆柱体
# T(r) = T_center - q'''*r²/(4*k)
q_triple = q_prime / (np.pi * R_f**2)  # 体积热源 (W/m³)

T_fuel = np.zeros_like(r)
T_clad = np.zeros_like(r)

# 包壳外表面温度（由对流换热决定）
T_co = T_cool + q_prime / (2 * np.pi * R_co * h_cool)

# 包壳内表面温度
T_ci = T_co + q_prime * np.log(R_co/R_ci) / (2 * np.pi * k_c)

# 燃料表面温度
T_fs = T_ci + q_prime * g_gap / (2 * np.pi * R_f * 0.5)  # 间隙传热

# 燃料中心温度
T_center = T_fs + q_triple * R_f**2 / (4 * k_f)

# 填充温度数组
for i, ri in enumerate(r):
    if ri <= R_f:
        T_fuel[i] = T_center - q_triple * ri**2 / (4 * k_f)
    elif ri <= R_ci:
        T_fuel[i] = np.nan  # 间隙
    else:
        T_clad[i] = T_ci + q_prime * np.log(ri/R_ci) / (2 * np.pi * k_c)

T_total = np.where(np.isnan(T_fuel), T_clad, T_fuel)

# 计算热应力（厚壁圆筒理论）
# 燃料芯块热应力
sigma_rf = np.zeros_like(r)  # 径向应力
sigma_tf = np.zeros_like(r)  # 环向应力

for i, ri in enumerate(r):
    if ri <= R_f:
        # 有内热源的圆盘热应力
        sigma_rf[i] = alpha_f * E_f * q_triple / (16 * k_f * (1-nu_f)) * (R_f**2 - ri**2)
        sigma_tf[i] = alpha_f * E_f * q_triple / (16 * k_f * (1-nu_f)) * (R_f**2 - 3*ri**2)

# 包壳热应力（简化模型）
sigma_rc = np.zeros_like(r)
sigma_tc = np.zeros_like(r)

for i, ri in enumerate(r):
    if ri >= R_ci:
        delta_T_clad = T_clad[i] - T_cool
        sigma_tc[i] = alpha_c * E_c * delta_T_clad / (1 - nu_c)

# 外压引起的应力（拉梅公式）
sigma_p_outer = -p_cool * R_co**2 / (R_co**2 - R_ci**2) * (1 - R_ci**2/r**2)
sigma_p_outer = np.where(r >= R_ci, sigma_p_outer, 0)

# 组合应力
sigma_r_total = sigma_rf + sigma_rc + sigma_p_outer
sigma_t_total = sigma_tf + sigma_tc + sigma_p_outer

# von Mises等效应力
sigma_vm = np.sqrt(0.5 * ((sigma_r_total - sigma_t_total)**2 + sigma_t_total**2 + sigma_r_total**2))

print(f"燃料中心温度: {T_center:.1f}°C")
print(f"燃料表面温度: {T_fs:.1f}°C")
print(f"包壳内表面温度: {T_ci:.1f}°C")
print(f"包壳外表面温度: {T_co:.1f}°C")
print(f"最大von Mises应力: {np.max(sigma_vm)/1e6:.2f} MPa")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 燃料棒截面示意图
ax1 = axes[0, 0]
# 燃料芯块
fuel_circle = plt.Circle((0, 0), R_f*1000, fill=True, color='orange', alpha=0.7, label='Fuel')
ax1.add_patch(fuel_circle)
# 包壳内表面
clad_inner = plt.Circle((0, 0), R_ci*1000, fill=False, color='gray', linewidth=2, linestyle='--')
ax1.add_patch(clad_inner)
# 包壳外表面
clad_outer = plt.Circle((0, 0), R_co*1000, fill=True, color='silver', alpha=0.5, label='Cladding')
ax1.add_patch(clad_outer)

ax1.set_xlim(-6, 6)
ax1.set_ylim(-6, 6)
ax1.set_aspect('equal')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('y (mm)', fontsize=11)
ax1.set_title('Fuel Rod Cross Section', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)

# 2. 径向温度分布
ax2 = axes[0, 1]
ax2.plot(r*1000, T_total, 'b-', linewidth=2, label='Temperature')
ax2.axvline(x=R_f*1000, color='r', linestyle='--', label='Fuel surface')
ax2.axvline(x=R_ci*1000, color='g', linestyle='--', label='Clad inner')
ax2.axvline(x=R_co*1000, color='m', linestyle='--', label='Clad outer')
ax2.set_xlabel('Radius (mm)', fontsize=11)
ax2.set_ylabel('Temperature (°C)', fontsize=11)
ax2.set_title('Radial Temperature Distribution', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 3. 径向应力分布
ax3 = axes[1, 0]
ax3.plot(r*1000, sigma_r_total/1e6, 'r-', linewidth=2, label='Radial stress')
ax3.plot(r*1000, sigma_t_total/1e6, 'b-', linewidth=2, label='Hoop stress')
ax3.axvline(x=R_f*1000, color='gray', linestyle='--', alpha=0.5)
ax3.axvline(x=R_ci*1000, color='gray', linestyle='--', alpha=0.5)
ax3.axvline(x=R_co*1000, color='gray', linestyle='--', alpha=0.5)
ax3.set_xlabel('Radius (mm)', fontsize=11)
ax3.set_ylabel('Stress (MPa)', fontsize=11)
ax3.set_title('Radial Stress Distribution', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 4. von Mises等效应力
ax4 = axes[1, 1]
ax4.plot(r*1000, sigma_vm/1e6, 'g-', linewidth=2)
ax4.axvline(x=R_f*1000, color='gray', linestyle='--', alpha=0.5)
ax4.axvline(x=R_ci*1000, color='gray', linestyle='--', alpha=0.5)
ax4.axvline(x=R_co*1000, color='gray', linestyle='--', alpha=0.5)
ax4.set_xlabel('Radius (mm)', fontsize=11)
ax4.set_ylabel('von Mises Stress (MPa)', fontsize=11)
ax4.set_title('von Mises Equivalent Stress', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case3_fuel_rod.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例3完成，图片已保存")

# 安全裕度评估
sigma_yield_clad = 200e6  # 包壳屈服强度 (Pa)
safety_margin = sigma_yield_clad / np.max(sigma_vm)
print(f"\n安全评估:")
print(f"包壳屈服强度: {sigma_yield_clad/1e6:.1f} MPa")
print(f"最大等效应力: {np.max(sigma_vm)/1e6:.2f} MPa")
print(f"安全裕度: {safety_margin:.2f}")

4.4 案例4：电子器件散热变形分析

问题描述：现代电子器件功率密度不断提高，芯片在工作时产生大量热量。温度升高导致PCB板和封装材料的热膨胀，可能引起焊点失效、翘曲变形等问题。准确预测电子器件的热-结构响应对于可靠性设计至关重要。

器件几何：

• 芯片尺寸：$10\times 10\times 0.5$ mm³
• 基板尺寸：$50\times 50\times 1.6$ mm³
• 焊球直径：$0.4$ mm
• 焊球数量：$16\times 16=256$
• 焊球间距：$1.0$ mm

材料参数：

• 芯片（硅）：
  • 弹性模量：$E=130$ GPa
  • 热膨胀系数：$\alpha =2.6\times 10^{-6}$ /K
  • 导热系数：$k=150$ W/(m·K)
• 基板（FR4）：
  • 弹性模量：$E=22$ GPa
  • 热膨胀系数：$\alpha =14\times 10^{-6}$ /K
  • 导热系数：$k=0.3$ W/(m·K)
• 焊球（SAC305）：
  • 弹性模量：$E=50$ GPa
  • 热膨胀系数：$\alpha =21\times 10^{-6}$ /K

工况条件：

• 芯片功耗：$P=15$ W
• 环境温度：$T_{amb}=25$°C
• 对流换热系数：$h=50$ W/(m²·K)（自然对流）
• 散热器热阻：$R_{hs}=2$ K/W

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例4：电子器件散热变形分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Circle, FancyBboxPatch
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题027'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "="*60)
print("案例4：电子器件散热变形分析")
print("="*60)

# 器件几何参数
L_chip = 10e-3  # 芯片边长 (m)
H_chip = 0.5e-3  # 芯片厚度 (m)
L_sub = 50e-3  # 基板边长 (m)
H_sub = 1.6e-3  # 基板厚度 (m)
d_ball = 0.4e-3  # 焊球直径 (m)
n_balls = 16  # 每边焊球数
P_ball = 1.0e-3  # 焊球间距 (m)

# 材料参数
# 芯片 (硅)
E_chip = 130e9  # 弹性模量 (Pa)
alpha_chip = 2.6e-6  # 热膨胀系数 (/K)
k_chip = 150  # 导热系数 (W/(m·K))
nu_chip = 0.28  # 泊松比

# 基板 (FR4)
E_sub = 22e9  # 弹性模量 (Pa)
alpha_sub = 14e-6  # 热膨胀系数 (/K)
k_sub = 0.3  # 导热系数 (W/(m·K))
nu_sub = 0.3  # 泊松比

# 焊球 (SAC305)
E_ball = 50e9  # 弹性模量 (Pa)
alpha_ball = 21e-6  # 热膨胀系数 (/K)
nu_ball = 0.35  # 泊松比

# 工况条件
P_chip = 15  # 芯片功耗 (W)
T_amb = 25  # 环境温度 (°C)
h_conv = 50  # 对流换热系数 (W/(m²·K))
R_hs = 2  # 散热器热阻 (K/W)

# 计算温度分布（简化热阻网络模型）
# 芯片结温
T_junction = T_amb + P_chip * R_hs

# 芯片表面温度分布（考虑热扩散）
n_grid = 50
x = np.linspace(-L_sub/2, L_sub/2, n_grid)
y = np.linspace(-L_sub/2, L_sub/2, n_grid)
X, Y = np.meshgrid(x, y)

# 温度分布：芯片区域温度高，向外扩散
T_board = np.zeros_like(X)
for i in range(n_grid):
    for j in range(n_grid):
        dist = np.sqrt(X[i,j]**2 + Y[i,j]**2)
        if dist < L_chip/2:
            # 芯片区域
            T_board[i,j] = T_junction - 10 * (dist / (L_chip/2))**2
        else:
            # 基板区域，指数衰减
            T_board[i,j] = T_junction * np.exp(-(dist - L_chip/2) / 0.01) + T_amb

# 计算热变形（简化模型）
# 基板翘曲主要由CTE失配引起
delta_alpha = alpha_sub - alpha_chip
delta_T_avg = np.mean(T_board) - T_amb

# 基板翘曲（简支梁模型）
# 假设基板中心固定，边缘自由
kappa = 6 * delta_alpha * delta_T_avg / H_sub  # 曲率

# 计算变形场
w_deform = np.zeros_like(X)
for i in range(n_grid):
    for j in range(n_grid):
        r = np.sqrt(X[i,j]**2 + Y[i,j]**2)
        w_deform[i,j] = 0.5 * kappa * r**2

# 焊球剪切应变（由CTE失配引起）
# 最外层焊球承受的剪切最大
r_max = (n_balls - 1) / 2 * P_ball
gamma_ball = delta_alpha * delta_T_avg * r_max / H_sub

# 焊球剪切应力
tau_ball = E_ball / (2 * (1 + nu_ball)) * gamma_ball

print(f"芯片结温: {T_junction:.1f}°C")
print(f"基板平均温升: {delta_T_avg:.1f}°C")
print(f"最大翘曲变形: {np.max(w_deform)*1e6:.2f} μm")
print(f"焊球最大剪切应变: {gamma_ball*100:.3f}%")
print(f"焊球最大剪切应力: {tau_ball/1e6:.2f} MPa")

# 绘制结果
fig = plt.figure(figsize=(16, 12))

# 1. 器件结构示意图
ax1 = fig.add_subplot(2, 2, 1)
# 基板
sub_rect = Rectangle((-L_sub/2*1000, -H_sub*1000), L_sub*1000, H_sub*1000, 
                      fill=True, color='green', alpha=0.5, label='Substrate')
ax1.add_patch(sub_rect)
# 芯片
chip_rect = Rectangle((-L_chip/2*1000, 0), L_chip*1000, H_chip*1000, 
                       fill=True, color='blue', alpha=0.7, label='Chip')
ax1.add_patch(chip_rect)
# 焊球
for i in range(n_balls):
    for j in range(n_balls):
        x_ball = (-(n_balls-1)/2 + i) * P_ball * 1000
        y_ball = (-(n_balls-1)/2 + j) * P_ball * 1000
        if abs(x_ball) <= L_chip/2*1000 and abs(y_ball) <= L_chip/2*1000:
            ball = Circle((x_ball, -d_ball/2*1000), d_ball/2*1000, 
                         fill=True, color='gray', alpha=0.8)
            ax1.add_patch(ball)

ax1.set_xlim(-30, 30)
ax1.set_ylim(-3, 2)
ax1.set_aspect('equal')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('z (mm)', fontsize=11)
ax1.set_title('Electronic Package Cross Section', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)

# 2. 基板温度分布
ax2 = fig.add_subplot(2, 2, 2)
contour2 = ax2.contourf(X*1000, Y*1000, T_board, levels=20, cmap='hot')
ax2.set_aspect('equal')
ax2.set_xlabel('x (mm)', fontsize=11)
ax2.set_ylabel('y (mm)', fontsize=11)
ax2.set_title('PCB Temperature Distribution (°C)', fontsize=12, fontweight='bold')
cbar2 = plt.colorbar(contour2, ax=ax2)
cbar2.set_label('Temperature (°C)', fontsize=10)

# 3. 基板翘曲变形（3D）
ax3 = fig.add_subplot(2, 2, 3, projection='3d')
surf3 = ax3.plot_surface(X*1000, Y*1000, w_deform*1e6, cmap='viridis', 
                         alpha=0.8, rstride=2, cstride=2)
ax3.set_xlabel('x (mm)', fontsize=10)
ax3.set_ylabel('y (mm)', fontsize=10)
ax3.set_zlabel('Deformation (μm)', fontsize=10)
ax3.set_title('PCB Warpage Deformation', fontsize=12, fontweight='bold')

# 4. 焊球应力分布
ax4 = fig.add_subplot(2, 2, 4)
# 计算每个焊球的剪切应力
ball_stresses = []
ball_positions = []
for i in range(n_balls):
    for j in range(n_balls):
        x_ball = (-(n_balls-1)/2 + i) * P_ball
        y_ball = (-(n_balls-1)/2 + j) * P_ball
        if abs(x_ball) <= L_chip/2 and abs(y_ball) <= L_chip/2:
            r_ball = np.sqrt(x_ball**2 + y_ball**2)
            tau = tau_ball * r_ball / r_max
            ball_stresses.append(tau)
            ball_positions.append((x_ball*1000, y_ball*1000))

ball_stresses = np.array(ball_stresses)
scatter4 = ax4.scatter([p[0] for p in ball_positions], 
                       [p[1] for p in ball_positions], 
                       c=ball_stresses/1e6, s=100, cmap='jet', vmin=0)
ax4.set_aspect('equal')
ax4.set_xlabel('x (mm)', fontsize=11)
ax4.set_ylabel('y (mm)', fontsize=11)
ax4.set_title('Solder Ball Shear Stress (MPa)', fontsize=12, fontweight='bold')
cbar4 = plt.colorbar(scatter4, ax=ax4)
cbar4.set_label('Shear Stress (MPa)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case4_electronics.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例4完成，图片已保存")

# 可靠性评估
# 焊球疲劳寿命估算（基于Coffin-Manson方程）
delta_gamma = gamma_ball  # 剪切应变幅
N_f_solder = 0.5 * (0.4 / delta_gamma)**(1/0.6)  # 疲劳寿命循环数
print(f"\n可靠性评估:")
print(f"焊球估算疲劳寿命: {N_f_solder:.0f} 热循环")
if N_f_solder > 1000:
    print("✓ 焊球设计满足可靠性要求")
else:
    print("⚠ 焊球疲劳寿命不足，建议优化设计")

4.5 案例5：汽车排气系统热应力分析

问题描述：汽车排气系统在工作时承受高温废气的热冲击和发动机振动载荷。排气歧管、催化转化器和消声器等部件的温度可达800°C以上，热应力和热疲劳是主要的失效模式。

排气歧管几何：

• 歧管长度：$L=0.4$ m
• 管径：$D=0.045$ m
• 壁厚：$t=0.003$ m
• 歧管分支数：$N_b=4$
• 法兰厚度：$t_f=0.015$ m

材料参数（不锈钢）：

• 弹性模量：$E=200$ GPa（随温度变化）
• 泊松比：$\nu =0.3$
• 热膨胀系数：$\alpha =16\times 10^{-6}$ /K
• 导热系数：$k=25$ W/(m·K)
• 密度：$\rho =7850$ kg/m³

工况条件：

• 废气温度：$T_{ex}=850$°C（峰值）
• 废气温度：$T_{ex}=600$°C（稳态）
• 废气压力：$p_{ex}=0.15$ MPa
• 环境温度：$T_{amb}=30$°C
• 对流换热系数（内侧）：$h_i=200$ W/(m²·K)
• 对流换热系数（外侧）：$h_o=30$ W/(m²·K)

Python仿真代码：

# -*- coding: utf-8 -*-
"""
主题027：热-流-固耦合分析
Python仿真代码 - 包含6个工程案例
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle, FancyBboxPatch
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题027'
os.makedirs(output_dir, exist_ok=True)

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

print("="*60)
print("主题027：热-流-固耦合分析")
print("="*60)

# ============================================================
# 案例1：热交换器管板应力分析
# ============================================================
print("\n" + "="*60)
print("案例1：热交换器管板应力分析")
print("="*60)

# 几何参数
D = 0.6  # 管板直径 (m)
t = 0.05  # 管板厚度 (m)
N_tubes = 100  # 管孔数量
d_tube = 0.025  # 管孔直径 (m)
P = 0.04  # 管间距 (m)

# 材料参数
E = 200e9  # 弹性模量 (Pa)
nu = 0.3  # 泊松比
alpha = 12e-6  # 热膨胀系数 (/K)
k = 45  # 导热系数 (W/(m·K))

# 工况条件
T_tube = 300  # 管程温度 (°C)
T_shell = 150  # 壳程温度 (°C)
p_tube = 5e6  # 管程压力 (Pa)
p_shell = 2e6  # 壳程压力 (Pa)
T_ref = 20  # 参考温度 (°C)

# 创建管板网格
n_r = 50
n_theta = 60
r = np.linspace(0, D/2, n_r)
theta = np.linspace(0, 2*np.pi, n_theta)
R, Theta = np.meshgrid(r, theta)
X = R * np.cos(Theta)
Y = R * np.sin(Theta)

# 生成管孔位置（三角形排列）
n_rings = int((D/2) / P) + 1
tube_positions = []
for i in range(n_rings):
    r_ring = i * P
    if r_ring > D/2 - d_tube:
        break
    n_tubes_ring = max(1, int(2 * np.pi * r_ring / P))
    for j in range(n_tubes_ring):
        angle = 2 * np.pi * j / n_tubes_ring
        if i % 2 == 1:
            angle += np.pi / n_tubes_ring
        x_tube = r_ring * np.cos(angle)
        y_tube = r_ring * np.sin(angle)
        if np.sqrt(x_tube**2 + y_tube**2) < D/2 - d_tube:
            tube_positions.append((x_tube, y_tube))

tube_positions = tube_positions[:N_tubes]
print(f"实际布置管孔数量: {len(tube_positions)}")

# 计算温度场
T_field = np.zeros_like(X)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        x, y = X[i,j], Y[i,j]
        min_dist = float('inf')
        for tx, ty in tube_positions:
            dist = np.sqrt((x-tx)**2 + (y-ty)**2)
            if dist < min_dist:
                min_dist = dist
        
        if min_dist < d_tube/2:
            T_field[i,j] = T_tube
        else:
            T_field[i,j] = T_shell + (T_tube - T_shell) * np.exp(-min_dist / 0.05)

# 计算热应力
sigma_th = E * alpha * (T_field - T_ref) / (1 - nu)
delta_p = p_tube - p_shell
sigma_p = delta_p * D / (4 * t)
K_t = 3.0
sigma_total = K_t * (sigma_th + sigma_p)
sigma_max = 300e6
sigma_total = np.clip(sigma_total, 0, sigma_max * 1.5)

print(f"管板最大热应力: {np.max(sigma_th)/1e6:.2f} MPa")
print(f"压力引起的应力: {sigma_p/1e6:.2f} MPa")
print(f"组合最大应力: {np.max(sigma_total)/1e6:.2f} MPa")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

ax1 = axes[0, 0]
circle = plt.Circle((0, 0), D/2, fill=False, color='black', linewidth=2)
ax1.add_patch(circle)
for tx, ty in tube_positions:
    hole = plt.Circle((tx, ty), d_tube/2, fill=True, color='lightblue', alpha=0.7)
    ax1.add_patch(hole)
ax1.set_xlim(-D/2*1.1, D/2*1.1)
ax1.set_ylim(-D/2*1.1, D/2*1.1)
ax1.set_aspect('equal')
ax1.set_xlabel('x (m)', fontsize=11)
ax1.set_ylabel('y (m)', fontsize=11)
ax1.set_title('Tube Sheet Geometry', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
contour2 = ax2.contourf(X, Y, T_field, levels=20, cmap='hot')
ax2.set_aspect('equal')
ax2.set_xlabel('x (m)', fontsize=11)
ax2.set_ylabel('y (m)', fontsize=11)
ax2.set_title('Temperature Distribution (°C)', fontsize=12, fontweight='bold')
cbar2 = plt.colorbar(contour2, ax=ax2)
cbar2.set_label('Temperature (°C)', fontsize=10)

ax3 = axes[1, 0]
contour3 = ax3.contourf(X, Y, sigma_th/1e6, levels=20, cmap='jet')
ax3.set_aspect('equal')
ax3.set_xlabel('x (m)', fontsize=11)
ax3.set_ylabel('y (m)', fontsize=11)
ax3.set_title('Thermal Stress (MPa)', fontsize=12, fontweight='bold')
cbar3 = plt.colorbar(contour3, ax=ax3)
cbar3.set_label('Stress (MPa)', fontsize=10)

ax4 = axes[1, 1]
contour4 = ax4.contourf(X, Y, sigma_total/1e6, levels=20, cmap='jet')
ax4.set_aspect('equal')
ax4.set_xlabel('x (m)', fontsize=11)
ax4.set_ylabel('y (m)', fontsize=11)
ax4.set_title('Total Stress (MPa)', fontsize=12, fontweight='bold')
cbar4 = plt.colorbar(contour4, ax=ax4)
cbar4.set_label('Stress (MPa)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_heat_exchanger.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例1完成，图片已保存")

safety_factor = sigma_max / np.max(sigma_total)
print(f"\n安全评估: 安全系数 = {safety_factor:.2f}")

# ============================================================
# 案例2：航空发动机叶片热疲劳分析
# ============================================================
print("\n" + "="*60)
print("案例2：航空发动机叶片热疲劳分析")
print("="*60)

L_blade = 0.12
C_root = 0.04
C_tip = 0.025
t_max = 0.006
N_c = 5

E_0 = 120e9
nu = 0.3
alpha = 14e-6
k = 15
rho = 8200

T_g = 1400
T_c = 600
omega = 15000 * 2 * np.pi / 60
h_g = 2000
h_c = 800

n_span = 50
n_chord = 60
span = np.linspace(0, L_blade, n_span)
chord = np.linspace(0, 1, n_chord)

def airfoil_shape(x_c, t_c, chord_len):
    y_t = 5 * t_c * (0.2969*np.sqrt(x_c) - 0.1260*x_c - 0.3516*x_c**2 + 
                     0.2843*x_c**3 - 0.1015*x_c**4)
    return y_t * chord_len

T_blade = np.zeros((n_span, n_chord))
for i, s in enumerate(span):
    C_local = C_root + (C_tip - C_root) * s / L_blade
    T_g_local = T_g - 50 * s / L_blade
    
    for j, x_c in enumerate(chord):
        y_t = airfoil_shape(x_c, 0.12, 1.0)
        Bi_g = h_g * t_max / k
        Bi_c = h_c * t_max / k
        T_surface = (Bi_g * T_g_local + Bi_c * T_c) / (Bi_g + Bi_c)
        T_blade[i, j] = T_surface - (T_surface - T_c) * y_t * 2

T_avg = np.mean(T_blade)
delta_T = T_blade - T_avg
E_T = E_0 * (1 - 0.0003 * (T_blade - 20))
E_T = np.clip(E_T, E_0 * 0.5, E_0)
sigma_th_blade = E_T * alpha * delta_T / (1 - nu)

radius = 0.2 + span
sigma_centrifugal = rho * omega**2 * radius[:, np.newaxis] * C_root / 2
sigma_centrifugal = np.tile(sigma_centrifugal, (1, n_chord))
sigma_total_blade = sigma_th_blade + sigma_centrifugal

delta_epsilon = alpha * (np.max(T_blade) - np.min(T_blade))
epsilon_f = 0.5
c = -0.6
N_f = (epsilon_f / delta_epsilon)**(1/abs(c))

print(f"叶片最高温度: {np.max(T_blade):.1f}°C")
print(f"叶片最低温度: {np.min(T_blade):.1f}°C")
print(f"最大热应力: {np.max(np.abs(sigma_th_blade))/1e6:.2f} MPa")
print(f"最大离心应力: {np.max(sigma_centrifugal)/1e6:.2f} MPa")
print(f"最大组合应力: {np.max(sigma_total_blade)/1e6:.2f} MPa")
print(f"估算热疲劳寿命: {N_f:.0f} 循环")

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

ax1 = axes[0, 0]
s_plot = [0, L_blade/2, L_blade]
for idx, s in enumerate(s_plot):
    C_local = C_root + (C_tip - C_root) * s / L_blade
    x_af = np.linspace(0, 1, 50) * C_local
    y_upper = airfoil_shape(np.linspace(0, 1, 50), 0.12, C_local)
    y_lower = -y_upper
    
    offset = idx * C_root * 1.5
    ax1.plot(x_af + offset, y_upper + s*1000, 'b-', linewidth=2)
    ax1.plot(x_af + offset, y_lower + s*1000, 'b-', linewidth=2)
    ax1.fill_between(x_af + offset, y_lower + s*1000, y_upper + s*1000, alpha=0.3)
    
    for k in range(N_c):
        x_cool = C_local * (0.2 + 0.6 * k / (N_c-1)) + offset
        y_cool = s * 1000
        circle = plt.Circle((x_cool, y_cool), 0.002, fill=True, color='cyan', alpha=0.8)
        ax1.add_patch(circle)

ax1.set_xlabel('Chord Position (mm)', fontsize=11)
ax1.set_ylabel('Span Position (mm)', fontsize=11)
ax1.set_title('Turbine Blade Geometry', fontsize=12, fontweight='bold')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
Span, Chord = np.meshgrid(span*1000, chord)
im2 = ax2.contourf(Span, Chord, T_blade.T, levels=20, cmap='hot')
ax2.set_xlabel('Span (mm)', fontsize=11)
ax2.set_ylabel('Normalized Chord', fontsize=11)
ax2.set_title('Temperature Distribution (°C)', fontsize=12, fontweight='bold')
cbar2 = plt.colorbar(im2, ax=ax2)
cbar2.set_label('Temperature (°C)', fontsize=10)

ax3 = axes[1, 0]
im3 = ax3.contourf(Span, Chord, sigma_th_blade.T/1e6, levels=20, cmap='RdBu_r')
ax3.set_xlabel('Span (mm)', fontsize=11)
ax3.set_ylabel('Normalized Chord', fontsize=11)
ax3.set_title('Thermal Stress (MPa)', fontsize=12, fontweight='bold')
cbar3 = plt.colorbar(im3, ax=ax3)
cbar3.set_label('Stress (MPa)', fontsize=10)

ax4 = axes[1, 1]
im4 = ax4.contourf(Span, Chord, sigma_total_blade.T/1e6, levels=20, cmap='jet')
ax4.set_xlabel('Span (mm)', fontsize=11)
ax4.set_ylabel('Normalized Chord', fontsize=11)
ax4.set_title('Total Stress (MPa)', fontsize=12, fontweight='bold')
cbar4 = plt.colorbar(im4, ax=ax4)
cbar4.set_label('Stress (MPa)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case2_turbine_blade.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例2完成，图片已保存")

# ============================================================
# 案例3：核反应堆燃料棒热流固耦合
# ============================================================
print("\n" + "="*60)
print("案例3：核反应堆燃料棒热流固耦合")
print("="*60)

D_f = 8.2e-3
D_c = 9.5e-3
t_c = 0.57e-3
L = 3.8
g_gap = 0.08e-3

R_f = D_f / 2
R_ci = D_c / 2 - t_c
R_co = D_c / 2

k_f = 2.5
alpha_f = 10e-6
E_f = 200e9
nu_f = 0.25

k_c = 13
alpha_c = 6.5e-6
E_c = 75e9
nu_c = 0.33

q_prime = 40e3
T_cool = 300
p_cool = 15.5e6
v_cool = 4.5
h_cool = 12000

n_r = 100
r = np.linspace(0, R_co, n_r)

q_triple = q_prime / (np.pi * R_f**2)

T_fuel = np.zeros_like(r)
T_clad = np.zeros_like(r)

T_co = T_cool + q_prime / (2 * np.pi * R_co * h_cool)
T_ci = T_co + q_prime * np.log(R_co/R_ci) / (2 * np.pi * k_c)
T_fs = T_ci + q_prime * g_gap / (2 * np.pi * R_f * 0.5)
T_center = T_fs + q_triple * R_f**2 / (4 * k_f)

for i, ri in enumerate(r):
    if ri <= R_f:
        T_fuel[i] = T_center - q_triple * ri**2 / (4 * k_f)
    elif ri <= R_ci:
        T_fuel[i] = np.nan
    else:
        T_clad[i] = T_ci + q_prime * np.log(ri/R_ci) / (2 * np.pi * k_c)

T_total = np.where(np.isnan(T_fuel), T_clad, T_fuel)

sigma_rf = np.zeros_like(r)
sigma_tf = np.zeros_like(r)

for i, ri in enumerate(r):
    if ri <= R_f:
        sigma_rf[i] = alpha_f * E_f * q_triple / (16 * k_f * (1-nu_f)) * (R_f**2 - ri**2)
        sigma_tf[i] = alpha_f * E_f * q_triple / (16 * k_f * (1-nu_f)) * (R_f**2 - 3*ri**2)

sigma_rc = np.zeros_like(r)
sigma_tc = np.zeros_like(r)

for i, ri in enumerate(r):
    if ri >= R_ci:
        delta_T_clad = T_clad[i] - T_cool
        sigma_tc[i] = alpha_c * E_c * delta_T_clad / (1 - nu_c)

sigma_p_outer = -p_cool * R_co**2 / (R_co**2 - R_ci**2) * (1 - R_ci**2/r**2)
sigma_p_outer = np.where(r >= R_ci, sigma_p_outer, 0)

sigma_r_total = sigma_rf + sigma_rc + sigma_p_outer
sigma_t_total = sigma_tf + sigma_tc + sigma_p_outer
sigma_vm = np.sqrt(0.5 * ((sigma_r_total - sigma_t_total)**2 + sigma_t_total**2 + sigma_r_total**2))

print(f"燃料中心温度: {T_center:.1f}°C")
print(f"燃料表面温度: {T_fs:.1f}°C")
print(f"包壳内表面温度: {T_ci:.1f}°C")
print(f"包壳外表面温度: {T_co:.1f}°C")
print(f"最大von Mises应力: {np.max(sigma_vm)/1e6:.2f} MPa")

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

ax1 = axes[0, 0]
fuel_circle = plt.Circle((0, 0), R_f*1000, fill=True, color='orange', alpha=0.7, label='Fuel')
ax1.add_patch(fuel_circle)
clad_inner = plt.Circle((0, 0), R_ci*1000, fill=False, color='gray', linewidth=2, linestyle='--')
ax1.add_patch(clad_inner)
clad_outer = plt.Circle((0, 0), R_co*1000, fill=True, color='silver', alpha=0.5, label='Cladding')
ax1.add_patch(clad_outer)
ax1.set_xlim(-6, 6)
ax1.set_ylim(-6, 6)
ax1.set_aspect('equal')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('y (mm)', fontsize=11)
ax1.set_title('Fuel Rod Cross Section', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
ax2.plot(r*1000, T_total, 'b-', linewidth=2, label='Temperature')
ax2.axvline(x=R_f*1000, color='r', linestyle='--', label='Fuel surface')
ax2.axvline(x=R_ci*1000, color='g', linestyle='--', label='Clad inner')
ax2.axvline(x=R_co*1000, color='m', linestyle='--', label='Clad outer')
ax2.set_xlabel('Radius (mm)', fontsize=11)
ax2.set_ylabel('Temperature (°C)', fontsize=11)
ax2.set_title('Radial Temperature Distribution', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
ax3.plot(r*1000, sigma_r_total/1e6, 'r-', linewidth=2, label='Radial stress')
ax3.plot(r*1000, sigma_t_total/1e6, 'b-', linewidth=2, label='Hoop stress')
ax3.axvline(x=R_f*1000, color='gray', linestyle='--', alpha=0.5)
ax3.axvline(x=R_ci*1000, color='gray', linestyle='--', alpha=0.5)
ax3.axvline(x=R_co*1000, color='gray', linestyle='--', alpha=0.5)
ax3.set_xlabel('Radius (mm)', fontsize=11)
ax3.set_ylabel('Stress (MPa)', fontsize=11)
ax3.set_title('Radial Stress Distribution', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
ax4.plot(r*1000, sigma_vm/1e6, 'g-', linewidth=2)
ax4.axvline(x=R_f*1000, color='gray', linestyle='--', alpha=0.5)
ax4.axvline(x=R_ci*1000, color='gray', linestyle='--', alpha=0.5)
ax4.axvline(x=R_co*1000, color='gray', linestyle='--', alpha=0.5)
ax4.set_xlabel('Radius (mm)', fontsize=11)
ax4.set_ylabel('von Mises Stress (MPa)', fontsize=11)
ax4.set_title('von Mises Equivalent Stress', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case3_fuel_rod.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例3完成，图片已保存")

sigma_yield_clad = 200e6
safety_margin = sigma_yield_clad / np.max(sigma_vm)
print(f"安全裕度: {safety_margin:.2f}")

# ============================================================
# 案例4：电子器件散热变形分析
# ============================================================
print("\n" + "="*60)
print("案例4：电子器件散热变形分析")
print("="*60)

L_chip = 10e-3
H_chip = 0.5e-3
L_sub = 50e-3
H_sub = 1.6e-3
d_ball = 0.4e-3
n_balls = 16
P_ball = 1.0e-3

E_chip = 130e9
alpha_chip = 2.6e-6
k_chip = 150
nu_chip = 0.28

E_sub = 22e9
alpha_sub = 14e-6
k_sub = 0.3
nu_sub = 0.3

E_ball = 50e9
alpha_ball = 21e-6
nu_ball = 0.35

P_chip = 15
T_amb = 25
h_conv = 50
R_hs = 2

T_junction = T_amb + P_chip * R_hs

n_grid = 50
x = np.linspace(-L_sub/2, L_sub/2, n_grid)
y = np.linspace(-L_sub/2, L_sub/2, n_grid)
X, Y = np.meshgrid(x, y)

T_board = np.zeros_like(X)
for i in range(n_grid):
    for j in range(n_grid):
        dist = np.sqrt(X[i,j]**2 + Y[i,j]**2)
        if dist < L_chip/2:
            T_board[i,j] = T_junction - 10 * (dist / (L_chip/2))**2
        else:
            T_board[i,j] = T_junction * np.exp(-(dist - L_chip/2) / 0.01) + T_amb

delta_alpha = alpha_sub - alpha_chip
delta_T_avg = np.mean(T_board) - T_amb
kappa = 6 * delta_alpha * delta_T_avg / H_sub

w_deform = np.zeros_like(X)
for i in range(n_grid):
    for j in range(n_grid):
        r = np.sqrt(X[i,j]**2 + Y[i,j]**2)
        w_deform[i,j] = 0.5 * kappa * r**2

r_max = (n_balls - 1) / 2 * P_ball
gamma_ball = delta_alpha * delta_T_avg * r_max / H_sub
tau_ball = E_ball / (2 * (1 + nu_ball)) * gamma_ball

print(f"芯片结温: {T_junction:.1f}°C")
print(f"基板平均温升: {delta_T_avg:.1f}°C")
print(f"最大翘曲变形: {np.max(w_deform)*1e6:.2f} μm")
print(f"焊球最大剪切应变: {gamma_ball*100:.3f}%")
print(f"焊球最大剪切应力: {tau_ball/1e6:.2f} MPa")

fig = plt.figure(figsize=(16, 12))

ax1 = fig.add_subplot(2, 2, 1)
sub_rect = Rectangle((-L_sub/2*1000, -H_sub*1000), L_sub*1000, H_sub*1000, 
                      fill=True, color='green', alpha=0.5, label='Substrate')
ax1.add_patch(sub_rect)
chip_rect = Rectangle((-L_chip/2*1000, 0), L_chip*1000, H_chip*1000, 
                       fill=True, color='blue', alpha=0.7, label='Chip')
ax1.add_patch(chip_rect)
for i in range(n_balls):
    for j in range(n_balls):
        x_ball = (-(n_balls-1)/2 + i) * P_ball * 1000
        y_ball = (-(n_balls-1)/2 + j) * P_ball * 1000
        if abs(x_ball) <= L_chip/2*1000 and abs(y_ball) <= L_chip/2*1000:
            ball = Circle((x_ball, -d_ball/2*1000), d_ball/2*1000, 
                         fill=True, color='gray', alpha=0.8)
            ax1.add_patch(ball)
ax1.set_xlim(-30, 30)
ax1.set_ylim(-3, 2)
ax1.set_aspect('equal')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('z (mm)', fontsize=11)
ax1.set_title('Electronic Package Cross Section', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)

ax2 = fig.add_subplot(2, 2, 2)
contour2 = ax2.contourf(X*1000, Y*1000, T_board, levels=20, cmap='hot')
ax2.set_aspect('equal')
ax2.set_xlabel('x (mm)', fontsize=11)
ax2.set_ylabel('y (mm)', fontsize=11)
ax2.set_title('PCB Temperature Distribution (°C)', fontsize=12, fontweight='bold')
cbar2 = plt.colorbar(contour2, ax=ax2)
cbar2.set_label('Temperature (°C)', fontsize=10)

ax3 = fig.add_subplot(2, 2, 3, projection='3d')
surf3 = ax3.plot_surface(X*1000, Y*1000, w_deform*1e6, cmap='viridis', 
                         alpha=0.8, rstride=2, cstride=2)
ax3.set_xlabel('x (mm)', fontsize=10)
ax3.set_ylabel('y (mm)', fontsize=10)
ax3.set_zlabel('Deformation (μm)', fontsize=10)
ax3.set_title('PCB Warpage Deformation', fontsize=12, fontweight='bold')

ax4 = fig.add_subplot(2, 2, 4)
ball_stresses = []
ball_positions = []
for i in range(n_balls):
    for j in range(n_balls):
        x_ball = (-(n_balls-1)/2 + i) * P_ball
        y_ball = (-(n_balls-1)/2 + j) * P_ball
        if abs(x_ball) <= L_chip/2 and abs(y_ball) <= L_chip/2:
            r_ball = np.sqrt(x_ball**2 + y_ball**2)
            tau = tau_ball * r_ball / r_max
            ball_stresses.append(tau)
            ball_positions.append((x_ball*1000, y_ball*1000))

ball_stresses = np.array(ball_stresses)
scatter4 = ax4.scatter([p[0] for p in ball_positions], 
                       [p[1] for p in ball_positions], 
                       c=ball_stresses/1e6, s=100, cmap='jet', vmin=0)
ax4.set_aspect('equal')
ax4.set_xlabel('x (mm)', fontsize=11)
ax4.set_ylabel('y (mm)', fontsize=11)
ax4.set_title('Solder Ball Shear Stress (MPa)', fontsize=12, fontweight='bold')
cbar4 = plt.colorbar(scatter4, ax=ax4)
cbar4.set_label('Shear Stress (MPa)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case4_electronics.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例4完成，图片已保存")

delta_gamma = gamma_ball
N_f_solder = 0.5 * (0.4 / delta_gamma)**(1/0.6)
print(f"焊球估算疲劳寿命: {N_f_solder:.0f} 热循环")

# ============================================================
# 案例5：汽车排气系统热应力分析
# ============================================================
print("\n" + "="*60)
print("案例5：汽车排气系统热应力分析")
print("="*60)

L_manifold = 0.4
D_pipe = 0.045
t_wall = 0.003
N_branches = 4
t_flange = 0.015

E_steel = 200e9
nu_steel = 0.3
alpha_steel = 16e-6
k_steel = 25
rho_steel = 7850

T_exhaust_peak = 850
T_exhaust_steady = 600
p_exhaust = 0.15e6
T_amb = 30
h_inner = 200
h_outer = 30

n_x = 50
n_theta = 36
x = np.linspace(0, L_manifold, n_x)
theta = np.linspace(0, 2*np.pi, n_theta)
X, Theta = np.meshgrid(x, theta)

r_inner = D_pipe / 2
r_outer = r_inner + t_wall

T_inner = T_exhaust_peak - (T_exhaust_peak - T_amb) / (1 + h_inner/h_outer + h_inner*t_wall/k_steel)
T_outer = T_inner - h_inner * (T_exhaust_peak - T_amb) * t_wall / k_steel

T_wall = np.zeros((n_theta, n_x))
for i in range(n_theta):
    for j in range(n_x):
        T_local = T_inner * np.exp(-x[j] / 0.5) + T_amb
        T_wall[i, j] = T_local

delta_T_wall = T_wall - T_amb
sigma_th_manifold = E_steel * alpha_steel * delta_T_wall / (1 - nu_steel)
sigma_p_hoop = p_exhaust * r_inner / t_wall
sigma_p_long = p_exhaust * r_inner / (2 * t_wall)
sigma_vm_manifold = np.sqrt(sigma_th_manifold**2 + sigma_p_hoop**2 - sigma_th_manifold*sigma_p_hoop)

print(f"歧管内壁温度: {T_inner:.1f}°C")
print(f"歧管外壁温度: {T_outer:.1f}°C")
print(f"最大热应力: {np.max(sigma_th_manifold)/1e6:.2f} MPa")
print(f"内压环向应力: {sigma_p_hoop/1e6:.2f} MPa")
print(f"最大von Mises应力: {np.max(sigma_vm_manifold)/1e6:.2f} MPa")

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

ax1 = axes[0, 0]
for i in range(N_branches):
    angle = i * 2 * np.pi / N_branches
    x_branch = np.linspace(0, L_manifold*1000, 20)
    y_branch = np.sin(angle) * 20
    ax1.plot(x_branch, [y_branch]*20, 'b-', linewidth=8, alpha=0.7)
    ax1.plot([0, 0], [0, y_branch], 'b-', linewidth=6, alpha=0.7)

flange = Rectangle((-10, -30), 20, 60, fill=True, color='gray', alpha=0.5)
ax1.add_patch(flange)
ax1.set_xlim(-20, L_manifold*1000 + 20)
ax1.set_ylim(-40, 40)
ax1.set_xlabel('Length (mm)', fontsize=11)
ax1.set_ylabel('Radial Position (mm)', fontsize=11)
ax1.set_title('Exhaust Manifold Geometry', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.set_aspect('equal')

ax2 = axes[0, 1]
contour2 = ax2.contourf(X*1000, Theta*180/np.pi, T_wall, levels=20, cmap='hot')
ax2.set_xlabel('Length (mm)', fontsize=11)
ax2.set_ylabel('Circumferential Angle (deg)', fontsize=11)
ax2.set_title('Wall Temperature Distribution (°C)', fontsize=12, fontweight='bold')
cbar2 = plt.colorbar(contour2, ax=ax2)
cbar2.set_label('Temperature (°C)', fontsize=10)

ax3 = axes[1, 0]
contour3 = ax3.contourf(X*1000, Theta*180/np.pi, sigma_th_manifold/1e6, levels=20, cmap='jet')
ax3.set_xlabel('Length (mm)', fontsize=11)
ax3.set_ylabel('Circumferential Angle (deg)', fontsize=11)
ax3.set_title('Thermal Stress (MPa)', fontsize=12, fontweight='bold')
cbar3 = plt.colorbar(contour3, ax=ax3)
cbar3.set_label('Stress (MPa)', fontsize=10)

ax4 = axes[1, 1]
contour4 = ax4.contourf(X*1000, Theta*180/np.pi, sigma_vm_manifold/1e6, levels=20, cmap='jet')
ax4.set_xlabel('Length (mm)', fontsize=11)
ax4.set_ylabel('Circumferential Angle (deg)', fontsize=11)
ax4.set_title('von Mises Stress (MPa)', fontsize=12, fontweight='bold')
cbar4 = plt.colorbar(contour4, ax=ax4)
cbar4.set_label('Stress (MPa)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case5_exhaust_system.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例5完成，图片已保存")

delta_T_cycle = T_exhaust_peak - T_exhaust_steady
C_fatigue = 5e12
m_fatigue = 3
N_fatigue = C_fatigue / (delta_T_cycle**m_fatigue)
print(f"估算热疲劳寿命: {N_fatigue:.0f} 循环")

# ============================================================
# 案例6：深海管道热屈曲分析
# ============================================================
print("\n" + "="*60)
print("案例6：深海管道热屈曲分析")
print("="*60)

D_o = 0.3239
t_pipe = 0.0127
L_pipe = 1000
H_buried = 1.5

D_i = D_o - 2 * t_pipe
r_o = D_o / 2
r_i = D_i / 2

E_pipe = 207e9
sigma_y = 450e6
alpha_pipe = 12e-6
k_pipe = 45
nu_pipe = 0.3

T_inlet = 120
T_seawater = 4
p_oper = 15e6
mu_soil = 0.6
q_soil = 50e3

x_pipe = np.linspace(0, L_pipe, 200)
h_conv = 100
U_overall = 1 / (1/h_conv + t_pipe/k_pipe)

T_pipe = T_seawater + (T_inlet - T_seawater) * np.exp(-x_pipe / 5000)

delta_T_pipe = T_pipe - T_seawater
sigma_thermal = E_pipe * alpha_pipe * delta_T_pipe

A_cross = np.pi * (r_o**2 - r_i**2)
I_moment = np.pi * (r_o**4 - r_i**4) / 4
EA = E_pipe * A_cross

k_soil = 100e6
P_cr_lateral = 2 * np.sqrt(k_soil * E_pipe * I_moment)
P_cr_upheaval = q_soil * L_pipe / 10
nu_effect = nu_pipe * p_oper * np.pi * r_i**2 / A_cross
P_actual = sigma_thermal * A_cross + nu_effect * A_cross

SF_lateral = P_cr_lateral / np.max(P_actual)
SF_upheaval = P_cr_upheaval / np.max(P_actual)

print(f"管道入口温度: {T_inlet:.1f}°C")
print(f"管道出口温度: {T_pipe[-1]:.1f}°C")
print(f"最大热应力: {np.max(sigma_thermal)/1e6:.2f} MPa")
print(f"侧向屈曲临界载荷: {P_cr_lateral/1e6:.2f} MN")
print(f"上浮屈曲临界载荷: {P_cr_upheaval/1e6:.2f} MN")
print(f"最大实际轴向力: {np.max(P_actual)/1e6:.2f} MN")
print(f"侧向屈曲安全系数: {SF_lateral:.2f}")
print(f"上浮屈曲安全系数: {SF_upheaval:.2f}")

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

ax1 = axes[0, 0]
seabed = Rectangle((-50, -2), L_pipe + 100, 2, fill=True, color='sandybrown', alpha=0.7)
ax1.add_patch(seabed)
pipe = Rectangle((0, 0), L_pipe, D_o*10, fill=True, color='steelblue', alpha=0.8)
ax1.add_patch(pipe)
water = Rectangle((-50, D_o*10), L_pipe + 100, 50, fill=True, color='lightblue', alpha=0.3)
ax1.add_patch(water)
ax1.set_xlim(-50, L_pipe + 50)
ax1.set_ylim(-5, 60)
ax1.set_xlabel('Pipeline Length (m)', fontsize=11)
ax1.set_ylabel('Elevation (m)', fontsize=11)
ax1.set_title('Subsea Pipeline System', fontsize=12, fontweight='bold')
ax1.text(L_pipe/2, 30, 'Seawater', ha='center', fontsize=10)
ax1.text(L_pipe/2, -1, 'Seabed', ha='center', fontsize=10)

ax2 = axes[0, 1]
ax2.plot(x_pipe, T_pipe, 'r-', linewidth=2, label='Pipe Temperature')
ax2.axhline(y=T_seawater, color='b', linestyle='--', label='Seawater Temperature')
ax2.fill_between(x_pipe, T_seawater, T_pipe, alpha=0.3, color='red')
ax2.set_xlabel('Pipeline Length (m)', fontsize=11)
ax2.set_ylabel('Temperature (°C)', fontsize=11)
ax2.set_title('Temperature Distribution Along Pipeline', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
ax3.plot(x_pipe, sigma_thermal/1e6, 'b-', linewidth=2)
ax3.axhline(y=sigma_y/1e6, color='r', linestyle='--', label='Yield Strength')
ax3.fill_between(x_pipe, 0, sigma_thermal/1e6, alpha=0.3, color='blue')
ax3.set_xlabel('Pipeline Length (m)', fontsize=11)
ax3.set_ylabel('Thermal Stress (MPa)', fontsize=11)
ax3.set_title('Thermal Stress Distribution', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
ax4.plot(x_pipe, P_actual/1e6, 'g-', linewidth=2, label='Actual Axial Force')
ax4.axhline(y=P_cr_lateral/1e6, color='r', linestyle='--', label='Lateral Buckling Limit')
ax4.axhline(y=P_cr_upheaval/1e6, color='m', linestyle='--', label='Upheaval Buckling Limit')
ax4.fill_between(x_pipe, 0, P_actual/1e6, alpha=0.3, color='green')
ax4.set_xlabel('Pipeline Length (m)', fontsize=11)
ax4.set_ylabel('Axial Force (MN)', fontsize=11)
ax4.set_title('Buckling Analysis', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case6_subsea_pipeline.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例6完成，图片已保存")

print(f"\n屈曲风险评估:")
if SF_lateral > 2.0 and SF_upheaval > 2.0:
    print("✓ 管道设计安全，无屈曲风险")
elif SF_lateral > 1.5 or SF_upheaval > 1.5:
    print("⚠ 管道存在潜在屈曲风险，建议增加约束")
else:
    print("✗ 管道屈曲风险高，需要重新设计")

print("\n" + "="*60)
print("所有案例仿真完成！")
print("="*60)

print("\n生成的图片文件:")
print("  - case1_heat_exchanger.png")
print("  - case2_turbine_blade.png")
print("  - case3_fuel_rod.png")
print("  - case4_electronics.png")
print("  - case5_exhaust_system.png")
print("  - case6_subsea_pipeline.png")