src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2008-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "COxidationMurphyShaddix.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 template<class CloudType>
  32 Foam::label Foam::COxidationMurphyShaddix<CloudType>::maxIters_ = 1000;
  33
  34 template<class CloudType>
  35 Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::tolerance_ = 1e-06;
  36
  37
  38 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  39
  40 template<class CloudType>
  41 Foam::COxidationMurphyShaddix<CloudType>::COxidationMurphyShaddix
  42 (
  43     const dictionary& dict,
  44     CloudType& owner
  45 )
  46 :
  47     SurfaceReactionModel<CloudType>
  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>
  81 Foam::COxidationMurphyShaddix<CloudType>::~COxidationMurphyShaddix()
  82 {}
  83
  84
  85 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  86
  87 template<class CloudType>
  88 bool Foam::COxidationMurphyShaddix<CloudType>::active() const
  89 {
  90     return true;
  91 }
  92
  93
  94 template<class CloudType>
  95 Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
  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 scalarField& dMassVolatile,
 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     // Molar emission rate of volatiles per unit surface area
 146     const scalar qVol = sum(dMassVolatile)/(WVol_*Ap);
 147
 148     // Total molar concentration of the carrier phase [kmol/m^3]
 149     const scalar C = pc/(specie::RR*Tc);
 150
 151     if (debug)
 152     {
 153         Pout<< "mass  = " << mass << nl
 154             << "fComb = " << fComb << nl
 155             << "Ap    = " << Ap << nl
 156             << "dt    = " << dt << nl
 157             << "C     = " << C << nl
 158             << endl;
 159     }
 160
 161     // Molar reaction rate per unit surface area [kmol/(m^2.s)]
 162     scalar qCsOld = 0;
 163     scalar qCs = 1;
 164
 165     const scalar qCsLim = mass*fComb/(WC_*Ap*dt);
 166
 167     if (debug)
 168     {
 169         Pout << "qCsLim = " << qCsLim << endl;
 170     }
 171
 172     label iter = 0;
 173     while ((mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
 174     {
 175         qCsOld = qCs;
 176         const scalar PO2Surface = ppO2*exp(-(qCs + qVol)*d/(2*C*D));
 177         qCs = A_*exp(-E_/(specie::RR*T))*pow(PO2Surface, n_);
 178         qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
 179         qCs = min(qCs, qCsLim);
 180
 181         if (debug)
 182         {
 183             Pout<< "iter = " << iter
 184                 << ", qCsOld = " << qCsOld
 185                 << ", qCs = " << qCs
 186                 << nl << endl;
 187         }
 188
 189         iter++;
 190     }
 191
 192     if (iter > maxIters_)
 193     {
 194         WarningIn
 195         (
 196             "scalar Foam::COxidationMurphyShaddix<CloudType>::calculate(...)"
 197         )   << "iter limit reached (" << maxIters_ << ")" << nl << endl;
 198     }
 199
 200     // Calculate the number of molar units reacted
 201     scalar dOmega = qCs*Ap*dt;
 202
 203     // Add to carrier phase mass transfer
 204     dMassSRCarrier[O2GlobalId_] += -dOmega*WO2_;
 205     dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + WO2_);
 206
 207     // Add to particle mass transfer
 208     dMassSolid[CsLocalId_] += dOmega*WC_;
 209
 210     const scalar HC =
 211         this->owner().composition().solids().properties()[CsLocalId_].Hf()
 212       + this->owner().composition().solids().properties()[CsLocalId_].cp()*T;
 213     const scalar HCO2 =
 214         this->owner().mcCarrierThermo().speciesData()[CO2GlobalId_].H(T);
 215     const scalar HO2 =
 216         this->owner().mcCarrierThermo().speciesData()[O2GlobalId_].H(T);
 217
 218     // Heat of reaction
 219     return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);
 220 }
 221
 222
 223 // ************************************************************************* //