singleStepCombustion.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-2020 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 \*---------------------------------------------------------------------------*/
25 
26 #include "singleStepCombustion.H"
27 #include "fvmSup.H"
28 
29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 namespace combustionModels
34 {
35 
36 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
37 
38 template<class ReactionThermo, class ThermoType>
40 {
41  const Reaction<ThermoType>& reaction = reaction_();
42  const scalar Wu = mixture_.specieThermos()[fuelIndex_].W();
43 
44  forAll(reaction.lhs(), i)
45  {
46  const label speciei = reaction.lhs()[i].index;
47  const scalar stoichCoeff = reaction.lhs()[i].stoichCoeff;
48  specieStoichCoeffs_[speciei] = -stoichCoeff;
49  qFuel_.value() += mixture_.specieThermos()[speciei].hc()*stoichCoeff/Wu;
50  }
51 
52  forAll(reaction.rhs(), i)
53  {
54  const label speciei = reaction.rhs()[i].index;
55  const scalar stoichCoeff = reaction.rhs()[i].stoichCoeff;
56  specieStoichCoeffs_[speciei] = stoichCoeff;
57  qFuel_.value() -= mixture_.specieThermos()[speciei].hc()*stoichCoeff/Wu;
58  specieProd_[speciei] = -1;
59  }
60 
61  Info << "Fuel heat of combustion :" << qFuel_.value() << endl;
62 }
63 
64 
65 template<class ReactionThermo, class ThermoType>
67 {
68  const label O2Index = mixture_.species()["O2"];
69  const scalar Wu = mixture_.specieThermos()[fuelIndex_].W();
70 
71  stoicRatio_ =
72  (mixture_.specieThermos()[inertIndex_].W()
73  * specieStoichCoeffs_[inertIndex_]
74  + mixture_.specieThermos()[O2Index].W()
75  * mag(specieStoichCoeffs_[O2Index]))
76  / (Wu*mag(specieStoichCoeffs_[fuelIndex_]));
77 
78  s_ =
79  (mixture_.specieThermos()[O2Index].W()
80  * mag(specieStoichCoeffs_[O2Index]))
81  / (Wu*mag(specieStoichCoeffs_[fuelIndex_]));
82 
83  Info << "stoichiometric air-fuel ratio :" << stoicRatio_.value() << endl;
84 
85  Info << "stoichiometric oxygen-fuel ratio :" << s_.value() << endl;
86 }
87 
88 
89 template<class ReactionThermo, class ThermoType>
91 {
92  const Reaction<ThermoType>& reaction = reaction_();
93 
94  scalar Wm = 0.0;
95  scalar totalMol = 0.0;
96  forAll(reaction.rhs(), i)
97  {
98  label speciei = reaction.rhs()[i].index;
99  totalMol += mag(specieStoichCoeffs_[speciei]);
100  }
101 
102  scalarList Xi(reaction.rhs().size());
103 
104  forAll(reaction.rhs(), i)
105  {
106  const label speciei = reaction.rhs()[i].index;
107  Xi[i] = mag(specieStoichCoeffs_[speciei])/totalMol;
108  Wm += Xi[i]*mixture_.specieThermos()[speciei].W();
109  }
110 
111  forAll(reaction.rhs(), i)
112  {
113  const label speciei = reaction.rhs()[i].index;
114  Yprod0_[speciei] = mixture_.specieThermos()[speciei].W()/Wm*Xi[i];
115  }
116 
117  Info << "Maximum products mass concentrations:" << nl;
118  forAll(Yprod0_, i)
119  {
120  if (Yprod0_[i] > 0)
121  {
122  Info<< " " << mixture_.species()[i] << ": " << Yprod0_[i] << nl;
123  }
124  }
125 
126  // Normalize the stoichiometric coeff to mass
127  forAll(specieStoichCoeffs_, i)
128  {
129  specieStoichCoeffs_[i] =
130  specieStoichCoeffs_[i]
131  * mixture_.specieThermos()[i].W()
132  / (mixture_.specieThermos()[fuelIndex_].W()
133  * mag(specieStoichCoeffs_[fuelIndex_]));
134  }
135 }
136 
137 
138 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
139 
140 template<class ReactionThermo, class ThermoType>
142 (
143  const word& modelType,
144  const ReactionThermo& thermo,
146  const word& combustionProperties
147 )
148 :
149  ThermoCombustion<ReactionThermo>(modelType, thermo, turb),
150  mixture_
151  (
152  dynamic_cast<const multiComponentMixture<ThermoType>&>(this->thermo())
153  ),
154  reaction_
155  (
157  (
158  mixture_.species(),
159  mixture_.specieThermos(),
160  this->subDict("reaction")
161  )
162  ),
163  stoicRatio_(dimensionedScalar("stoicRatio", dimless, 0)),
164  s_(dimensionedScalar("s", dimless, 0)),
165  qFuel_(dimensionedScalar("qFuel", sqr(dimVelocity), 0)),
166  specieStoichCoeffs_(mixture_.species().size(), 0.0),
167  Yprod0_(mixture_.species().size(), 0.0),
168  fres_(Yprod0_.size()),
169  inertIndex_(mixture_.species()[thermo.lookup("inertSpecie")]),
170  fuelIndex_(mixture_.species()[thermo.lookup("fuel")]),
171  specieProd_(Yprod0_.size(), 1),
172  wFuel_
173  (
174  IOobject
175  (
176  this->thermo().phasePropertyName("wFuel"),
177  this->mesh().time().timeName(),
178  this->mesh(),
181  ),
182  this->mesh(),
184  ),
185  semiImplicit_(readBool(this->coeffs_.lookup("semiImplicit")))
186 {
187  forAll(fres_, fresI)
188  {
189  IOobject header
190  (
191  "fres_" + mixture_.species()[fresI],
192  this->mesh().time().timeName(),
193  this->mesh()
194  );
195 
196  fres_.set
197  (
198  fresI,
199  new volScalarField
200  (
201  header,
202  this->mesh(),
203  dimensionedScalar("fres" + name(fresI), dimless, 0)
204  )
205  );
206  }
207 
208  calculateqFuel();
209 
210  massAndAirStoichRatios();
211 
212  calculateMaxProducts();
213 
214  if (semiImplicit_)
215  {
216  Info<< "Combustion mode: semi-implicit" << endl;
217  }
218  else
219  {
220  Info<< "Combustion mode: explicit" << endl;
221  }
222 }
223 
224 
225 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
226 
227 template<class ReactionThermo, class ThermoType>
229 {}
230 
231 
232 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
233 
234 template<class ReactionThermo, class ThermoType>
236 (
238 ) const
239 {
240  const label specieI = mixture_.species()[Y.member()];
241 
242  volScalarField wSpecie
243  (
244  wFuel_*specieStoichCoeffs()[specieI]
245  );
246 
247  if (semiImplicit_)
248  {
249  const label fNorm = specieProd()[specieI];
250  const volScalarField fres(this->fres(specieI));
251  wSpecie /= max(fNorm*(Y - fres), scalar(1e-2));
252 
253  return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
254  }
255  else
256  {
257  return wSpecie + fvm::Sp(0.0*wSpecie, Y);
258  }
259 }
260 
261 
262 template<class ReactionThermo, class ThermoType>
265 {
266  const label fuelI = fuelIndex();
267  volScalarField& YFuel =
268  const_cast<volScalarField&>(this->thermo().composition().Y(fuelI));
269 
270  return -qFuel()*(R(YFuel) & YFuel);
271 }
272 
273 
274 template<class ReactionThermo, class ThermoType>
276 {
278  {
279  return true;
280  }
281  else
282  {
283  return false;
284  }
285 }
286 
287 
288 template<class ReactionThermo, class ThermoType>
290 {
291  const Reaction<ThermoType>& reaction = reaction_();
292 
293  const label O2Index = mixture_.species()["O2"];
294  const volScalarField& YFuel = mixture_.Y()[fuelIndex_];
295  const volScalarField& YO2 = mixture_.Y()[O2Index];
296 
297  // reactants
298  forAll(reaction.lhs(), i)
299  {
300  const label speciei = reaction.lhs()[i].index;
301  if (speciei == fuelIndex_)
302  {
303  fres_[speciei] = max(YFuel - YO2/s_, scalar(0));
304  }
305  else if (speciei == O2Index)
306  {
307  fres_[speciei] = max(YO2 - YFuel*s_, scalar(0));
308  }
309  }
310 
311  // products
312  forAll(reaction.rhs(), i)
313  {
314  const label speciei = reaction.rhs()[i].index;
315  if (speciei != inertIndex_)
316  {
317  forAll(fres_[speciei], celli)
318  {
319  if (fres_[fuelIndex_][celli] > 0.0)
320  {
321  // rich mixture
322  fres_[speciei][celli] =
323  Yprod0_[speciei]
324  * (1.0 + YO2[celli]/s_.value() - YFuel[celli]);
325  }
326  else
327  {
328  // lean mixture
329  fres_[speciei][celli] =
330  Yprod0_[speciei]
331  * (
332  1.0
333  - YO2[celli]/s_.value()*stoicRatio_.value()
334  + YFuel[celli]*stoicRatio_.value()
335  );
336  }
337  }
338  }
339  }
340 }
341 
342 
343 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
344 
345 } // End namespace combustionModels
346 } // End namespace Foam
347 
348 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void calculateMaxProducts()
Calculate maximum products at stoichiometric mixture.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: ReactionI.H:57
bool readBool(Istream &)
Definition: boolIO.C:60
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
rhoReactionThermo & thermo
Definition: createFields.H:28
static word member(const word &name)
Return member (name without the extension)
Definition: IOobject.C:159
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
dynamicFvMesh & mesh
A class for handling words, derived from string.
Definition: word.H:59
Foam::multiComponentMixture.
CombustionModel< rhoReactionThermo > & reaction
word timeName
Definition: getTimeIndex.H:3
static const char nl
Definition: Ostream.H:260
Thermo model wrapper for combustion models.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual tmp< volScalarField > Qdot() const
Heat release rate [kg/m/s^3].
virtual bool read()
Update properties from given dictionary.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
#define R(A, B, C, D, E, F, K, M)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
PtrList< volScalarField > & Y
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
singleStepCombustion(const word &modelType, const ReactionThermo &thermo, const compressibleMomentumTransportModel &turb, const word &combustionProperties)
Construct from components.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
Definition: ReactionI.H:64
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
A class for managing temporary objects.
Definition: PtrList.H:53
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
void fresCorrect()
Calculates the residual for all components.
void massAndAirStoichRatios()
Calculate air/fuel and oxygen/fuel ratio.
Abstract base class for turbulence models (RAS, LES and laminar).
Calculate the matrix for implicit and explicit sources.
virtual tmp< fvScalarMatrix > R(volScalarField &Y) const
Fuel consumption rate matrix.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:812
const dimensionSet dimVelocity