简介这份资源面向学习计算流体力学与数值模拟的高校学生、科研人员及工程技术人员提供二维非定常Navier-Stokes方程的MATLAB求解示例帮助读者理解不可压缩流体在二维空间内运动的数值实现思路。压缩包共76个文件约463KB以44个m脚本为核心涵盖网格定义、形函数与高斯点计算、对流项与粘性项矩阵组装、压力修正及边界约束等模块另有29张png结果图与1个html页面用于展示速度场与流线演化便于对照代码理解每一步计算。资源已有1947人学习下载适合作为有限元或有限差分求解N-S方程的入门实践素材。读者可从中获得从方程离散、矩阵组装到时间推进与可视化输出的完整代码框架并借助示例图直观感受非定常流动的演化过程为空气动力学、水力学等方向的数值模拟打下基础。1. 二维非定常 Navier-Stokes为什么“能跑出漩涡”不等于“算对了”很多人第一次在 MATLAB 里跑二维非定常 Navier-Stokes 示例看到流场里卷出几个漂亮的涡就觉得大功告成。但真正决定结果可信度的往往不是那几个涡好不好看而是压力泊松方程有没有解对、对流项离散是否守恒、时间推进是否满足稳定性约束。二维 Naiver-Stoke 非定常 示例的价值恰恰在于它把“不可压”这个约束如何一步步变成可执行的矩阵运算摊开给你看速度场先被对流和扩散推动再被压力梯度强行拉回无散状态。搞懂这条链路你才知道代码里每一行 reshape、每一次 sparse 求解在干什么也才能在结果发散时判断是物理不稳定性还是数值翻车。这份 MATLAB 示例值得逐行拆因为它把教科书上跳过的投影步细节补上了。2. 投影法不可压约束怎么变成一次泊松求解2.1 从动量方程到中间速度不可压 Navier-Stokes 的难点不在动量方程本身而在连续性方程没有时间导数——速度场必须时刻满足散度为零压力没有独立的演化方程它只是一个拉格朗日乘子。常见的实现思路是 Chorin 投影法先忽略压力梯度用当前压力场把速度往前推半步得到一个“中间速度” u*这个中间速度一般不满足无散条件再用压力泊松方程求出压力修正量把 u* 投影回无散空间。写成半离散形式动量方程的时间推进是% u,v: 当前时刻速度含边界 % p: 当前时刻压力 % dt: 时间步长, nu: 运动粘度 % 对流项用迎风/中心差分扩散项用中心差分 u_star u dt * (-convect_u(u,v) nu * laplacian_u(u,v)); v_star v dt * (-convect_v(u,v) nu * laplacian_v(u,v)); % 注意这里没有压力梯度项压力在投影步才出现这段代码里convect_u和laplacian_u分别是对流和扩散算子的离散。关键点在于中间速度 u* 的边界条件必须和最终速度一致否则投影步会把边界误差放大成全局散度。常见的实现思路是把壁面法向速度设为零切向速度按无滑移或自由滑移处理u* 在壁面处直接继承这个值。2.2 压力泊松方程怎么组装对 u* 取散度并要求修正后的速度散度为零就得到压力泊松方程∇²p (1/dt) ∇·u*离散后这是一个稀疏线性系统。二维结构化网格上用五点差分格式组装拉普拉斯矩阵右端项是中间速度的离散散度。MATLAB 里通常用sparse组装再用\或pcg求解。这里有个容易被忽略的细节纯 Neumann 边界条件下拉普拉斯矩阵是奇异的压力只能确定到相差一个常数必须固定一个参考点或者对右端项做零均值化否则求解器会给出警告甚至发散。% 组装压力泊松方程 A*p rhs % Nx,Ny: 网格点数, dx,dy: 网格间距 N Nx * Ny; % 五点差分拉普拉斯内部点 A spalloc(N, N, 5*N); rhs zeros(N,1); for j 2:Ny-1 for i 2:Nx-1 idx (j-1)*Nx i; A(idx, idx) -2/dx^2 - 2/dy^2; A(idx, idx-1) 1/dx^2; A(idx, idx1) 1/dx^2; A(idx, idx-Nx) 1/dy^2; A(idx, idxNx) 1/dy^2; % 右端项中间速度的离散散度 rhs(idx) (u_star(idx1)-u_star(idx-1))/(2*dx) ... (v_star(idxNx)-v_star(idx-Nx))/(2*dy); rhs(idx) rhs(idx) / dt; end end % 固定参考点消除奇异性 A(1,:) 0; A(1,1) 1; rhs(1) 0; p A \ rhs;这段代码里A(1,:)0; A(1,1)1就是把第一个网格点的压力钉死为零这是处理纯 Neumann 问题的标准做法。如果不做这一步A\rhs可能返回一个巨大偏移的压力场速度修正后虽然散度接近零但压力梯度会把流场整体推飞。参数dx、dy必须和速度离散用的网格间距一致否则散度算子和梯度算子不匹配投影后的残余散度会随时间累积。2.3 速度修正与散度检查求出压力后用压力梯度修正中间速度u u* - dt * ∇p这一步的离散必须和压力泊松方程里的散度算子互为伴随否则投影不精确。常见的实现思路是散度用中心差分梯度也用中心差分保证离散意义下的∇·∇p和拉普拉斯矩阵一致。修正完一定要算一次最大散度作为数值健康的指标。% 用压力梯度修正速度 for j 2:Ny-1 for i 2:Nx-1 idx (j-1)*Nx i; u(idx) u_star(idx) - dt * (p(idx1)-p(idx-1))/(2*dx); v(idx) v_star(idx) - dt * (p(idxNx)-p(idx-Nx))/(2*dy); end end % 检查投影后的最大散度 div_max 0; for j 2:Ny-1 for i 2:Nx-1 idx (j-1)*Nx i; div (u(idx1)-u(idx-1))/(2*dx) (v(idxNx)-v(idx-Nx))/(2*dy); div_max max(div_max, abs(div)); end end fprintf(max divergence after projection: %.3e\n, div_max);div_max应该和机器精度同一量级如果它在 1e-3 以上说明压力求解没收敛或者边界处理有误。这个检查比看流线图有用得多流线图只能告诉你“像不像”散度告诉你“对不对”。3. 对流项离散为什么中心差分会在高雷诺数下翻车3.1 中心差分的色散与振荡对流项 u·∇u 是非线性的也是整个求解里最容易出问题的地方。中心差分精度高但在高雷诺数下会产生非物理振荡原因是中心差分对高频模态的耗散为零能量会堆积在网格尺度上。二维非定常示例里如果直接用中心差分跑一个剪切层你会看到速度场出现棋盘格一样的抖动时间步越小抖动越明显——这不是物理不稳定是离散格式的色散误差。常见的实现思路是改用迎风差分或者加人工粘性。一阶迎风稳定但数值耗散大会把涡抹平二阶迎风或者 QUICK 格式折中一些。MATLAB 示例里为了代码简洁往往用一阶迎风代价是有效雷诺数被数值粘性拉低。% 一阶迎风对流项u 方向为例 % 根据 u 的符号选择上游差分 for j 2:Ny-1 for i 2:Nx-1 idx (j-1)*Nx i; if u(idx) 0 du_dx (u(idx) - u(idx-1))/dx; else du_dx (u(idx1) - u(idx))/dx; end % v 方向同理此处省略 convect_u(idx) u(idx)*du_dx v(idx)*du_dy; end end迎风差分的逻辑是信息沿速度方向传播所以取上游点的值。u(idx)0时上游在左侧用u(idx)-u(idx-1)否则上游在右侧。这个分支在向量化时可以用max/min或者逻辑索引批量处理但逐点循环更容易看清机制。注意迎风差分引入的数值粘性量级是|u|*dx/2如果物理粘度nu比这个小很多你算出来的雷诺数其实是被数值粘性主导的涡的衰减会比理论快。3.2 时间推进显式格式的稳定性红线非定常问题的时间推进常用显式格式比如前向欧拉或者低阶 Runge-Kutta。显式格式的稳定性受 CFL 条件约束对流要求dt dx/|u|扩散要求dt dx^2/(4*nu)。二维问题里扩散约束通常更苛刻网格一加密dt要按平方缩小这就是为什么高分辨率非定常模拟的计算量涨得飞快。% 估算稳定时间步长 u_max max(abs(u(:))); v_max max(abs(v(:))); dt_conv min(dx/max(u_max,eps), dy/max(v_max,eps)); dt_diff 1 / (2*nu*(1/dx^2 1/dy^2)); dt 0.5 * min(dt_conv, dt_diff); % 安全系数 0.50.5这个安全系数是血泪经验理论上 CFL1 是边界但非线性项和边界条件会让实际稳定域比线性分析小取一半留余量。如果跑着跑着速度突然指数增长先检查dt是不是超了再看对流项离散有没有在高速区产生负值。3.3 边界条件无滑移壁面怎么落地二维方腔驱动流是这类示例的经典算例三面壁面无滑移顶盖以恒定速度拖动。实现时速度边界直接赋值压力边界用 Neumann 条件法向梯度为零。容易翻车的地方是角点顶盖和侧壁的交点速度既受顶盖拖动又受侧壁约束处理不一致会让角点附近出现压力奇异迭代久了散度检查会报警。常见的实现思路是角点速度取平均或者直接归零压力在角点用单侧差分。4. 避坑与排查非定常求解里最常见的四类翻车4.1 现象压力场出现棋盘格振荡原因压力和速度用了同位网格中心差分的散度和梯度算子不满足 inf-sup 条件压力出现奇偶失联。解决改用交错网格MAC 网格或者在同位网格上加 Rhie-Chow 插值。MATLAB 示例如果为了代码简洁用同位网格压力振荡是预期内的不要以为是求解器坏了。4.2 现象散度检查始终在 1e-2 量级降不下去原因压力泊松方程的边界条件和速度边界不匹配或者拉普拉斯矩阵没有处理奇异性。解决先确认纯 Neumann 问题固定了参考点再检查壁面处法向速度是否严格为零。如果壁面法向速度有 O(dt) 的残余散度就会卡在这个量级。4.3 现象涡的强度随时间单调衰减比物理预期快原因一阶迎风差分的数值粘性太大有效粘度远高于设定值。解决加密网格数值粘性随 dx 线性下降或者换二阶迎风/中心差分加小量人工粘性。判断方法是算一下|u|*dx/2和nu的比值如果远大于 1说明数值粘性主导。4.4 现象跑几百步后速度突然爆炸原因显式时间推进超了稳定域或者对流项在局部高速区产生了负值触发了非线性不稳定。解决把dt减半重跑观察爆炸步数是否推后如果是就是稳定性问题。同时检查对流项离散在速度过零处有没有符号判断错误。5. 从示例到可信结果离散一致性比代码长度更重要这份二维非定常 Navier-Stokes 示例最值得深挖的地方不是它用了多少行代码而是它把投影法的每一步都暴露出来让你能逐项验证离散一致性。散度算子和梯度算子是否互为伴随、压力泊松方程的边界是否和速度边界匹配、对流项离散的数值粘性是否淹没了物理粘度——这三个问题决定了你跑出来的涡是物理的还是数值的。常见的实现思路是在示例基础上加一个收敛性检查固定物理参数把网格加密一倍、时间步减半看同一时刻的涡量场是否收敛如果两次结果差异很大说明离散误差还没进入渐近区这时候谈“结果对不对”没有意义。我一般会这么理解这份资源的价值它是一块可以拆开看的底盘而不是一辆能直接上路的车。你可以把迎风差分换成中心差分加滤波把显式推进换成半隐式把方腔换成圆柱绕流每换一处都回头看一眼散度检查和稳定性约束。从那以后我每次改完离散格式都强制先跑一遍散度检查和网格收敛再去看流场图。希望帮到你。本文还有配套的精品资源点击获取
