OpenFOAM-6的twoPhaseEulerFoam学习

[Samuel-Tu]

最近结合东岳老师的论文，初步学习了下OpenFOAM-6中twoPhaseEulerFoam源代码中的动量方程（依据pU/UEqns.H），得到了一些收获，也有很多疑惑。现在分享给大家，也希望大家能帮忙解惑。

pU/UEqns.H中的动量方程形式是东岳老师论文中提到的全守恒形式，因为密度没有被除，见下图1。 对应的代码及相关注释如下。其中，g、压力梯度项都放入了压力修正里面去了，而虚拟质量力和拽力保留了下来。

 U1Eqn =
        (
            fvm::ddt(alpha1, rho1, U1)//图1等号左边第一项  
          + fvm::div(alphaRhoPhi1, U1)//图1等号左边第二项  
          - fvm::Sp(contErr1, U1)//图1等号左边第三项，查阅coutErr.H后发现多了一项fvOption
          + MRF.DDt(alpha1*rho1 + Vm, U1)//不知道这项对应什么
          + phase1.turbulence().divDevRhoReff(U1)//可能对应图1等号左边第三项和第四项??
         ==
          - Vm
           *(
                fvm::ddt(U1)
              + fvm::div(phi1, U1)
              - fvm::Sp(fvc::div(phi1), U1)
              - DDtU2
            )//应该是虚拟质量力，论文中忽略了这一项，和东岳流体中的虚拟质量力形式吻合
          + fvOptions(alpha1, rho1, U1)//fvOption，不知道对应什么
        );
        U1Eqn.relax();
        U1Eqn += fvm::Sp(Kd, U1);//对应拽力

如果将代码中的MRF.DDT项、fvOPtion项和虚拟质量项全部去掉，则与图1中的全守恒公式完全对应。但对于divDevRhoReff（U1）项我通过Doxygen查看了相关代码（不知道是否正确），如下：

- fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))
       - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)

其中laplacian项对应与图1全守恒公式左边第五项，但nuEff()*dev2项似乎应该等于 ，与论文中的形式 差了 这部分差的值放到那里去了呢？
未完待续。。。。

[东岳]

是的，后面的kI放在了压力里面

[Samuel-Tu]

@东岳 感谢解惑(｡･ω･｡)ﾉ♡

[Samuel-Tu]

继续更新。。。
组成UEqn之后，#include EEqns.H，看了下内容是和thermo有关的，怀疑是组成能量方程的，由于我的方向不涉及热交换，所以忽略了这一步。
继续接着看 #include pU/pEqn.H ，由于相1和相2的方程类似，因此只以相1为例：
首先定语rAU1为矩阵主对角线的倒数，即1/A，代码如下：

1.0 /(
        U1Eqn.A()
      + max(phase1.residualAlpha() - alpha1, scalar(0))
       *rho1/runTime.deltaT()

但不知道max(phase1.residualAlpha()-alpha1 , scalar(0)) * rho2/runTime.deltaT()的作用是什么。
不过分析来看，如果alpha1大于residualAlpha()的话，这一项就没有作用， 那么就rAU1=1/A，刚好就是符合我的认知。而如果alpha1小于residualAlpha()的话，就多了一项。
我怀疑这是为了防止alpha1过小的一种数值处理，但是不知到原理是什么。希望大家能够帮忙解惑

[Samuel-Tu]

接着上面，发现

surfaceScalarField alpharAUf1
(
    fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1)
);

alpharAUf1=alpha1/A1，这是从体心插值到了面心，后面泊松方程用RC插值到面心的时候需要用。另外这里也比较了alpha1和phase1.residualAlpha()的大小。如果alpha1小于phase1.residualAlpha()，则取phase1.residualAlpha()。这似乎是东岳老师论文中提到了防止相体积分数过小的一种处理方法1：直接指定一个小值。
接着定义了HbyA1=H1/A1，代码如下：

HbyA1 =
        rAU1
       *(
            U1Eqn.H()
          + max(phase1.residualAlpha() - alpha1, scalar(0))
           *rho1*U1.oldTime()/runTime.deltaT()
        );

这一项是非对角线元素除以主对角线元素，符合预期，但同时这里也用了一个alpha1和residualAlpha()的比较，可能也是防止alpha1过小的数值处理。
接着定义了ghSnGradRho:

    surfaceScalarField ghSnGradRho
    (
        "ghSnGradRho",
        ghf*fvc::snGrad(rho)*mesh.magSf()
    );

查阅相关ghf源码，发现：

     Info<< "Calculating field g.h\n" << endl;
     dimensionedScalar ghRef(- mag(g)*hRef);
     volScalarField gh("gh", (g & mesh.C()) - ghRef);
     surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);

说明ghSnGradRho指

此项在泊松方程中的RC插值到面上有用。
接下来组合成了和g有关的泊松方程中的项phig1：

    surfaceScalarField phig1
    (
        alpharAUf1
       *(
           ghSnGradRho
         - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
        )
    );

对应于
不过符号我还没有检查，初步推导的时候感觉似乎符号有点问题。

[Samuel-Tu]

接下来定义了HbyA1、g、湍流、壁面润滑力和升力造成的通量phiHbyA1:

    surfaceScalarField phiHbyA1
    (
        IOobject::groupName("phiHbyA", phase1.name()),
        fvc::flux(HbyA1)
      + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
       *(
            MRF.absolute(phi1.oldTime())
          - fvc::flux(U1.oldTime())
        )/runTime.deltaT()
      - phiF1()
      - phig1
    );

对应公式为 phiF1对应升力、湍流和壁面润滑力的作用，phig1上条评论已经分析了。代码中还多了一个修正项phiCorrCoeff1项，不知道这项是干嘛的，希望大家解答。。
注意到上述代码还缺少了显式的曳力项，因此在组成最后的压力泊松方程等号左边的项时，作如下处理：

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        alphaf1*(phiHbyA1 + rAUKd1*MRF.absolute(phi2))
      + alphaf2*(phiHbyA2 + rAUKd2*MRF.absolute(phi1))
    );
    MRF.makeRelative(phiHbyA);

对应公式
其中rAUKd1为

    surfaceScalarField rAUKd1(fvc::interpolate(rAU1*Kd));

[Samuel-Tu]

@东岳 在 OpenFOAM-6的twoPhaseEulerFoam学习 中说：

是的，后面的kI放在了压力里面

最开始我以为是把kI放入了压力泊松方程中，但是我pEqn里一直都没有出现kI，后来我在xiaopingqiu的个人网站上看见了一种说法：动量方程中的压力是雷诺时均压力和雷诺应力的各向同性分量（2k/3)之和。
这种观点可以在 Pope 2000 书第 88 页找到依据。在 “ The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM® and Matlab®” 这本书的第 699 页，也提到 k 是被放到压力项里去了，目的在于使动量方程中只含有 νt 这一个跟湍流有关的未知量——来源于xiaopingqiu

[Dingcy]

我有一个疑问，在源代码中pUf/UEqns.H代表了什么呢？因为我在twoPhaseEulerFoam.C文件中看见了如下代码

if (faceMomentum)
            {
                #include "pUf/UEqns.H"
                #include "EEqns.H"
                #include "pUf/pEqn.H"
                #include "pUf/DDtU.H"
            }
            else
            {
                #include "pU/UEqns.H"
                #include "EEqns.H"
                #include "pU/pEqn.H"
                #include "pU/DDtU.H"
            }

很明显，在计算过程中进行了选择，代码中的faceMomentum代表了什么？我在算例文件中的设置中没发现进行选择的选项？

[Zhy2022]

@Samuel-Tu kI项是雷诺应力项对照柯西应力项中正应力部分写出来的，而p正好是柯西应力中的正应力部分，所以两者放在一起