关于MULES的疑问



  • Dear FOAMers:

    大家好,本人最近在看MULES算法,有些疑问,希望大家能为我解惑,可能比较长,首先表示感谢:papa:

    我使用的是interFOAM求解器,以OF6为例,为了简便,假设MULESCorr为false,在alphaEqn.H中的求解过程为:

        for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
        {
            #include "alphaSuSp.H"               //声明源项,全部为zeroField类型
    
            surfaceScalarField phir(phic*mixture.nHatf());      //nHatf为自由表面的单位法向量
    
            tmp<surfaceScalarField> talphaPhi1Un        //组建alpha通量,采用时间项若采用Euler,则cnCoeff= 1.0
            (
                fvc::flux
                (
                    phiCN(),
                    cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
                    alphaScheme
                )
              + fvc::flux  //压缩项
                (
                   -fvc::flux(-phir, alpha2, alpharScheme),
                    alpha1,
                    alpharScheme
                )
            );
    
            /** 省略一部分 */
    
            else
            {
                alphaPhi10 = talphaPhi1Un;
    
                MULES::explicitSolve                     
                (
                    geometricOneField(),                                          //此处不太理解,请看后文
                    alpha1,
                    phiCN,
                    alphaPhi10,
                    Sp,
                    (Su + divU*min(alpha1(), scalar(1)))(),
                    oneField(),
                    zeroField()
                );
            }
    	
    	//#include "setAlphaMR.H"
    
            alpha2 = 1.0 - alpha1;
    
            mixture.correct();
        }
    

    观察MULES::explicitSolve函数,其传入8个参数,查阅其代码,有:

    void Foam::MULES::explicitSolve
    (
        const RhoType& rho,         //对应前面的geometricOneField()
        volScalarField& psi,           //alpha1
        const surfaceScalarField& phi,    
        surfaceScalarField& phiPsi,
        const SpType& Sp,      //显式计算的源项?
        const SuType& Su,     //隐式计算的源项?
        const PsiMaxType& psiMax,
        const PsiMinType& psiMin
    )
    {
        const fvMesh& mesh = psi.mesh();
    
        psi.correctBoundaryConditions();
    
        if (fv::localEulerDdt::enabled(mesh))
        {
            const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
            limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false);
            explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
        }
        else
        {
            const scalar rDeltaT = 1.0/mesh.time().deltaTValue();  //时间步的倒数
            limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false);  //这个函数非常长,没仔细看不知道什么意思,但估计是保证alpha的结果在0-1之间?
            explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su); //最后更新值还是调用的这个六参数的函数
        }
    }
    

    关于上述六参数的函数:

    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); //这里将通量除以了体积并幅值给psiIf
    
       //省略一小部分
    
        else
        {
            psiIf =
            (
                rho.oldTime().field()*psi0*rDeltaT
              + Su.field()
              - psiIf
            )/(rho.field()*rDeltaT - Sp.field());
        }
    
        psi.correctBoundaryConditions();
    }
    

    我试图把上述else中的代码变成公式:
    0_1542966236253_14b0163c-8163-4d9b-bd90-88c93ac0f089-image.png

    此公式与离散相方程后的结果类似,只是多出了rho,查阅第一段代码可以发现,实际应用的时候传入的参数为geometricOneField(),即为1.0,因此这其实就是一个系数,但是我不明白这样定义的意义在哪里?是有其他的公式会用到这个参数吗?为什么起了密度这个名字?这个地方一开始把我搞晕了好久。

    此外,在interFOAM中,Su和Sp均为0,我查看代码发现其一个表示显式源项,另一个表示隐式源项,如果我想在相方程中加入显式源项,请问该如何添加呢?直接调用函数的时候把我的源项表示的场传进去(给Su)就行了吗?

    还有,上述代码中的注释都是我自己的理解,发帖时写上去的,请大神们看看有没有哪里理解错误的,请批评指正,谢谢!

    不知道我说明白了没,希望大家为我解惑,谢谢!当然,如果有其他关于MULES的问题也可以在这里一起讨论!


  • 管理员

    因此这其实就是一个系数,但是我不明白这样定义的意义在哪里?是有其他的公式会用到这个参数吗?为什么起了密度这个名字?这个地方一开始把我搞晕了好久。

    在某些情况下需要传入密度,这里只不过假定不可压缩,为geometricOneField()

    Su和Sp均为0,我查看代码发现其一个表示显式源项,另一个表示隐式源项,如果我想在相方程中加入显式源项,请问该如何添加呢?直接调用函数的时候把我的源项表示的场传进去(给Su)就行了吗?

    理论上是可以的,数值稳定性问题需要进一步研究

    你的理解是正确的。图中公式也是正确的。另外,limit函数是用来求$\lambda$的值,可以参考Boris and book 1973年关于Flux correct transport的文章,对应文章中的$\lambda$



  • @东岳 感谢东岳老师的回复,关于源项的问题,之前我是这么解决的:

    if (MULESCorr)
        6     {
        7         fvScalarMatrix alpha1Eqn
        8         (
        9             fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
       10           + fv::gaussConvectionScheme<scalar>
       11             (
       12                 mesh,
       13                 phi,
       14                 upwind<scalar>(mesh, phi)
       15             ).fvmDiv(phi, alpha1)
    
                     //在这个位置添加显式的源项表达式
    
       16         );
       17 
       18         solve(alpha1Eqn);
       19 
       20         Info<< "Phase-1 volume fraction = "
       21             << alpha1.weightedAverage(mesh.Vsc()).value()
       22             << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
       23             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
       24             << endl;
    

    之前在上述代码注释处添加了源项(那时候看不懂代码,就发现这里是相方程。。。),这个源项比较简单,就是一个时间的简单函数,且之前的结果看来稳定性还算可以,不过当时写程序的时候并没有在后面的MULES::correct中加入源项表达式,目前没有评估出这个差异会有多大。

    此外,东岳老师,我还有个小疑问,这个bool变量MULESCorr决定的两种求解方式的区别大吗?您有做过类似的测试吗?谢谢!


  • 讲师



  • @yhdthu

    感谢,没想到知乎上也有大神研究这个!


Log in to reply
 


CFD中文网 | 东岳流体学术 | 东岳流体商业 | 吉ICP备20003622号-1