RaviPetersen.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) 2013-2016 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 "RaviPetersen.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 namespace laminarFlameSpeedModels
34 {
35  defineTypeNameAndDebug(RaviPetersen, 0);
36 
38  (
39  laminarFlameSpeed,
40  RaviPetersen,
41  dictionary
42  );
43 }
44 }
45 
46 
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 
49 Foam::laminarFlameSpeedModels::RaviPetersen::RaviPetersen
50 (
51  const dictionary& dict,
52  const psiuReactionThermo& ct
53 )
54 :
55  laminarFlameSpeed(dict, ct),
56  coeffsDict_(dict.subDict(typeName + "Coeffs").subDict(fuel_)),
57  pPoints_(coeffsDict_.lookup("pPoints")),
58  EqRPoints_(coeffsDict_.lookup("EqRPoints")),
59  alpha_(coeffsDict_.lookup("alpha")),
60  beta_(coeffsDict_.lookup("beta")),
61  TRef_(readScalar(coeffsDict_.lookup("TRef")))
62 {
63  checkPointsMonotonicity("equivalenceRatio", EqRPoints_);
64  checkPointsMonotonicity("pressure", pPoints_);
65  checkCoefficientArrayShape("alpha", alpha_);
66  checkCoefficientArrayShape("beta", beta_);
67 }
68 
69 
70 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
71 
73 {}
74 
75 
76 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
77 
78 void Foam::laminarFlameSpeedModels::RaviPetersen::checkPointsMonotonicity
79 (
80  const word& name,
81  const List<scalar>& x
82 ) const
83 {
84  for (label i = 1; i < x.size(); i ++)
85  {
86  if (x[i] <= x[i-1])
87  {
89  (
90  coeffsDict_
91  ) << "Data points for the " << name
92  << " do not increase monotonically" << endl
93  << exit(FatalIOError);
94  }
95  }
96 }
97 
98 
99 void Foam::laminarFlameSpeedModels::RaviPetersen::checkCoefficientArrayShape
100 (
101  const word& name,
102  const List<List<List<scalar>>>& x
103 ) const
104 {
105  bool ok = true;
106 
107  ok &= x.size() == EqRPoints_.size() - 1;
108 
109  forAll(x, i)
110  {
111  ok &= x[i].size() == pPoints_.size();
112 
113  forAll(x[i], j)
114  {
115  ok &= x[i][j].size() == x[i][0].size();
116  }
117  }
118 
119  if (!ok)
120  {
122  (
123  coeffsDict_
124  ) << "Inconsistent size of " << name << " coefficients array" << endl
125  << exit(FatalIOError);
126  }
127 }
128 
129 
130 inline bool Foam::laminarFlameSpeedModels::RaviPetersen::interval
131 (
132  const List<scalar>& xPoints,
133  const scalar x,
134  label& xIndex,
135  scalar& xXi,
136  scalar& xLim
137 ) const
138 {
139  if (x < xPoints.first())
140  {
141  xIndex = 0;
142  xXi = 0.0;
143  xLim = xPoints.first();
144  return false;
145  }
146 
147  else if (x > xPoints.last())
148  {
149  xIndex = xPoints.size() - 2;
150  xXi = 1.0;
151  xLim = xPoints.last();
152  return false;
153  }
154 
155  for (xIndex = 0; x > xPoints[xIndex+1]; xIndex ++)
156  {
157  // increment xIndex until xPoints[xIndex] < x < xPoints[xIndex+1]
158  }
159 
160  xXi = (x - xPoints[xIndex])/(xPoints[xIndex+1] - xPoints[xIndex]);
161  xLim = x;
162 
163  return true;
164 }
165 
166 
167 inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::polynomial
168 (
169  const List<scalar>& coeffs,
170  const scalar x
171 ) const
172 {
173  scalar xPow = 1.0;
174  scalar y = 0.0;
175  forAll(coeffs, i)
176  {
177  y += coeffs[i]*xPow;
178  xPow *= x;
179  }
180  return y;
181 }
182 
183 
184 inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::dPolynomial
185 (
186  const List<scalar>& coeffs,
187  const scalar x
188 ) const
189 {
190  scalar xPow = 1.0;
191  scalar y = 0.0;
192  for (label i = 1; i < coeffs.size(); i++)
193  {
194  y += i*coeffs[i]*xPow;
195  xPow *= x;
196  }
197  return y;
198 }
199 
200 
201 inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::THatPowB
202 (
203  const label EqRIndex,
204  const label pIndex,
205  const scalar EqR,
206  const scalar Tu
207 ) const
208 {
209  return pow
210  (
211  Tu/TRef_,
212  polynomial(beta_[EqRIndex][pIndex],EqR)
213  );
214 }
215 
216 
217 inline Foam::scalar
218 Foam::laminarFlameSpeedModels::RaviPetersen::correlationInRange
219 (
220  const label EqRIndex,
221  const label pIndex,
222  const scalar EqR,
223  const scalar Tu
224 ) const
225 {
226  // standard correlation
227  return
228  polynomial(alpha_[EqRIndex][pIndex],EqR)
229  *THatPowB(EqRIndex, pIndex, EqR, Tu);
230 }
231 
232 
233 inline Foam::scalar
234 Foam::laminarFlameSpeedModels::RaviPetersen::correlationOutOfRange
235 (
236  const label EqRIndex,
237  const label pIndex,
238  const scalar EqR,
239  const scalar EqRLim,
240  const scalar Tu
241 ) const
242 {
243  scalar A = polynomial(alpha_[EqRIndex][pIndex], EqRLim);
244  scalar dA = dPolynomial(alpha_[EqRIndex][pIndex], EqRLim);
245  scalar dB = dPolynomial(beta_[EqRIndex][pIndex], EqRLim);
246  scalar TB = THatPowB(EqRIndex, pIndex, EqRLim, Tu);
247 
248  // linear extrapolation from the bounds of the correlation
249  return max(TB*(A + (dA + A*log(Tu/TRef_)*dB)*(EqR - EqRLim)), 0.0);
250 }
251 
252 
253 inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::speed
254 (
255  const scalar EqR,
256  const scalar p,
257  const scalar Tu
258 ) const
259 {
260  scalar Su = 0, s;
261 
262  label EqRIndex, pIndex;
263  scalar EqRXi, pXi;
264  scalar EqRLim, pLim;
265  bool EqRInRange;
266 
267  EqRInRange = interval(EqRPoints_, EqR, EqRIndex, EqRXi, EqRLim);
268 
269  interval(pPoints_, p, pIndex, pXi, pLim);
270 
271  for (label pI = 0; pI < 2; pI ++)
272  {
273  if (EqRInRange)
274  {
275  s = correlationInRange(EqRIndex, pIndex + pI, EqR, Tu);
276  }
277  else
278  {
279  s = correlationOutOfRange(EqRIndex, pIndex + pI, EqR, EqRLim, Tu);
280  }
281 
282  Su += (1 - pXi)*s;
283  pXi = 1 - pXi;
284  }
285 
286  return Su;
287 }
288 
289 
290 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
291 
294 {
295  const volScalarField& p = psiuReactionThermo_.p();
296  const volScalarField& Tu = psiuReactionThermo_.Tu();
297 
298  volScalarField EqR
299  (
300  IOobject
301  (
302  "EqR",
303  p.time().timeName(),
304  p.db(),
307  false
308  ),
309  p.mesh(),
310  dimensionedScalar("EqR", dimless, 0.0)
311  );
312 
313  if (psiuReactionThermo_.composition().contains("ft"))
314  {
315  const volScalarField& ft = psiuReactionThermo_.composition().Y("ft");
316 
317  EqR =
319  (
320  psiuReactionThermo_.lookup("stoichiometricAirFuelMassRatio")
321  )*ft/max(1 - ft, SMALL);
322  }
323  else
324  {
325  EqR = equivalenceRatio_;
326  }
327 
329  (
330  new volScalarField
331  (
332  IOobject
333  (
334  "Su0",
335  p.time().timeName(),
336  p.db(),
339  false
340  ),
341  p.mesh(),
342  dimensionedScalar("Su0", dimVelocity, 0.0)
343  )
344  );
345 
346  volScalarField& Su0 = tSu0.ref();
347 
348  forAll(Su0, celli)
349  {
350  Su0[celli] = speed(EqR[celli], p[celli], Tu[celli]);
351  }
352 
353  return tSu0;
354 }
355 
356 
357 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
dimensionedScalar log(const dimensionedScalar &ds)
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:221
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:633
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:715
T & first()
Return the first element of the list.
Definition: UListI.H:114
Macros for easy insertion into run-time selection tables.
scalar y
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
A class for handling words, derived from string.
Definition: word.H:59
const Time & time() const
Return time.
Definition: IOobject.C:227
addToRunTimeSelectionTable(laminarFlameSpeed, constant, dictionary)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
Foam::psiuReactionThermo.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
Abstract class for laminar flame speed.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:331
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const Mesh & mesh() const
Return mesh.
A class for managing temporary objects.
Definition: PtrList.H:54
T & last()
Return the last element of the list.
Definition: UListI.H:128
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
tmp< volScalarField > operator()() const
Return the laminar flame speed [m/s].
Definition: RaviPetersen.C:293
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
Namespace for OpenFOAM.
IOerror FatalIOError
const dimensionSet dimVelocity