applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Application
  26     adiabaticFlameT
  27
  28 Description
  29     Calculates the adiabatic flame temperature for a given fuel over a
  30     range of unburnt temperatures and equivalence ratios.
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "argList.H"
  35 #include "objectRegistry.H"
  36 #include "Time.H"
  37 #include "dictionary.H"
  38 #include "IFstream.H"
  39 #include "OSspecific.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
  52 int main(int argc, char *argv[])
  53 {
  54     argList::validArgs.clear();
  55     argList::validArgs.append("controlFile");
  56     argList args(argc, argv);
  57
  58     fileName controlFileName(args.additionalArgs()[0]);
  59
  60     // Construct control dictionary
  61     IFstream controlFile(controlFileName);
  62
  63     // Check controlFile stream is OK
  64     if (!controlFile.good())
  65     {
  66         FatalErrorIn(args.executable())
  67             << "Cannot read file " << controlFileName
  68             << exit(FatalError);
  69     }
  70
  71     dictionary control(controlFile);
  72
  73
  74     scalar T0(readScalar(control.lookup("T0")));
  75     word fuelName(control.lookup("fuel"));
  76     scalar n(readScalar(control.lookup("n")));
  77     scalar m(readScalar(control.lookup("m")));
  78
  79
  80     Info<< nl << "Reading Burcat data dictionary" << endl;
  81
  82     fileName BurcatCpDataFileName(findEtcFile("thermoData/BurcatCpData"));
  83
  84     // Construct control dictionary
  85     IFstream BurcatCpDataFile(BurcatCpDataFileName);
  86
  87     // Check BurcatCpData stream is OK
  88     if (!BurcatCpDataFile.good())
  89     {
  90         FatalErrorIn(args.executable())
  91             << "Cannot read file " << BurcatCpDataFileName
  92             << exit(FatalError);
  93     }
  94
  95     dictionary CpData(BurcatCpDataFile);
  96
  97
  98     scalar stoicO2 = n + m/4.0;
  99     scalar stoicN2 = (0.79/0.21)*(n + m/4.0);
 100     scalar stoicCO2 = n;
 101     scalar stoicH2O = m/2.0;
 102
 103     thermo fuel
 104     (
 105         "fuel",
 106         thermo(CpData.lookup(fuelName))
 107     );
 108
 109     thermo oxidant
 110     (
 111         "oxidant",
 112         stoicO2*thermo(CpData.lookup("O2"))
 113       + stoicN2*thermo(CpData.lookup("N2"))
 114     );
 115
 116     dimensionedScalar stoichiometricAirFuelMassRatio
 117     (
 118         "stoichiometricAirFuelMassRatio",
 119         dimless,
 120         (oxidant.W()*oxidant.nMoles())/fuel.W()
 121     );
 122
 123     Info<< "stoichiometricAirFuelMassRatio "
 124         << stoichiometricAirFuelMassRatio << ';' << endl;
 125
 126     for (int i=0; i<300; i++)
 127     {
 128         scalar equiv = (i + 1)*0.01;
 129         scalar ft = 1/(1 + stoichiometricAirFuelMassRatio.value()/equiv);
 130
 131         Info<< "phi = " << equiv << nl
 132             << "ft = " << ft << endl;
 133
 134         scalar o2 = (1.0/equiv)*stoicO2;
 135         scalar n2 = (0.79/0.21)*o2;
 136         scalar fres = max(1.0 - 1.0/equiv, 0.0);
 137         scalar ores = max(1.0/equiv - 1.0, 0.0);
 138         scalar fburnt = 1.0 - fres;
 139
 140         thermo fuel
 141         (
 142             "fuel",
 143             thermo(CpData.lookup(fuelName))
 144         );
 145         Info<< "fuel " << fuel << ';' << endl;
 146
 147         thermo oxidant
 148         (
 149             "oxidant",
 150             o2*thermo(CpData.lookup("O2"))
 151           + n2*thermo(CpData.lookup("N2"))
 152         );
 153         Info<< "oxidant " << (1/oxidant.nMoles())*oxidant << ';' << endl;
 154
 155         thermo reactants
 156         (
 157             "reactants",
 158             fuel + oxidant
 159         );
 160         Info<< "reactants " << (1/reactants.nMoles())*reactants << ';' << endl;
 161
 162         thermo burntProducts
 163         (
 164             "burntProducts",
 165           + (n2 - (0.79/0.21)*ores*stoicO2)*thermo(CpData.lookup("N2"))
 166           + fburnt*stoicCO2*thermo(CpData.lookup("CO2"))
 167           + fburnt*stoicH2O*thermo(CpData.lookup("H2O"))
 168         );
 169         Info<< "burntProducts "
 170             << (1/burntProducts.nMoles())*burntProducts << ';' << endl;
 171
 172         thermo products
 173         (
 174             "products",
 175             fres*fuel
 176           + n2*thermo(CpData.lookup("N2"))
 177           + fburnt*stoicCO2*thermo(CpData.lookup("CO2"))
 178           + fburnt*stoicH2O*thermo(CpData.lookup("H2O"))
 179           + ores*stoicO2*thermo(CpData.lookup("O2"))
 180         );
 181
 182         Info<< "products " << (1/products.nMoles())*products << ';' << endl;
 183
 184         scalar Tad = products.TH(reactants.H(T0), 1000.0);
 185         Info<< "Tad = " << Tad << nl << endl;
 186     }
 187
 188     Info<< nl << "end" << endl;
 189
 190     return 0;
 191 }
 192
 193
 194 // ************************************************************************* //