twoPhaseEulerFoam中的dilatation rates

dzw05

最近在推导体积分数alpha的方程。理论上来说compressibleInterFoam和twoPhaseEulerFoam中的体积分数方程应该是相同的（不考虑相间的质量传递时），即“东岳流体”中的公式：

在compressibleInterFoam中，dilatation项定义为：

dgdt =
            (
                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
            );

这个和东岳流体中的定义是一致的。
而在twoPhaseEulerFoam中，dilatation项的定义为：

            fluid.dgdt() =
            (
                alpha1*(pEqnComp2 & p_rgh)
              - alpha2*(pEqnComp1 & p_rgh)
            );

这个显然和公式推导出的不一致，但是在后续的体积分数方程计算时，我并没有看到程序对这项有什么特殊处理，这是否是错误的？
另外还想问问各位大神，pEqnComp2 & p_rgh是代表什么意思？pEqnComp2 是fvScalarMatrix类型，而p_rgh是volScalarField类型，中间的操作符&表示什么，是OpenFoam进行算符重构了吗？

李东岳

dgdt要和MULES求解的相方程匹配。单看twoPhaseEulerFoam，在OpenFOAM-2.2.x中定义为：

dgdt =
            (
                pos(alpha2)*(pEqnComp2 & p)/rho2
              - pos(alpha1)*(pEqnComp1 & p)/rho1
            );

在OpenFOAM-2.3.x中定义为：

fluid.dgdt() =
            (
                alpha1*(pEqnComp2 & p_rgh)
              - alpha2*(pEqnComp1 & p_rgh)
            );

我并没有看到程序对这项有什么特殊处理

你应该看的是OpenFOAM-2.3.x以后的版本。你可以看一下alphaEqn的处理。

&是进行重构了。矩阵稀疏×场，在OpenFOAM中很常见，类似一种反离散 ，还原回去。

dzw05

@李东岳多谢东岳。还有几处不太明白。
（1）OpenFoam对 &的重构是在哪定义的呢？具体到当前的这个&，类似于反离散？具体返回什么值呢？
（2）我看的是OpenFOAM-4.0.x。无论dgdt的定义如何，最后应该都能回归到体积分数的求解公式，也就是：
OpenFOAM-2.3中twoPhaseEulerFoam里dgdt的定义，与OpenFOAM4.0.x里compressibleInterFoam的定义是相同的，最后在显式和隐式源项相加后还能恢复到公式。
但是OpenFOAM4.0.x里twoPhaseEulerFoam的定义中，显式和隐式源项相加后并不能得到体积分数公式。在twoPhaseSystem中，其源项定义为：

 forAll(dgdt_, celli)
        {
            if (dgdt_[celli] > 0.0)
            {
                Sp[celli] -= dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
                Su[celli] += dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
            }
            else if (dgdt_[celli] < 0.0)
            {
                Sp[celli] += dgdt_[celli]/max(alpha1[celli], 1e-4);
            }
        }

也就是得到，alpha*Sp+Su=dgdt + Su这个结果，这个公式显然与体积分数公式不同。体积分数公式中右端项为
而现在得到的结果是：。

李东岳

是的，你说的对。我又看了下，在23x以后，dgdt中pEqnComp的定义为：

pEqnComp1 =
            (
                fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
              - (fvOptions(alpha1, rho1)&rho1)
              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
            )/rho1
          + (alpha1*psi1/rho1)*correction(fvm::ddt(p));

        pEqnComp2 =
            (
                fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
              - (fvOptions(alpha2, rho2)&rho2)
              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
            )/rho2
          + (alpha2*psi2/rho2)*correction(fvm::ddt(p));

在这俩个定义中，和分别包含了alpha，因此再乘以时候，就是alpha1*alpha2。

dzw05

@李东岳多谢多谢，是我自己看程序不仔细，pEqnComp1确实已经包含了alpha1。

MrZZQi

请问 pEqnComp1 & p_rgh 这一项为什么是$ D\rho /Dt$ ？而且在新旧两个版本的pEqncomp1表达式中都是这样的。我搜索到了一些资料，https://www.cfd-online.com/Forums/openfoam/92947-openfoam-syntax.html https://www.cfd-online.com/Forums/openfoam-solving/70965-formulation-compressibleinterfoam-4.html 但都没有解释好这个问题。请问有可以参考的资料吗？谢谢!

李东岳

@MrZZQi http://dyfluid.com/reactingTwoPhaseEulerFoam.html

MrZZQi

@东岳感谢回复，我没说清楚，我对方程本身没有疑问，只是对之前讨论中您提到的

&是进行重构了。矩阵稀疏×场，在OpenFOAM中很常见，类似一种反离散

这句话不是很理解，所以不知道为什么pEqnComp1 & p_rgh这段代码对应的是$D \rho /D t$。请问有关于这种反离散的资料吗？

李东岳

你知道对$\int u\rd V$离散后是$[a u]$, 对整个网格场进行离散就是一个列向量$[a_0 u_0,a_1u_1,...,a_nu_n]^T$，其中$a$就是离散稀疏（行向量），$u$就是求解的量（列向量）。所以如果给定$a$，再给定$u$，进行车中相乘，就是$[a_0 u_0,a_1u_1,...,a_nu_n]^T$，类似一种反离散。没什么资料

MrZZQi

@东岳非常感谢！那我试着推推看这里

CFD中文网

twoPhaseEulerFoam中的dilatation rates