INSE学习笔记——5. 投影法
投影法也称为时间分裂法，分步法，按照思路分为两类：
- 第一类基于 N-S 方程组的半离散格式，基于 Helmholtz 分解定理；
- 第二类基于 N-S 方程组的全离散格式，采用近似算子分裂/近似矩阵分裂。
下文只包括第一类投影法。
Helmholtz 分解定理
对于任意光滑的向量场\(V\)可以唯一分解为一个无源场\(V^d\)和一个无旋场\(\nabla \phi\) 之和。
\[ V = V^d + \nabla \phi \]
满足\(\nabla \cdot V^d = 0\), 和\(\nabla \times \nabla \phi=0\)。
(这一说法是不严谨的，因为唯一性需要向量场满足一定的边界条件，暂时略过边界的考虑)
投影法简介
考虑无量纲化的二维不可压缩方程组
\[ \begin{aligned} \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot \nabla) \mathbf{u} &= -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u}\\ \nabla \cdot \mathbf{u} &= 0\\ \end{aligned} \]
1.无压力增量投影法
我们介绍最早由 Chorin 和 Temam 提出的投影法，利用显式欧拉离散得到：
\[ \begin{aligned} \frac{\mathbf{u}^{n+1}-\mathbf{u}^n}{\Delta t}+ \left[(\mathbf{u}^n\cdot \nabla) \mathbf{u}^n - \frac{1}{Re} \nabla^2 \mathbf{u}^n \right] &= -\nabla p^{n+1} \\ \nabla \cdot \mathbf{u}^{n+1} &= 0\\ \end{aligned} \]
在其中引入预测速度\(\mathbf{u}^*\)，把动量方程拆成两个部分：
\[ \begin{aligned} \frac{\mathbf{u}^{*}-\mathbf{u}^n}{\Delta t}+ \left[(\mathbf{u}^n\cdot \nabla) \mathbf{u}^n - \frac{1}{Re} \nabla^2 \mathbf{u}^n \right] = 0 \\ \frac{\mathbf{u}^{n+1}-\mathbf{u}^*}{\Delta t} = -\nabla p^{n+1} \end{aligned} \]
对第二个式子取散度，得到
\[ \nabla^2 p^{n+1} = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^* \]
因此，得到计算流程如下：
- 从\(t^n\)时刻的第 n 步 \(\mathbf{u}^n\) 出发，
- 基于\(\mathbf{u}^n\)，计算预测速度\(\mathbf{u}^*\)满足 \[ \begin{aligned} \frac{\mathbf{u}^{*}-\mathbf{u}^n}{\Delta t}+ \left[(\mathbf{u}^n\cdot \nabla) \mathbf{u}^n - \frac{1}{Re} \nabla^2 \mathbf{u}^n \right] = 0 \end{aligned} \]
- 基于\(\mathbf{u}^*\)，计算压强\(p^{n+1}\)满足 \[ \nabla^2 p^{n+1} = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^* \]
- 得到第 n+1 步的速度\(\mathbf{u}^{n+1}\) \[ \mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t \nabla p^{n+1} \]
最后一步体现了Helmholtz 分解定理和投影：
\[ \mathbf{u}^{*} = \mathbf{u}^{n+1} + \Delta t \nabla p^{n+1} \]
- 在利用 \(\mathbf{u}^n\) 计算预测速度 \(\mathbf{u}^*\) 时，直接忽视了压强场的约束作用；
- 预测速度\(\mathbf{u}^*\)根据 Helmholtz 分解定理被分解成无源场 \(\mathbf{u}^{n+1}\) 和无旋场 \(\Delta t \nabla p^{n+1}\) ；
- 通过预测速度 \(\mathbf{u}^*\) 和压强 \(p^{n+1}\) 计算 \(\mathbf{u}^{n+1}\) 就是通过压强的约束作用，把速度重新投影到散度为零的子空间内：\(\mathbf{u}^{n+1}=P(\mathbf{u}^*)\)，最终的速度严格满足\(\nabla \cdot \mathbf{u}^{n+1}=0\)。
这种最早的投影法通常称作无压力增量投影法：因为计算中间层的预测速度\(\mathbf{u}^*\)时，没有考虑压强梯度的影响。可以改成隐式的时间离散格式，可以在计算预测速度时考虑压强的影响，因此精度有限。有一系列衍生的投影法。
注意，实践中可能不是直接离散\(\nabla^2 p^{n+1} = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^*\)，而是离散下面两个式子，把离散后的结果进行合并：
\[ \begin{aligned} \frac{\mathbf{u}^{n+1}-\mathbf{u}^*}{\Delta t} = -\nabla p^{n+1}\\ \nabla \cdot \mathbf{u}^{n+1} =0 \end{aligned} \]
这样可以保证离散的 \(\nabla \cdot \mathbf{u}^{n+1} =0\) 严格成立。
2.压力增量投影法
首先，时间上不使用显式 Euler，而改用二阶向后差分格式
\[ \frac{\partial \mathbf{u}^{n+1}}{\partial t} \approx \frac{3\mathbf{u}^{n+1}-4\mathbf{u}^n + \mathbf{u}^{n-1}}{2\Delta t} \]
动量方程离散为
\[ \begin{align} \frac{3\mathbf{u}^{n+1}-4\mathbf{u}^n + \mathbf{u}^{n-1}}{2\Delta t} + (\mathbf{u}^{n+1}\cdot \nabla) \mathbf{u}^{n+1} &= -\nabla p^{n+1} + \frac{1}{Re} \nabla^2 \mathbf{u}^{n+1}\\ \end{align} \]
(不等价地)拆为下面两个式子
\[ \begin{align} \frac{3\mathbf{u}^{*}-4\mathbf{u}^n + \mathbf{u}^{n-1}}{2\Delta t} + (\mathbf{u}^{*}\cdot \nabla) \mathbf{u}^{*} &= -\nabla p^{n} + \frac{1}{Re} \nabla^2 \mathbf{u}^{*}\\ \frac{3(\mathbf{u}^{n+1}-\mathbf{u}^*)}{2\Delta t} &= -\nabla (p^{n+1}-p^n) \end{align} \]
第一步，预测步计算中间速度\(\mathbf{u}^*\)，考虑压强梯度的影响(需要隐式求解\(\mathbf{u}^*\))
\[ \begin{align} \frac{3\mathbf{u}^{*}-4\mathbf{u}^n + \mathbf{u}^{n-1}}{2\Delta t} + (\mathbf{u}^*\cdot \nabla) \mathbf{u}^* &= -\nabla p^{n} + \frac{1}{Re} \nabla^2 \mathbf{u}^* \end{align} \]
第二步，求解压强增量\(q = p^{n+1}-p^n\)，对第二个式子求散度得
\[ \begin{align} \nabla^2 q &=\frac{3}{2\Delta t} \nabla \cdot \mathbf{u}^* \end{align} \]
第三步，更新\(\mathbf{u}^{n+1}\), \(p^{n+1}\)
\[ \begin{align} \mathbf{u}^{n+1} &= \mathbf{u}^* - \frac{2\Delta t}{3} \nabla q\\ p^{n+1} &= p^n + q \end{align} \]
3.旋度型压力增量投影法
第一步，预测步, 计算中间速度\(\mathbf{u}^*\)，考虑压强梯度的影响(需要隐式求解\(\mathbf{u}^*\))
\[ \begin{align} \frac{3\mathbf{u}^{*}-4\mathbf{u}^n + \mathbf{u}^{n-1}}{2\Delta t} + (\mathbf{u}^*\cdot \nabla) \mathbf{u}^* &= -\nabla p^{n} + \frac{1}{Re} \nabla^2 \mathbf{u}^*\\ \end{align} \]
第二步，从求解压强增量\(q=p^{n+1}-p^n\) 变成求解新的变量\(\phi\)
\[ \phi = p^{n+1}-p^n + \frac{1}{Re} \nabla \cdot \mathbf{u}^* \]
满足如下方程组
\[ \begin{align} \nabla^2 \phi &=\frac{3}{2\Delta t} \nabla \cdot \mathbf{u}^* \end{align} \]
第三步，更新\(\mathbf{u}^{n+1}\), \(p^{n+1}\)
\[ \begin{align} \mathbf{u}^{n+1} &= \mathbf{u}^* - \frac{2\Delta t}{3} \nabla \phi\\ p^{n+1} &= p^n + \phi - \frac{1}{Re} \nabla \cdot \mathbf{u}^* \end{align} \]
这一种和压力增量型投影法相比，略有不同，可以视作对压力增量型投影方法的改进，因为如果不可压方程组对速度有第一类边界条件时，对\(\phi\)的求解比对\(q\)的求解可以给出更合理的边界条件，略。