计算气固颗粒反应流动，稠密颗粒流动的基础上添加了气固反应，组分输运方程的代码如下：
{
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));
}
到目前为止，我也没看出来有啥问题。所以，想请教一下吧里的大佬，能否帮我解答一下疑惑？如何解决这个问题，请各位大佬不吝赐教！！！
