各位老师们好，目前我在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;
想请教一下大家有没有什么更改的建议。
