reactingTwoPhaseEulerFoam求解器的相变求解过程

[zwl]

从reactingTwoPhaseEulerFoam求解器的介绍和给出的算例可以看出，该可压缩求解器可通过组分输运方程考虑两种组分以上的多相流问题。通过看该求解器代码，我发现质量源项只出现在了组分输运方程中，而相方程中没有源项，不论是蒸发还是沸腾液相和气相都有质量交换，那么该求解器是通过什么方式考虑相间质量传递的呢？如有熟悉该求解的高手，闲暇之余，还请指点一二，捅破这层窗户纸~,先行谢过！

[李东岳]

相方程直接对应压力方程，因此在压力房产里面

[zwl]

@东岳 非常感谢东岳老师的指点，去看了压力方程后确实发现在压力方程里考虑了质量传递，但我有一点不明白，相方程是先求解的，相方程里不给质量源相，那体积分数分布怎么保证是对的？对不起，也许这是简单的问题，但我对方程的求解细节不太清楚，还请解答一下。谢谢。

[李东岳]

Here:

if (phase1_.divU().valid() && phase2_.divU().valid())
    {
        tdgdt =
        (
            alpha2()
           *phase1_.divU()()()
          - alpha1()
           *phase2_.divU()()()
        );
    }

// Set the phase dilatation rates
            if (pEqnComp1.valid())
            {
                phase1.divU(-pEqnComp1 & p_rgh);
            }
            if (pEqnComp2.valid())
            {
                phase2.divU(-pEqnComp2 & p_rgh);
            }

[zwl]

@东岳 非常感谢，之前一直以为这些值都是在类之间传递，赋值大部分都是在类内部完成，却忽略了压力方程内进行赋值进而用来求解相方程，技巧性确实不错。

[zwl]

@东岳 今天仔细推导了可压缩两相方程，基本都能与代码都能对的上，pEqnComp代表的就是可压缩项与质量源相的和，但是为什么要返给相方程pEqnComp与p_rgh的积？这跟您之前分析的其它可压缩求解器也不一样？能分析下为什么这么处理吗？公式中没有质量源相，相方程可压缩部分与压力的关系，这样处理的目的和原理是什么？谢谢。

[李东岳]

为什么要返给相方程pEqnComp与p_rgh的积？

pEqnComp是矩阵系数，与p_rgh的积才是真正的volScalarField

[zwl]

@东岳 那么从物理角度应该怎么理解呢？dgdt*p_rgh什么物理意义？应该不只是体标量场这么简单吧。非常感谢东岳老师耐心的回复。

[李东岳]

那个是dgdt，不是dgdt*p_rgh哟

[zwl]

@东岳 对，您说的是对的。那我这样表达下我的意思吧：
reactingTwoPhaseEulerFaom的dgdt这里叫dgdt1；compressibleInterFoam里的dgdt这里叫dgdt2；我根据代码解读看到dgdt1=alpha1* alpha2* dgdt2* p_rgh，而且我感觉dgdt1应该定义成alpha1* alpha2* dgdt2才对，膨胀率乘以压力的处理和意义是公式里没有的，所以我才一直有疑惑。不知道我说明白了没有。

[zwl]

@zwl 在 reactingTwoPhaseEulerFoam求解器的相变求解过程 中说：

dgdt1=alpha1alpha2dgdt2p_rgh

乘号漏掉了：dgdt1=alpha1 * alpha2 * dgdt2*p_rgh

[李东岳]

我根据代码解读看到dgdt1=alpha1* alpha2* dgdt2* p_rgh

是dgdt1=alpha1* alpha2* dgdt2吧，就像我上面说的那个

那个是dgdt，不是dgdt*p_rgh哟

膨胀率乘以压力的处理和意义是公式里没有的

是没有啊，那个是dgdt，不是dgdt乘以p_rgh..

[zwl]

@东岳 代码如下：

    tmp<fvScalarMatrix> pEqnComp1;
    tmp<fvScalarMatrix> pEqnComp2;

    // Construct the compressibility parts of the pressure equation
    if (pimple.transonic())
    {
        if (phase1.compressible())
        {
            surfaceScalarField phid1
            (
                IOobject::groupName("phid", phase1.name()),
                fvc::interpolate(psi1)*phi1
            );

            pEqnComp1 =
                (
                    phase1.continuityError()
                  - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
                )/rho1           //相方程可压缩部分 
              + correction
                (
                    (alpha1/rho1)*
                    (
                        psi1*fvm::ddt(p_rgh)
                      + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
                    )
                );

            deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
            pEqnComp1.ref().relax();
        }

        if (phase1.compressible())
        {
            phase1.divU(-pEqnComp1 & p_rgh);//divU(divU_)函数返回自身参数
        }

 
    // Construct the dilatation rate source term
    tmp<volScalarField::Internal> tdgdt;
 
    if (phase1_.divU().valid() && phase2_.divU().valid())
    {
        tdgdt =
        (
            alpha2()
           *phase1_.divU()()()
          - alpha1()
           *phase2_.divU()()()
        );
    }

如此看来是不是dgdt1=alpha1* alpha2* dgdt2* p_rgh？难道我还有遗漏？

[李东岳]

没理解，就像我之前说的

pEqnComp是矩阵系数，与p_rgh的积才是真正的volScalarField

代码中的dgdt2* p_rgh实际上是公式中的dgdt2，pEqnComp1是矩阵系数并不是值，与p_rgh的积才是dgdt2

[zwl]

@东岳 懂了，谢谢。