COxidationMurphyShaddix.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) 2011-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 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 template<class CloudType>
33 
34 template<class CloudType>
36 
37 
38 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
39 
40 template<class CloudType>
42 (
43  const dictionary& dict,
44  CloudType& owner
45 )
46 :
47  SurfaceReactionModel<CloudType>(dict, owner, typeName),
48  D0_(readScalar(this->coeffDict().lookup("D0"))),
49  rho0_(readScalar(this->coeffDict().lookup("rho0"))),
50  T0_(readScalar(this->coeffDict().lookup("T0"))),
51  Dn_(readScalar(this->coeffDict().lookup("Dn"))),
52  A_(readScalar(this->coeffDict().lookup("A"))),
53  E_(readScalar(this->coeffDict().lookup("E"))),
54  n_(readScalar(this->coeffDict().lookup("n"))),
55  WVol_(readScalar(this->coeffDict().lookup("WVol"))),
56  CsLocalId_(-1),
57  O2GlobalId_(owner.composition().carrierId("O2")),
58  CO2GlobalId_(owner.composition().carrierId("CO2")),
59  WC_(0.0),
60  WO2_(0.0),
61  HcCO2_(0.0)
62 {
63  // Determine Cs ids
64  label idSolid = owner.composition().idSolid();
65  CsLocalId_ = owner.composition().localId(idSolid, "C");
66 
67  // Set local copies of thermo properties
68  WO2_ = owner.thermo().carrier().W(O2GlobalId_);
69  const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
70  WC_ = WCO2 - WO2_;
71 
72  HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
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  D0_(srm.D0_),
88  rho0_(srm.rho0_),
89  T0_(srm.T0_),
90  Dn_(srm.Dn_),
91  A_(srm.A_),
92  E_(srm.E_),
93  n_(srm.n_),
94  WVol_(srm.WVol_),
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>
108 {}
109 
110 
111 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
112 
113 template<class CloudType>
115 (
116  const scalar dt,
117  const label celli,
118  const scalar d,
119  const scalar T,
120  const scalar Tc,
121  const scalar pc,
122  const scalar rhoc,
123  const scalar mass,
124  const scalarField& YGas,
125  const scalarField& YLiquid,
126  const scalarField& YSolid,
127  const scalarField& YMixture,
128  const scalar N,
129  scalarField& dMassGas,
130  scalarField& dMassLiquid,
131  scalarField& dMassSolid,
132  scalarField& dMassSRCarrier
133 ) const
134 {
135  // Fraction of remaining combustible material
136  const label idSolid = CloudType::parcelType::SLD;
137  const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
138 
139  // Surface combustion until combustible fraction is consumed
140  if (fComb < small)
141  {
142  return 0.0;
143  }
144 
145  const SLGThermo& thermo = this->owner().thermo();
146 
147  // Cell carrier phase O2 species density [kg/m^3]
148  const scalar rhoO2 = rhoc*thermo.carrier().Y(O2GlobalId_)[celli];
149 
150  if (rhoO2 < small)
151  {
152  return 0.0;
153  }
154 
155  // Particle surface area [m^2]
156  const scalar Ap = constant::mathematical::pi*sqr(d);
157 
158  // Calculate diffision constant at continuous phase temperature
159  // and density [m^2/s]
160  const scalar D = D0_*(rho0_/rhoc)*pow(Tc/T0_, Dn_);
161 
162  // Far field partial pressure O2 [Pa]
163  const scalar ppO2 = rhoO2/WO2_*RR*Tc;
164 
165  // Total molar concentration of the carrier phase [kmol/m^3]
166  const scalar C = pc/(RR*Tc);
167 
168  if (debug)
169  {
170  Pout<< "mass = " << mass << nl
171  << "fComb = " << fComb << nl
172  << "Ap = " << Ap << nl
173  << "dt = " << dt << nl
174  << "C = " << C << nl
175  << endl;
176  }
177 
178  // Molar reaction rate per unit surface area [kmol/(m^2.s)]
179  scalar qCsOld = 0;
180  scalar qCs = 1;
181 
182  const scalar qCsLim = mass*fComb/(WC_*Ap*dt);
183 
184  if (debug)
185  {
186  Pout<< "qCsLim = " << qCsLim << endl;
187  }
188 
189  label iter = 0;
190  while ((mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
191  {
192  qCsOld = qCs;
193  const scalar PO2Surface = ppO2*exp(-(qCs + N)*d/(2*C*D));
194  qCs = A_*exp(-E_/(RR*T))*pow(PO2Surface, n_);
195  qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
196  qCs = min(qCs, qCsLim);
197 
198  if (debug)
199  {
200  Pout<< "iter = " << iter
201  << ", qCsOld = " << qCsOld
202  << ", qCs = " << qCs
203  << nl << endl;
204  }
205 
206  iter++;
207  }
208 
209  if (iter > maxIters_)
210  {
212  << "iter limit reached (" << maxIters_ << ")" << nl << endl;
213  }
214 
215  // Calculate the number of molar units reacted
216  scalar dOmega = qCs*Ap*dt;
217 
218  // Add to carrier phase mass transfer
219  dMassSRCarrier[O2GlobalId_] += -dOmega*WO2_;
220  dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + WO2_);
221 
222  // Add to particle mass transfer
223  dMassSolid[CsLocalId_] += dOmega*WC_;
224 
225  const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
226 
227  // carrier sensible enthalpy exchange handled via change in mass
228 
229  // Heat of reaction [J]
230  return dOmega*(WC_*HsC - (WC_ + WO2_)*HcCO2_);
231 }
232 
233 
234 // ************************************************************************* //
#define readScalar
Definition: doubleScalar.C:38
dictionary dict
Graphite solid properties.
Definition: C.H:48
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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
const PtrList< solidProperties > & properties() const
Return the solidProperties properties.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
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.
Limited to C(s) + O2 -> CO2.
stressControl lookup("compactNormalStress") >> compactNormalStress
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:98
dimensionedScalar exp(const dimensionedScalar &ds)
COxidationMurphyShaddix(const dictionary &dict, CloudType &owner)
Construct from dictionary.
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package, and provides:
Definition: SLGThermo.H:62
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)
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const scalar RR
Universal gas constant (default in [J/(kmol K)])
virtual ~COxidationMurphyShaddix()
Destructor.
Templated surface reaction model class.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
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.