简介这份MATLAB示例面向流体力学学习者与数值计算入门者围绕二维非定常Navier-Stokes方程展开帮助读者理解不可压缩流体在二维空间内的运动求解流程。资源包共76个文件以44个m脚本为核心涵盖网格定义、边界条件设置、速度场与压力初始化、时间推进、空间离散、压力修正及可视化等环节另有29个png结果图、1个html说明页和1个txt许可文件压缩包约463KB结构紧凑便于按模块查阅。已有1947人学习下载说明其在教学与自学场景中具有一定参考价值。读者可从中获得有限差分与时间迭代的实现思路、压力泊松方程修正的代码组织方式以及quiver、contourf等可视化手段的调用范例适合用于课程设计、算法验证与编程练习也能加深对流体力学基本概念的理解。1. 非定常二维 Navier-Stokes为什么“能跑出漩涡”不等于“算对了”很多人第一次在 MATLAB 里跑二维非定常 Navier-Stokes 示例看到流场里卷出几个漂亮的涡就觉得大功告成。但真正决定结果可信度的往往不是那几个涡而是压力泊松方程怎么解、对流项怎么离散、时间推进用显式还是半隐式。二维非定常 N-S 方程的核心难点在于速度场必须每一步都满足不可压约束而压力本身没有独立的演化方程它只是让散度归零的拉格朗日乘子。这份 MATLAB 示例的价值就在于把“投影法”这条主线完整摊开从动量方程预测速度到解压力泊松方程修正速度再到时间步进与边界处理。搞懂它你才能判断一个示例是教学玩具还是可复用的求解骨架。2. 投影法拆解压力为什么是“事后补”的2.1 不可压约束与压力泊松方程的来源二维非定常不可压 N-S 方程写成无量纲形式是两条动量方程描述速度随时间的演化连续方程要求速度散度为零。麻烦在于动量方程里含压力梯度但没有任何一条方程直接告诉压力怎么随时间变化。常见的实现思路是把动量方程先“半离散”用当前时刻的速度算一个中间速度这个中间速度只满足动量方程、不满足散度为零然后用压力梯度去修正它让修正后的速度散度归零。把修正关系代入连续方程就得到压力泊松方程。这一步是整个投影法的枢纽——压力不是物理上独立演化的量而是为了让速度场满足约束而反解出来的。理解这一点后面所有离散选择都有了依据。2.2 中间速度的预测步与修正步预测步的离散形式取决于对流项和粘性项怎么处理。显式处理对流、隐式处理粘性是最常见的折中对流项用当前速度直接算粘性项用隐式格式避免时间步受粘性限制。修正步则是用解出的压力梯度去减中间速度。这里有个容易忽略的细节修正步之后的速度在离散意义下散度是否严格为零取决于压力泊松方程的离散算子和散度算子是否互为伴随。如果两者用的差分模板不一致修正后的速度仍会残留散度长时间推进就会累积成质量不守恒。示例里通常用交错网格或同位网格加 Rhie-Chow 插值来保证这一点前者实现直观但索引繁琐后者索引简单但需要额外修正项。2.3 时间推进格式与 CFL 约束非定常问题的“非定常”三个字意味着时间精度和稳定性同样重要。显式对流项带来 CFL 约束时间步长必须小于网格尺度除以当地速度的量级否则对流项会发散。粘性项若也显式还会多一个扩散数约束在细网格上这个约束往往比 CFL 更苛刻。示例里常见的做法是对流显式、粘性隐式这样时间步主要由 CFL 决定而不是被粘性拖死。但隐式粘性意味着每一步要解一个线性系统MATLAB 里通常用稀疏矩阵加反斜杠求解或者用预条件迭代。时间推进用一阶欧拉还是二阶 Adams-Bashforth直接影响涡的耗散速度一阶格式数值耗散大涡会“糊”得快看起来稳定但物理上偏阻尼。判断一个示例是否可信先看它有没有在文档里交代时间格式的阶数。% 预测步显式对流 隐式粘性示意非完整可运行代码 % u, v 为当前速度dt 为时间步Re 为雷诺数dx, dy 为网格间距 [du_dx, du_dy] gradient(u, dx, dy); [dv_dx, dv_dy] gradient(v, dx, dy); conv_u u .* du_dx v .* du_dy; % 对流项显式 conv_v u .* dv_dx v .* dv_dy; % 粘性项组装为稀疏线性系统 A*u_star rhs隐式求解 A speye(Nx*Ny) - dt/Re * laplacian2d(Nx, Ny, dx, dy); u_star A \ (u(:) - dt * conv_u(:)); v_star A \ (v(:) - dt * conv_v(:));上面这段示意里laplacian2d负责组装二维五点拉普拉斯稀疏矩阵A的对角占优程度由dt/Re决定u_star、v_star是只满足动量方程的中间速度。注意对流项用的是当前时刻速度这是显式处理的标志粘性项进了左端矩阵所以时间步不受扩散数限制。实际示例里还要处理边界条件对矩阵的修改否则边界上的隐式求解会引入非物理反射。3. 离散细节网格、边界与散度残差3.1 交错网格与同位网格的取舍交错网格把压力放在单元中心、速度分量放在面上好处是压力梯度和散度天然使用相邻面速度离散散度严格为零压力泊松方程不会出现棋盘格振荡。代价是索引复杂每个速度分量有自己的网格尺寸写起来容易错位。同位网格把所有量放在同一套节点上索引统一、代码短但直接离散会遭遇压力-速度解耦出现棋盘格压力。常见的补救是 Rhie-Chow 动量插值在面速度里加一个压力梯度的修正项。示例如果用的是同位网格却没提 Rhie-Chow那压力场大概率有高频振荡只是被可视化平滑掉了。判断方法很简单把压力场画出来看有没有相邻节点正负跳变。3.2 边界条件的实现方式二维非定常示例里边界条件往往比内点离散更容易翻车。入口给定速度、出口给定压力或零梯度、壁面用无滑移这些在连续方程里很清晰落到离散上却要区分“强加”和“弱加”。强加是把边界节点值直接写入线性系统弱加是通过 ghost 节点或通量修正间接施加。投影法里尤其要注意压力泊松方程的边界条件必须和动量方程的边界条件相容否则修正步会在边界附近产生虚假的散度源。一个血泪经验是出口用零梯度压力时如果速度出口也强加零梯度两者叠加会让出口附近的质量通量不守恒涡在出口处被“憋”住。稳妥做法是出口压力固定、速度用对流流出条件。3.3 散度残差作为可信度指标示例跑完最该看的不是流线图而是每一步修正后的最大散度残差。理论上投影法修正后散度应为零实际受线性求解器容差和离散相容性影响会有一个小残差。如果这个残差随时间单调增长说明压力泊松方程没解到位或者散度算子与梯度算子不伴随。把残差曲线画出来比任何漂亮涡图都更能说明求解器是否健康。常见的实现思路是每步记录max(abs(divergence))和 CFL 数一起输出到日志长期跑的时候一眼就能看出哪一步开始失稳。% 每步修正后计算离散散度残差作为健康检查 div (u(2:end-1,2:end-1) - u(1:end-2,2:end-1)) / dx ... (v(2:end-1,2:end-1) - v(2:end-1,1:end-2)) / dy; max_div max(abs(div(:))); fprintf(step %d, CFL%.3f, max|div|%.3e\n, step, cfl, max_div);这段检查代码里div只算内点避免边界强加值干扰判断max_div的量级应该和线性求解器容差同阶如果突然跳升两三个数量级基本就是压力方程没收敛或边界条件打架。把它和 CFL 一起打印是排查非定常推进最省事的习惯。4. 避坑与排查非定常推进里最容易翻车的几处4.1 现象跑几十步后速度突然爆掉。原因CFL 按初始速度估的涡加速后超限。解决每步用当前最大速度重算 CFL超阈值就自动减时间步或报警别用固定步长硬扛。4.2 现象压力场出现棋盘格。原因同位网格直接离散缺 Rhie-Chow 修正。解决要么换交错网格要么在面速度里补压力梯度修正项并检查修正项系数是否与动量方程一致。4.3 现象出口附近涡被“吸”住不出去。原因出口速度强加零梯度与压力边界不相容。解决出口压力固定速度用一阶或二阶对流流出条件让涡自然流出。4.4 现象长时间跑总质量缓慢漂移。原因散度算子与梯度算子不伴随修正步残留散度累积。解决检查两者差分模板是否匹配必要时在修正步后加一次散度清理或改用严格相容的离散对。4.5 现象涡耗散得比物理预期快。原因一阶时间格式数值耗散大或对流项迎风离散过重。解决换二阶 Adams-Bashforth 或 Runge-Kutta对流项用中心差分加少量人工粘性而不是全迎风。5. 从示例到可复用骨架设计取舍与验证习惯把这份二维非定常 N-S 示例当教学玩具还是当骨架差别在于你有没有追问每一步的离散选择。投影法之所以流行是因为它把“满足约束”和“推进物理”解耦代价是压力泊松方程每步都要解计算量随网格规模上升。显式对流换来实现简单和并行友好代价是 CFL 约束隐式粘性换来大时间步代价是每步线性求解。交错网格换来严格无散度代价是索引复杂同位网格换来代码短代价是 Rhie-Chow 修正和棋盘格风险。这些取舍没有绝对优劣只有和你的雷诺数范围、网格规模、是否长期推进相匹配。验证一个非定常求解器我一般会走三步先用泰勒-格林涡的解析衰减对比数值衰减率看时间格式阶数是否对得上再跑方腔驱动流看涡心位置随雷诺数的迁移是否合理最后盯散度残差和动能随时间的变化动能应该单调衰减无外力时残差应该稳定在小量级。这三步里任何一步翻车都说明某个离散环节有问题而不是“网格不够细”能糊弄过去的。从那以后我每次拿到一个非定常 N-S 示例都强制先看它有没有输出散度残差和 CFL再看时间格式阶数写没写清楚最后才去看涡图画得好不好看。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/5 6:56:22 · 阅读 1.2k 非定常二维Navier-Stokes方程MATLAB实现：投影法、压力泊松方程与CFL约束
简介这份MATLAB示例面向流体力学学习者与数值计算入门者围绕二维非定常Navier-Stokes方程展开帮助读者理解不可压缩流体在二维空间内的运动求解流程。资源包共76个文件以44个m脚本为核心涵盖网格定义、边界条件设置、速度场与压力初始化、时间推进、空间离散、压力修正及可视化等环节另有29个png结果图、1个html说明页和1个txt许可文件压缩包约463KB结构紧凑便于按模块查阅。已有1947人学习下载说明其在教学与自学场景中具有一定参考价值。读者可从中获得有限差分与时间迭代的实现思路、压力泊松方程修正的代码组织方式以及quiver、contourf等可视化手段的调用范例适合用于课程设计、算法验证与编程练习也能加深对流体力学基本概念的理解。1. 非定常二维 Navier-Stokes为什么“能跑出漩涡”不等于“算对了”很多人第一次在 MATLAB 里跑二维非定常 Navier-Stokes 示例看到流场里卷出几个漂亮的涡就觉得大功告成。但真正决定结果可信度的往往不是那几个涡而是压力泊松方程怎么解、对流项怎么离散、时间推进用显式还是半隐式。二维非定常 N-S 方程的核心难点在于速度场必须每一步都满足不可压约束而压力本身没有独立的演化方程它只是让散度归零的拉格朗日乘子。这份 MATLAB 示例的价值就在于把“投影法”这条主线完整摊开从动量方程预测速度到解压力泊松方程修正速度再到时间步进与边界处理。搞懂它你才能判断一个示例是教学玩具还是可复用的求解骨架。2. 投影法拆解压力为什么是“事后补”的2.1 不可压约束与压力泊松方程的来源二维非定常不可压 N-S 方程写成无量纲形式是两条动量方程描述速度随时间的演化连续方程要求速度散度为零。麻烦在于动量方程里含压力梯度但没有任何一条方程直接告诉压力怎么随时间变化。常见的实现思路是把动量方程先“半离散”用当前时刻的速度算一个中间速度这个中间速度只满足动量方程、不满足散度为零然后用压力梯度去修正它让修正后的速度散度归零。把修正关系代入连续方程就得到压力泊松方程。这一步是整个投影法的枢纽——压力不是物理上独立演化的量而是为了让速度场满足约束而反解出来的。理解这一点后面所有离散选择都有了依据。2.2 中间速度的预测步与修正步预测步的离散形式取决于对流项和粘性项怎么处理。显式处理对流、隐式处理粘性是最常见的折中对流项用当前速度直接算粘性项用隐式格式避免时间步受粘性限制。修正步则是用解出的压力梯度去减中间速度。这里有个容易忽略的细节修正步之后的速度在离散意义下散度是否严格为零取决于压力泊松方程的离散算子和散度算子是否互为伴随。如果两者用的差分模板不一致修正后的速度仍会残留散度长时间推进就会累积成质量不守恒。示例里通常用交错网格或同位网格加 Rhie-Chow 插值来保证这一点前者实现直观但索引繁琐后者索引简单但需要额外修正项。2.3 时间推进格式与 CFL 约束非定常问题的“非定常”三个字意味着时间精度和稳定性同样重要。显式对流项带来 CFL 约束时间步长必须小于网格尺度除以当地速度的量级否则对流项会发散。粘性项若也显式还会多一个扩散数约束在细网格上这个约束往往比 CFL 更苛刻。示例里常见的做法是对流显式、粘性隐式这样时间步主要由 CFL 决定而不是被粘性拖死。但隐式粘性意味着每一步要解一个线性系统MATLAB 里通常用稀疏矩阵加反斜杠求解或者用预条件迭代。时间推进用一阶欧拉还是二阶 Adams-Bashforth直接影响涡的耗散速度一阶格式数值耗散大涡会“糊”得快看起来稳定但物理上偏阻尼。判断一个示例是否可信先看它有没有在文档里交代时间格式的阶数。% 预测步显式对流 隐式粘性示意非完整可运行代码 % u, v 为当前速度dt 为时间步Re 为雷诺数dx, dy 为网格间距 [du_dx, du_dy] gradient(u, dx, dy); [dv_dx, dv_dy] gradient(v, dx, dy); conv_u u .* du_dx v .* du_dy; % 对流项显式 conv_v u .* dv_dx v .* dv_dy; % 粘性项组装为稀疏线性系统 A*u_star rhs隐式求解 A speye(Nx*Ny) - dt/Re * laplacian2d(Nx, Ny, dx, dy); u_star A \ (u(:) - dt * conv_u(:)); v_star A \ (v(:) - dt * conv_v(:));上面这段示意里laplacian2d负责组装二维五点拉普拉斯稀疏矩阵A的对角占优程度由dt/Re决定u_star、v_star是只满足动量方程的中间速度。注意对流项用的是当前时刻速度这是显式处理的标志粘性项进了左端矩阵所以时间步不受扩散数限制。实际示例里还要处理边界条件对矩阵的修改否则边界上的隐式求解会引入非物理反射。3. 离散细节网格、边界与散度残差3.1 交错网格与同位网格的取舍交错网格把压力放在单元中心、速度分量放在面上好处是压力梯度和散度天然使用相邻面速度离散散度严格为零压力泊松方程不会出现棋盘格振荡。代价是索引复杂每个速度分量有自己的网格尺寸写起来容易错位。同位网格把所有量放在同一套节点上索引统一、代码短但直接离散会遭遇压力-速度解耦出现棋盘格压力。常见的补救是 Rhie-Chow 动量插值在面速度里加一个压力梯度的修正项。示例如果用的是同位网格却没提 Rhie-Chow那压力场大概率有高频振荡只是被可视化平滑掉了。判断方法很简单把压力场画出来看有没有相邻节点正负跳变。3.2 边界条件的实现方式二维非定常示例里边界条件往往比内点离散更容易翻车。入口给定速度、出口给定压力或零梯度、壁面用无滑移这些在连续方程里很清晰落到离散上却要区分“强加”和“弱加”。强加是把边界节点值直接写入线性系统弱加是通过 ghost 节点或通量修正间接施加。投影法里尤其要注意压力泊松方程的边界条件必须和动量方程的边界条件相容否则修正步会在边界附近产生虚假的散度源。一个血泪经验是出口用零梯度压力时如果速度出口也强加零梯度两者叠加会让出口附近的质量通量不守恒涡在出口处被“憋”住。稳妥做法是出口压力固定、速度用对流流出条件。3.3 散度残差作为可信度指标示例跑完最该看的不是流线图而是每一步修正后的最大散度残差。理论上投影法修正后散度应为零实际受线性求解器容差和离散相容性影响会有一个小残差。如果这个残差随时间单调增长说明压力泊松方程没解到位或者散度算子与梯度算子不伴随。把残差曲线画出来比任何漂亮涡图都更能说明求解器是否健康。常见的实现思路是每步记录max(abs(divergence))和 CFL 数一起输出到日志长期跑的时候一眼就能看出哪一步开始失稳。% 每步修正后计算离散散度残差作为健康检查 div (u(2:end-1,2:end-1) - u(1:end-2,2:end-1)) / dx ... (v(2:end-1,2:end-1) - v(2:end-1,1:end-2)) / dy; max_div max(abs(div(:))); fprintf(step %d, CFL%.3f, max|div|%.3e\n, step, cfl, max_div);这段检查代码里div只算内点避免边界强加值干扰判断max_div的量级应该和线性求解器容差同阶如果突然跳升两三个数量级基本就是压力方程没收敛或边界条件打架。把它和 CFL 一起打印是排查非定常推进最省事的习惯。4. 避坑与排查非定常推进里最容易翻车的几处4.1 现象跑几十步后速度突然爆掉。原因CFL 按初始速度估的涡加速后超限。解决每步用当前最大速度重算 CFL超阈值就自动减时间步或报警别用固定步长硬扛。4.2 现象压力场出现棋盘格。原因同位网格直接离散缺 Rhie-Chow 修正。解决要么换交错网格要么在面速度里补压力梯度修正项并检查修正项系数是否与动量方程一致。4.3 现象出口附近涡被“吸”住不出去。原因出口速度强加零梯度与压力边界不相容。解决出口压力固定速度用一阶或二阶对流流出条件让涡自然流出。4.4 现象长时间跑总质量缓慢漂移。原因散度算子与梯度算子不伴随修正步残留散度累积。解决检查两者差分模板是否匹配必要时在修正步后加一次散度清理或改用严格相容的离散对。4.5 现象涡耗散得比物理预期快。原因一阶时间格式数值耗散大或对流项迎风离散过重。解决换二阶 Adams-Bashforth 或 Runge-Kutta对流项用中心差分加少量人工粘性而不是全迎风。5. 从示例到可复用骨架设计取舍与验证习惯把这份二维非定常 N-S 示例当教学玩具还是当骨架差别在于你有没有追问每一步的离散选择。投影法之所以流行是因为它把“满足约束”和“推进物理”解耦代价是压力泊松方程每步都要解计算量随网格规模上升。显式对流换来实现简单和并行友好代价是 CFL 约束隐式粘性换来大时间步代价是每步线性求解。交错网格换来严格无散度代价是索引复杂同位网格换来代码短代价是 Rhie-Chow 修正和棋盘格风险。这些取舍没有绝对优劣只有和你的雷诺数范围、网格规模、是否长期推进相匹配。验证一个非定常求解器我一般会走三步先用泰勒-格林涡的解析衰减对比数值衰减率看时间格式阶数是否对得上再跑方腔驱动流看涡心位置随雷诺数的迁移是否合理最后盯散度残差和动能随时间的变化动能应该单调衰减无外力时残差应该稳定在小量级。这三步里任何一步翻车都说明某个离散环节有问题而不是“网格不够细”能糊弄过去的。从那以后我每次拿到一个非定常 N-S 示例都强制先看它有没有输出散度残差和 CFL再看时间格式阶数写没写清楚最后才去看涡图画得好不好看。希望帮到你。本文还有配套的精品资源点击获取