CFD Online Discussion Forums - Solving for an additional species CO in coalChemistryFoam

Hi,

Let me frist give you a brief background:

I have added multiple reactions for char combustion chemistry and modeified coalChemistryFoam. It was simple as I added following reactions in chemkin folder in chem.inp file:
ELEMENTS
O H C N
END
SPECIES
N2 O2 CH4 H2 CO2 H2O CO
END
REACTIONS
CH4 + 2 O2 => CO2 + 2 H2O 7.000E+12 0.00 20000
H2 + 0.5 O2 => H2O 1.500E+11 0.00 20000
CO + H2O => CO2 + H2 5.000E+12 0.00 20000
END

The fvSolution file already has details for CO. The amount of CO comes from Char combustion model and I have modified following file:

~/src/lagrangian/mycoalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C

Instead of CO2, CO is formed due to reaction of C with O2, H2O and CO2.

Questions:

1. I dont see that there is CO in the gas phase. Is there an issue in coupling for the two phases or there is an issue in the set-up of the case?

2. If anyone has implemented multiple-step char combustion chemistry in openfoam, can you please check the following file and let me know if you find anything.

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "COxidationMultiStepChemistry.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CloudType>
Foam::COxidationMultiStepChemistry<CloudType>::
COxidationMultiStepChemistry
(
const dictionary& dict,
CloudType& owner
)
:
SurfaceReactionModel<CloudType>
(
dict,
owner,
typeName
),
Sb_(dimensionedScalar(this->coeffDict().lookup("Sb")).value()),
C1_(dimensionedScalar(this->coeffDict().lookup("C1")).value()),
Ka_(dimensionedScalar(this->coeffDict().lookup("Ka")).value()),
Ea_(dimensionedScalar(this->coeffDict().lookup("Ea")).value()),
Kb_(dimensionedScalar(this->coeffDict().lookup("Kb")).value()),
Eb_(dimensionedScalar(this->coeffDict().lookup("Eb")).value()),
Kc_(dimensionedScalar(this->coeffDict().lookup("Kc")).value()),
Ec_(dimensionedScalar(this->coeffDict().lookup("Ec")).value()),
CsLocalId_(-1),
O2GlobalId_(owner.composition().globalCarrierId("O 2")),
CO2GlobalId_(owner.composition().globalCarrierId(" CO2")),
COGlobalId_(owner.composition().globalCarrierId("C O")),
H2OGlobalId_(owner.composition().globalCarrierId(" H2O")),
H2GlobalId_(owner.composition().globalCarrierId("H 2")),
WC_(0.0),
WO2_(0.0),
WH2_(0.0),
WCO2_(0.0),
WCO_(0.0),
WH2O_(0.0)
{
// Determine Cs ids
label idSolid = owner.composition().idSolid();
CsLocalId_ = owner.composition().localId(idSolid, "C");
// Set local copies of thermo properties
WO2_ = owner.mcCarrierThermo().speciesData()[O2GlobalId_].W();
WH2_ = owner.mcCarrierThermo().speciesData()[H2GlobalId_].W();
WCO2_ = owner.mcCarrierThermo().speciesData()[CO2GlobalId_].W();
WH2O_ = owner.mcCarrierThermo().speciesData()[H2OGlobalId_].W();
WC_ = WCO2_ - WO2_;
WCO_ = owner.mcCarrierThermo().speciesData()[COGlobalId_].W();
if (Sb_ < 0)
{
FatalErrorIn
(
"COxidationMultiStepChemistry"
"("
"const dictionary&, "
"CloudType&"
")"
) << "Stoichiometry of reaction, Sb, must be greater than zero" << nl
<< exit(FatalError);
}
}

// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class CloudType>
Foam::COxidationMultiStepChemistry<CloudType>::
~COxidationMultiStepChemistry()
{}

// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
bool Foam::COxidationMultiStepChemistry<CloudType>::act ive() const
{
return true;
}

template<class CloudType>
Foam::scalar Foam::COxidationMultiStepChemistry<CloudType>::cal culate
(
const scalar dt,
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 scalarField& dMassVolatile,
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;
}
// Local mass fraction of O2 in the carrier phase
const scalar YO2 = this->owner().mcCarrierThermo().Y(O2GlobalId_)[cellI];
const scalar YCO2 = this->owner().mcCarrierThermo().Y(CO2GlobalId_)[cellI];
const scalar YH2O = this->owner().mcCarrierThermo().Y(H2OGlobalId_)[cellI];
// Diffusion rate coefficient
const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
// Kinetic rate
const scalar Rka = max(Tc*(-0.0168+0.0000132*Tc) , 0) ;
const scalar Rkb = Kb_*exp(-Eb_/(specie::RR*Tc));
const scalar Rkc = Kc_*exp(-Ec_/(specie::RR*Tc));
// Particle surface area
const scalar Ap = mathematicalConstant:pi*sqr(d);
// Change in C mass [kg]
scalar dmCa = Ap*rhoc*specie::RR*Tc*YO2/WO2_*D0*Rka/(D0 + Rka);
scalar dmCb = Ap*rhoc*specie::RR*Tc*YCO2/WCO2_*D0*Rkb/(D0 + Rkb);
scalar dmCc = Ap*rhoc*specie::RR*Tc*YH2O/WH2O_*D0*Rkc/(D0 + Rkc);
scalar dmC = dmCa + dmCb +dmCc;
// Limit mass transfer by availability of C
dmC = min(mass*fComb, dmC);
// Change in O2 mass [kg]
const scalar dmO2 = dmCa/WC_*Sb_*WO2_;
// Change in CO2 mass [kg]
const scalar dmCOa = dmC+dmO2;
const scalar dmCO2 = dmCb/WC_*WCO2_;
const scalar dmCOb = dmCb + dmCO2;
const scalar dmH2O = dmCc/WC_*WH2O_;
const scalar dmCOc = dmCc/WC_*WCO_;
const scalar dmH2 = dmCc + dmH2O - dmCOc;
const scalar dmCO = dmCOa + dmCOb + dmCOc;
// Update local particle C mass
dMassSolid[CsLocalId_] += dmC;
// Update carrier O2 and CO2 mass
dMassSRCarrier[O2GlobalId_] -= dmO2;
dMassSRCarrier[CO2GlobalId_] -= dmCO2;
dMassSRCarrier[COGlobalId_] += dmCO;
dMassSRCarrier[H2GlobalId_] += dmH2;
dMassSRCarrier[H2OGlobalId_] -= dmH2O;
const scalar HC =
this->owner().composition().solids().properties()[CsLocalId_].Hf()
+ this->owner().composition().solids().properties()[CsLocalId_].cp()*T;
const scalar HcCO2 =
this->owner().mcCarrierThermo().speciesData()[CO2GlobalId_].H(T);
const scalar HcCO =
this->owner().mcCarrierThermo().speciesData()[COGlobalId_].H(T);
const scalar HcH2 =
this->owner().mcCarrierThermo().speciesData()[H2GlobalId_].H(T);
const scalar HcH2O =
this->owner().mcCarrierThermo().speciesData()[H2OGlobalId_].H(T);
const scalar HcO2 =
this->owner().mcCarrierThermo().speciesData()[O2GlobalId_].H(T);
// Heat of reaction [J]
return (dmC*HC + dmO2*HcO2 + HcCO2*dmCO2 + HcH2O*dmH2O - dmCO*HcCO - dmH2*HcH2);
}

Many Thanks.
NiA