1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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
29 Calculates the adiabatic flame temperature for a given fuel over a
30 range of unburnt temperatures and equivalence ratios.
32 \*---------------------------------------------------------------------------*/
36 #include "dictionary.H"
38 #include "OSspecific.H"
40 #include "specieThermo.H"
41 #include "janafThermo.H"
42 #include "perfectGas.H"
46 typedef specieThermo<janafThermo<perfectGas> > thermo;
49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51 int main(int argc, char *argv[])
53 argList::validArgs.clear();
54 argList::validArgs.append("controlFile");
55 argList args(argc, argv);
57 fileName controlFileName(args.additionalArgs()[0]);
59 // Construct control dictionary
60 IFstream controlFile(controlFileName);
62 // Check controlFile stream is OK
63 if (!controlFile.good())
65 FatalErrorIn(args.executable())
66 << "Cannot read file " << controlFileName
70 dictionary control(controlFile);
73 scalar T0(readScalar(control.lookup("T0")));
74 word fuelName(control.lookup("fuel"));
75 scalar n(readScalar(control.lookup("n")));
76 scalar m(readScalar(control.lookup("m")));
79 Info<< nl << "Reading Burcat data dictionary" << endl;
81 fileName BurcatCpDataFileName(dotFoam("thermoData/BurcatCpData"));
83 // Construct control dictionary
84 IFstream BurcatCpDataFile(BurcatCpDataFileName);
86 // Check BurcatCpData stream is OK
87 if (!BurcatCpDataFile.good())
89 FatalErrorIn(args.executable())
90 << "Cannot read file " << BurcatCpDataFileName
94 dictionary CpData(BurcatCpDataFile);
97 scalar stoicO2 = n + m/4.0;
98 scalar stoicN2 = (0.79/0.21)*(n + m/4.0);
100 scalar stoicH2O = m/2.0;
105 thermo(CpData.lookup(fuelName))
111 stoicO2*thermo(CpData.lookup("O2"))
112 + stoicN2*thermo(CpData.lookup("N2"))
115 dimensionedScalar stoichiometricAirFuelMassRatio
117 "stoichiometricAirFuelMassRatio",
119 (oxidant.W()*oxidant.nMoles())/fuel.W()
122 Info<< "stoichiometricAirFuelMassRatio "
123 << stoichiometricAirFuelMassRatio << ';' << endl;
125 for (int i=0; i<300; i++)
127 scalar equiv = (i + 1)*0.01;
128 scalar ft = 1/(1 + stoichiometricAirFuelMassRatio.value()/equiv);
130 Info<< "phi = " << equiv << nl
131 << "ft = " << ft << endl;
133 scalar o2 = (1.0/equiv)*stoicO2;
134 scalar n2 = (0.79/0.21)*o2;
135 scalar fres = max(1.0 - 1.0/equiv, 0.0);
136 scalar ores = max(1.0/equiv - 1.0, 0.0);
137 scalar fburnt = 1.0 - fres;
142 thermo(CpData.lookup(fuelName))
144 Info<< "fuel " << fuel << ';' << endl;
149 o2*thermo(CpData.lookup("O2"))
150 + n2*thermo(CpData.lookup("N2"))
152 Info<< "oxidant " << (1/oxidant.nMoles())*oxidant << ';' << endl;
159 Info<< "reactants " << (1/reactants.nMoles())*reactants << ';' << endl;
164 + (n2 - (0.79/0.21)*ores*stoicO2)*thermo(CpData.lookup("N2"))
165 + fburnt*stoicCO2*thermo(CpData.lookup("CO2"))
166 + fburnt*stoicH2O*thermo(CpData.lookup("H2O"))
168 Info<< "burntProducts "
169 << (1/burntProducts.nMoles())*burntProducts << ';' << endl;
175 + n2*thermo(CpData.lookup("N2"))
176 + fburnt*stoicCO2*thermo(CpData.lookup("CO2"))
177 + fburnt*stoicH2O*thermo(CpData.lookup("H2O"))
178 + ores*stoicO2*thermo(CpData.lookup("O2"))
181 Info<< "products " << (1/products.nMoles())*products << ';' << endl;
183 scalar Tad = products.TH(reactants.H(T0), 1000.0);
184 Info<< "Tad = " << Tad << nl << endl;
187 Info<< nl << "end" << endl;
193 // ************************************************************************* //