Code of the Week #4: solve(UEqn == -fvc::grad(p))

rheoworks · 2018 年10 月 23 日 14:52

抱歉，本周的每周一行来晚了。

今天，我们要讲的是

solve(UEqn == -fvc::grad(p));

这一行程序来自于 icoFoam.C ，从字面上理解是求解方程组。

首先，来看一下括号里面各项的定义：

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

以类 fvVectorMatrix 定义了一个以 UEqn 为名的变量，括号内的各项已经分别在 code of the week #1: \nabla(UU) div(phi, U), #2: \dfrac{\partial U}{\partial t} ddt(U), #3 \nabla\cdot(\nu\nabla U) laplacian(nu, U) 分别进行了解读，其分别代表了方程等号左边的各项

\begin{equation} \frac{\partial U}{\partial t} + \nabla\cdot(UU) - \nabla\cdot(\nu\nabla U) = -\nabla p \end{equation} \tag{1}

对方程左边部分离散化后可以得到以下矩阵方程。

\begin{equation} A_\mathrm{P} \mathbf U_\mathrm{P}^r + \sum A_\mathrm{N} \mathbf U_\mathrm{N}^n - E_\mathrm{P}^n=0\end{equation} \tag{2}

其中 U^r_P 指待求网格的速度， U^n_N 指相邻网格已知速度值 (前一时间步)， E_P^n 为根据该点已知速度计算所得的数。这一方程加上 \nabla p 就成为动量预测方程。

\begin{equation} A_\mathrm{P} \mathbf U_\mathrm{P}^r + \sum A_\mathrm{N} \mathbf U_\mathrm{N}^n - E_\mathrm{P}^n=-\nabla p \end{equation} \tag{3}

可简写为

A U = B

该行程序双等号右边的 -fvc::grad(p) 是动量方程右边的部分，而使用 fvc:: 表示利用前一时间步的值来获得此梯度。

solve(UEqn == -fvc::grad(p)) 这一行就是求解 U 方程获取预测速度，但不能满足连续性方程。因此，需要进行修正。

如果我们定义 \mathbf U_\mathrm{P}^{n+1} 为 PISO 循环后最终收敛的速度，那么，进行修正后的 \mathbf U_\mathrm{P}^{n+1} 就应当满足动量方程和连续性方程。

速度 \mathbf U_\mathrm{P}^{n+1} 符合动量方程，则有

\begin{equation} A_\mathrm{P} \mathbf U_\mathrm{P}^{n+1} + \sum A_\mathrm{N}\mathbf U_\mathrm{N}^{n+1} - E_\mathrm{P}^{n+1} = - \nabla p^{n+1} \end{equation} \tag{4}

其中， p^{n+1} 为 PISO 循环内最后求得的压力。

\begin{equation} \mathbf U_\mathrm{P}^{n+1} = \mathbf{U}^{*,n+1} - \frac{1}{A_\mathrm{P}}\nabla p^{n+1} \end{equation} \tag{5}

式中， \mathbf{U}^{*,n+1} 为预测速度值。速度满足连续性方程，只需要对 \mathbf U_\mathrm{P}^{n+1} 进行求散度操作。

\begin{equation} \nabla \cdot \mathbf U_\mathrm{P}^{n+1}=0 \end{equation}\tag{6}

离散后得

\begin{equation} \sum \left(\mathbf U_\mathrm{P}^{n+1} \cdot \mathbf S_f \right)=0 \end{equation}\tag{7}

而此时，该网格面上的中心值也未求得，因此，需要另辟他法，Rhie-Chow 插值建议这样进行获取面上的 \mathbf U_{\mathrm{P,}f}^{n+1} ：

\begin{equation} \mathbf U_{\mathrm{P},f}^{n+1} = \mathbf{U}^{*,n+1}_f - \frac{1}{A_{\mathrm{P},f}} \nabla_f p^{n+1} \end{equation}\tag{8}

将获取的面心上的速度 \mathbf U_{\mathrm{P,}f}^{n+1} 代入连续方程，

\begin{equation} \sum \left[ \left( \mathbf{U}^{*,n+1}_f - \frac{1}{A_{\mathrm{P},f}}\nabla_f p^{n+1} \right) \cdot \mathbf S_f \right] = 0 \end{equation}\tag{9}

写成张量形式，有

\begin{equation} \nabla \cdot (\mathbf{U}^{*,n+1}) = \nabla \cdot\left(\frac{1}{A_{\mathrm{P},f}} \nabla p^{n+1}\right) \end{equation} \tag{10}

PISO 循环

方程的解 \mathbf U_\mathrm{P}^{n+1} 和压力 p^{n+1} 应当满足连续性方程和动量方程。而到目前为止，仅有通过动量预测方程求出来的 \mathbf U_\mathrm{P}^r 。动量预测方程和替代速度连续性方程无法自行求解。Issa 于 1986 年提出 PISO (Pressure Implicit with Splitting of Operators) 技术来解决这个问题，通过对当前时间步内的多次修正来获得最终的 \mathbf U_\mathrm{P}^{n+1} 和 p^{n+1} 。

icoFoam 使用 PISO 算法解 N-S 方程。其基本方法为

根据初始条件（压力场、速度场等）求解预测速度场 \mathbf U^r
根据预测速度 \mathbf U^r 或第 4 步求得的速度 \mathbf U^{n+1} ，计算出无压力梯度项的速度场 \mathbf U^n
根据无压力梯度项的速度场 \mathbf U^n ，求解压力（泊松）方程 (9) 或 (10)，得到压力场 p^n
根据压力场 p^n ，使用方程 (8) 修正预测速度 \mathbf U^{n+1} 以满足连续性方程，获得一个新的速度 \mathbf U^{n+1}
接第 2 步

Jackslow · 2019 年3 月 28 日 03:36

有个疑问，根据楼主对（2）的解读，

LHS的后两项都是已知速度值，结合（3）的解读“

“为什么同是已知值，一个要用fvc一个要用fvm?

rheoworks · 2019 年3 月 30 日 00:24

非常感谢 @Jackslow 读得这么仔细。

你的意思是不是说，在 UEqn 的离散中都采用了 fvm::，而在式 (3) 中的 \nabla p 又采用了 fvc::，两个似乎不能对应。

是这样的，UEqn 代表式 (1) 的左边三项，而式 (2) 是式是经过有限体积离散后写成的矩阵方程，中间做了处理，写成了待求网格 (第 1 项)，相邻网格 (第 2 项)，和其他值 (第 3 项)，各项与式 (1) 各项没有直接的对应关系。式 (1) LHS 各项都是待求速度。

PISO 的思想是求得一个预测速度，通过修正获得满足连续性方程的准确速度值。

zhaxi666 · 2019 年4 月 2 日 07:16

写得太棒了。
\rho

dengbin · 2019 年4 月 27 日 19:09

前辈，请问按照Rhie-Chow的插值建议，我们最后求得的岂不是最后是面心速度值？而不是体心速度值了？

chnrdu · 2019 年4 月 28 日 00:13

利用 Rhie-Chow 插值获得未知量的面心表达式，再代回离散方程，求体心的值。

dengbin · 2019 年4 月 30 日 08:31

前辈，对于面心表达式，回代离散方程的这个过程，可以给更多一点说明吗？谢谢

chnrdu · 2019 年5 月 1 日 00:46

首先，方程 (7) 是连续性方程的离散形式，不知道面心值，因此，需要使用 Rhie-Chow 插值方法 (8)。将这个方程代回方程 (7)，获得方程(9) ，方程 (9) 就是利用 Rhie-Chow 插值获得的连续性方程的离散形式。由于式 (9) 是离散形式，需要将其写成向量形式 (10)，从而可用 OpenFOAM 的方式来写方程。式 (10) 也就是通常所说的压力 (泊松) 方程，这个方程利用无压力项的动量方程 (2) 求得的速度 \boldsymbol{U}_p^* ，来求解压力。最后利用式 (8) 修正，获得面上的通量， Code of the Week #18: phi -= pEqn.flux()。而后利用方程 (8) 最后得到满足连续性方程的体心值 Code of the Week #19: fvc::reconstruct。

piteqiu · 2019 年7 月 23 日 05:44

还有个问题：式8,9,10中需要用到 ,此时只有单元中心的已知，还是不知道面上的值。面上的值是用单元值线性插值得到的吗？

chnrdu · 2019 年7 月 23 日 07:11

是否使用线性插值，根据 system/fvSchemes 里面的

interpolationSchemes
{
    interpolate(U)    linear;  // 确定 U 的面值的插值方式
}

来确定。

其他的插值方式查询，可以将 interpolate(U) linear 中的 linear 换成 banana，求解器遇到不能解释的，就会提醒哪些可用。