All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
equilibriumFlameT.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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  equilibriumFlameT
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 "dictionary.H"
36 #include "IFstream.H"
37 #include "etcFiles.H"
38 #include "IOmanip.H"
39 #include "dimensionedTypes.H"
40 
41 #include "specie.H"
42 #include "perfectGas.H"
43 #include "thermo.H"
44 #include "janafThermo.H"
45 #include "absoluteEnthalpy.H"
46 
47 using namespace Foam;
48 
50  thermo;
51 
52 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
53 
54 int main(int argc, char *argv[])
55 {
56  #include "removeCaseOptions.H"
57 
58  argList::validArgs.append("properties dictionary");
59  argList args(argc, argv);
60 
61  const fileName propertiesDictName = args[1];
62 
63  // Construct control dictionary
64  IFstream propertiesDict(propertiesDictName);
65 
66  // Check propertiesDict stream is OK
67  if (!propertiesDict.good())
68  {
70  << "Cannot read file " << propertiesDictName
71  << abort(FatalError);
72  }
73 
74  dictionary control(propertiesDict);
75 
76 
77  scalar P(control.lookup<scalar>("P"));
78  const word fuelName(control.lookup("fuel"));
79  scalar n(control.lookup<scalar>("n"));
80  scalar m(control.lookup<scalar>("m"));
81 
82 
83  Info<< nl << "Reading thermodynamic data dictionary" << endl;
84 
85  fileName thermoDataFileName(findEtcFile("thermoData/thermoData"));
86 
87  // Construct control dictionary
88  IFstream thermoDataFile(thermoDataFileName);
89 
90  // Check thermoData stream is OK
91  if (!thermoDataFile.good())
92  {
94  << "Cannot read file " << thermoDataFileName
95  << abort(FatalError);
96  }
97 
98  dictionary thermoData(thermoDataFile);
99 
100 
101  Info<< nl << "Reading thermodynamic data for relevant species"
102  << nl << endl;
103 
104  // Reactants (mole-based)
105  thermo FUEL(thermoData.subDict(fuelName)); FUEL *= FUEL.W();
106 
107  // Oxidant (mole-based)
108  thermo O2(thermoData.subDict("O2")); O2 *= O2.W();
109  thermo N2(thermoData.subDict("N2")); N2 *= N2.W();
110 
111  // Intermediates (mole-based)
112  thermo H2(thermoData.subDict("H2")); H2 *= H2.W();
113 
114  // Products (mole-based)
115  thermo CO2(thermoData.subDict("CO2")); CO2 *= CO2.W();
116  thermo H2O(thermoData.subDict("H2O")); H2O *= H2O.W();
117  thermo CO(thermoData.subDict("CO")); CO *= CO.W();
118 
119 
120  // Product dissociation reactions
121 
122  thermo CO2BreakUp
123  (
124  CO2 == CO + 0.5*O2
125  );
126 
127  thermo H2OBreakUp
128  (
129  H2O == H2 + 0.5*O2
130  );
131 
132 
133  // Stoichiometric number of moles of species for one mole of fuel
134  scalar stoicO2 = n + m/4.0;
135  scalar stoicN2 = (0.79/0.21)*(n + m/4.0);
136  scalar stoicCO2 = n;
137  scalar stoicH2O = m/2.0;
138 
139  // Oxidant
140  thermo oxidant
141  (
142  "oxidant",
143  stoicO2*O2
144  + stoicN2*N2
145  );
146 
147  dimensionedScalar stoichiometricAirFuelMassRatio
148  (
149  "stoichiometricAirFuelMassRatio",
150  dimless,
151  oxidant.Y()/FUEL.W()
152  );
153 
154  Info<< "stoichiometricAirFuelMassRatio "
155  << stoichiometricAirFuelMassRatio << ';' << endl;
156 
157  Info<< "Equilibrium flame temperature data ("
158  << P/1e5 << " bar)" << nl << nl
159  << setw(3) << "Phi"
160  << setw(12) << "ft"
161  << setw(7) << "T0"
162  << setw(12) << "Tad"
163  << setw(12) << "Teq"
164  << setw(12) << "Terror"
165  << setw(20) << "O2res (mole frac)" << nl
166  << endl;
167 
168 
169  // Loop over equivalence ratios
170  for (int i=0; i<16; i++)
171  {
172  scalar equiv = 0.6 + i*0.05;
173  scalar ft = 1/(1 + stoichiometricAirFuelMassRatio.value()/equiv);
174 
175  // Loop over initial temperatures
176  for (int j=0; j<28; j++)
177  {
178  scalar T0 = 300.0 + j*100.0;
179 
180  // Number of moles of species for one mole of fuel
181  scalar o2 = (1.0/equiv)*stoicO2;
182  scalar n2 = (0.79/0.21)*o2;
183  scalar fres = max(1.0 - 1.0/equiv, 0.0);
184  scalar fburnt = 1.0 - fres;
185 
186  // Initial guess for number of moles of product species
187  // ignoring product dissociation
188  scalar oresInit = max(1.0/equiv - 1.0, 0.0)*stoicO2;
189  scalar co2Init = fburnt*stoicCO2;
190  scalar h2oInit = fburnt*stoicH2O;
191 
192  scalar ores = oresInit;
193  scalar co2 = co2Init;
194  scalar h2o = h2oInit;
195 
196  scalar co = 0.0;
197  scalar h2 = 0.0;
198 
199  // Total number of moles in system
200  scalar N = fres + n2 + co2 + h2o + ores;
201 
202 
203  // Initial guess for adiabatic flame temperature
204  scalar adiabaticFlameTemperature =
205  T0
206  + (fburnt/(1.0 + o2 + n2))/(1.0/(1.0 + (1.0 + 0.79/0.21)*stoicO2))
207  *2000.0;
208 
209  scalar equilibriumFlameTemperature = adiabaticFlameTemperature;
210 
211 
212  // Iteration loop for adiabatic flame temperature
213  for (int j=0; j<20; j++)
214  {
215  if (j > 0)
216  {
217  co = co2*
218  min
219  (
220  CO2BreakUp.Kn(P, equilibriumFlameTemperature, N)
221  /::sqrt(max(ores, 0.001)),
222  1.0
223  );
224 
225  h2 = h2o*
226  min
227  (
228  H2OBreakUp.Kn(P, equilibriumFlameTemperature, N)
229  /::sqrt(max(ores, 0.001)),
230  1.0
231  );
232 
233  co2 = co2Init - co;
234  h2o = h2oInit - h2;
235  ores = oresInit + 0.5*co + 0.5*h2;
236  }
237 
238  thermo reactants
239  (
240  FUEL + o2*O2 + n2*N2
241  );
242 
243  thermo products
244  (
245  fres*FUEL + ores*O2 + n2*N2
246  + co2*CO2 + h2o*H2O + co*CO + h2*H2
247  );
248 
249 
250  scalar equilibriumFlameTemperatureNew =
251  products.THa(reactants.Ha(P, T0), P, adiabaticFlameTemperature);
252 
253  if (j==0)
254  {
255  adiabaticFlameTemperature = equilibriumFlameTemperatureNew;
256  }
257  else
258  {
259  equilibriumFlameTemperature = 0.5*
260  (
261  equilibriumFlameTemperature
262  + equilibriumFlameTemperatureNew
263  );
264  }
265  }
266 
267  Info<< setw(3) << equiv
268  << setw(12) << ft
269  << setw(7) << T0
270  << setw(12) << adiabaticFlameTemperature
271  << setw(12) << equilibriumFlameTemperature
272  << setw(12)
273  << adiabaticFlameTemperature - equilibriumFlameTemperature
274  << setw(12) << ores/N
275  << endl;
276  }
277  }
278 
279  Info<< nl << "end" << endl;
280 
281  return 0;
282 }
283 
284 
285 // ************************************************************************* //
fluidReactionThermo & thermo
Definition: createFields.H:28
A class for handling file names.
Definition: fileName.H:79
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:153
Thermodynamics mapping class to expose the absolute enthalpy functions.
const dimensionSet dimless
A class for handling words, derived from string.
Definition: word.H:59
Functions to search &#39;etc&#39; directories for configuration files etc.
Extract command arguments and options from the supplied argc and argv parameters. ...
Definition: argList.H:102
errorManip< error > abort(error &err)
Definition: errorManip.H:131
fileName findEtcFile(const fileName &, bool mandatory=false)
Search for a file using findEtcFiles.
Definition: etcFiles.C:252
Istream and Ostream manipulators taking arguments.
static const char nl
Definition: Ostream.H:260
Input from file stream.
Definition: IFstream.H:81
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
Definition: thermo.H:52
messageStream Info
label N
Definition: createFields.H:9
label n
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Foam::argList args(argc, argv)
Namespace for OpenFOAM.
scalar T0
Definition: createFields.H:22