简介：这份资源面向学习数值计算、偏微分方程求解的本科生与工程技术人员，聚焦用Matlab实现有限差分法求解二维热传导问题。内容围绕热传导方程∂u/∂t=α(∂²u/∂x²+∂²u/∂y²)展开，涵盖网格划分、边界条件设定、温度场初始化、时间步长更新，以及用追赶法处理三对角矩阵以提升求解效率，并配有结果可视化，帮助读者理解离散化与数值稳定性等关键环节。资源包共4个文件，以3个m脚本为主，分别承担主程序、追赶法求解与算例演示，另含1张jpg示意图辅助理解，整体约27KB，轻量便于直接运行调试。目前已有4518人学习下载，适合希望掌握有限差分编程实现、对照代码复现二维温度场演化过程的读者参考。
1. 从一张温度云图说起：为什么有限差分法解二维热传导值得你花时间
如果你做过热仿真，大概率遇到过这样的场景：一块金属板左侧被火焰加热，右侧用冷却水冲刷，你想知道板上每一点的温度随时间怎么变。商业软件当然能算，但当你需要反复调参数、嵌入自己的控制逻辑、或者只是想搞明白求解器内部到底在干什么时，自己写一个有限差分法求解器反而是最快、最可控的路径。Matlab 实现有限差分法解二维热传导问题，本质上就是把这个物理过程翻译成矩阵运算：把连续的温度场离散成网格点，用差分公式近似偏导数，再让时间一步步往前推。它适合三类人：正在学数值传热的学生、需要快速验证热设计方案的工程师、以及想把求解器嵌入更大优化循环的开发者。这篇文章不讲空泛的公式推导，而是从网格划分、差分格式、边界条件、稳定性判据一路写到可运行的代码和踩坑记录，让你能直接复现并改造成自己的版本。
2. 二维热传导方程的离散化：从偏微分方程到矩阵方程
2.1 控制方程与定解条件
二维热传导方程的标准形式是：
ρ c ∂T/∂t = ∂/∂x(k ∂T/∂x) + ∂/∂y(k ∂T/∂y) + Q
其中 ρ 是密度，c 是比热容，k 是导热系数，Q 是内热源。对于各向同性、物性为常数的材料，方程简化为：
∂T/∂t = α (∂²T/∂x² + ∂²T/∂y²) + Q/(ρc)
这里 α = k/(ρc) 是热扩散率，单位是 m²/s。要解这个方程，还需要初始条件 T(x,y,0) = T₀(x,y) 和边界条件。常见的边界条件有三类：第一类给定边界温度（Dirichlet），第二类给定热流密度（Neumann），第三类给定对流换热系数和环境温度（Robin）。在 Matlab 里，我一般把温度场存成一个二维矩阵 T，行对应 y 方向，列对应 x 方向，这样后面用矩阵索引操作最顺手。
2.2 显式差分格式的推导与稳定性条件
显式格式是最直观的：时间导数用前向差分，空间二阶导数用中心差分。对内部节点 (i,j)：
T(i,j,n+1) = T(i,j,n) + α Δt [ (T(i+1,j,n) - 2T(i,j,n) + T(i-1,j,n))/Δx² + (T(i,j+1,n) - 2T(i,j,n) + T(i,j-1,n))/Δy² ]
这个式子可以直接写成矩阵运算，但前提是满足稳定性条件。对于二维显式格式，稳定性要求：
α Δt (1/Δx² + 1/Δy²) ≤ 1/2
如果 Δx = Δy = h，则 Δt ≤ h²/(4α)。这是硬约束，不满足就会振荡发散。我见过太多人在这里翻车：网格加密一倍，时间步长忘了缩小四倍，结果算了几步温度就飞到 10 的几十次方。所以实际写代码时，我会先根据网格间距自动计算最大允许时间步长，再乘一个 0.8 到 0.9 的安全系数。
2.3 边界条件的离散处理
边界条件是有限差分法里最容易出错的地方。以左侧边界 x=0 为例，如果是第一类边界条件，直接赋值 T(1,:) = T_left 即可。如果是第二类边界条件，比如给定热流 q，则用一阶向前差分近似：
-k (T(2,j) - T(1,j))/Δx = q
整理得 T(1,j) = T(2,j) + q Δx / k。注意符号：热流流入为正时，边界温度比内部高。如果是第三类边界条件，对流换热：
-k (T(2,j) - T(1,j))/Δx = h (T_env - T(1,j))
整理得 T(1,j) = (k T(2,j) + h Δx T_env) / (k + h Δx)。这些式子看起来简单，但边界上的节点在循环里要单独处理，不能和内部节点混在一起。我习惯把边界更新写成一个独立的函数，输入当前温度场和边界参数，输出更新后的边界值，这样主循环干净，调试也方便。
3. 用 Matlab 把求解器跑起来：网格、时间步与主循环
3.1 网格生成与物性参数设置
先定义几何和物性。假设一块 0.1 m × 0.1 m 的方形板，材料是铝，导热系数 k=237 W/(m·K)，密度 ρ=2700 kg/m³，比热容 c=900 J/(kg·K)，热扩散率 α = k/(ρc) ≈ 9.75e-5 m²/s。网格取 50×50，则 Δx = Δy = 0.1/49 ≈ 2.04e-3 m。根据稳定性条件，Δt ≤ h²/(4α) ≈ (2.04e-3)²/(4×9.75e-5) ≈ 0.0107 s。取 Δt = 0.008 s，安全系数约 0.75。
% 几何与物性参数 Lx = 0.1; Ly = 0.1; % 板尺寸，单位 m k = 237; rho = 2700; c = 900; % 铝的物性 alpha = k / (rho * c); % 热扩散率 % 网格参数 Nx = 50; Ny = 50; % 网格点数 dx = Lx / (Nx - 1); dy = Ly / (Ny - 1); x = linspace(0, Lx, Nx); y = linspace(0, Ly, Ny); [X, Y] = meshgrid(x, y); % 时间步长 dt_max = 1 / (2 * alpha * (1/dx^2 + 1/dy^2)); dt = 0.8 * dt_max; % 安全系数 0.8 t_end = 2.0; % 模拟总时长，单位 s Nt = ceil(t_end / dt);这段代码里，dt_max 是根据稳定性条件算出来的理论上限，实际用的 dt 打了八折。Nt 是总步数，后面主循环用。注意 meshgrid 生成的 X 和 Y 是矩阵，方便后面画图或设置空间变化的初始条件。
3.2 初始条件与边界条件的代码实现
初始温度设成 25°C 均匀，左边界突然升到 200°C，右边界保持 25°C，上下边界绝热。绝热边界用第二类边界条件，热流为零。
% 初始温度场 T = 25 * ones(Ny, Nx); % 单位 °C T_left = 200; % 左边界温度 T_right = 25; % 右边界温度 % 预分配温度场历史（可选，用于动画） T_history = zeros(Ny, Nx, min(Nt, 100)); save_idx = round(linspace(1, Nt, min(Nt, 100))); % 主循环 for n = 1:Nt T_old = T; % 内部节点更新（显式差分） for i = 2:Ny-1 for j = 2:Nx-1 T(i,j) = T_old(i,j) + alpha * dt * ... ((T_old(i+1,j) - 2*T_old(i,j) + T_old(i-1,j))/dx^2 + ... (T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1))/dy^2); end end % 边界条件 T(:,1) = T_left; % 左边界 Dirichlet T(:,end) = T_right; % 右边界 Dirichlet T(1,:) = T(2,:); % 下边界绝热（零热流） T(end,:) = T(end-1,:); % 上边界绝热 % 保存快照 if ismember(n, save_idx) idx = find(save_idx == n, 1); T_history(:,:,idx) = T; end end这段代码里，内部节点用双重循环更新，虽然 Matlab 里可以用向量化加速，但双重循环更直观，也方便你改成非均匀网格或变物性。边界条件在每次迭代后强制赋值，绝热边界用相邻内部节点的温度直接复制，等价于零梯度。注意 T(1,:) 对应 y=0 的下边界，T(end,:) 对应 y=Ly 的上边界，索引别搞反。
3.3 结果可视化与温度云图
算完之后，用 contourf 画温度云图，或者用 imagesc 看温度分布。
figure; contourf(X, Y, T, 20, 'LineColor', 'none'); colorbar; colormap('jet'); xlabel('x (m)'); ylabel('y (m)'); title(['温度分布，t = ', num2str(t_end), ' s']); axis equal tight; % 如果保存了历史，可以画动画 figure; for idx = 1:size(T_history, 3) contourf(X, Y, T_history(:,:,idx), 20, 'LineColor', 'none'); colorbar; caxis([25 200]); title(['t = ', num2str(save_idx(idx)*dt, '%.2f'), ' s']); drawnow; endcontourf 的第四个参数是等值线数量，20 条一般够看。caxis 固定色标范围，动画才不会闪。如果网格更密，比如 200×200，双重循环会明显变慢，这时候可以把内部节点更新改成向量化形式，用 2D 卷积或者稀疏矩阵乘法，速度能快一个数量级。但向量化代码可读性差，我一般先用循环跑通，确认物理结果对了再优化。
4. 避坑指南：显式格式解二维热传导的五个血泪教训
4.1 时间步长超限导致温度振荡发散
现象：算了几十步后，温度场出现棋盘状振荡，数值迅速增大到 1e10 以上。原因：显式格式的稳定性条件被违反，通常是网格加密后忘了缩小时间步长。解决：在代码开头根据 dx、dy 和 alpha 自动计算 dt_max，并打印出来。如果手动设 dt，加一个断言检查 dt <= dt_max。我习惯在循环里每 100 步检查 max(abs(T(:))) 是否超过初始最大值的 10 倍，超过就报错退出，避免浪费时间。
4.2 边界条件索引搞反导致热流方向错误
现象：左边界加热，但温度云图显示热量往左边界外跑，或者边界温度始终不变。原因：Matlab 矩阵索引 (i,j) 中 i 对应行（y 方向），j 对应列（x 方向），但很多人按数学坐标 (x,y) 去写，结果 T(1,:) 被当成左边界，实际是下边界。解决：在代码里明确注释每个边界对应的物理位置，画图时用 meshgrid 生成的 X、Y 验证。我一般会在设置边界后立刻画一张初始温度场，确认边界位置和预期一致再往下算。
4.3 绝热边界用错差分格式引入虚假热流
现象：上下边界绝热，但算出来的温度场在边界附近出现异常梯度，总热量不守恒。原因：绝热边界应该用零热流条件，即 ∂T/∂n = 0，离散后是 T(1,:) = T(2,:)。但有人用一阶向前差分写成 (T(2,:) - T(1,:))/dy = 0，然后解出 T(1,:) = T(2,:)，这没错；但如果用中心差分或者搞错符号，就会引入非零热流。解决：绝热边界直接用相邻节点赋值，不要引入额外的热流项。如果边界是第三类条件，确保 h 和 T_env 的单位一致，h 的单位是 W/(m²·K)。
4.4 物性参数单位不统一导致结果量级错误
现象：温度变化极快或极慢，和理论估计差几个数量级。原因：长度用 mm、导热系数用 W/(m·K)、比热容用 kJ/(kg·K)，单位没统一到 SI。解决：所有输入参数统一换算成 SI 单位：长度 m，温度 K 或 °C（差分方程里用 °C 没问题，因为只涉及温差），时间 s，功率 W。热扩散率 α 的单位是 m²/s，算出来之后检查一下数量级，金属一般在 1e-5 到 1e-4 之间，塑料在 1e-7 左右。
4.5 双重循环效率过低导致大规模网格算不动
现象：网格 200×200 时，算 1000 步要等好几分钟。原因：Matlab 的 for 循环在旧版本里很慢，即使新版本有 JIT 加速，双重循环仍然不如向量化。解决：把内部节点更新写成向量化形式，用切片操作：
T(2:end-1, 2:end-1) = T_old(2:end-1, 2:end-1) + alpha * dt * ... ((T_old(3:end, 2:end-1) - 2*T_old(2:end-1, 2:end-1) + T_old(1:end-2, 2:end-1))/dx^2 + ... (T_old(2:end-1, 3:end) - 2*T_old(2:end-1, 2:end-1) + T_old(2:end-1, 1:end-2))/dy^2);这样一行顶双重循环，速度快很多。但要注意切片索引的边界，别越界。向量化之后，边界条件仍然单独赋值，顺序不能乱。
5. 从能算到算得准：验证、加速与进阶技巧
5.1 用一维解析解验证二维求解器
验证求解器最可靠的方法是和解析解对比。考虑一个一维瞬态热传导问题：初始温度均匀 T₀，左边界突然升到 T_s，右边界保持 T₀，无内热源。解析解是无穷级数，但取前几项就够。在二维求解器里，把上下边界设成绝热，左右边界设成 Dirichlet，且初始温度只沿 x 方向变化，那么温度场应该只沿 x 方向变化，退化成 一维问题。把数值解和解析解画在同一张图上，如果最大误差在 1% 以内，说明差分格式和边界条件都对了。我一般会取 x 方向中间一条线，对比不同时刻的温度分布，这样最直观。
5.2 隐式格式与 Crank-Nicolson 的取舍
显式格式受稳定性限制，时间步长小，算长时间瞬态问题时步数太多。隐式格式无条件稳定，可以用大时间步，但每步要解线性方程组。在 Matlab 里，二维隐式格式的系数矩阵是块三对角矩阵，可以用稀疏矩阵和反斜杠直接求解。对于 50×50 网格，隐式格式每步求解大约 0.01 秒，而显式格式每步 0.001 秒但需要 10 倍步数，总时间差不多。但如果网格到 200×200，隐式格式的优势就出来了。Crank-Nicolson 格式精度更高，但可能引入非物理振荡，尤其是边界条件突变时。我一般先用显式格式快速验证物理模型，确认无误后再换隐式格式做长时间模拟。
5.3 变物性与内热源的处理
实际工程中，导热系数可能随温度变化，或者有内热源。变物性处理起来很简单：在每个时间步更新 k 和 α，然后重新计算稳定性条件。如果 α 变大，dt 要相应减小。内热源 Q 直接加到更新公式里：T(i,j) = T_old(i,j) + dt * (α * 拉普拉斯项 + Q/(ρc))。注意 Q 的单位是 W/m³，如果给的是热流密度 W/m²，要除以体积或者换算成源项。我见过有人把 Q 直接当成温度增量加上去，结果量级完全不对。
5.4 用稀疏矩阵加速隐式求解
隐式格式的核心是求解 A * T_new = b，其中 A 是稀疏矩阵。对于二维问题，可以用 kron 函数构造：
% 构造二维拉普拉斯矩阵（隐式格式用） e = ones(Nx, 1); A1 = spdiags([e -2*e e], -1:1, Nx, Nx); I = speye(Nx); A = kron(I, A1) + kron(A1, I); % 注意 kron 的顺序然后每个时间步更新右端项 b，用 T_new = A \ b 求解。稀疏矩阵求解在 Matlab 里非常快，200×200 网格也能秒级完成。但要注意 kron 的顺序：如果 T 按行优先存储，kron(I, A1) 和 kron(A1, I) 的位置要对调。我一般会先用小网格测试，确认矩阵和手动循环的结果一致，再放大网格。
5.5 一个容易忽略的细节：温度单位与差分符号
最后说一个玄学问题：温度用 °C 还是 K。差分方程里只有温差，所以 °C 和 K 等价。但如果你用 °C 算热辐射或者温度相关物性，记得换算成 K。另外，热流方向：傅里叶定律是 q = -k ∇T，负号不能丢。在边界条件里，如果热流流入计算域，q 为正，边界温度比内部高；如果热流流出，q 为负。我习惯在代码里统一用“流入为正”的约定，然后在注释里写清楚，避免符号混乱。
写到这里，这个求解器已经能跑、能验证、能加速了。我自己的习惯是：每改一个参数，先跑 10 步看温度场是否合理，再跑完整模拟。如果结果不对，先检查边界条件索引和单位，再检查稳定性条件。希望帮到你。
本文还有配套的精品资源，点击获取
