applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 Application
  25     adiabaticFlameT
  26
  27 Description
  28     Calculates the equilibrium flame temperature for a given fuel and
  29     pressure for a range of unburnt gas temperatures and equivalence
  30     ratios; the effects of dissociation on O2, H2O and CO2 are included.
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "argList.H"
  35 #include "Time.H"
  36 #include "dictionary.H"
  37 #include "IFstream.H"
  38 #include "OSspecific.H"
  39 #include "IOmanip.H"
  40
  41 #include "specieThermo.H"
  42 #include "janafThermo.H"
  43 #include "perfectGas.H"
  44
  45 using namespace Foam;
  46
  47 typedef specieThermo<janafThermo<perfectGas> > thermo;
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 int main(int argc, char *argv[])
  52 {
  53     argList::validArgs.append("controlFile");
  54     argList args(argc, argv);
  55
  56     const fileName controlFileName = args[1];
  57
  58     // Construct control dictionary
  59     IFstream controlFile(controlFileName);
  60
  61     // Check controlFile stream is OK
  62     if (!controlFile.good())
  63     {
  64         FatalErrorIn(args.executable())
  65             << "Cannot read file " << controlFileName
  66             << abort(FatalError);
  67     }
  68
  69     dictionary control(controlFile);
  70
  71
  72     scalar P(readScalar(control.lookup("P")));
  73     const word fuelName(control.lookup("fuel"));
  74     scalar n(readScalar(control.lookup("n")));
  75     scalar m(readScalar(control.lookup("m")));
  76
  77
  78     Info<< nl << "Reading Burcat data dictionary" << endl;
  79
  80     fileName BurcatCpDataFileName(findEtcFile("thermoData/BurcatCpData"));
  81
  82     // Construct control dictionary
  83     IFstream BurcatCpDataFile(BurcatCpDataFileName);
  84
  85     // Check BurcatCpData stream is OK
  86     if (!BurcatCpDataFile.good())
  87     {
  88         FatalErrorIn(args.executable())
  89             << "Cannot read file " << BurcatCpDataFileName
  90             << abort(FatalError);
  91     }
  92
  93     dictionary thermoData(BurcatCpDataFile);
  94
  95
  96     Info<< nl << "Reading Burcat data for relevant species" << nl << endl;
  97
  98     // Reactants
  99     thermo FUEL(thermoData.lookup(fuelName));
 100     thermo O2(thermoData.lookup("O2"));
 101     thermo N2(thermoData.lookup("N2"));
 102
 103     // Products
 104     thermo CO2(thermoData.lookup("CO2"));
 105     thermo H2O(thermoData.lookup("H2O"));
 106
 107     // Product fragments
 108     thermo CO(thermoData.lookup("CO"));
 109     thermo H2(thermoData.lookup("H2"));
 110
 111
 112     // Product dissociation reactions
 113
 114     thermo CO2BreakUp
 115     (
 116         CO2 == CO + 0.5* O2
 117     );
 118
 119     thermo H2OBreakUp
 120     (
 121         H2O == H2 + 0.5*O2
 122     );
 123
 124
 125     // Stoiciometric number of moles of species for one mole of fuel
 126     scalar stoicO2 = n + m/4.0;
 127     scalar stoicN2 = (0.79/0.21)*(n + m/4.0);
 128     scalar stoicCO2 = n;
 129     scalar stoicH2O = m/2.0;
 130
 131     // Oxidant
 132     thermo oxidant
 133     (
 134         "oxidant",
 135         stoicO2*O2
 136       + stoicN2*N2
 137     );
 138
 139     dimensionedScalar stoichiometricAirFuelMassRatio
 140     (
 141         "stoichiometricAirFuelMassRatio",
 142         dimless,
 143         (oxidant.W()*oxidant.nMoles())/FUEL.W()
 144     );
 145
 146     Info<< "stoichiometricAirFuelMassRatio "
 147         << stoichiometricAirFuelMassRatio << ';' << endl;
 148
 149     Info<< "Equilibrium flame temperature data ("
 150         << P/1e5 << " bar)" << nl << nl
 151         << setw(3) << "Phi"
 152         << setw(12) << "ft"
 153         << setw(7) << "T0"
 154         << setw(12) << "Tad"
 155         << setw(12) << "Teq"
 156         << setw(12) << "Terror"
 157         << setw(20) << "O2res (mole frac)" << nl
 158         << endl;
 159
 160
 161     // Loop over equivalence ratios
 162     for (int i=0; i<16; i++)
 163     {
 164         scalar equiv = 0.6 + i*0.05;
 165         scalar ft = 1/(1 + stoichiometricAirFuelMassRatio.value()/equiv);
 166
 167     // Loop over initial temperatures
 168     for (int j=0; j<28; j++)
 169     {
 170         scalar T0 = 300.0 + j*100.0;
 171
 172         // Number of moles of species for one mole of fuel
 173         scalar o2 = (1.0/equiv)*stoicO2;
 174         scalar n2 = (0.79/0.21)*o2;
 175         scalar fres = max(1.0 - 1.0/equiv, 0.0);
 176         scalar fburnt = 1.0 - fres;
 177
 178         // Initial guess for number of moles of product species
 179         // ignoring product dissociation
 180         scalar oresInit = max(1.0/equiv - 1.0, 0.0)*stoicO2;
 181         scalar co2Init = fburnt*stoicCO2;
 182         scalar h2oInit = fburnt*stoicH2O;
 183
 184         scalar ores = oresInit;
 185         scalar co2 = co2Init;
 186         scalar h2o = h2oInit;
 187
 188         scalar co = 0.0;
 189         scalar h2 = 0.0;
 190
 191         // Total number of moles in system
 192         scalar N = fres + n2 + co2 + h2o + ores;
 193
 194
 195         // Initial guess for adiabatic flame temperature
 196         scalar adiabaticFlameTemperature =
 197             T0
 198           + (fburnt/(1.0 + o2 + n2))/(1.0/(1.0 + (1.0 + 0.79/0.21)*stoicO2))
 199            *2000.0;
 200
 201         scalar equilibriumFlameTemperature = adiabaticFlameTemperature;
 202
 203
 204         // Iteration loop for adiabatic flame temperature
 205         for (int j=0; j<20; j++)
 206         {
 207
 208             if (j > 0)
 209             {
 210                 co = co2*
 211                     min
 212                     (
 213                         CO2BreakUp.Kn(equilibriumFlameTemperature, P, N)
 214                        /::sqrt(max(ores, 0.001)),
 215                         1.0
 216                     );
 217
 218                 h2 = h2o*
 219                     min
 220                     (
 221                         H2OBreakUp.Kn(equilibriumFlameTemperature, P, N)
 222                        /::sqrt(max(ores, 0.001)),
 223                         1.0
 224                     );
 225
 226                 co2 = co2Init - co;
 227                 h2o = h2oInit - h2;
 228                 ores = oresInit + 0.5*co + 0.5*h2;
 229             }
 230
 231             thermo reactants
 232             (
 233                 FUEL + o2*O2 + n2*N2
 234             );
 235
 236             thermo products
 237             (
 238                 fres*FUEL + ores*O2 + n2*N2
 239               + co2*CO2 + h2o*H2O + co*CO + h2*H2
 240             );
 241
 242
 243             scalar equilibriumFlameTemperatureNew =
 244                 products.TH(reactants.H(T0), adiabaticFlameTemperature);
 245
 246             if (j==0)
 247             {
 248                 adiabaticFlameTemperature = equilibriumFlameTemperatureNew;
 249             }
 250             else
 251             {
 252                 equilibriumFlameTemperature = 0.5*
 253                 (
 254                     equilibriumFlameTemperature
 255                   + equilibriumFlameTemperatureNew
 256                 );
 257             }
 258         }
 259
 260         Info<< setw(3) << equiv
 261             << setw(12) << ft
 262             << setw(7) << T0
 263             << setw(12) << adiabaticFlameTemperature
 264             << setw(12) << equilibriumFlameTemperature
 265             << setw(12)
 266             << adiabaticFlameTemperature - equilibriumFlameTemperature
 267             << setw(12) << ores/N
 268             << endl;
 269     }
 270     }
 271
 272     Info<< nl << "end" << endl;
 273
 274     return 0;
 275 }
 276
 277
 278 // ************************************************************************* //