MPPICInterFoam 中如何通过alphac得到alpha1和alpha2?

Jasper 0

大家好，本人最近在研究MPPICInterFoam中的代码，在createFields中，有如下代码。对alphac进行定义之后，最后两行通过 alphac=1-theta 对alphac进行了更新，但是我好奇alpha1和alpha2是如何更新的呢？请大佬帮忙解答一下，谢谢。

 volScalarField alphac
 (
     IOobject
     (
         "alphac",
         runTime.timeName(),
         mesh,
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
     dimensionedScalar(dimless, Zero),
     zeroGradientFvPatchScalarField::typeName
 );
 alphac.oldTime();
 
 volScalarField alphacRho(alphac*rho);
 alphacRho.oldTime();
 
 Info<< "Constructing kinematicCloud " << endl;
 basicKinematicMPPICCloud kinematicCloud
 (
     "kinematicCloud",
     rho,
     U,
     mu,
     g
 );
 
 // Particle fraction upper limit
 scalar alphacMin
 (
     1.0
   - (
         kinematicCloud.particleProperties().subDict("constantProperties")
        .get<scalar>("alphaMax")
     )
 );
 
 // Update alphac from the particle locations
 alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
 alphac.correctBoundaryConditions();

星星星星晴

@jasper-0 在 MPPICInterFoam 中如何通过alphac得到alpha1和alpha2? 中说：

MPPICInterFoam

你这个看起来不是openfoam.org 原生的solver吧？
如果是下面的这个solver的话，alpha1 和 alpha2 会自动更新

https://github.com/TonkomoLLC/MPPICInterDyMFoam/blob/master/MPPICInterDyMFoam/createFields.H

volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());

Jasper 0

@星星星星晴哦我跟您的是一样的，我只是选取了一部分代码。您给出的这两行在alphac的定义上面，我认为作为代码层面，是不是不会返回去求alpha1和2了？

星星星星晴

@jasper-0 不好意思，这部分我不懂，我不做VOF的，只说是在createfield.H里面，alpha1和alpha2是这么定义的，而且会随着mixture.alpha1 和mixture.alpha2 的变化更新，至于具体怎么算的，你这个还要看你的solver 还有后续的过程，这只是初始化，肯定后续有计算和更新的

Jasper 0

@星星星星晴谢谢大佬了，我再继续研究一下

Tens

@jasper-0 我个人理解是这样的，createFields.H只是创建初始场，所以没必要这么准确，而在循环中的alpha1是在alphaEqn.H里直接求，alpha2=1-alpha1。

Jasper 0

@tens 谢谢您的回复，通过alphaEqn.H 发现有alphac和alpha1的关系

Tens

@jasper-0

Jasper 0

@tens 在 MPPICInterFoam 中如何通过alphac得到alpha1和alpha2? 中说：

@jasper-0

通过alphaEqn.H中的206-226行有如下代码，这几行代码应该是求解出来了乘以alphac之后的alpha1，那alpha2怎么是1-alpha1呢？不是应该alphac-alpha1吗？还望高手解答一下啊

else
        {
            alphaPhi10 = talphaPhi1Un;

            MULES::explicitSolve
            (
                alphac,
                alpha1,
                phiCN,
                alphaPhi10,
                Sp,
                (Su + divU*min(alpha1(), scalar(1)))(),
                oneField(),
                zeroField()
            );
        }

        alpha2 = 1.0 - alpha1;

        mixture.correct();
    }

Tens

@jasper-0
我的理解，alphac是颗粒的体积分数，1-alphac是连续相的体积分数，alpha1是连续相中液相占连续相的体积分数，alpha2是连续相中气相占连续相的体积分数。

Tens

@tens 睡觉前发现说反了，alphac是连续相的，1-alphac是液滴的，总之我理解的是alpha1和alpha2是占连续相的而不是总的

Jasper 0

@tens 对，所以我好奇怎么在程序中alpha2 = 1.0 - alpha1？而不是alpha2 = alphac - alpha1？

zhouxu

@tens 膜拜大佬的作息

Tens

@jasper-0 难道讲的不好理解吗
因为alpha1和alpha2是占连续相的而不是总的，连续相中只有液相或者气相，所以除了液相(alpha1)就是气相(alpha2=1-alpha1)，这个1不是总体积，而是总的连续相。
换个意思说，如果全部是占总体积的分数，颗粒为1-alphac，连续液相为alpha1*alphac，连续气相为alpha2 * alphac

Jasper 0

@tens 感谢大佬的解答！也就是说通过以下代码求得的alpha1是占连续相的体积分数？那最后结果跑出来后要想得到总的alpha1是不是要在paraview里面做计算alpha1*alphac？

 MULES::explicitSolve

            (

                alphac,

                alpha1,

                phiCN,

                alphaPhi10,

                Sp,

                (Su + divU*min(alpha1(), scalar(1)))(),

                oneField(),

                zeroField()

            );

        }

Tens

@jasper-0 是的吧

Jasper 0

@tens 多谢大佬！

Jasper 0

@jasper-0 在 MPPICInterFoam 中如何通过alphac得到alpha1和alpha2? 中说：

@tens 感谢大佬的解答！也就是说通过以下代码求得的alpha1是占连续相的体积分数？那最后结果跑出来后要想得到总的alpha1是不是要在paraview里面做计算alpha1*alphac？
 MULES::explicitSolve

            (

                alphac,

                alpha1,

                phiCN,

                alphaPhi10,

                Sp,

                (Su + divU*min(alpha1(), scalar(1)))(),

                oneField(),

                zeroField()

            );

        }

我发现interFoam中alpha1（即alpha1）的求解对应下面这个公式，

而MPPICInterFoam中通过以上代码求得的alphan为

公式右边第二项是用alpha占连续相的体积分数算出来的，
而和都是总的相分数，

所以此方程右边第二项是不是应该乘以个alphac呢？即

还望大佬解答一下，万分感谢！

Tens

@jasper-0 这个公式是哪里看到的，字母对应的意义不理解，在程序中相体积分数应该是通过控制方程求的吧

fvScalarMatrix alpha1Eqn
        (
            (
                LTS
              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alphac, alpha1)
              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            )
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
         - fvm::Sp(fvc::ddt(alphac) + fvc::div(phiCN), alpha1)
         ==
            Su + fvm::Sp(Sp + divU, alpha1)
        );

        alpha1Eqn.solve();

Jasper 0

@tens 这一步应该是用upwind格式求出来alpha1之后然后下面再通过MULES::correct进行修正。如果不进行修正的话，就是我之前提到的MULES::explicitSolve，在MULESTemplates.C里面定义了这个函数：

  template<class RdeltaTType, class RhoType, class SpType, class SuType>
 
  void Foam::MULES::explicitSolve
 
  (
 
      const RdeltaTType& rDeltaT,
 
      const RhoType& rho,
 
      volScalarField& psi,
 
      const surfaceScalarField& phiPsi,
 
      const SpType& Sp,
 
      const SuType& Su
 
  )
 
  {
 
      Info<< "MULES: Solving for " << psi.name() << endl;
 
   
 
      const fvMesh& mesh = psi.mesh();
 
   
 
      scalarField& psiIf = psi;
 
      const scalarField& psi0 = psi.oldTime();
 
   
 
      psiIf = 0.0;
 
      fvc::surfaceIntegrate(psiIf, phiPsi);
 
   
 
      if (mesh.moving())
 
      {
 
          psiIf =
 
          (
 
              mesh.Vsc0()().field()*rho.oldTime().field()
 
             *psi0*rDeltaT/mesh.Vsc()().field()
 
            + Su.field()
 
            - psiIf
 
          )/(rho.field()*rDeltaT - Sp.field());
 
      }
 
      else
 
      {
 
          psiIf =
 
          (
 
              rho.oldTime().field()*psi0*rDeltaT
 
            + Su.field()
 
            - psiIf
 
          )/(rho.field()*rDeltaT - Sp.field());
 
      }
 
   
 
      psi.correctBoundaryConditions();
 
  }

其中

 psiIf =
     (
         rho.oldTime().field()*psi0*rDeltaT
       + Su.field()
       - psiIf
     )/(rho.field()*rDeltaT - Sp.field());

就是我上面提到的方程
alphan应该就是最终求得的新的时间步的体积分数

Tens

@jasper-0

懂了，现在就是说是红色部分的问题对吧，我看了代码中这部分是通过下面代码求的

fvc::surfaceIntegrate(psiIf, phiPsi)

而里面的phiPsi则对应的这一部分

通过追踪

alphaPhi10 = talphaPhi1Un;

//又查到

tmp<surfaceScalarField> talphaPhi1Un
        (
            fvc::flux
            (
                phiCN(),
                cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
                alphaScheme
            )
          + fvc::flux
            (
               -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );

而其中的phiCN则

    tmp<surfaceScalarField> phiCN(alphaPhic);

    // Calculate the Crank-Nicolson off-centred volumetric flux
    if (ocCoeff > 0)
    {
        phiCN = cnCoeff*alphaPhic + (1.0 - cnCoeff)*alphaPhic.oldTime();
    }

如果不采用CN格式phiCN就等于alphaPhic，又查到

 alphacf = fvc::interpolate(alphac);
alphaPhic = alphacf*phi;

所以本质上图一中红色部分还是通过alphac*alpha1得到的

Jasper 0

@tens 懂了懂了，大佬牛！

Tens

@jasper-0
互相交流，以前从没这么认真的看过这么底层的东西，这次也学到了很多

Jasper 0

@tens

zhe

@jasper-0 你好，看到你这个问题却是很有意思。解释的code与公式对应起来了。但我相问一下，那么最终为什么interFoam和MPPICInterFoam的公式还是不一样啊？

Jasper 0

@zhe MPPICInterFoam 就是多了个alphac啊，考虑拉格朗日粒子对体积分数的影响

zhe

@jasper-0 换句话说，alphac在interFoam里就是1，所以可以直接省略。谢谢。

Jasper 0

@zhe 对

zhe

我测试了一下大家说的，在MPPICInterFoam里，加入粒子。但是alpha.water和alphac是不符合逻辑的。在最底部分，应该是只有水，那么alpha的值就该都是0.38。为什么在alpha.water却是0.62？

Jasper 0

@zhe alpha.water 是占连续相的体积分数，如果没有空气只有水的话，理论上应该alpha.water都是1吧，你看看是不是你alpha.water图例显示的问题？

zhe

@jasper-0 理论上，是应该这样的。我有两个问题：第一，哪里提到alpha1和aloha2只是占连续相的比值？我们希望它只是占连续相比值，但code里并没有提过吧。我是怕先入为主了，因为确实对不上。第二，你对比的alpha1的求解公式，帮助很大。但从物理意义上来说，alpha1*alphac就是考虑了粒子影响后的求解方式？它俩相乘是啥意义啊？

Jasper 0

@tens 在 MPPICInterFoam 中如何通过alphac得到alpha1和alpha2? 中说：

@jasper-0

懂了，现在就是说是红色部分的问题对吧，我看了代码中这部分是通过下面代码求的

fvc::surfaceIntegrate(psiIf, phiPsi)

而里面的phiPsi则对应的这一部分

通过追踪

alphaPhi10 = talphaPhi1Un;

//又查到

tmp<surfaceScalarField> talphaPhi1Un
        (
            fvc::flux
            (
                phiCN(),
                cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
                alphaScheme
            )
          + fvc::flux
            (
               -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );

而其中的phiCN则

    tmp<surfaceScalarField> phiCN(alphaPhic);

    // Calculate the Crank-Nicolson off-centred volumetric flux
    if (ocCoeff > 0)
    {
        phiCN = cnCoeff*alphaPhic + (1.0 - cnCoeff)*alphaPhic.oldTime();
    }

如果不采用CN格式phiCN就等于alphaPhic，又查到

 alphacf = fvc::interpolate(alphac);
alphaPhic = alphacf*phi;

所以本质上图一中红色部分还是通过alphac*alpha1得到的

@zhe 看这条消息，code中就是按这个公式求的，alphan*alphac是真实的相体积分数，因此求得的alphan就是占连续相的体积分数

zhe

@jasper-0 不好意思，最近没看到回复。如果像你说的alpha1*alphac就是真实的体积分数，那么alpha1才应该是整个cell里的alpha1啊，否则不需要再去乘以alphac.
我之前怀疑过，MPPICInterFoam并不是像我们想象的那样，alpha1+alpha2+alpha_p = alphac + alpha_p = 1。而是按照interFoam的alpha1 + alpha2 = 1进行了两相流的计算。随后在动量方程和连续方程中，乘以alphac（这个值与lapha1和alpha2无关，而是用1-alpha_p来的）的效果，以保证真实流过的体积是考虑过粒子占比后的，也能解释方程右边加入的momentum transfer的力。所以可以说OpenFOAM目前做的三相流耦合，还是秉承了粒子在流体中不占体积来的（只是有相互作用力）。
以上如果有不对的，希望可以多交流，感谢！

Jasper 0

@zhe 在 MPPICInterFoam 中如何通过alphac得到alpha1和alpha2? 中说：

那么alpha1才应该是整个cell里的alpha1啊

没看懂这是什么意思。举个例子，比如alphac=1-theta=0.8, alpha1=0.5，那么phase1所占整个cell的体积分数就是0.8*0.5=0.4，phase2也是0.4，0.4+0.4+0.2=1，在MPPICInterFoam中粒子是考虑了体积分数的

zhe

@jasper-0 MPPICInterFoam中考虑了体积分数是为了保证质量守恒。你举的例子，alpha1=0.5本身就是占据整个cell的体积分数。而乘以alphac后得到的是实际占据cell的体积分数（因为cell中被粒子占据了一部分）。所以，整个过程中，alpha1一直都是0.5（alpha2=1-alpha1=1-0.5=0.5)这样是符合算法的。而上述例子中最后的0.4+0.4+0.2=1。可是OpenFOAM对于alpha2的求解是通过1-alpha1. 0.4+0.4不等于1。
考虑一个简单的例子，只有粒子和水。alphac=1-theta=0.8。那么如果考虑了粒子的占比，那么水在这里就该是alpha1=alphac=0.8.那么在根据alpha1alphac=0.64?这个值是什么？
然而回到我说的没有考虑粒子占比，对于alpha1的值，一直是占据cell的，那么alpha1就该是等于1（alpha2=1-alpha1=1-1=0. 符合假设前提）。而在实际运动过程中，由于粒子占比，alpha1实际的流通量是alpha1alphac=1*0.8=0.8. 这里考虑了alphac的影响，但是跟alpha1本身的值没有任何关系。
我做了个最简单的验证，上下一共5个网格，最下面的放了一颗粒子。下面3个网格是水，上面是空气。得到的数据是

希望能帮助大家。

Tens

@zhe alpha1一直都不是占cell的体积分数，只是占连续相的体积分数，不是alpha1+alpha2=alphac，而是再次将alphac看做单位1，alpha1+alpha2=1
对于你的例子，只有粒子和水，那就是alphac=0.8，alpha1=1

Tens

@zhe 对alpha1和alpha2来说，分母不是Vcell，而是Vc

zhe

@tens 你的意思我明白，但是将alphac看做单位1是混淆了概念。alphac就是alphac,它是个值，并是有时不是1。这个操作在OpenFOAM中没有。在概念上来说，你说的通，我是认同的。但是在存在MPPICInterFoam中，这种把alphac当作1,是不存在的。alpha1一直以来都是占cell的体积分数，你所说的把alphac当作1，也不是不对，而是当作1，其实就是cell本身。

这里有人会迷糊，说那么alphac到底干啥用？它就是用来体现粒子占据一定体积后，对流体总体积造成了影响，不是来梳理alphac, alpha1, alpha2的关系的。

zhe

@tens 之前你说不是alpha1+alpha2=alphac. 而这个图里又表示V1+V2=Vc?

Tens

@zhe 体积和占比并不是一个概念

zhe

@tens 不好意思，我没看到。V1和V2是体积，因为咱们一直在说的是体积分数的问题。你的意思我是懂的，是很正统的三相流耦合。但在OpenFOAM的MPPICInterFoam中，不是通过真正的三相体积分数为1来处理的。并且OpenFOAM对粒子处理的本质就是忽略粒子的形态，体积，只考虑粒子对流体的作用力。

Jasper 0

Tens

@zhe 对于欧拉场来讲，是单纯的两相流，因为考虑了颗粒体积，所以连续相的体积只有alphac*Vcell，即Vc，所以连续相的计算是在这个基础上进行的，这也是为什么控制方程中都要乘alphac，可能结合单纯的MPPICFoam和InterFoam更好理解，只是将InterFoam的处理运用到MPPICFoam的连续相中

zhe

@tens 没错，意见一致。欧拉场中计算的所有都是在单纯两相流中乘以alphac的作用。所以在计算alpha1和alpha2时，OpenFOAM仍然按照单纯两相流计算，从interFoam中照搬来的。

Tens