MATLAB流体与传热仿真程序解析：从离散格式到收敛验证

简介：MATLAB流体计算与传热程序包，面向需要开展CFD数值模拟或热分析的研究者、工程师及高年级本科生，可解决流动场求解、温度分布预测等多类典型问题。压缩包内含11个.m文件，整体约10KB，全部为可直接运行的MATLAB脚本，兼顾流体流动模拟与热传导、对流等传热问题；文件虽少，但覆盖了有限差分、有限体积法以及TDMA、SOR迭代等核心算法，既有基础算例，也有二维扩展，适合配合《数值传热学》或《计算流体力学导论》教材同步学习。目前已有275人学习下载，常用于课程设计、课题验证或作为自定义求解器的起点。通过阅读示例脚本，可快速理解离散化思路与边界条件设置；修改物性参数即可复现不同工况下流场与温度场的变化，并借助MATLAB自带绘图工具可视化计算结果。对希望掌握流体计算与传热具体编程实现的读者而言，这套代码提供了清晰的对照模板，是入门数值模拟的实用起点。

1. 这份 MATLAB 流体与传热程序，到底能拿来干什么

接手一个名为“fluid heat_matlab”的压缩包，里面装着几十个 .m 文件，第一反应通常是兴奋，第二反应往往是翻车：变量名是拼音缩写、没有子函数调用说明、跑出来的云图怎么看都像一张彩虹噪点图。这类 MATLAB 编写的流体计算和传热程序，绝大多数是教学型或预研型的二维简化模型，用有限差分或有限体积法把连续性方程、动量方程和能量方程离散成代数方程组，再用矩阵迭代求解。它解决的核心问题是：在不需要商业 CFD 软件授权、不想被 Fluent 的网格划分折腾半天时，快速验证一个流动与传热耦合的物理过程到底能不能收敛、趋势对不对。

适合读这篇文章的人，是做课程设计、毕业设计预研、实验台参数估算的本科生和研究生，以及刚转行到仿真方向的工程师。你能从一份别人写好的 MATLAB 程序里拿到的不只是结果云图，更是一套“如何把物理方程变成可迭代的代数系统”的完整思路。我接下来会按从方程到代码、从流动到传热、从踩坑到验证的顺序，把这类程序的阅读、修改和复现方法拆开讲清楚。

2. 从控制方程到差分方程：读懂程序前先看清离散框架

2.1 控制方程与差分格式的选用逻辑

流体计算和传热程序的数学起点是三个守恒方程：质量守恒、动量守恒和能量守恒。在二维不可压缩、常物性的简化假设下，方程组可以写成 N-S 方程与对流-扩散型能量方程。绝大多数教学程序不会直接求解原始变量形式的 N-S 方程，而是做两步简化：要么引入流函数-涡量法消去压力梯度项，要么采用 SIMPLE 算法在交错网格上做压力修正。选择哪种方案，直接决定你后续读代码时看哪些数组、调哪些参数。

差分格式的选择看 Peclet 数和雷诺数。中心差分在网格 Peclet 数小于 2 时精度高，超过后会产生非物理振荡；迎风格式无条件稳定但数值耗散大。程序里如果出现类似D = Dc*(T(i+1,j) - 2*T(i,j) + T(i-1,j))/(dx^2)这样的行，就是典型中心差分，而if(u(i,j)>0) T_left = T(i-1,j); else T_left = T(i+1,j); end则是迎风逻辑。读这份程序时，第一件事不是看每一行，而是用grep或 MATLAB 的查找功能搜if分支和差分系数，确认作者混用了哪些格式，这决定了你在后续修改时改哪里。

2.2 一维非稳态传热的最小可复现算例：隐式格式怎么写

在直接啃二维流动耦合代码之前，我强烈建议先在一维问题上跑通隐式离散。大多数 fluid heat 程序包里都会有一段一维热传导的启动代码，它虽然简单，却是理解整个程序矩阵组装方式的钥匙。隐式格式的离散方程长这样：

% 一维非稳态导热，隐式格式 % 参数定义：L=1m 杆长，N=50 网格，alpha=1e-4 热扩散系数 % dt=0.1s，t_end=10s，两端边界分别为 0 和 100 摄氏度 L = 1.0; N = 50; dx = L/(N-1); alpha = 1e-4; dt = 0.1; t_end = 10.0; x = linspace(0, L, N)'; T = zeros(N, 1); T(end) = 100; % 右端固定为100摄氏度 % 组装系数矩阵 A 和右端项 b，采用隐式（后向欧拉） % 每步求解 A * T_new = T_old + 边界贡献 A = sparse(N, N); b = zeros(N, 1); r = alpha * dt / dx^2; % 网格傅里叶数（扩散稳定性参数） % 内部节点方程: -r*T(i-1) + (1+2r)*T(i) - r*T(i+1) = T_old(i) for i = 2:N-1 A(i, i-1) = -r; A(i, i) = 1 + 2*r; A(i, i+1) = -r; end % 边界节点：左端0摄氏度（第一类边界），右端100摄氏度 A(1,1) = 1; b(1) = 0; A(N,N) = 1; b(N) = 100; for t = 0:dt:t_end T = A \ b; % 直接用左除求解稀疏三对角系统 % 每步更新右端项：内部节点取上一时刻的 T b(2:N-1) = T(2:N-1); end plot(x, T); title('隐式格式温度分布');

这段代码逻辑上分两步：先组装一次系数矩阵，因为网格和边界条件不变时不需每步重建；后进入时间循环，每步左除求解。关键参数是r = alpha*dt/dx^2，它叫网格傅里叶数。显式格式要求r <= 0.5才能稳定，隐式格式原则上对r没有上限，但太大的dt会让时间方向截断误差变大，温度曲线会偏离真实解。我一般会先用dt = dx^2/(2*alpha)算一个参考值，再逐步放大，观察曲线是否产生阶梯状畸变。

2.3 三对角追赶法：为什么老程序里总有个 tdma 函数

几乎所有老牌 MATLAB 流体传热程序里都会有一个tdma.m或tridiag.m的独立文件，上面代码里我用了A \ b直接求解，但在几十年前的 MATLAB 版本里，稀疏矩阵解法的效率远不如今天，作者们习惯手写 Thomas 算法。这个算法只有三步：消去下对角线、更新右端项、回代求未知量。对 50 个网格而言两者都能瞬解，但如果你的研究把网格加到了 5000 个，A \ b仍然很快，因为A声明为sparse后 MATLAB 会自动走稀疏 LU 分解；反倒是手写 tdma 在 FORTRAN 风格迁移的代码里经常出现索引 bug。

当你在压缩包里看到 tdma 函数时，先别急着删掉换成左除。它往往还承担着另一个任务：把三对角求解嵌入迭代循环里，每次只更新一条网格线，这叫线松弛。在二维问题里，交替方向隐式（ADI）格式会把 x 方向和 y 方向各扫一遍，每步都解一个三对角系统。这种写法在今天依然有价值，因为它的内存占用恒定，不需要存储整个二维稀疏矩阵。读这种代码时，重点看它a/b/c/d四个输入向量是不是按网格线顺序填充的，边界点是否被重复计算成了未知量。

3. 让流体动起来：从流函数-涡量法看懂求解器骨架

3.1 为什么教学程序偏爱流函数-涡量法

压缩包里的流动求解部分，十有八九是流函数-涡量法，而不是 SIMPLE 算法。原因很实际：二维不可压缩流场用流函数-涡量表述后，连续性方程自动满足，压力项通过交叉求导消掉，最后只剩两个输运方程——涡量输运方程和流函数泊松方程。MATLAB 写起来非常紧凑：涡量场omega和流函数场psi都是按网格索引访问的矩阵，迭代逻辑直观，画等值线也方便。

代价是边界条件处理比原始变量法更绕：固壁上的涡量不能直接给定，要由流函数和法向速度梯度推导，常见公式是omega_wall = (psi_wall - psi_inner) * 2 / dn^2。程序里如果有一段omega(1,:) = 2*(psi(1,:) - psi(2,:))/dy^2;之类的赋值，那就是在更新边界涡量。这个公式的系数 2 和分母的dn^2是根据一阶单边差商加泰勒展开推导出来的，改网格间距时必须同步修改，否则边界层的涡量会产生几个量级的误差，流场速度分布整个错掉。

3.2 方腔流求解核心代码：迭代骨架与松弛参数

顶盖驱动方腔流是这类程序的标配测试算例。下面这段代码展示了流函数-涡量法的主迭代结构，包含内迭代和外迭代两层，适合接在传热方程之前做速度场基础。

% 顶盖驱动方腔流，流函数-涡量法核心迭代骨架 % 参数设置 Nx = 51; Ny = 51; % 计算网格（含边界） Lx = 1.0; Ly = 1.0; % 方腔尺寸 dx = Lx/(Nx-1); dy = Ly/(Ny-1); Re = 100; % 雷诺数，越小越容易收敛 U_top = 1.0; % 顶盖移动速度 omega = zeros(Nx, Ny); % 涡量场（内节点有效） psi = zeros(Nx, Ny); % 流函数场 dt = 0.005; % 虚拟时间步，受CFL条件约束 omega_tol = 1e-6; % 涡量内迭代收敛判据 max_iter = 8000; % 流函数泊松方程的松弛因子：0~2，超松弛可取1.7附近 % 涡量输运推进采用显式格式 for iter = 1:max_iter % 1) 由涡量求流函数：SOR迭代解 laplacian(psi) = -omega for isor = 1:20 psi_old = psi; for i = 2:Nx-1 for j = 2:Ny-1 psi(i,j) = 0.25 * (psi(i+1,j) + psi(i-1,j) + ... psi(i,j+1) + psi(i,j-1) + ... dx*dx * omega(i,j)); end end % 对内部点再施加几何平均超松弛 psi = psi_old * (1 - 1.7) + psi * 1.7; end % 2) 由流函数更新各点速度（二阶中心差分） u = zeros(Nx, Ny); v = zeros(Nx, Ny); u(2:Nx-1, 2:Ny-1) = (psi(2:Nx-1, 3:Ny) - psi(2:Nx-1, 1:Ny-2)) / (2*dy); v(2:Nx-1, 2:Ny-1) = -(psi(3:Nx, 2:Ny-1) - psi(1:Nx-2, 2:Ny-1)) / (2*dx); % 3) 更新边界涡量 omega(1,:) = 2*(psi(2,:) - psi(1,:)) / dx^2; omega(Nx,:) = 2*(psi(Nx-1,:) - psi(Nx,:)) / dx^2; omega(:,1) = 2*(psi(:,2) - psi(:,1)) / dy^2; omega(:,Ny) = -2*(psi(:,Ny-1) - psi(:,Ny)) / dy^2; % 注意顶盖是移动壁面，涡量有额外项: -2*U_top/dx omega(Nx,:) = omega(Nx,:) - 2*U_top/dx; % 4) 涡量输运方程显式推进（对流+扩散） omega_new = omega; for i = 2:Nx-1 for j = 2:Ny-1 conv = -u(i,j)*(omega(i+1,j)-omega(i-1,j))/(2*dx) ... -v(i,j)*(omega(i,j+1)-omega(i,j-1))/(2*dy); diff = (omega(i+1,j)+omega(i-1,j)-2*omega(i,j))/dx^2 ... + (omega(i,j+1)+omega(i,j-1)-2*omega(i,j))/dy^2; omega_new(i,j) = omega(i,j) + dt*(conv + diff/Re); end end omega = omega_new; % 5) 判断涡量残差是否达到容差 if max(max(abs(omega-omega_new))) < omega_tol break; end end

代码里最关键的是两处数值处理。第一处是边界涡量的方向符号：顶盖移动壁面那行加了一个- 2*U_top/dx，符号反了会导致顶盖下方形成错误方向的回流涡。第二处是 SOR 超松弛因子的位置：我把它放在内层 20 次迭代之后，而不是每次更新后立刻松弛，这是为了避免松弛过程与泊松方程收敛互相干扰。流函数-涡量法的迭代速度完全取决于这两个参数配合，Re超过 500 时不缩小dt，程序大概率当场发散，现象是流函数等值线出现锯齿状断裂。

3.3 判断收敛的标准：不要只看残差曲线

程序跑完输出一张图，残差曲线一路下滑，看上去很完美，但很多时候这是假收敛。我的经验是至少看三个指标：涡量残差、流函数残差、以及计算的 Nusselt 数是否随迭代稳定在某个值到三位小数。很多老程序只输出涡量残差，因为作者当年只关心流场；你现在要耦合传热的话，必须额外监测热通量或壁面平均 Nusselt 数。这个值如果每次迭代都在第二或第三位小数震荡，说明流场还存在低频振荡，直接去算温度场会把温度云图弄出带状纹理。

4. 传热方程的耦合：让流体把热量捎带走

4.1 能量方程离散：中心差分与迎风的临界判断

流动解收敛后，温度场在同一个流场上求解。能量方程形式上是一个纯对流-扩散方程，它的关键无量纲数是 Peclet 数Pe = u*dx/alpha。当Pe > 2时，中心差分会产生空间振荡，解出来的温度会在高温热源旁边出现负值。这个负温度在物理上不成立，但程序不会报错，它只是安静地画出蓝绿色斑块。

处理办法是把它改成迎风差分。迎风逻辑在 MATLAB 里最简洁的写法不是 if 判断，而是用max(u,0)和min(u,0)把对流项拆成正负两部分。我在下面代码里给出这种写法，它能顺手完成向量化，比循环快一个数量级以上。

4.2 在已收敛流场上求解温度场的实用代码

% 在同一流场上求解稳态温度场（迎风差分，标量温度 T） % 前置：需要 u, v 速度场已由流函数-涡量法收敛 Nx = size(u,1); Ny = size(u,2); alpha = 0.01; % 热扩散系数 dx = Lx/(Nx-1); dy = Ly/(Ny-1); T = zeros(Nx, Ny); % 初始温度 T(:,1) = 1.0; % 左边界恒温1，其余边界绝热或等温视问题而定 relax_T = 0.6; % 温度场欠松弛因子，防止发散 T_tol = 1e-8; max_T_iter = 20000; % 边界已知时只需迭代内部节点 for iter = 1:max_T_iter T_old = T; % 中心扩散项（二阶中心差分） d2Tdx2 = (T(3:Nx,2:Ny-1) + T(1:Nx-2,2:Ny-1) - 2*T(2:Nx-1,2:Ny-1))/dx^2; d2Tdy2 = (T(2:Nx-1,3:Ny) + T(2:Nx-1,1:Ny-2) - 2*T(2:Nx-1,2:Ny-1))/dy^2; % 对流项：一阶迎风 conv_x = max(u(2:Nx-1,2:Ny-1), 0) .* (T(3:Nx,2:Ny-1) - T(2:Nx-1,2:Ny-1))/dx ... + min(u(2:Nx-1,2:Ny-1), 0) .* (T(2:Nx-1,2:Ny-1) - T(1:Nx-2,2:Ny-1))/dx; conv_y = max(v(2:Nx-1,2:Ny-1), 0) .* (T(2:Nx-1,3:Ny) - T(2:Nx-1,2:Ny-1))/dy ... + min(v(2:Nx-1,2:Ny-1), 0) .* (T(2:Nx-1,2:Ny-1) - T(2:Nx-1,1:Ny-2))/dy; L_T = alpha*(d2Tdx2 + d2Tdy2) - (conv_x + conv_y); % 点隐式更新（Jacobi 欠松弛） T(2:Nx-1, 2:Ny-1) = T(2:Nx-1, 2:Ny-1) + relax_T * (L_T - T(2:Nx-1, 2:Ny-1)); % 边界条件：Dirichlet 左边界固定，其余设为第二类绝热（法向梯度为0） T(1,:) = 1.0; T(Nx,:) = T(Nx-1,:); T(:,Ny) = T(:,Ny-1); if max(max(abs(T - T_old))) < T_tol break; end end

这段代码的思想叫分离求解：在同一个速度场上迭代温度，每步更新后重新施加强边界条件。迎风项用min/max函数做方向分裂，同样是Pe > 2，中心差分会震荡，迎风则只会增加一点数值耗散，不会出现负温度。relax_T是欠松弛因子，取值在 0.5~0.8 之间比较稳妥，取 1.0 的话在局部速度大的区域会出现温度场通道间振荡。最终的温度分布你需要根据物理直觉做检查——在方腔流中，热量应该被顶盖的剪切运动带向下游，形成明显的倾斜等温线，而不是一圈一圈的同心圆。

4.3 耦合迭代的必要性：什么时候必须回流场

上面代码假设流场已收敛且不随温度改变，这叫单向耦合，适用于常物性流体。如果程序里出现了mu = mu0*(1 - beta*(T-T0))这一类行，说明作者考虑了浮力或变粘度，流场和温度场必须双向耦合。强耦合的迭代策略是：每做 10~20 步温度迭代，就用新温度修正物性，再重新跑一轮流函数-涡量迭代。这种嵌套循环在 MATLAB 里很耗时间，我通常会在温度场迭代前先判断整个温差的相对变化，如果相对变化小于千分之一，就断开物性更新，让程序自动退化为单向耦合，以节省后半段的计算时间。

5. 避坑指南：MATLAB 流体传热程序最常见的五个翻车点

5.1 现象：全程矩阵运算仍慢得离谱，网格加到 200 以上每步耗时好几秒

原因：老程序大量使用zeros(Nx,Ny)预分配，但在迭代循环体内部反复创建临时数组，或对整矩阵做feval和匿名函数递归调用，触发 MATLAB 的 JIT 编译降级。解决：第一件事打开 MATLAB 的代码分析器（profile on然后跑一小段），看哪一行消耗占比最高。常见做法是不改算法，只把循环内的omega = function_handle(omega, params)展开成显式数组运算法，速度能提升 5~10 倍。

5.2 现象：流函数等值线图左右不对称，明显该是方腔涡却画成了螺旋

原因：边界涡量公式的符号或者分母长度用混了。尤其是顶盖移动壁面，dx和dy在网格等距时差异不大，但一旦程序支持非均匀网格，dx在近壁处和中心处完全不同，而作者只写了一个统一常数。解决：把边界赋值那一行改用2*(psi临近点 - psi壁面) / (壁面到临近点实际间距)^2，并且打印一遍边界涡量的最大绝对值，正常方腔流应在 10 到 100 量级，如果出现 1e4 以上的尖峰，先去查符号。

5.3 现象：温度场在网格加密后反而出现波浪形振荡

原因：中心差分在对流占优区域失效，Pe 数大于 2。解决：检查加密后的网格间距dx，计算局部最大速度乘dx/alpha。如果超过 2，把那一行的差分改成迎风。不要全局迎风，全局迎风会把整个温度场耗散得过于平滑，等温线细节全部丢失。

5.4 现象：迭代残差曲线一开始下降，到某个值后横盘不动，甚至缓慢上升

原因：松弛因子选取不当。SOR 因子的最优值接近 2/(1+pi/N)，网格越密越接近 2，但超过临界值就震荡；温度场的欠松弛因子取 1.0 也会导致同样的横盘。解决：每次都做一次松弛因子扫描实验，在同一网格下分别取omega_sor = 1.5, 1.7, 1.85，跑 200 步看残差，选收敛最快且不震荡的值。这是这类程序调试中最常见的“玄学”环节，本质上不是玄学，是特征值估计没做。

5.5 现象：加入传热模块后原流动程序的计算结果变了，但温度没变

原因：传热模块引用了全局变量u、v，而主程序在每次迭代后修改了速度场的数组名或索引偏移，导致传热模块读取的是一份过期数据。解决：用 MATLAB 的结构体把速度场、网格参数、物性参数打包成一个params结构体传入子函数，避免全局变量名冲突。我见过最离谱的一次，是传热子函数里的v和主程序控制流程里的逻辑标志v（valid 的意思）重名，程序不报错，但温度场永远只有边界那一条线是热的。这种血泪教训，值得你在压缩包代码里提前用搜索功能排查一遍所有单字母变量。

6. 验证与改造技巧：怎么让这份程序从“能跑”变成“可信”

6.1 网格无关性验证的实操步骤

拿到任何程序，第一件事不是改物理参数，而是做网格无关性验证。选三个网格系列，比如 41×41、81×81、161×161，固定雷诺数和边界温度，分别计算壁面平均 Nusselt 数Nu = q_w * L / (k * (T_hot - T_cold))，把三个值列成表格对比。如果相邻两档的差异小于 1%，就用中间档做正式计算；如果差异大于 5%，说明当前网格完全不够，程序给出的流场细节全是伪结构。

6.2 验证解的小技巧：拿二维方腔流文献值做基准

顶盖驱动方腔流在 Re=100、Re=400、Re=1000 时的方腔中心涡量极值、涡心位置都有公开发表的基准数据。我一般在程序里加一行断言：如果中心涡量极值与 Ghia 等基准值差超过 2%，就打印警告。这比肉眼看云图可靠得多。另一个验证思路是把它和 Fluent 在同一网格下的结果做横向对比，但注意 Fluent 默认用有限体积法和交错网格，和流函数-涡量法的对比误差在 1%~3% 以内就算正常，超过 5% 优先怀疑自己的边界条件写错。

6.3 从脚本进化成工具箱的三个改造方向

第一，把所有物理参数集中到文件头部一个struct里，删除散落在各子函数里的magic number。第二，把迭代解函数从脚本改成带输入输出参数的函数式，方便你批量扫描不同 Re、Pr 数。第三，把每步迭代的残差、最大速度点、平均 Nu 存入一个 history 矩阵，跑完直接plot出收敛曲线。我个人的血泪经验是，凡是超过 100 行的 MATLAB 流体程序，都不要用脚本硬写；哪怕只是把变量名从a1改成u_eddy，后续调试时间能省一半。希望这篇文章能帮你把压缩包里那份程序真正读透，也希望你在改完代码后，能感觉到数值仿真从“黑匣子”变成“看得见每一步的透明过程”。

本文还有配套的精品资源，点击获取