MATLAB求解瞬态热传导方程：四种数值格式详解与实现

说到热传导方程，很多用过商业仿真软件的同学第一反应都是“这不是ANSYS或者COMSOL里点几下就完事了吗，干嘛还要用MATLAB自己折腾？”这个问题我被问过很多次。商业软件确实方便，但等你真遇到非标边界条件、突然想加一个体积热源、或者想把导热和流动耦合起来算的时候，你就会发现，不懂底层的数值格式，你连调试都不知道从哪下手。用MATLAB从头写一个求解器，目的从来不是为了跟商业软件比规模，而是把“温度到底是怎么一步步算出来”这件事彻底弄明白。

今天这算例的定位很明确：用一维和二维的瞬态导热问题当抓手，把显式、隐式、Crank-Nicolson和ADI这几种常用计算格式挨个过一遍。一维问题帮你建立“离散化+时间推进”的基本概念，二维问题帮你理解为什么直接照搬一维思路会出事、以及ADI格式是怎么救场的。适合正在学数值计算的学生、刚接触热仿真的工程师，以及所有想在MATLAB里快速验证自己导热模型的同学。看完这篇文章，你至少能独立搭建一个一维和一个二维的瞬态导热求解器，而且会清楚地知道每种格式的脾气和踩坑点。

1. 控制方程与计算格式的取舍逻辑

1.1 先从傅里叶定律说起

热传导的物理本质是热量从高温区域向低温区域传递，这个过程的宏观规律由傅里叶定律描述：热流密度正比于温度梯度，比例系数就是导热系数k。把傅里叶定律和能量守恒结合，就能得到我们熟悉的瞬态热传导方程：

[
\frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} \right)
]

其中 (\alpha = k/(\rho c_p)) 称为热扩散系数，单位是 m²/s。这个参数决定了温度扰动在材料内部传播的快慢。钢材的α大约在 (1 \times 10^{-5}) 到 (2 \times 10^{-5}) m²/s 之间，铝的α大约是钢的5到10倍，所以铝材温度均匀化速度远快于钢材，这个差异在仿真里非常关键。

方程本身是一类经典的抛物线型偏微分方程（PDE），它的特点很鲜明：信息以无限大速度衰减传播，也就是说某一点的温度变化几乎会同时影响全场，只是影响幅度随距离衰减。这个特性决定了数值格式的选择和边界条件的处理方式。

1.2 离散化是数值求解的必经之路

解析解只在极少数简化条件下存在，比如无限大平面、恒定边界条件、规则几何。工程上绝大多数问题都要靠数值方法，核心思想就一句话：把连续的求解域离散成有限个网格点，把偏导数替换成差分近似，把PDE转化成一堆代数方程，然后用计算机迭代求解。

在MATLAB里有两条路可以走：一是直接调用PDE工具箱，二是自己写差分格式。我的建议是，学习阶段一定要自己手写一遍差分格式。原因很简单：PDE工具箱封装了太多内部细节，用完你会觉得“哦跑出来了”，但对差分格式、稳定性、边界条件植入这些核心概念还是云里雾里。自己写代码虽然繁琐，但每一步都能看到温度场的演变，遇到问题也知道改哪里。

空间域的离散通常统一用二阶中心差分：(\partial^2 T/\partial x^2 \approx (T_{i+1} - 2T_i + T_{i-1})/\Delta x^2)。这个近似精度是二阶，也就是网格加密一倍，误差缩小到四分之一。时间域的离散可以有不同的玩法，这就是我们要讨论的“计算格式”问题。

1.3 三种时间推进格式的性格对比

把空间离散后的方程写成半离散形式后，时间推进上有三种常见格式：

显式格式（前向欧拉）：下一时刻的温度直接用当前时刻的相邻点温度计算。公式是
[
T_i^{n+1} = T_i^n + F_o (T_{i+1}^n - 2T_i^n + T_{i-1}^n)
]
其中 (F_o = \alpha \Delta t / \Delta x^2) 是傅里叶数。这个格式编程极其简单，一个for循环就能搞定，但代价是要受严格稳定性限制。一维情况下要求 (F_o \le 0.5)，二维情况下更苛刻：(F_o \le 0.25)。一旦时间步长超过这个限制，温度场就会开始抖动、锯齿化，然后直接发散。

隐式格式（后向欧拉）：空间二阶导数取下一时刻的值，也就是
[
T_i^{n+1} - F_o(T_{i+1}^{n+1} - 2T_i^{n+1} + T_{i-1}^{n+1}) = T_i^n
]
所有未知量耦合在一起，需要联立求解三对角方程组。每步计算量比显式大，但好处是无条件稳定，想用多大的时间步长都行。不过这里有个陷阱：稳定性不代表准确性，步长过大的时候结果同样会偏离真实解，只是不会发散而已。

Crank-Nicolson格式：显式和隐式各取一半，像一个双方各让一步的调解方案。时间上的截断误差是二阶，精度比前两个都高，同时无条件稳定。代价也是要解三对角方程组，复杂性跟隐式差不多，但计算精度收益很高。因此工程模拟中最常推荐的其实就是Crank-Nicolson。

2. 一维算例：金属杆的温度驰豫过程

2.1 问题设定与物理意义

先做一个最简单的一维算例：一根长度 (L=1) m 的金属杆，初始温度均匀为100°C。从 (t=0) 开始，左端温度突变为0°C并保持，右端保持绝热（即温度梯度为零），观察杆内温度如何随时间变化。

这个设定在工程里很有代表性——想象一个金属工件一端接触冷水进行淬火，另一端包着保温材料。热扩散系数 (\alpha = 10^{-4}) m²/s，大概介于钢和铝之间。

空间网格取 (N_x = 100)，则 (\Delta x = 0.01) m。时间步长设为多少，取决于用哪种格式。

2.2 显式格式实现与稳定性管理

取傅里叶数 (F_o = 0.25)，满足稳定性要求 (F_o \le 0.5)。于是时间步长：
[
\Delta t = F_o \frac{\Delta x^2}{\alpha} = 0.25 \times \frac{0.0001}{10^{-4}} = 0.25 \text{ s}
]

跑1000步，对应模拟250秒的物理时间。MATLAB实现如下：

matlab复制% 参数设置
L = 1;                 % 杆长，单位m
Nx = 100;              % 空间网格数
dx = L / Nx;           % 空间步长
alpha = 1e-4;          % 热扩散系数，m^2/s
Fo = 0.25;             % 傅里叶数
dt = Fo * dx^2 / alpha;  % 时间步长 = 0.25s
Nt = 1000;             % 时间步数

% 初始条件
T = 100 * ones(Nx+1, 1);
T(1) = 0;              % 左端固定为0

% 时间推进
for n = 1:Nt
    Tn = T;
    % 内部节点更新
    T(2:end-1) = Tn(2:end-1) + Fo * (Tn(3:end) - 2*Tn(2:end-1) + Tn(1:end-2));
    % 右端绝热边界：T(Nx+1) = T(Nx)  =>  温度梯度为零
    T(end) = T(end-1);
end

% 画图
x = 0:dx:L;
plot(x, T, 'LineWidth', 1.5);
xlabel('位置 (m)');
ylabel('温度 (°C)');
title('一维导热显式格式结果 t=250s');

注意绝热边界的实现方式：令最后两个网格点的温度相等，本质上是边界处的梯度近似为零，这是Neumann边界条件最简单的处理手法，一阶精度就已经够用。

如果 (F_o) 取0.6会怎样？试试就知道，算个几十步温度场就开始锯齿状振荡，然后数值爆炸。这个现象我建议每个读者都亲手试一次，感受一下“看着它发散”的酸爽。

2.3 隐式格式与Crank-Nicolson实现

隐式和Crank-Nicolson都需要解三对角方程组。MATLAB里直接构造稀疏矩阵然后用反斜杠求解最方便，不需要手写Thomas算法：

matlab复制% 隐式格式
Fo = 5;                % 比显式大得多，试试看
dt = Fo * dx^2 / alpha;
Nt = 200;              % 总物理时间约100s

A = zeros(Nx-1, Nx-1);
for i = 1:Nx-1
    A(i,i) = 1 + 2*Fo;
    if i > 1,     A(i,i-1) = -Fo; end
    if i < Nx-1,  A(i,i+1) = -Fo; end
end

% 时间推进
T = 100 * ones(Nx+1, 1);
T(1) = 0;

for n = 1:Nt
    b = T(2:end-1);
    b(1)   = b(1) + Fo * T(1);   % 左端已知温度贡献
    b(end) = b(end);              % 右端绝热，暂无需修改
    T(2:end-1) = A \ b;
    T(end) = T(end-1);
end

Crank-Nicolson的矩阵形式稍微不同，构造矩阵时系数变成 (1 + F_o)，对角线的负系数是 (-F_o/2)。其他部分类似，我就不重复贴全代码了。

用大步长跑出来以后你会发现，结果确实不会发散，但温度曲线在早期可能有一点点“僵硬感”——因为大步长牺牲了时间分辨率。这正是我反复强调的：无条件稳定不等于无条件精确。

2.4 三种格式的对比结论

我自己做项目的习惯是：quick check用显式，正式结果用Crank-Nicolson。隐式用得相对少，因为它的时间一阶精度在长时程模拟里会造成比较明显的数值耗散，温度场会比真实情况“钝化”不少。

3. 二维算例：用ADI格式搞定矩形板的散热

3.1 二维问题为什么不能直接硬来

二维瞬态热传导方程：
[
\frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right)
]

如果照搬一维的全隐式思路，每个时间步需要解一个 ((N_x-1)\times(N_y-1)) 大小的块三对角方程组。假设 (N_x=N_y=100)，那就是一个10000×10000的矩阵，虽然稀疏，但求解代价比一维提升了不止一个量级。显式格式更不用提，稳定性条件变成 (F_o = \alpha\Delta t(1/\Delta x^2 + 1/\Delta y^2) \le 0.5)，网格哪怕稍微加密一点，时间步长就被压得很小，计算量剧增。

ADI格式（Alternating Direction Implicit，交替方向隐式）的巧妙之处在于：把二维问题拆成两个一维问题，每一时间步分两个半步推进。第一个半步在x方向隐式、y方向显式；第二个半步在y方向隐式、x方向显式。这样每个半步只需要解若干个（对每一行或每一列）独立的三对角方程组，计算量大幅下降，同时保持无条件稳定。

3.2 ADI格式的数学推导

取一个标准的Peaceman-Rachford格式。已知第n步的温度场 (T^n)：

第一步（x方向隐式）：
[
\frac{T^{n+1/2} - T^n}{\Delta t/2} = \alpha \left( \frac{\partial^2 T^{n+1/2}}{\partial x^2} + \frac{\partial^2 T^n}{\partial y^2} \right)
]

第二步（y方向隐式）：
[
\frac{T^{n+1} - T^{n+1/2}}{\Delta t/2} = \alpha \left( \frac{\partial^2 T^{n+1/2}}{\partial x^2} + \frac{\partial^2 T^{n+1}}{\partial y^2} \right)
]

离散后每个方向都变成三对角方程组。x方向的半步，对每一行 (j) 求解一个沿x方向的三对角系统；y方向的半步，对每一列 (i) 求解一个沿y方向的三对角系统。实现起来其实就是两层循环。

3.3 MATLAB代码实现

算例设定：一块 (1\text{m} \times 1\text{m}) 的正方形金属板，四周边界固定为0°C。初始时刻整个板内部温度为100°C，然后让板自然冷却。热扩散系数仍取 (10^{-4}) m²/s。网格取 (N_x = N_y = 100)。

matlab复制% 二维ADI求解热传导方程
Lx = 1; Ly = 1;
Nx = 100; Ny = 100;
dx = Lx / Nx;
dy = Ly / Ny;
alpha = 1e-4;

% 用较大的时间步长测试ADI的稳定性
dt = 1;                % 1秒
Nt = 200;              % 总共200秒

% 傅里叶数（用于参考）
Fo_x = alpha * dt / dx^2;
Fo_y = alpha * dt / dy^2;

% 初始条件
T = 100 * ones(Nx+1, Ny+1);
% 边界固定为0（Dirichlet条件）
T([1 end], :) = 0;
T(:, [1 end]) = 0;

% 构造三对角矩阵（x和y方向各自一组）
ex = ones(Nx-1, 1);
Ax = spdiags([-Fo_x*ex/2, (1+Fo_x)*ex, -Fo_x*ex/2], -1:1, Nx-1, Nx-1);

ey = ones(Ny-1, 1);
Ay = spdiags([-Fo_y*ey/2, (1+Fo_y)*ey, -Fo_y*ey/2], -1:1, Ny-1, Ny-1);

% 时间推进
for n = 1:Nt
    % 第一步：x方向隐式，y方向显式
    % 对每个内部行 j = 2:Ny，构造右端项并求解
    for j = 2:Ny
        rhs = T(2:end-1, j) + Fo_y/2 * (T(2:end-1, j+1) - 2*T(2:end-1, j) + T(2:end-1, j-1));
        % 考虑到边界贡献：如果 x方向边界不为0，需要调整rhs
        % 本例中x方向边界为0，因此无需额外调整
        T(2:end-1, j) = Ax \ rhs;
    end

    % 第二步：y方向隐式，x方向显式
    for i = 2:Nx
        rhs = T(i, 2:end-1)' + Fo_x/2 * (T(i+1, 2:end-1) - 2*T(i, 2:end-1) + T(i-1, 2:end-1))';
        % 同理，y方向边界为0，无需调整
        T(i, 2:end-1) = (Ay \ rhs)';
    end
end

% 可视化
[X, Y] = meshgrid(0:dx:Lx, 0:dy:Ly);
surf(X, Y, T');
xlabel('x (m)');
ylabel('y (m)');
zlabel('温度 (°C)');
title('二维热传导ADI结果 t=200s');

这里有个细节：当边界温度不是零时，右端项需要显式地加上边界贡献。上面代码因为边界是0，所以省掉了这部分。实际工程问题中经常遇到非零边界或者对流边界，一定要记得把它们作为源项加到rhs里，否则边界条件就丢了。

还有一点值得注意：上面用了稀疏矩阵 spdiags 而不是全矩阵。当网格数放大到 (200 \times 200) 甚至 (500 \times 500) 时，全矩阵存储的内存开销是灾难性的，用稀疏矩阵能大幅降低内存占用和求解时间。

3.4 结果分析与格式验证

跑出来的温度场是一个从边界向内逐渐冷却的对称分布，中心温度下降得最慢，这是直觉上完全合理的——中心点距离所有边界最远，热量往外传递的路径最长。

验证ADI格式正确性的一个简单办法是拿二维问题的稳态解对比。当 (t \to \infty) 时，板内温度应趋于四周边界的平均值，也就是0°C。算到几百秒后看板中心温度是否趋近于0，可以直观判断程序对错。更严格的方法是收敛性检验：固定物理时间，不断加密网格和时间步长，观察解是否逐渐趋近同一个值。

我实际跑的时候，dt=1 对应的傅里叶数 (F_o = 0.1)（两个方向都是），这个步长如果换显式格式确实也能跑，但代价是每步更新所有内部节点，在二维网格里每步要处理约一万个点的计算。ADI的优势在更大步长下体现更明显：把 dt 加到10，仍然稳定，而同样的步长显式格式早就炸了。

4. 实际调试中的坑与排查技巧

4.1 显式格式发散：先怀疑稳定性

显式格式跑着跑着温度场出现棋盘状锯齿，或者干脆NaN了，90%的情况是傅里叶数超限。不要上来就怀疑程序逻辑，先把时间步长缩小到原来的三分之一，如果能稳定，说明就是稳定性问题。一维问题 (F_o \le 0.5)，二维问题 (F_o \le 0.25)，这是显式格式的硬性约束，没有商量余地。

另外注意热物性参数的单位。很多人用着用着发现结果怪怪的，翻回去一看，导热系数用的是W/(m·K)，密度用了kg/m³，比热用了J/(kg·K)，算出来的α看着没问题，但单位一混，数量级就偏了。

4.2 边界条件处理：最常见也最容易错

Dirichlet边界（给定温度）实现最简单，直接在每步更新后把边界点的值重新赋值为给定温度就行。Neumann边界（给定热流或绝热）需要从方程中消去边界外的虚拟点，否则会引入额外的误差甚至不稳定。最简单稳妥的办法是用一阶单边差分：绝热边界就是让边界点的温度等于相邻内部点的温度，这个做法虽然精度只有一阶，但胜在稳定可靠。如果对精度有要求，可以用二阶的单边差分，或者在边界外设置虚拟节点并采用中心差分。

对流边界（Robin条件）处理起来最麻烦，它把温度和热流耦合在一起。我的建议是：初学阶段先把Dirichlet和Neumann吃透，对流边界弄清楚之后再碰，不然后期调试时边界和内部方程互相干扰，问题很难定位。

4.3 程序对不对？先找片区域跟解析解对比

代码写完以后，第一件事永远是验证。最可靠的验证方式是找一个有解析解的简单问题来测。一维问题里，四周恒温的无限大平板冷却问题，或者半无限大物体表面温度跳变问题，都有现成的解析解（用误差函数表达）。拿数值解跟解析解在不同时刻对一下，误差在可接受范围内，程序基本就算验证通过了。

还有一个极其实用的技巧：检查能量守恒。对绝热边界问题，整个物体内部的平均温度应该不随时间变化。如果边界是绝热但平均温度在慢慢往下掉，说明边界条件的地方有热量泄漏。这个检查方法不用解析解也能做，而且能快速定位边界实现的问题。

4.4 MATLAB性能优化的几个小技巧

写MATLAB数值仿真跟写通用脚本不太一样。第一，尽量向量化循环。像一维显式格式里 T(2:end-1) = Tn(2:end-1) + Fo * (Tn(3:end) - 2*Tn(2:end-1) + Tn(1:end-2)) 这种写法，比 for i = 2:Nx 逐点更新快得多。第二，在循环外预先构造好系数矩阵，不要在时间步内重复组装。时间循环里只做矩阵求解，性能会提升很大。第三，稀疏矩阵一定要用 spdiags 构造，别用稀疏全矩阵然后 sparse() 转换，内存效率差别巨大。这些习惯养成以后，跑几百乘几百的网格都不在话下。

最后再分享两个小心得

第一，学数值格式最忌讳的就是“只知道代码不知道原理”。每种格式都是前人踩过无数坑之后设计出来的，它解决的问题、付出的代价，都写在数学推导里。建议把显式、隐式、Crank-Nicolson这三种格式的推导过程亲手推一遍，推完你会发现自己能一眼看出别人的代码在用什么格式、大概会有多准。

第二，做热仿真时永远记得问自己一个问题：这个结果符合物理直觉吗？如果一块板在冷却过程中出现了局部温度升高的现象，而你又没有设置内热源，那大概率是数值问题而不是物理问题。带着这个质疑去检查代码，通常能迅速缩小排查范围。

搞懂这些格式之后，你会发现再去看商业软件里那些看似神秘的时间步长设置、松弛因子、收敛阈值，背后其实都是这些最基本的数值方法在起作用。工具在变，但底层的逻辑几十年都没变过。