模型背景
我正在使用 OpenFOAM 的 COxidationKineticDiffusionLimitedRate 模型：
配置如下：（Sb, C1, C2, E），温度约在 1250 K；
从源码可知 dmC 的计算包含 Ap·rhoc、RR·Tc·YO₂/WO₂、D₀·Rₖ/(D₀+Rₖ)，其组合单位似乎为 J/m，但源码最终把输出 dmC 当作 碳质量消耗 kg 来处理；
#include "COxidationKineticDiffusionLimitedRate.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CloudType>
Foam::COxidationKineticDiffusionLimitedRate<CloudType>::
COxidationKineticDiffusionLimitedRate
(
const dictionary& dict,
CloudType& owner
)
:
SurfaceReactionModel<CloudType>(dict, owner, typeName),
Sb_(this->coeffDict().getScalar("Sb")),
C1_(this->coeffDict().getScalar("C1")),
C2_(this->coeffDict().getScalar("C2")),
E_(this->coeffDict().getScalar("E")),
CsLocalId_(-1),
O2GlobalId_(owner.composition().carrierId("O2")),
CO2GlobalId_(owner.composition().carrierId("CO2")),
WC_(0.0),
WO2_(0.0),
HcCO2_(0.0)
{
// Determine Cs ids
label idSolid = owner.composition().idSolid();
CsLocalId_ = owner.composition().localId(idSolid, "C");
// Set local copies of thermo properties
WO2_ = owner.thermo().carrier().W(O2GlobalId_);
const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
WC_ = WCO2 - WO2_;
HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
Info<< " C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
}
template<class CloudType>
Foam::COxidationKineticDiffusionLimitedRate<CloudType>::
COxidationKineticDiffusionLimitedRate
(
const COxidationKineticDiffusionLimitedRate<CloudType>& srm
)
:
SurfaceReactionModel<CloudType>(srm),
Sb_(srm.Sb_),
C1_(srm.C1_),
C2_(srm.C2_),
E_(srm.E_),
CsLocalId_(srm.CsLocalId_),
O2GlobalId_(srm.O2GlobalId_),
CO2GlobalId_(srm.CO2GlobalId_),
WC_(srm.WC_),
WO2_(srm.WO2_),
HcCO2_(srm.HcCO2_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
(
const scalar dt,
const scalar Re,
const scalar nu,
const label celli,
const scalar d,
const scalar T,
const scalar Tc,
const scalar pc,
const scalar rhoc,
const scalar mass,
const scalarField& YGas,
const scalarField& YLiquid,
const scalarField& YSolid,
const scalarField& YMixture,
const scalar N,
scalarField& dMassGas,
scalarField& dMassLiquid,
scalarField& dMassSolid,
scalarField& dMassSRCarrier
) const
{
// Fraction of remaining combustible material
const label idSolid = CloudType::parcelType::SLD;
const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
// Surface combustion active combustible fraction is consumed
if (fComb < SMALL)
{
return 0.0;
}
const SLGThermo& thermo = this->owner().thermo();
// Local mass fraction of O2 in the carrier phase
const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[celli];
// Diffusion rate coefficient
const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
// Kinetic rate
const scalar Rk = C2_*exp(-E_/(RR*Tc));
// Particle surface area
const scalar Ap = constant::mathematical::pi*sqr(d);
// Change in C mass [kg]
scalar dmC = Ap*rhoc*RR*Tc*YO2/WO2_*D0*Rk/(D0 + Rk)*dt;
// Limit mass transfer by availability of C
dmC = min(mass*fComb, dmC);
// Molar consumption
const scalar dOmega = dmC/WC_;
// Change in O2 mass [kg]
const scalar dmO2 = dOmega*Sb_*WO2_;
// Mass of newly created CO2 [kg]
const scalar dmCO2 = dOmega*(WC_ + Sb_*WO2_);
// Update local particle C mass
dMassSolid[CsLocalId_] += dOmega*WC_;
// Update carrier O2 and CO2 mass
dMassSRCarrier[O2GlobalId_] -= dmO2;
dMassSRCarrier[CO2GlobalId_] += dmCO2;
const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
// carrier sensible enthalpy exchange handled via change in mass
// Heat of reaction [J]
return dmC*HsC - dmCO2*HcCO2_;
}
COxidationKineticDiffusionLimitedRate 是 OpenFOAM 中用于描述 焦炭表面碳 (C(s)) 与氧气 (O₂) 反应生成 CO₂ 的动力学–扩散受限表面反应模型。
该模型只支持单一路径：C(s) + Sb·O₂ → CO₂；无法处理生成 CO 或 H₂ 的气化反应（如 char + CO₂ 或 char + H₂O）。(OpenFOAM)
社区用户曾指出，若使用 tutorial 中默认参数（Sb=1、C1=5e‑12、C2=0.002、E=7.9e7），在合理温度（如 Tc = 1250 K）下计算 Rk → 0，从而导致 dmC/dt → 0，看似模型不起作用。(CFD Online)
核心关注点
1. dmC 的单位与物理意义：
公式如下（简化版本）：
dmC = Ap * rhoc * RR * Tc * YO2 / WO2_ * (D0 * Rk / (D0 + Rk)) * dt;
各项单位：
Ap * rhoc → kg/m
RR*Tc*YO2/WO2_ → J/kg 或 J·kmol/kg²
D0 * Rk / (D0 + Rk) → 1/s
dt → s
组合后似乎是 J/m（能量密度单位），而非质量单位（kg），令人难以直接信服 dmC 实为 kg 碳质量损耗。
社区质疑这一点，但也指出源码后续使用 dmC/WC 计算摩尔变化，再映射为气相 O₂ 与 CO₂ 质量，最终模型确实是按质量守恒执行。(CFD Online)
2. 模型常参数来源与适用性
默认参数值（tutorial 案例中提供）在合理温度下几乎不启动：
当 Tc = 1250 K，Rk ≈ 0，导致 dmC/dt = 0，即不发生碳消耗；
需要极高温度（如 Tc → 50000 K）才会有所反应，这显然不现实。(CFD Online)
参数来源据称为 Baum & Street (1971) 研究，但社区提问第二点是：
确切常数来源不明；
模型中使用 Tc（载体温度）而非 Tp（碳颗粒表面温度），与文献描述不一致。(CFD Online)
总结
我对以下三点非常困惑，想向大家询问：
单位体系是否正确？
虽然中间结果为 J/m，但 dmC 最终确实代表 kg 碳损耗？有哪些地方负责单位换算以确保质量守恒？
为何采用载体温度 Tc 而非颗粒表面温度 Tp？
文献模型多数以 Tp（carbon particle temperature）为基础，但代码却用 Tc，是否存在误用或简化？
默认参数值是否合理
常数如 C2=0.002, E=7.9e7 是否适用？
