简介这是一套面向流体力学数值计算学习者与Matlab开发者的二维稳态Navier-Stokes方程有限元求解程序重点解决复杂几何与非均匀边界条件下速度—压力场的离散求解问题适合具备有限元基础和Matlab编程经验的读者研读。压缩包为zip格式共45个文件全部为.m源文件体积约19KB整体轻量便于逐文件阅读和二次修改。内容围绕几何网格生成、基函数与高斯积分、弱形式变分、矩阵与右端项组装、Dirichlet/Robin/应力边界处理、Newton初值构造、压力定点约束以及线性系统求解等环节展开并配有解析解与误差计算脚本可用于验证收敛性和分析流场结果。目前已有312人学习。读者可借此梳理二维Navier-Stokes方程的有限元实现流程理解边界条件与压力规范化的处理思路并在此基础上迁移到热传导、扩散等同类偏微分方程问题。1. 二维稳态 Navier-Stokes 有限元程序先判断它能不能替你省掉三天调试把二维稳态 Navier-Stokes 方程从公式搬到 MATLAB 里最容易翻车的往往不是弱形式而是出口边界和压力零空间。这个有限元求解程序面向不可压缩稳态流动核心是 Taylor-Hood 单元、牛顿迭代、边界条件装配和后处理适合做微通道、管道流、方腔驱动这类二维算例。第一次搭流体有限元的人可以照着它跑通一条完整链路已经写过的人重点看单元选型、续算参数和排查表能少走不少重复调试的弯路。2. 弱形式到 Taylor-Hood二维 N-S 方程离散时先锁死这三个量二维不可压缩稳态 Navier-Stokes 方程写成[ -\mu \Delta \mathbf{u}\rho(\mathbf{u}\cdot\nabla)\mathbf{u}\nabla p\mathbf{f}, \qquad \nabla\cdot\mathbf{u}0 ]其中 (\mathbf{u}(u_x,u_y)) 是速度(p) 是压力(\rho) 是密度(\mu) 是动力粘度。有限元程序要做的不是把这条方程原样塞进矩阵而是先做弱形式、选单元对、处理压力零空间再把非线性项送进牛顿迭代。下面三个量如果一开始没锁死后面调残差会非常痛苦速度-压力单元对、压力基准点、自由度编号。2.1 弱形式与压力零空间对动量方程乘速度测试函数 (\mathbf{v})对连续性方程乘压力测试函数 (q)分部积分后得到弱形式求 ((\mathbf{u},p))使得对任意 ((\mathbf{v},q)) 有[ \mu\int_\Omega \nabla\mathbf{u}:\nabla\mathbf{v},d\Omega \rho\int_\Omega (\mathbf{u}\cdot\nabla)\mathbf{u}\cdot\mathbf{v},d\Omega -\int_\Omega p,\nabla\cdot\mathbf{v},d\Omega \int_\Omega \mathbf{f}\cdot\mathbf{v},d\Omega ]以及[ -\int_\Omega q,\nabla\cdot\mathbf{u},d\Omega0 ]入口和壁面用 Dirichlet 条件直接约束速度出口常见做法是只给压力基准不强制速度让出口法向应力自然释放。压力只以梯度形式出现所以 (p) 和 (pC) 对应同一个速度场离散后切线矩阵会天然缺一个约束。数值上必须固定一个压力自由度或者显式做零均值约束否则线性系统会多一个零特征值。我一般先把物性参数和收敛参数收进一个结构体避免后面改雷诺数时到处翻脚本par.rho 1.0; % 密度长度、速度、粘度先统一量纲 par.mu 0.02; % 动力粘度单位与rho、U、H保持一致 par.U 1.0; % 特征速度后面用来算Re和入口剖面 par.H 1.0; % 特征高度 par.Re par.rho * par.U * par.H / par.mu; par.tol 1e-8; % 相对残差阈值按初始残差归一 par.maxit 30; % 牛顿步上限 par.damp0 1.0; % 首步阻尼发散时改到0.3~0.7par.Re不是独立参数改mu或U都会改变它所以入口流量、物性和网格特征长度要一起核对。par.tol建议用相对残差绝对残差会随网格量和方程缩放变化拿绝对残差判断收敛很容易误判。par.damp0不是物理量只是牛顿步的初始缩放低雷诺数下取 1 通常没问题高雷诺数下先降到 0.5 更稳妥。提示如果压力场整体漂移、速度场却看起来正常先检查压力基准点不要急着怀疑对流项。2.2 Taylor-Hood 单元与自由度编号二维三角形网格上速度用二次单元 P2压力用线性单元 P1这就是常说的 Taylor-Hood 单元对。每个三角形速度有 6 个节点压力有 3 个节点速度节点包括三个顶点和三条边的中点。P2-P1 满足离散的 LBB 条件压力不会出现棋盘振荡。等阶 P1-P1 虽然实现简单但压力容易抖动除非额外加压力稳定项否则不建议拿它做第一个版本。自由度编号要一次定清楚。常见做法是先把 P1 顶点编号再给每条边追加一个中点编号速度节点总数就是nP1 nEdge。每个速度节点有两个分量压力节点仍用原来的 P1 顶点。下面这段是编号骨架% mesh.nodes: P1顶点坐标mesh.tri: P1三角形连接 % mesh.edges: 每行两个端点边中点按全局唯一边表生成 nP1 size(mesh.nodes, 1); nEdge size(mesh.edges, 1); mesh.nVel nP1 nEdge; % P2速度节点数 mesh.nPress nP1; % P1压力节点数 % 第1列ux第2列uy不要把所有ux排完再排uy后处理会很难受 mesh.uDof reshape(1:2*mesh.nVel, mesh.nVel, 2); mesh.pDof 2*mesh.nVel (1:mesh.nPress);mesh.uDof(i,1)是第 i 个速度节点的 (u_x) 自由度mesh.uDof(i,2)是 (u_y)。压力自由度整体偏移到速度后面后面组装 15×15 的单元矩阵时前 12 个是速度后 3 个是压力索引不容易乱。边中点必须全局唯一如果每个单元各建一个中点相邻单元的速度自由度不共享连续性方程会漏质量残差看着能降流量却对不上。压力基准点通常选出口中心或出口下角固定p0。只固定一个点就够不要给整条出口压力边界都赋零那会把出口流场压死回流区容易失真。2.3 单元矩阵组装与稀疏矩阵单元组装分两步先算局部残差和局部切线矩阵再按全局自由度索引累加到稀疏矩阵。二维 P2-P1 三角形单元每个单元有 (6\times2315) 个局部自由度所以局部矩阵是 15×15。下面这段是全局组装骨架function [R, J] assemble_ns(mesh, u, p, par) % 返回全局残差R和切线矩阵J尚未施加Dirichlet边界 nU 2 * mesh.nVel; nP mesh.nPress; R zeros(nU nP, 1); rows []; cols []; vals []; for e 1:size(mesh.tri, 1) [Re, Je, idx] local_ns_element(mesh, e, u, p, par); R(idx) R(idx) Re; [rr, cc, vv] find(Je); rows [rows; idx(rr)]; cols [cols; idx(cc)]; vals [vals; vv]; end J sparse(rows, cols, vals, nU nP, nU nP); endidx的长度必须是 15顺序建议固定为[u1x u1y ... u6x u6y p1 p2 p3]。local_ns_element里用高斯积分算形函数和梯度粘性项用 (\nabla N_i\cdot\nabla N_j)压力项用 (N_i \nabla\cdot \mathbf{v})连续性项用 (q_j \nabla\cdot \mathbf{u})。对流项的非线性在于当前速度 (\mathbf{u}_h) 也参与积分牛顿切线里会出现两部分贡献测试函数方向的对流导数和试探函数方向的对流导数漏掉任何一项牛顿迭代都可能只降几步就卡住。sparse组装时不要用三重循环直接给J(i,j)J(i,j)...赋值二维网格稍微加密就会慢得明显。先把rows/cols/vals攒起来最后一次构造稀疏矩阵内存和速度都更好控制。局部矩阵Je里如果有零元素find会自动跳过不影响结果。3. 牛顿迭代与续算把非线性残差从 1e-2 压到 1e-8 的具体写法稳态 N-S 的离散方程是一组非线性代数方程牛顿法每步解一个线性系统 (J\delta-R)再更新 (u\leftarrow u\alpha\delta)。真正难的不是写这个公式而是保证残差和切线矩阵同源、初值别太离谱、步长别一上来就拉满。高雷诺数算例如果直接零初值开牛顿第一步残差变大是常见现象不要把它当成程序写错了。3.1 牛顿残差和切线矩阵要同源残差和切线矩阵必须来自同一套形函数、同一组高斯点、同一个边界条件处理函数。常见错误是残差里用了当前速度切线矩阵里却用了上一次迭代的速度或者残差加了体积力切线矩阵漏了对应的线性项。这样牛顿法可能仍然会降但收敛阶会从二次掉到线性最后卡在 (10^{-4}) 附近。function [u, p, hist] newton_ns(mesh, par, u, p) if nargin 3 [u, p] stokes_initial(mesh, par); % 默认先解Stokes end for it 1:par.maxit [R, J] assemble_ns(mesh, u, p, par); [R, J] apply_bc(R, J, mesh, par); % 边界行必须在R和J上同步 if it 1 par.res0 max(1, norm(R, inf)); end res norm(R, inf) / par.res0; hist(it) res; if res par.tol break; end du -(J \ R); alpha line_search(u, p, du, R, mesh, par); u u alpha * du(1:2*mesh.nVel); p p alpha * du(2*mesh.nVel1:end); end endR是离散残差J是对应切线矩阵du是牛顿增量alpha是线搜索得到的步长。par.res0用第一步残差做归一避免网格加密后绝对残差量级变化导致阈值失效。apply_bc一定要同时改R和J只在切线矩阵里置单位行、残差里不改解出来的增量不会满足 Dirichlet 条件边界节点会慢慢漂掉。残差历史hist最好画出来看。正常牛顿应该在前几步快速下降接近收敛时下降更快如果残差在前两步下降、第三步突然抬头通常是线搜索步长太大或者对流项切线不一致。3.2 用 Stokes 解做初值再按 Re 续算零初值在低雷诺数下能用到了 Re100 以上就容易把牛顿法带偏。更常见的做法是先解 Stokes 方程把粘性项和压力项的解作为初始场再按雷诺数逐步加对流项。每一步用上一步的解热启动牛顿步数通常能压到 5 步以内。ReList [0 5 20 50 100 200]; [u, p] stokes_initial(mesh, par); % Re0对应的线性解 for k 1:numel(ReList) par.Re ReList(k); par.mu par.rho * par.U * par.H / max(par.Re, eps); [u, p, hist] newton_ns(mesh, par, u, p); fprintf(Re%g, res%.3e, it%d\n, par.Re, hist(end), numel(hist)); endReList不要一开始就写得太密先看每一档的牛顿步数。如果某一档从 5 步突然涨到 20 步说明步长偏大下一档可以把增量减半。max(par.Re, eps)是为了避免 Re0 时除以零Re0 那一步实际上就是 Stokes 解不需要再算粘性系数。热启动时u和p的尺寸必须和当前网格一致换网格后要重新插值不能直接拿旧网格的解往新自由度上套。3.3 线搜索与收敛判据线搜索的目标是让残差范数下降不是让速度增量最小。下面这段是简化版回溯线搜索先试完整牛顿步残差没有足够下降就减半最多试 8 次。1e-4是下降判据里的一个小量实际用的时候不用太纠结关键是别让残差在迭代中反复抬头。function alpha line_search(u, p, du, R, mesh, par) alpha par.damp0; for k 1:8 uT u alpha * du(1:2*mesh.nVel); pT p alpha * du(2*mesh.nVel1:end); RT residual_only(mesh, uT, pT, par); if norm(RT, inf) (1 - 1e-4 * alpha) * norm(R, inf) return; end alpha 0.5 * alpha; end alpha min(alpha, 0.1); % 兜底下一步仍失败就缩小Re步长 endpar.damp0是初始步长低雷诺数取 1高雷诺数取 0.5 或 0.3。residual_only只算残差、不算切线矩阵线搜索里反复调用时开销小很多。如果线搜索连续失败不要一直加迭代次数先缩小雷诺数增量或者把网格长宽比降下来后者对切线矩阵条件数的影响经常被低估。收敛判据建议同时看残差和增量相对残差小于1e-8且速度增量的无穷范数小于1e-10再停。只看残差有时会在压力基准附近出现假收敛尤其是出口附近网格很稀的时候。4. 边界条件与网格入口、壁面、出口的约束顺序和参数表二维稳态流动的边界条件比三维好调但坑并不少。入口给什么剖面、壁面约束几个分量、出口要不要固定压力、压力基准点放哪里这几个决定会直接影响回流区和质量守恒。网格也不是越密越好出口距离不够、壁面首层太厚、单元长宽比太大都会让牛顿法在最后几步反复。4.1 入口给抛物型速度别只给平均速度槽道充分发展流的入口速度是抛物型分布壁面为零中心最大。若入口给成均匀速度入口附近会立刻出现一个很大的横向速度梯度高雷诺数本文还有配套的精品资源点击获取
