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

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

今天,我们要讲的是

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 方程。其基本方法为

  1. 根据初始条件(压力场、速度场等)求解预测速度场 \mathbf U^r

  2. 根据预测速度 \mathbf U^r 或第 4 步求得的速度 \mathbf U^{n+1} ,计算出无压力梯度项的速度场 \mathbf U^n

  3. 根据无压力梯度项的速度场 \mathbf U^n ,求解压力(泊松)方程 (9) 或 (10),得到压力场 p^n

  4. 根据压力场 p^n ,使用方程 (8) 修正预测速度 \mathbf U^{n+1} 以满足连续性方程,获得一个新的速度 \mathbf U^{n+1}

  5. 接第 2 步

2 Likes

有个疑问,根据楼主对(2)的解读,

LHS的后两项都是已知速度值,结合(3)的解读“

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

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

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

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

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

1 Like

写得太棒了。
\rho

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

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

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

首先,方程 (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

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

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

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

来确定。

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