热冲击分析 — 急速温度变化引起的裂纹发生与比奥数的作用

好问题。比奥数是物体内部热传导阻力与表面热传递阻力的比值。公式为

这在很大程度上取决于材料的脆性。例如，陶瓷（Al₂O₃，k=30 W/(m·K)）这种脆性材料，在比奥数超过1的条件下，比如高温气体（h=1000 W/(m²·K)）作用于薄壁部件（Lc=0.01m）时，

本质区别在于时间导数项的有无。热冲击问题完全由时间变化主导，需使用非定常（过渡）傅立叶方程：

需要使用考虑热应变的本构关系（广义胡克定律）。对于各向同性材料，热应力σ_th与以下关系相关：

为了捕捉热冲击的急速变化，必须采用无条件稳定的隐式法（如后向差分法）。显式法对时间步长Δt有严格的稳定性条件（CFL条件），含有微小单元的网格会导致计算成本膨胀。隐式法在每个时间步需求解线性方程组

热冲击初期（温度急变瞬间）应采用自动步长控制（Automatic Time Stepping）。例如Ansys中用"DELTIM"命令设定初始、最小、最大步长。经验上，热传播过特征长度Lc需要的时间的1/10以下应作为最小步长。若热扩散率为α=k/(ρc_p)，则

理想情况是网格一致，但不是强制要求。多数求解器都有在不同网格间映射（插值）结果的功能。不过精度损失确实可能发生。特别是在热冲击产生陡峭温度梯度的区域，若结构网格过粗，会低估热应变。对策是在温度梯度大的地方（表面和裂纹尖端）充分细化两种网格，并提高单元阶数。Abaqus中可用"*TIE"或"*COUPLING"控制映射。评估映射误差的方法是，对比热分析网格节点温度与映射后结构网格节点温度，运行验证分析。

实务上以下6步较为可行。
1. **前提梳理**：从规格书或试验确定最大温度差ΔT_max、热传递系数h、曝露时间t_exp。计算比奥数以掌握问题规模。
2. **几何和网格**：温度梯度大的表面层需特别细化网格。边界层网格有效。厚度方向至少5个单元以上。
3. **材料设置**：尽可能输入温度依存性（k, c_p, α, E, ν）。若缺乏高温侧材料数据，参考JAHM（高温材料数据库）等。
4. **边界条件**：重现热冲击时，表面热传递系数用阶跃函数或陡峭函数给定。考虑对流和辐射的耦合。
5. **求解执行**：过渡热传导→静结构的顺序进行单向耦合分析。时间步长自动控制。
6. **后处理**：提取最大热应力发生时刻和位置、温度履历、热应力随时间变化等。

至少要确认以下4点。
**1. 能量守恒**：求解器输出的热能守恒（输入热量 vs 内能增加+散失热量）误差在1%以内否。Ansys用"PRENERGY"命令确认。
**2. 网格依赖性**：最细网格（基准网格）结果与粗一级网格结果的最大热应力相差5%以内否。特别要改变表面网格尺寸验证。
**3. 时间步长依赖性**：最小时间步长减半后，温度、应力峰值是否不变。
**4. 简化理论值对比**：平板表面与中心的温度差ΔT_s-c可用简化式（恒热流近似）概算

裂纹由拉应力引起，故第一主应力（σ1）最为关键。评估步骤如下：
1. **最大拉应力的位置和时刻**：过渡分析中，σ1是否超过材料拉伸强度（如陶瓷按JIS R 1606四点弯曲强度）。超过则裂纹风险高。
2. **应力梯度**：表面拉、内部压的分布。拉应力区域的深度是裂纹侵入深度的参考。
3. **热应力集中系数**：几何棱角或孔洞边缘处，理论热应力乘以形状系数（Kt）。用FEM最大应力除以平板理论公称热应力，求得实际Kt。
4. **断裂力学参数**：假设初始缺陷，用J积分等评估其位置的应力强度因子KI，对比材料断裂韧性KIC（ASTM E399规范等）。

取决于用途。**Ansys Mechanical**工作流完善，特别是"Steady-State Thermal"→"Transient Thermal"→"Static Structural"的单向耦合操作简便。用APDL脚本可实现高度控制。**Abaqus/Standard**在*COUPLED TEMPERATURE-DISPLACEMENT分析步中强势，能做完全的热-结构双向耦合分析（考虑热对应力的反馈）。若问题涉及塑性发热或摩擦热，Abaqus更有优势。**COMSOL Multiphysics**的耦合设置用GUI直观，可把热传导和固体力学作为一个"物理场"同步求解。原型设计和参数扫描较快。

方法基本相同，但输入格式有差。
- **Ansys**：Engineering Data中以温度为独立变量，表格式输入。k、c_p、α、E、ν等单独定义。APDL中用"MPTEMP"和"MPDATA"命令。
- **Abaqus**：材料属性中用"Depvar"，或编用户子程序"USDFLD"或"UMAT"，支持更复杂的依存关系（如退火效应）。
- **COMSOL**：材料库丰富，许多材料的温度依存数据预装。自定义函数（解析函数、插值函数）的GUI设置简单。
通用技巧是：给定100°C间隔的数据点5个以上，精度会显著提升。特别陶瓷的比热c_p温度依存性大，需注意。

可以，但前处理和配置耗时较多。**Code_Aster**（Salome-Meca环境）功能强大，支持过渡热传导（THER_NON_LINE）和热应力（MECA_STATIQUE）的耦合分析。可处理材料温度依存和非线性。问题在于所有配置都需写在命令文件（.comm）里，学习曲线陡峭。**CalculiX**亦然。商用软件的GUI负责做网格质量检查和求解器稳定化的默认设置，用户需自觉意识到这些细节。小规模检讨和教学用可，但生产设计要快速得可信结果，商用软件的完成度和技术支持仍有优势。

主要有3个原因。
1. **陡峭的边界条件**：热传递系数h用阶跃函数给定，初期温度变化率接近无穷，数值不稳定。对策：给现实的上升时间（如1毫秒），用线性或指数函数。Ansys中用"Function"工具。
2. **材料特性不连续**：温度依存数据表中相邻温度点材料值差异过大（如100°C处k=50、101°C处k=10）。数据平滑化或用更细的温度间隔。
3. **不当的时间步长**：隐式法也会在极大步长下产生非物理振荡。把自动步控的"最大步长"限在热时常（Lc²/α）量级以下。

很可能是**约束条件**或**观察时刻**出错。首先检查热膨胀变形能否自由进行。若物体全部固定，温度上升区域会全体压缩。热冲击中，冷的内部充当"拘束材"，压制加热表面的膨胀，故表面被拉伸。要重现此机制，需正确施加对称边界条件，仅约束刚体运动。其次，确认看的是过渡分析的哪个时刻。热冲击直后，表面受热膨胀，一时可能出现表面压应力。热传向内部后应力平衡改变。必须用动画查看应力随时间的变化过程。

这是存在**几何奇点**（尖锐的角、裂纹尖端）时的典型症状。线性弹性分析中，若曲率半径趋零，理论应力趋无穷。FEM网格越细，越接近无穷大，不收敛。对策两个：
1. **修正为现实的几何形状**：设计上不存在完全尖角。引入制造倒角（如R0.1mm），应力集中因子就落到有限值。
2. **采用断裂力学方法**：若问题前提是裂纹已存在，则用应力强度因子等参数评估。