简介这是一份基于MATLAB的对流换热数值计算项目面向学习计算传热学、需要完成课程设计或入门CFD的工科学生与工程师。资源共3个文件包含可运行的.m脚本、PDF说明文档和DOCX计算说明书压缩包仅746KB下载与解压都很轻量。内容从对流换热基本理论出发系统覆盖有限体积法离散、边界条件如定温、定热流等施加、离散方程组的迭代求解以及温度场后处理等关键环节可帮助读者理解如何将纳维-斯托克斯方程和能量方程转化为可执行的MATLAB代码并掌握GMRES等求解器的实际用法。附带的说明文档还梳理了数学模型构建、离散过程与验证分析思路可作为独立完成同类数值实验的参考模板。目前已有358人学习下载适合配合教材或课程进行实践训练。1. 为什么“对流换热数值计算”值得用 MATLAB 从头做一遍对流换热数值计算听起来像是 Fluent 的活但真正会搜到这个标题的人手里的需求往往很小要么是课程设计要一个二维方腔自然对流的可运行程序要么是做热分析却拿不到商业软件授权要么已经用热分析软件算完想用 MATLAB 独立验证壁面换热系数。反直觉的一点是网格只有 41×41 就能把经典方腔算例跑出有意义的 Nusselt 数而整套流程在 MATLAB 里只需要流函数、涡量和温度三个二维矩阵。相比 Ansys 瞬态热分析它没有网格划分和求解器配置的门槛相比手推理论解它又能给出全场分布。这套方法适合能写循环、但没接触过 CFD 工具链的工程师和研究生。下面从控制方程讲到数据导出把这条最简路径完整走通。2. 对流换热控制方程先化成 MATLAB 能算的流函数-涡量格式2.1 能量方程与动量方程怎么放进同一套网格先写二维不可压缩流体的控制方程。常压空气自然对流可以当作常物性处理动量方程里的浮力用 Boussinesq 近似即除浮力项之外密度为常数。能量方程写成∂T/∂t u∂T/∂x v∂T/∂y α(∂²T/∂x² ∂²T/∂y²)涡量方程写成∂ω/∂t u∂ω/∂x v∂ω/∂y ν∇²ω gβ ∂T/∂x强制对流没有最后这个浮力源项自然对流的全部驱动力都在这里。直接求解原始变量 u、v、p 会遇到压力修正问题MATLAB 里最常见的做法是换成流函数-涡量法u ∂ψ/∂yv -∂ψ/∂x涡量 ω -∇²ψ。这样连续性方程自动满足压力项也在涡量方程里消掉剩下两个输运方程加一个泊松方程。网格上每个节点只需要存 ψ、ω、T 三个数省掉了 SIMPLE 那套压力迭代逻辑。无量纲化时自然对流的控制参数是瑞利数 Ra gβΔTL³/(να) 和普朗特数 Pr ν/α。如果用 α/L 作为特征速度控制方程整理成下面这组标准化形式∂T/∂t u∂T/∂x v∂T/∂y ∇²T∂ω/∂t u∂ω/∂x v∂ω/∂y Pr∇²ω Ra·Pr·∂T/∂x∇²ψ -ω注意看能量方程里已经没有 α涡量方程里也没有 g 和 β所有物性都收进 Ra 和 Pr。这就是为什么搜到的 MATLAB 程序大概率没有输入空气的导热系数也不需要任何单位换算。2.2 无量纲数先定网格和迭代策略计算前先把这几组无量纲数的角色搞清楚否则网格和松弛因子只能盲试。下表是计算中最常用的四组数无量纲数定义在程序中的作用Prν/α决定显式时间步上限Pr 越大扩散步越受限RagβΔTL³/(να)决定浮力强度和壁面热边界层厚度GrRa/Pr判断自然对流起始的辅助参数NuhL/k程序输出量最终要拿它和文献值比对Ra 的物理含义是浮力与粘性耗散的比值。Ra 升高壁面附近热边界层变薄厚度大致正比于 Ra^(-1/4)。41×41 网格计算 Ra1e4 够用Ra 提到 1e6 还用这套网格左壁热边界层里就只剩一个网格点温度梯度被严重低估Nu 一定会偏低。工程上判断网格够不够不用看流量残差直接对比两套网格算出的平均 Nu 差即可。2.3 为什么用流函数-涡量法而不是 SIMPLER很多从 Fluent 转过来的用户第一反应是找压力-速度耦合算法。二维方腔自然对流这个经典对标算例里流函数-涡量比 SIMPLER 合适得多原因有三条。第一边界简单。方腔四周都是固定壁面ψ 在边界恒为 0不需要进出口压力边界。第二未知量少。SIMPLER 每步要解 u、v、p、T 四个场ψ-ω 解法只要解三个流函数泊松方程用 SOR 迭代控制在 60 次以内即可。第三涡量边界条件有明确的物理含义壁面剪切层的强度直接从 ψ 的二阶导数得到。代价是三维问题里没有统一的流函数ψ-ω 法在三维很难扩展。所以这个方案锁在二维热对流计算、换热分析这类场景属于验证级工具而非生产级 CFD。选型时记住做方案对比、课程设计、热分析校核用它真要预测复杂几何的流场再回到 Fluent 也不迟。3. 用 MATLAB 算二维方腔自然对流的可运行代码3.1 参数设置和网格初始化的最小写法二维方腔自然对流是热对流数值计算领域最经典的验证题左壁高温 Th1右壁低温 Tc0上下绝热腔内流体是空气Pr0.71。计算域是边长为 1 的正方形均匀网格。网格从 41×41 起步这个规模在 MATLAB 里每次主迭代只需要几秒。% 无量纲控制参数 Ra 1e4; % 瑞利数先取经典对标值 Pr 0.71; % 空气普朗特数 nx 41; ny 41; % 网格节点数 L 1.0; h L/(nx-1); % 方腔边长和网格步长 dt 0.5*h^2/(1Pr); % 显式稳定时间步 % 流函数、涡量、温度场 psi zeros(ny, nx); % 流函数四周壁面恒为 0 omega zeros(ny, nx); % 涡量 T zeros(ny, nx); T(:,1) 1.0; % 左壁高温 T(:,nx) 0.0; % 右壁低温dt 的写法是有讲究的。显式格式扩散项的时间步上限是 Δt ≤ 0.5h²/max(1,Pr)写成 0.5*h^2/(1Pr) 比随便给一个 1e-3 更稳。后面涡量方程多一个浮力源项这个 dt 对 Ra1e4 是安全的。注意在这个无量纲程序里不要引入真实空气的比热容和导热系数出现这些量通常意味着单位换算已经出错。3.2 温度场、涡量场、流函数的主迭代骨架稳态解通过时间推进获得。每个时间步做四件事解流函数泊松方程更新壁面涡量显式推进内部温度与涡量最后复位温度边界。顺序固定后整个主循环就是这个骨架maxIter 50000; wSOR 1.6; % 流函数方程的松弛因子 for iter 1:maxIter % 第一步用当前涡量解流函数泊松方程 % 采用逐点 SOR内迭代 60 次已经足够 for sub 1:60 for j 2:ny-1 for i 2:nx-1 psi_s 0.25*(psi(j,i1)psi(j,i-1)psi(j1,i)psi(j-1,i) ... h^2*omega(j,i)); psi(j,i) (1-wSOR)*psi(j,i) wSOR*psi_s; end end end % 第二步用新 psi 更新壁面涡量 for j 2:ny-1 omega(j,1) -2*psi(j,2)/h^2; % 左壁 omega(j,nx) -2*psi(j,nx-1)/h^2; % 右壁 end for i 2:nx-1 omega(1,i) -2*psi(2,i)/h^2; % 下壁 omega(ny,i) -2*psi(ny-1,i)/h^2; % 上壁 end % 第三步显式推进内部温度与涡量 [T, omega] updateConvection(T, omega, psi, h, dt, Ra, Pr); % 第四步温度边界复位绝热壁用零梯度 T(1,:) T(2,:); % 下壁绝热 T(ny,:) T(ny-1,:); % 上壁绝热 T(:,1) 1.0; T(:,nx) 0; end壁面涡量那两行是很多初学脚本算不收敛的重灾区。以左壁为例边界 ψ0靠近壁面的内点 ψ(j,2) 做泰勒展开可得 ω_b -∂²ψ/∂x²|_{wall} ≈ -2ψ(j,2)/h²。如果这里写成 omega(j,1)0相当于把壁面当成了自由滑移边界方腔里根本不会形成旋转涡胞Nu 会算出纯导热的 1.0。SOR 松弛因子取 1.6 是因为流函数方程是正定泊松方程超松弛能明显减少子迭代次数温度场和涡量场是显式推进别再对它们做松弛。3.3 内部节点温度与涡量的离散格式内部节点更新是传热计算精度的主要来源。扩散项用中心差分对流项用一阶迎风保证系数矩阵对角占优。这段写成独立函数方便之后换二阶迎风或添加内热源。function [T, omega] updateConvection(T, omega, psi, h, dt, Ra, Pr) [ny, nx] size(T); T_new T; omega_new omega; for j 2:ny-1 for i 2:nx-1 u (psi(j1,i) - psi(j-1,i))/(2*h); v -(psi(j,i1) - psi(j,i-1))/(2*h); % 能量方程对流通量一阶迎风 if u 0 convx_T u*(T(j,i)-T(j,i-1))/h; else convx_T u*(T(j,i1)-T(j,i))/h; end if v 0 convy_T v*(T(j,i)-T(j-1,i))/h; else convy_T v*(T(j1,i)-T(j,i))/h; end % 扩散项中心差分 lapT (T(j,i1)T(j,i-1)T(j1,i)T(j-1,i) ... - 4*T(j,i))/h^2; T_new(j,i) T(j,i) dt*(lapT - convx_T - convy_T); % 涡量输运多一项浮力源项 Ra*Pr*dT/dx if u 0 convx_w u*(omega(j,i)-omega(j,i-1))/h; else convx_w u*(omega(j,i1)-omega(j,i))/h; end if v 0 convy_w v*(omega(j,i)-omega(j-1,i))/h; else convy_w v*(omega(j1,i)-omega(j,i))/h; end lapOmega (omega(j,i1)omega(j,i-1)omega(j1,i) ... omega(j-1,i) - 4*omega(j,i))/h^2; dTdx (T(j,i1)-T(j,i-1))/(2*h); omega_new(j,i) omega(j,i) ... dt*(Pr*lapOmega - convx_w - convy_w ... Ra*Pr*dTdx); end end T T_new; omega omega_new; end浮力项 dTdx 的符号决定流动方向。左壁热、右壁冷腔内左侧流体受热上升温度沿 x 方向下降∂T/∂x0产生正涡量最终形成逆时针主涡。如果符号写反涡的方向会反转温度场仍然能迭代出结果但流线显示为顺时针回流这个特征在可视化阶段一眼就能抓出来。一阶迎风的数值粘性在 Ra1e4 时对平均 Nu 的影响大约在 2%~5%验证性分析完全可接受程序跑通后再换 QUICK 或二阶迎风也不迟。4. 对流换热迭代的收敛参数和发散排查4.1 显式格式的三个必调参数显式时间推进的程序里三个参数直接决定能不能收敛。第一个是时间步 dt按 0.5h²/(1Pr) 取。超过这个值扩散项立即出现交替方向的高频振荡表现为残差先降后跳温度场出现棋盘状分布。第二个是流函数方程的 SOR 松弛因子取 1.3~1.6。取 1.8 以上在网格数增大后可能直接发散取 1.0 只能保证收敛但速度慢一倍。第三个是收敛判据不要只看流函数残差Nusselt 数连续多步不变才算稳态。下表是这三个参数的常用区间和典型症状参数常见区间超出区间的典型表现时间步 dt≤ 0.5h²/(1Pr)温度场棋盘振荡NaNSOR 松弛因子1.3~1.6大于 1.8 时 psi 迭代发散Nu 收敛门限1e-5 ~ 1e-6门限过大时结果未达稳态收敛检查我一般放在主循环中间每 500 步看一次平均 Nu 的变化% 前 10000 步跳过之后计算平均 Nu % 平均 Nu 计算函数在后文 5.1 中给出 if mod(iter, 500) 0 iter 10000 Nu_avg_current calcNuLocal(T, h); Nu_avg_current mean(Nu_avg_current(2:end-1)); if abs(Nu_avg_current - Nu_avg_prev) 1e-5 break; end Nu_avg_prev Nu_avg_current; end提示用 Nu 的步间差值做收敛判据能避开一个经典坑——流函数场已经不再变化但温度场还在非常缓慢地向冷壁传热。这种情况在 Ra 较高时很常见。4.2 角点温度冲突与绝热壁的常见写法方腔四个角点是边界条件的冲突区。左壁定温 T1、下壁绝热左下角同时属于这两条边。按绝热处理角点温度会被拖向内部温度造成左壁底部温度梯度偏小按左壁处理绝热壁的零梯度条件又被破坏。平均 Nu 对左壁全高度积分时角点误差会被带进去。常见做法是对角点做折中T(1,1) 0.5*(T(2,1) T(1,2)); % 左下角折中 T(1,nx) 0.5*(T(2,nx) T(1,nx-1)); % 右下角折中 T(ny,1) 0.5*(T(ny-1,1) T(ny,2)); T(ny,nx) 0.5*(T(ny-1,nx) T(ny,nx-1));另一个常见错误是绝热壁写成 T(1,:)T(1,:)。这一行在循环里永远不改变程序也不报错但下壁温度被冻结在初始值 0 上下半腔的浮力驱动就消失了。正确的零梯度写法是让边界节点等于相邻内点即主循环里的 T(1,:)T(2,:)。这两种写法的差别在收敛后看下壁附近等温线会非常明显错误写法下等温线会垂直扎进下壁。4.3 高瑞利数发散与数值粘性的取舍Ra 升到 1e5 以上一阶迎风加 41×41 网格的组合会出现两个症状一是近壁速度峰值被磨平二是主涡中心区出现非物理的波动。这两个症状都来自数值粘性而非浮力模型错误。应对策略是网格翻倍到 81×81同时把时间步按 h² 关系缩小。不要单独靠调小 dt 继续跑旧网格那样只会把已经被污染的温度场算得更光滑Nu 仍然偏低只是不再发散而已。Ra 到 1e6 量级时壁面热边界层厚度大约 0.02L41 个节点完全无法解析显式时间推进也开始力不从心推荐换成交替方向隐式格式处理扩散项。另外注意强制对流和自然对流的参数不能混用。用 Fluent 算管流时入口 Re5000 是强制对流主导方腔自然对流没有来流速度控制参数是 Ra。拿 Re 去估算这类问题的无量纲时间步会得到完全错误的量级这是从热分析软件切到 MATLAB 脚本时最容易被忽略的坑。5. 用 Nusselt 数验证热分析结果并导出 MATLAB 全场数据5.1 计算左壁局部与平均 Nusselt 数换热强度的最终输出是 Nusselt 数。左壁是高温壁局部 Nu 定义是 Nu_local -∂T/∂x|_{wall}无量纲温度下不需要导热系数和特征长度。用温度场向腔内取前向差分function Nu_local calcNuLocal(T, h) [ny, nx] size(T); Nu_local zeros(ny, 1); for j 2:ny-1 % 左壁面温度梯度用内侧一点做前向差分 Nu_local(j) -(T(j,1)-T(j,2))/h; end end Nu_local calcNuLocal(T, h); Nu_avg mean(Nu_local(2:end-1,1)); % 平均 Nusselt 数这个差分方向与壁面涡量一致都是朝计算域内部取点。Ra1e4、Pr0.71 时41×41 网格算出的平均 Nu 应在 2.2~2.3 之间对应经典基准解 2.24 附近。如果接近 1.0温度场接近纯导热优先查浮力项符号和壁面涡量如果明显大于 3多半是网格太粗导致壁面梯度被放大。5.2 温度场与流线同时画验证热羽流形态温度场和流函数场同时画出来是判断热对流是否真正开启的最直接方法。温度场用颜色云图速度场用流线叠加。x 0:h:1; y 0:h:1; figure; contourf(x, y, T, 20); colorbar; colormap(hot); hold on; % 流函数求速度场 u zeros(ny,nx); v zeros(ny,nx); for j 2:ny-1 for i 2:nx-1 u(j,i) (psi(j1,i)-psi(j-1,i))/(2*h); v(j,i) -(psi(j,i1)-psi(j,i-1))/(2*h); end end streamslice(x, y, u, v, 2);健康的温度场左壁根部有一个明显的高温薄层顶部向右侧弯折右侧冷壁底部有对称的低温薄层流线应在方腔中心形成一个逆时针大涡。如果等温线平直、没有边界层弯曲直接回主循环检查涡量边界条件。streamslice 的密度参数设为 2流线太密反而看不清主涡结构。工程汇报时温度云图加流线图就是热分析报告里需要的那两张图。5.3 把温度场和壁面 Nu 导出为 CSV验证通过后把完整温度场和左壁局部 Nu 导出来后续可以交给 Python 的 pandas 或 Excel 做进一步处理。用 writematrix 两行写完T_export [0, y; x, T]; % 首行和首列存放坐标 writematrix(T_export, cavity_T_field.csv); Nu_export [y(2:ny-1), Nu_local(2:ny-1)]; writematrix(Nu_export, cavity_Nu_left.csv);导出前确认主循环已经按 5.1 的收敛门限退出否则存下来的只是中间帧。快速校验的另一个技巧是改变 Ra 再算一组Ra1e3 与 Ra1e4 的平均 Nu 比值应接近 (1e4/1e3)^0.25这个指数对应层流自然对流关联式 Nu C·Ra^n 中的 n偏离超过 10% 就回头检查边界条件处理。跑出这部分结果后把温度场 CSV 与 Nu 分布 CSV 一起归档整套 MATLAB 对流换热脚本就能独立承担换热分析校核任务并和 Ansys 瞬态热分析结果互相印证。本文还有配套的精品资源点击获取