简介这是面向流体力学初学者与MATLAB用户的二维非定常Navier-Stokes方程求解示例聚焦不可压缩流体在二维空间内的运动模拟适合学习有限差分、时间推进及压力修正等核心数值方法。资源共76个文件含44个m脚本、29个png结果图、1个txt说明及1个html演示页覆盖网格定义、边界条件、初始化、数值求解与可视化全流程压缩包约463KB。文件组织清晰主程序与各功能模块分离便于按步骤理解代码逻辑。目前已有1941人浏览学习。通过该示例可掌握使用quiver绘制速度矢量场、contourf显示流线等可视化技巧同时加深对N-S方程中压力梯度项、粘性扩散项及连续性约束的实际认识为后续研究飞行器空气动力学、水力学等工程问题打下基础。1. 二维非定常 N-S 方程这个 MATLAB 示例能让你把“流动仿真”从玄学变成看得见的曲线做过流体仿真的同学都有这种体验学了一堆 N-S 方程推导真要自己在 MATLAB 里跑一个非定常流场却不知道第一行代码该写什么。网上的开源求解器动不动就是 C 几万行装依赖就得折腾半天教科书附录里的示例代码又往往只给定常版一碰到时间项就集体沉默。这份“二维 Navier-Stokes非定常MATLAB 示例”把整个问题压缩到了一个可读、可改、能跑出流场图的工程里它用流函数涡量法离散不可压缩流动的控制方程通过时间步进模拟顶盖驱动方腔流这一类经典算例输出速度场、涡量场和流函数分布。你不用理解全部数值方法细节照着改雷诺数、改网格密度就能看到涡旋从生成到发展的全过程。适合正在选做课程设计、刚接手流体仿真任务、或者想找一个干净底子做二次开发的 MATLAB 使用者。我把它完整拆了一遍从公式到代码行逐一讲清楚。2. 数学模型与离散方案为什么选流函数涡量法而不是直接解压力2.1 非定常不可压缩 N-S 方程的完整形式在二维直角坐标系下不可压缩非定常 Navier-Stokes 方程组的完整形式包含三个分量方程两个动量方程和一个连续性方程。速度分量 u 和 v 分别沿 x 和 y 方向p 为压力ν 为运动黏度% 动量方程x 方向 % du/dt u*du/dx v*du/dy -1/rho * dp/dx nu * (d2u/dx2 d2u/dy2) % 动量方程y 方向 % dv/dt u*dv/dx v*dv/dy -1/rho * dp/dy nu * (d2v/dx2 d2v/dy2) % 连续性方程 % du/dx dv/dy 0直接数值求解这套原始变量u,v,p的非定常方程会遇到一个经典麻烦连续性方程没有显式的时间演化形式压力场的求解需要反复迭代修正程序写出来又复杂又容易发散。这个示例采用的处理方式是引入流函数 ψ 和涡量 ω把问题从“三个方程耦合求解”降维成“两个方程分步求解”稳定性也更好控制。这是上世纪七八十年代发展起来、至今仍广泛用于教学和工程评估的经典方案。2.2 非定常涡量输运方程与运动学关系的离散通过取旋度消去压力梯度项得到非定常涡量输运方程% 涡量输运方程非定常核心 % dω/dt u*dω/dx v*dω/dy nu * (d2ω/dx2 d2ω/dy2) % 其中 ω dv/dx - du/dy % 而流函数满足泊松方程d2ψ/dx2 d2ψ/dy2 -ω % 速度可由流函数恢复u dψ/dyv -dψ/dx方程里每一项的物理含义值得说清楚左边第一项是涡量随时间的变化率即非定常项这是本示例区别于定常求解器的核心第二、三项是涡量的对流输运它让流场中的涡旋结构随时间移动和变形右边是分子黏性造成的涡量扩散。整个时间推进的循环就是先根据当前时刻的速度场计算涡量的对流和扩散增量更新得到下一时刻的涡量场再由更新后的涡量场求解泊松方程得到流函数场最后从流函数差分出速度场进入下一步。每一步之间环环相扣是理解这段代码的主线。空间离散采用二阶中心差分时间推进采用显式 Euler——虽然精度只到一阶但胜在实现简单、代码短符合示例定位。在实际做项目时我会把时间推进换成 Adams-Bashforth 二阶格式这个在后面“改造”小节给出具体做法。对于显式格式需要注意步长限制时间步增量不能超过网格尺寸和对流速度共同决定的 CFL 上限也不能超过黏性扩散的稳定极限否则算到一半流场就炸了。3. 程序模块拆解从参数初始化到流场可视化的完整链路3.1 工程目录与文件组织这个示例代码设计得比较轻量整体目录结构大致如下NS2D_Unsteady/ ├── main.m % 主程序参数设置、时间循环、后处理 ├── initField.m % 初始化流函数、涡量、速度场 ├── solvePoisson.m % 泊松方程求解流函数 ├── vorticityUpdate.m % 涡量输运方程的时间推进 ├── boundaryCondition.m % 边界涡量与流函数的约束 └── plotResult.m % 可视化输出主程序负责装配整个计算流程各子函数按职责拆分。这种结构最大的优点是把“数值方法”和“数据流程”分开了改边界条件只需要动 boundaryCondition.m改时间推进格式只需要动 vorticityUpdate.m不需要在一大段顺序代码里寻找修改点。如果只是课程作业把所有子函数并进一个 main.m 也能跑但二次开发时你会后悔没有遵守这个模块边界。3.2 参数初始化与网格生成主程序开头是参数设置段。以顶盖驱动方腔流为例方腔边长为 L网格等分为 nx×ny 个单元雷诺数 Re 是核心无量纲参数% 顶盖驱动方腔流参数设置 L 1.0; % 方腔边长 nx 64; ny 64; % 网格点数 Re 100; % 雷诺数定义Re U_top * L / nu nu U_top * L / Re; % 运动黏度 U_top 1.0; % 顶盖运动速度 dt 0.001; % 时间步长 Tmax 20; % 总模拟时长 dx L / (nx-1); % 空间步长 dy L / (ny-1); % 空间步长 x linspace(0, L, nx); y linspace(0, L, ny);雷诺数直接控制流态Re100 时流场有一个稳定的主涡Re1000 时角落会出现次级涡Re5000 以上流动变得非定常甚至混乱。这个参数是研究流动状态最重要的旋钮改它比改网格更直观地影响流场结构。初学者很容易把这里的运动黏度计算漏掉直接写 nu0.01 之类随意数值这会让“雷诺数”名存实亡后续和文献结果对比就失去了依据。3.3 涡量输运方程的时间推进实现涡量更新的核心代码不长但集中了这个示例的数值方法精华。显式时间推进的意思很直接用当前时刻的涡量空间导数去估计时间变化率然后外推一个时间步长。这里给出 vorticityUpdate.m 的关键循环体% 涡量输运方程显式时间推进以二维数组方式实现循环速度远快于双重 for function omega_new vorticityUpdate(omega, psi, dx, dy, dt, nu, nx, ny) % 由流函数计算速度场u dpsi/dyv -dpsi/dx中心差分 u zeros(nx, ny); v zeros(nx, ny); u(2:end-1, 2:end-1) (psi(2:end-1, 3:end) - psi(2:end-1, 1:end-2)) / (2*dy); v(2:end-1, 2:end-1) -(psi(3:end, 2:end-1) - psi(1:end-2, 2:end-1)) / (2*dx); % 涡量的空间导数中心差分 dOmega_dx (omega(2:end-1, 3:end) - omega(2:end-1, 1:end-2)) / (2*dx); dOmega_dy (omega(3:end, 2:end-1) - omega(1:end-2, 2:end-1)) / (2*dy); % 拉普拉斯算子二阶中心差分 lapOmega (omega(2:end-1, 3:end) - 2*omega(2:end-1, 2:end-1) omega(2:end-1, 1:end-2)) / dx^2 ... (omega(3:end, 2:end-1) - 2*omega(2:end-1, 2:end-1) omega(1:end-2, 2:end-1)) / dy^2; % 时间推进dω/dt -u*dω/dx - v*dω/dy nu * lapOmega omega_new zeros(nx, ny); omega_new(2:end-1, 2:end-1) omega(2:end-1, 2:end-1) dt * ( ... -u(2:end-1, 2:end-1) .* dOmega_dx ... -v(2:end-1, 2:end-1) .* dOmega_dy ... nu * lapOmega); end参数含义对应地解释一下u 和 v 是当前时刻速度场由流函数 ψ 差分得到dOmega_dx、dOmega_dy 是对流项的一阶导数lapOmega 是黏性扩散项的二阶导数。时间步进公式中对流项是“涡量被流场携带运动”的数学表达扩散项是“涡量从高浓度向低浓度抹平”的数学表达。代码里使用 MATLAB 矩阵整体运算而不是 for 循环嵌套是一个很实用的性能技巧同样的网格规模向量化写法一般比嵌套循环快一个数量级时间循环总共要跑上万步这个差距直接决定你是喝一杯咖啡还是等一晚。3.4 泊松方程求解与边界条件每个时间步内拿到新的涡量场后需要求解流函数泊松方程。示例直接调用 MATLAB 的反斜杠算子求解稀疏线性方程组% 泊松方程d2(psi)/dx2 d2(psi)/dy2 -omega % 组装五点差分格式的系数矩阵 A N (nx-2) * (ny-2); A sparse(N, N); rhs zeros(N, 1); for j 2:ny-1 for i 2:nx-1 idx (j-2)*(nx-2) (i-1); A(idx, idx) -2/dx^2 - 2/dy^2; if i 2, A(idx, idx-1) 1/dx^2; end if i nx-1, A(idx, idx1) 1/dx^2; end if j 2, A(idx, idx-(nx-2)) 1/dy^2; end if j ny-1, A(idx, idx(nx-2)) 1/dy^2; end rhs(idx) -omega(i, j); end end psi_vec A \ rhs; % 稀疏矩阵求解 psi zeros(nx, ny); psi(2:end-1, 2:end-1) reshape(psi_vec, nx-2, ny-2);A 矩阵每一行对应方腔内部一个网格点的离散方程对角线是 -2/dx²-2/dy²四邻域各有权重。方程右侧 -ω 就是涡量对流的“源项”物理意义是涡旋结构决定了流动的形态。用反斜杠求解的好处是 MATLAB 自动选用稀疏直接法或迭代法对 64×64 这类网格在单时间步内达到毫秒级延迟。需要知道的一个局限是矩阵组装用了两层 for 循环只在初始化时执行一次所以即使写得粗糙也不会成为性能瓶颈真正决定程序快慢的是时间循环里的涡量更新段。边界条件是这个示例中最容易写错、也最能体现数值功底的地方。顶盖驱动流对 ψ 和 ω 的约束如下% 边界条件设置 % 顶盖y Lu U_topv 0dpsi/dy U_top % 底面、左壁、右壁无滑移psi 0即不穿透 % 边界涡量由流函数导数反推 % 左壁omega(1,j) -2 * psi(2,j) / dx^2 % 右壁omega(nx,j) -2 * psi(nx-1,j) / dx^2 % 底面omega(i,1) -2 * psi(i,2) / dy^2 % 顶盖omega(i,ny) -2 * psi(i,ny-1) / dy^2 - 2*U_top/dy边界涡量的公式来自“无滑移边界下单侧导数近似”壁面涡量正比于紧邻壁面的流函数值。顶盖那一项里多出来的 -2U_top/dy是因为顶盖在运动剪切层内部顶着有限速度差。很多人在这一步忽略顶盖项结果就是流场整体形状对但顶盖附近的涡量峰值明显偏小与文献基准解一比就露馅了。3.5 可视化输出示例自带绘图段输出三个子图涡量云图、流函数等值线、速度矢量图。绘图的代码比较常规% 可视化当前时刻流场 figure(1) subplot(1, 3, 1) contourf(x, y, omega, 20); colorbar; title(涡量场) subplot(1, 3, 2) contour(x, y, psi, 20); title(流函数) subplot(1, 3, 3) quiver(x(1:4:end), y(1:4:end), u(1:4:end, 1:4:end), v(1:4:end, 1:4:end)) title(速度场) drawnow可视化不仅可以用来“看个大概”还有实际价值观察涡量云图是否对称、流函数等值线是否光滑能快速判断边界条件有没有写对。绘图本身只消耗很少的时间比例建议每几百步输出一次而不要每一步都画图——否则 matplotlib 式的地帧重绘会让时间循环显著变慢。实际做项目时我会把特定时刻的完整流场数据存成 .mat 文件后处理单独做避免反复重算。4. 运行指南与算例验证把命令跑通再把结果跑对4.1 从零到第一个流场图的完整步骤拿到源码后的标准操作流程如下在 MATLAB 里打开 main.m直接 F5 运行。无需安装任何工具箱稀疏矩阵反斜杠、contourf、quiver 这些都是 MATLAB 基础功能。文件里如果包含子函数把它们全部加入当前工作目录或统一放在一个文件夹中。运行过程中命令行窗口会输出当前时间步和模拟时间几秒后弹出流场图。% 在 MATLAB 命令行中运行 main % 或者直接点击编辑器中的“运行”按钮建议第一次跑用默认参数64×64 网格、Re100、总模拟时长 Tmax20。Re100 属于层流稳态工况流场在 5~10 个无量纲时间单位后就基本稳定你会看到一个位于方腔中央偏右的大涡。这时把 Re 提高到 1000缩短 Tmax 到 10 左右观察涡心偏移和角涡出现的规律。如果出现 NaN 或者数值发散优先减小 dt。4.2 如何验证结果与经典基准解对比跑通不等于跑对。几乎所有流体数值代码都需要和经典基准解对照否则改了一个参数后结果逐渐失真你却浑然不知。有两类最常用的对比方式。第一类是流场特征对照方腔中心剖面上 u 速度沿竖直方向的分布曲线与文献中的基准数据对比Re100 时曲线应当重合第二类是涡量极值对照顶盖处的涡量最大值随雷诺数的变化规律。后者尤其适合检查边界涡量公式是否实现正确。% 提取竖直中线上的 u 速度分布用于结果对比 midLine round(nx/2); uVertical u(midLine, :); plot(uVertical, y); % 应近似关于 y 对称的 S 形曲线 hold on; % 此处可加载文献基准数据 plot(ref_u, ref_y, ro)需要特别说明的是不要期望曲线与文献基准解完全重合。网格分辨率 64×64 加上一阶时间精度中心线上尾部的误差在 1% 到 3% 之间属于正常范围误差偏大且曲线形状不一才需要考虑数值问题。4.3 参数扫描方法与实验设计研究流动的一个常见需求是“看不同雷诺数下的流场演化”。写一个外层循环自动扫描 Re 值并把各个算例的结果保存成独立文件即可。另一个实用的实验设计是“固定 Re、加密网格”观察结果是否收敛——如果 Re100 下 32×32 和 128×128 网格得到的速度曲线差异很大说明你的网格还没达到网格无关性后续所有结果都要打一个问号。这个“网格无关性验证”是论文和报告中必须呈现的核心内容示例代码同样支持这种用法直接改 nx、ny 重新运行即可。5. 常见问题排查关于发散、边界涡量和可视化异常的排雷日志5.1 时间步进几十步后流场出现 NaN现象程序运行 20 到 50 个时间步后涡量场突然出现 NaN随后整个流场崩溃绘图窗口空白或显示错误。原因分析显式格式的时间步长超出了稳定极限。涡量输运方程包含对流项和扩散项对流项的稳定性要求 dt ≤ dx/|u|max扩散项要求 dt ≤ dx²/(4ν)。实际程序中报错前的几万步通常都是“潜伏期”网格越细黏度越低潜伏期越短。解决办法也很直接把 dt 从 0.001 调整到 0.0005 或更低或者更系统地写一个自动计算最大稳定步长的函数。有经验之后我习惯把 dt 取三个限制量中的最小值再乘 0.5 的保险系数彻底告别这类问题。5.2 边界附近涡量云图出现锯齿状等值线现象涡量云图整体轮廓正常但紧贴壁面处等值线不再光滑出现锯齿纹理。原因分析边界涡量公式中的单侧差分只有一阶精度内部点使用二阶中心差分这种精度不匹配会在边界引起虚假振荡用高分辨率云图看得特别清楚。解决办法用二阶精度的单侧公式替代一阶公式。例如左壁的涡量计算公式从 -2ψ(2,j)/dx² 改写为 (-3ψ(1,j) 4ψ(2,j) - ψ(3,j)) / (2dx²)。改完边界涡量等值线明显光滑。这里需要提醒一个反直觉的事实边界涡量的精度直接决定整体求解精度不是边界内部各占一半而是边界负责给内部提供“锚点”锚歪了内部再准也无济于事。5.3 速度矢量图看起来方向和预期相反现象顶盖驱动流中靠近顶盖的水平速度方向与顶盖运动方向相反流场存在一个明显的大涡但旋转方向反了。原因分析流函数与速度的关系符号约定不一致常见的是 ψ 的 y 向导数符号弄反。解决办法检查速度恢复公式正确关系是 u ∂ψ/∂y、v -∂ψ/∂x。这里有一个查验方法方腔内一个大涡对应 ψ 呈山峰状分布沿 x 方向看 ψ 上升、沿 y 方向看 ψ 下降就能推断出正确的速度分量方向。这类错误与雷洛数、网格无关是纯粹的符号约定问题思想方法最有效。5.4 高位 Re 下程序越跑越慢、内存占用飙升现象把 Re 从 100 改到 10000 后程序快要完成时突然变得极慢甚至报告内存不足。原因分析高频流动出现小尺度涡结构在固定网格下 CFL 要求时间步大幅缩小总步数翻倍。而 MATLAB 每次循环如果意外产生了完整矩阵的拼接或动态增长会出现内存碎片化。解决办法提前用 zeros 初始化所有大型数组如果确实需要高频非定常计算则要考虑加密网格、减少总模拟时长、或将结果降采样保存。需要说明的是显式方法的扩散项对时间步长的限制与 dt ≤ dx²/(2ν) 相关即细网格、低黏度下时间步会急剧压缩这是显式方法的固有局限改用隐式处理扩散项可以绕过这个坎但实现复杂度也上来了。6. 二次开发实战把“能跑的示例”改造成“能研究问题的小型求解器”先说第一个建议动手改之前先跑一遍原版代码并对照标准算例验证确认你手上的版本是“基线正确”的。任何后续改造都要基于这个基线否则会遭遇“改了十处代码和跑出一个不知道哪里的错”的混乱局面。基线验证通过后按由浅到深的顺序选择改造点第一步把时间推进从显式 Euler 升级为二阶 Adams-Bashforth。做法是保存前一步的涡量增量当前增量与上一步增量做加权组合ω(n1) ω(n) dt×(1.5×F(n) - 0.5×F(n-1))其中 F 表示涡量输运方程的右端。这个改动只涉及循环内的几行代码但时间精度从一阶升到二阶长期演化更稳而且为后续加源项提供了干净入口。这是“性价比”最高的一处改造。第二步改造几何模型。原示例默认方腔流修改壁面形状需要更换网格生成和边界条件两处。比如改成圆柱绕流抛弃笛卡尔网格是更好的选择整体框架才换得更彻底但改成长宽比不同的矩形通道流则只需调整 L 和 nx、ny几乎零成本。这类保守改造恰好适合练手。第三步把压力场恢复出来。流函数-涡量法算完不直接给压力场研究受力就需要恢复压力。做法是对动量方程取散度得到压力泊松方程边界用诺伊曼条件这与原示例的泊松求解套路一致可以直接复用稀疏矩阵框架。恢复压力后可以计算圆柱或翼型的升阻力系数后续再研究降噪、受力分析等工程问题才有根基。第四步把算例参数文件独立出来。将 Re、网格尺寸、求解时长、输出频率全部放进一个结构体配合循环脚本自动批量跑参数工况批量后处理时效率极高。我现在做多工况对比都是这套流程先把单例结果验证到和文献数据对得上再放开参数扫描最后统一出图和分析。从那以后我每次改造求解器都强制走一遍“基线复现到确认没问题再动刀”的流程改一个参数就重新观察一次流场演化杜绝原地大改。希望你也能把求解器的每一次改动变成一次可以回滚、可以对照的实验记录这比灵光一现的“改完就跑”高效得多。希望帮到你。本文还有配套的精品资源点击获取
简介这是面向流体力学初学者与MATLAB用户的二维非定常Navier-Stokes方程求解示例聚焦不可压缩流体在二维空间内的运动模拟适合学习有限差分、时间推进及压力修正等核心数值方法。资源共76个文件含44个m脚本、29个png结果图、1个txt说明及1个html演示页覆盖网格定义、边界条件、初始化、数值求解与可视化全流程压缩包约463KB。文件组织清晰主程序与各功能模块分离便于按步骤理解代码逻辑。目前已有1941人浏览学习。通过该示例可掌握使用quiver绘制速度矢量场、contourf显示流线等可视化技巧同时加深对N-S方程中压力梯度项、粘性扩散项及连续性约束的实际认识为后续研究飞行器空气动力学、水力学等工程问题打下基础。1. 二维非定常 N-S 方程这个 MATLAB 示例能让你把“流动仿真”从玄学变成看得见的曲线做过流体仿真的同学都有这种体验学了一堆 N-S 方程推导真要自己在 MATLAB 里跑一个非定常流场却不知道第一行代码该写什么。网上的开源求解器动不动就是 C 几万行装依赖就得折腾半天教科书附录里的示例代码又往往只给定常版一碰到时间项就集体沉默。这份“二维 Navier-Stokes非定常MATLAB 示例”把整个问题压缩到了一个可读、可改、能跑出流场图的工程里它用流函数涡量法离散不可压缩流动的控制方程通过时间步进模拟顶盖驱动方腔流这一类经典算例输出速度场、涡量场和流函数分布。你不用理解全部数值方法细节照着改雷诺数、改网格密度就能看到涡旋从生成到发展的全过程。适合正在选做课程设计、刚接手流体仿真任务、或者想找一个干净底子做二次开发的 MATLAB 使用者。我把它完整拆了一遍从公式到代码行逐一讲清楚。2. 数学模型与离散方案为什么选流函数涡量法而不是直接解压力2.1 非定常不可压缩 N-S 方程的完整形式在二维直角坐标系下不可压缩非定常 Navier-Stokes 方程组的完整形式包含三个分量方程两个动量方程和一个连续性方程。速度分量 u 和 v 分别沿 x 和 y 方向p 为压力ν 为运动黏度% 动量方程x 方向 % du/dt u*du/dx v*du/dy -1/rho * dp/dx nu * (d2u/dx2 d2u/dy2) % 动量方程y 方向 % dv/dt u*dv/dx v*dv/dy -1/rho * dp/dy nu * (d2v/dx2 d2v/dy2) % 连续性方程 % du/dx dv/dy 0直接数值求解这套原始变量u,v,p的非定常方程会遇到一个经典麻烦连续性方程没有显式的时间演化形式压力场的求解需要反复迭代修正程序写出来又复杂又容易发散。这个示例采用的处理方式是引入流函数 ψ 和涡量 ω把问题从“三个方程耦合求解”降维成“两个方程分步求解”稳定性也更好控制。这是上世纪七八十年代发展起来、至今仍广泛用于教学和工程评估的经典方案。2.2 非定常涡量输运方程与运动学关系的离散通过取旋度消去压力梯度项得到非定常涡量输运方程% 涡量输运方程非定常核心 % dω/dt u*dω/dx v*dω/dy nu * (d2ω/dx2 d2ω/dy2) % 其中 ω dv/dx - du/dy % 而流函数满足泊松方程d2ψ/dx2 d2ψ/dy2 -ω % 速度可由流函数恢复u dψ/dyv -dψ/dx方程里每一项的物理含义值得说清楚左边第一项是涡量随时间的变化率即非定常项这是本示例区别于定常求解器的核心第二、三项是涡量的对流输运它让流场中的涡旋结构随时间移动和变形右边是分子黏性造成的涡量扩散。整个时间推进的循环就是先根据当前时刻的速度场计算涡量的对流和扩散增量更新得到下一时刻的涡量场再由更新后的涡量场求解泊松方程得到流函数场最后从流函数差分出速度场进入下一步。每一步之间环环相扣是理解这段代码的主线。空间离散采用二阶中心差分时间推进采用显式 Euler——虽然精度只到一阶但胜在实现简单、代码短符合示例定位。在实际做项目时我会把时间推进换成 Adams-Bashforth 二阶格式这个在后面“改造”小节给出具体做法。对于显式格式需要注意步长限制时间步增量不能超过网格尺寸和对流速度共同决定的 CFL 上限也不能超过黏性扩散的稳定极限否则算到一半流场就炸了。3. 程序模块拆解从参数初始化到流场可视化的完整链路3.1 工程目录与文件组织这个示例代码设计得比较轻量整体目录结构大致如下NS2D_Unsteady/ ├── main.m % 主程序参数设置、时间循环、后处理 ├── initField.m % 初始化流函数、涡量、速度场 ├── solvePoisson.m % 泊松方程求解流函数 ├── vorticityUpdate.m % 涡量输运方程的时间推进 ├── boundaryCondition.m % 边界涡量与流函数的约束 └── plotResult.m % 可视化输出主程序负责装配整个计算流程各子函数按职责拆分。这种结构最大的优点是把“数值方法”和“数据流程”分开了改边界条件只需要动 boundaryCondition.m改时间推进格式只需要动 vorticityUpdate.m不需要在一大段顺序代码里寻找修改点。如果只是课程作业把所有子函数并进一个 main.m 也能跑但二次开发时你会后悔没有遵守这个模块边界。3.2 参数初始化与网格生成主程序开头是参数设置段。以顶盖驱动方腔流为例方腔边长为 L网格等分为 nx×ny 个单元雷诺数 Re 是核心无量纲参数% 顶盖驱动方腔流参数设置 L 1.0; % 方腔边长 nx 64; ny 64; % 网格点数 Re 100; % 雷诺数定义Re U_top * L / nu nu U_top * L / Re; % 运动黏度 U_top 1.0; % 顶盖运动速度 dt 0.001; % 时间步长 Tmax 20; % 总模拟时长 dx L / (nx-1); % 空间步长 dy L / (ny-1); % 空间步长 x linspace(0, L, nx); y linspace(0, L, ny);雷诺数直接控制流态Re100 时流场有一个稳定的主涡Re1000 时角落会出现次级涡Re5000 以上流动变得非定常甚至混乱。这个参数是研究流动状态最重要的旋钮改它比改网格更直观地影响流场结构。初学者很容易把这里的运动黏度计算漏掉直接写 nu0.01 之类随意数值这会让“雷诺数”名存实亡后续和文献结果对比就失去了依据。3.3 涡量输运方程的时间推进实现涡量更新的核心代码不长但集中了这个示例的数值方法精华。显式时间推进的意思很直接用当前时刻的涡量空间导数去估计时间变化率然后外推一个时间步长。这里给出 vorticityUpdate.m 的关键循环体% 涡量输运方程显式时间推进以二维数组方式实现循环速度远快于双重 for function omega_new vorticityUpdate(omega, psi, dx, dy, dt, nu, nx, ny) % 由流函数计算速度场u dpsi/dyv -dpsi/dx中心差分 u zeros(nx, ny); v zeros(nx, ny); u(2:end-1, 2:end-1) (psi(2:end-1, 3:end) - psi(2:end-1, 1:end-2)) / (2*dy); v(2:end-1, 2:end-1) -(psi(3:end, 2:end-1) - psi(1:end-2, 2:end-1)) / (2*dx); % 涡量的空间导数中心差分 dOmega_dx (omega(2:end-1, 3:end) - omega(2:end-1, 1:end-2)) / (2*dx); dOmega_dy (omega(3:end, 2:end-1) - omega(1:end-2, 2:end-1)) / (2*dy); % 拉普拉斯算子二阶中心差分 lapOmega (omega(2:end-1, 3:end) - 2*omega(2:end-1, 2:end-1) omega(2:end-1, 1:end-2)) / dx^2 ... (omega(3:end, 2:end-1) - 2*omega(2:end-1, 2:end-1) omega(1:end-2, 2:end-1)) / dy^2; % 时间推进dω/dt -u*dω/dx - v*dω/dy nu * lapOmega omega_new zeros(nx, ny); omega_new(2:end-1, 2:end-1) omega(2:end-1, 2:end-1) dt * ( ... -u(2:end-1, 2:end-1) .* dOmega_dx ... -v(2:end-1, 2:end-1) .* dOmega_dy ... nu * lapOmega); end参数含义对应地解释一下u 和 v 是当前时刻速度场由流函数 ψ 差分得到dOmega_dx、dOmega_dy 是对流项的一阶导数lapOmega 是黏性扩散项的二阶导数。时间步进公式中对流项是“涡量被流场携带运动”的数学表达扩散项是“涡量从高浓度向低浓度抹平”的数学表达。代码里使用 MATLAB 矩阵整体运算而不是 for 循环嵌套是一个很实用的性能技巧同样的网格规模向量化写法一般比嵌套循环快一个数量级时间循环总共要跑上万步这个差距直接决定你是喝一杯咖啡还是等一晚。3.4 泊松方程求解与边界条件每个时间步内拿到新的涡量场后需要求解流函数泊松方程。示例直接调用 MATLAB 的反斜杠算子求解稀疏线性方程组% 泊松方程d2(psi)/dx2 d2(psi)/dy2 -omega % 组装五点差分格式的系数矩阵 A N (nx-2) * (ny-2); A sparse(N, N); rhs zeros(N, 1); for j 2:ny-1 for i 2:nx-1 idx (j-2)*(nx-2) (i-1); A(idx, idx) -2/dx^2 - 2/dy^2; if i 2, A(idx, idx-1) 1/dx^2; end if i nx-1, A(idx, idx1) 1/dx^2; end if j 2, A(idx, idx-(nx-2)) 1/dy^2; end if j ny-1, A(idx, idx(nx-2)) 1/dy^2; end rhs(idx) -omega(i, j); end end psi_vec A \ rhs; % 稀疏矩阵求解 psi zeros(nx, ny); psi(2:end-1, 2:end-1) reshape(psi_vec, nx-2, ny-2);A 矩阵每一行对应方腔内部一个网格点的离散方程对角线是 -2/dx²-2/dy²四邻域各有权重。方程右侧 -ω 就是涡量对流的“源项”物理意义是涡旋结构决定了流动的形态。用反斜杠求解的好处是 MATLAB 自动选用稀疏直接法或迭代法对 64×64 这类网格在单时间步内达到毫秒级延迟。需要知道的一个局限是矩阵组装用了两层 for 循环只在初始化时执行一次所以即使写得粗糙也不会成为性能瓶颈真正决定程序快慢的是时间循环里的涡量更新段。边界条件是这个示例中最容易写错、也最能体现数值功底的地方。顶盖驱动流对 ψ 和 ω 的约束如下% 边界条件设置 % 顶盖y Lu U_topv 0dpsi/dy U_top % 底面、左壁、右壁无滑移psi 0即不穿透 % 边界涡量由流函数导数反推 % 左壁omega(1,j) -2 * psi(2,j) / dx^2 % 右壁omega(nx,j) -2 * psi(nx-1,j) / dx^2 % 底面omega(i,1) -2 * psi(i,2) / dy^2 % 顶盖omega(i,ny) -2 * psi(i,ny-1) / dy^2 - 2*U_top/dy边界涡量的公式来自“无滑移边界下单侧导数近似”壁面涡量正比于紧邻壁面的流函数值。顶盖那一项里多出来的 -2U_top/dy是因为顶盖在运动剪切层内部顶着有限速度差。很多人在这一步忽略顶盖项结果就是流场整体形状对但顶盖附近的涡量峰值明显偏小与文献基准解一比就露馅了。3.5 可视化输出示例自带绘图段输出三个子图涡量云图、流函数等值线、速度矢量图。绘图的代码比较常规% 可视化当前时刻流场 figure(1) subplot(1, 3, 1) contourf(x, y, omega, 20); colorbar; title(涡量场) subplot(1, 3, 2) contour(x, y, psi, 20); title(流函数) subplot(1, 3, 3) quiver(x(1:4:end), y(1:4:end), u(1:4:end, 1:4:end), v(1:4:end, 1:4:end)) title(速度场) drawnow可视化不仅可以用来“看个大概”还有实际价值观察涡量云图是否对称、流函数等值线是否光滑能快速判断边界条件有没有写对。绘图本身只消耗很少的时间比例建议每几百步输出一次而不要每一步都画图——否则 matplotlib 式的地帧重绘会让时间循环显著变慢。实际做项目时我会把特定时刻的完整流场数据存成 .mat 文件后处理单独做避免反复重算。4. 运行指南与算例验证把命令跑通再把结果跑对4.1 从零到第一个流场图的完整步骤拿到源码后的标准操作流程如下在 MATLAB 里打开 main.m直接 F5 运行。无需安装任何工具箱稀疏矩阵反斜杠、contourf、quiver 这些都是 MATLAB 基础功能。文件里如果包含子函数把它们全部加入当前工作目录或统一放在一个文件夹中。运行过程中命令行窗口会输出当前时间步和模拟时间几秒后弹出流场图。% 在 MATLAB 命令行中运行 main % 或者直接点击编辑器中的“运行”按钮建议第一次跑用默认参数64×64 网格、Re100、总模拟时长 Tmax20。Re100 属于层流稳态工况流场在 5~10 个无量纲时间单位后就基本稳定你会看到一个位于方腔中央偏右的大涡。这时把 Re 提高到 1000缩短 Tmax 到 10 左右观察涡心偏移和角涡出现的规律。如果出现 NaN 或者数值发散优先减小 dt。4.2 如何验证结果与经典基准解对比跑通不等于跑对。几乎所有流体数值代码都需要和经典基准解对照否则改了一个参数后结果逐渐失真你却浑然不知。有两类最常用的对比方式。第一类是流场特征对照方腔中心剖面上 u 速度沿竖直方向的分布曲线与文献中的基准数据对比Re100 时曲线应当重合第二类是涡量极值对照顶盖处的涡量最大值随雷诺数的变化规律。后者尤其适合检查边界涡量公式是否实现正确。% 提取竖直中线上的 u 速度分布用于结果对比 midLine round(nx/2); uVertical u(midLine, :); plot(uVertical, y); % 应近似关于 y 对称的 S 形曲线 hold on; % 此处可加载文献基准数据 plot(ref_u, ref_y, ro)需要特别说明的是不要期望曲线与文献基准解完全重合。网格分辨率 64×64 加上一阶时间精度中心线上尾部的误差在 1% 到 3% 之间属于正常范围误差偏大且曲线形状不一才需要考虑数值问题。4.3 参数扫描方法与实验设计研究流动的一个常见需求是“看不同雷诺数下的流场演化”。写一个外层循环自动扫描 Re 值并把各个算例的结果保存成独立文件即可。另一个实用的实验设计是“固定 Re、加密网格”观察结果是否收敛——如果 Re100 下 32×32 和 128×128 网格得到的速度曲线差异很大说明你的网格还没达到网格无关性后续所有结果都要打一个问号。这个“网格无关性验证”是论文和报告中必须呈现的核心内容示例代码同样支持这种用法直接改 nx、ny 重新运行即可。5. 常见问题排查关于发散、边界涡量和可视化异常的排雷日志5.1 时间步进几十步后流场出现 NaN现象程序运行 20 到 50 个时间步后涡量场突然出现 NaN随后整个流场崩溃绘图窗口空白或显示错误。原因分析显式格式的时间步长超出了稳定极限。涡量输运方程包含对流项和扩散项对流项的稳定性要求 dt ≤ dx/|u|max扩散项要求 dt ≤ dx²/(4ν)。实际程序中报错前的几万步通常都是“潜伏期”网格越细黏度越低潜伏期越短。解决办法也很直接把 dt 从 0.001 调整到 0.0005 或更低或者更系统地写一个自动计算最大稳定步长的函数。有经验之后我习惯把 dt 取三个限制量中的最小值再乘 0.5 的保险系数彻底告别这类问题。5.2 边界附近涡量云图出现锯齿状等值线现象涡量云图整体轮廓正常但紧贴壁面处等值线不再光滑出现锯齿纹理。原因分析边界涡量公式中的单侧差分只有一阶精度内部点使用二阶中心差分这种精度不匹配会在边界引起虚假振荡用高分辨率云图看得特别清楚。解决办法用二阶精度的单侧公式替代一阶公式。例如左壁的涡量计算公式从 -2ψ(2,j)/dx² 改写为 (-3ψ(1,j) 4ψ(2,j) - ψ(3,j)) / (2dx²)。改完边界涡量等值线明显光滑。这里需要提醒一个反直觉的事实边界涡量的精度直接决定整体求解精度不是边界内部各占一半而是边界负责给内部提供“锚点”锚歪了内部再准也无济于事。5.3 速度矢量图看起来方向和预期相反现象顶盖驱动流中靠近顶盖的水平速度方向与顶盖运动方向相反流场存在一个明显的大涡但旋转方向反了。原因分析流函数与速度的关系符号约定不一致常见的是 ψ 的 y 向导数符号弄反。解决办法检查速度恢复公式正确关系是 u ∂ψ/∂y、v -∂ψ/∂x。这里有一个查验方法方腔内一个大涡对应 ψ 呈山峰状分布沿 x 方向看 ψ 上升、沿 y 方向看 ψ 下降就能推断出正确的速度分量方向。这类错误与雷洛数、网格无关是纯粹的符号约定问题思想方法最有效。5.4 高位 Re 下程序越跑越慢、内存占用飙升现象把 Re 从 100 改到 10000 后程序快要完成时突然变得极慢甚至报告内存不足。原因分析高频流动出现小尺度涡结构在固定网格下 CFL 要求时间步大幅缩小总步数翻倍。而 MATLAB 每次循环如果意外产生了完整矩阵的拼接或动态增长会出现内存碎片化。解决办法提前用 zeros 初始化所有大型数组如果确实需要高频非定常计算则要考虑加密网格、减少总模拟时长、或将结果降采样保存。需要说明的是显式方法的扩散项对时间步长的限制与 dt ≤ dx²/(2ν) 相关即细网格、低黏度下时间步会急剧压缩这是显式方法的固有局限改用隐式处理扩散项可以绕过这个坎但实现复杂度也上来了。6. 二次开发实战把“能跑的示例”改造成“能研究问题的小型求解器”先说第一个建议动手改之前先跑一遍原版代码并对照标准算例验证确认你手上的版本是“基线正确”的。任何后续改造都要基于这个基线否则会遭遇“改了十处代码和跑出一个不知道哪里的错”的混乱局面。基线验证通过后按由浅到深的顺序选择改造点第一步把时间推进从显式 Euler 升级为二阶 Adams-Bashforth。做法是保存前一步的涡量增量当前增量与上一步增量做加权组合ω(n1) ω(n) dt×(1.5×F(n) - 0.5×F(n-1))其中 F 表示涡量输运方程的右端。这个改动只涉及循环内的几行代码但时间精度从一阶升到二阶长期演化更稳而且为后续加源项提供了干净入口。这是“性价比”最高的一处改造。第二步改造几何模型。原示例默认方腔流修改壁面形状需要更换网格生成和边界条件两处。比如改成圆柱绕流抛弃笛卡尔网格是更好的选择整体框架才换得更彻底但改成长宽比不同的矩形通道流则只需调整 L 和 nx、ny几乎零成本。这类保守改造恰好适合练手。第三步把压力场恢复出来。流函数-涡量法算完不直接给压力场研究受力就需要恢复压力。做法是对动量方程取散度得到压力泊松方程边界用诺伊曼条件这与原示例的泊松求解套路一致可以直接复用稀疏矩阵框架。恢复压力后可以计算圆柱或翼型的升阻力系数后续再研究降噪、受力分析等工程问题才有根基。第四步把算例参数文件独立出来。将 Re、网格尺寸、求解时长、输出频率全部放进一个结构体配合循环脚本自动批量跑参数工况批量后处理时效率极高。我现在做多工况对比都是这套流程先把单例结果验证到和文献数据对得上再放开参数扫描最后统一出图和分析。从那以后我每次改造求解器都强制走一遍“基线复现到确认没问题再动刀”的流程改一个参数就重新观察一次流场演化杜绝原地大改。希望你也能把求解器的每一次改动变成一次可以回滚、可以对照的实验记录这比灵光一现的“改完就跑”高效得多。希望帮到你。本文还有配套的精品资源点击获取
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
