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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "COxidationMurphyShaddix.H"
  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>
  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 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 // ************************************************************************* //