61int main(
int argc,
char *argv[])
65 "Calculate the equilibrium flame temperature for a given fuel and"
66 " pressure for a range of unburnt gas temperatures and equivalence"
68 "Includes the effects of dissociation on O2, H2O and CO2."
71 argList::noParallel();
72 argList::noFunctionObjects();
73 argList::addArgument(
"controlFile");
80 IFstream controlFile(controlFileName);
83 if (!controlFile.good())
86 <<
"Cannot read file " << controlFileName
93 const scalar P(control.get<scalar>(
"P"));
94 const word fuelName(control.get<
word>(
"fuel"));
95 const scalar
n(control.get<scalar>(
"n"));
96 const scalar m(control.get<scalar>(
"m"));
99 Info<<
nl <<
"Reading thermodynamic data dictionary" <<
endl;
104 IFstream thermoDataFile(thermoDataFileName);
107 if (!thermoDataFile.good())
110 <<
"Cannot read file " << thermoDataFileName
111 <<
abort(FatalError);
117 Info<<
nl <<
"Reading thermodynamic data for relevant species"
121 thermo FUEL(thermoData.subDict(fuelName)); FUEL *= FUEL.W();
124 thermo O2(thermoData.subDict(
"O2")); O2 *= O2.W();
128 thermo H2(thermoData.subDict(
"H2")); H2 *= H2.W();
131 thermo CO2(thermoData.subDict(
"CO2")); CO2 *= CO2.W();
133 thermo CO(thermoData.subDict(
"CO")); CO *= CO.W();
150 scalar stoicO2 =
n + m/4.0;
151 scalar stoicN2 = (0.79/0.21)*(
n + m/4.0);
153 scalar stoicH2O = m/2.0;
165 "stoichiometricAirFuelMassRatio",
170 Info<<
"stoichiometricAirFuelMassRatio "
171 << stoichiometricAirFuelMassRatio <<
';' <<
endl;
173 Info<<
"Equilibrium flame temperature data ("
174 << P/1e5 <<
" bar)" <<
nl <<
nl
180 <<
setw(12) <<
"Terror"
181 <<
setw(20) <<
"O2res (mole frac)" <<
nl
186 for (
int i=0; i<16; i++)
188 scalar equiv = 0.6 + i*0.05;
189 scalar ft = 1/(1 + stoichiometricAirFuelMassRatio.value()/equiv);
192 for (
int j=0; j<28; j++)
194 scalar
T0 = 300.0 + j*100.0;
197 scalar o2 = (1.0/equiv)*stoicO2;
198 scalar n2 = (0.79/0.21)*o2;
199 scalar fres =
max(1.0 - 1.0/equiv, 0.0);
200 scalar fburnt = 1.0 - fres;
204 scalar oresInit =
max(1.0/equiv - 1.0, 0.0)*stoicO2;
205 scalar co2Init = fburnt*stoicCO2;
206 scalar h2oInit = fburnt*stoicH2O;
208 scalar ores = oresInit;
209 scalar co2 = co2Init;
210 scalar h2o = h2oInit;
216 scalar
N = fres + n2 + co2 + h2o + ores;
220 scalar adiabaticFlameTemperature =
222 + (fburnt/(1.0 + o2 + n2))/(1.0/(1.0 + (1.0 + 0.79/0.21)*stoicO2))
225 scalar equilibriumFlameTemperature = adiabaticFlameTemperature;
229 for (
int j=0; j<20; j++)
236 CO2BreakUp.Kn(P, equilibriumFlameTemperature,
N)
237 /::sqrt(
max(ores, 0.001)),
244 H2OBreakUp.Kn(P, equilibriumFlameTemperature,
N)
245 /::sqrt(
max(ores, 0.001)),
251 ores = oresInit + 0.5*co + 0.5*h2;
261 fres*FUEL + ores*O2 + n2*
N2
262 + co2*CO2 + h2o*
H2O + co*CO + h2*H2
266 scalar equilibriumFlameTemperatureNew =
267 products.THa(reactants.Ha(P,
T0), P, adiabaticFlameTemperature);
271 adiabaticFlameTemperature = equilibriumFlameTemperatureNew;
275 equilibriumFlameTemperature = 0.5*
277 equilibriumFlameTemperature
278 + equilibriumFlameTemperatureNew
286 <<
setw(12) << adiabaticFlameTemperature
287 <<
setw(12) << equilibriumFlameTemperature
289 << adiabaticFlameTemperature - equilibriumFlameTemperature
290 <<
setw(12) << ores/
N
Istream and Ostream manipulators taking arguments.
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
Input from file stream, using an ISstream.
Thermodynamics mapping class to expose the absolute enthalpy functions.
Extract command arguments and options from the supplied argc and argv parameters.
T get(const label index) const
Get a value from the argument at index.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
A class for handling file names.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
scalar W() const
Molecular weight [kg/kmol].
A class for handling words, derived from Foam::string.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Functions to search 'etc' directories for configuration files etc.
messageStream Info
Information stream (stdout output on master, null elsewhere)
fileName findEtcFile(const fileName &name, const bool mandatory=false, unsigned short location=0777)
Search for a single FILE within the etc directories.
Omanip< int > setw(const int i)
Ostream & endl(Ostream &os)
Add newline and flush stream.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
constexpr char nl
The newline '\n' character (0x0a)
Foam::argList args(argc, argv)
const Vector< label > N(dict.get< Vector< label > >("N"))