of中的多相流求解器

[mohui]

在of中我有个疑问，关于多相流求解器。目前来看，两相流求解器，对于求解alpha方程都是第二相用1-alpha1得到，这样子虽然保证体积守恒怎么保证质量守恒，特别是有相变的求解器（interphasechangeFoam）。另外of中多相流求解器中的alpha方程是套了循环，也就是说每一相都求解了，那么最后是怎么保证每一个cell中所有的相的体积分数总和为1?

[CFD中文网]

$\alpha_1=1-\alpha_2$是用来保证alpha有界的，因为：
\begin{equation}
0<\alpha<1
\end{equation}
守恒体现在连续性方程中：
\begin{equation}
\frac{\partial \rho\alpha}{\partial t}+\nabla\cdot(\rho\mathbf{U}\alpha)=0
\end{equation}

of中多相流求解器中的alpha方程是套了循环

具体指的是什么？MULES么？

[mohui]

@cfd-china OF中的多相流求解器是求解n相，所以需要套循环，alpha求解是用mules求解，然后我仔细看了多相流求解器，发现除了mules求解每一相之外，还有mules limitsum修正，这个不知道我理解对不对。理由是，我做了个测试，在compressibleinterFoam求解器中，对应的第二相也用mules求解，但是无法保证总的体积分数为1，算到一定时间算崩了。所以我推断多相流求解中的mules limitsum是保证总的体积分数守恒，关于这个理解不知道对不对，求解答。谢谢

[队长别开枪]

VOF方程求解分为几何重构和纯代数方法两类，OpenFOAM中的MULES求解方法属于纯代数方法，具有效率高、易于实现、适用于任意多面体网格等优点，类似的还有CICSAM方法，代数类方法的求解重点在于计算单元表面体积流率和克服数值弥散。据我所知，在interFoam和interPhaseChangeFoam中都是只求解主相alpha1而且都只支持两相流动。

[mohui]

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

[队长别开枪]

@mohui 这个有点复杂，等我今天开完组会我们一起讨论一下。

[mohui]

@队长别开枪 恩，我做了测试，如果第二相采用mules求解会导致alphasum无界。

[李东岳]

@mohui
如果alpha1和alpha2采用完全相同的数值处理，并且边界条件给定的互相对应，求解alpha1和alpha2的方程和求解alpha1方程并alpha2=1-alpha1应该是相同的。

可以贴出来你做的代码改动，你用什么算例测试的？

[mohui]

@李东岳 ，首先我是看了您写的关于compressibleInterFoam求解器解析，然后推导完之后发现了个小小的错误。dgdt里面是没有alpha的。
然后这个是我推导的可压缩两相方程：
，
根据这个方程结合源码，这是修改compressibleInterFoam中的alpha方程的代码：

for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
    {
        volScalarField::Internal Sp
        (
            IOobject
            (
                "Sp",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
        );

        volScalarField::Internal Su
        (
            IOobject
            (
                "Su",
                runTime.timeName(),
                mesh
            ),
           
            divU*min(alpha1, scalar(1))
        );

        forAll(dgdt, celli)
        {
            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
            {
                Sp[celli] -= dgdt[celli]*alpha1[celli];
                Su[celli] += dgdt[celli]*alpha1[celli];
            }
            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
            {
                Sp[celli] += dgdt[celli]*(max(1.0 - alpha1[celli],scalar(0)));
            }
        }


        surfaceScalarField alphaPhi1
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );

        MULES::explicitSolve
        (
            geometricOneField(),
            alpha1,
            phi,
            alphaPhi1,
            Sp,
            Su,
            1,
            0
        );

        //slove for alpha2
        volScalarField::Internal Sp1
        (
            IOobject
            (
                "Sp1",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedScalar("Sp1", dgdt.dimensions(), 0.0)
        );

        volScalarField::Internal Su1
        (
            IOobject
            (
                "Su1",
                runTime.timeName(),
                mesh
            ),
            
            divU*min(alpha2, scalar(1))
        );
       
        forAll(dgdt, celli)
        {
            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
            {
                Sp1[celli] -= dgdt[celli]*alpha1[celli];
                
            }
            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
            {
                Sp1[celli] += dgdt[celli]*alpha2[celli];
                Su1[celli] -= dgdt[celli]*(max(1.0 - alpha1[celli],scalar(0)));
            }
        }
         surfaceScalarField alphaPhi2
        (
            fvc::flux
            (
                phi,
                alpha2,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(phir, alpha1, alpharScheme),
                alpha2,
                alpharScheme
            )
        );
         MULES::explicitSolve
        (
            geometricOneField(),
            alpha2,
            phi,
            alphaPhi2,
            Sp1,
            Su1,
            1,
            0
        );

[李东岳]

可以把我那个错误得地方指出来么我好改一下。

我看了下你的方程和贴出来的代码没看出问题。能否把alpha1方程去掉，单独求解alpha2? 另外你的错误结果是什么？alpha1+alpha2>1么？还是结果严重不符合物理。

[mohui]

@李东岳 ，我说的错误是这个，里面是没有alpha1，alpha2的。对，算的结果是alpha1+alpha2最大值为2.。 左边是原程序，右边是我改的程序。两者还是有点区别的。

[队长别开枪]

@mohui Hello, 不好意思现在才回复你。OpenFOAM的MULES求解器我不大熟悉，我自己本身是做多面体网格上的PLIC算法的，最后评估算法优劣的参数中就有一个mass/volume conservation error，这个值越小越好，如果算法能够保证这个参数恒等于零，那么分别求解alpha1和alpha2没有什么问题，但现实中算法误差和数值误差的存在，这个值不可能恒等于零，所以一般次相alpha2直接用1-alpha1计算，如果这个时候分别计算alpha1和alpha2，我觉得会叠加误差，以至于出现你遇到的alpha1+alpha2>1的情况，更新密度粘度场后误差再传进动量方程和压力泊松方程，最后造成求解发散。这些是我的看法。

[mohui]

@队长别开枪 恩，确实是这样子的。

[mohui]

@队长别开枪 这个是compressiblemultiphaseFoam的关于求解alpha的源码，从源码看，除了用explicitsolve，在这前还使用了limit以及limitsum 限制，我在想是不是这样子才会使得求解每一相同时还能保证alpha之和为1.

void Foam::multiphaseMixtureThermo::solveAlphas
(
    const scalar cAlpha
)
{
    static label nSolves=-1;
    nSolves++;

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

    surfaceScalarField phic(mag(phi_/mesh_.magSf()));
    phic = min(cAlpha*phic, max(phic));

    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
    int phasei = 0;

    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
    {
        phaseModel& alpha = phase();

        phiAlphaCorrs.set
        (
            phasei,
            new surfaceScalarField
            (
                fvc::flux
                (
                    phi_,
                    alpha,
                    alphaScheme
                )
            )
        );

        surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];

        forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
        {
            phaseModel& alpha2 = phase2();

            if (&alpha2 == &alpha) continue;

            surfaceScalarField phir(phic*nHatf(alpha, alpha2));

            phiAlphaCorr += fvc::flux
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha,
                alpharScheme
            );
        }

        MULES::limit
        (
            1.0/mesh_.time().deltaT().value(),
            geometricOneField(),
            alpha,
            phi_,
            phiAlphaCorr,
            zeroField(),
            zeroField(),
            1,
            0,
            3,
            true
        );

        phasei++;
    }

    MULES::limitSum(phiAlphaCorrs);

    rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);

    volScalarField sumAlpha
    (
        IOobject
        (
            "sumAlpha",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
        dimensionedScalar("sumAlpha", dimless, 0)
    );


    volScalarField divU(fvc::div(fvc::absolute(phi_, U_)));


    phasei = 0;

    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
    {
        phaseModel& alpha = phase();

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

        volScalarField::DimensionedInternalField Sp
        (
            IOobject
            (
                "Sp",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
            dimensionedScalar("Sp", alpha.dgdt().dimensions(), 0.0)
        );

        volScalarField::DimensionedInternalField Su
        (
            IOobject
            (
                "Su",
                mesh_.time().timeName(),
                mesh_
            ),
            // Divergence term is handled explicitly to be
            // consistent with the explicit transport solution
            divU*min(alpha, scalar(1))
        );

        {
            const scalarField& dgdt = alpha.dgdt();

            forAll(dgdt, celli)
            {
                if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
                {
                    Sp[celli] += dgdt[celli]*alpha[celli];
                    Su[celli] -= dgdt[celli]*alpha[celli];
                }
                else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
                {
                    Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
                }
            }
        }

        forAllConstIter(PtrDictionary<phaseModel>, phases_, phase2)
        {
            const phaseModel& alpha2 = phase2();

            if (&alpha2 == &alpha) continue;

            const scalarField& dgdt2 = alpha2.dgdt();

            forAll(dgdt2, celli)
            {
                if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
                {
                    Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
                    Su[celli] += dgdt2[celli]*alpha[celli];
                }
                else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
                {
                    Sp[celli] += dgdt2[celli]*alpha2[celli];
                }
            }
        }

        MULES::explicitSolve
        (
            geometricOneField(),
            alpha,
            phiAlpha,
            Sp,
            Su
        );

[mohui]

后来我又做了测试，把limit以及limitsum注释掉了，重新编译运算，发现alpha的总和仍然是守恒的，按道理我注释掉了这两个，说白了也是只有explicitsolve求解，运算的时候结果很好的保证了有界性，还请各位大神能否解释这是为什么。

[mohui]

不好意思，昨天编译的仓促，实际上并没有真正的编译我所改掉的程序，所以运行的还是原来的程序。直到刚才我才正确重新编译，运算发现马上浮点溢出。说明limit和limtsum确实有作用。

[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}