欧拉及拉格朗日下的动量方程压力项竟然不一样

[李东岳]

算例一直发散，整不明白。今天看了半天代码，跟方程对不上，发现差了一个相分数。很奇怪。于是翻了几个文献，发现

%(#ff0000)[Eulerian-Eulerian和Eulerian-Lagrangian下的动量方程压力项竟然不一样]

Eulerian-Lagrangian压力项不需要乘以相分数

主要看压力项，取自三篇不同的SCI，没有相分数！才发现

Eulerian-Eulerian压力项需要乘以相分数

.

[李东岳]

上面应该是某文章里面有个笔误的Eulerian-Lagraigian方程，应该没有相分数。
有误，见下文。

植入代码的时候要小心，如果方程有笔误的话，怎么可能收敛呢？不过我自己的文章也有很多笔误，没办法，方程太多..

[李东岳]

昨天在植入算法的时候（参考多相流的数学模型里面的方程），发现最后推导出来的压力泊松方程里面带一个${\alpha }_{\mathrm{c}}^2$ 。之前从来没见过带平方的压力泊松方程于是回头仔细仔细再仔细的看方程，看是不是哪里错了。

到头来，发现最开始的连续相动量方程文献里面竟然不一样，就是上文说的，压力梯度项相差一个相分数${\alpha }_{\mathrm{c}}^2$。此处不赘述。昨晚上搞完之后就出去喝酒了。现整理思路，总结一下。

欧拉欧拉以及欧拉拉格朗日下连续相方程的压力梯度项区别

现在把俩种形式的方程写一下，欧拉欧拉下连续相方程A：

$$\begin{array}{}\text{(1)}& \frac{\mathrm{\partial }\left({\alpha }_{\mathrm{c}}{\mathbf{U}}_{\mathrm{c}}\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_{\mathrm{c}}\left({\mathbf{U}}_{\mathrm{c}}\otimes {\mathbf{U}}_{\mathrm{c}}\right)\right)-\mathrm{\nabla }\cdot \left({\alpha }_{\mathrm{c}}{\mathbf{R}}_{\mathrm{c}}\right)=-{\alpha }_{\mathrm{c}}\mathrm{\nabla }\frac{p_{\mathrm{c}}}{{\rho }_{\mathrm{c}}}+{\alpha }_{\mathrm{c}}\mathbf{g}-{\mathbf{M}}_{\mathrm{d}}.\end{array}$$

欧拉拉格朗日下连续相方程B：

$$\begin{array}{}\text{(2)}& \frac{\mathrm{\partial }\left({\alpha }_{\mathrm{c}}{\mathbf{U}}_{\mathrm{c}}\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_{\mathrm{c}}\left({\mathbf{U}}_{\mathrm{c}}\otimes {\mathbf{U}}_{\mathrm{c}}\right)\right)-\mathrm{\nabla }\cdot \left({\alpha }_{\mathrm{c}}{\mathbf{R}}_{\mathrm{c}}\right)=-\mathrm{\nabla }\frac{p_{\mathrm{c}}}{{\rho }_{\mathrm{c}}}+{\alpha }_{\mathrm{c}}\mathbf{g}-{\mathbf{M}}_{\mathrm{d}}.\end{array}$$

方程$\text{(1)}$和$\text{(2)}$的区别在于一个$\alpha$。文中符号含义请参考多相流的数学模型，此略。

这个区别大么？

目前并不清楚。但很明确的是，从方程$\text{(1)}$最后推导出来的压力泊松方程为：
$$\begin{array}{}\text{(3)}& \frac{\mathrm{\partial }{\alpha }_{\mathrm{c}}}{\mathrm{\partial }t}+{\mathrm{\nabla }}_{\mathbf{x}}\cdot \left({\alpha }_{\mathrm{c}}(\mathbf{P}\mathbf{r}\mathbf{e}+\mathbf{H}\mathbf{b}\mathbf{y}\mathbf{A})\right)={\mathrm{\nabla }}_{\mathbf{x}}\cdot \left(\frac{{\alpha }_{\mathrm{c}}^2}{A_{\mathrm{P}}}\mathrm{\nabla }\frac{p_{\mathrm{c}}}{{\rho }_{\mathrm{c}}}\right).\end{array}$$
从方程$\text{(2)}$最后推导出来的压力泊松方程为：

$$\begin{array}{}\text{(4)}& \frac{\mathrm{\partial }{\alpha }_{\mathrm{c}}}{\mathrm{\partial }t}+{\mathrm{\nabla }}_{\mathbf{x}}\cdot \left({\alpha }_{\mathrm{c}}(\mathbf{P}\mathbf{r}\mathbf{e}+\mathbf{H}\mathbf{b}\mathbf{y}\mathbf{A})\right)={\mathrm{\nabla }}_{\mathbf{x}}\cdot \left(\frac{{\alpha }_{\mathrm{c}}}{A_{\mathrm{P}}}\mathrm{\nabla }\frac{p_{\mathrm{c}}}{{\rho }_{\mathrm{c}}}\right).\end{array}$$

在写代码的时候，分别对应俩种形式：

fvScalarMatrix pEqn
        (
            fvm::laplacian(alphacf*rAUcf, p)
         ==
            fvc::ddt(alphac) + fvc::div(alphacf*phiHbyA)
        );

fvScalarMatrix pEqn
        (
            fvm::laplacian(alphacf的平方*rAUcf, p)
         ==
            fvc::ddt(alphac) + fvc::div(alphacf*phiHbyA)
        );

为什么？

并不知晓。昨天把本帖随手发了个朋友圈，下面这个朋友的评论非常在理

自古高手在朋友圈啊 but why? 接下来，只能从物理的角度来推方程了。

首先，我们要考虑颗粒的大小：

• 大颗粒：颗粒体积相对网格单元体积较大，或者数量众多，导致不能忽略，即${\alpha }_{\mathrm{d}}\ne 0$
• 小颗粒：反之，${\alpha }_{\mathrm{d}}\approx 0$

然后，除了颗粒的受力，除了曳力（阻力）外，我们主要考虑颗粒所受的压力梯度，如下图中的$\mathrm{\nabla }p$。

最后，我们的压力梯度交换项就是图中的${\mathbf{M}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}}$。${\mathbf{M}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}}$的值就是导致带相分数或者不带相分数的区别。

考虑大颗粒小颗粒俩种情况，如下图：

• 对于小颗粒，${\alpha }_{\mathrm{d}}\approx 0$导致${\mathbf{M}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}}\approx 0$
• 对于大颗粒或者颗粒群，${\mathbf{M}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}}\approx -{\alpha }_{\mathrm{d}}\mathrm{\nabla }p$（图中忘记了一个负号）

好了，重新考虑方程$\text{(2)}$中右边的压力项有（${\rho }_{\mathrm{c}}=1$）：
$$\begin{array}{}\text{(5)}& -\mathrm{\nabla }{p}_{\mathrm{c}}-{\mathbf{M}}_{\mathrm{d}}=-\mathrm{\nabla }{p}_{\mathrm{c}}+{\alpha }_{\mathrm{d}}\mathrm{\nabla }{p}_{\mathrm{c}}=-{\alpha }_{\mathrm{d}}\mathrm{\nabla }{p}_{\mathrm{c}}\end{array}$$
这样，方程$\text{(1)}$和方程eqref{B}就相符了。

结论

方程$\text{(1)}$和方程$\text{(2)}$的区别就在于考没考虑粒子的压力梯度项：

• %(#ff0000)[考虑粒子的压力梯度，则有相分数，对应方程]$\text{(1)}$
• %(#ff0000)[没考虑粒子的压力梯度，则没有相分数，对应方程]$\text{(2)}$

我在好多论文里面看到直接给出方程$\text{(1)}$和方程$\text{(2)}$，没有一点解释。是太过基础了么？总之，现在我明白了。

（虽然英文图，但我自己画的）

[李东岳]

找到了个文章，有一点讨论：

The multiphase particle-in-cell (MP-PIC) method for dense particulate flows

1996年的MPPIC算法鼻祖文章，单篇被引300+，感兴趣的可以看看。

[linhan.ge]

这篇参考文献有相关讨论：
Zhou, Z.Y., Kuang, S.B., Chu, K.W., Yu, A.B., 2010. Discrete particle simulation of particle–fluid flow: model formulations and their applicability. Journal of Fluid Mechanics 661, 482-510.

[linhan.ge]

@东岳 方程(3)和(4)是不是颠倒了？(1)推出(4)，(2)推出(3)。

[李东岳]

@linhan-ge 是的，已更正，目前正确了。

[李东岳]

PRESSURE GRADIENT TERM这个图中的$F=\frac{m}{V}\mathrm{\nabla }p$，其中的$V$应该为$\rho$