4. 二维稳态导热有限差分法

二维稳态导热问题的控制方程是拉普拉斯方程（无内热源）或泊松方程（有内热源）。有限差分法的核心思路是：把连续的温度场 T(x,y) 离散成有限个节点上的温度值，用差商代替导数，把偏微分方程转化为一组代数方程，再求解这组方程得到各节点温度。本节围绕“网格怎么划、节点方程怎么写、方程组怎么组装和求解”这条主线展开，并以一个矩形平板导热的算例贯穿说明。

一、网格划分与节点编号

设矩形区域宽 W、高 H，沿 x 方向划分 N_x 个区间、沿 y 方向划分 N_y 个区间，则节点数为 (N_x+1)×(N_y+1)，步长分别为 Δx = W/N_x、Δy = H/N_y。为简化书写，通常取均匀网格 Δx = Δy = h。节点用 (i, j) 表示，i 为 x 方向列号（0～N_x），j 为 y 方向行号（0～N_y）。

编号方式直接影响方程组系数矩阵的结构。常用做法是“按行扫描”：先固定 j，令 i 从 0 到 N_x 依次编号，再换下一行。这样同一行的相邻节点编号相差 1，上下相邻行对应节点编号相差 (N_x+1)。对于二维问题，每个内部节点的方程只涉及自身和上下左右四个邻居，因此系数矩阵是稀疏的“五对角”带状矩阵，这一结构对后续选择求解器很关键。

二、内部节点与边界节点的差分格式

对内部节点 (i, j)，用中心差分近似二阶导数：

∂²T/∂x² ≈ (T[i+1,j] − 2T[i,j] + T[i−1,j]) / h²
∂²T/∂y² ≈ (T[i,j+1] − 2T[i,j] + T[i,j−1]) / h²

代入 Δx = Δy = h 的拉普拉斯方程，整理得内部节点的五点格式：

T[i,j] = (T[i+1,j] + T[i−1,j] + T[i,j+1] + T[i,j−1]) / 4

即内部节点温度等于四邻节点温度的算术平均。若存在内热源 q_v（泊松方程），则多出一项 h²q_v/k，形式为 T[i,j] = (四邻之和 + h²q_v/k)/4。

边界节点的处理取决于边界条件类型。第一类边界（给定温度）最简单：边界节点温度已知，不参与未知量求解，只作为相邻内部节点方程中的已知项。第二类边界（给定热流密度 q″）需要引入虚拟节点或采用单侧差分。以右边界 x = W 处给定热流 q″ 为例，用虚拟节点 (N_x+1, j) 表示边界外侧，由能量守恒可得：

T[Nx,j] = (2T[Nx−1,j] + T[Nx,j+1] + T[Nx,j−1] + 2h·q″/k) / 4

第三类边界（对流换热，换热系数 h_c、环境温度 T_∞）的边界节点方程需同时包含导热和对流项，以右边界为例：

T[Nx,j] = (2T[Nx−1,j] + T[Nx,j+1] + T[Nx,j−1] + 2h·h_c·T∞/k) / (4 + 2h·h_c/k)

注意这里的 h 是网格步长，h_c 是对流换热系数，两者符号相近但含义完全不同，编程时务必用不同变量名区分，这是初学者最容易出错的地方之一。

三、线性方程组组装

把所有未知节点（内部节点 + 第二、三类边界节点）按编号顺序排列成向量 T，每个节点写一个方程，最终得到 A·T = b。组装时逐节点处理：若某邻居是已知温度的第一类边界节点，把该项移到右端项 b 中；若邻居是未知节点，则保留在系数矩阵 A 中。以内部节点为例，方程写成标准形式：

4·T[i,j] − T[i−1,j] − T[i+1,j] − T[i,j−1] − T[i,j+1] = 0

若左邻节点是已知温度 T_L，则方程变为 4·T[i,j] − T[i+1,j] − T[i,j−1] − T[i,j+1] = T_L。这样 A 的对角元为 4（边界节点可能不同），非对角元为 −1，矩阵对称正定，这是五点格式的重要性质。

以一个具体算例说明：取 1 m × 1 m 方形平板，k = 1 W/(m·K)，划分 4×4 网格（h = 0.25 m），共 25 个节点。设左边界 T = 100 ℃，右边界 T = 0 ℃，上下边界绝热（q″ = 0）。未知量为 3 列 × 5 行 = 15 个内部及绝热边界节点。绝热边界节点方程中缺少外侧邻居项，系数相应减少，例如上边界节点：

3·T[i,0] − T[i−1,0] − T[i+1,0] − T[i,1] = 0   （i 为内部列）

角点 (0,0) 同时受两个绝热边界影响，方程变为 2·T[0,0] − T[1,0] − T[0,1] = 0。组装完成后，A 为 15×15 的稀疏矩阵，b 中只有靠近左边界 T = 100 ℃ 的节点对应非零项。

四、直接求解与结果校核

对于中小规模网格（节点数几百到几千），直接法（高斯消元、LU 分解）稳定可靠。Python 中可用 numpy.linalg.solve 直接求解，或用 scipy.sparse.linalg.spsolve 利用稀疏性加速。求解后应做两项校核：一是把解代回原方程检查残差 ‖A·T − b‖ 是否接近机器精度；二是检查物理合理性，例如温度应介于边界温度之间、等温线不应出现异常振荡。

上述算例若用 4×4 网格求解，中心区域温度大致在 50 ℃ 附近，靠近左边界处温度梯度较大，靠近右边界处较平缓，符合一维导热的定性规律。若把网格加密到 20×20，中心温度会略有变化并趋于收敛值，这正是网格无关性验证的意义。实际工程中建议至少做一次网格加密对比，确认关键节点温度变化小于允许误差后再采用该网格结果。