equilibriumFlameT.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 equilibriumFlameT
29
30Group
31 grpThermophysicalUtilities
32
33Description
34 Calculate the equilibrium flame temperature for a given fuel and
35 pressure for a range of unburnt gas temperatures and equivalence
36 ratios. Includes the effects of dissociation on O2, H2O and CO2.
37
38\*---------------------------------------------------------------------------*/
39
40#include "argList.H"
41#include "Time.H"
42#include "dictionary.H"
43#include "IFstream.H"
44#include "OSspecific.H"
45#include "etcFiles.H"
46#include "IOmanip.H"
47
48#include "specie.H"
49#include "perfectGas.H"
50#include "thermo.H"
51#include "janafThermo.H"
52#include "absoluteEnthalpy.H"
53
54using namespace Foam;
55
57 thermo;
58
59// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
60
61int main(int argc, char *argv[])
62{
63 argList::addNote
64 (
65 "Calculate the equilibrium flame temperature for a given fuel and"
66 " pressure for a range of unburnt gas temperatures and equivalence"
67 " ratios.\n"
68 "Includes the effects of dissociation on O2, H2O and CO2."
69 );
70
71 argList::noParallel();
72 argList::noFunctionObjects(); // Never use function objects
73 argList::addArgument("controlFile");
74
75 argList args(argc, argv);
76
77 const auto controlFileName = args.get<fileName>(1);
78
79 // Construct control dictionary
80 IFstream controlFile(controlFileName);
81
82 // Check controlFile stream is OK
83 if (!controlFile.good())
84 {
86 << "Cannot read file " << controlFileName
87 << abort(FatalError);
88 }
89
90 dictionary control(controlFile);
91
92
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"));
97
98
99 Info<< nl << "Reading thermodynamic data dictionary" << endl;
100
101 fileName thermoDataFileName(findEtcFile("thermoData/thermoData"));
102
103 // Construct control dictionary
104 IFstream thermoDataFile(thermoDataFileName);
105
106 // Check thermoData stream is OK
107 if (!thermoDataFile.good())
108 {
110 << "Cannot read file " << thermoDataFileName
111 << abort(FatalError);
112 }
113
114 dictionary thermoData(thermoDataFile);
115
116
117 Info<< nl << "Reading thermodynamic data for relevant species"
118 << nl << endl;
119
120 // Reactants (mole-based)
121 thermo FUEL(thermoData.subDict(fuelName)); FUEL *= FUEL.W();
122
123 // Oxidant (mole-based)
124 thermo O2(thermoData.subDict("O2")); O2 *= O2.W();
125 thermo N2(thermoData.subDict("N2")); N2 *= N2.W();
126
127 // Intermediates (mole-based)
128 thermo H2(thermoData.subDict("H2")); H2 *= H2.W();
129
130 // Products (mole-based)
131 thermo CO2(thermoData.subDict("CO2")); CO2 *= CO2.W();
132 thermo H2O(thermoData.subDict("H2O")); H2O *= H2O.W();
133 thermo CO(thermoData.subDict("CO")); CO *= CO.W();
134
135
136 // Product dissociation reactions
137
138 thermo CO2BreakUp
139 (
140 CO2 == CO + 0.5*O2
141 );
142
143 thermo H2OBreakUp
144 (
145 H2O == H2 + 0.5*O2
146 );
147
148
149 // Stoiciometric number of moles of species for one mole of fuel
150 scalar stoicO2 = n + m/4.0;
151 scalar stoicN2 = (0.79/0.21)*(n + m/4.0);
152 scalar stoicCO2 = n;
153 scalar stoicH2O = m/2.0;
154
155 // Oxidant
156 thermo oxidant
157 (
158 "oxidant",
159 stoicO2*O2
160 + stoicN2*N2
161 );
162
163 dimensionedScalar stoichiometricAirFuelMassRatio
164 (
165 "stoichiometricAirFuelMassRatio",
166 dimless,
167 oxidant.Y()/FUEL.W()
168 );
169
170 Info<< "stoichiometricAirFuelMassRatio "
171 << stoichiometricAirFuelMassRatio << ';' << endl;
172
173 Info<< "Equilibrium flame temperature data ("
174 << P/1e5 << " bar)" << nl << nl
175 << setw(3) << "Phi"
176 << setw(12) << "ft"
177 << setw(7) << "T0"
178 << setw(12) << "Tad"
179 << setw(12) << "Teq"
180 << setw(12) << "Terror"
181 << setw(20) << "O2res (mole frac)" << nl
182 << endl;
183
184
185 // Loop over equivalence ratios
186 for (int i=0; i<16; i++)
187 {
188 scalar equiv = 0.6 + i*0.05;
189 scalar ft = 1/(1 + stoichiometricAirFuelMassRatio.value()/equiv);
190
191 // Loop over initial temperatures
192 for (int j=0; j<28; j++)
193 {
194 scalar T0 = 300.0 + j*100.0;
195
196 // Number of moles of species for one mole of fuel
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;
201
202 // Initial guess for number of moles of product species
203 // ignoring product dissociation
204 scalar oresInit = max(1.0/equiv - 1.0, 0.0)*stoicO2;
205 scalar co2Init = fburnt*stoicCO2;
206 scalar h2oInit = fburnt*stoicH2O;
207
208 scalar ores = oresInit;
209 scalar co2 = co2Init;
210 scalar h2o = h2oInit;
211
212 scalar co = 0.0;
213 scalar h2 = 0.0;
214
215 // Total number of moles in system
216 scalar N = fres + n2 + co2 + h2o + ores;
217
218
219 // Initial guess for adiabatic flame temperature
220 scalar adiabaticFlameTemperature =
221 T0
222 + (fburnt/(1.0 + o2 + n2))/(1.0/(1.0 + (1.0 + 0.79/0.21)*stoicO2))
223 *2000.0;
224
225 scalar equilibriumFlameTemperature = adiabaticFlameTemperature;
226
227
228 // Iteration loop for adiabatic flame temperature
229 for (int j=0; j<20; j++)
230 {
231 if (j > 0)
232 {
233 co = co2*
234 min
235 (
236 CO2BreakUp.Kn(P, equilibriumFlameTemperature, N)
237 /::sqrt(max(ores, 0.001)),
238 1.0
239 );
240
241 h2 = h2o*
242 min
243 (
244 H2OBreakUp.Kn(P, equilibriumFlameTemperature, N)
245 /::sqrt(max(ores, 0.001)),
246 1.0
247 );
248
249 co2 = co2Init - co;
250 h2o = h2oInit - h2;
251 ores = oresInit + 0.5*co + 0.5*h2;
252 }
253
254 thermo reactants
255 (
256 FUEL + o2*O2 + n2*N2
257 );
258
259 thermo products
260 (
261 fres*FUEL + ores*O2 + n2*N2
262 + co2*CO2 + h2o*H2O + co*CO + h2*H2
263 );
264
265
266 scalar equilibriumFlameTemperatureNew =
267 products.THa(reactants.Ha(P, T0), P, adiabaticFlameTemperature);
268
269 if (j==0)
270 {
271 adiabaticFlameTemperature = equilibriumFlameTemperatureNew;
272 }
273 else
274 {
275 equilibriumFlameTemperature = 0.5*
276 (
277 equilibriumFlameTemperature
278 + equilibriumFlameTemperatureNew
279 );
280 }
281 }
282
283 Info<< setw(3) << equiv
284 << setw(12) << ft
285 << setw(7) << T0
286 << setw(12) << adiabaticFlameTemperature
287 << setw(12) << equilibriumFlameTemperature
288 << setw(12)
289 << adiabaticFlameTemperature - equilibriumFlameTemperature
290 << setw(12) << ores/N
291 << endl;
292 }
293 }
294
295 Info<< nl << "End" << nl << endl;
296
297 return 0;
298}
299
300
301// ************************************************************************* //
Istream and Ostream manipulators taking arguments.
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
label n
Y[inertIndex] max(0.0)
water
Definition: H2O.H:62
Input from file stream, using an ISstream.
Definition: IFstream.H:57
Liquid N2.
Definition: N2.H:63
Thermodynamics mapping class to expose the absolute enthalpy functions.
Extract command arguments and options from the supplied argc and argv parameters.
Definition: argList.H:124
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
A class for handling file names.
Definition: fileName.H:76
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.
Definition: word.H:68
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Functions to search 'etc' directories for configuration files etc.
Namespace for OpenFOAM.
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.
Definition: etcFiles.C:446
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
scalar T0
Definition: createFields.H:22
const Vector< label > N(dict.get< Vector< label > >("N"))