半无限固体非稳态热传导

当然，现实中并不存在“无限厚的物体”。但是，当热量突然施加到厚壁或地面时，在热量到达背面之前的初始阶段，背面“等同于不存在”，因此可以将其视为半无限体来处理。

举一些实际工作中常见的例子：

• 焊接刚输入热量后：电弧焊或激光焊接时，热量进入母材并扩散的初始阶段。HAZ（热影响区）的温度分布正可以用半无限体的解来估算。
• 激光加工的表面温度：从激光照射瞬间开始数毫秒内的表面温升。这直接关系到加工深度的控制。
• 地中的冻结深度：冬季地表温度急剧下降时，估算冻结会进展到多深。这成为水管埋设深度的设计依据。
• 铸造/锻造的急冷：高温金属用水或油淬火时，表面附近的冷却速度。用于判断是否会发生马氏体相变。

半无限固体（$x \geq 0$）的一维非稳态热传导方程如下：

这里 $\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$ 这一个变量。这被称为玻尔兹曼变换。

意味着“位置 $x$ 远一倍的点，其温度变化在时间 $t$ 变为4倍时是相同的”。如果 $\eta$ 相同，则温度分布的形状相同——这称为自相似性。可以想象成，随着时间的推移，温度分布按 $\sqrt{t}$ 的比例向深处“延伸”。

将 $\theta(\eta) = (T - T_s)/(T_i - T_s)$ 无量纲化后，控制方程变为：

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

代回原始温度，就得到了半无限固体温度分布最重要的公式：

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

记住以下这些值会比较方便：

也就是说，在 $\eta = 2$（即 $x = 4\sqrt{\alpha t}$）以深的地方，几乎没有温度变化。这成为渗透深度的依据。

正是如此。热渗透深度（thermal penetration depth）是温度变化几乎到达的深度的参考值：

例如，钢（$\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$ 秒左右的时间内，半无限体近似是有效的。

表面（$x = 0$）的热流通过对 erfc 解对 $x$ 求导得到：

它与 $1/\sqrt{t}$ 成正比，因此初期需要极大的热流，并随时间衰减。$t \to 0$ 时 $q \to \infty$，这对应于初始时刻表面温度梯度为无穷大。实际上，任何热源都有有限的输出功率，所以在极初期，恒定热流条件比恒定温度条件更符合现实。

半无限固体典型的边界条件有三种模式：

Case 2 的表面温度为 $T(0,t) = T_i + \frac{2q_0''}{k}\sqrt{\frac{\alpha t}{\pi}}$，与 $\sqrt{t}$ 成正比上升。在激光加工中，由于输出功率恒定，Case 2 往往更接近实际情况。

可以用傅里叶数 $\mathrm{Fo} = \alpha t / L^2$ 来判定（$L$ 是板厚）。粗略地说：

• $\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\}$ 形式很像吧？区别在于是一阶时间导数还是二阶。

绝对不能用均匀网格。erfc 的温度分布在表面附近梯度很陡，越往深处越平坦。因此，必须采用从表面开始按指数规律增大尺寸的不等分割网格（偏置网格）。

在 Ansys 中，将 Bias Factor 设置为 10~20 左右，会自动生成合适的不等分割网格。在 Abaqus 中，使用 edge seed 的 single bias 功能。

很敏锐的观察。在 $t = 0$ 附近，表面热流按 $1/\sqrt{t}$ 发散，因此初始时间步长也必须设置得非常小。具体来说：