of中的多相流求解器



  • @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的通量计算。



  • @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:



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



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



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



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

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

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



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



  • @李东岳

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

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

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

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

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


  • 管理员

    @maoyanjun_dut

    确实!好眼力。
    @李东岳



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



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



  • @李东岳 好的👌,到时候跟前辈讨论学习一下~


登录后回复
 

与 CFD 中国 的连接断开,我们正在尝试重连,请耐心等待