FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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.
29 
30  Calculates the equilibrium flame temperature for a given fuel and
31  pressure for a range of unburnt gas temperatures and equivalence
32  ratios; the effects of dissociation on O2, H2O and CO2 are included.
33 
34 Usage
35 
36  - equilibriumFlameT [OPTIONS] <controlFile>
37 
38  @param <controlFile> \n
39  @todo Detailed description of argument.
40 
41  @param -case <dir>\n
42  Case directory.
43 
44  @param -parallel \n
45  Run in parallel.
46 
47  @param -help \n
48  Display help message.
49 
50  @param -doc \n
51  Display Doxygen API documentation page for this application.
52 
53  @param -srcDoc \n
54  Display Doxygen source documentation page for this application.
55 
56 \*---------------------------------------------------------------------------*/
57 
58 #include <OpenFOAM/argList.H>
59 #include <OpenFOAM/Time.H>
60 #include <OpenFOAM/dictionary.H>
61 #include <OpenFOAM/IFstream.H>
62 #include <OpenFOAM/OSspecific.H>
63 #include <OpenFOAM/IOmanip.H>
64 
65 #include <specie/specieThermo.H>
66 #include <specie/janafThermo.H>
67 #include <specie/perfectGas.H>
68 
69 using namespace Foam;
70 
72 
73 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
74 
75 int main(int argc, char *argv[])
76 {
78  argList::validArgs.append("controlFile");
79  argList args(argc, argv);
80 
81  fileName controlFileName(args.additionalArgs()[0]);
82 
83  // Construct control dictionary
84  IFstream controlFile(controlFileName);
85 
86  // Check controlFile stream is OK
87  if (!controlFile.good())
88  {
90  << "Cannot read file " << controlFileName
91  << abort(FatalError);
92  }
93 
94  dictionary control(controlFile);
95 
96 
97  scalar P(readScalar(control.lookup("P")));
98  word fuel(control.lookup("fuel"));
99  scalar n(readScalar(control.lookup("n")));
100  scalar m(readScalar(control.lookup("m")));
101 
102 
103  Info<< nl << "Reading Burcat data dictionary" << endl;
104 
105  fileName BurcatCpDataFileName(findEtcFile("thermoData/BurcatCpData"));
106 
107  // Construct control dictionary
108  IFstream BurcatCpDataFile(BurcatCpDataFileName);
109 
110  // Check BurcatCpData stream is OK
111  if (!BurcatCpDataFile.good())
112  {
114  << "Cannot read file " << BurcatCpDataFileName
115  << abort(FatalError);
116  }
117 
118  dictionary thermoData(BurcatCpDataFile);
119 
120 
121  Info<< nl << "Reading Burcat data for relevant species" << nl << endl;
122 
123  // Reactants
124  thermo FUEL(thermoData.lookup(fuel));
125  thermo O2(thermoData.lookup("O2"));
126  thermo N2(thermoData.lookup("N2"));
127 
128  // Products
129  thermo CO2(thermoData.lookup("CO2"));
130  thermo H2O(thermoData.lookup("H2O"));
131 
132  // Product fragments
133  thermo CO(thermoData.lookup("CO"));
134  thermo H2(thermoData.lookup("H2"));
135 
136 
137  // Product dissociation reactions
138 
139  thermo CO2BreakUp
140  (
141  CO2 == CO + 0.5* O2
142  );
143 
144  thermo H2OBreakUp
145  (
146  H2O == H2 + 0.5*O2
147  );
148 
149 
150  // Stoiciometric number of moles of species for one mole of fuel
151  scalar stoicO2 = n + m/4.0;
152  scalar stoicN2 = (0.79/0.21)*(n + m/4.0);
153  scalar stoicCO2 = n;
154  scalar stoicH2O = m/2.0;
155 
156  // Oxidant
157  thermo oxidant
158  (
159  "oxidant",
160  stoicO2*O2
161  + stoicN2*N2
162  );
163 
164  dimensionedScalar stoichiometricAirFuelMassRatio
165  (
166  "stoichiometricAirFuelMassRatio",
167  dimless,
168  (oxidant.W()*oxidant.nMoles())/FUEL.W()
169  );
170 
171  Info<< "stoichiometricAirFuelMassRatio "
172  << stoichiometricAirFuelMassRatio << ';' << endl;
173 
174  Info<< "Equilibrium flame temperature data ("
175  << P/1e5 << " bar)" << nl << nl
176  << setw(3) << "Phi"
177  << setw(12) << "ft"
178  << setw(7) << "T0"
179  << setw(12) << "Tad"
180  << setw(12) << "Teq"
181  << setw(12) << "Terror"
182  << setw(20) << "O2res (mole frac)" << nl
183  << endl;
184 
185 
186  // Loop over equivalence ratios
187  for (int i=0; i<16; i++)
188  {
189  scalar equiv = 0.6 + i*0.05;
190  scalar ft = 1/(1 + stoichiometricAirFuelMassRatio.value()/equiv);
191 
192  // Loop over initial temperatures
193  for (int j=0; j<28; j++)
194  {
195  scalar T0 = 300.0 + j*100.0;
196 
197  // Number of moles of species for one mole of fuel
198  scalar o2 = (1.0/equiv)*stoicO2;
199  scalar n2 = (0.79/0.21)*o2;
200  scalar fres = max(1.0 - 1.0/equiv, 0.0);
201  scalar fburnt = 1.0 - fres;
202 
203  // Initial guess for number of moles of product species
204  // ignoring product dissociation
205  scalar oresInit = max(1.0/equiv - 1.0, 0.0)*stoicO2;
206  scalar co2Init = fburnt*stoicCO2;
207  scalar h2oInit = fburnt*stoicH2O;
208 
209  scalar ores = oresInit;
210  scalar co2 = co2Init;
211  scalar h2o = h2oInit;
212 
213  scalar co = 0.0;
214  scalar h2 = 0.0;
215 
216  // Total number of moles in system
217  scalar N = fres + n2 + co2 + h2o + ores;
218 
219 
220  // Initial guess for adiabatic flame temperature
221  scalar adiabaticFlameTemperature =
222  T0
223  + (fburnt/(1.0 + o2 + n2))/(1.0/(1.0 + (1.0 + 0.79/0.21)*stoicO2))
224  *2000.0;
225 
226  scalar equilibriumFlameTemperature = adiabaticFlameTemperature;
227 
228 
229  // Iteration loop for adiabatic flame temperature
230  for (int j=0; j<20; j++)
231  {
232 
233  if (j > 0)
234  {
235  co = co2*
236  min
237  (
238  CO2BreakUp.Kn(equilibriumFlameTemperature, P, N)
239  /::sqrt(max(ores, 0.001)),
240  1.0
241  );
242 
243  h2 = h2o*
244  min
245  (
246  H2OBreakUp.Kn(equilibriumFlameTemperature, P, N)
247  /::sqrt(max(ores, 0.001)),
248  1.0
249  );
250 
251  co2 = co2Init - co;
252  h2o = h2oInit - h2;
253  ores = oresInit + 0.5*co + 0.5*h2;
254  }
255 
256  thermo reactants
257  (
258  FUEL + o2*O2 + n2*N2
259  );
260 
261  thermo products
262  (
263  fres*FUEL + ores*O2 + n2*N2
264  + co2*CO2 + h2o*H2O + co*CO + h2*H2
265  );
266 
267 
268  scalar equilibriumFlameTemperatureNew =
269  products.TH(reactants.H(T0), adiabaticFlameTemperature);
270 
271  if (j==0)
272  {
273  adiabaticFlameTemperature = equilibriumFlameTemperatureNew;
274  }
275  else
276  {
277  equilibriumFlameTemperature = 0.5*
278  (
279  equilibriumFlameTemperature
280  + equilibriumFlameTemperatureNew
281  );
282  }
283  }
284 
285  Info<< setw(3) << equiv
286  << setw(12) << ft
287  << setw(7) << T0
288  << setw(12) << adiabaticFlameTemperature
289  << setw(12) << equilibriumFlameTemperature
290  << setw(12)
291  << adiabaticFlameTemperature - equilibriumFlameTemperature
292  << setw(12) << ores/N
293  << endl;
294  }
295  }
296 
297  Info<< nl << "end" << endl;
298 
299  return 0;
300 }
301 
302 
303 // ************************ vim: set sw=4 sts=4 et: ************************ //