简介面向流体力学、传热学与化学工程等领域的Matlab有限体积法求解器用于瞬态对流扩散偏微分方程的数值模拟。资源以有限体积法为核心提供主程序与启动模块、网格生成、物理属性设置等完整骨架并配有PDE图示和边界条件说明可在Matlab中直接运行和修改参数帮助理解守恒离散、通量计算等数值核心。压缩包共343个文件以m源码为主300个另有png示意图、md文档、mlx实时脚本、mat数据及ipynb对照实验并附IAPWS-IF97物性函数涵盖球坐标一维扩散解析解对比、顶盖驱动腔流、相变焓法等多个经典算例整体仅999KB结构清晰、便于按需检索代码注释与示例文档也方便二次开发。已有156人学习浏览适合需要快速搭建PDE求解框架或系统学习有限体积法的研究生、工程师与科研爱好者。1. 有限体积求解器并没有那么难FVM 解瞬态对流扩散的清晰骨架同样是求解瞬态对流扩散偏微分方程我在 MATLAB 里从零写过有限差分、也封装过有限体积。最后留下来反复用的反而是这套不到 200 行核心逻辑的 FVM 求解器。它的主语不是“方法论文档”而是 FVTdemo.m 这种能直接跑出 diff_pde.jpg 和 diff_pde_3d.jpg 的演示脚本还有一个 IAPWS_IF97.m 帮你把水蒸气物性算进去。比起把每个偏微分方程都写成专用程序这套解算器只做一件事把守恒形式离散到控制体上剩下的系数、边界、时间推进全部作为参数暴露出来。适合想弄懂偏微分方程离散过程、又不想一开始就抱大型工具箱的 MATLAB 使用者。2. FVM 离散从守恒方程到三对角矩阵2.1 为什么要把方程改写成守恒积分式有限体积法处理瞬态对流扩散问题的第一步是把偏微分方程改写成控制体上的积分守恒形式∂/∂t ∫V ρφ dV ∮S (ρuφ − Γ∇φ)·n dS ∫V S dV左边第二项是通过控制体表面的净通量其中 ρuφ 是对流通量−Γ∇φ 是扩散通量。这个形式的价值在于散度定理被用在了离散层面相邻控制体共享同一个界面通量从一侧流出就必然从另一侧流入因此全局守恒是自动满足的。我在用有限差分处理强对流问题时经常遇到数值振荡而写成这种通量平衡后即使网格粗糙也没有“虚拟源项”出现。相比有限元法要处理形函数和弱形式有限差分法用点值直接做差分近似FVM 的优势在于通量可以从控制体的真实几何边界上走出来。对传热或流体这类需要严格守恒的物理场FVM 不会像中心差分那样在间断附近产生振荡也不会像有限元法那样在细网格下调试迎风稳定化参数时让人头疼。这在处理强对流占主导的瞬态问题时尤其明显。2.2 一维均匀网格的离散系数组装以一个一维控制体 i 为例相邻控制体 w 和 e 布置在面上。采用迎风格式处理对流项扩散项用中心差分时间上用隐式欧拉离散后的代数方程为aP φP aE φE aW φW aP0 φP0 Su其中系数定义为aE De max(0, −Fe)aW Dw max(0, Fw)aP0 ρAΔx/Δt。F 表示面质量流量D 表示扩散传导率。如果流量从西向东F0迎风会保留西侧贡献反过来也一样。这样做的好处是系数矩阵保持对角占优隐式求解不会出现无界振荡。在 MATLAB 中定义偏微分方程时不必像符号工具那样显式写出方程名只需把对流速度、扩散系数和源项数组填好就行。下面的代码用spdiags直接组装一维问题的三对角矩阵后续只改Gamma和u就能切换纯扩散与对流扩散% 一维瞬态对流扩散 FVM 矩阵组装隐式欧拉 N 60; L 1.0; dx L/N; u 0.3; Gamma 0.01; rho 1.0; A 1.0; dt 0.001; % 时间步长决定 aP0 的大小 D Gamma*A/dx; % 扩散传导率 F rho*u*A; % 界面质量流量 % 迎风离散系数均匀网格 aW D max(0, F); % 西侧系数 aE D max(0, -F); % 东侧系数 aP0 rho*A*dx/dt; % 时间项系数 % 三对角矩阵主体 diag_main (aW aE aP0) * ones(N,1); diag_u -aE * ones(N-1,1); % 上对角线 diag_l -aW * ones(N-1,1); % 下对角线 M spdiags([diag_l, diag_main, diag_u], [-1 0 1], N, N); % 边界条件左、右为第一类Dirichlet M(1,:) 0; M(1,1) 1; M(N,:) 0; M(N,N) 1;这个矩阵M的每一行对应一个控制体的系数平衡。spdiags的第二个参数[-1 0 1]指定三条对角线的偏移量主对角线放aW aE aP0上对角线放-aE下对角线放-aW。显式边界处理之所以直接整行清零再置对角线为 1是为了让边界节点的值不参与内部通量计算。若换成第二类或第三类边界只需改写对应行的通量系数而不必动内部循环。2.3 时间推进方案的选择为什么默认推荐隐式表格列出三种常见时间格式它们会直接影响步长上限和耗散误差时间格式稳定条件精度特点显式欧拉CFL uΔt/Δx ≤ 1 且扩散数 ΓΔt/(Δx²) ≤ 0.5一阶每个时间步只做一次矩阵向量乘但步长受限隐式欧拉无条件稳定一阶三对角求解耗散略大Crank-Nicolson无条件稳定高振荡下可能产生伪振荡二阶对时间项做梯形平均适合平滑初值我在写 FVTdemo.m 内部的循环时默认走隐式欧拉因为瞬态对流扩散问题在时间方向上是刚性的网格加密之后显式格式的步长会被压到不可接受的程度。隐式格式每步都在解稀疏线性系统MATLAB 的三对角求解器开销很小多花的时间换来的是可以放心的dt。如果你只需要观察长时间行为甚至可以把dt放大一个数量级让时间项系数变小趋近稳态解。3. FVTdemo.m 实战从参数表到可视化3.1 先读发布文档再读代码资源里的 FVTdemo.html 不是手工写的而是 MATLAB 的publish功能生成的代码报告。浏览器打开后能在同一页面看到源码、输出图和运行注释。建议的运行顺序是先打开 FVTdemo.html 观察主程序的调用流程再用 MATLAB 的edit FVTdemo.m对照源码。很多初学 FVM 的人直接打开 .m 文件看到一长串系数组装就失去耐心其实那份 html 已经把“输入参数 - 网格生成 - 时间循环 - 绘图”的流水线标出来了。FVTdemo.m 的开头一般在定义物理参数然后调用meshgrid或linspace生成网格。拿这份代码改案例时我习惯把参数集中在脚本顶部像查表一样替换。表 3.1 是主程序里最常见的参数组参数名常见初始值含义调整建议Lx, Ly1.0计算域尺寸影响网格分辨率取值要配合 Nx/NyNx, Ny40控制体数量加密后必须复查 CFL 条件u0.2对流速度大于 1 可能让前沿变陡需缩短 dtGamma0.01扩散系数代表热扩散率、分子扩散系数dt0.0005时间步长隐式可放大但会损失瞬态精度nStep2000时间步总数由总模拟时间 T_end nStep*dt 决定bcLeftDirichlet左边界类型常用 Dirichlet/Neumann决定矩阵首行系数最容易混淆的是把Nx当作网格节点数。FVM 的控制体中心和节点不同Nx是控制体数量若边界条件是 Dirichlet实际存储数组长度仍是Nx因为边界值被强制写到首末控制体上。使用时注意网格坐标在中心点采样而不是边缘。3.2 时间循环和残差监视主程序推开核心循环后结构通常是预分配结果数组再按预设频率保存剖面。下面的框架可以直接嵌进自己的脚本% 从 FVTdemo.m 中抽出的核心时间循环框架 phi zeros(Nx, Ny); % 初始浓度或温度场 phi_store zeros(Nx, Ny, nStep/nOut 1); phi_store(:,:,1) phi; for k 1:nStep % 组装右端项时间项phi_old 源项 rhs aP0 * phi; % 上一时刻的 phi顺序为列向量 rhs apply_source(rhs, x, y); % 自定义源项可以写成匿名函数 % 求解稀疏线性系统M 来自上一节的系数矩阵 phi_vec M \ rhs(:); % 边界值覆盖第一类边界 phi_vec(idx_left) phiL; phi_vec(idx_right) phiR; % 重塑回二维网格 phi reshape(phi_vec, Nx, Ny); % 每 nOut 步保存一次方便后续画动画 if mod(k, nOut) 0 phi_store(:,:, k/nOut1) phi; end end这里有个细节M \ rhs用的是 MATLAB 内置稀疏直接解法。非稳态问题一般只需一次 LU 分解之后每个时间步只做一次回代因此不必在循环里反复分解。如果问题变成非线性比如源项依赖phi^2就得把M拆成可更新的部分在循环内部重算对角线并用迭代求解器bicgstab代替直接解。3.3 可视化与数据导出示例运行结束后FVTdemo 会生成 diff_pde.jpg 和 diff_pde_3d.jpg这正是 MATLAB 画图指令contourf和surf的输出。下面的代码把某个时刻的剖面导出为 CSV方便在 Python 或 Excel 里复核% 导出中心线剖面数据 xmid x; % 取网格中心坐标 phi_mid squeeze(phi_store(:, Ny/2, end)); % 最后时刻中心线 T_out table(xmid(:), phi_mid(:), VariableNames, {x, phi}); writetable(T_out, fvm_profiles.csv);CSV 文件名不要带中文MATLAB 的writetable在部分 Linux 版本下对中文路径处理不够稳定。squeeze用于去掉二维数组里的单一维度避免导出后列数不整齐。这个 CSV 可以直接用readtable(fvm_profiles.csv)读回也可以扔给任何第三方库做差分对比。4. 不止一维球坐标、焓法和顶盖驱动空腔的扩展用法4.1 球坐标扩散和 FVTool、FiPy 的解析解对比资源里的diffusion1Dspherical_analytic_vs_FVTool_vs_Fipy.ipynb是一个 Jupyter Notebook核心任务是把一维球形扩散的数值解和解析解放在一起比较。球坐标下的径向扩散方程为 ∂T/∂t α/r² ∂/∂r(r² ∂T/∂r)直接套用笛卡尔坐标系数会出错。需要把控制体积改为球壳体积面面积是 4πr² 而不是常量。% 球坐标系下的控制体积几何修正 dr 0.005; r dr/2:dr:1; N length(r); rmid r(:); % 左界面和右界面半径 rL rmid - dr/2; rR rmid dr/2; % 体积: 球壳 (4/3)π(rR^3 - rL^3) V (4*pi/3) .* (rR.^3 - rL.^3); % 界面面积 AL 4*pi .* rL.^2; AR 4*pi .* rR.^2; % 扩散传导率取左右界面的几何平均值 D alpha .* (AL AR) ./ 2;注意在这个离散式里如果保留源项请先把它乘以控制体体积在球坐标中边界处的面积为零首末控制体的几何修正对结果影响很大。我在写这个示例时常常忽略rL(1)0处的面积导致边界通量丢失发现后改为显式设置对称边界条件数值解才与解析解重合。解析解来自球坐标系分离变量解通常表现为无限级数。Jupyter Notebook 里用 FiPy 做交叉验证的好处是两边用的是完全不同的离散实现若数值解和解析解、FiPy 解的偏差都在 1e-3 以内可以确认自己的 MATLAB FVM 主程序没有结构性错误。比较时建议固定同一个 CFL别单纯比计算时间。4.2 焓法处理相变IAPWS_IF97 的调用位置相变问题不能用单一温度方程直接求解因为潜热会让焓曲线出现拐点。phaseChangeEnthalpyMethodExample.m用的是焓法把能量守恒写成含焓的形式每步先更新焓再由 h-T 映射表读出温度。IAPWS_IF97.m 在这里不是求解器而是物性查询函数供应传热计算所需的密度、比焓和比热。我自己常用的调用模式是先查压力对应的饱和温度再返回焓% 访问 IAPWS_IF97 的几个典型入口 p 1e5; % 1 bar Tsat IAPWS_IF97(Tsat_p, p); % 饱和温度 h_liq IAPWS_IF97(h_pT, p, Tsat-1); % 过冷液焓 h_vap IAPWS_IF97(h_pT, p, Tsat1); % 过热蒸汽焓 latent h_vap - h_liq; % 汽化潜热IAPWS_IF97的接口有多个入口Tsat_p输入压力输出饱和温度h_pT输入压力和温度输出焓值。用的时候要留意单位IF97 公式的压力基准为 Pa温度基准为 K焓基准为 J/kg有的版本会将 MPa 混入导致焓值差几个数量级。焓法循环里每步都调用物性函数会明显变慢我一般先建一个焓-温度查找表循环里只查表并线性插值。示例文件求解对象附加依赖典型收敛判据diffusion1Dspherical_analytic_vs_...球坐标一维扩散FiPy / FVToolmaxphaseChangeEnthalpyMethodExample.m一维相变 Stefan 问题IAPWS_IF97.m液固界面位置误差 2%SteadyLidDrivenCavityExample.m不可压稳态顶盖驱动空腔无中心速度残差 1e-6FVTdemo.m瞬态对流扩散无总浓度相对变化 1e-54.3 顶盖驱动空腔从标量输运到矢量流场SteadyLidDrivenCavityExample.m呈现的是不可压 Navier-Stokes 的基准算例。它和瞬态对流扩散共享 FVM 思想但多了压力-速度耦合常见做法是 SIMPLE 算法先假定压力场解动量方程得到速度场再修正压力让速度满足连续性方程。顶盖驱动空腔的边界条件很直观左、右、下边界为无滑移上边界给水平速度 U1。我在实际跑这个例子时只需要把上一步的标量扩散求解器嵌套进两个循环外层迭代压力修正内层解两个方向的动量方程。边界条件不再是简单的 Dirichlet而是速度分量的逐一指定。如果发现压力振荡通常是没有在网格上使用交错布置或压力参考点选得不对。这个示例的意义是告诉你FVM 的核心离散一旦写对从标量方程到矢量方程只是加装耦合器的过程。5. 收敛性检查与常见坑CFL、残差和 CSV 导出5.1 显示 CFL 检查和网格独立性瞬态对流扩散最容易出的问题是时间步长太大前沿出现“越界”振荡。即使在隐式格式里太大 dt 也会造成物理上的失真。我通常在每个脚本开头计算一次 CFL 数% 显式上限参考隐式可放宽但不建议盲目放大 cfl u * dt / dx; if cfl 1 warning(CFL%.2f 超过1建议减小 dt 或加密 dx, cfl); end这里 dx 取所有方向的最小网格间距u 取最大速度模。CFL 不满足时隐式格式虽然不会发散但数值解会把对流前沿抹得更平看起来像人为增大了扩散系数。网格独立性检查比收敛判据更实用把网格数从 40 翻倍到 80如果同一时刻的剖面最大变化小于 2%可以进入参数研究如果还在大幅移动则继续加密并同步调整 dt。5.2 质量守恒残差监视有限体积法最好的自检指标是总质量守恒。在每个时间步记录 sum(phi.*V) 随时间的变化理想情况下应保持一个常数或等于边界注入总量。下面的代码片段可以加到循环末尾% 统计总质量/总能量输出到命令行 total sum(phi(:) .* V(:)); % 一维时 V 是控制体积 fprintf(%4d total %.8f\n, k, total);若总质量在某一步突然跳变先去查边界条件如果边界是 Neumann应确认通量项没有在矩阵行里被意外清零。这个检查比盯着等值线图直观得多。5.3 把结果发给其他工具分析画完 3D 图后需要把几个时刻的数据连同坐标一起导出。用writetable导出 CSV 已经很常见如果你的同事用 Python 做后处理建议保存 numpy 能读的纯数字矩阵不要混入文本表头。% 保存二进制数据避免 CSV 在大网格下读写太慢 save(phi_final.mat, phi, x, dt, nStep);或者用dlmwrite(phi_final.csv, phi, precision, %.6e)只写数字矩阵。导出的 CSV 可直接用readmatrix读回 MATLAB这样后续做 FFT 仿真或频谱分析时不用重新跑一遍数值解。读入 Python 时用np.loadtxt完全避开 MATLAB 表格表头的解析问题。这也是我在多语言对照验证时最常用的导出方式。本文还有配套的精品资源点击获取 用MATLAB从零实现有限体积法：瞬态对流扩散求解骨架
简介面向流体力学、传热学与化学工程等领域的Matlab有限体积法求解器用于瞬态对流扩散偏微分方程的数值模拟。资源以有限体积法为核心提供主程序与启动模块、网格生成、物理属性设置等完整骨架并配有PDE图示和边界条件说明可在Matlab中直接运行和修改参数帮助理解守恒离散、通量计算等数值核心。压缩包共343个文件以m源码为主300个另有png示意图、md文档、mlx实时脚本、mat数据及ipynb对照实验并附IAPWS-IF97物性函数涵盖球坐标一维扩散解析解对比、顶盖驱动腔流、相变焓法等多个经典算例整体仅999KB结构清晰、便于按需检索代码注释与示例文档也方便二次开发。已有156人学习浏览适合需要快速搭建PDE求解框架或系统学习有限体积法的研究生、工程师与科研爱好者。1. 有限体积求解器并没有那么难FVM 解瞬态对流扩散的清晰骨架同样是求解瞬态对流扩散偏微分方程我在 MATLAB 里从零写过有限差分、也封装过有限体积。最后留下来反复用的反而是这套不到 200 行核心逻辑的 FVM 求解器。它的主语不是“方法论文档”而是 FVTdemo.m 这种能直接跑出 diff_pde.jpg 和 diff_pde_3d.jpg 的演示脚本还有一个 IAPWS_IF97.m 帮你把水蒸气物性算进去。比起把每个偏微分方程都写成专用程序这套解算器只做一件事把守恒形式离散到控制体上剩下的系数、边界、时间推进全部作为参数暴露出来。适合想弄懂偏微分方程离散过程、又不想一开始就抱大型工具箱的 MATLAB 使用者。2. FVM 离散从守恒方程到三对角矩阵2.1 为什么要把方程改写成守恒积分式有限体积法处理瞬态对流扩散问题的第一步是把偏微分方程改写成控制体上的积分守恒形式∂/∂t ∫V ρφ dV ∮S (ρuφ − Γ∇φ)·n dS ∫V S dV左边第二项是通过控制体表面的净通量其中 ρuφ 是对流通量−Γ∇φ 是扩散通量。这个形式的价值在于散度定理被用在了离散层面相邻控制体共享同一个界面通量从一侧流出就必然从另一侧流入因此全局守恒是自动满足的。我在用有限差分处理强对流问题时经常遇到数值振荡而写成这种通量平衡后即使网格粗糙也没有“虚拟源项”出现。相比有限元法要处理形函数和弱形式有限差分法用点值直接做差分近似FVM 的优势在于通量可以从控制体的真实几何边界上走出来。对传热或流体这类需要严格守恒的物理场FVM 不会像中心差分那样在间断附近产生振荡也不会像有限元法那样在细网格下调试迎风稳定化参数时让人头疼。这在处理强对流占主导的瞬态问题时尤其明显。2.2 一维均匀网格的离散系数组装以一个一维控制体 i 为例相邻控制体 w 和 e 布置在面上。采用迎风格式处理对流项扩散项用中心差分时间上用隐式欧拉离散后的代数方程为aP φP aE φE aW φW aP0 φP0 Su其中系数定义为aE De max(0, −Fe)aW Dw max(0, Fw)aP0 ρAΔx/Δt。F 表示面质量流量D 表示扩散传导率。如果流量从西向东F0迎风会保留西侧贡献反过来也一样。这样做的好处是系数矩阵保持对角占优隐式求解不会出现无界振荡。在 MATLAB 中定义偏微分方程时不必像符号工具那样显式写出方程名只需把对流速度、扩散系数和源项数组填好就行。下面的代码用spdiags直接组装一维问题的三对角矩阵后续只改Gamma和u就能切换纯扩散与对流扩散% 一维瞬态对流扩散 FVM 矩阵组装隐式欧拉 N 60; L 1.0; dx L/N; u 0.3; Gamma 0.01; rho 1.0; A 1.0; dt 0.001; % 时间步长决定 aP0 的大小 D Gamma*A/dx; % 扩散传导率 F rho*u*A; % 界面质量流量 % 迎风离散系数均匀网格 aW D max(0, F); % 西侧系数 aE D max(0, -F); % 东侧系数 aP0 rho*A*dx/dt; % 时间项系数 % 三对角矩阵主体 diag_main (aW aE aP0) * ones(N,1); diag_u -aE * ones(N-1,1); % 上对角线 diag_l -aW * ones(N-1,1); % 下对角线 M spdiags([diag_l, diag_main, diag_u], [-1 0 1], N, N); % 边界条件左、右为第一类Dirichlet M(1,:) 0; M(1,1) 1; M(N,:) 0; M(N,N) 1;这个矩阵M的每一行对应一个控制体的系数平衡。spdiags的第二个参数[-1 0 1]指定三条对角线的偏移量主对角线放aW aE aP0上对角线放-aE下对角线放-aW。显式边界处理之所以直接整行清零再置对角线为 1是为了让边界节点的值不参与内部通量计算。若换成第二类或第三类边界只需改写对应行的通量系数而不必动内部循环。2.3 时间推进方案的选择为什么默认推荐隐式表格列出三种常见时间格式它们会直接影响步长上限和耗散误差时间格式稳定条件精度特点显式欧拉CFL uΔt/Δx ≤ 1 且扩散数 ΓΔt/(Δx²) ≤ 0.5一阶每个时间步只做一次矩阵向量乘但步长受限隐式欧拉无条件稳定一阶三对角求解耗散略大Crank-Nicolson无条件稳定高振荡下可能产生伪振荡二阶对时间项做梯形平均适合平滑初值我在写 FVTdemo.m 内部的循环时默认走隐式欧拉因为瞬态对流扩散问题在时间方向上是刚性的网格加密之后显式格式的步长会被压到不可接受的程度。隐式格式每步都在解稀疏线性系统MATLAB 的三对角求解器开销很小多花的时间换来的是可以放心的dt。如果你只需要观察长时间行为甚至可以把dt放大一个数量级让时间项系数变小趋近稳态解。3. FVTdemo.m 实战从参数表到可视化3.1 先读发布文档再读代码资源里的 FVTdemo.html 不是手工写的而是 MATLAB 的publish功能生成的代码报告。浏览器打开后能在同一页面看到源码、输出图和运行注释。建议的运行顺序是先打开 FVTdemo.html 观察主程序的调用流程再用 MATLAB 的edit FVTdemo.m对照源码。很多初学 FVM 的人直接打开 .m 文件看到一长串系数组装就失去耐心其实那份 html 已经把“输入参数 - 网格生成 - 时间循环 - 绘图”的流水线标出来了。FVTdemo.m 的开头一般在定义物理参数然后调用meshgrid或linspace生成网格。拿这份代码改案例时我习惯把参数集中在脚本顶部像查表一样替换。表 3.1 是主程序里最常见的参数组参数名常见初始值含义调整建议Lx, Ly1.0计算域尺寸影响网格分辨率取值要配合 Nx/NyNx, Ny40控制体数量加密后必须复查 CFL 条件u0.2对流速度大于 1 可能让前沿变陡需缩短 dtGamma0.01扩散系数代表热扩散率、分子扩散系数dt0.0005时间步长隐式可放大但会损失瞬态精度nStep2000时间步总数由总模拟时间 T_end nStep*dt 决定bcLeftDirichlet左边界类型常用 Dirichlet/Neumann决定矩阵首行系数最容易混淆的是把Nx当作网格节点数。FVM 的控制体中心和节点不同Nx是控制体数量若边界条件是 Dirichlet实际存储数组长度仍是Nx因为边界值被强制写到首末控制体上。使用时注意网格坐标在中心点采样而不是边缘。3.2 时间循环和残差监视主程序推开核心循环后结构通常是预分配结果数组再按预设频率保存剖面。下面的框架可以直接嵌进自己的脚本% 从 FVTdemo.m 中抽出的核心时间循环框架 phi zeros(Nx, Ny); % 初始浓度或温度场 phi_store zeros(Nx, Ny, nStep/nOut 1); phi_store(:,:,1) phi; for k 1:nStep % 组装右端项时间项phi_old 源项 rhs aP0 * phi; % 上一时刻的 phi顺序为列向量 rhs apply_source(rhs, x, y); % 自定义源项可以写成匿名函数 % 求解稀疏线性系统M 来自上一节的系数矩阵 phi_vec M \ rhs(:); % 边界值覆盖第一类边界 phi_vec(idx_left) phiL; phi_vec(idx_right) phiR; % 重塑回二维网格 phi reshape(phi_vec, Nx, Ny); % 每 nOut 步保存一次方便后续画动画 if mod(k, nOut) 0 phi_store(:,:, k/nOut1) phi; end end这里有个细节M \ rhs用的是 MATLAB 内置稀疏直接解法。非稳态问题一般只需一次 LU 分解之后每个时间步只做一次回代因此不必在循环里反复分解。如果问题变成非线性比如源项依赖phi^2就得把M拆成可更新的部分在循环内部重算对角线并用迭代求解器bicgstab代替直接解。3.3 可视化与数据导出示例运行结束后FVTdemo 会生成 diff_pde.jpg 和 diff_pde_3d.jpg这正是 MATLAB 画图指令contourf和surf的输出。下面的代码把某个时刻的剖面导出为 CSV方便在 Python 或 Excel 里复核% 导出中心线剖面数据 xmid x; % 取网格中心坐标 phi_mid squeeze(phi_store(:, Ny/2, end)); % 最后时刻中心线 T_out table(xmid(:), phi_mid(:), VariableNames, {x, phi}); writetable(T_out, fvm_profiles.csv);CSV 文件名不要带中文MATLAB 的writetable在部分 Linux 版本下对中文路径处理不够稳定。squeeze用于去掉二维数组里的单一维度避免导出后列数不整齐。这个 CSV 可以直接用readtable(fvm_profiles.csv)读回也可以扔给任何第三方库做差分对比。4. 不止一维球坐标、焓法和顶盖驱动空腔的扩展用法4.1 球坐标扩散和 FVTool、FiPy 的解析解对比资源里的diffusion1Dspherical_analytic_vs_FVTool_vs_Fipy.ipynb是一个 Jupyter Notebook核心任务是把一维球形扩散的数值解和解析解放在一起比较。球坐标下的径向扩散方程为 ∂T/∂t α/r² ∂/∂r(r² ∂T/∂r)直接套用笛卡尔坐标系数会出错。需要把控制体积改为球壳体积面面积是 4πr² 而不是常量。% 球坐标系下的控制体积几何修正 dr 0.005; r dr/2:dr:1; N length(r); rmid r(:); % 左界面和右界面半径 rL rmid - dr/2; rR rmid dr/2; % 体积: 球壳 (4/3)π(rR^3 - rL^3) V (4*pi/3) .* (rR.^3 - rL.^3); % 界面面积 AL 4*pi .* rL.^2; AR 4*pi .* rR.^2; % 扩散传导率取左右界面的几何平均值 D alpha .* (AL AR) ./ 2;注意在这个离散式里如果保留源项请先把它乘以控制体体积在球坐标中边界处的面积为零首末控制体的几何修正对结果影响很大。我在写这个示例时常常忽略rL(1)0处的面积导致边界通量丢失发现后改为显式设置对称边界条件数值解才与解析解重合。解析解来自球坐标系分离变量解通常表现为无限级数。Jupyter Notebook 里用 FiPy 做交叉验证的好处是两边用的是完全不同的离散实现若数值解和解析解、FiPy 解的偏差都在 1e-3 以内可以确认自己的 MATLAB FVM 主程序没有结构性错误。比较时建议固定同一个 CFL别单纯比计算时间。4.2 焓法处理相变IAPWS_IF97 的调用位置相变问题不能用单一温度方程直接求解因为潜热会让焓曲线出现拐点。phaseChangeEnthalpyMethodExample.m用的是焓法把能量守恒写成含焓的形式每步先更新焓再由 h-T 映射表读出温度。IAPWS_IF97.m 在这里不是求解器而是物性查询函数供应传热计算所需的密度、比焓和比热。我自己常用的调用模式是先查压力对应的饱和温度再返回焓% 访问 IAPWS_IF97 的几个典型入口 p 1e5; % 1 bar Tsat IAPWS_IF97(Tsat_p, p); % 饱和温度 h_liq IAPWS_IF97(h_pT, p, Tsat-1); % 过冷液焓 h_vap IAPWS_IF97(h_pT, p, Tsat1); % 过热蒸汽焓 latent h_vap - h_liq; % 汽化潜热IAPWS_IF97的接口有多个入口Tsat_p输入压力输出饱和温度h_pT输入压力和温度输出焓值。用的时候要留意单位IF97 公式的压力基准为 Pa温度基准为 K焓基准为 J/kg有的版本会将 MPa 混入导致焓值差几个数量级。焓法循环里每步都调用物性函数会明显变慢我一般先建一个焓-温度查找表循环里只查表并线性插值。示例文件求解对象附加依赖典型收敛判据diffusion1Dspherical_analytic_vs_...球坐标一维扩散FiPy / FVToolmaxphaseChangeEnthalpyMethodExample.m一维相变 Stefan 问题IAPWS_IF97.m液固界面位置误差 2%SteadyLidDrivenCavityExample.m不可压稳态顶盖驱动空腔无中心速度残差 1e-6FVTdemo.m瞬态对流扩散无总浓度相对变化 1e-54.3 顶盖驱动空腔从标量输运到矢量流场SteadyLidDrivenCavityExample.m呈现的是不可压 Navier-Stokes 的基准算例。它和瞬态对流扩散共享 FVM 思想但多了压力-速度耦合常见做法是 SIMPLE 算法先假定压力场解动量方程得到速度场再修正压力让速度满足连续性方程。顶盖驱动空腔的边界条件很直观左、右、下边界为无滑移上边界给水平速度 U1。我在实际跑这个例子时只需要把上一步的标量扩散求解器嵌套进两个循环外层迭代压力修正内层解两个方向的动量方程。边界条件不再是简单的 Dirichlet而是速度分量的逐一指定。如果发现压力振荡通常是没有在网格上使用交错布置或压力参考点选得不对。这个示例的意义是告诉你FVM 的核心离散一旦写对从标量方程到矢量方程只是加装耦合器的过程。5. 收敛性检查与常见坑CFL、残差和 CSV 导出5.1 显示 CFL 检查和网格独立性瞬态对流扩散最容易出的问题是时间步长太大前沿出现“越界”振荡。即使在隐式格式里太大 dt 也会造成物理上的失真。我通常在每个脚本开头计算一次 CFL 数% 显式上限参考隐式可放宽但不建议盲目放大 cfl u * dt / dx; if cfl 1 warning(CFL%.2f 超过1建议减小 dt 或加密 dx, cfl); end这里 dx 取所有方向的最小网格间距u 取最大速度模。CFL 不满足时隐式格式虽然不会发散但数值解会把对流前沿抹得更平看起来像人为增大了扩散系数。网格独立性检查比收敛判据更实用把网格数从 40 翻倍到 80如果同一时刻的剖面最大变化小于 2%可以进入参数研究如果还在大幅移动则继续加密并同步调整 dt。5.2 质量守恒残差监视有限体积法最好的自检指标是总质量守恒。在每个时间步记录 sum(phi.*V) 随时间的变化理想情况下应保持一个常数或等于边界注入总量。下面的代码片段可以加到循环末尾% 统计总质量/总能量输出到命令行 total sum(phi(:) .* V(:)); % 一维时 V 是控制体积 fprintf(%4d total %.8f\n, k, total);若总质量在某一步突然跳变先去查边界条件如果边界是 Neumann应确认通量项没有在矩阵行里被意外清零。这个检查比盯着等值线图直观得多。5.3 把结果发给其他工具分析画完 3D 图后需要把几个时刻的数据连同坐标一起导出。用writetable导出 CSV 已经很常见如果你的同事用 Python 做后处理建议保存 numpy 能读的纯数字矩阵不要混入文本表头。% 保存二进制数据避免 CSV 在大网格下读写太慢 save(phi_final.mat, phi, x, dt, nStep);或者用dlmwrite(phi_final.csv, phi, precision, %.6e)只写数字矩阵。导出的 CSV 可直接用readmatrix读回 MATLAB这样后续做 FFT 仿真或频谱分析时不用重新跑一遍数值解。读入 Python 时用np.loadtxt完全避开 MATLAB 表格表头的解析问题。这也是我在多语言对照验证时最常用的导出方式。本文还有配套的精品资源点击获取