求解相方程越界问题



  • 各位前辈们好!

    我目前尝试在compressibleInterFoam里面添加空化模型,前期经过热心前辈的指点,现在终于可以跑了,在此表示特别感谢!但是跑起来发现求解相方程alpha.water竟然越界了!!

    Courant Number mean: 0.000805796 max: 0.474105
    Time = 5e-07
    
    PIMPLE: iteration 1
    MULES: Solving for alpha.water
    Liquid phase volume fraction = 1  Min(alpha.water) = 1  Max(alpha.water) = 1
      Min(alpha.vapor) = -2.22045e-16  Max(alpha.vapor) = 2.22045e-16
    
    diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
    smoothSolver:  Solving for T, Initial residual = 0.97202, Final residual = 9.69394e-06, No Iterations 100
    min(T) 293.15
    GAMGPCG:  Solving for p_rgh, Initial residual = 1, Final residual = 0.000777804, No Iterations 1
    max(U) 10.3973
    min(p_rgh) 7567.84
    GAMGPCG:  Solving for p_rgh, Initial residual = 0.040294, Final residual = 4.4669e-05, No Iterations 1
    max(U) 10.3357
    min(p_rgh) 17859.6
    PIMPLE: iteration 2
    MULES: Solving for alpha.water
    Liquid phase volume fraction = 1  Min(alpha.water) = 1  Max(alpha.water) = 1
      Min(alpha.vapor) = -4.44089e-16  Max(alpha.vapor) = 4.44089e-16
    

    查找了好多论坛里面好多前辈讨论的帖子,自己分析是alpha方程里面源项的处理有问题,alpha方程如下:![D5TYQKE_0X{1C46}Z](LLAB.png](/assets/uploads/files/1593657820208-d5tyqke_0x-1c46-z-llab.png)

    我的alpha方程中源项离散代码如下:

    	Pair<tmp<volScalarField>> vDotAlphal =
            mixture->vDotAlphal();
        const volScalarField& vDotcAlphal = vDotAlphal[0]();
        const volScalarField& vDotvAlphal = vDotAlphal[1]();
        const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
    
        for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
        {
            volScalarField::Internal Sp
            (
                IOobject
                (
                    "Sp",
                    runTime.timeName(),
                    mesh
                ),
    		    vDotvmcAlphal  //这里参考的是interPhaseChangeFoam的alpha方程MULES:explicitSolve中Sp的设置
            );
    
            volScalarField::Internal Su
            (
                IOobject
                (
                    "Su",
                    runTime.timeName(),
                    mesh
                ),
                // Divergence term is handled explicitly to be
                // consistent with the explicit transport solution
                divU*min(alpha1, scalar(1))+vDotcAlphal   //这里参考的是interPhaseChangeFoam的alpha方程MULES:explicitSolve中Su的设置
            );
    
            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]*(1.0 - alpha1[celli]);
                }
            }
    
    
            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
            );
    

    等式右端前两项的离散参考的是东岳老师compressibleInterFoam解析中的方法,在此不做任何改动。最后一项空化源项的离散参考interPhaseChangeFoam中alpha方程中显式MULES中Sp和Su的处理,

    由于alpha.water求解的大于1,推断源项离散出了问题,查找资料发现源项离散要保证Sp<0,,,在这里我添加了vDotvmcAlphal (即蒸发率-冷凝率),这一项在不可压求解器的显式MULES中直接设置为了Sp,不可压求解器相方程对应代码如下:

    else
            {
                MULES::explicitSolve
                (
                    geometricOneField(),
                    alpha1,
                    phi,
                    talphaPhiCorr.ref(),
                    vDotvmcAlphal,   //Sp部分
                    (divU*alpha1 + vDotcAlphal)(),//Su部分
                    1,
                    0
                );
    

    按照道理我直接将这一项加到可压缩求解器中的Sp里面,逻辑上应该没有问题,,但确实alpha求解出现了越界问题

    麻烦各位前辈能不能指点一下我这个思路存在哪些问题呢?或者能够给我提供一点建议呢?十分感谢



  • @小考拉 alpha方程如下:61c83717-3135-4a66-9cdd-c1af5ba33159-image.png


Log in to reply
 

CFD中文网 2016 - 2020 | 京ICP备15017992号-2