基于MATLAB的有限体积法对流换热数值求解与实现

简介：这是一份面向热工、能源与航空航天领域学习者的对流换热数值计算MATLAB项目资料，以有限体积法为主线，解决从物理模型建立、偏微分方程离散到计算求解全流程的实际问题，适合本科高年级及工程师快速上手。压缩包内共3个文件：PDF说明对理论框架进行梳理，DOCX计算说明书详述建模与边界条件设定等步骤，MATLAB脚本则直接演示流动与温度场的数值实现，整体仅746KB，轻量易用。资料已有358人学习下载，可配合课程设计、毕业设计或工程自查使用。内容围绕纳维-斯托克斯方程与能量方程展开，重点覆盖Dirichlet、Neumann等边界条件的施加、迭代求解器的选用以及温度场/速度场的可视化验证，能让读者结合代码和文档快速跑通算例，理解对流系数与温度分布的规律，省去从零搭建程序的繁琐过程，同时也为后续开展更复杂换热模拟打下基础。

1. 为什么对流换热要在 MATLAB 里做数值求解

打开压缩包时我在想，一个用 MATLAB 写的“对流换热数值计算”，能比商业软件多讲出什么？拆完heat_convection.m和说明文档后，结论比较明确：这套资料的价值不在于算出了多复杂的几何，而在于把有限体积离散的每一步——界面插值、系数组装、边界条件施加——都压缩到了可以直接追踪的矩阵运算里。对流换热在工程里无处不在，轴承冷却、电子散热、室内自然对流都属于这类问题；能用 MATLAB 把最小可行的求解器写通，理解层次跟只点软件的流形完全不同。对正在做课程设计、毕业设计或准备仿真二次开发的人来说，这份代码是一个很合适的底座。它解决的问题很朴素：给定已知或已解出的流场，求温度场分布和壁面换热系数。

2. 从控制方程到有限体积离散：对流项是误差的主要来源

温度场由对流和扩散共同决定。在不可压缩流动中，能量方程的守恒形式为：

ρ c_p (∂T/∂t + ∇·(uT)) = ∇·(k∇T) + S

u 是速度矢量。如果速度场已经由流场计算给出（比如用 SIMPLE 算法解出的稳态流动），那 T 的方程就是一个线性对流扩散方程，这也是heat_convection.m的主线思路：先把流动当已知，再解温度。理解了这一点再去读代码，就能意识到压力耦合跟温度是分开处理的，丢掉了 N-S 方程里非线性的麻烦，方便先验证热求解部分是否正确。

2.1 有限体积离散：守恒是基本原则

有限体积法不会把偏微分方程直接差分化，而是对每个控制体做积分。对任意控制体 P，时间项和源项乘以体积，界面上的通量写成年对面上的流量 F 和扩散导 D 的组合。离散后的代数方程是：

a_P T_P = a_E T_E + a_W T_W + a_N T_N + a_S T_S + b

其中 a_E、a_W 等系数由扩散导与对流流量的某种组合决定。界面上的物理量无法直接用网格节点值表示，需要做插值，这就是所谓“格式”问题。

2.2 界面插值格式对比：Pé 数决定稳定性

中心差分把界面温度取为两侧节点平均值，精度为二阶，但对流占主导时会导致负系数，迭代求解出现振荡。迎风差分根据流动方向取上游节点值，虽然只有一阶精度，却能保证系数满足对角占优，用迭代法更容易收敛。

这里 Pé = ρ c_p |u| Δx / k，当 Pé 大于 2 时，中心差分会在界面附近产生非物理振荡。工程计算宁可牺牲一阶精度，也要保证对角占优，这也是heat_convection.m采用迎风型系数的基础。

2.3 一维迎风离散的 MATLAB 片段与系数含义

为了说明系数是怎么组装的，下面给出一维迎风组装循环。

% 一维对流扩散稳态：设置边界条件后滚动中间节点 rho = 1.2; cp = 1005; k = 0.026; % 空气物性，SI 单位 L = 1; Nx = 20; dx = L/(Nx-1); u = 0.1 * ones(Nx,1); % 已知速度场 F = rho*cp*u; % 对流流量（1D 简化） D = k/dx; % 扩散导 aP = zeros(Nx,1); aE = zeros(Nx,1); aW = zeros(Nx,1); b = zeros(Nx,1); for i = 2:Nx-1 Fe = F(i); Fw = F(i); % 面值，线性插值后更准确 aE(i) = D + max(-Fe, 0); % 东侧系数 aW(i) = D + max( Fw, 0); % 西侧系数 aP(i) = aE(i) + aW(i); % 对角系数等于邻居之和 b(i) = 0; % 无内热源 end

逻辑说明：当 F > 0 时流动方向是从西到东，上游在西侧，所以西侧系数带上完整对流项，东侧只保留扩散；max函数把方向信息压缩进去，避免写 if-else 分支。按这组物性算，Pé 约为 92.7，中心差分早已无法收敛，迎风此刻仍能给出物理上可接受的单调温度分布。这段逻辑在heat_convection.m里被扩展成二维，系数从数组变成稀疏矩阵 A，边界条件也相应改成绝热或恒温约束。

3. heat_convection.m 的实现：矩阵组装、边界条件与求解器

3.1 程序骨架和网格定义

读heat_convection.m之前，先看说明文档里的流程图。整体流程是：建立矩形网格 → 给定速度场 → 计算每个控制体四边界面上的流量与扩散导 → 组装系数矩阵 A 和右侧向量 b → 施加边界条件 → 求解线性方程组 → 后处理出温度云图和对流换热系数。

网格是结构化矩形网格，节点按列优先编号。因为只做换热部分的计算，速度不参与能量方程内部迭代，传热问题被控制在一个线性方程组里，比完整流固耦合小得多。计算说明书里建议网格从此小到大递进：先用 20×20 跑通，再逐步加密，避免一开始就在大网格上调不出收敛行为。

3.2 二维组装循环：稀疏矩阵是唯一合理的写法

% 二维 FVM 能量方程组装：等距网格，迎风，Dirichlet 边界用大系数法 N = Nx*Ny; A = sparse(N,N); b = zeros(N,1); tol = 1e30; % 大系数，用于固定壁温约束 for j = 2:Ny-1 for i = 2:Nx-1 idx = j + (i-1)*Ny; De = k*dy/dx; Dw = De; Dn = k*dx/dy; Ds = Dn; Fe = rho*cp*u_face_e(i,j)*dy; Fw = rho*cp*u_face_w(i,j)*dy; Fn = rho*cp*v_face_n(i,j)*dx; Fs = rho*cp*v_face_s(i,j)*dx; aE = De + max(-Fe,0); aW = Dw + max( Fw,0); aN = Dn + max(-Fn,0); aS = Ds + max( Fs,0); aP = aE + aW + aN + aS; A(idx,idx) = aP; A(idx,idx+Ny) = -aE; % 东邻居编号差 Ny A(idx,idx-Ny) = -aW; A(idx,idx+1) = -aN; A(idx,idx-1) = -aS; b(idx) = S_rate*dx*dy; % 内热源项 end end

逻辑说明：界面流量 Fe 用速度场在界面上的值乘以界面面积 dy，再乘 ρc_p 变成热容流率。若界面速度为零，Fe=0，方程退化为纯导热，系数就是扩散导，这保证代码能同时覆盖对流和导热两种工况。稀疏矩阵 A 的索引按列优先编号保持物理邻居关系：西邻居编号减 Ny，东邻居编号加 Ny；上下邻居在内存上相邻，增量为 ±1。对流量项用max(-Fe,0)而不是abs(Fe)，是因为迎风逻辑只关心“从哪个方向进入控制体”，符号本身已经包含方向信息。

3.3 边界条件三种施加方式与代码对应

绝热边界是 Neumann 的特例，令 q=0 即可。边界条件施加完毕后，用 MATLAB 内置稀疏直接求解器：

T = A\b; % 直接求解线性方程组 T = reshape(T, Ny, Nx); % 转成物理网格，方便画图

这段求解方式在节点规模小于一万时效率可观。网格加大后需要切到bicgstab(A,b,tol,200)或 GMRES，配合对角占优的迎风矩阵，收敛速度会比默认直接求解更稳定。遇到“矩阵奇异”报错时，先检查是否所有固定壁温边界都加了tol大系数，再检查稀疏矩阵尺寸是否为 Nx*Ny。

4. 验证与参数调试：解析解对标、Pé 数与松弛因子

4.1 用充分发展流场的解析解做对标

验证案例选用二维平行通道：入口给定抛物线速度剖面，壁面恒温 Tw，流体中心温度 Tc。在热充分发展段，无量纲温度剖面接近抛物线分布，可以用解析解做基准。

在heat_convection.m外层包一个测试脚本，计算 L2 相对误差：

% 验收脚本：计算无量纲 L2 误差 T_num = reshape(T, Ny, Nx)'; y = linspace(0, H, Ny); T_exact = (Tw - Tc)*(y/H) .* (1 - y/H) + Tc; L2 = norm(T_num(:)-T_exact(:), 2) / norm(T_exact(:)-Tc, 2); fprintf('L2 relative error = %.3e at %d x %d grid\n', L2, Nx, Ny);

一阶迎风的典型收敛趋势是：网格尺寸减半，误差近似减半。网格加密结果：

这组测试的意义不在误差绝对值，而在于确认程序没有出现中心差分负系数导致的振荡。如果看到 T 在空间上呈锯齿波动，第一检查max(-Fe,0)方向是否写反，第二检查是否漏乘了 ρc_p 流量项。

4.2 松弛因子与残差：怎么判断迭代收敛

直接使用A\b时不需要松弛因子，但扩展到瞬态或流固耦合后，常用 SOR 代替直接求解。SOR 的迭代格式是：

T^(k+1) = T^(k) + ω (T_new - T^(k))

ω 过小收敛慢，过大会发散。工程经验一般先取 0.7 试算，观察残差序列变化：

r = norm(A*T - b, inf); fprintf('residual at step %d: %.2e\n', step, r);

直接求解下残差应当降到 1e-12 以下；若停留在 1e-3 量级，多半是 Nx*Ny 与稀疏矩阵尺寸不匹配，或边界条件没有覆盖全部边界节点。我的常用检查是打印full(A)的若干子块，快速确认组装是否对称、边界行是否只剩对角。

4.3 三类排错：发散、奇异、残差下不去

时间紧的情况下，先用 20×20 网格跑通，再调整到 100×100。真正要看的不是云图漂不漂亮，而是壁面换热系数 h = q_wall / (T_wall - T_bulk) 是否随网格逼近稳定值，这直接引出网格无关性验证。

5. 网格无关性验证与 VTK 导出：把结果落成可复核的文件

5.1 网格无关性验证的计算方法

网格无关性不能只加密一次。分别在 Nx×Ny、2Nx×2Ny、4Nx×4Ny 三套网格上计算目标量，比如出口截面平均温度 T_bulk 或壁面平均 Nu，然后用 Richardson 外推估计收敛比：

R = (f2 - f1) / (f3 - f2)

R 接近 1 表示单调收敛；若 |R| 明显小于 1 且符号振荡，就要回头检查边界条件。网格收敛指数按常用近似公式：

GCI = 3 * |ε| / (r^p - 1)

网格加密比 r=2，迎风离散的理论精度 p=1。三套网格结果合并写成记录表，比单张温度云图更有说服力。

5.2 把温度场导出成 VTK 格式

MATLAB 的surf图适合自己看，但汇报或跟 CFD 结果对比时，导成 VTK 更通用。最小结构化网格导出函数：

function write_vtk(filename, X, Y, T) [Ny,Nx] = size(T); fid = fopen(filename, 'w'); fprintf(fid, '# vtk DataFile Version 3.0\n'); fprintf(fid, 'temperature field\nASCII\nDATASET STRUCTURED_GRID\n'); fprintf(fid, 'DIMENSIONS %d %d 1\n', Nx, Ny); fprintf(fid, 'POINTS %d float\n', Nx*Ny); fprintf(fid, '%g %g 0\n', [X(:)'; Y(:)']); fprintf(fid, 'POINT_DATA %d\nSCALARS T float 1\nLOOKUP_TABLE default\n', Nx*Ny); fprintf(fid, '%g\n', T(:)); fclose(fid); end

把write_vtk挂到主程序尾部，每次跑完自动生成 VTK 文件，配合 ParaView 做切面和剖面叠加，可以将温度场和速度矢量放到同一坐标系下复核，整条调试链路就闭环了。

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