简介面向流体力学、有限元数值计算与MATLAB编程学习者一款基于MATLAB实现二维稳态Navier-Stokes方程有限元求解的程序包覆盖从网格生成、基函数与变分形式构造到线性系统组装、边界条件处理及后处理的完整数值流程。程序全程采用自编MATLAB代码而非直接调用PDE工具箱便于初学者理解每个步骤的底层实现也为针对特定流动问题修改算法留出空间。压缩包共45个文件全部为.m源文件体积仅19KB结构紧凑其中包含均匀网格生成、局部基函数、矩阵与向量组装、Dirichlet/应力/Robin边界条件处理、Newton迭代求解、H_s误差计算等模块可直接运行配套测试脚本观察速度与压力分布。已有312人学习浏览。借助源码中清晰的分步实现与测试算例读者不仅能掌握速度-压力耦合问题中有限元离散与求解的要点还可将这套框架迁移至热传导、扩散方程等其他偏微分方程场景是开展教学实验或科研预研的轻量参考。1. 这个程序到底在算什么二维稳态N-S方程有限元求解的本质你手里这个标题——二维稳态Navier-Stoks方程有限元求解程序其实少了一个e标准写法是Navier-Stokes。但这不妨碍它指向一个非常具体的工程诉求在二维计算域里给定不可压缩流体的边界速度用有限元法算出稳态的速度场和压力场。它不追踪时间不算湍流只解决边界条件给定了稳定流场长什么样这一件事。适合用来做格式验证、写论文里的CFD算例或者给非稳态求解器铺路。它的核心代码并不长真正难的是理解弱形式、选对速度压力空间、调准迭代参数——这三个地方任何一个翻车结果都会很离谱。2. 先推弱形式为什么有限元能求解N-S方程2.1 强形式到弱形式把光滑性要求降下来二维稳态不可压缩N-S方程的强形式是三条方程质量守恒∂u/∂x ∂v/∂y 0 动量守恒x u∂u/∂x v∂u/∂y -∂p/∂x ν(∂²u/∂x² ∂²u/∂y²) 动量守恒y u∂v/∂x v∂v/∂y -∂p/∂y ν(∂²v/∂x² ∂²v/∂y²)这里的u和v是速度分量p是压力ν是运动黏度。强形式要求u和v至少有两阶连续导数真实流场往往达不到这个光滑度尤其在角点、边界层处。有限元的思路是乘上测试函数在区域内积分再对扩散项做一次分部积分把二阶导数转嫁给测试函数。这样弱形式只要求速度的一阶导数和测试函数的一阶导数存在。最终弱形式可以写成找速度u和压力p使得对所有测试函数v和q都有ν∫(∇u:∇v)dx ∫((u·∇)u)·v dx ∫p(∇·v)dx ∫q(∇·u)dx 0其中第一项代表黏性扩散第二项是对流项第三项是压力梯度贡献第四项是连续性约束。这个式子就是后续编程时写进变分形式的原材料。注意这里压力项取了正号这是通过把动量方程中压力梯度用分部积分得到的符号约定不同不改变物理结果但会改变鞍点矩阵的对称性。2.2 混合有限元与LBB条件为什么P1/P1会翻车观察弱形式会发现速度u和压力p同时出现在一个鞍点系统里不像泊松方程那样只有一个未知量。这意味着速度空间和压力空间不能随意选。经典的LBB条件也叫inf-sup条件要求离散后的速度压力和连续情形一致否则压力场会出现棋盘状振荡。我见过许多从标量有限元转过来的人第一反应是速度和压力都用P1线性元结果算出的压力云图呈红黑相间的棋盘看起来诡异又没法用。原因就是P1/P1不满足LBB条件。教科书级别的稳定做法是Taylor-Hood元速度用P2三角形节点和边中点共6节点压力用P1仅三角形顶点。在二维三角剖分上这个选择既稳定又精度足够。所以代码里速度空间写成VectorFunctionSpace(mesh, P, 2)压力空间写成FunctionSpace(mesh, P, 1)是两个不同的度。如果你非要P1/P1也可以但必须加稳定化项比如PSPG或罚项。新手不要一开始碰稳定化那属于调参玄学范畴老老实实P2/P1绝大多数算例都能正常跑。2.3 最小可运行程序FEniCS实现方腔驱动流下面用FEniCS实现一个经典算例——方腔驱动流。计算域是[0,1]×[0,1]的正方形顶盖以水平速度1向右滑动其余三边无滑移壁面。这个算例是N-S求解器的hello world几何简单但流动结构丰富适合验证程序正确性。from dolfin import * # 1. 网格与物理参数 mesh UnitSquareMesh(32, 32) # 32 x 32 三角形网格 Re 100.0 nu Constant(1.0 / Re) # 无量纲黏度 # 2. Taylor-Hood 混合空间速度 P2压力 P1 V VectorFunctionSpace(mesh, P, 2) Q FunctionSpace(mesh, P, 1) W V * Q # 3. 未知函数和测试函数 w Function(W) # 速度压力组合变量 u, p split(w) v, q TestFunctions(W) # 4. 边界条件无滑移 顶盖拖动 压力锁定 noslip DirichletBC(W.sub(0), Constant((0, 0)), on_boundary !near(x[1], 1.0)) lid DirichletBC(W.sub(0), Constant((1, 0)), near(x[1], 1.0)) # 固定一个压力点消掉常数自由度 pin DirichletBC(W.sub(1), Constant(0), x[0] DOLFIN_EPS x[1] DOLFIN_EPS, methodpointwise) bcs [noslip, lid, pin] # 5. 非线性变分形式 F (nu * inner(grad(u), grad(v)) * dx inner(dot(grad(u), u), v) * dx p * div(v) * dx q * div(u) * dx) # 6. 牛顿求解 J derivative(F, w) problem NonlinearVariationalProblem(F, w, bcs, J) solver NonlinearVariationalSolver(problem) solver.parameters[newton_solver][maximum_iterations] 20 solver.parameters[newton_solver][absolute_tolerance] 1e-8 solver.solve() # 7. 查看中心点压力值验证压力唯一性 print(p at center:, p((0.5, 0.5)))代码逻辑说明第1步用UnitSquareMesh生成结构化三角网格Re100时黏度取1/Re。第2步构造Taylor-Hood空间。第4步的边界条件区分顶盖和其余边界near(x[1], 1.0)是浮点安全写法避免x[1] 1.0的判断失误。压力固定点选在左下角附近DOLFIN_EPS是FEniCS内置的浮点零阈值。参数说明maximum_iterations20是牛顿迭代次数的上限一般Re100在10次内收敛absolute_tolerance1e-8是残差绝对容忍度如果只关心定性流场设1e-6也可以但后处理速度曲线时会看到毛刺。这段代码在Re100下从零初值通常能直接收敛如果发散不要硬扛先看第4章的续延策略。3. 网格与边界条件误差藏在三角剖分的细节里3.1 三角剖分怎么影响精度误差积分的按阶收敛有限元程序本质是在每个三角形上做积分和组装。UnitSquareMesh(32, 32)把每个正方形切分成两个直角三角形生成的一致剖分。这种结构化网格对方腔算例足够但对复杂几何就得用非结构三角网格。网格质量最关键的指标是三角形最小角长宽比过大的薄三角形会让Jacobian矩阵病态直接拉低精度。在二维三角剖分中有限元的误差理论有一个常用的表述如果精确解足够光滑P2速度场的H^1范数误差是O(h²)P1压力场的L^2误差也是O(h²)但压力H^1误差只有O(h)。这里h表示最大三角形边长。当你把网格加密一倍速度误差应该缩小大约4倍压力L^2误差缩小4倍。做数值实验时这个按阶收敛的规律能帮你判断程序写对了没有。常见做法是画一张网格加密表h1/16、1/32、1/64分别计算速度误差再比较误差比。如果你的误差不随加密而下降说明边界条件、积分阶次或空间选择有问题不是网格不够细。3.2 边界条件的数学与实现无滑移、顶盖驱动和压力锁定方腔驱动流的三类边界条件都滑到Dirichlet类型顶盖速度给定为(1,0)其余边速度为零。代码里的noslip边界表达式on_boundary !near(x[1], 1.0)必须排除顶盖否则顶盖同时被强制为零驱动就没了。这是最常见的低级错误之一后面避坑章详谈。压力在Dirichlet边界上不需要给定条件但纯Neumann问题里的压力常数自由度会让线性系统奇异。FEniCS里要么固定一个压力点要么施加零均值约束。固定点的位置有讲究不要放在顶盖与边墙的角点因为那个地方速度不连续压力有奇异性放在左下角附近相对安全。如果你用零均值约束则在变分形式里额外加一项∫p dx 0然后用Lagrange乘子处理但FEniCS中点固定更省事。注意角点奇异性顶盖与左右墙交界处速度从0突变为1数学上不可微。任何一阶精度以上的有限元在这里都会出现压力峰值增大现象。不要把这个当程序bug加密网格后峰值会继续增长这是真实解的弱奇异特性。处理技巧是后处理时避开角点或者接受离开角点两个网格单元后结果可信。3.3 网格收敛性验证从粗网格到细网格第一次跑通程序后不要直接看云图先做网格收敛性验证。以Re100为例分别在16、32、64网格密度下求解然后提取y0.5水平中线上u速度的分布。如果三条曲线基本重合说明网格够了如果明显分离继续加密。提取中心线数据的代码可以这样写import numpy as np def probe_u(u_func, y00.5, n33): xs np.linspace(0, 1, n) values [] for x in xs: # FEniCS中Function可直接用Point调用 ui u_func(x, y0) # 返回(u, v) values.append(ui[0]) return xs, np.array(values) xs, u_line probe_u(u, 0.5, 33)这个函数在网格加密后重复调用把速度曲线画在一起对比。Ghia等人1982年在方腔驱动流基准解中给出过Re100、400、1000等工况的中线速度数据把曲线和文献点放在同一张图里偏差肉眼可见。代码逻辑说明u_func是FEniCS中的Function对象可以直接传入坐标点计算插值结果。ui[0]提取x方向速度分量。采样点数n不用太多33个点足够画出一条光滑曲线。加密网格后如果曲线抖动剧烈不是采样问题而是网格太粗导致的数值振荡。4. 稳态求解的迭代策略与参数调优4.1 Picard迭代与Newton迭代什么时候选哪个稳态N-S方程的对流项是非线性的非线性迭代是求解器的核心。最常见的两种线性化方法是Picard迭代和Newton迭代。Picard的思路是用上一步速度结果替代对流项里的速度系数把非线性问题改成每步一个Stokes-like线性问题。它的好处是对初值不敏感从零初值也能收敛但收敛速率是线性的需要较多迭代步。Newton迭代则对弱形式做Frechet求导得到Jacobian矩阵每次迭代都修正非线性残差。它收敛快二阶速率但初值必须足够接近真解否则残差振荡甚至直接发散。解决初值依赖性的标准手段是续延法从一个小雷诺数比如Re1开始求解以这个解作为Re10的初值一步步推到目标雷诺数。FEniCS中的NonlinearVariationalProblem默认使用Newton代码里通过J derivative(F, w)提供了Jacobian。如果你的程序在高Re下发散不要慌先换成续延策略。续延实现很直接因为nu是Constant对象可以在循环中修改它的值nu_vals [1.0, 0.1, 0.01] # 对应 Re1, 10, 100 for val in nu_vals: nu.assign(1.0 / val) problem.solve()逻辑说明nu.assign()更新黏度值FEniCS会把新值传播到所有用到nu的表达式里。上一个雷诺数收敛的w自动成为当前雷诺数的初始猜测。这个循环是续延法的核心比手动调初值可靠得多。4.2 三个必调参数松弛因子、容忍度、最大迭代次数如果你是手写Picard循环或者自研MATLAB程序需要关注三个参数。第一是松弛因子记为omega。每步Picard求解得到新速度u_new后不直接替换而是做线性混合w.vector()[:] omega * u_new.vector() (1 - omega) * w.vector()omega取0.7到0.9比较常见。omega越小越稳但收敛越慢omega太大容易振荡。对稳态N-S我一般从0.7起步如果残差单调下降就逐步提到0.9。第二是容忍度。绝对容忍度设1e-8是稳妥值但注意它只针对绝对残差向量范数。如果你只关心工程意义上的速度分布1e-6就能出图了但要做网格收敛阶验证至少要1e-8否则误差都落在迭代未收敛的残差里。第三是最大迭代次数。Newton设20次足够Picard可能需要50次。如果达到上限仍未收敛说明初值不合适不要盲目加迭代次数先看雷诺数和网格是否匹配。4.3 线性求解器与预处理别忘了鞍点系统的特殊性每次非线性迭代内部都要解一个大型线性方程组。N-S离散出的矩阵是鞍点结构不是普通的正定矩阵用CG这类针对正定矩阵的方法会直接失败。推荐的做法是小网格用直接求解器MUMPS大网格用GMRES配合块预处理。FEniCS中可以指定solver.parameters[linear_solver] mumps如果网格超过10万节点MUMPS内存会膨胀这时改用solver.parameters[linear_solver] gmres solver.parameters[preconditioner] amg逻辑说明amg是代数多重网格预条件对速度块有效但压力块还需要额外处理。FEniCS内置的amg配合hypre在大多数二维算例上表现不错。如果你看到GMRES残差卡在1e-3下不去八成是压力常数自由度没锁住或者预条件子没匹配好。这个坑在第5章还会展开。预条件器的选择有一些经验成分但有一条铁律不要用ilu做N-S的预条件至少对几百阶的矩阵不管用。二维方腔在64×64网格下自由度数大约1万多ilu填充率失控运行时间会翻几倍。5. 避坑二维稳态N-S有限元求解的5个常见翻车现场5.1 现象压力云图呈棋盘状原因速度空间和压力空间都用了P1不满足LBB条件。解决改用Taylor-Hood元也就是速度P2、压力P1。检查空间定义V VectorFunctionSpace(mesh, P, 2) Q FunctionSpace(mesh, P, 1)如果已经用了P2/P1仍出现棋盘压力那可能是网格太粗尤其是Re较高时。先加密网格再看压力云图是否改善。5.2 现象线性求解器报错崩溃残差卡住不降原因压力常数自由度没有固定。压力在方程中只以梯度出现整体加一个常数仍满足方程系统奇异。解决在代码里加pin DirichletBC(W.sub(1), Constant(0), x[0] DOLFIN_EPS x[1] DOLFIN_EPS, methodpointwise)。如果用了MUMPS奇异性会直接报错用GMRES则表现为残差长期不降。加了pin之后两者都会恢复正常。5.3 现象牛顿迭代发散残差振荡或者出现NaN原因雷诺数过高零初值离真解太远或角点奇异性导致局部残差过大。解决不要硬扛用续延法。从Re1逐步增加到目标值。如果即使Re很小也发散检查边界条件是否把顶盖速度意外设成了零。5.4 现象流动没有被驱动起来整个流场接近静止原因无滑移边界条件覆盖了顶盖。常见写法是noslip DirichletBC(W.sub(0), Constant((0,0)), on_boundary)这样顶盖也被强制为零。解决排除顶盖区域用!near(x[1], 1.0)。注意这个Bug不会报错只会得到一个平庸解非常容易忽略。5.5 现象加密网格后速度曲线出现显著偏差原因粗网格下流场结构没被解析特别是高雷诺数下边界层和涡心位置偏移。解决对网格做收敛性检查。可以对比16、32、64网格中心线速度曲线如果两条加密后的曲线重合就认为网格无关。如果一直不重合检查求解器容忍度是否足够迭代是否真正收敛。6. 验证程序用基准解和网格收敛阶给求解器上保险程序跑通后最重要的不是急着出图而是验证它算得对。方向有两个一是和公开基准解对比二是验证误差按理论阶次收敛。基准解对比最简单的做法是使用方腔驱动流文献数据。Ghia等在1982年发表的Re100、400、1000方腔流数据是经典参考。你把程序中y0.5中心线上的u速度和文献离散点画在一个坐标系如果趋势一致且数值差距在网格加密后不断缩小就可以认为求解器正确。注意文献数据是特定网格加密到极致的参考解粗网格有偏差是正常的。第二个验证是误差阶次检验。如果你有解析解比如用人为构造的方法制造解技术直接计算L2误差。但方腔流没有解析解可以用网格收敛阶近似取三套网格h、h/2、h/4比较相邻网格上解的差别。P2/P1元理论预期速度L2误差阶是2压力L2误差阶也是2。如果实际误差比明显偏离比如只有阶1.3那说明程序里有隐藏问题比如积分阶次不够或边界条件施加错误。下面这段代码可以帮你做收敛阶粗测import numpy as np def compute_velocity_l2_diff(u_fine, u_coarse): # 把粗网格解插值到细网格自由度上计算L2范数差 V_fine u_fine.function_space() P interpolate(u_coarse, V_fine) diff u_fine - P return np.sqrt(assemble(inner(diff, diff) * dx)) # 假设有16和32网格的解误差比约为4说明二阶 err1 compute_velocity_l2_diff(u32, u16) err2 compute_velocity_l2_diff(u64, u32) ratio err1 / err2 print(收敛阶约为:, np.log2(ratio) if ratio 0 else 误差不降)这段代码利用了FEniCS的interpolate把粗网格解转换到细网格空间再计算差函数的L2范数。若ratio接近4说明误差以O(h²)下降若接近2只是一阶说明压力或空间选取有问题。说说我自己的教训早先我算Re1000时为了省时间只用了32×32网格结果涡心位置和文献差了10%以上。返工到64×64才基本吻合白白浪费一个晚上。所以先用粗网格检查程序逻辑再用中等网格看趋势最后加密到目标网格这个过程不要省。验证不是走形式它是给自己上的保险。希望帮到你。本文还有配套的精品资源点击获取