二维非定常纳维-斯托克斯流体有限元求解器

资 源 简 介

该项目包含名为UNSTEADY_NAVIER_STOKES的MATLAB脚本，旨在通过有限元方法求解二维纳维-斯托克斯（Navier-Stokes）方程。虽然源描述中提及稳态，但根据脚本名称及常规数值计算逻辑，该工具主要架构用于处理非定常流体流动问题，或通过时间演化方法求解稳态解。其空间离散化完全采用标准的伽辽金（Galerkin）方法。在具体的有限元网格与插值函数选择上，项目采用了混合单元策略：对于压力场，使用双线性四边形单元（Bilinear Quadrangular Elements）进行近似；对于速度场，则使用高一阶的双二次四边形单元（Biquadratic Quadrangular Elements）进行近似。这种这种速度与压力的特定插值组合（也就是经典的Taylor-Hood单元配置）严格满足LBB（Ladyzhenskaya-Babuska-Brezzi）条件。满足LBB条件是至关重要的，因为它确保了离散化后得到的代数方程组具有唯一解，并保证了数值计算的稳定性，防止在压力解中出现虚假的数值震荡。程序的主要功能包括构建网格、组装质量矩阵与刚度矩阵、处理非线性对流项、应用边界条件以及求解最终的线性系统。

详 情 说 明

2D Unsteady Navier-Stokes Finite Element Solver

项目简介

本项目是一个基于MATLAB开发的二维非定常纳维-斯托克斯（Navier-Stokes）方程有限元求解器。该程序专为模拟不可压缩流体的非定常流动而设计，默认配置为经典的顶盖驱动方腔流（Lid-driven Cavity）物理模型。

求解器核心采用标准的伽辽金（Galerkin）有限元方法，通过混合单元策略保证了数值计算的稳定性。代码完全独立，内部集成了网格生成、矩阵组装、边界条件处理、非线性迭代及时间步进求解的全流程，无需依赖外部网格文件或第三方工具箱。

功能特性

• 物理模型：针对不可压缩牛顿流体的二维非定常流动进行求解。
• 空间离散：采用Q2-Q1 Taylor-Hood混合单元配置。

* 速度场 ($u, v$)：使用9节点双二次四边形单元（Q2）近似，精度高。 * 压力场 ($p$)：使用4节点双线性四边形单元（Q1）近似。 * 稳定性：该组合严格满足LBB（Ladyzhenskaya-Babuska-Brezzi）条件，有效防止压力场出现伪震荡。

• 时间离散：采用隐式时间步进方案（Implicit Time Stepping），保证了时间推进的数值稳定性。
• 非线性求解：内置Picard迭代（定点迭代）算法，用于处理纳维-斯托克斯方程中的非线性对流项，在每个时间步内通过迭代收敛至真实解。
• 全耦合求解：采用单块（Monolithic）方法构建全局线性系统，同时求解速度和压力自由度。
• 自包含网格生成：脚本内部通过算法自动生成结构化四边形网格及对应的节点拓扑关系。

系统要求

• MATLAB R2016b 或更高版本。
• 无需额外的工具箱（Toolbox），仅依赖MATLAB基础函数（如 mesh, sparse, ` 等）。

使用方法

直接在MATLAB环境中运行主脚本（即包含 main 函数的文件）。程序将依次执行以下步骤：

1. 初始化物理参数（雷诺数、时间步长等）。
2. 生成有限元网格。
3. 组装系统矩阵。
4. 进入时间步循环，实时求解流场并在图形窗口中动态绘制结果。

代码实现与算法详解

本项目核心逻辑在 main 函数中实现，其具体流程如下：

1. 参数定义与初始化

程序首先定义了仿真的关键物理与数值参数，默认设置为：

• 雷诺数 (Re)：100
• 时间步长 (dt)：0.05
• 网格密度：8x8 个单元（指宏观单元数，实际自由度更高）
• 迭代控制：设定了非线性Picard迭代的最大次数和收敛容差。

2. 网格生成与拓扑构建

不同于读取外部网格，代码内部实现了结构化网格生成器：

• 节点生成：在 $[0,1] times [0,1]$ 区域内生成节点。由于采用混合单元，速度节点密度大约是压力节点的4倍。
• 连通性矩阵：

* 构建 elem_u：记录每个Q2单元的9个速度节点索引，遵循角点-边中点-中心点的节点排序。 * 构建 elem_p：记录每个Q1单元的4个压力节点索引（仅角点）。

3. 自由度（DOF）管理

系统构建了全局解向量 U_global，其结构为堆叠形式：$[u, v, p]^T$。

• 前 $N$ 个分量为水平速度 $u$。
• 中间 $N$ 个分量为垂直速度 $v$。
• 最后 $M$ 个分量为压力 $p$。

4. 边界条件定义

代码硬编码了顶盖驱动方腔流的边界条件：

• 速度边界：

* 上边界（$y=1$）：设定 $u=1, v=0$。 * 左、右、下边界：设定无滑移条件 $u=0, v=0$。

• 压力边界：为了消除不可压缩流动中压力场常数模式导致的奇异性，固定了第一个压力节点的值为0。

5. 矩阵组装（恒定部分）

在时间循环开始前，预先计算不随时间变化的矩阵以提高效率：

• M (质量矩阵)：源于时间导数项 $partial u / partial t$。
• K (刚度矩阵/拉普拉斯矩阵)：源于粘性扩散项 $nabla^2 u$。
• Dx, Dy (梯度/散度算子)：源于压力梯度项 $nabla p$ 和连续性方程 $nabla cdot u = 0$ 中的混合积分。
• 矩阵组装使用高斯积分法（Gaussian Quadrature），在子函数 assemble_constant_matrices 中实现。

6. 非定常与非线性求解循环

程序核心是一个双层循环结构：

• 外层循环（时间步）：推进物理时间。
• 内层循环（Picard迭代）：解决当前时间步内的非线性耦合。

1. 更新对流矩阵：基于上一次迭代的速度场 $u_k$，重新组装对流矩阵 $C(u_k)$。 2. 构建线性系统： * 动量方程部分包含：惯性项（$M/dt$）、粘性项（$K/Re$）、对流项（$C$）和压力梯度项（$Dx^T, Dy^T$）。 * 连续性方程部分包含：散度约束（$-Dx, -Dy$）。 3. 应用边界条件：通过修改线性系统矩阵（置1法/替换法），将狄利克雷边界条件强行施加到线性方程组中。 4. 求解：利用MATLAB的稀疏矩阵直接求解器（`）计算新的解向量。 5. 收敛判定：比较两次迭代解的范数，若残差小于设定容差则进入下一时间步。

7. 可视化

在时间步进过程中，程序会定期调用绘图函数，将计算得到的全局解向量还原到网格节点上，直观展示流场的演化过程。