半无限固体的非定常热传导

半无限固体的非定常热传导的理论基础

概要和适用场景

• 焊接入热的初期阶段：电弧焊或激光焊时，母材受到的热量在初期扩散。HAZ（热影响区）的温度分布可以用半无限体解来估算
• 激光加工的表面温度：激光照射瞬间到数毫秒间的表面温度上升。这与加工深度的制御直接相关
• 地中的冻结深度：冬季地表温度骤降时，冻结能进展到多深。水管埋设深度的设计依据
• 铸造锻造的快速冷却：高温金属用水或油急冷时，表面附近的冷却速度。判断是否发生马氏体变态的关键

支配方程式和相似变量

这里 $\alpha = k/(\rho c_p)$ 是热扩散率 [m²/s]。钢铁约 $1.2 \times 10^{-5}$，铝约 $9.7 \times 10^{-5}$ m²/s。这个值越大，温度变化传播得越快。

• 初始条件：$T(x, 0) = T_i$（全体均匀温度）
• 表面边界条件：$T(0, t) = T_s$（$t > 0$ 时表面温度变为常数）
• 远方条件：$T(\infty, t) = T_i$（充分远处保持初始温度）

还有一个关键技巧是引入相似变量：

这个变量变换可以把偏微分方程转化为常微分方程。$x$ 和 $t$ 两个独立变量会并入 $\eta$ 一个变量。这叫做波尔兹曼变换。

分析解的导出

边界条件是 $\theta(0) = 0$, $\theta(\infty) = 1$。这个方程用误差函数 erf 可以求解：

代入原温度后，得到半无限固体温度分布的最重要的公式：

其中 $\operatorname{erfc}(\eta) = 1 - \operatorname{erf}(\eta)$ 是互补误差函数。

也就是说 $\eta = 2$（$x = 4\sqrt{\alpha t}$）以深处基本没有温度变化。这就是浸透深度的理论根据。

浸透深度和表面热流束

比如钢铁（$\alpha \approx 1.2 \times 10^{-5}$ m²/s）在焊接中瞬间受热的情况：

• $t = 1$ 秒后：$\delta \approx 4\sqrt{1.2 \times 10^{-5}} \approx 14$ mm
• $t = 10$ 秒后：$\delta \approx 44$ mm
• $t = 100$ 秒后：$\delta \approx 139$ mm

如果板的厚度远大于这个 $\delta$，就可以放心用半无限体近似。比如20 mm厚的钢板，焊接后 $t \lesssim 2$ 秒左右，半无限体近似都是有效的。

它与 $1/\sqrt{t}$ 成正比，所以初期需要剧烈的热流，之后逐渐减弱。当 $t \to 0$ 时，$q \to \infty$ —— 这对应表面的温度梯度无限陡峭。实际的热源出力有限，所以非常初期用一定温度条件不如一定热流条件更符合实际。

边界条件的3种模式

Case 2 的表面温度为 $T(0,t) = T_i + \frac{2q_0''}{k}\sqrt{\frac{\alpha t}{\pi}}$，与 $\sqrt{t}$ 成正比而上升。激光加工通常出力恒定，所以 Case 2 更接近实际。

半无限体近似的妥当性判定

• $\mathrm{Fo} < 0.05$：半无限体近似的误差不足1%。可以放心使用
• $0.05 < \mathrm{Fo} < 0.2$：裏面的影响开始出现。要注意精度
• $\mathrm{Fo} > 0.2$：已经不能算半无限体了。应该用有限厚度的傅里叶级数解

比如厚度50 mm的钢板（$\alpha = 1.2 \times 10^{-5}$）中 $\mathrm{Fo} = 0.05$ 对应的时间是 $t = 0.05 \times 0.05^2 / (1.2 \times 10^{-5}) \approx 10.4$ 秒。也就是焊接后10秒以内，用半无限体来计算没有问题。

半无限固体的非定常热传导的数值计算方法

分析解无法使用的场景

• 温度依赖物性：钢的热传导率从0°C的约60 W/(m·K)降到800°C的约30 W/(m·K)。变成非线性就没有分析解了
• 二维三维效应：如焊接焊道宽度有限，需要考虑横向热扩散
• 随时间变化的边界条件：如焊接枪的移动、激光的脉冲照射
• 相变：融熔池形成和凝固（Stefan 问题）
• 复合材料异种材料接合：界面热阻存在的情况

用有限元法离散化半无限固体

1. 确保足够大的分析区域：浸透深度 $\delta = 4\sqrt{\alpha t_{\max}}$ 的2〜3倍的区域，在远方边界设置 $T = T_i$（Dirichlet条件）。最简单可靠
2. 使用无限单元：如Ansys的INFIN111或Abaqus的CINPE4等，内置远方温度衰减的特殊单元在边界配置。能大幅缩小分析区域

热传导的FEM表述中，Galerkin法得到以下离散系统：

其中 $[C]$ 是热容量矩阵，$[K]$ 是热传导矩阵，$\{F\}$ 是热负荷向量。这与结构分析的 $[M]\{\ddot{u}\} + [K]\{u\} = \{F\}$ 形式相似，区别只是一阶还是二阶时间导数。

网格策略

在Ansys中设置Bias Factor 为 10〜20，就能自动得到不错的非均匀网格。Abaqus用 edge seed 的 single bias 功能。

时间积分方案

• 初期时间步长：$\Delta t_1 = \Delta x_{\min}^2 / (6\alpha)$ 左右（表面最小单元的热扩散时间的1/6）
• 时间步长的增长：每个步长乘以1.2〜1.5倍（几何级数增大）
• 自动时间步长制御：交给商用求解器的自动制御最安全。启用Abaqus的自动时间步长或Ansys的AUTO TIME STEPPING

陰解法（backward Euler、Crank-Nicolson）没有稳定性约束，但为了精度初期还是需要细小的步长。

无限单元的活用

用无限单元能把分析区域缩到浸透深度的1.5倍，大幅削减单元数。但无限单元的方向配置错误会导致无意义的结果，必须用理论解来验证。

半无限固体的非定常热传导的实际应用

焊接HAZ的温度预测

线热源入热量 $Q$ [J/m] 在厚板表面移动的情况，从罗森塔尔移动点热源解推导出的冷却速度近似式：

简化版是，对于板厚 $d$ > 浸透深度的情况：

比如入热 $Q = 1.5$ kJ/mm、无预热（$T_0 = 20$°C）、板厚25 mm的软钢，$t_{8/5} \approx 12$ 秒。小于6秒会有氢脆风险，大于30秒韧性下降。这直接关系焊接施工条件的设计。

激光加工的表面温度

例如，SUS304（$k = 16$ W/(m·K)、$\alpha = 4.2 \times 10^{-6}$ m²/s）被功率密度 $I_0 = 10^8$ W/m² 的激光照射 1 ms 时：

$T(0, 1\text{ms}) = 20 + \frac{2 \times 10^8}{16}\sqrt{\frac{4.2 \times 10^{-6} \times 10^{-3}}{\pi}} \approx 20 + 1,450 = 1,470$°C

超过熔点1400°C，开始融解。调控照射时间就能制御加工深度。但一旦开始融解，潜热效应会减缓实际温度上升，精密预测需要相变模型（Stefan问题）。

地中的冻结深度推测

$\operatorname{erfc}(\eta_f) = 5/20 = 0.25$ → $\eta_f \approx 0.81$

土壤的热扩散率约 $\alpha \approx 5 \times 10^{-7}$ m²/s，经过一周（$t = 6 \times 10^5$ 秒）：

$x_f = 2 \times 0.81 \times \sqrt{5 \times 10^{-7} \times 6 \times 10^5} \approx 0.89$ m

也就是约90 cm处被冻结。防冻需要埋设在这个深度以下。日本寒冷地带的水道管埋设深度1.0〜1.2 m 正是基于这个计算。但实际土壤含水冻结伴随潜热放出，冻结进行会比 erfc 解慢。更精密的计算用诺伊曼解（Stefan问题的分析解）。

接触温度的计算

这里 $e = \sqrt{k \rho c_p}$ 是热浸透率（thermal effusivity），单位J/(m²·K·s^{1/2})。

比如200°C的钢（$e \approx 14,000$）和20°C的铜（$e \approx 37,000$）接触：

$T_c = (14000 \times 200 + 37000 \times 20)/(14000 + 37000) \approx 69$°C

接触面温度偏向铜一侧——热浸透率大的物体有更强的"推送温度"的力量。冬天不锈钢的扶手摸起来冷，木制扶手温暖，原因就在于扶手材料的热浸透率不同。CAE中的接触热问题也用同样的原理。

半无限固体的非定常热传导的软件对比

商业工具中的设置方法

1. 用 Case 1（一定温度边界）的 FEM 结果和 erfc 解比较
2. 确认表面温度的时间历史和深度方向的温度分布都一致
3. 表面热流束 $q_s'' = k(T_s - T_i)/\sqrt{\pi\alpha t}$ 的比较（检验数值微分精度）
4. 网格收敛性：3层以上网格密度的解是否收敛

NAFEMS 的 T2 基准问题（一维非定常热传导）正是这个。

开源代码中的实现

分析解例（Python）：

import numpy as np
from scipy.special import erfc

alpha = 1.2e-5   # 钢的热扩散率 [m²/s]
Ti, Ts = 20, 800  # 初始温度, 表面温度 [°C]
t = 5.0           # 时间 [s]
x = np.linspace(0, 0.05, 100)  # 深度 [m]

T = Ti + (Ts - Ti) * erfc(x / (2 * np.sqrt(alpha * t)))
delta = 4 * np.sqrt(alpha * t)  # 浸透深度
print(f"浸透深度: {delta*1000:.1f} mm")

用 FEniCS 解同样的问题，一维网格 + Dirichlet BC + backward Euler，几十行代码就行。可以用来检验商业求解器的结果。

半无限固体的非定常热传导的前沿研究

非傅里叶热传导

这种情况需要Cattaneo-Vernotte 方程（双曲型热传导方程）：

$\tau_q$ 是热缓和时间（金属约 $10^{-11}$〜$10^{-13}$ 秒）。这个方程是波动方程和扩散方程的混合体，出现有限速度的"热波"传播。飞秒激光加工、半导体热点分析、极低温物理中很重要。

逆问题的应用

• 铸造的界面热传达系数同定：从铸型内的温度历史反算熔液-铸型间的传热系数
• 火箭发动机热负荷推测：从壁面温度计测推算表面热流
• 急冷焙烧的冷却曲线分析：从内部温度推测表面冷却速度

因为有分析解，逆问题的核函数（格林函数）已知，配合正则化法（Tikhonov正则化等）就能得到稳定的逆分析。结合机械学习（PINN: Physics-Informed Neural Network）的最新方法研究也在进行。

半无限固体的非定常热传导的半无限固体的非定常热传导常见问题与调试

常见的失败和对策

验证用基准问题

1. NAFEMS T2：一维非定常热传导。一定温度边界，32秒后深度0.01 m处的温度验证。正解用erfc解计算
2. Case 1自作基准：钢（$k=50$, $\rho=7850$, $c_p=500$）中 $T_i=20$, $T_s=500$°C。$t=10$ s后的温度分布与erfc解比较。$x=5$ mm处正解约 $T \approx 434$°C
3. Case 2（一定热流）：同物性下 $q_0''=10^6$ W/m²。表面温度的时间历史与分析解比较
4. 表面热流验证：Case 1 的 FEM 表面热流输出与 $k(T_s-T_i)/\sqrt{\pi\alpha t}$ 比较。数值微分的精度检查

所有情况 2%以内的误差就说明网格和时间步长足够了。