请教贴：拉格朗日耦合进密度基求解器，能量方程报错

[Rachel0096]

各位老师们好，目前我在OpenFOAM-v1912中尝试把一个压力基的拉格朗日两相反应求解器（基于coalChemstryFoam）耦合进detonationFoam（密度基）中。目前编译没有问题，可以顺利通过，但是在测试算例时运行会报错（使用的是NS_Sutherland模型）：

=== Start Sensible Enthalpy Transport ====


--> FOAM FATAL ERROR:
incompatible fields for operation
    [rhoE] - [h]

    From function void Foam::checkMethod(const Foam::fvMatrix<Type>&, const Foam::fvMatrix<Type>&, const char*) [with Type = double]
    in file /share/home/zhoufan/OpenFOAM/OpenFOAM-v1912/src/finiteVolume/lnInclude/fvMatrix.C at line 1337.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  void Foam::checkMethod<double>(Foam::fvMatrix<double> const&, Foam::fvMatrix<double> const&, char const*) at ??:?
#3  Foam::tmp<Foam::fvMatrix<double> > Foam::operator-<double>(Foam::tmp<Foam::fvMatrix<double> > const&, Foam::tmp<Foam::fvMatrix<double> > const&) at ??:?
#4  ? at ??:?
#5  __libc_start_main in /lib64/libc.so.6
#6  ? at ??:?
Aborted (core dumped)

这里显示动量方程和组分方程求解都没有问题，在求解能量方程时遇到了赋值不统一的问题。
这是原本detonationFoam中rhoEEqn的代码：

surfaceScalarField sigmaDotU
(
    "sigmaDotU",
    (
        fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
      + fvc::dotInterpolate(mesh.Sf(), tauMC)
    )
  & fvc::interpolate(U)    
);

solve
(
    fvm::ddt(rhoE)
  + fvc::div(rhoEPhi)
  - fvc::div(sigmaDotU)
  ==
    reaction->Qdot()
);

e = rhoE/rho - 0.5*magSqr(U);
e.correctBoundaryConditions();
thermo.correct();
rhoE.boundaryFieldRef() == rho.boundaryField()*(e.boundaryField()+0.5*magSqr(U.boundaryField()));


solve
(
    fvm::ddt(rho, e) - fvc::ddt(rho, e)
  + thermophysicalTransport->divq(e)
);
thermo.correct();
rhoE = rho*(e + 0.5*magSqr(U));

p.ref() = rho()/psi();
p.correctBoundaryConditions();
rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();

Info<< "min/max(p) = "<< min(p).value() << ", " << max(p).value() << endl;
Info<< "min/max(T) = "<< min(T).value() << ", " << max(T).value() << endl;

这是我修改添加颗粒源项之后的代码：

clock_t t_begin_h = std::clock();
Info <<  "=== Start Sensible Enthalpy Transport ====" << endl;
auto dt = runTime.deltaT();

surfaceScalarField sigmaDotU
(
    "sigmaDotU",
    (
        fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
      + fvc::dotInterpolate(mesh.Sf(), tauMC)
    )
  & fvc::interpolate(U)    
);

solve
(
    fvm::ddt(rhoE)
  + fvc::div(rhoEPhi)
  - fvc::div(sigmaDotU)
  - fvm::laplacian(turbulence->mut()/Prt + turbulence->alpha(), e)
  ==
    combustion->Qdot()
  + AlParcels.Sh(e)
  + radiation->Sh(thermo, e)
);

e = rhoE/rho - 0.5*magSqr(U);
e.correctBoundaryConditions();
thermo.correct();
rhoE.boundaryFieldRef() == rho.boundaryField()*(e.boundaryField()+0.5*magSqr(U.boundaryField()));

solve
(
    fvm::ddt(rho, e) - fvc::ddt(rho, e)
);

thermo.correct();
radiation->correct(); //add by vv
rhoE = rho*(e + 0.5*magSqr(U));

p.ref() = rho()/psi();
p.correctBoundaryConditions();
rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();

Info<< "min/max(p) = "<< min(p).value() << ", " << max(p).value() << endl;
Info<< "min/max(T) = "<< min(T).value() << ", " << max(T).value() << endl;

clock_t t_end_h = std::clock();
double elapsed_secs_h = double(t_end_h - t_begin_h) / CLOCKS_PER_SEC;

Info <<  "=== Sensible Enthalpy Transport [Done]: "
     << elapsed_secs_h
     << " s ====\n"
     << endl;

想请教一下大家有没有什么更改的建议。