瞬态导热的数值解法
理论：θ法的三种选择
概述
求解瞬态导热时，该选显式法还是隐式法？
各自弱点不同。显式法每步计算量小，但时间步长有严格上限，超过就迅速发散。后向欧拉法无条件稳定但只有一阶精度，大步长时误差大。Crank-Nicolson法为二阶精度很准确，但遇到突变时解会振荡，可能给出物理上不可能的温度。我们用同一个平板急冷问题把三者并列比较。
θ法公式
θ = 0 为显式法，θ = 1 为后向欧拉法，θ = 1/2 为Crank-Nicolson法。
稳定条件与放大因子
显式法需要r ≤ 1/2。Crank-Nicolson法步长大时，短波长分量的放大因子接近−1，即突变分量几乎不衰减，每步变号残留下来，这就是振荡的根源。
克兰克与尼科尔森
Crank-Nicolson法由英国的约翰·克兰克和菲莉丝·尼科尔森于1947年发表，作为数值求解热传导型方程的实用方法。那时电子计算机尚未普及，这是在减少手摇计算机工作量的同时提高精度的办法。几乎同时，显式法的不稳定性开始用冯·诺伊曼提出的稳定性分析来解释。至今这两种思路仍是许多热分析软件时间积分的基础。
算例：100℃平板急冷
厚1，热扩散率1，初始100℃，两面保持0℃，20等分（Δx = 0.05），比较t = 0.1时的中心温度。精确值为47.449℃。
显式法
| r | 步数 | 中心温度 |
|---|---|---|
| 0.40 | 100 | 47.22 °C |
| 0.50 | 80 | 46.87 °C |
| 0.52 | 77 | 164.6 ℃（发散） |
后向欧拉法与Crank-Nicolson法
| r | 步数 | 后向欧拉：中心（误差） | CN：中心（误差） | CN：第一步后最低温度 |
|---|---|---|---|---|
| 0.4 | 100 | 47.67 (+0.22) | 47.45 (−0.003) | 70.8 °C |
| 2 | 20 | 48.55 (+1.10) | 47.44 (−0.011) | 23.6 °C |
| 10 | 4 | 52.28 (+4.83) | 47.14 (−0.31) | −28.4 °C |
Crank-Nicolson法中心温度几乎正确，表面附近却出现了低于0℃的温度。
急冷初期，表面附近温度急剧变化。这种突变由短波长分量构成，Crank-Nicolson法在大步长时几乎不衰减它们，只是变号残留。所以第一步就出现−28℃，这在两面0℃、初始100℃的问题中绝不可能出现。像中心这样变化平滑的地方很准确，只看平均值就察觉不到。在物性随温度变化或凝固这类随温度改变状态的问题中，这种不可能的值会把计算带偏。后向欧拉法强烈衰减这些分量，不会振荡，但整体被抹平，误差较大。实际中常用前几步后向欧拉法、再切换到Crank-Nicolson法的做法。
选择与注意
- 显式法遵守r ≤ 1/2（二维1/4、三维1/6），由最小单元决定。
- 隐式法也要确认改变时间步长结果不变。
- 有突变的问题用Crank-Nicolson法时，开头用后向欧拉法。
- 检查是否出现物理上不可能的值（超出初值和边界值范围）。
- 非线性问题可考虑自动调整时间步长。
“不该凝固的地方凝固了”
某铸造热分析中，熔融金属刚注入型腔，与型壁交界的几个点瞬间低于固相线温度，由此扩展出奇怪的凝固。原因是大时间步长下使用Crank-Nicolson法产生的振荡：交界处温度突变使温度暂时向下摆动，凝固判定对此作出了反应。把前几步改为后向欧拉法后，不自然的凝固消失了。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 显式法r > 1/2 | 发散 | 遵守稳定条件 |
| 隐式法步长过大 | 误差大 | 改变步长确认 |
| 突变问题直接用CN | 不可能的温度 | 开头用后向欧拉 |
| 只检查平均值 | 漏掉振荡 | 也检查最小最大值 |
| 忽略最小单元 | 局部发散 | 按最小尺寸定步长 |
我想了解相关内容。