简介这是一套使用Matlab实现二维稳态Navier-Stokes方程有限元求解的完整源代码面向流体计算与有限元数值方法的学习者和科研人员针对复杂几何与非均匀边界条件问题可帮助理解不可压缩流动从偏微分方程到离散代数系统的完整求解流程。压缩包共45个m文件整体仅19KB全部为Matlab脚本和函数文件涵盖主程序与子函数未附带数据或文档覆盖均匀网格生成、拉格朗日基函数与高斯积分、Dirichlet/Robin/应力边界条件处理、矩阵组装、牛顿迭代求解及误差计算等关键环节。目前已有311人学习下载。源码模块划分清晰既适合逐行研读以掌握有限元方法在Navier-Stokes方程中的落地实现也便于替换源项、边界条件或几何域后推广至热传导、扩散等其他偏微分方程是一份紧凑实用的流体计算参考程序掌握这一过程也有助于后续研究扩展。 前阵子我把二维稳态Navier–Stokes方程这里先提一句项目名里的“Stoks”是笔误全名应该是Navier–Stokes的有限元求解程序从头到尾完整实现了一遍。这个东西解决的核心问题很直观在给定计算域、边界条件和雷诺数的情况下用有限元方法求出二维流场的速度分量u(x,y)、v(x,y)和压力p(x,y)。网上写NS方程理论的帖子不少但能把“理论方程”一路推到“可运行代码”的完整记录并不多。这篇主要面向正在学有限元或计算流体力学的同学也适合想自己动手写二维流场求解器的工程师我会把方案选型、离散细节、代码实现和调试心得一锅端出来尽量讲清楚每个选择背后的原因。1. 为什么用有限元解NS方程——先从方案选型说起1.1 几种主流离散方法对照拿到二维稳态NS方程这个问题第一个要拍板的事情不是写代码而是选数值离散方法。不少初学者上来就问为什么不直接写个有限差分为什么不学OpenFOAM用有限体积这里面的取舍直接决定了程序的数据结构、边界处理方式和可扩展性。方法典型应用优势劣势有限差分 FDM结构化网格上的简单流动编程简单、概念直观复杂几何边界难处理守恒性依赖网格正交性有限体积 FVM工程CFD主流Fluent、OpenFOAM天然保证通量守恒适合复杂几何高阶精度构造繁琐数学理论分析不如FEM清晰有限元 FEM固体力学、流体力学、多物理场耦合复杂几何适应性强可灵活选择插值阶次理论误差估计完备编程复杂度较高不可压缩问题需要满足LBB条件我最终选有限元核心原因是想把NS方程的数学结构吃透。有限元方法能把连续方程从“微积分问题”干净地转化为“代数问题”而且对非规则边界特别友好。二维问题规模不大就算用直接求解器也撑得住所以不用担心线性系统太大导致内存爆炸——这是我敢用FEM做NS求解的重要原因。1.2 把连续方程“翻译”成代数系统——弱形式这件事很多同学卡在有限元的第一道坎就是弱形式。别怕弱形式的本质就是把偏微分方程两边乘上一个“任意测试函数”再在计算域Ω上积分最后利用分部积分把导数“甩”一部分给测试函数。这么做的直接收益有两个一是降低了对解函数光滑度的要求二是边界条件可以自然地进到变分公式里。对于二维稳态不可压缩N-S方程无量纲形式写作连续方程∇·u 0动量方程(u·∇)u −∇p (1/Re)∇²u其中u(u,v)是速度场p是压力Re是雷诺数。把这个方程组分别乘以速度测试函数v和压力测试函数q然后分部积分得到弱形式。这一步做完一个偏微分方程组就变成了一个积分等式接下来才能把计算域切成有限单元在单元上用基函数做插值。2. 核心环节一二维稳态NS方程的离散与线性化2.1 不要踩等阶插值的坑为什么必须用Taylor-Hood元在有限元里选择速度场和压力场的插值空间时有个非常反直觉的规定速度和压力的插值阶次不能随便一样。如果用一阶线性元同时插值速度与压力解出来的压力场常常会出现棋盘状振荡一团黑一团白完全不是真实物理。这个问题的深层原因是不可压缩约束带来的“inf-sup条件”也叫LBB条件。简单理解速度场和压力场的工作空间不匹配会造成压力模式“锁不住”产生伪振荡。可靠的做法是采用Taylor-Hood元即速度场用二阶P2单元压力场用一阶P1连续单元。二维三角形单元下速度节点共有6个3个顶点3个边中点压力节点只有3个3个顶点这套组合能满足LBB条件流体计算领域经典得不能再经典。实际编码时要注意不要单独存储压力场的“单元形函数数组”因为P2速度单元和P1压力单元使用同一套几何网格自由度编号是分开的组装矩阵时要维护好“速度自由度列表”和“压力自由度列表”两份索引。我第一次写的时候就因为索引错位结果压力场呈现出一种奇怪的锯齿形态排查了半天才发现是自由度映射出了问题。2.2 非线性迭代Picard还是NewtonNS方程里的对流项(u·∇)u是非线性的离散后得到的是一组非线性代数方程必须用迭代法求解。常见的两种选择是Picard迭代和Newton迭代。Picard迭代的思路是把对流项中的速度场用上一轮已知迭代值然后求解一个线性化的动量方程相当于在“冻结”系数。它的收敛是线性的速度偏慢但胜在稳健从零初值出发也不容易发散。Newton迭代则是对整个残量向量做一阶泰勒展开求解一个带Jacobian矩阵的线性系统。它收敛快理想情况下是二阶收敛但Jacobian矩阵的推导和组装比Picard复杂一些。还有一个麻烦是它对初始猜测极其敏感初值给得太差会直接发散。我的实际做法是分段组合前5到8次迭代先用Picard把流场带出个大概后续再切换到Newton迭代去冲刺。这样兼顾了稳健性和收敛速度。在代码层面我用一个迭代控制参数控制切换时机比如当Picard的相对残差降到1e-2时切Newton。程序里预留两个函数assemble_picard_system() 和 assemble_newton_system()共用一套网格和自由度映射切换成本很低。雷诺数一旦超过1000Newton迭代基本离不开阻尼。我常用的阻尼方式是限制每个迭代步的速度增量幅度比如把增量缩放因子设为0.5到0.8实测对稳定收敛帮助极大。2.3 边界条件与载荷施加边界条件的正确处理是NS方程有限元程序能不能算得动、算得准的关键。二维稳态问题常见的三类边界条件入口边界给定速度分布Dirichlet条件比如u(U0,0)即抛物入流或均匀入流。壁面边界无滑移条件uv0。出口边界如果给出准确的速度分布可以继续用Dirichlet更自然的方式是所谓“零法向应力”或自然边界条件对应弱形式中边界积分自动消失的情形此时不需要在边界上强制任何速度分量。有一个容易忽略的点压力场的参考点。不可压缩NS方程中的压力是定义在一个常数差范围内的如果不固定压力基准线性系统是奇异的直接求解器会报错。常见做法是固定计算域内某个节点的压力为零或者在压力自由度上施加一个平均值约束。我在代码里用的是“固定一个内部节点压力为0”的方案简单粗暴稳定。3. 实操过程从网格到求解器3.1 网格生成别在第一步偷懒程序跑出来效果不理想很多时候不是求解器的问题而是网格不合格。二维方腔流这类规则计算域可以用Gmsh生成三角形网格也可以自己写脚本生成四边形网格再剖分三角形。一个通用经验角点附近和速度梯度大的区域必须局部加密。我的做法是用Gmsh的geo脚本控制网格尺寸。这里给个可直接套用的方案计算域单位正方形[0,1]×[0,1]最大网格尺寸0.02对应每方向约50个单元中等规模近壁面加密设置 boundary layer第一层厚度0.002层数10增长因子1.2物理模型二维、稳态、不可压缩流体网格生成后导出时选择med或msh格式。msh文件的节点坐标和单元连接信息在Gmsh里用gmsh -tol的基本窗口操作就能导出为MSH 2.2版本并解析节点列表和单元列表。需要注意检查网格质量最小三角形质量指数最好大于0.3否则单元大钝角会导致刚度矩阵条件数变差直接求解器虽然不会崩但解出来会有意想不到的振荡。3.2 单元组装与全局系统单元层面的计算是有限元程序的核心循环。对着每个三角形单元需要完成取出该单元6个速度节点和3个压力节点的坐标在参考单元上建立P2和P1的形函数及导数用Gauss积分计算单元刚度矩阵将对流项、扩散项、压力梯度和连续方程阻力项四块按自由度顺序组装成局部矩阵把局部矩阵通过编号映射散到全局稀疏矩阵中。全局系统的结构是块矩阵[ A B^T ; B 0 ] [ u ; p ] [ f ; 0 ]其中A包含扩散项和对流项的贡献B是速度散度压力梯度耦合项右下角是零块。这个块结构对应的是鞍点问题也解释了为什么最好用能处理非正定系统的直接求解器。数值积分时三角形单元我用7点Gauss积分足以精确积分P2元对应的四阶多项式。不要为了图快降到“1点积分”或“3点积分”那会直接损精度高Re下还会出现不可控的单元寄生振荡。3.3 线性求解器选择二维问题我强烈建议用直接求解器系统自由度数通常在1万到10万之间直接求解完全扛得住。用UMFPACK或者PARDISO这种稀疏LU分解求解器一次分解之后解多右端向量非常快这对非线性迭代特别友好——每一轮只需回代一次。我的推荐组合小规模自由度2万UMFPACK / Suitesparse自带的lu。中大规模自由度5万以上且迭代轮次多PARDISO如果你用的是MKL。如果想扩展到大三维问题再考虑迭代法GMRESILU预条件但二维稳态问题真没必要自虐。组装时务必用CSR或者COO稀疏格式。我一开始图方便用了稠密矩阵存全部自由度1万自由度时内存直接爆掉。换成稀疏存储后内存占用从GB级降到几十MB。3.4 后处理与基准验证程序写完后第一件事不是看漂亮云图而是跑基准题。二维方腔流lid-driven cavity flow是NS方程有限元代码绕不开的“试金石”上边界以恒定速度向右运动其余三个边界固定Re100Re400Re1000是三个经典验证工况。我把程序算出的几何中心涡心位置和速度分量拿出来和Ghia等在1982年发表的基准值对照Re400时涡心约在(0.5547, 0.6055)我的程序在网格规模80×80时误差小于1%。这里提醒一句初学时跑方腔流特别容易被“压差云图看着很漂亮”骗到正确做法是导出中心线上的速度分量做定量对比否则某些振荡问题会被视觉掩盖。后处理我用的是Paraview的VTK格式输出。程序里将每个单元上的速度分量和压力写成一个legacy VTK文件Paraview里直接绘制速度流线、压力云图和涡量场。流线用Tecplot或Paraview的stream tracer都很方便。4. 容易踩的坑与调试心得4.1 压力场出现棋盘振荡现象压力云图呈现红白交替格状分布严重时噪声淹没物理信号。原因最常见的是速度压力插值空间不满足LBB条件。如果你用的是P1/P1速度和压力都线性出现伪振荡是必然。另一个原因是压力没有设置参考值整体未约束导致求解失败或解异常。处理首先升级到P2/P1插值。其次检查压力自由度是否至少固定了一个节点。最后如果压力仍然抖观察网格是否有极小角度的长条单元这类单元会放大压力梯度插值误差。4.2 迭代发散先怀疑初值现象Newton迭代到第三步残差突然飙升到1e10。原因NS方程对流项强非线性零初值出发时速度场毫无物理基础Jacobian矩阵在迭代中途把解推出稳定域。处理两条路一条是先Picard迭代另一条是给速度场一个合理初值比如入口速度的线性缩小版本。配合阻尼更新Newton迭代就稳定很多。实测如果Re1000阻尼系数0.6Picard先走10步Newton再走5步一般就能让残差掉到1e-8量级。4.3 高雷诺数下数值振荡是对流项在捣乱现象Re超过500以后速度场在某些区域出现类似“波纹”的振荡网格加密后反而更明显。原因对流占优对流项系数远大于扩散项时标准Galerkin有限元缺少迎风效应离散后会在网格尺度上产生虚假数值振荡。处理这类问题有两条技术路线。一条是加密网格直到网格Peclet数Pe ρU h/ μ降到2以下但Re高时网格需求爆炸。另一条是引入稳定化项例如SUPG或GLS稳定方法在对流方向添加人工耗散在不显著牺牲精度的情况下压制数值振荡。二维稳态程序里SUPG实现相对简单属于性价比极高的投资。4.4 收敛判据不能只看残差绝对值现象程序报告“收敛”但速度场完全不对。原因残差绝对值取决于方程无量纲化方式有时候物理尺度很小残差降几个量级也未必代表解真正收敛。处理我用相对残差加解变化量双判据。具体是每一轮迭代计算残差向量的L2范数并除以初始残差的L2范数当该比值小于1e-8同时速度解两次迭代的最大变化量小于1e-8才判定收敛。另外务必监控连续方程残差。因为压力只通过连续方程间接出现在求解系统中连续方程残差不过关压力场通常已经出了问题。5. 一个小建议把程序按模块拆开最后分享一点维护层面的心得。二维稳态NS方程的有限元程序虽然逻辑复杂但功能模块非常清晰我在实际编码中严格分成五个部分网格读取、自由度和索引管理、单元刚度矩阵组装、非线性迭代与线性求解、后处理输出。每个模块之间用接口传递数据避免母函数里堆砌一大堆全局数组。一开始可能是为了跑通一个案例但你会发现当你把程序的模块边界划清楚后换边界条件、换稳定化方法、甚至换成三维扩展都只是替换其中一个模块的事。我自己后来在这个框架上加了自然对流和温度场耦合只多写了半天的代码量收益非常可观。从零把NS方程有限元求解器搭出来的过程确实磨人但它会逼着你把从弱形式到稀疏求解的每个环节都搞清楚这种理解深度是直接调用成熟软件很难获得的。本文还有配套的精品资源点击获取
