气固颗粒反应组分输运方程问题

[书生kao]

计算气固颗粒反应流动，稠密颗粒流动的基础上添加了气固反应，组分输运方程的代码如下：

{
    combustion->correct();
    Qdot = combustion->Qdot();
    volScalarField Yt(0.0*Y[0]);

    const volScalarField muEff(turbulence->muEff());

    forAll(Y, i)
    {
        if (i != inertIndex && composition.active(i))
        {
            volScalarField& Yi = Y[i];

            fvScalarMatrix YiEqn
            (
                fvm::ddt(alphac, rhoc, Yi)
              + mvConvection->fvmDiv(alphaRhoPhic, Yi)
              - fvm::laplacian
                (
                    fvc::interpolate(alphac)
                    *fvc::interpolate(muEff),
                    Yi
                )
              ==
                vanadiumParcels.SYi(i, Yi)
              + combustion->R(Yi)
              + fvOptions(rhoc, Yi)
            );

            YiEqn.relax();

            fvOptions.constrain(YiEqn);

            YiEqn.solve("Yi");

            fvOptions.correct(Yi);

            Yi.max(0.0);
            Yt += Yi;
        }
    }

    Y[inertIndex] = scalar(1) - Yt;
    Y[inertIndex].max(0.0);
}

编译是没有问题的，目前遇到的问题是各组分的质量分数不对。
反应是生成Cl2，消耗O2，inertSpecie为N2，计算过程氯气的质量分数不断增加（未越界），而N2的质量分数减少。实际上即便完全反应，N2的质量分数下降应该不大，计算结果不符合实际现象。
个人理解问题应该是vanadiumParcels.SYi(i, Yi)，可以理解为组分输运方程的源项，其代码为

template<class CloudType>
inline Foam::tmp<Foam::fvScalarMatrix> Foam::ReactingCloud<CloudType>::SYi
(
    const label i,
    volScalarField& Yi
) const
{
    if (this->solution().coupled())
    {
        if (this->solution().semiImplicit("Yi"))
        {
            tmp<volScalarField> trhoTrans
            (
                volScalarField::New
                (
                    this->name() + ":rhoTrans",
                    this->mesh(),
                    dimensionedScalar(dimMass/dimTime/dimVolume, 0)
                )
            );

            volScalarField& sourceField = trhoTrans.ref();

            sourceField.primitiveFieldRef() =
                rhoTrans_[i]/(this->db().time().deltaTValue()*this->mesh().V());

            const dimensionedScalar Yismall("Yismall", dimless, small);

            return
                fvm::Sp(neg(sourceField)*sourceField/(Yi + Yismall), Yi)
              + pos0(sourceField)*sourceField;
        }
        else
        {
            tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
            fvScalarMatrix& fvm = tfvm.ref();

            fvm.source() = -rhoTrans_[i]/this->db().time().deltaTValue();

            return tfvm;
        }
    }

    return tmp<fvScalarMatrix>(new fvScalarMatrix(Yi, dimMass/dimTime));
}

到目前为止，我也没看出来有啥问题。所以，想请教一下吧里的大佬，能否帮我解答一下疑惑？如何解决这个问题，请各位大佬不吝赐教！！！