MATLAB实现对流换热数值计算：方腔自然对流模拟

简介面向工程热物理、航空航天、能源系统等领域研究者的MATLAB对流换热数值计算入门项目以有限体积法为核心系统讲解纳维-斯托克斯方程与能量方程的求解思路适合需要开展换热模拟、温度场分析的中高级学习者。资源压缩包共3个文件PDF说明文档梳理理论背景Word计算说明书详解建模与验证过程.m脚本为可直接运行的MATLAB源码整体仅746KB便于下载与快速部署。目前已有358人学习下载。内容覆盖对流换热理论、数学模型构建、有限体积离散、Dirichlet/Neumann/Robin边界条件施加等关键环节代码可按步骤复现并可视化温度场与速度矢量帮助读者从方程推导走向数值实验验证。通过该案例既能掌握完整换热数值计算流程也可提升MATLAB矩阵运算与迭代求解的编程实践能力对撰写课程报告或完成工程项目具有直接参考价值。1. 对流换热数值计算用 MATLAB 把热分析跑成可复现脚本对流换热数值计算的第一步很多工程师会直接开一套商业 CFD 软件但在电子散热、电池热管理和小尺寸风道这类二维简化模型里MATLAB 反而是更快的路径物理假设明确、参数可批量扫描、算法对审查完全透明。这个压缩包对应的就是一套标准解法——用有限差分网格求解二维不可压缩对流换热的温度场和速度场输出壁面换热系数与平均努塞尔数解决手算传热只能给平均值、给不出局部温度和流线分布的问题。适合计算传热学的学生、刚转热分析方向的仿真工程师以及要把传热学方程固化成公司内部脚本的研发人员。2. 对流换热控制方程与无量纲化N-S 方程怎么变成代码参数2.1 二维层流对流换热的四组控制方程写求解器之前先把物理模型钉死在纸面上。二维不可压缩层流对流换热控制方程由连续方程、x 方向动量方程、y 方向动量方程和能量方程组成。连续方程是∂u/∂x ∂v/∂y 0x 方向动量方程∂u/∂t u∂u/∂x v∂u/∂y −(1/ρ)∂p/∂x ν(∂²u/∂x² ∂²u/∂y²)y 方向动量方程∂v/∂t u∂v/∂x v∂v/∂y −(1/ρ)∂p/∂y ν(∂²v/∂x² ∂²v/∂y²) gβ(T−T_c)能量方程∂T/∂t u∂T/∂x v∂T/∂y α(∂²T/∂x² ∂²T/∂y²)这是原始变量形式u、v、p、T 直接作为未知量。y 方向动量方程里最后一项 gβ(T−Tc) 是浮力项也是自然对流的驱动源强迫对流场景去掉这一项换成入口速度边界条件即可。ν 是运动粘度α 是热扩散率β 是体积膨胀系数。方程本身不稀奇工程里真正容易踩坑的是浮力项的处理方式。如果温差大、物性变化超过 10%Boussinesq 假设就不够用了那时需要改成物性随温度变化的分段线性模型方程的刚度和数值稳定性都会明显变差。2.2 Boussinesq 假设的适用边界Boussinesq 假设的出发点是密度变化只在浮力项里出现其他位置按常数处理等价于把动量方程里的密度波动压缩到浮力项 gβ(T−Tc) 中。空气在 300 K 附近、温差 30 K 的工况这个假设的误差在 5% 以内温差超过 100 K或者介质处于超临界状态就不建议继续用。判断指标可以量化工程上一般要求 βΔT 0.1。β 对理想气体就是 1/T所以 300 K 环境下 ΔT 不应超过 30 K。超过这个范围要么换变物性模型要么在代码里把密度和导热系数按温度查表更新。很多从教材抄来的代码默认用户工况满足假设实际项目里温差一大结果直接偏掉。2.3 无量纲特征数Pr、Ra、Gr 与 Nu 的工程含义代码里不能直接带温度、速度算那样对流项和扩散项的数量级差距太大会让迭代极难收敛。需要先无量纲化这是热分析代码最关键的预处理也是最容易做错的地方。取特征长度 L、特征温差 ΔT Th − Tc参考速度取 α/L无量纲量定义为X x/L Y y/L U uL/α V vL/α τ tα/L² Θ (T − Tc)/(Th − Tc)无量纲化之后三个特征数自然浮现特征数定义物理含义Pr 普朗特数ν/α动量扩散与热量扩散能力之比空气约 0.71水约 7Ra 瑞利数gβΔTL³/(να)浮力驱动对流与耗散和热扩散的竞争Nu 努塞尔数hL/k对流换热强度与纯导热的比值Gr 与 Ra 的关系是 Gr Ra/Pr。电子散热报告习惯用 Ra暖通教材习惯用 Gr代码里应统一用 Ra输出不再乘 Pr。2.4 边界条件的无量纲写法以方腔自然对流为例左壁面高温 Th右壁面低温 Tc上下壁面绝热。无量纲边界条件写成左壁Θ 1u v 0右壁Θ 0u v 0上下壁∂Θ/∂n 0u v 0壁面涡量不能直接给数值要根据流函数反算。左壁的 ω_left −(8ψ₂ − ψ₃)/(2h²)这个式子的推导体现在第 3 章离散过程里。3. 有限差分离散与 MATLAB 求解流程从网格到压力修正3.1 原始变量法与涡量-流函数法的取舍MATLAB 里求解二维对流换热有原始变量法和涡量-流函数法两条路径。原始变量法直接解 u、v、p、T 四个未知场难点在压力-速度耦合要用 SIMPLE 算法或投影法反复迭代压力。涡量-流函数法消去压力变量未知量缩减为涡量 ω、流函数 ψ 和温度 T方程从四个变三个。工程选择的逻辑很明确如果只做二维自然对流和强迫对流的教学验证涡量-流函数法在 MATLAB 里最容易收敛、代码最短如果最终目标是三维模型或者湍流过渡原始变量法不可绕开。多数课程毕设和产品初版的热分析脚本涡量-流函数法够用。对热分析本身来说温度场其实是最容易求的难在流场。涡量-流函数法的求解链条包含三个环节涡量输运方程、流函数泊松方程、温度输运方程。三者每个时间步交替更新任何一步的离散符号错误整个温度场都会一起出错。3.2 空间离散格式中心差分与迎风差分二阶导数项的扩散项直接采用中心差分∂²φ/∂x² ≈ (φ(i1,j) − 2φ(i,j) φ(i−1,j)) / h²一阶对流项在教学代码里默认用中心差分u∂φ/∂x ≈ u(i,j)·(φ(i,j1) − φ(i,j−1)) / (2h)但网格稀或流速大时中心差分会产生非物理振荡。原因是它的截断误差里藏着一个负扩散项。把对流项换成迎风差分按速度方向选取上游节点u∂φ/∂x ≈ max(u,0)·(φ(i,j) − φ(i,j−1))/h min(u,0)·(φ(i,j1) − φ(i,j))/h迎风差分只有一阶精度数值扩散明显但保证解的单调性。调试阶段先用迎风差分把流场跑稳定再换回中心差分提升精度是我常用的顺序。提示涡量-流函数法解出的压力场只能是后处理推导值。如果后续衔接结构热应力分析需要精确压力场建议尽早切换到原始变量法。3.3 涡量输运方程与流函数泊松方程的离散方腔自然对流的无量纲涡量输运方程为∂ω/∂τ U∂ω/∂X V∂ω/∂Y Pr(∂²ω/∂X² ∂²ω/∂Y²) Pr·Ra·∂Θ/∂X时间项用前向欧拉显式格式。每步内顺序固定先更新涡量和温度再由涡量求流函数最后反算速度。流函数泊松方程用 SOR 逐次超松弛求解SOR 格式为ψ(i,j) (1−γ)·ψ(i,j) γ/4·[ψ(i1,j) ψ(i−1,j) ψ(i,j1) ψ(i,j−1) h²·ω(i,j)]松弛因子 γ 取 1.51.8。瑞利数越高γ 越要往 1.5 靠强浮力环境下过大的松弛因子直接引发震荡。3.4 时间步长选取与收敛判据显式推进格式的稳定性由扩散项决定。方腔自然对流的动量扩散系数是 Pr热扩散系数是 1。二维显式格式的稳定条件近似为Δτ ≤ h² / (2·(Pr 1))以 81×81 网格、Pr 0.71 为例h 0.0125Δτ 约为 4.5e-5。这个步长下收敛稳态需要上万次迭代MATLAB 嵌套循环会比较慢。工程上常见做法是在核心循环内尽量向量化外部迭代保持不变。实际代码里我给 dt 乘一个 0.10.5 的安全系数每 200 步监控一次残差。收敛判据看温度场和涡量场的相对变化量同时小于 1e-6这比单看温度场要严格因为温度收敛不代表流场已经稳定。4. 方腔自然对流 MATLAB 算例代码、参数与热分析结果提取4.1 可运行的完整 MATLAB 脚本下面这段是涡量-流函数法求解方腔自然对流的经典教学实现不依赖任何工具箱纯 MATLAB 脚本直接运行。% 方腔自然对流涡量-流函数法 % Pr 默认 0.71空气Ra 可调 clear; clc; N 81; % 网格数 N x N L 1.0; h L/(N-1); % 网格间距 Pr 0.71; Ra 1e4; % 瑞利数 T zeros(N,N); T(:,1) 1.0; % 左壁无量纲温度 1 T(:,N) 0.0; % 右壁无量纲温度 0 psi zeros(N,N); % 流函数边界取 0 omega zeros(N,N); % 涡量 gamma 1.7; % SOR 松弛因子 dt 0.2 * h^2 / (1 Pr); % 显式差分稳定时间步 maxiter 20000; tol 1e-6; for step 1:maxiter T_old T; omega_old omega; omega_new omega; T_new T; % 1) 内部点更新涡量和温度 for i 2:N-1 for j 2:N-1 % 由流函数差分离散计算速度 u (psi(i,j1) - psi(i,j-1)) / (2*h); % u dpsi/dy v -(psi(i1,j) - psi(i-1,j)) / (2*h); % v -dpsi/dx % 涡量拉普拉斯项 lap_w (omega(i1,j) omega(i-1,j) ... omega(i,j1) omega(i,j-1) - 4*omega(i,j)) / h^2; % 涡量对流项中心差分 conv_w u*(omega(i,j1)-omega(i,j-1))/(2*h) ... v*(omega(i1,j)-omega(i-1,j))/(2*h); % 浮力源项Pr*Ra*dTheta/dX buoy Pr*Ra*(T(i1,j)-T(i-1,j))/(2*h); omega_new(i,j) omega(i,j) dt*(Pr*lap_w - conv_w buoy); % 能量方程 lap_T (T(i1,j)T(i-1,j)T(i,j1)T(i,j-1)-4*T(i,j))/h^2; conv_T u*(T(i,j1)-T(i,j-1))/(2*h) v*(T(i1,j)-T(i-1,j))/(2*h); T_new(i,j) T(i,j) dt*(lap_T - conv_T); end end omega omega_new; T T_new; % 2) 上下壁绝热温度零梯度镜像 for j 2:N-1 T(1,j) T(2,j); T(N,j) T(N-1,j); end % 3) SOR 求解流函数泊松方程 for iter_sor 1:30 for i 2:N-1 for j 2:N-1 psi(i,j) (1-gamma)*psi(i,j) gamma/4 * ... (psi(i1,j)psi(i-1,j)psi(i,j1)psi(i,j-1)h^2*omega(i,j)); end end end % 4) 用新流函数更新壁面涡量 for i 2:N-1 omega(i,1) -(8*psi(i,2) - psi(i,3)) / (2*h^2); omega(i,N) -(8*psi(i,N-1) - psi(i,N-2)) / (2*h^2); end for j 2:N-1 omega(1,j) -(8*psi(2,j) - psi(3,j)) / (2*h^2); omega(N,j) -(8*psi(N-1,j) - psi(N-2,j)) / (2*h^2); end % 5) 周期性残差监控 if mod(step, 200) 0 err_T max(max(abs(T - T_old))); err_w max(max(abs(omega - omega_old))); if err_T tol err_w tol fprintf(收敛于第 %d 步err_T%.2eerr_w%.2e\n, ... step, err_T, err_w); break; end end end代码整体逻辑分五段内部点显式推涡量和温度、绝热边界镜像、SOR 解流函数、由流函数更新壁面涡量、周期性残差判断。u 的符号取决于 ψ 对 Y 的偏导v 是 ψ 对 X 偏导的负值这个约定必须和浮力源项里的 ∂Θ/∂X 一致否则流场会逆向旋转。4.2 代码的关键参数说明调参时优先动下面几个量注意它们之间互相牵连参数典型取值调参注意Ra1e31e6超过 1e5 时 81×81 网格明显不够需加密并减小 dtPr0.71空气、7水换介质后 dt 要按 Pr 重算N41、81、161网格每加密两倍计算量增约 8 倍gamma1.51.8高 Ra 时降到 1.5 以下防振荡dth²/(2·(Pr1)) 再乘安全系数不建议加大显式格式没有精度换步长空间壁面涡量更新公式是最容易抄错的位置。以左壁为例从流函数在壁面附近做泰勒展开并代入无滑移条件 ψ_壁 0、∂ψ/∂n 0得到 ω_wall −∂²ψ/∂n² ≈ −(8ψ₂ − ψ₃)/(2h²)。注意 ψ₂ 是距离壁面第一个网格点不是第二个。如果用了原始变量法和投影法壁面涡量不需要显式计算。4.3 平均努塞尔数与温度场可视化收敛后提取工程报告用到的关键指标平均努塞尔数。左壁热壁的无量纲温度梯度为 −∂Θ/∂X用一阶差分近似% 热壁面局部努塞尔数i 2:N-1 避开角点 Nu_local zeros(N,1); for i 2:N-1 Nu_local(i) -(T(i,2) - T(i,1)) / h; end Nu_avg mean(Nu_local(2:N-1)); fprintf(平均努塞尔数 Nu %.4f\n, Nu_avg); % 温度场云图与流函数等值线 figure; contourf(linspace(0,1,N), linspace(0,1,N), T, 40); colorbar; title([温度场 Ra, num2str(Ra)]); figure; contour(linspace(0,1,N), linspace(0,1,N), psi, 30); title(流函数等值线);Ra 1e4、Pr 0.71、81×81 网格的经典参考值是 Nu ≈ 2.24。跑出来落在 2.22.3 区间说明代码正确。局部努塞尔数在热壁面底部最大顶部最小绘制成曲线能看到边界层发展的全过程。4.4 最常见的不收敛与错误模式三个高频问题按出现概率排序。第一个是流场方向反转。涡量和浮力项的符号配合出错直接表现是左壁附近流线向下走而不是向上。检查浮力项 Pr·Ra·(T(i1,j)−T(i−1,j))/(2h) 的差分方向i 增加的方向必须与重力作用方向相反。第二个是温度场震荡发散。通常不是时间步长问题而是对流项中心差分在高 Ra 下产生波动。把能量方程的对流项换成迎风差分验证稳定则说明需要加密网格。第三个是流函数突然出现 NaN。壁面涡量公式符号反了或 SOR 松弛因子大于 1.8 是常见原因逐项检查即可。5. 网格无关性验证与 Nu 结果的工程校准5.1 网格无关性的验证操作网格无关性验证不是可选项是热分析报告里审计必查的内容。具体操作是同一 Ra、同一 Pr 下分别以 N41、81、161 跑三遍记录平均 Nu。相邻网格的 Nu 相对偏差小于 1% 视为达到网格无关。% 将第 4 章脚本封装成 run_cavity(Ra, N) 后批量扫描 Ra_set [1e3, 1e4, 1e5]; N_set [41, 81, 161]; Nu_table zeros(3,3); for k 1:3 for n 1:3 Nu_table(k,n) run_cavity(Ra_set(k), N_set(n)); fprintf(Ra%.0e N%d Nu%.4f\n, ... Ra_set(k), N_set(n), Nu_table(k,n)); end end工程中更隐蔽的问题是Ra 提高后网格不变收敛也没有报错但 Nu 比文献值偏低超过 10%。原因是对流项的数值扩散在稀网格上把壁面温度梯度磨平了。加密后若 Nu 上升说明原解偏保守设计余量被算大了。5.2 用文献参考值校准你的求解器热分析代码本身不产生结论结论要放在基准解上校准。方腔自然对流在 Pr 0.71 下的经典基准数据de Vahl Davis可以直接做回归测试Ra参考平均 Nu可接受误差区间1e31.118±1%1e42.243±1%1e54.519±2%1e68.800±3%161×161 网格在 Ra 1e5 时 Nu 低于 4.3先查近壁边界层的网格分辨率再查壁面涡量公式系数。二维方腔问题没有隐藏变量差值超出误差区间就一定是数值格式的问题。校准通过之后再把边界换成实际产品的热源加载方式做进一步的瞬态热分析推演。5.3 从温度场提取工程用的换热系数把无量纲 Nu 还原成有量纲的对流换热系数 h公式为h Nu·k / Lk 是流体导热系数L 是特征长度。举例空气 k 0.026 W/(m·K)方腔边长 0.1 mRa1e4 时 Nu2.24得到 h ≈ 0.58 W/(m²·K)。这就是设计报告里要配的那个关键参数。MATLAB 换算时注意单位统一。长度用米、温度用开尔文k 的单位必须是 W/(m·K)。工程物性表给的若是摄氏温度下的导热系数按 273.15 K 换算后再代入。导出温度场时用 40 层等高线图观察近壁区的温度梯度分布若热壁面附近的等温线间距明显大于冷壁面说明网格在边界层方向分布不足加密时优先加密近壁区而不是全场均匀加密。本文还有配套的精品资源点击获取