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-2011 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
  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>(dict, owner, typeName),
  48     D0_(readScalar(this->coeffDict().lookup("D0"))),
  49     rho0_(readScalar(this->coeffDict().lookup("rho0"))),
  50     T0_(readScalar(this->coeffDict().lookup("T0"))),
  51     Dn_(readScalar(this->coeffDict().lookup("Dn"))),
  52     A_(readScalar(this->coeffDict().lookup("A"))),
  53     E_(readScalar(this->coeffDict().lookup("E"))),
  54     n_(readScalar(this->coeffDict().lookup("n"))),
  55     WVol_(readScalar(this->coeffDict().lookup("WVol"))),
  56     CsLocalId_(-1),
  57     O2GlobalId_(owner.composition().globalCarrierId("O2")),
  58     CO2GlobalId_(owner.composition().globalCarrierId("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.thermo().carrier().W(O2GlobalId_);
  69     const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
  70     WC_ = WCO2 - WO2_;
  71
  72     HcCO2_ = owner.thermo().carrier().Hc(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>
  81 Foam::COxidationMurphyShaddix<CloudType>::COxidationMurphyShaddix
  82 (
  83     const COxidationMurphyShaddix<CloudType>& srm
  84 )
  85 :
  86     SurfaceReactionModel<CloudType>(srm),
  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 {}
 101
 102
 103 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 104
 105 template<class CloudType>
 106 Foam::COxidationMurphyShaddix<CloudType>::~COxidationMurphyShaddix()
 107 {}
 108
 109
 110 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 111
 112 template<class CloudType>
 113 Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
 114 (
 115     const scalar dt,
 116     const label cellI,
 117     const scalar d,
 118     const scalar T,
 119     const scalar Tc,
 120     const scalar pc,
 121     const scalar rhoc,
 122     const scalar mass,
 123     const scalarField& YGas,
 124     const scalarField& YLiquid,
 125     const scalarField& YSolid,
 126     const scalarField& YMixture,
 127     const scalar N,
 128     scalarField& dMassGas,
 129     scalarField& dMassLiquid,
 130     scalarField& dMassSolid,
 131     scalarField& dMassSRCarrier
 132 ) const
 133 {
 134     // Fraction of remaining combustible material
 135     const label idSolid = CloudType::parcelType::SLD;
 136     const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
 137
 138     // Surface combustion until combustible fraction is consumed
 139     if (fComb < SMALL)
 140     {
 141         return 0.0;
 142     }
 143
 144     const SLGThermo& thermo = this->owner().thermo();
 145
 146     // Cell carrier phase O2 species density [kg/m^3]
 147     const scalar rhoO2 = rhoc*thermo.carrier().Y(O2GlobalId_)[cellI];
 148
 149     if (rhoO2 < SMALL)
 150     {
 151         return 0.0;
 152     }
 153
 154     // Particle surface area [m^2]
 155     const scalar Ap = constant::mathematical::pi*sqr(d);
 156
 157     // Calculate diffision constant at continuous phase temperature
 158     // and density [m^2/s]
 159     const scalar D = D0_*(rho0_/rhoc)*pow(Tc/T0_, Dn_);
 160
 161     // Far field partial pressure O2 [Pa]
 162     const scalar ppO2 = rhoO2/WO2_*specie::RR*Tc;
 163
 164     // Total molar concentration of the carrier phase [kmol/m^3]
 165     const scalar C = pc/(specie::RR*Tc);
 166
 167     if (debug)
 168     {
 169         Pout<< "mass  = " << mass << nl
 170             << "fComb = " << fComb << nl
 171             << "Ap    = " << Ap << nl
 172             << "dt    = " << dt << nl
 173             << "C     = " << C << nl
 174             << endl;
 175     }
 176
 177     // Molar reaction rate per unit surface area [kmol/(m^2.s)]
 178     scalar qCsOld = 0;
 179     scalar qCs = 1;
 180
 181     const scalar qCsLim = mass*fComb/(WC_*Ap*dt);
 182
 183     if (debug)
 184     {
 185         Pout<< "qCsLim = " << qCsLim << endl;
 186     }
 187
 188     label iter = 0;
 189     while ((mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
 190     {
 191         qCsOld = qCs;
 192         const scalar PO2Surface = ppO2*exp(-(qCs + N)*d/(2*C*D));
 193         qCs = A_*exp(-E_/(specie::RR*T))*pow(PO2Surface, n_);
 194         qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
 195         qCs = min(qCs, qCsLim);
 196
 197         if (debug)
 198         {
 199             Pout<< "iter = " << iter
 200                 << ", qCsOld = " << qCsOld
 201                 << ", qCs = " << qCs
 202                 << nl << endl;
 203         }
 204
 205         iter++;
 206     }
 207
 208     if (iter > maxIters_)
 209     {
 210         WarningIn
 211         (
 212             "scalar Foam::COxidationMurphyShaddix<CloudType>::calculate(...)"
 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 // ************************************************************************* //