企业官网建设流程全解析

简介：这份资源是面向流体力学与数值计算学习者的二维稳态Navier-Stokes方程有限元求解程序，基于Matlab实现，适合已具备偏微分方程与有限元基础、希望深入理解流体数值求解流程的高年级本科生、研究生及工程技术人员。程序围绕动量方程与连续性方程展开，涵盖几何离散、函数空间定义、变分形式推导、矩阵组装、线性系统求解及速度压力后处理等完整环节，并包含Dirichlet、Robin、应力边界处理与Newton迭代初始化等模块。压缩包共45个文件，全部为m脚本，约19KB，按功能可大致分为网格生成、基函数与高斯积分、矩阵与向量组装、边界条件处理、误差计算及求解测试等类别，结构紧凑便于逐文件研读。目前已有312人学习下载。通过阅读与调试源码，读者可掌握有限元法求解流体问题的实现细节，并在此基础上迁移至热传导、扩散等方程，提升Matlab编程与算法定制能力。

1. 二维稳态Navier-Stokes方程有限元求解：从“跑不动”到“看得见流场”的完整程序拆解

很多人第一次接触流体计算，都是被那个非线性项坑到怀疑人生。二维稳态Navier-Stokes方程本身不算长，但真拿有限元去离散，迎风、迭代、压力震荡、网格依赖，随便一个都能让程序卡在“残差降不下去”的状态。这份资源给的是一个完整的二维稳态N-S方程有限元求解程序，用Matlab写成，主打不可压缩牛顿流体在规则/非规则区域内的稳态流场求解。它适合两类人：一类是正在做流体力学课程设计、需要快速拿到速度场和压力场结果的学生；另一类是想把有限元从结构力学扩展到流体、需要一个能改能跑的代码底子的开发者。程序本身不依赖商业工具箱，核心就是刚度矩阵组装、非线性迭代和线性方程组求解三件事，但把这三件事串起来，就是一个能出图的求解器。

2. 程序架构与目录结构：先看清模块边界再动手

2.1 整体分层：从网格到后处理的五段式

这个程序的结构是典型的“脚本驱动 + 函数模块”混合风格。主脚本负责参数设定和流程调度，底层函数各管一件事。常见做法是分成五层：网格生成层、单元计算层、全局组装层、非线性求解层、后处理层。这样分的好处是，你改边界条件不用动组装逻辑，换迭代方法不用碰网格数据。

目录大致长这样：

ns2d_solver/ ├── main.m % 主入口，参数、调用、出图 ├── mesh/ │ ├── generate_rect.m % 矩形区域结构化网格 │ └── generate_obstacle.m % 带圆柱/方柱绕流的非结构网格 ├── fem/ │ ├── shape_Q1Q1.m % Q1-Q1 单元形函数及导数 │ ├── local_stiff.m % 单元刚度矩阵（粘性项） │ ├── local_convect.m % 单元对流矩阵（非线性项） │ └── assemble.m % 全局稀疏矩阵组装 ├── solver/ │ ├── newton_iter.m % Newton 迭代主循环 │ └── linear_solve.m % 直接法/迭代法切换 ├── post/ │ ├── plot_velocity.m % 速度矢量图、流线图 │ └── plot_pressure.m % 压力云图 └── data/ └── bc_config.mat % 边界条件配置

这个结构不是强制的，但如果你拿到的是一个散乱脚本，建议先按这个思路拆一遍。拆完你会发现，调试时能定位到具体层，而不是在一千行里找那个写错的节点编号。

2.2 网格与自由度编号：最容易埋雷的地方

二维稳态N-S用Q1-Q1单元（双线性速度和压力）会面临LBB条件不满足的问题，程序里通常用两种方式绕开：一是压力做稳定化处理，二是用Q2-Q1（速度二次、压力线性）。这份资源用的是Q1-Q1加压力稳定化，代码里体现在局部矩阵多了一个稳定化项。

自由度编号规则是：先排所有速度节点，再排压力节点。假设网格有N个节点，速度自由度是2N（u和v），压力自由度是N，总自由度3N。组装时全局矩阵是3N×3N的稀疏矩阵。这里有个血泪经验：节点编号和自由度编号的映射关系一定要在组装前算好，不要每次循环里现算，否则Matlab的循环开销会让你等到怀疑人生。

% 自由度映射：节点i的速度u、v和压力p在全局向量中的位置 nNodes = size(nodes,1); uDof = 1:2:2*nNodes; % u 自由度 vDof = 2:2:2*nNodes; % v 自由度 pDof = 2*nNodes+1:3*nNodes; % p 自由度 % 单元自由度索引（以Q1单元4节点为例） elemDof = zeros(12,1); for k = 1:4 nid = elemNodes(k); elemDof(2*k-1) = uDof(nid); elemDof(2*k) = vDof(nid); end elemDof(9:12) = pDof(elemNodes);

这段代码的关键是elemDof这个12×1的索引向量，它把局部单元矩阵的12个自由度映射到全局。参数说明：nodes是N×2的节点坐标矩阵，elemNodes是当前单元的4个节点编号。注意压力自由度排在最后，这样速度块和压力块在全局矩阵里是分开的，方便后面做分块预处理。

2.3 边界条件处理：Dirichlet和Neumann的落位

入口给速度（Dirichlet），出口给压力（Neumann），壁面给无滑移。程序里用一个大数组bcFlag标记每个自由度的类型：0表示自由，1表示Dirichlet，2表示Neumann。组装完全局矩阵后，对Dirichlet自由度做行替换——把对应行改成单位行，右端项改成给定值。这是最直接的做法，缺点是会破坏矩阵对称性，但对直接求解器影响不大。

% Dirichlet 边界条件施加 for i = 1:length(dirichletDof) dof = dirichletDof(i); val = dirichletVal(i); A(dof,:) = 0; A(dof,dof) = 1; b(dof) = val; end

逻辑说明：dirichletDof是预先收集好的Dirichlet自由度编号，dirichletVal是对应值。行替换后，该自由度对应的方程变成“u_i = val”，其他方程里这个自由度的影响通过列保留。注意这里没有做列消去，所以右端项不需要修正，但矩阵不再对称。如果你要用共轭梯度法，得改成对称消去法。

3. 核心模块拆解：从单元矩阵到Newton迭代的完整链路

3.1 单元刚度矩阵与对流矩阵的组装

粘性项对应的是拉普拉斯算子，单元刚度矩阵就是形函数导数乘积的积分。Q1单元在参考正方形上做2×2高斯积分，映射到物理单元。对流项是非线性的，因为速度本身是未知量，所以单元对流矩阵依赖于当前迭代步的速度值。

function [Ke, Ce] = local_matrices(uv, elemNodes, nodes, Re) % uv: 当前单元的速度值 [u1 v1 u2 v2 u3 v3 u4 v4] % 返回粘性刚度 Ke 和对流矩阵 Ce（12x12） [gaussPts, gaussWts] = gauss2x2(); Ke = zeros(12,12); Ce = zeros(12,12); for gp = 1:4 [N, dNdx, dNdy, detJ] = shape_Q1Q1(gaussPts(gp,:), elemNodes, nodes); w = gaussWts(gp) * detJ; % 粘性项: (1/Re) * grad(N) . grad(N) Ke = Ke + w * (1/Re) * (dNdx'*dNdx + dNdy'*dNdy); % 对流项: (u_h . grad) u_h uh = N' * uv(1:2:end); vh = N' * uv(2:2:end); Ce = Ce + w * N' * (uh*dNdx + vh*dNdy); end end

参数说明：Re是雷诺数，uv是当前单元8个速度自由度值，gaussPts和gaussWts是2×2高斯点和权重。shape_Q1Q1返回形函数矩阵N（4×1）和导数矩阵dNdx、dNdy（4×1），以及雅可比行列式detJ。注意Ke和Ce都是12×12，但只有速度块非零，压力块在稳定化项里单独加。

3.2 压力稳定化：让Q1-Q1能用的关键一步

Q1-Q1不满足LBB条件，直接解会出现压力棋盘震荡。程序里用的是压力投影稳定化（PSPG）的简化版，在连续性方程里加一个与动量方程残差相关的项。代码上表现为压力块多了一个稳定化矩阵。

% 压力稳定化项（简化 PSPG） tau = h_elem^2 / (2*Re); % 稳定化参数，h_elem 为单元尺寸 for gp = 1:4 [N, dNdx, dNdy, detJ] = shape_Q1Q1(gaussPts(gp,:), elemNodes, nodes); w = gaussWts(gp) * detJ; % 稳定化矩阵加到压力-速度耦合块 Kstab = Kstab + w * tau * (dNdx*dNdx' + dNdy*dNdy'); end

逻辑说明：tau是稳定化参数，常见取法是单元尺寸平方除以2倍雷诺数。这个项加到压力方程里，相当于给压力场加了一点扩散，抑制震荡。注意tau不能太大，否则会过度抹平压力梯度，一般取0.1到1.0之间的系数乘上去。

3.3 Newton迭代：非线性项怎么收敛

稳态N-S是非线性方程组，程序用Newton法迭代。每一步需要组装残差向量和雅可比矩阵。雅可比矩阵由三部分组成：粘性刚度、对流矩阵、以及对流项对速度的导数（也就是“几何”部分）。

% Newton 迭代主循环 maxIter = 20; tol = 1e-6; for iter = 1:maxIter [R, J] = assemble_residual_jacobian(uv, p, mesh, Re); % 施加边界条件 [J, R] = apply_bc(J, R, bcFlag, bcVal); % 求解线性系统 delta = J \ (-R); % 更新解 uv = uv + delta(1:2*nNodes); p = p + delta(2*nNodes+1:end); % 检查收敛 if norm(delta) < tol fprintf('Newton 收敛于第 %d 步，残差 %.2e\n', iter, norm(delta)); break; end end

参数说明：R是残差向量，J是雅可比矩阵，delta是增量。收敛判据用的是增量范数，也可以改用残差范数。注意Newton法对初值敏感，程序里通常用Stokes解（忽略对流项）作为初值，或者用低雷诺数解逐步延拓到高雷诺数。

3.4 线性求解器选择：直接法还是迭代法

Matlab里最省事的是\，对3N×3N的稀疏矩阵，N几千以内直接法还能扛。N上万就得换迭代法，比如GMRES加ILU预处理。程序里留了一个linear_solve.m做切换。

function x = linear_solve(A, b, method) switch method case 'direct' x = A \ b; case 'gmres' setup.type = 'ilutp'; setup.droptol = 1e-3; [L, U] = ilu(A, setup); [x, flag, relres] = gmres(A, b, [], 1e-8, 200, L, U); if flag ~= 0 warning('GMRES 未收敛，相对残差 %.2e', relres); end end end

逻辑说明：direct直接调用Matlab的稀疏直接求解器，适合中小规模。gmres用不完全LU分解做预处理，适合大规模。注意droptol控制ILU的稀疏性，太小会导致L、U接近稠密，太大则预处理效果差。一般从1e-3开始试。

4. 避坑与排查：那些让残差降不下去的常见问题

4.1 压力棋盘震荡：现象、原因、解决

现象：压力云图出现黑白相间的棋盘格，速度场看起来正常但压力完全不可用。原因：Q1-Q1单元不满足LBB条件，压力自由度比速度自由度约束不足。解决：加压力稳定化项，或者换Q2-Q1单元。如果已经加了稳定化还有震荡，检查tau是不是太小，或者网格是不是太粗。

4.2 Newton迭代发散：残差越迭代越大

现象：前几步残差下降，然后突然反弹，最后NaN。原因：初值太远，或者雷诺数太高，非线性项占主导。解决：先用Stokes解做初值，然后用延拓法——从Re=1开始，每步乘1.5，用上一步的解做下一步初值。如果还不行，检查对流矩阵的符号是不是写反了。

4.3 出口边界条件设置错误：流场被“堵住”

现象：出口附近速度异常大或异常小，流线不平行。原因：出口给了速度Dirichlet而不是压力Neumann。不可压缩流出口应该让流体自由流出，给压力参考值（通常为0），速度由连续性方程决定。解决：检查bcFlag里出口节点的标记，压力自由度给Neumann，速度自由度给自由。

4.4 网格质量导致雅可比行列式为负

现象：组装时报错“detJ <= 0”，或者结果完全乱掉。原因：网格单元扭曲太严重，或者节点编号顺序不对（比如逆时针写成顺时针）。解决：检查网格生成函数，确保每个单元节点按逆时针排列。如果是非结构网格，用patch画出来看一眼，肉眼检查有没有翻转单元。

4.5 线性求解器不收敛：GMRES残差卡住

现象：GMRES迭代到最大步数还没到容差，或者直接报错。原因：矩阵条件数太大，或者ILU预处理不够。解决：先换直接法验证问题本身有没有解。如果直接法能解，说明是预处理问题，调小droptol或者换ilu的选项。如果直接法也解不了，检查雅可比矩阵是不是奇异——通常是边界条件没施加够，导致压力参考值缺失。

5. 二次开发与改造：从能跑到好用

5.1 换网格：从矩形到圆柱绕流

程序默认给的是矩形区域，想算圆柱绕流，改generate_obstacle.m。核心是生成非结构网格，然后在圆柱表面标记壁面边界。常见做法是用delaunay三角剖分，再删掉圆柱内部的单元。改完后注意检查圆柱表面的节点是不是都在边界上，否则无滑移条件施加不上。

% 生成带圆柱的网格（简化示意） [x, y] = meshgrid(linspace(0, 10, 80), linspace(0, 5, 40)); nodes = [x(:), y(:)]; % 删除圆柱内部节点 r = sqrt((nodes(:,1)-3).^2 + (nodes(:,2)-2.5).^2); nodes(r < 0.5, :) = []; % 三角剖分 elem = delaunay(nodes(:,1), nodes(:,2)); % 标记圆柱表面节点为壁面 wallNodes = find(abs(r - 0.5) < 0.05);

参数说明：圆柱中心在(3, 2.5)，半径0.5。wallNodes是壁面节点索引，后面在bcFlag里把这些节点的速度自由度设为1（Dirichlet，值为0）。注意三角剖分后要删掉圆柱内部的单元，否则会在固体区域算出虚假流动。

5.2 换迭代方法：从Newton到Picard

Newton法收敛快但初值敏感，Picard法收敛慢但稳健。改法很简单，把newton_iter.m里的雅可比矩阵换成只含粘性项和对流项，不含对流项对速度的导数。这样每步只解一个线性化的问题，相当于不动点迭代。

% Picard 迭代（对比 Newton） for iter = 1:maxIter [K, C] = assemble_stokes_convect(uv, mesh, Re); A = K + C; b = -C*uv + bcRHS; % 注意符号 uv_new = A \ b; if norm(uv_new - uv) < tol break; end uv = uv_new; end

逻辑说明：Picard法把对流项当成已知量，用上一步的速度算。收敛速度是线性的，但不容易发散。实际用的时候可以先用Picard跑几步，再用Newton加速。

5.3 加时间项：从稳态到非稳态

稳态程序改成非稳态，核心是加质量矩阵和时间离散。最简单的是后向欧拉，每步解一个稳态问题，但右端项多一个M*(u_n - u_{n-1})/dt。

% 后向欧拉时间推进 M = assemble_mass(mesh); % 质量矩阵 dt = 0.01; for t = 1:nSteps % 右端项加时间项 b_time = M * (uv - uv_old) / dt; [R, J] = assemble_residual_jacobian(uv, p, mesh, Re); R = R + b_time; delta = J \ (-R); uv_old = uv; uv = uv + delta(1:2*nNodes); p = p + delta(2*nNodes+1:end); end

参数说明：M是质量矩阵，dt是时间步长。注意时间步长受CFL条件限制，一般取网格尺寸除以最大速度的0.5倍以下。加时间项后，程序从稳态求解器变成非稳态求解器，可以算涡街脱落这类瞬态问题。

5.4 验证方法：用解析解和基准算例

改完程序怎么知道对不对？两个办法。一是用解析解，比如Poiseuille流，矩形区域，入口给抛物线速度，出口给压力0，壁面无滑移，稳态解就是抛物线。算出来和解析解对比，误差应该在网格尺寸的平方量级。二是用基准算例，比如方腔驱动流，Re=100时中心涡位置有标准值，对比一下就能判断程序对不对。

% Poiseuille 流验证 % 解析解: u(y) = 4*Umax*y*(H-y)/H^2, v=0, p = -8*Umax*mu*x/H^2 % 算完后取出口截面速度，和解析解对比 u_exact = 4*Umax*y_out*(H-y_out)/H^2; error = norm(u_num - u_exact) / norm(u_exact); fprintf('相对误差: %.2e\n', error);

逻辑说明：y_out是出口截面的y坐标，u_num是数值解。如果误差在1e-3到1e-2之间，说明程序基本正确。如果误差很大，先检查边界条件，再检查网格收敛性。

从那以后我每次改完求解器，都强制走一遍Poiseuille验证，不管多急都先跑这个。因为流体计算里，边界条件写错一个符号，结果可能看起来“像那么回事”，但一对比解析解就露馅。希望帮到你。

本文还有配套的精品资源，点击获取