src/thermophysicalModels/reactionThermo/mixtures/singleStepReactingMixture/singleStepReactingMixture.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2010-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 "singleStepReactingMixture.H"
  27 #include "fvMesh.H"
  28
  29 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
  30
  31 template<class ThermoType>
  32 void Foam::singleStepReactingMixture<ThermoType>::calculateqFuel()
  33 {
  34     const Reaction<ThermoType>& reaction = this->operator[](0);
  35     const  scalar Wu = this->speciesData()[fuelIndex_].W();
  36
  37     forAll(reaction.lhs(), i)
  38     {
  39         const label specieI = reaction.lhs()[i].index;
  40         const scalar stoichCoeff = reaction.lhs()[i].stoichCoeff;
  41         specieStoichCoeffs_[specieI] = -stoichCoeff;
  42         qFuel_.value() += this->speciesData()[specieI].hc()*stoichCoeff/Wu;
  43     }
  44
  45     forAll(reaction.rhs(), i)
  46     {
  47         const label specieI = reaction.rhs()[i].index;
  48         const scalar stoichCoeff = reaction.rhs()[i].stoichCoeff;
  49         specieStoichCoeffs_[specieI] = stoichCoeff;
  50         qFuel_.value() -= this->speciesData()[specieI].hc()*stoichCoeff/Wu;
  51         specieProd_[specieI] = -1;
  52     }
  53
  54     Info << "Fuel heat of combustion :" << qFuel_.value() << endl;
  55 }
  56
  57
  58 template<class ThermoType>
  59 void Foam::singleStepReactingMixture<ThermoType>::massAndAirStoichRatios()
  60 {
  61     const label O2Index = this->species()["O2"];
  62     const scalar Wu = this->speciesData()[fuelIndex_].W();
  63
  64     stoicRatio_ =
  65        (this->speciesData()[inertIndex_].W()
  66       * specieStoichCoeffs_[inertIndex_]
  67       + this->speciesData()[O2Index].W()
  68       * mag(specieStoichCoeffs_[O2Index]))
  69       / (Wu*mag(specieStoichCoeffs_[fuelIndex_]));
  70
  71     s_ =
  72         (this->speciesData()[O2Index].W()
  73       * mag(specieStoichCoeffs_[O2Index]))
  74       / (Wu*mag(specieStoichCoeffs_[fuelIndex_]));
  75
  76     Info << "stoichiometric air-fuel ratio :" << stoicRatio_.value() << endl;
  77
  78     Info << "stoichiometric oxygen-fuel ratio :" << s_.value() << endl;
  79 }
  80
  81
  82 template<class ThermoType>
  83 void Foam::singleStepReactingMixture<ThermoType>::calculateMaxProducts()
  84 {
  85     const Reaction<ThermoType>& reaction = this->operator[](0);
  86
  87     scalar Wm = 0.0;
  88     scalar totalMol = 0.0;
  89     forAll(reaction.rhs(), i)
  90     {
  91         label specieI = reaction.rhs()[i].index;
  92         totalMol += mag(specieStoichCoeffs_[specieI]);
  93     }
  94
  95     scalarList Xi(reaction.rhs().size());
  96
  97     forAll(reaction.rhs(), i)
  98     {
  99         const label specieI = reaction.rhs()[i].index;
 100         Xi[i] = mag(specieStoichCoeffs_[specieI])/totalMol;
 101
 102         Wm += Xi[i]*this->speciesData()[specieI].W();
 103     }
 104
 105     forAll(reaction.rhs(), i)
 106     {
 107         const label specieI = reaction.rhs()[i].index;
 108         Yprod0_[specieI] =  this->speciesData()[specieI].W()/Wm*Xi[i];
 109     }
 110
 111     // Normalize the stoichiometric coeff to mass
 112     forAll(specieStoichCoeffs_, i)
 113     {
 114         specieStoichCoeffs_[i] =
 115             specieStoichCoeffs_[i]
 116           * this->speciesData()[i].W()
 117           / (this->speciesData()[fuelIndex_].W()
 118           * mag(specieStoichCoeffs_[fuelIndex_]));
 119     }
 120 }
 121
 122
 123 template<class ThermoType>
 124 void Foam::singleStepReactingMixture<ThermoType>::fresCorrect()
 125 {
 126     const Reaction<ThermoType>& reaction = this->operator[](0);
 127
 128     label O2Index = this->species()["O2"];
 129     const volScalarField& YFuel = this->Y()[fuelIndex_];
 130     const volScalarField& YO2 = this->Y()[O2Index];
 131
 132     // reactants
 133     forAll(reaction.lhs(), i)
 134     {
 135         const label specieI = reaction.lhs()[i].index;
 136         if (specieI == fuelIndex_)
 137         {
 138             fres_[specieI] =  max(YFuel - YO2/s_, scalar(0.0));
 139         }
 140         else if (specieI == O2Index)
 141         {
 142             fres_[specieI] =  max(YO2 - YFuel*s_, scalar(0.0));
 143         }
 144     }
 145
 146
 147     // products
 148     forAll(reaction.rhs(), i)
 149     {
 150         const label specieI = reaction.rhs()[i].index;
 151         if (specieI != inertIndex_)
 152         {
 153             forAll(fres_[specieI], cellI)
 154             {
 155                 if (fres_[fuelIndex_][cellI] > 0.0)
 156                 {
 157                     // rich mixture
 158                     fres_[specieI][cellI] =
 159                         Yprod0_[specieI]
 160                       * (1.0 + YO2[cellI]/s_.value() - YFuel[cellI]);
 161                 }
 162                 else
 163                 {
 164                     // lean mixture
 165                     fres_[specieI][cellI] =
 166                         Yprod0_[specieI]
 167                       * (
 168                             1.0
 169                           - YO2[cellI]/s_.value()*stoicRatio_.value()
 170                           + YFuel[cellI]*stoicRatio_.value()
 171                         );
 172                 }
 173             }
 174         }
 175     }
 176 }
 177
 178
 179 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 180
 181 template<class ThermoType>
 182 Foam::singleStepReactingMixture<ThermoType>::singleStepReactingMixture
 183 (
 184     const dictionary& thermoDict,
 185     const fvMesh& mesh
 186 )
 187 :
 188     reactingMixture<ThermoType>(thermoDict, mesh),
 189     stoicRatio_(dimensionedScalar("stoicRatio", dimless, 0.0)),
 190     s_(dimensionedScalar("s", dimless, 0.0)),
 191     qFuel_(dimensionedScalar("qFuel", sqr(dimVelocity), 0.0)),
 192     specieStoichCoeffs_(this->species_.size(), 0.0),
 193     Yprod0_(this->species_.size(), 0.0),
 194     fres_(Yprod0_.size()),
 195     inertIndex_(this->species()[thermoDict.lookup("inertSpecie")]),
 196     fuelIndex_(this->species()[thermoDict.lookup("fuel")]),
 197     specieProd_(Yprod0_.size(), 1)
 198 {
 199     if (this->size() == 1)
 200     {
 201         forAll(fres_, fresI)
 202         {
 203             IOobject header
 204             (
 205                 "fres_" + this->species()[fresI],
 206                 mesh.time().timeName(),
 207                 mesh,
 208                 IOobject::NO_READ,
 209                 IOobject::NO_WRITE
 210             );
 211
 212             fres_.set
 213             (
 214                 fresI,
 215                 new volScalarField
 216                 (
 217                     header,
 218                     mesh,
 219                     dimensionedScalar("fres" + name(fresI), dimless, 0.0)
 220                 )
 221             );
 222         }
 223
 224         calculateqFuel();
 225
 226         massAndAirStoichRatios();
 227
 228         calculateMaxProducts();
 229
 230         autoPtr<chemistryReader<ThermoType> >::clear();
 231     }
 232     else
 233     {
 234         FatalErrorIn
 235         (
 236             "singleStepReactingMixture::<ThermoType>::singleStepReactingMixture"
 237             "("
 238                 "const dictionary&, "
 239                 "const fvMesh&"
 240             ")"
 241         )   << "Only one reaction required for single step reaction"
 242             << exit(FatalError);
 243     }
 244 }
 245
 246
 247 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 248
 249 template<class ThermoType>
 250 void Foam::singleStepReactingMixture<ThermoType>::read
 251 (
 252     const dictionary& thermoDict
 253 )
 254 {}
 255
 256
 257 // ************************************************************************* //