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-2021 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>
32 Foam::label Foam::COxidationMurphyShaddix<CloudType>::maxIters_ = 1000;
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_(this->coeffDict().template lookup<scalar>("D0")),
49  rho0_(this->coeffDict().template lookup<scalar>("rho0")),
50  T0_(this->coeffDict().template lookup<scalar>("T0")),
51  Dn_(this->coeffDict().template lookup<scalar>("Dn")),
52  A_(this->coeffDict().template lookup<scalar>("A")),
53  E_(this->coeffDict().template lookup<scalar>("E")),
54  n_(this->coeffDict().template lookup<scalar>("n")),
55  WVol_(this->coeffDict().template lookup<scalar>("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.composition().carrier().Wi(O2GlobalId_);
69  const scalar WCO2 = owner.composition().carrier().Wi(CO2GlobalId_);
70  WC_ = WCO2 - WO2_;
71 
72  HcCO2_ = owner.composition().carrier().Hf(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 parcelThermo& thermo = this->owner().thermo();
146  const basicSpecieMixture& carrier = this->owner().composition().carrier();
147 
148  // Cell carrier phase O2 species density [kg/m^3]
149  const scalar rhoO2 = rhoc*carrier.Y(O2GlobalId_)[celli];
150 
151  if (rhoO2 < small)
152  {
153  return 0.0;
154  }
155 
156  // Particle surface area [m^2]
157  const scalar Ap = constant::mathematical::pi*sqr(d);
158 
159  // Calculate diffusion constant at continuous phase temperature
160  // and density [m^2/s]
161  const scalar D = D0_*(rho0_/rhoc)*pow(Tc/T0_, Dn_);
162 
163  // Far field partial pressure O2 [Pa]
164  const scalar ppO2 = rhoO2/WO2_*RR*Tc;
165 
166  // Total molar concentration of the carrier phase [kmol/m^3]
167  const scalar C = pc/(RR*Tc);
168 
169  if (debug)
170  {
171  Pout<< "mass = " << mass << nl
172  << "fComb = " << fComb << nl
173  << "Ap = " << Ap << nl
174  << "dt = " << dt << nl
175  << "C = " << C << nl
176  << endl;
177  }
178 
179  // Molar reaction rate per unit surface area [kmol/m^2/s]
180  scalar qCsOld = 0;
181  scalar qCs = 1;
182 
183  const scalar qCsLim = mass*fComb/(WC_*Ap*dt);
184 
185  if (debug)
186  {
187  Pout<< "qCsLim = " << qCsLim << endl;
188  }
189 
190  label iter = 0;
191  while ((mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
192  {
193  qCsOld = qCs;
194  const scalar PO2Surface = ppO2*exp(-(qCs + N)*d/(2*C*D));
195  qCs = A_*exp(-E_/(RR*T))*pow(PO2Surface, n_);
196  qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
197  qCs = min(qCs, qCsLim);
198 
199  if (debug)
200  {
201  Pout<< "iter = " << iter
202  << ", qCsOld = " << qCsOld
203  << ", qCs = " << qCs
204  << nl << endl;
205  }
206 
207  iter++;
208  }
209 
210  if (iter > maxIters_)
211  {
213  << "iter limit reached (" << maxIters_ << ")" << nl << endl;
214  }
215 
216  // Calculate the number of molar units reacted
217  scalar dOmega = qCs*Ap*dt;
218 
219  // Add to carrier phase mass transfer
220  dMassSRCarrier[O2GlobalId_] += -dOmega*WO2_;
221  dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + WO2_);
222 
223  // Add to particle mass transfer
224  dMassSolid[CsLocalId_] += dOmega*WC_;
225 
226  const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
227 
228  // carrier sensible enthalpy exchange handled via change in mass
229 
230  // Heat of reaction [J]
231  return dOmega*(WC_*HsC - (WC_ + WO2_)*HcCO2_);
232 }
233 
234 
235 // ************************************************************************* //
dictionary dict
Graphite solid properties.
Definition: C.H:48
fluidReactionThermo & thermo
Definition: createFields.H:28
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
const PtrList< solidProperties > & properties() const
Return the solidProperties properties.
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package, and provides:
Definition: parcelThermo.H:58
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Specialisation of basicMixture for a mixture consisting of a number for molecular species...
Limited to C(s) + O2 -> CO2.
const dimensionedScalar RR
Universal gas constant: default SI units: [J/kmol/K].
dimensionedScalar exp(const dimensionedScalar &ds)
COxidationMurphyShaddix(const dictionary &dict, CloudType &owner)
Construct from dictionary.
const solidMixtureProperties & solids() const
Return reference to the global (additional) solids.
Definition: parcelThermo.C:95
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
volScalarField rhoc(IOobject(rhocValue.name(), runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, rhocValue)
static const char nl
Definition: Ostream.H:260
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 > &)
virtual ~COxidationMurphyShaddix()
Destructor.
PtrList< volScalarField > & Y()
Return the mass-fraction fields.
Templated surface reaction model class.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
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.