Stokes 问题
关键词： 稳态 线性 Stokes方程 Taylor-Hood单元
在流体力学领域，对于 Stokes 方程的混合有限元方法的研究一直是一个热点问题。通常的有限元法的求解困难在于：Stokes 方程要求有限元空间的组合必须满足 Ladyzhenskaya–Babuška–Brezzi（LBB）（或 inf–sup）相容性条件。正是这一条件的限制排除了传统的等阶插值有限元空间的使用。
求解 Stokes 方程使用最广泛的有限元族之一是 Taylor–Hood 族。 它由用于速度场的连续 $P_q (q \ge 2)$ 拉格朗日单元和用于压力场的连续 $P_{q−1}$ 拉格朗日单元组成。对于 $q = 2$，一维和二维的 Taylor–Hood 单元如下图所示。
本节以二维 Stokes 问题和 q9-q4 四边形单元为例，讲解如何使用 FEtch 系统提供的 Taylor–Hood 单元对流体力学问题进行全耦合求解。
控制方程
微分形式
Stokes 方程描述了稳定的不可压缩的低雷诺数流体。 对于域 $\Omega\subset\mathbb{R}^d\ (1\le d \le 3)$，Stokes 方程为：
\[-\mu\Delta \boldsymbol{u} + \nabla p = \boldsymbol{f} \ \left(\text{in } \Omega\right) \tag{1a}\label{eq1}\] \[\nabla \cdot \boldsymbol{u} = 0 \ \left(\text{in } \Omega\right) \tag{1b}\]其中，$\boldsymbol{u}:\Omega \to \mathbb{R}^d$ 是速度场，$p:\Omega \to \mathbb{R}$ 是压力场，$\boldsymbol{f}:\Omega \to \mathbb{R}^d$ 是源项， $\mu$ 是黏性系数。
通常仅考虑第一类边界条件，即在边界上给定速度和压力
\[\boldsymbol{u}=\bar{\boldsymbol{u}},\ p=\bar{p}\ \left(\text{on } \Gamma\right)\]弱形式
根据变分原理，Stokes 方程的弱形式如下
\[\int_{\Omega}\mu \nabla \boldsymbol{u} \cdot \nabla \delta \boldsymbol{u} d{\Omega}-\int_{\Omega}p(\nabla \cdot \delta \boldsymbol{u}) d{\Omega}-\int_{\Omega}(\nabla \cdot \boldsymbol{u}) \delta p d{\Omega}=\int_{\Omega} \boldsymbol{f} \cdot \delta u d{\Omega}\]在二维直角坐标系下，令 $\boldsymbol{u}=(u,v)$，$\boldsymbol{f}=(f_x,f_y)$，上式可写成
\[\int_{\Omega}{\mu \left (\frac{\partial u}{\partial x}\frac{\partial \delta u}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial \delta u}{\partial y}+\frac{\partial v}{\partial x}\frac{\partial \delta v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial \delta v}{\partial y}\right)d{\Omega}}\\-\int_{\Omega}{\left (p(\frac{\partial \delta u}{\partial x}+\frac{\partial \delta v}{\partial y})+(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y})\delta p\right)d{\Omega}}=\int_{\Omega}{\left({f_{x}}\delta u+{f_{y}}\delta v\right)d{\Omega}} \tag{2}\label{eq2}\]脚本编写
一个具体问题的有限元语言表述，通常需要包含 5 种类型的文件：PDE、FBC、MDI、GCN 和 NFE。由于方程 $\eqref{eq2}$ 不存在边界积分项，所以这里不需要写 FBC 文件。下面依次给出其他 4 个文件的详细介绍。
PDE 文件
在这个例子中，采用全耦合的方法，将流速和压力视为同一个物理场。根据有限元语言的语法规则，将方程 $\eqref{eq2}$ 的体积分项写入 PDE 文件如下。
| uvp.pde | 说明 |
|---|---|
| defi disp u v p coor x y shap %1 %2 shap %1 4 p gaus %3 func evv mate emu fu fv 1.0; 0.0; 0.0 func evv=+[u/x]+[v/y] stif dist= +[u/x;u/x]*emu+[u/y;u/y]*emu +[v/x;v/x]*emu+[v/y;v/y]*emu -[p;evv]-[evv;p] load=+[u]*fu+[v]*fv end |
系统关键字（信息定义） 定义未知量 定义坐标分量 定义采用的形函数 特别定义 p 变量的形函数由 4 个节点插值 定义采用的积分方式 声明自定义函数 定义用到的材料参数和默认参数值 （空行） 系统关键字（向量函数表达式） 自定义函数的具体表达式 （空行） 系统关键字（单元刚度矩阵表达式定义） 刚度矩阵项，对应式 $\eqref{eq2}$ 左侧 （空行） 荷载向量项，对应式 $\eqref{eq2}$ 右侧 （空行） 系统关键字（文件结束） |
这里将 shap 行和 gaus 行的参数设为待定，参数值 %1、%2 和 %3 将从 MDI 文件中取得。
注意，脚本文本中存在两个 shap 行。多出来的第二个 shap 行是 Taylor-Hood 单元所特有的，用于对未知量 p 的形函数进行特别定义，指定其为 4 节点的线性插值。
MDI 文件
MDI 文件负责组织 PDE 和 FBC 文件，声明程序具体使用的坐标系统、单元类型和积分方法。这里文件命名为 stokes.mdi，具体内容如下。
| stokes.mdi | 说明 |
|---|---|
| 2dxy #a 0 3 u v p_ pde uvp q9g3 # |
在二维直角坐标系下求解 a 场，无初值，3 个自由度，名称为 u v p，p 采用低阶插值 调用 uvp.pde，采用 9 节点四边形单元，3 点高斯积分 （文件结束符，可留空） |
特别注意，自由度 p 的右侧加了下划线 “_”，代表它采用的形函数的阶次要比当前的单元阶次低一阶。同时，下划线 “_” 的使用一定要和 PDE 文件中的 shap 行匹配好，否则会出现致命错误。
MDI 文件声明的单元类型和积分方式将以参数的形式传递给 PDE 文件中的的 %1、%2 和 %3。如果后续想要使用其他的单元类型或者积分方法，只需修改 MDI 文件即可。
GCN 文件
GCN 文件给出计算流程信息，系统规定其文件名字要与 MDI 文件一致。作为线性稳态问题，本算例的计算流程非常简单，只需填写少量代码，具体内容如下。
| stokes.gcn | 说明 |
|---|---|
| defi a ell #startumf a #solvumf a |
算法引入段 a 场，采用的算法模板文件为 ell.nfe (空行) 开始命令流段 初始化 a 场 求解 a 场，使用 umf 求解器 |
注意，这里推荐使用 umf 求解器，它可以很好地保证计算的稳定性和精度。如果采用 sin 求解器，由于存在压力项矩阵的对角位置为零的情况，计算通常是不收敛的。此时需要对控制方程添加稳定项，进而需要改写 PDE 文件。感兴趣的高级用户可自行尝试实现。
NFE 文件
GCN 文件使用了 ell 算法，系统会基于模板文件自动生成以下的 stokesa.nfe 文件。这里需要修改的是 gidpost 输出语句：
- 速度（矢量）和压力（标量）在物理本质上是两个场，在将结果保存到 GiD 后处理文件时需要进行区分，即采用两次 gidpost 语句，分开输出。
修改完成后的文件内容如下。
| stokesa.nfe | 说明 |
|---|---|
| defi stif s load f equation matrix = [s] forc = [f] solution u gidpost("u v",1.0) u(1),u(2) gidpost("p",1.0) u(3) end |
DEFI段 定义单元刚度矩阵 s 定义荷载向量 f （空行） EQUATION段 总刚矩阵公式 荷载向量公式 （空行） SOLUTION段 定义求解的场向量 将流速结果写入 GiD 后处理文件 将压力结果写入 GiD 后处理文件 （空行） 结束标记 |
算例
问题描述
考虑一个二维的 $1\ \mathrm{m}×1\ \mathrm{m}$ 的方腔，内部流体的初始速度和压力都是零，现在顶部一无限长的滑板以 $1 \mathrm{m/s}$ 的速度向右拖动顶部流体，求稳定后方腔内部的流场速度和压力分布情况。
设流体的黏性系数 $\mu =1.0\times 10^{-3}\ \mathrm{Pa\cdot s}$，体积力 $f_x=f_y=0.0$。
边界条件为：上面边界上的流体速度 $u=1.0$，$v=0$。其他三边上都为固壁边界条件 $u=0$，$v=0$。在 $(0,0,0)$ 点处给定参考压力 $p=0$。
前处理
填写材料参数文件
系统自动生成了 stokes.mat 文件，用于设置材料参数。参数初始值全部为 PDE 文件给出的默认值。可根据具体问题进行修改。对于本算例，修改完成后的文件内容如下。
| stokes.mat | 说明 |
|---|---|
| 1 4 aeq9g3 num emu fu fv 1 1.0 0.0 0.0 # |
9 节点四边形面单元 材料参数 参数值 分隔符 |
赋单元材料
施加第一类边界条件
其中：
u-I 为边界类型，-1 表示约束，1 表示自由
u-D 为相应的约束值或施加节点力
由于有不同的边界条件交接，角点位置上的边界条件可能出现歧义。为此，继续明确定义角点处的边界条件。
网格剖分
根据所开发程序的计算要求，网格剖分采用 9 节点等参单元。共使用了 400 个单元，1681 个节点。
计算结果
水平流速分量的分布
竖直流速分量的分布
流速大小的分布
压力的分布
如果继续加密网格进行计算分析，可以进一步提高计算精度，逼近真解，从而验证 Taylor-Hood 单元和计算程序的有效性。这里不再给出具体过程，留给用户自行测试。
感谢您的支持，我们会继续努力的!
打开支付宝或微信扫一扫，即可进行扫码打赏哦
