COxidationIntrinsicRate.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2014-2018 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  << "Stoichiometry of reaction, Sb, must be greater than zero" << nl
71  << exit(FatalError);
72  }
73 
74  const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
75  const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
76  Info<< " C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
77 }
78 
79 
80 template<class CloudType>
82 (
84 )
85 :
87  Sb_(srm.Sb_),
88  C1_(srm.C1_),
89  rMean_(srm.rMean_),
90  theta_(srm.theta_),
91  Ai_(srm.Ai_),
92  Ei_(srm.Ei_),
93  Ag_(srm.Ag_),
94  tau_(srm.tau_),
95  CsLocalId_(srm.CsLocalId_),
96  O2GlobalId_(srm.O2GlobalId_),
97  CO2GlobalId_(srm.CO2GlobalId_),
98  WC_(srm.WC_),
99  WO2_(srm.WO2_),
100  HcCO2_(srm.HcCO2_)
101 {}
102 
103 
104 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
105 
106 template<class CloudType>
109 {}
110 
111 
112 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
113 
114 template<class CloudType>
116 (
117  const scalar dt,
118  const label celli,
119  const scalar d,
120  const scalar T,
121  const scalar Tc,
122  const scalar pc,
123  const scalar rhoc,
124  const scalar mass,
125  const scalarField& YGas,
126  const scalarField& YLiquid,
127  const scalarField& YSolid,
128  const scalarField& YMixture,
129  const scalar N,
130  scalarField& dMassGas,
131  scalarField& dMassLiquid,
132  scalarField& dMassSolid,
133  scalarField& dMassSRCarrier
134 ) const
135 {
136  // Fraction of remaining combustible material
137  const label idSolid = CloudType::parcelType::SLD;
138  const scalar Ychar = YMixture[idSolid]*YSolid[CsLocalId_];
139 
140  // Surface combustion until combustible fraction is consumed
141  if (Ychar < small)
142  {
143  return 0.0;
144  }
145 
146  const SLGThermo& thermo = this->owner().thermo();
147 
148  // Local mass fraction of O2 in the carrier phase []
149  const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[celli];
150 
151  // Quick exit if oxidant not present
152  if (YO2 < rootVSmall)
153  {
154  return 0.0;
155  }
156 
157  // Diffusion rate coefficient [m2/s]
158  const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
159 
160  // Apparent density of pyrolysis char [kg/m3]
161  const scalar rhop = 6.0*mass/(constant::mathematical::pi*pow3(d));
162 
163  // Knusden diffusion coefficient [m2/s]
164  const scalar Dkn = 97.0*rMean_*sqrt(T/WO2_);
165 
166  // Effective diffusion [m2/s]
167  const scalar De = theta_/sqr(tau_)/(1.0/Dkn + 1/D0);
168 
169  // Cell carrier phase O2 species density [kg/m^3]
170  const scalar rhoO2 = rhoc*YO2;
171 
172  // Partial pressure O2 [Pa]
173  const scalar ppO2 = rhoO2/WO2_*RR*Tc;
174 
175  // Intrinsic reactivity [1/s]
176  const scalar ki = Ai_*exp(-Ei_/RR/T);
177 
178  // Thiele modulus []
179  const scalar phi =
180  max(0.5*d*sqrt(Sb_*rhop*Ag_*ki*ppO2/(De*rhoO2)), rootVSmall);
181 
182  // Effectiveness factor []
183  const scalar eta = max(3.0*sqr(phi)*(phi/tanh(phi) - 1.0), 0.0);
184 
185  // Chemical rate [kmol/m2/s]
186  const scalar R = eta*d/6.0*rhop*Ag_*ki;
187 
188  // Particle surface area [m2]
189  const scalar Ap = constant::mathematical::pi*sqr(d);
190 
191  // Change in C mass [kg]
192  scalar dmC = Ap*rhoc*RR*Tc*YO2/WO2_*D0*R/(D0 + R)*dt;
193 
194  // Limit mass transfer by availability of C
195  dmC = min(mass*Ychar, dmC);
196 
197  // Molar consumption [kmol]
198  const scalar dOmega = dmC/WC_;
199 
200  // Change in O2 mass [kg]
201  const scalar dmO2 = dOmega*Sb_*WO2_;
202 
203  // Mass of newly created CO2 [kg]
204  const scalar dmCO2 = dOmega*(WC_ + Sb_*WO2_);
205 
206  // Update local particle C mass
207  dMassSolid[CsLocalId_] += dOmega*WC_;
208 
209  // Update carrier O2 and CO2 mass
210  dMassSRCarrier[O2GlobalId_] -= dmO2;
211  dMassSRCarrier[CO2GlobalId_] += dmCO2;
212 
213  const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
214 
215  // carrier sensible enthalpy exchange handled via change in mass
216 
217  // Heat of reaction [J]
218  return dmC*HsC - dmCO2*HcCO2_;
219 }
220 
221 
222 // ************************************************************************* //
Collection of constants.
dimensionedScalar tanh(const dimensionedScalar &ds)
dictionary dict
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
surfaceScalarField & phi
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const PtrList< solidProperties > & properties() const
Return the solidProperties properties.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
rhoReactionThermo & thermo
Definition: createFields.H:28
const solidMixtureProperties & solids() const
Return reference to the global (additional) solids.
Definition: SLGThermo.C:134
PtrList< volScalarField > & Y()
Return the mass-fraction fields.
virtual ~COxidationIntrinsicRate()
Destructor.
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensionedScalar exp(const dimensionedScalar &ds)
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package, and provides:
Definition: SLGThermo.H:62
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
static const char nl
Definition: Ostream.H:265
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const basicSpecieMixture & carrier() const
Return reference to the gaseous components.
Definition: SLGThermo.C:108
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
COxidationIntrinsicRate(const dictionary &dict, CloudType &owner)
Construct from dictionary.
#define R(A, B, C, D, E, F, K, M)
Intrinsic char surface reaction mndel.
messageStream Info
const scalar RR
Universal gas constant (default in [J/(kmol K)])
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.
Templated surface reaction model class.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69