of中的多相流求解器

mohui

@队长别开枪，最近做了一些测试，发现多相流求解器为什么求解每一相的同时保证alpha的总和有界。这里最关键的地方代码是limitsum通量限制器。这个是limitsum的源码，

void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
{
    forAll(phiPsiCorrs[0], facei)
    {
        scalar sumPos = 0;
        scalar sumNeg = 0;

        for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
        {
            if (phiPsiCorrs[phasei][facei] > 0)
            {
                sumPos += phiPsiCorrs[phasei][facei];
            }
            else
            {
                sumNeg += phiPsiCorrs[phasei][facei];
            }
        }

        scalar sum = sumPos + sumNeg;

        if (sum > 0 && sumPos > VSMALL)
        {
            scalar lambda = -sumNeg/sumPos;

            for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
            {
                if (phiPsiCorrs[phasei][facei] > 0)
                {
                    phiPsiCorrs[phasei][facei] *= lambda;
                }
            }
        }
        else if (sum < 0 && sumNeg < -VSMALL)
        {
            scalar lambda = -sumPos/sumNeg;

            for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
            {
                if (phiPsiCorrs[phasei][facei] < 0)
                {
                    phiPsiCorrs[phasei][facei] *= lambda;
                }
            }
        }
    }
}

看完代码后，可以比较清晰的看出，它其实做的工作量就是统计每一个面上的每一相的alpha的通量，然后进行了修正，得到了修正系数再乘上了alpha的通量，输出。

surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
        phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha)

这两行代码是指，最后这个修正的通量加上原来的通量，作为alpha的通量计算。

队长别开枪

@mohui OpenFOAM的另一个非官方版本（foam-extend）里alpha方程的源代码比官方的看起来简洁一点，比如求解器interFoam，extend版本（4.0）中的alphaEqn.H为

{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

    surfaceScalarField phic = mag(phi/mesh.magSf());
    phic = min(interface.cAlpha()*phic, max(phic));
    surfaceScalarField phir = phic*interface.nHatf();

    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
        surfaceScalarField phiAlpha =
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
                alpha1,
                alpharScheme
            );

        MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);

        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
    }

    Info<< "Liquid phase volume fraction = "
        << alpha1.weightedAverage(mesh.V()).value()
        << "  Min(alpha1) = " << min(alpha1).value()
        << "  Max(alpha1) = " << max(alpha1).value()
        << endl;
}

官方版本（4.x）中的为

{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

    tmp<fv::ddtScheme<scalar>> ddtAlpha
    (
        fv::ddtScheme<scalar>::New
        (
            mesh,
            mesh.ddtScheme("ddt(alpha)")
        )
    );

    // Set the off-centering coefficient according to ddt scheme
    scalar ocCoeff = 0;
    if
    (
        isType<fv::EulerDdtScheme<scalar>>(ddtAlpha())
     || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha())
    )
    {
        ocCoeff = 0;
    }
    else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()))
    {
        if (nAlphaSubCycles > 1)
        {
            FatalErrorInFunction
                << "Sub-cycling is not supported "
                   "with the CrankNicolson ddt scheme"
                << exit(FatalError);
        }

        ocCoeff =
            refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
           .ocCoeff();
    }
    else
    {
        FatalErrorInFunction
            << "Only Euler and CrankNicolson ddt schemes are supported"
            << exit(FatalError);
    }

    scalar cnCoeff = 1.0/(1.0 + ocCoeff);

    // Standard face-flux compression coefficient
    surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));

    // Add the optional isotropic compression contribution
    if (icAlpha > 0)
    {
        phic *= (1.0 - icAlpha);
        phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
    }

    surfaceScalarField::Boundary& phicBf =
        phic.boundaryFieldRef();

    // Do not compress interface at non-coupled boundary faces
    // (inlets, outlets etc.)
    forAll(phic.boundaryField(), patchi)
    {
        fvsPatchScalarField& phicp = phicBf[patchi];

        if (!phicp.coupled())
        {
            phicp == 0;
        }
    }

    tmp<surfaceScalarField> phiCN(phi);

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

    if (MULESCorr)
    {
        fvScalarMatrix alpha1Eqn
        (
            (
                LTS
              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            )
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
        );

        alpha1Eqn.solve();

        Info<< "Phase-1 volume fraction = "
            << alpha1.weightedAverage(mesh.Vsc()).value()
            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
            << endl;

        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
        alphaPhi = talphaPhiUD();

        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
        {
            Info<< "Applying the previous iteration compression flux" << endl;
            MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);

            alphaPhi += talphaPhiCorr0();
        }

        // Cache the upwind-flux
        talphaPhiCorr0 = talphaPhiUD;

        alpha2 = 1.0 - alpha1;

        mixture.correct();
    }


    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
        surfaceScalarField phir(phic*mixture.nHatf());

        tmp<surfaceScalarField> talphaPhiUn
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
               -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );

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

        if (MULESCorr)
        {
            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
            volScalarField alpha10("alpha10", alpha1);

            MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr.ref(), 1, 0);

            // Under-relax the correction for all but the 1st corrector
            if (aCorr == 0)
            {
                alphaPhi += talphaPhiCorr();
            }
            else
            {
                alpha1 = 0.5*alpha1 + 0.5*alpha10;
                alphaPhi += 0.5*talphaPhiCorr();
            }
        }
        else
        {
            alphaPhi = talphaPhiUn;

            MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
        }

        alpha2 = 1.0 - alpha1;

        mixture.correct();
    }

    if (alphaApplyPrevCorr && MULESCorr)
    {
        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
    }

    if
    (
        word(mesh.ddtScheme("ddt(rho,U)"))
     == fv::EulerDdtScheme<vector>::typeName
    )
    {
        rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
    }
    else
    {
        if (ocCoeff > 0)
        {
            // Calculate the end-of-time-step alpha flux
            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
        }

        // Calculate the end-of-time-step mass flux
        rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
    }

    Info<< "Phase-1 volume fraction = "
        << alpha1.weightedAverage(mesh.Vsc()).value()
        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
        << endl;
}

具体区别我不大清楚，但看起来还是有区别的。

李东岳

@mohui
有点迷惑。
如果你要验证求解一个相方程和两个相方程的区别，你就不需要研究MULES::limiter()，
如果你要验证MULES的限制器，你可以直接采用求解单alpha方程来仔细研究通量限制器。
哪个是研究的重点？:confused:

李东岳

@mohui 在 of中的多相流求解器中说：

我说的错误是这个，里面是没有alpha1，alpha2的

是说 compreisblInterFoam解析中下面这个方程有误？

\begin{equation}
\mathrm{dgdt}= \left( \frac{\alpha_2}{\rho_2}\frac{\mathrm{D}\rho_2}{\mathrm{D}t} - \frac{\alpha_1}{\rho_1}\frac{\mathrm{D}\rho_1}{\mathrm{D}t} \right)
\end{equation}

mohui

@李东岳对，这个我推导里面是没有alpha的，源程序里也是没有的。

mohui

@队长别开枪这两个区别在于后者采用了CrankNicolson格式。现在我最新给的代码是解释多相流求解器如何保证alpha之和有界的，你可以看下

李东岳

@mohui 非常感谢！原文已更新:cheeky:

有关MULES的验证（暂时说不上研究，因为基金会已经植入了），我这里分为俩部分验证，第一部分验证已经完毕并且投稿，第二部分验证还在进行中，不过进度很慢毕竟还有别的课题在进行，我只能建议你参考Flux Correct Transport部分数值算法内容，MULES的思想来源于此，只不过基金会将此方法用于多相流

有关求解一个方程或者两个方程的验证，应该不是很难，并且结果应该一致。在你之前发的相图上，我没有看出太大差别，是否有曲线图。另外，是否对比过1D Riemann问题的解析解？

mohui

@李东岳东岳老师，你可能还对我提的问题存在误区，我并不是对mules做验证，仅是分析OF中多相流求解器和两相流求解器对与两相流求解的区别。以我目前的测试检验，认为多相流求解器每一相求解并能够同时保证总的alpha之和等于的1的关键在于limitsum这个限制器，关于这个限制器的代码我也上传了，通过代码可以分析的出来，就是在mules求解之前对alphaphi进行了修正，以保证流入流出单元的alphaphi总和为0。其次，就如我前面所言，如果单纯的mules求解每一相是无法保证alpha1+alpha2=1的，虽然我给的图上没有很大区别，但是我输出alpha1+alpha2最大值为2，这是不对的。以上是我这几天的发现。

maoyanjun_dut

@李东岳

@李东岳在 of中的多相流求解器中说：

@mohui 非常感谢！原文已更新:cheeky:

有关MULES的验证（暂时说不上研究，因为基金会已经植入了），我这里分为俩部分验证，第一部分验证已经完毕并且投稿，第二部分验证还在进行中，不过进度很慢毕竟还有别的课题在进行，我只能建议你参考Flux Correct Transport部分数值算法内容，MULES的思想来源于此，只不过基金会将此方法用于多相流

有关求解一个方程或者两个方程的验证，应该不是很难，并且结果应该一致。在你之前发的相图上，我没有看出太大差别，是否有曲线图。另外，是否对比过1D Riemann问题的解析解？

岳哥，公式20更新了，公式19却没更新。建议更新一下，否则看着怪怪的。

赵一铭

@maoyanjun_dut

确实！好眼力。
@李东岳

yhdthu

@李东岳前辈，compressibleInterFoam解析的(19)式忘记改了，括号内分子没有alpha:laughing:

李东岳

谢谢，刚改了，interPhaseChangeFoam一直没时间写，估计得10月份开会回来之后能写...
:big_mouth:

yhdthu

@李东岳好的👌，到时候跟前辈讨论学习一下～

Jacobian

@mohui 没太看懂limitSum为何能保证体积分数之和为1，可否讲解一下？

dzw05

@mohui 其实也没太懂为什么limitsum可以保证alpha总和有界。limitsum限制的是所谓的修正通量，即“高阶通量”减去“迎风通量”的部分，至于为什么这部分通量总和等于零就可以保证alpha总和有界，这个还是没看懂。

小龙

@李东岳 @mohui @赵一铭我最近再做空化模拟，OF中只有多相流VOF模型，VOF做空化不太适合，想在OF里编译一下mixture模型，前辈们觉得实现的可行性高吗。实现的难度大吗？我初步的想法是基于interPhaseChangeFoam来修改。前辈们有什么建议？

李东岳

如果单纯的mules求解每一相是无法保证alpha1+alpha2=1的

现在我就想确定，alpha2通过对应方程推导，用MULES求解不可行的吗？

MULES的出发点是保证变量的有界。所以如果你用MULES求解alpha1，那么理论上alpha1是有界的。alpha2的求解没有调用MULES，他是通过1-alpha1算出来的，如果alpha1有界，alpha2也有界。

如果你用算法同时求解alpha1和alpha2，如何处理耦合？这俩个变量是耦合在一起的，就像速度和压力。你不能单独的去分离求解，如果分离求解就就需要迭代。迭代就导致计算速度变慢。因此现存大厂据我所知都是只求解alpha1，然后alpha2=1-alpha1。

我在去你年底验证了MULES，理论上是可以保证有界，但是真实计算的时候，还是会越界。尤其是物理模型比较复杂的时候。个人觉得这方面内容搞出来，绝对是个好文。

hehe

此回复已被删除！

CFD中文网

of中的多相流求解器