SCOPELaminarFlameSpeed.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) 2011-2012 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 "IFstream.H"
27 #include "SCOPELaminarFlameSpeed.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 namespace laminarFlameSpeedModels
35 {
36  defineTypeNameAndDebug(SCOPE, 0);
37 
39  (
40  laminarFlameSpeed,
41  SCOPE,
42  dictionary
43  );
44 }
45 }
46 
47 
48 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
49 
50 Foam::laminarFlameSpeedModels::SCOPE::polynomial::polynomial
51 (
52  const dictionary& polyDict
53 )
54 :
55  FixedList<scalar, 7>(polyDict.lookup("coefficients")),
56  ll(readScalar(polyDict.lookup("lowerLimit"))),
57  ul(readScalar(polyDict.lookup("upperLimit"))),
58  llv(polyPhi(ll, *this)),
59  ulv(polyPhi(ul, *this)),
60  lu(0)
61 {}
62 
63 
64 Foam::laminarFlameSpeedModels::SCOPE::SCOPE
65 (
66  const dictionary& dict,
67  const psiuReactionThermo& ct
68 )
69 :
70  laminarFlameSpeed(dict, ct),
71 
72  coeffsDict_
73  (
74  dictionary
75  (
76  IFstream
77  (
78  fileName
79  (
80  dict.lookup("fuelFile")
81  )
82  )()
83  ).subDict(typeName + "Coeffs")
84  ),
85  LFL_(readScalar(coeffsDict_.lookup("lowerFlamabilityLimit"))),
86  UFL_(readScalar(coeffsDict_.lookup("upperFlamabilityLimit"))),
87  SuPolyL_(coeffsDict_.subDict("lowerSuPolynomial")),
88  SuPolyU_(coeffsDict_.subDict("upperSuPolynomial")),
89  Texp_(readScalar(coeffsDict_.lookup("Texp"))),
90  pexp_(readScalar(coeffsDict_.lookup("pexp"))),
91  MaPolyL_(coeffsDict_.subDict("lowerMaPolynomial")),
92  MaPolyU_(coeffsDict_.subDict("upperMaPolynomial"))
93 {
94  SuPolyL_.ll = max(SuPolyL_.ll, LFL_) + SMALL;
95  SuPolyU_.ul = min(SuPolyU_.ul, UFL_) - SMALL;
96 
97  SuPolyL_.lu = 0.5*(SuPolyL_.ul + SuPolyU_.ll);
98  SuPolyU_.lu = SuPolyL_.lu - SMALL;
99 
100  MaPolyL_.lu = 0.5*(MaPolyL_.ul + MaPolyU_.ll);
101  MaPolyU_.lu = MaPolyL_.lu - SMALL;
102 
103  if (debug)
104  {
105  Info<< "phi Su (T = Tref, p = pref)" << endl;
106  label n = 200;
107  for (int i=0; i<n; i++)
108  {
109  scalar phi = (2.0*i)/n;
110  Info<< phi << token::TAB << SuRef(phi) << endl;
111  }
112  }
113 }
114 
115 
116 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
117 
119 {}
120 
121 
122 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
123 
124 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::polyPhi
125 (
126  scalar phi,
127  const polynomial& a
128 )
129 {
130  scalar x = phi - 1.0;
131 
132  return
133  a[0]
134  *(
135  scalar(1)
136  + x*(a[1] + x*(a[2] + x*(a[3] + x*(a[4] + x*(a[5] + x*a[6])))))
137  );
138 }
139 
140 
141 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::SuRef
142 (
143  scalar phi
144 ) const
145 {
146  if (phi < LFL_ || phi > UFL_)
147  {
148  // Return 0 beyond the flamibility limits
149  return scalar(0);
150  }
151  else if (phi < SuPolyL_.ll)
152  {
153  // Use linear interpolation between the low end of the
154  // lower polynomial and the lower flamibility limit
155  return SuPolyL_.llv*(phi - LFL_)/(SuPolyL_.ll - LFL_);
156  }
157  else if (phi > SuPolyU_.ul)
158  {
159  // Use linear interpolation between the upper end of the
160  // upper polynomial and the upper flamibility limit
161  return SuPolyU_.ulv*(UFL_ - phi)/(UFL_ - SuPolyU_.ul);
162  }
163  else if (phi < SuPolyL_.lu)
164  {
165  // Evaluate the lower polynomial
166  return polyPhi(phi, SuPolyL_);
167  }
168  else if (phi > SuPolyU_.lu)
169  {
170  // Evaluate the upper polynomial
171  return polyPhi(phi, SuPolyU_);
172  }
173  else
174  {
175  FatalErrorIn("laminarFlameSpeedModels::SCOPE::SuRef(scalar phi)")
176  << "phi = " << phi
177  << " cannot be handled by SCOPE function with the "
178  "given coefficients"
179  << exit(FatalError);
180 
181  return scalar(0);
182  }
183 }
184 
185 
187 (
188  scalar phi
189 ) const
190 {
191  if (phi < MaPolyL_.ll)
192  {
193  // Beyond the lower limit assume Ma is constant
194  return MaPolyL_.llv;
195  }
196  else if (phi > MaPolyU_.ul)
197  {
198  // Beyond the upper limit assume Ma is constant
199  return MaPolyU_.ulv;
200  }
201  else if (phi < SuPolyL_.lu)
202  {
203  // Evaluate the lower polynomial
204  return polyPhi(phi, MaPolyL_);
205  }
206  else if (phi > SuPolyU_.lu)
207  {
208  // Evaluate the upper polynomial
209  return polyPhi(phi, MaPolyU_);
210  }
211  else
212  {
213  FatalErrorIn("laminarFlameSpeedModels::SCOPE::Ma(scalar phi)")
214  << "phi = " << phi
215  << " cannot be handled by SCOPE function with the "
216  "given coefficients"
217  << exit(FatalError);
218 
219  return scalar(0);
220  }
221 }
222 
223 
224 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
225 (
226  scalar p,
227  scalar Tu,
228  scalar phi
229 ) const
230 {
231  static const scalar Tref = 300.0;
232  static const scalar pRef = 1.013e5;
233 
234  return SuRef(phi)*pow((Tu/Tref), Texp_)*pow((p/pRef), pexp_);
235 }
236 
237 
238 Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
239 (
240  const volScalarField& p,
241  const volScalarField& Tu,
242  scalar phi
243 ) const
244 {
245  tmp<volScalarField> tSu0
246  (
247  new volScalarField
248  (
249  IOobject
250  (
251  "Su0",
252  p.time().timeName(),
253  p.db(),
256  ),
257  p.mesh(),
258  dimensionedScalar("Su0", dimVelocity, 0.0)
259  )
260  );
261 
262  volScalarField& Su0 = tSu0();
263 
264  forAll(Su0, celli)
265  {
266  Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi);
267  }
268 
269  forAll(Su0.boundaryField(), patchi)
270  {
271  scalarField& Su0p = Su0.boundaryField()[patchi];
272  const scalarField& pp = p.boundaryField()[patchi];
273  const scalarField& Tup = Tu.boundaryField()[patchi];
274 
275  forAll(Su0p, facei)
276  {
277  Su0p[facei] = Su0pTphi(pp[facei], Tup[facei], phi);
278  }
279  }
280 
281  return tSu0;
282 }
283 
284 
285 Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
286 (
287  const volScalarField& p,
288  const volScalarField& Tu,
289  const volScalarField& phi
290 ) const
291 {
292  tmp<volScalarField> tSu0
293  (
294  new volScalarField
295  (
296  IOobject
297  (
298  "Su0",
299  p.time().timeName(),
300  p.db(),
303  ),
304  p.mesh(),
305  dimensionedScalar("Su0", dimVelocity, 0.0)
306  )
307  );
308 
309  volScalarField& Su0 = tSu0();
310 
311  forAll(Su0, celli)
312  {
313  Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi[celli]);
314  }
315 
316  forAll(Su0.boundaryField(), patchi)
317  {
318  scalarField& Su0p = Su0.boundaryField()[patchi];
319  const scalarField& pp = p.boundaryField()[patchi];
320  const scalarField& Tup = Tu.boundaryField()[patchi];
321  const scalarField& phip = phi.boundaryField()[patchi];
322 
323  forAll(Su0p, facei)
324  {
325  Su0p[facei] =
326  Su0pTphi
327  (
328  pp[facei],
329  Tup[facei],
330  phip[facei]
331  );
332  }
333  }
334 
335  return tSu0;
336 }
337 
338 
340 (
341  const volScalarField& phi
342 ) const
343 {
344  tmp<volScalarField> tMa
345  (
346  new volScalarField
347  (
348  IOobject
349  (
350  "Ma",
351  phi.time().timeName(),
352  phi.db(),
355  ),
356  phi.mesh(),
357  dimensionedScalar("Ma", dimless, 0.0)
358  )
359  );
360 
361  volScalarField& ma = tMa();
362 
363  forAll(ma, celli)
364  {
365  ma[celli] = Ma(phi[celli]);
366  }
367 
368  forAll(ma.boundaryField(), patchi)
369  {
370  scalarField& map = ma.boundaryField()[patchi];
371  const scalarField& phip = phi.boundaryField()[patchi];
372 
373  forAll(map, facei)
374  {
375  map[facei] = Ma(phip[facei]);
376  }
377  }
378 
379  return tMa;
380 }
381 
382 
385 {
387  {
388  const volScalarField& ft = psiuReactionThermo_.composition().Y("ft");
389 
390  return Ma
391  (
393  (
394  psiuReactionThermo_.lookup("stoichiometricAirFuelMassRatio")
395  )*ft/(scalar(1) - ft)
396  );
397  }
398  else
399  {
400  const fvMesh& mesh = psiuReactionThermo_.p().mesh();
401 
402  return tmp<volScalarField>
403  (
404  new volScalarField
405  (
406  IOobject
407  (
408  "Ma",
409  mesh.time().timeName(),
410  mesh,
413  ),
414  mesh,
416  )
417  );
418  }
419 }
420 
421 
424 {
426  {
427  const volScalarField& ft = psiuReactionThermo_.composition().Y("ft");
428 
429  return Su0pTphi
430  (
434  (
435  psiuReactionThermo_.lookup("stoichiometricAirFuelMassRatio")
436  )*ft/(scalar(1) - ft)
437  );
438  }
439  else
440  {
441  return Su0pTphi
442  (
446  );
447  }
448 }
449 
450 
451 // ************************************************************************* //
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
#define readScalar
Definition: doubleScalar.C:38
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
virtual volScalarField & p()
Pressure [Pa].
Definition: basicThermo.C:466
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
tmp< volScalarField > operator()() const
Return the laminar flame speed [m/s].
messageStream Info
const psiuReactionThermo & psiuReactionThermo_
const Mesh & mesh() const
Return mesh.
dynamicFvMesh & mesh
virtual basicSpecieMixture & composition()=0
Return the composition of the multi-component mixture.
Namespace for OpenFOAM.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
label n
scalar equivalenceRatio_
Equivalence ratio of a homogeneous mixture.
dictionary dict
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
tmp< volScalarField > Ma() const
Return the Markstein number.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
label patchi
Macros for easy insertion into run-time selection tables.
surfaceScalarField & phi
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
PtrList< volScalarField > & Y()
Return the mass-fraction fields.
addToRunTimeSelectionTable(laminarFlameSpeed, constant, dictionary)
virtual const volScalarField & Tu() const =0
Unburnt gas temperature [K].
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
error FatalError
bool contains(const word &specieName) const
Does the mixture include this specie?
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const dimensionSet dimVelocity
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118