COxidationIntrinsicRate.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 #include "mathematicalConstants.H"
28 
29 using namespace Foam::constant;
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 template<class CloudType>
35 (
36  const dictionary& dict,
37  CloudType& owner
38 )
39 :
40  SurfaceReactionModel<CloudType>(dict, owner, typeName),
41  Sb_(readScalar(this->coeffDict().lookup("Sb"))),
42  C1_(readScalar(this->coeffDict().lookup("C1"))),
43  rMean_(readScalar(this->coeffDict().lookup("rMean"))),
44  theta_(readScalar(this->coeffDict().lookup("theta"))),
45  Ai_(readScalar(this->coeffDict().lookup("Ai"))),
46  Ei_(readScalar(this->coeffDict().lookup("Ei"))),
47  Ag_(readScalar(this->coeffDict().lookup("Ag"))),
48  tau_(this->coeffDict().lookupOrDefault("tau", sqrt(2.0))),
49  CsLocalId_(-1),
50  O2GlobalId_(owner.composition().carrierId("O2")),
51  CO2GlobalId_(owner.composition().carrierId("CO2")),
52  WC_(0.0),
53  WO2_(0.0),
54  HcCO2_(0.0)
55 {
56  // Determine Cs ids
57  label idSolid = owner.composition().idSolid();
58  CsLocalId_ = owner.composition().localId(idSolid, "C");
59 
60  // Set local copies of thermo properties
61  WO2_ = owner.thermo().carrier().W(O2GlobalId_);
62  const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
63  WC_ = WCO2 - WO2_;
64 
65  HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
66 
67  if (Sb_ < 0)
68  {
70  (
71  "COxidationIntrinsicRate<CloudType>"
72  "("
73  "const dictionary&, "
74  "CloudType&"
75  ")"
76  ) << "Stoichiometry of reaction, Sb, must be greater than zero" << nl
77  << exit(FatalError);
78  }
79 
80  const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
81  const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
82  Info<< " C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
83 }
84 
85 
86 template<class CloudType>
88 (
90 )
91 :
93  Sb_(srm.Sb_),
94  C1_(srm.C1_),
95  rMean_(srm.rMean_),
96  theta_(srm.theta_),
97  Ai_(srm.Ai_),
98  Ei_(srm.Ei_),
99  Ag_(srm.Ag_),
100  tau_(srm.tau_),
101  CsLocalId_(srm.CsLocalId_),
102  O2GlobalId_(srm.O2GlobalId_),
103  CO2GlobalId_(srm.CO2GlobalId_),
104  WC_(srm.WC_),
105  WO2_(srm.WO2_),
106  HcCO2_(srm.HcCO2_)
107 {}
108 
109 
110 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
111 
112 template<class CloudType>
115 {}
116 
117 
118 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
119 
120 template<class CloudType>
122 (
123  const scalar dt,
124  const label cellI,
125  const scalar d,
126  const scalar T,
127  const scalar Tc,
128  const scalar pc,
129  const scalar rhoc,
130  const scalar mass,
131  const scalarField& YGas,
132  const scalarField& YLiquid,
133  const scalarField& YSolid,
134  const scalarField& YMixture,
135  const scalar N,
136  scalarField& dMassGas,
137  scalarField& dMassLiquid,
138  scalarField& dMassSolid,
139  scalarField& dMassSRCarrier
140 ) const
141 {
142  // Fraction of remaining combustible material
143  const label idSolid = CloudType::parcelType::SLD;
144  const scalar Ychar = YMixture[idSolid]*YSolid[CsLocalId_];
145 
146  // Surface combustion until combustible fraction is consumed
147  if (Ychar < SMALL)
148  {
149  return 0.0;
150  }
151 
152  const SLGThermo& thermo = this->owner().thermo();
153 
154  // Local mass fraction of O2 in the carrier phase []
155  const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
156 
157  // Quick exit if oxidant not present
158  if (YO2 < ROOTVSMALL)
159  {
160  return 0.0;
161  }
162 
163  // Diffusion rate coefficient [m2/s]
164  const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
165 
166  // Apparent density of pyrolysis char [kg/m3]
167  const scalar rhop = 6.0*mass/(constant::mathematical::pi*pow3(d));
168 
169  // Knusden diffusion coefficient [m2/s]
170  const scalar Dkn = 97.0*rMean_*sqrt(T/WO2_);
171 
172  // Effective diffusion [m2/s]
173  const scalar De = theta_/sqr(tau_)/(1.0/Dkn + 1/D0);
174 
175  // Cell carrier phase O2 species density [kg/m^3]
176  const scalar rhoO2 = rhoc*YO2;
177 
178  // Partial pressure O2 [Pa]
179  const scalar ppO2 = rhoO2/WO2_*RR*Tc;
180 
181  // Intrinsic reactivity [1/s]
182  const scalar ki = Ai_*exp(-Ei_/RR/T);
183 
184  // Thiele modulus []
185  const scalar phi =
186  max(0.5*d*sqrt(Sb_*rhop*Ag_*ki*ppO2/(De*rhoO2)), ROOTVSMALL);
187 
188  // Effectiveness factor []
189  const scalar eta = max(3.0*sqr(phi)*(phi/tanh(phi) - 1.0), 0.0);
190 
191  // Chemical rate [kmol/m2/s]
192  const scalar R = eta*d/6.0*rhop*Ag_*ki;
193 
194  // Particle surface area [m2]
195  const scalar Ap = constant::mathematical::pi*sqr(d);
196 
197  // Change in C mass [kg]
198  scalar dmC = Ap*rhoc*RR*Tc*YO2/WO2_*D0*R/(D0 + R)*dt;
199 
200  // Limit mass transfer by availability of C
201  dmC = min(mass*Ychar, dmC);
202 
203  // Molar consumption [kmol]
204  const scalar dOmega = dmC/WC_;
205 
206  // Change in O2 mass [kg]
207  const scalar dmO2 = dOmega*Sb_*WO2_;
208 
209  // Mass of newly created CO2 [kg]
210  const scalar dmCO2 = dOmega*(WC_ + Sb_*WO2_);
211 
212  // Update local particle C mass
213  dMassSolid[CsLocalId_] += dOmega*WC_;
214 
215  // Update carrier O2 and CO2 mass
216  dMassSRCarrier[O2GlobalId_] -= dmO2;
217  dMassSRCarrier[CO2GlobalId_] += dmCO2;
218 
219  const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
220 
221  // carrier sensible enthalpy exchange handled via change in mass
222 
223  // Heat of reaction [J]
224  return dmC*HsC - dmCO2*HcCO2_;
225 }
226 
227 
228 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensionedScalar pow3(const dimensionedScalar &ds)
const basicSpecieMixture & carrier() const
Return reference to the gaseous components.
Definition: SLGThermo.C:108
COxidationIntrinsicRate(const dictionary &dict, CloudType &owner)
Construct from dictionary.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Collection of constants.
virtual scalar calculate(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 scalar N, scalarField &dMassGas, scalarField &dMassLiquid, scalarField &dMassSolid, scalarField &dMassSRCarrier) const
Update surface reactions.
#define R(A, B, C, D, E, F, K, M)
Intrinsic char surface reaction mndel.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
dimensionedScalar exp(const dimensionedScalar &ds)
const PtrList< solidProperties > & properties() const
Return the solidProperties properties.
messageStream Info
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package, and provides:
Definition: SLGThermo.H:62
Templated surface reaction model class.
dictionary dict
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
const solidMixtureProperties & solids() const
Return reference to the global (additional) solids.
Definition: SLGThermo.C:140
dimensionedScalar tanh(const dimensionedScalar &ds)
surfaceScalarField & phi
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
const scalar RR
Universal gas constant (default in [J/(kmol K)])
PtrList< volScalarField > & Y()
Return the mass-fraction fields.
error FatalError
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
virtual ~COxidationIntrinsicRate()
Destructor.