求解相方程越界问题

[小考拉]

各位前辈们好！

我目前尝试在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方程如下：

我的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方程如下：