MULES方法求解alpha方程出错，有谁知道怎样修改吗？

[mark]

[mark]

{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");
 
    surfaceScalarField phir("phir", phic*interface.nHatf());
 
    Pair<tmp<volScalarField> > vDotAlphal =
        twoPhaseProperties->vDotAlphal();
    const volScalarField& vDotcAlphal = vDotAlphal[0]();
    const volScalarField& vDotvAlphal = vDotAlphal[1]();
    const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
 
    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
    {
        volScalarField::DimensionedInternalField Sp
        (
            IOobject
            (
                "Sp",
                runTime.timeName(),
                mesh
            ),
            vDotvmcAlphal
        );
 
        volScalarField::DimensionedInternalField Su
        (
            IOobject
            (
                "Su",
                runTime.timeName(),
                mesh
            ),
            // Divergence term is handled explicitly to be
            // consistent with the explicit transport solution
            vDotcAlphal //divU*min(alpha1, scalar(1))
        );
 
        forAll(dgdt, celli)
        {
            if (dgdt[celli] >= 0.0 )
            {
                Su[celli] += (dgdt[celli]*alpha1[celli]*alpha2[celli]);
            }
            else
            {
                Sp[celli] += (dgdt[celli]*alpha2[celli]);
            }
            if (divU[celli] >= 0.0 )
            {
                Su[celli] += alpha1[celli]*divU[celli];
            }
            else
            {
                Sp[celli] += divU[celli];
            }
        }
 
 
        surfaceScalarField phiAlpha1
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );
 
        MULES::explicitSolve
        (
            geometricOneField(),
            alpha1,
            phi,
            phiAlpha1,
            Sp,
            Su,
            1,
            0
        );
 
        surfaceScalarField rho1f(fvc::interpolate(rho1));
        surfaceScalarField rho2f(fvc::interpolate(rho2));
        rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
 
        alpha1 = min(max(alpha1, scalar(0)), scalar(1));
        alpha2 = scalar(1) - alpha1;
    }
 
    Info<< "Liquid phase volume fraction = "
        << alpha1.weightedAverage(mesh.V()).value()
        << "  Min(" << alpha1.name() <<") = " << min(alpha1).value()
        << "  Min(" << alpha2.name() <<") = " << min(alpha2).value()
        << endl;
}

[mark]

[mark]

[mark]

1 , sp su 该怎样修改？

[mark]

2 sp su 怎样赋值给每个网格，而不是整体的源项？