COxidationMurphyShaddix.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
30 #include "mathematicalConstants.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 template<class CloudType>
36 
37 template<class CloudType>
39 
40 
41 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
42 
43 template<class CloudType>
45 (
46  const dictionary& dict,
47  CloudType& owner
48 )
49 :
50  SurfaceReactionModel<CloudType>(dict, owner, typeName),
51  D0_(this->coeffDict().getScalar("D0")),
52  rho0_(this->coeffDict().getScalar("rho0")),
53  T0_(this->coeffDict().getScalar("T0")),
54  Dn_(this->coeffDict().getScalar("Dn")),
55  A_(this->coeffDict().getScalar("A")),
56  E_(this->coeffDict().getScalar("E")),
57  n_(this->coeffDict().getScalar("n")),
58  WVol_(this->coeffDict().getScalar("WVol")),
59  CsLocalId_(-1),
60  O2GlobalId_(owner.composition().carrierId("O2")),
61  CO2GlobalId_(owner.composition().carrierId("CO2")),
62  WC_(0.0),
63  WO2_(0.0),
64  HcCO2_(0.0)
65 {
66  // Determine Cs ids
67  label idSolid = owner.composition().idSolid();
68  CsLocalId_ = owner.composition().localId(idSolid, "C");
69 
70  // Set local copies of thermo properties
71  WO2_ = owner.thermo().carrier().W(O2GlobalId_);
72  const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
73  WC_ = WCO2 - WO2_;
74 
75  HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
76 
77  const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
78  const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
79  Info<< " C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
80 }
81 
82 
83 template<class CloudType>
85 (
87 )
88 :
90  D0_(srm.D0_),
91  rho0_(srm.rho0_),
92  T0_(srm.T0_),
93  Dn_(srm.Dn_),
94  A_(srm.A_),
95  E_(srm.E_),
96  n_(srm.n_),
97  WVol_(srm.WVol_),
98  CsLocalId_(srm.CsLocalId_),
99  O2GlobalId_(srm.O2GlobalId_),
100  CO2GlobalId_(srm.CO2GlobalId_),
101  WC_(srm.WC_),
102  WO2_(srm.WO2_),
103  HcCO2_(srm.HcCO2_)
104 {}
105 
106 
107 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
108 
109 template<class CloudType>
111 (
112  const scalar dt,
113  const scalar Re,
114  const scalar nu,
115  const label celli,
116  const scalar d,
117  const scalar T,
118  const scalar Tc,
119  const scalar pc,
120  const scalar rhoc,
121  const scalar mass,
122  const scalarField& YGas,
123  const scalarField& YLiquid,
124  const scalarField& YSolid,
125  const scalarField& YMixture,
126  const scalar N,
127  scalarField& dMassGas,
128  scalarField& dMassLiquid,
129  scalarField& dMassSolid,
130  scalarField& dMassSRCarrier
131 ) const
132 {
133  // Fraction of remaining combustible material
134  const label idSolid = CloudType::parcelType::SLD;
135  const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
136 
137  // Surface combustion until combustible fraction is consumed
138  if (fComb < SMALL)
139  {
140  return 0.0;
141  }
142 
143  const SLGThermo& thermo = this->owner().thermo();
144 
145  // Cell carrier phase O2 species density [kg/m^3]
146  const scalar rhoO2 = rhoc*thermo.carrier().Y(O2GlobalId_)[celli];
147 
148  if (rhoO2 < SMALL)
149  {
150  return 0.0;
151  }
152 
153  // Particle surface area [m^2]
154  const scalar Ap = constant::mathematical::pi*sqr(d);
155 
156  // Calculate diffusion constant at continuous phase temperature
157  // and density [m^2/s]
158  const scalar D = D0_*(rho0_/rhoc)*pow(Tc/T0_, Dn_);
159 
160  // Far field partial pressure O2 [Pa]
161  const scalar ppO2 = rhoO2/WO2_*RR*Tc;
162 
163  // Total molar concentration of the carrier phase [kmol/m^3]
164  const scalar C = pc/(RR*Tc);
165 
166  if (debug)
167  {
168  Pout<< "mass = " << mass << nl
169  << "fComb = " << fComb << nl
170  << "Ap = " << Ap << nl
171  << "dt = " << dt << nl
172  << "C = " << C << nl
173  << endl;
174  }
175 
176  // Molar reaction rate per unit surface area [kmol/(m^2.s)]
177  scalar qCsOld = 0;
178  scalar qCs = 1;
179 
180  const scalar qCsLim = mass*fComb/(WC_*Ap*dt);
181 
182  if (debug)
183  {
184  Pout<< "qCsLim = " << qCsLim << endl;
185  }
186 
187  label iter = 0;
188  while ((mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
189  {
190  qCsOld = qCs;
191  const scalar PO2Surface = ppO2*exp(-(qCs + N)*d/(2*C*D));
192  qCs = A_*exp(-E_/(RR*T))*pow(PO2Surface, n_);
193  qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
194  qCs = min(qCs, qCsLim);
195 
196  if (debug)
197  {
198  Pout<< "iter = " << iter
199  << ", qCsOld = " << qCsOld
200  << ", qCs = " << qCs
201  << nl << endl;
202  }
203 
204  iter++;
205  }
206 
207  if (iter > maxIters_)
208  {
210  << "iter limit reached (" << maxIters_ << ")" << nl << endl;
211  }
212 
213  // Calculate the number of molar units reacted
214  scalar dOmega = qCs*Ap*dt;
215 
216  // Add to carrier phase mass transfer
217  dMassSRCarrier[O2GlobalId_] += -dOmega*WO2_;
218  dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + WO2_);
219 
220  // Add to particle mass transfer
221  dMassSolid[CsLocalId_] += dOmega*WC_;
222 
223  const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
224 
225  // carrier sensible enthalpy exchange handled via change in mass
226 
227  // Heat of reaction [J]
228  return dOmega*(WC_*HsC - (WC_ + WO2_)*HcCO2_);
229 }
230 
231 
232 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::constant::thermodynamic::RR
const scalar RR
Universal gas constant: default in [J/(kmol K)].
Definition: thermodynamicConstants.C:46
COxidationMurphyShaddix.H
Foam::SurfaceReactionModel
Templated surface reaction model class.
Definition: ReactingMultiphaseCloud.H:61
mathematicalConstants.H
Foam::SLGThermo
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package,...
Definition: SLGThermo.H:64
thermo
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::exp
dimensionedScalar exp(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:261
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::COxidationMurphyShaddix::calculate
virtual scalar calculate(const scalar dt, const scalar Re, const scalar nu, 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.
Definition: COxidationMurphyShaddix.C:111
nu
volScalarField & nu
Definition: readMechanicalProperties.H:176
Foam::Field< scalar >
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::DSMCCloud
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:71
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
T
const volScalarField & T
Definition: createFieldRefs.H:2
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::COxidationMurphyShaddix::COxidationMurphyShaddix
COxidationMurphyShaddix(const dictionary &dict, CloudType &owner)
Construct from dictionary.
Definition: COxidationMurphyShaddix.C:45
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::Re
scalarField Re(const UList< complex > &cf)
Extract real component.
Definition: complexField.C:159
D
const dimensionedScalar & D
Definition: solveBulkSurfactant.H:4
N
const Vector< label > N(dict.get< Vector< label >>("N"))
Foam::COxidationMurphyShaddix
Limited to C(s) + O2 -> CO2.
Definition: COxidationMurphyShaddix.H:54
Foam::C
Graphite solid properties.
Definition: C.H:50
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328