简介二维非定常 Navier-Stokes 方程的 MATLAB 实现示例面向流体力学初学者、CFD 课程学生以及需要开展流场仿真的工程技术人员。资源以有限差分法与压力修正算法为主线覆盖网格生成、边界条件施加、时间推进、泊松方程求解等典型环节并通过多种辅助函数构建出清晰的模块化程序框架便于逐段对照学习。压缩包共七十六个文件包括四十四份 MATLAB 源程序、二十九张过程结果图、一份许可文本和一份超文本说明文档整包大小约四百六十三 KB轻巧易取。二十九张图像展示不同时刻的速度场与压力场演化超文本文档辅助解释代码结构帮助读者快速掌握非定常流动模拟的关键步骤。已有约一千九百四十一人学习浏览适合用来巩固有限差分、压力修正等数值方法并可在修改边界条件或参数后拓展至更多流动场景。1. 一个能跑的MATLAB二维非定常Navier-Stokes示例从哪来、能算啥、坑在哪如果你在找的是那种“粘过去就能出图”的MATLAB二维非定常Navier-Stokes示例那我先把结论说清楚这份资源跑通之后你能得到的是方腔顶盖驱动流从启动到定常的完整演化过程以及一个自带解析解的衰减涡算例。它不算新用的也是很多人觉得“过时”的有限差分加投影分步法但恰恰是这套朴素的组合把非定常不可压缩流动里最难啃的压力耦合问题拆成了三步适合拿来当课程设计、毕业论文前期的验证基准。适合谁呢正在学计算流体力学、需要在MATLAB里从零搭NS求解器、或者是想快速验证自己离散格式有没有写错的人。我见过不少人下载后跑出NaN就放弃了绝大多数不是代码问题而是CFL条件、边界索引、迭代容差这三个地方没对齐。这份示例真正值钱的不是那几行MATLAB代码而是它替你标好了所有需要对齐的坐标。2. 从连续方程到离散格式交错网格与CFL条件怎么定2.1 为什么选“不可压缩二维显式时间推进”这个组合这份MATLAB示例求解的是二维不可压缩Navier-Stokes方程连续方程是速度散度为零动量方程是经典的对流项、压力梯度项和粘性项。二维意味着只有u、v两个速度分量外加一个压力p。不可压缩的难点不在于未知数多而在于压力没有独立的演化方程它是以“散度为零”这个约束形式出现。直接同时解三个场会遇到鞍点问题常规做法都是把压力处理成通过泊松方程求解的拉格朗日乘子。这套示例选的是显式时间推进也就是每个新时刻的速度都由上一时刻的已知量直接算出来。显式的好处是每步代价低、代码和直觉一致坏处是时间步长受稳定性约束。二维低雷诺数流场显式完全够用如果你以后要跑高雷诺数湍流再考虑隐式或压力投影的隐式版本也不迟。方程形式如下非定常体现在左边第一项时间导数$$\frac{\partial u}{\partial t} u\frac{\partial u}{\partial x} v\frac{\partial u}{\partial y} -\frac{\partial p}{\partial x} \frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} \frac{\partial^2 u}{\partial y^2}\right)$$代码里不会直接写成这种符号形式而是拆成有限差分。这里的关键是数值上满足连续方程的速度场和满足动量方程的速度场往往不一致所以需要投影分步法在每一步里“修正”一次。2.2 交错网格为什么压力和速度不能住在同一个格点上如果让u、v、p都定义在同一套网格节点上压力泊松方程的中心差分会允许一种棋盘状的振荡压力场存在正负压力交替排列但梯度为零导致动量方程根本感受不到它这种“伪模态”在迭代解泊松方程时特别容易被放大。使用交错网格是常见做法也就是把速度分量放在网格单元边界上压力放在单元中心。这份示例里典型的数组尺寸是这样分配的如果你看到u、v、p三个矩阵尺寸完全一样那就要警惕棋盘振荡% 计算域与网格参数 Lx 1; % x方向长度 Ly 1; % y方向长度 nx 64; % x方向控制体数量 ny 64; % y方向控制体数量 dx Lx / nx; dy Ly / ny; % 交错网格数组 % p 定义在单元中心尺寸为 nx x ny % u 定义在竖直单元边上尺寸为 (nx1) x ny % v 定义在水平单元边上尺寸为 nx x (ny1) u zeros(nx1, ny); % 水平速度分量 v zeros(nx, ny1); % 竖直速度分量 p zeros(nx, ny); % 压力这段初始化把网格数、计算域尺寸和三个场的存放位置全部框死了。网格坐标写法是x_c (0.5:nx)*dx意思是压力节点位于每个单元的中心速度分量相应地错开半个网格。索引错位是初学交错网格最容易翻车的地方建议在写迭代代码前先画一张小网格把 u(i,j)、v(i,j)、p(i,j) 各自对应的物理位置标出来再动手写循环。参数 Lx、Ly、nx、ny 决定了空间分辨率Re 决定流动特性这四个参数在后面的算例里会反复调整。2.3 dt 不是随便取的CFL 与粘性稳定性的合并约束显式时间推进有两个稳定性约束必须同时满足。第一个是对流CFL条件信息在一个时间步内不能跨过一个网格单元量化后是u*dt/dx和v*dt/dy都必须小于某个安全系数通常取0.25到0.5。第二个是显式扩散项的稳定性约束这个容易被忽略粘性项的显式处理要求2*dt/(Re*dx^2)小于1也就是dt Re*dx^2/2。很多人只盯着CFL结果低雷诺数、细网格时照样发散。综合判断时我一般会写成一个动态计算dt的代码而不是手工试% 动态时间步长同时满足对流与粘性稳定性 umax max(abs(u(:))); vmax max(abs(v(:))); % 对流限制 dt_cfl min(dx/umax, dy/vmax) * 0.25; % 粘性限制显式扩散稳定约束 dt_diff 0.5 * Re * min(dx^2, dy^2); % 实际采用两者中较小的 dt min(dt_cfl, dt_diff);这里的0.25是安全系数如果流场中出现局部高速区这个系数要再降低。dt_diff的公式来自一维扩散方程的von Neumann稳定性分析二维时取两个方向里更严格的那个即可。算完之后最好把dt打印出来看一眼如果比初始猜测值小一个数量级那是正常的别硬顶着大时间步跑会直接爆掉。3. 投影法三步走压力泊松方程的MATLAB求解细节3.1 第一步忽略压力梯度推进中间速度场投影法的第一步是先去掉动量方程里的压力梯度项用当前时刻的u、v计算一个新的中间速度场通常记为u星。这一步物理上不完整因为还没满足连续方程但它已经包含了对流和粘性的贡献。动量方程去掉压力项后u的预测方程为$$u^* u^n dt \left[ -\left(u\frac{\partial u}{\partial x} v\frac{\partial u}{\partial y}\right) \frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} \frac{\partial^2 u}{\partial y^2}\right) \right]$$MATLAB里的核心循环长这样% 预测步计算中间速度 ustar、vstar % 以 u 分量为例v 分量对称处理 % 粘性项二阶中心差分 uxx (u(3:nx1, :) - 2*u(2:nx, :) u(1:nx-1, :)) / dx^2; uyy (u(:, 2:ny) - 2*u(:, 1:ny-1) u(:, 1:ny-1)) / dy^2; % 注意边界索引 visc_u (uxx uyy) / Re; % 对流项一阶迎风先判断速度方向 % u*w 的乘积在交错网格上需要插值到同一位置 % 这里用简单的一阶迎风示意完整代码里会做半格插值 conv_u zeros(nx1, ny); for i 2:nx for j 2:ny-1 ue 0.5*(u(i1,j) u(i,j)); uw 0.5*(u(i,j) u(i-1,j)); un 0.5*(u(i,j1) u(i,j)); us 0.5*(u(i,j) u(i,j-1)); conv_u(i,j) ue*max(ue,0) - uw*max(uw,0) ... un*max(un,0) - us*max(us,0); end end % 中间速度 ustar(2:nx, 2:ny-1) u(2:nx, 2:ny-1) dt*(-conv_u(2:nx,2:ny-1) visc_u);这段代码里最需要注意的就是交错网格上的乘积插值对流项u*du/dx里的速度相乘必须在同一位置取值否则会产生数值不稳定。注释里“半格插值”的意思是相邻两个速度的平均。一阶迎风是低雷诺数算例的稳妥选择虽然耗散偏大但不会像中心差分那样在间断处产生振荡。如果你想要二阶精度可以把迎风换成二阶迎风或者Quick格式代价是实现复杂度上升。3.2 第二步解压力泊松方程这一步的坑最多中间速度u星不满足散度为零我们要找一个压力场用它的梯度去修正速度使得修正后的速度场无散度。推导过程很经典对修正公式两边取散度要求div(u^{n1})0就得到压力泊松方程$$\nabla^2 p \frac{1}{dt} \nabla \cdot u^*$$这个方程在代码里不会去组装大型稀疏矩阵而是用Jacobi迭代逐点更新。原因有两个一是迭代法在64乘64网格上收敛很快代码也直观二是避免在MATLAB里维护稀疏矩阵存储。代价是迭代次数和收敛容差需要手工控制。核心代码% 计算散度 RHS div_u zeros(nx, ny); div_u(2:nx-1, 2:ny-1) (ustar(3:nx, 2:ny-1) - ustar(2:nx-1, 2:ny-1)) / dx ... (vstar(2:nx-1, 3:ny) - vstar(2:nx-1, 2:ny-1)) / dy; rhs div_u / dt; % Jacobi 迭代求解压力泊松方程 p zeros(nx, ny); tol 1e-6; maxiter 2000; for iter 1:maxiter p_old p; % 内部节点五点差分 p(2:nx-1, 2:ny-1) 0.25 * ( ... p_old(1:nx-2, 2:ny-1) p_old(3:nx, 2:ny-1) ... p_old(2:nx-1, 1:ny-2) p_old(2:nx-1, 3:ny) - ... dx^2 * rhs(2:nx-1, 2:ny-1) ); % 边界Neumann 零梯度条件 p(1,:) p(2,:); p(nx,:) p(nx-1,:); p(:,1) p(:,2); p(:,ny) p(:,ny-1); if max(abs(p(:) - p_old(:))) tol break; end endrhs的计算是整个投影法的心脏这里散度差分必须和后面速度修正里的压力梯度差分使用完全相同的网格位置。如果散度用(u(i1)-u(i))/dx梯度也必须用对应的交错位置混用会导致修正后仍然有残余散度。Jacobi迭代的好处是天然满足Neumann边界条件边界处直接让边界压力等于相邻内点压力。tol取1e-6在大多数算例够了如果你观察流场有缓慢的周期波动可以调到1e-8。3.3 第三步修正速度让散度严格为零得到压力场后修正步用一个简单的梯度减法完成$$u^{n1} u^* - dt \frac{\partial p}{\partial x}, \quad v^{n1} v^* - dt \frac{\partial p}{\partial y}$$这一步做完之后理论上速度场就满足连续方程了。MATLAB实现很短% 修正步压力梯度修正速度场 u(2:nx, 2:ny-1) ustar(2:nx, 2:ny-1) - dt / dx * (p(2:nx, 2:ny-1) - p(1:nx-1, 2:ny-1)); v(2:nx-1, 2:ny) vstar(2:nx-1, 2:ny) - dt / dy * (p(2:nx-1, 2:ny) - p(2:nx-1, 1:ny-1)); % 主时间循环骨架 for n 1:nt % 1. 边界条件施加固壁无滑移或解析解边界 apply_boundary_conditions(u, v); % 2. 预测步得到 ustar, vstar % 3. 求解压力泊松方程 % 4. 修正速度场 % 5. 记录输出速度范数、残差、动画 end每步循环顺序不能乱先施加边界再预测、投影、修正。如果边界条件在预测步之后才施加边界上会积累误差。非定常的意义在于主循环里的每一步都对应真实的物理时间每一个中间场都是物理上可解释的流动状态。衰减涡算例验证的正是这一点同一时刻的数值解必须无限逼近解析解而不是只求“最终定常结果正确”。4. 边界条件与两个标准算例把代码真正跑起来4.1 顶盖驱动流非定常推进到定常的经典路径顶盖驱动流是所有NS代码的标配验证算例。计算域是单位方腔顶盖沿x方向以恒定速度运动其余三壁固定。雷诺数的定义是Re U0 * L / nu其中U0是顶盖速度L是方腔边长。这份示例里就取U0等于1、L等于1所以Re 1/nu代码里直接通过Re控制粘性系数大小。边界条件的实现是一行行硬编码的% 顶盖驱动流边界 % 下壁、左壁、右壁无滑移 u0, v0 u(:, 1) 0; % 下壁 u(:, ny) 1; % 顶盖水平速度设为 1 v(:, 1) 0; % 下壁 v(:, ny) 0; % 顶盖竖直速度为零 % 左壁与右壁 u(1, :) 0; u(nx1, :) 0; v(1, :) 0; v(nx1, :) 0;跑顶盖驱动流时有个容易误解的点题目叫非定常但顶盖驱动流的经典验证目标是“定常解”。启动瞬间整个流场是一潭死水的顶盖运动后涡逐渐发展经历一个过渡过程最终达到定常状态。所以这份示例既能展示非定常演化又能通过定常终态验证收敛性。判断定常别用“跑够多少步”而是看速度场的变化量% 定常收敛判据连续两步速度场变化小于阈值 diff_u max(abs(u(:) - u_prev(:))); diff_v max(abs(v(:) - v_prev(:))); if max(diff_u, diff_v) 1e-8 break; endRe100时粗网格跑到定常大概需要几千步细网格上万步。如果Re提高到1000流场会呈现出更复杂的二次涡结构收敛也更慢这也是后续练手的好对象。4.2 衰减涡算例用解析解自校验时间精度顶盖驱动流只能验证定常结果对不对验证非定常时间积分的正确性需要另一个算例衰减涡。它有一个解析解在方腔内用正弦函数构造一个涡其强度随时间指数衰减。这个算例的最大价值是能把每个时刻的数值解和精确解对比直接算出误差阶。衰减涡的解析速度场为% 衰减涡初始场正弦涡 % 解析解u sin(pi*x) * cos(pi*y) * exp(-2*pi^2*t/Re) % v -cos(pi*x) * sin(pi*y) * exp(-2*pi^2*t/Re) x_c ((1:nx)-0.5) * dx; y_c ((1:ny)-0.5) * dy; [Xc, Yc] meshgrid(x_c, y_c); % 初始时刻 t0 u sin(pi*Xc) .* cos(pi*Yc); v -cos(pi*Xc) .* sin(pi*Yc);注意初始场给的是单元中心位置而u、v数组在交错网格上分别位于不同位置所以要对x、y坐标做适当的偏移这又是一个索引细节。误差计算在任意时刻做% 在给定时刻对比数值解与解析解 t n * dt; u_exact sin(pi*Xc) .* cos(pi*Yc) .* exp(-2*pi^2*t/Re); err_u sqrt(sum((u(2:nx1, :) - u_exact(:,:)).^2) / numel(u_exact));这里我习惯用相对L2误差来判断而不是看最大误差因为最大误差对个别点太敏感。跑衰减涡时边界条件不能用固壁无滑移而应该把解析解直接赋在边界上这样误差只来自内部离散不会被边界处理污染。4.3 后处理看见流场才算真正跑通跑完计算后至少看三个图速度矢量图、压力云图、流线图。MATLAB里没有特别多技巧核心命令就这几个% 视觉化流场 figure quiver(x_u, y_u, u, v, 2); % 速度矢量2是缩放比例 hold on contourf(x_c, y_c, p, 20); % 压力云图 axis equal title(二维非定常方腔流瞬时流场); % 导出数据便于后续分析 save(cavity_result.mat, u, v, p, t, Re, nx, ny, dx, dy, dt);速度矢量图能直观看到主涡和角涡压力云图能看出高低压区分布流线图适合展示涡心位置。保存为后缀为mat的文件后之后再想画中心线速度剖面、算涡量、做动画都不用重跑一遍计算这点很重要因为在细网格上重跑一次的代价不小。5. 避坑指南发散、不收敛与伪非定常的五个典型现场5.1 前几步就爆NaN十有八九是dt超限现象程序启动后前几个时间步速度场就出现Inf或NaN绘图时报错整个计算直接终止。原因时间步太长超过了上一章说的对流或粘性稳定性边界。尤其是低雷诺数加细网格时粘性约束项Re*dx^2/2往往比CFL约束更严格只按CFL取dt必爆。另一个常见原因是初始速度场里存在不合理的局部高速动态dt计算时被低估。解决把dt计算改成动态的每步都重新扫描全场最大速度并同时套用对流和粘性两个约束。如果仍然爆炸把安全系数从0.25降到0.1试一次。这个问题的特征是“换个dt就从爆变稳”不用怀疑代码逻辑。5.2 压力场棋盘振荡检查你的网格是不是“真交错”现象速度场看起来大体合理但压力云图呈密集的红黑棋盘格交替数值上相邻格点压力正负交替出现。原因压力棋盘振荡是泊松方程零能模态的表现通常因为u、v、p数组尺寸相同也就是名义上“交错”实际没交错或者是散度和梯度算子用的差分方向不一致导致压力修正没有作用到速度的正确位置。解决先确认数组尺寸满足size(p)[nx,ny]、size(u)[nx1,ny]、size(v)[nx,ny1]。再检查修正步里的压力梯度索引u(2:nx,...)对应p(2:nx,...)-p(1:nx-1,...)如果写成了p(3:nx1,...)-p(2:nx,...)压力梯度会整体平移一格相当于用错了离散模板。这两点检查完棋盘振荡基本消失。5.3 顶盖驱动流左上角压力爆表角点奇点处理不当现象顶盖驱动流跑到几百步后左上角和右上角附近的压力异常大速度等值线在角点处剧烈扭曲。原因顶盖水平速度是1左壁水平速度是0这两个边界条件在角点直接冲突形成了数学上的奇点。网格越细角点处的速度梯度和压力越尖锐。解决在角点处对顶盖速度做一个平滑渐变常见做法是把角点速度设为零或者用一段余弦渐变过渡。代码里就是一行的事u(nx1, ny) 0; % 右上角速度置零消除奇点 u(1, ny) 0; % 左上角速度置零加这一行之后全局压力最大值会明显回落。这个改动只影响角点附近一个网格的精度对整体流场影响很小。5.4 衰减涡误差不降反升边界条件施加错了现象网格从32加密到64再到128误差几乎不变甚至在增大完全不符合二阶格式应有的表现。原因衰减涡算例的边界上如果把固壁条件当边界用边界附近的数值解会和解析解完全脱节误差由边界主导加密网格也救不回来。另一个原因是时间步长没同步缩小空间二阶但时间一阶时间误差盖过了空间误差。解决衰减涡的所有边界都直接赋解析解的值用当前位置和当前时刻的精确公式计算边界上的u、v。同时把dt减半和网格加密两件事分开做先固定dt只加密网格测空间精度再固定网格只减dt测时间精度就能定位到底是哪一阶拖了后腿。5.5 看着在动其实是伪非定常收敛判据太乐观现象顶盖驱动流跑了预设步数流场每一时刻都在缓慢变化看起来好像有非定常特征但速度变化量一直降不下去始终在某个水平抖动。原因两种情况。一种是收敛阈值设得太宽比如1e-4就停流场还没真正平衡另一种是压力泊松迭代的tol太松每步的残余散度没有被完全压掉导致速度场被反复注入微小扰动。解决把定常判据改成连续两步速度场最大差小于1e-8并让压力泊松迭代的tol也降到1e-8再看曲线是否平滑收敛。如果残差在某个数值附近震荡不降切割压力泊松迭代次数通常能直接见效。记住一个习惯判断定常要看历史曲线不要当场拍脑袋。6. 进阶从复现走向验证的快速检查法拿到这份MATLAB示例后我建议你先别急着改物理参数而是花半天时间做三件验证都很便宜。第一件是衰减涡的时间收敛阶测试。固定网格在64乘64从dt0.01开始逐步减半记录每一档的L2误差看误差是否大约按dt的第一次方下降。如果是说明时间推进和投影法配合正确如果误差不随dt变化说明时间误差不是主导问题可能在边界或初始场。dt速度L2误差比上一档变化的倍数0.012.1e-3-0.0051.05e-3约1/20.00255.2e-4约1/2第二件是顶盖驱动流的速度剖面验证。在定常状态下取出方腔水平中线上u速度随y的分布和同雷诺数的公开基准数据对比。Re100时中心处回流的负速度峰值大约在-0.2附近如果你的结果差出10%以上优先怀疑网格分辨率或角点平滑的影响。第三件是质量守恒检查。计算任意时刻全场散度的平方和这个数应当在压力泊松迭代的容差量级如果明显偏大说明投影法第三步的梯度索引写错了。性能优化方面整套代码最大的瓶颈在压力泊松方程的Jacobi迭代可以用MATLAB矩阵操作代替双重for循环加速比通常在5倍以上。如果网格超过128乘128把迭代改成预条件共轭梯度法会更划算但要注意边界条件的组装方式不同。这份示例包自带的算例控制在几十秒内跑完没必要为了性能牺牲代码可读性。从那以后我每次拿到任何NS求解代码都会强制自己先跑一遍衰减涡误差测试再谈别的。误差阶不对后面的流场图再漂亮也是废的。希望这份MATLAB二维非定常示例能帮你少走这几步弯路早日把注意力放到真正要做的流动问题上。本文还有配套的精品资源点击获取 
长河君 · 2026/10/23 11:20:14 · 阅读 1.2k MATLAB求解二维非定常Navier-Stokes：投影法与交错网格
简介二维非定常 Navier-Stokes 方程的 MATLAB 实现示例面向流体力学初学者、CFD 课程学生以及需要开展流场仿真的工程技术人员。资源以有限差分法与压力修正算法为主线覆盖网格生成、边界条件施加、时间推进、泊松方程求解等典型环节并通过多种辅助函数构建出清晰的模块化程序框架便于逐段对照学习。压缩包共七十六个文件包括四十四份 MATLAB 源程序、二十九张过程结果图、一份许可文本和一份超文本说明文档整包大小约四百六十三 KB轻巧易取。二十九张图像展示不同时刻的速度场与压力场演化超文本文档辅助解释代码结构帮助读者快速掌握非定常流动模拟的关键步骤。已有约一千九百四十一人学习浏览适合用来巩固有限差分、压力修正等数值方法并可在修改边界条件或参数后拓展至更多流动场景。1. 一个能跑的MATLAB二维非定常Navier-Stokes示例从哪来、能算啥、坑在哪如果你在找的是那种“粘过去就能出图”的MATLAB二维非定常Navier-Stokes示例那我先把结论说清楚这份资源跑通之后你能得到的是方腔顶盖驱动流从启动到定常的完整演化过程以及一个自带解析解的衰减涡算例。它不算新用的也是很多人觉得“过时”的有限差分加投影分步法但恰恰是这套朴素的组合把非定常不可压缩流动里最难啃的压力耦合问题拆成了三步适合拿来当课程设计、毕业论文前期的验证基准。适合谁呢正在学计算流体力学、需要在MATLAB里从零搭NS求解器、或者是想快速验证自己离散格式有没有写错的人。我见过不少人下载后跑出NaN就放弃了绝大多数不是代码问题而是CFL条件、边界索引、迭代容差这三个地方没对齐。这份示例真正值钱的不是那几行MATLAB代码而是它替你标好了所有需要对齐的坐标。2. 从连续方程到离散格式交错网格与CFL条件怎么定2.1 为什么选“不可压缩二维显式时间推进”这个组合这份MATLAB示例求解的是二维不可压缩Navier-Stokes方程连续方程是速度散度为零动量方程是经典的对流项、压力梯度项和粘性项。二维意味着只有u、v两个速度分量外加一个压力p。不可压缩的难点不在于未知数多而在于压力没有独立的演化方程它是以“散度为零”这个约束形式出现。直接同时解三个场会遇到鞍点问题常规做法都是把压力处理成通过泊松方程求解的拉格朗日乘子。这套示例选的是显式时间推进也就是每个新时刻的速度都由上一时刻的已知量直接算出来。显式的好处是每步代价低、代码和直觉一致坏处是时间步长受稳定性约束。二维低雷诺数流场显式完全够用如果你以后要跑高雷诺数湍流再考虑隐式或压力投影的隐式版本也不迟。方程形式如下非定常体现在左边第一项时间导数$$\frac{\partial u}{\partial t} u\frac{\partial u}{\partial x} v\frac{\partial u}{\partial y} -\frac{\partial p}{\partial x} \frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} \frac{\partial^2 u}{\partial y^2}\right)$$代码里不会直接写成这种符号形式而是拆成有限差分。这里的关键是数值上满足连续方程的速度场和满足动量方程的速度场往往不一致所以需要投影分步法在每一步里“修正”一次。2.2 交错网格为什么压力和速度不能住在同一个格点上如果让u、v、p都定义在同一套网格节点上压力泊松方程的中心差分会允许一种棋盘状的振荡压力场存在正负压力交替排列但梯度为零导致动量方程根本感受不到它这种“伪模态”在迭代解泊松方程时特别容易被放大。使用交错网格是常见做法也就是把速度分量放在网格单元边界上压力放在单元中心。这份示例里典型的数组尺寸是这样分配的如果你看到u、v、p三个矩阵尺寸完全一样那就要警惕棋盘振荡% 计算域与网格参数 Lx 1; % x方向长度 Ly 1; % y方向长度 nx 64; % x方向控制体数量 ny 64; % y方向控制体数量 dx Lx / nx; dy Ly / ny; % 交错网格数组 % p 定义在单元中心尺寸为 nx x ny % u 定义在竖直单元边上尺寸为 (nx1) x ny % v 定义在水平单元边上尺寸为 nx x (ny1) u zeros(nx1, ny); % 水平速度分量 v zeros(nx, ny1); % 竖直速度分量 p zeros(nx, ny); % 压力这段初始化把网格数、计算域尺寸和三个场的存放位置全部框死了。网格坐标写法是x_c (0.5:nx)*dx意思是压力节点位于每个单元的中心速度分量相应地错开半个网格。索引错位是初学交错网格最容易翻车的地方建议在写迭代代码前先画一张小网格把 u(i,j)、v(i,j)、p(i,j) 各自对应的物理位置标出来再动手写循环。参数 Lx、Ly、nx、ny 决定了空间分辨率Re 决定流动特性这四个参数在后面的算例里会反复调整。2.3 dt 不是随便取的CFL 与粘性稳定性的合并约束显式时间推进有两个稳定性约束必须同时满足。第一个是对流CFL条件信息在一个时间步内不能跨过一个网格单元量化后是u*dt/dx和v*dt/dy都必须小于某个安全系数通常取0.25到0.5。第二个是显式扩散项的稳定性约束这个容易被忽略粘性项的显式处理要求2*dt/(Re*dx^2)小于1也就是dt Re*dx^2/2。很多人只盯着CFL结果低雷诺数、细网格时照样发散。综合判断时我一般会写成一个动态计算dt的代码而不是手工试% 动态时间步长同时满足对流与粘性稳定性 umax max(abs(u(:))); vmax max(abs(v(:))); % 对流限制 dt_cfl min(dx/umax, dy/vmax) * 0.25; % 粘性限制显式扩散稳定约束 dt_diff 0.5 * Re * min(dx^2, dy^2); % 实际采用两者中较小的 dt min(dt_cfl, dt_diff);这里的0.25是安全系数如果流场中出现局部高速区这个系数要再降低。dt_diff的公式来自一维扩散方程的von Neumann稳定性分析二维时取两个方向里更严格的那个即可。算完之后最好把dt打印出来看一眼如果比初始猜测值小一个数量级那是正常的别硬顶着大时间步跑会直接爆掉。3. 投影法三步走压力泊松方程的MATLAB求解细节3.1 第一步忽略压力梯度推进中间速度场投影法的第一步是先去掉动量方程里的压力梯度项用当前时刻的u、v计算一个新的中间速度场通常记为u星。这一步物理上不完整因为还没满足连续方程但它已经包含了对流和粘性的贡献。动量方程去掉压力项后u的预测方程为$$u^* u^n dt \left[ -\left(u\frac{\partial u}{\partial x} v\frac{\partial u}{\partial y}\right) \frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} \frac{\partial^2 u}{\partial y^2}\right) \right]$$MATLAB里的核心循环长这样% 预测步计算中间速度 ustar、vstar % 以 u 分量为例v 分量对称处理 % 粘性项二阶中心差分 uxx (u(3:nx1, :) - 2*u(2:nx, :) u(1:nx-1, :)) / dx^2; uyy (u(:, 2:ny) - 2*u(:, 1:ny-1) u(:, 1:ny-1)) / dy^2; % 注意边界索引 visc_u (uxx uyy) / Re; % 对流项一阶迎风先判断速度方向 % u*w 的乘积在交错网格上需要插值到同一位置 % 这里用简单的一阶迎风示意完整代码里会做半格插值 conv_u zeros(nx1, ny); for i 2:nx for j 2:ny-1 ue 0.5*(u(i1,j) u(i,j)); uw 0.5*(u(i,j) u(i-1,j)); un 0.5*(u(i,j1) u(i,j)); us 0.5*(u(i,j) u(i,j-1)); conv_u(i,j) ue*max(ue,0) - uw*max(uw,0) ... un*max(un,0) - us*max(us,0); end end % 中间速度 ustar(2:nx, 2:ny-1) u(2:nx, 2:ny-1) dt*(-conv_u(2:nx,2:ny-1) visc_u);这段代码里最需要注意的就是交错网格上的乘积插值对流项u*du/dx里的速度相乘必须在同一位置取值否则会产生数值不稳定。注释里“半格插值”的意思是相邻两个速度的平均。一阶迎风是低雷诺数算例的稳妥选择虽然耗散偏大但不会像中心差分那样在间断处产生振荡。如果你想要二阶精度可以把迎风换成二阶迎风或者Quick格式代价是实现复杂度上升。3.2 第二步解压力泊松方程这一步的坑最多中间速度u星不满足散度为零我们要找一个压力场用它的梯度去修正速度使得修正后的速度场无散度。推导过程很经典对修正公式两边取散度要求div(u^{n1})0就得到压力泊松方程$$\nabla^2 p \frac{1}{dt} \nabla \cdot u^*$$这个方程在代码里不会去组装大型稀疏矩阵而是用Jacobi迭代逐点更新。原因有两个一是迭代法在64乘64网格上收敛很快代码也直观二是避免在MATLAB里维护稀疏矩阵存储。代价是迭代次数和收敛容差需要手工控制。核心代码% 计算散度 RHS div_u zeros(nx, ny); div_u(2:nx-1, 2:ny-1) (ustar(3:nx, 2:ny-1) - ustar(2:nx-1, 2:ny-1)) / dx ... (vstar(2:nx-1, 3:ny) - vstar(2:nx-1, 2:ny-1)) / dy; rhs div_u / dt; % Jacobi 迭代求解压力泊松方程 p zeros(nx, ny); tol 1e-6; maxiter 2000; for iter 1:maxiter p_old p; % 内部节点五点差分 p(2:nx-1, 2:ny-1) 0.25 * ( ... p_old(1:nx-2, 2:ny-1) p_old(3:nx, 2:ny-1) ... p_old(2:nx-1, 1:ny-2) p_old(2:nx-1, 3:ny) - ... dx^2 * rhs(2:nx-1, 2:ny-1) ); % 边界Neumann 零梯度条件 p(1,:) p(2,:); p(nx,:) p(nx-1,:); p(:,1) p(:,2); p(:,ny) p(:,ny-1); if max(abs(p(:) - p_old(:))) tol break; end endrhs的计算是整个投影法的心脏这里散度差分必须和后面速度修正里的压力梯度差分使用完全相同的网格位置。如果散度用(u(i1)-u(i))/dx梯度也必须用对应的交错位置混用会导致修正后仍然有残余散度。Jacobi迭代的好处是天然满足Neumann边界条件边界处直接让边界压力等于相邻内点压力。tol取1e-6在大多数算例够了如果你观察流场有缓慢的周期波动可以调到1e-8。3.3 第三步修正速度让散度严格为零得到压力场后修正步用一个简单的梯度减法完成$$u^{n1} u^* - dt \frac{\partial p}{\partial x}, \quad v^{n1} v^* - dt \frac{\partial p}{\partial y}$$这一步做完之后理论上速度场就满足连续方程了。MATLAB实现很短% 修正步压力梯度修正速度场 u(2:nx, 2:ny-1) ustar(2:nx, 2:ny-1) - dt / dx * (p(2:nx, 2:ny-1) - p(1:nx-1, 2:ny-1)); v(2:nx-1, 2:ny) vstar(2:nx-1, 2:ny) - dt / dy * (p(2:nx-1, 2:ny) - p(2:nx-1, 1:ny-1)); % 主时间循环骨架 for n 1:nt % 1. 边界条件施加固壁无滑移或解析解边界 apply_boundary_conditions(u, v); % 2. 预测步得到 ustar, vstar % 3. 求解压力泊松方程 % 4. 修正速度场 % 5. 记录输出速度范数、残差、动画 end每步循环顺序不能乱先施加边界再预测、投影、修正。如果边界条件在预测步之后才施加边界上会积累误差。非定常的意义在于主循环里的每一步都对应真实的物理时间每一个中间场都是物理上可解释的流动状态。衰减涡算例验证的正是这一点同一时刻的数值解必须无限逼近解析解而不是只求“最终定常结果正确”。4. 边界条件与两个标准算例把代码真正跑起来4.1 顶盖驱动流非定常推进到定常的经典路径顶盖驱动流是所有NS代码的标配验证算例。计算域是单位方腔顶盖沿x方向以恒定速度运动其余三壁固定。雷诺数的定义是Re U0 * L / nu其中U0是顶盖速度L是方腔边长。这份示例里就取U0等于1、L等于1所以Re 1/nu代码里直接通过Re控制粘性系数大小。边界条件的实现是一行行硬编码的% 顶盖驱动流边界 % 下壁、左壁、右壁无滑移 u0, v0 u(:, 1) 0; % 下壁 u(:, ny) 1; % 顶盖水平速度设为 1 v(:, 1) 0; % 下壁 v(:, ny) 0; % 顶盖竖直速度为零 % 左壁与右壁 u(1, :) 0; u(nx1, :) 0; v(1, :) 0; v(nx1, :) 0;跑顶盖驱动流时有个容易误解的点题目叫非定常但顶盖驱动流的经典验证目标是“定常解”。启动瞬间整个流场是一潭死水的顶盖运动后涡逐渐发展经历一个过渡过程最终达到定常状态。所以这份示例既能展示非定常演化又能通过定常终态验证收敛性。判断定常别用“跑够多少步”而是看速度场的变化量% 定常收敛判据连续两步速度场变化小于阈值 diff_u max(abs(u(:) - u_prev(:))); diff_v max(abs(v(:) - v_prev(:))); if max(diff_u, diff_v) 1e-8 break; endRe100时粗网格跑到定常大概需要几千步细网格上万步。如果Re提高到1000流场会呈现出更复杂的二次涡结构收敛也更慢这也是后续练手的好对象。4.2 衰减涡算例用解析解自校验时间精度顶盖驱动流只能验证定常结果对不对验证非定常时间积分的正确性需要另一个算例衰减涡。它有一个解析解在方腔内用正弦函数构造一个涡其强度随时间指数衰减。这个算例的最大价值是能把每个时刻的数值解和精确解对比直接算出误差阶。衰减涡的解析速度场为% 衰减涡初始场正弦涡 % 解析解u sin(pi*x) * cos(pi*y) * exp(-2*pi^2*t/Re) % v -cos(pi*x) * sin(pi*y) * exp(-2*pi^2*t/Re) x_c ((1:nx)-0.5) * dx; y_c ((1:ny)-0.5) * dy; [Xc, Yc] meshgrid(x_c, y_c); % 初始时刻 t0 u sin(pi*Xc) .* cos(pi*Yc); v -cos(pi*Xc) .* sin(pi*Yc);注意初始场给的是单元中心位置而u、v数组在交错网格上分别位于不同位置所以要对x、y坐标做适当的偏移这又是一个索引细节。误差计算在任意时刻做% 在给定时刻对比数值解与解析解 t n * dt; u_exact sin(pi*Xc) .* cos(pi*Yc) .* exp(-2*pi^2*t/Re); err_u sqrt(sum((u(2:nx1, :) - u_exact(:,:)).^2) / numel(u_exact));这里我习惯用相对L2误差来判断而不是看最大误差因为最大误差对个别点太敏感。跑衰减涡时边界条件不能用固壁无滑移而应该把解析解直接赋在边界上这样误差只来自内部离散不会被边界处理污染。4.3 后处理看见流场才算真正跑通跑完计算后至少看三个图速度矢量图、压力云图、流线图。MATLAB里没有特别多技巧核心命令就这几个% 视觉化流场 figure quiver(x_u, y_u, u, v, 2); % 速度矢量2是缩放比例 hold on contourf(x_c, y_c, p, 20); % 压力云图 axis equal title(二维非定常方腔流瞬时流场); % 导出数据便于后续分析 save(cavity_result.mat, u, v, p, t, Re, nx, ny, dx, dy, dt);速度矢量图能直观看到主涡和角涡压力云图能看出高低压区分布流线图适合展示涡心位置。保存为后缀为mat的文件后之后再想画中心线速度剖面、算涡量、做动画都不用重跑一遍计算这点很重要因为在细网格上重跑一次的代价不小。5. 避坑指南发散、不收敛与伪非定常的五个典型现场5.1 前几步就爆NaN十有八九是dt超限现象程序启动后前几个时间步速度场就出现Inf或NaN绘图时报错整个计算直接终止。原因时间步太长超过了上一章说的对流或粘性稳定性边界。尤其是低雷诺数加细网格时粘性约束项Re*dx^2/2往往比CFL约束更严格只按CFL取dt必爆。另一个常见原因是初始速度场里存在不合理的局部高速动态dt计算时被低估。解决把dt计算改成动态的每步都重新扫描全场最大速度并同时套用对流和粘性两个约束。如果仍然爆炸把安全系数从0.25降到0.1试一次。这个问题的特征是“换个dt就从爆变稳”不用怀疑代码逻辑。5.2 压力场棋盘振荡检查你的网格是不是“真交错”现象速度场看起来大体合理但压力云图呈密集的红黑棋盘格交替数值上相邻格点压力正负交替出现。原因压力棋盘振荡是泊松方程零能模态的表现通常因为u、v、p数组尺寸相同也就是名义上“交错”实际没交错或者是散度和梯度算子用的差分方向不一致导致压力修正没有作用到速度的正确位置。解决先确认数组尺寸满足size(p)[nx,ny]、size(u)[nx1,ny]、size(v)[nx,ny1]。再检查修正步里的压力梯度索引u(2:nx,...)对应p(2:nx,...)-p(1:nx-1,...)如果写成了p(3:nx1,...)-p(2:nx,...)压力梯度会整体平移一格相当于用错了离散模板。这两点检查完棋盘振荡基本消失。5.3 顶盖驱动流左上角压力爆表角点奇点处理不当现象顶盖驱动流跑到几百步后左上角和右上角附近的压力异常大速度等值线在角点处剧烈扭曲。原因顶盖水平速度是1左壁水平速度是0这两个边界条件在角点直接冲突形成了数学上的奇点。网格越细角点处的速度梯度和压力越尖锐。解决在角点处对顶盖速度做一个平滑渐变常见做法是把角点速度设为零或者用一段余弦渐变过渡。代码里就是一行的事u(nx1, ny) 0; % 右上角速度置零消除奇点 u(1, ny) 0; % 左上角速度置零加这一行之后全局压力最大值会明显回落。这个改动只影响角点附近一个网格的精度对整体流场影响很小。5.4 衰减涡误差不降反升边界条件施加错了现象网格从32加密到64再到128误差几乎不变甚至在增大完全不符合二阶格式应有的表现。原因衰减涡算例的边界上如果把固壁条件当边界用边界附近的数值解会和解析解完全脱节误差由边界主导加密网格也救不回来。另一个原因是时间步长没同步缩小空间二阶但时间一阶时间误差盖过了空间误差。解决衰减涡的所有边界都直接赋解析解的值用当前位置和当前时刻的精确公式计算边界上的u、v。同时把dt减半和网格加密两件事分开做先固定dt只加密网格测空间精度再固定网格只减dt测时间精度就能定位到底是哪一阶拖了后腿。5.5 看着在动其实是伪非定常收敛判据太乐观现象顶盖驱动流跑了预设步数流场每一时刻都在缓慢变化看起来好像有非定常特征但速度变化量一直降不下去始终在某个水平抖动。原因两种情况。一种是收敛阈值设得太宽比如1e-4就停流场还没真正平衡另一种是压力泊松迭代的tol太松每步的残余散度没有被完全压掉导致速度场被反复注入微小扰动。解决把定常判据改成连续两步速度场最大差小于1e-8并让压力泊松迭代的tol也降到1e-8再看曲线是否平滑收敛。如果残差在某个数值附近震荡不降切割压力泊松迭代次数通常能直接见效。记住一个习惯判断定常要看历史曲线不要当场拍脑袋。6. 进阶从复现走向验证的快速检查法拿到这份MATLAB示例后我建议你先别急着改物理参数而是花半天时间做三件验证都很便宜。第一件是衰减涡的时间收敛阶测试。固定网格在64乘64从dt0.01开始逐步减半记录每一档的L2误差看误差是否大约按dt的第一次方下降。如果是说明时间推进和投影法配合正确如果误差不随dt变化说明时间误差不是主导问题可能在边界或初始场。dt速度L2误差比上一档变化的倍数0.012.1e-3-0.0051.05e-3约1/20.00255.2e-4约1/2第二件是顶盖驱动流的速度剖面验证。在定常状态下取出方腔水平中线上u速度随y的分布和同雷诺数的公开基准数据对比。Re100时中心处回流的负速度峰值大约在-0.2附近如果你的结果差出10%以上优先怀疑网格分辨率或角点平滑的影响。第三件是质量守恒检查。计算任意时刻全场散度的平方和这个数应当在压力泊松迭代的容差量级如果明显偏大说明投影法第三步的梯度索引写错了。性能优化方面整套代码最大的瓶颈在压力泊松方程的Jacobi迭代可以用MATLAB矩阵操作代替双重for循环加速比通常在5倍以上。如果网格超过128乘128把迭代改成预条件共轭梯度法会更划算但要注意边界条件的组装方式不同。这份示例包自带的算例控制在几十秒内跑完没必要为了性能牺牲代码可读性。从那以后我每次拿到任何NS求解代码都会强制自己先跑一遍衰减涡误差测试再谈别的。误差阶不对后面的流场图再漂亮也是废的。希望这份MATLAB二维非定常示例能帮你少走这几步弯路早日把注意力放到真正要做的流动问题上。本文还有配套的精品资源点击获取