of中的多相流求解器



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


  • 管理员

    $\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么?



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


  • 自由表面模型副教授 OpenFOAM讲师

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



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


  • 自由表面模型副教授 OpenFOAM讲师

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



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


  • 网格教授 OpenFOAM教授 管理员

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

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



  • @李东岳 ,首先我是看了您写的关于compressibleInterFoam求解器解析,然后推导完之后发现了个小小的错误。dgdt里面是没有alpha的。
    然后这个是我推导的可压缩两相方程:
    0_1493962493553_2017-05-05_13-33-57.bmp
    根据这个方程结合源码,这是修改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
            );
    

  • 网格教授 OpenFOAM教授 管理员

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

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



  • @李东岳 0_1493990651369_2017-05-05_21-16-12.bmp ,我说的错误是这个,里面是没有alpha1,alpha2的。对,算的结果是alpha1+alpha2最大值为2.。0_1493994510812_2017-05-05_22-06-09_副本.jpg 左边是原程序,右边是我改的程序。两者还是有点区别的。


  • 自由表面模型副教授 OpenFOAM讲师

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



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



  • @队长别开枪 这个是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
            );


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



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



  • @队长别开枪 ,最近做了一些测试,发现多相流求解器为什么求解每一相的同时保证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的通量计算。


  • 自由表面模型副教授 OpenFOAM讲师

    @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;
    }
    

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


  • 网格教授 OpenFOAM教授 管理员

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


  • 网格教授 OpenFOAM教授 管理员

    @mohuiof中的多相流求解器 中说:

    我说的错误是这个,里面是没有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}