FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 2008-2010 OpenCFD Ltd.
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 
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 :
48  (
49  dict,
50  owner,
51  typeName
52  ),
53  D0_(dimensionedScalar(this->coeffDict().lookup("D0")).value()),
54  rho0_(dimensionedScalar(this->coeffDict().lookup("rho0")).value()),
55  T0_(dimensionedScalar(this->coeffDict().lookup("T0")).value()),
56  Dn_(dimensionedScalar(this->coeffDict().lookup("Dn")).value()),
57  A_(dimensionedScalar(this->coeffDict().lookup("A")).value()),
58  E_(dimensionedScalar(this->coeffDict().lookup("E")).value()),
59  n_(dimensionedScalar(this->coeffDict().lookup("n")).value()),
60  WVol_(dimensionedScalar(this->coeffDict().lookup("WVol")).value()),
61  CsLocalId_(-1),
62  O2GlobalId_(owner.composition().globalCarrierId("O2")),
63  CO2GlobalId_(owner.composition().globalCarrierId("CO2")),
64  WC_(0.0),
65  WO2_(0.0)
66 {
67  // Determine Cs ids
68  label idSolid = owner.composition().idSolid();
69  CsLocalId_ = owner.composition().localId(idSolid, "C");
70 
71  // Set local copies of thermo properties
72  WO2_ = owner.mcCarrierThermo().speciesData()[O2GlobalId_].W();
73  scalar WCO2 = owner.mcCarrierThermo().speciesData()[CO2GlobalId_].W();
74  WC_ = WCO2 - WO2_;
75 }
76 
77 
78 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
79 
80 template<class CloudType>
82 {}
83 
84 
85 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
86 
87 template<class CloudType>
89 {
90  return true;
91 }
92 
93 
94 template<class CloudType>
96 (
97  const scalar dt,
98  const label cellI,
99  const scalar d,
100  const scalar T,
101  const scalar Tc,
102  const scalar pc,
103  const scalar rhoc,
104  const scalar mass,
105  const scalarField& YGas,
106  const scalarField& YLiquid,
107  const scalarField& YSolid,
108  const scalarField& YMixture,
109  const scalar N,
110  scalarField& dMassGas,
111  scalarField& dMassLiquid,
112  scalarField& dMassSolid,
113  scalarField& dMassSRCarrier
114 ) const
115 {
116  // Fraction of remaining combustible material
117  const label idSolid = CloudType::parcelType::SLD;
118  const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
119 
120  // Surface combustion until combustible fraction is consumed
121  if (fComb < SMALL)
122  {
123  return 0.0;
124  }
125 
126  // Cell carrier phase O2 species density [kg/m^3]
127  const scalar rhoO2 =
128  rhoc*this->owner().mcCarrierThermo().Y(O2GlobalId_)[cellI];
129 
130  if (rhoO2 < SMALL)
131  {
132  return 0.0;
133  }
134 
135  // Particle surface area [m^2]
136  const scalar Ap = mathematicalConstant::pi*sqr(d);
137 
138  // Calculate diffision constant at continuous phase temperature
139  // and density [m^2/s]
140  const scalar D = D0_*(rho0_/rhoc)*pow(Tc/T0_, Dn_);
141 
142  // Far field partial pressure O2 [Pa]
143  const scalar ppO2 = rhoO2/WO2_*specie::RR*Tc;
144 
145  // Total molar concentration of the carrier phase [kmol/m^3]
146  const scalar C = pc/(specie::RR*Tc);
147 
148  if (debug)
149  {
150  Pout<< "mass = " << mass << nl
151  << "fComb = " << fComb << nl
152  << "Ap = " << Ap << nl
153  << "dt = " << dt << nl
154  << "C = " << C << nl
155  << endl;
156  }
157 
158  // Molar reaction rate per unit surface area [kmol/(m^2.s)]
159  scalar qCsOld = 0;
160  scalar qCs = 1;
161 
162  const scalar qCsLim = mass*fComb/(WC_*Ap*dt);
163 
164  if (debug)
165  {
166  Pout << "qCsLim = " << qCsLim << endl;
167  }
168 
169  label iter = 0;
170  while ((mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
171  {
172  qCsOld = qCs;
173  const scalar PO2Surface = ppO2*exp(-(qCs + N)*d/(2*C*D));
174  qCs = A_*exp(-E_/(specie::RR*T))*pow(PO2Surface, n_);
175  qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
176  qCs = min(qCs, qCsLim);
177 
178  if (debug)
179  {
180  Pout<< "iter = " << iter
181  << ", qCsOld = " << qCsOld
182  << ", qCs = " << qCs
183  << nl << endl;
184  }
185 
186  iter++;
187  }
188 
189  if (iter > maxIters_)
190  {
191  WarningIn
192  (
193  "scalar Foam::COxidationMurphyShaddix<CloudType>::calculate(...)"
194  ) << "iter limit reached (" << maxIters_ << ")" << nl << endl;
195  }
196 
197  // Calculate the number of molar units reacted
198  scalar dOmega = qCs*Ap*dt;
199 
200  // Add to carrier phase mass transfer
201  dMassSRCarrier[O2GlobalId_] += -dOmega*WO2_;
202  dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + WO2_);
203 
204  // Add to particle mass transfer
205  dMassSolid[CsLocalId_] += dOmega*WC_;
206 
207  const scalar HC =
208  this->owner().composition().solids().properties()[CsLocalId_].Hf()
209  + this->owner().composition().solids().properties()[CsLocalId_].cp()*T;
210  const scalar HCO2 =
211  this->owner().mcCarrierThermo().speciesData()[CO2GlobalId_].H(T);
212  const scalar HO2 =
213  this->owner().mcCarrierThermo().speciesData()[O2GlobalId_].H(T);
214 
215  // Heat of reaction
216  return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);
217 }
218 
219 
220 // ************************ vim: set sw=4 sts=4 et: ************************ //