RaviPetersen.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) 2013-2023 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 {
36 
38  (
42  );
43 }
44 }
45 
46 
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 
50 (
51  const dictionary& dict,
52  const psiuMulticomponentThermo& ct
53 )
54 :
56  coeffsDict_(dict.optionalSubDict(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_(coeffsDict_.lookup<scalar>("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 = psiuMulticomponentThermo_.p();
296  const volScalarField& Tu = psiuMulticomponentThermo_.Tu();
297 
298  volScalarField EqR
299  (
300  IOobject
301  (
302  "EqR",
303  p.time().name(),
304  p.db(),
307  false
308  ),
309  p.mesh(),
311  );
312 
313  if (psiuMulticomponentThermo_.containsSpecie("ft"))
314  {
315  const volScalarField& ft = psiuMulticomponentThermo_.Y("ft");
316 
317  EqR =
319  (
320  "stoichiometricAirFuelMassRatio",
321  dimless,
322  psiuMulticomponentThermo_.properties()
323  )*ft/max(1 - ft, small);
324  }
325  else
326  {
327  EqR = equivalenceRatio_;
328  }
329 
331  (
333  (
334  "Su0",
335  p.mesh(),
337  )
338  );
339 
340  volScalarField& Su0 = tSu0.ref();
341 
342  forAll(Su0, celli)
343  {
344  Su0[celli] = speed(EqR[celli], p[celli], Tu[celli]);
345  }
346 
347  return tSu0;
348 }
349 
350 
351 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
scalar y
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
Generic GeometricField class.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &, const HashPtrTable< Source > &=HashPtrTable< Source >())
Return a temporary field constructed from name,.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
Laminar flame speed obtained from Ravi and Petersen's correlation.
Definition: RaviPetersen.H:64
tmp< volScalarField > operator()() const
Return the laminar flame speed [m/s].
Definition: RaviPetersen.C:293
RaviPetersen(const dictionary &, const psiuMulticomponentThermo &)
Construct from dictionary and psiuMulticomponentThermo.
Definition: RaviPetersen.C:50
Abstract class for laminar flame speed.
Base-class for combustion fluid thermodynamic properties based on compressibility.
A class for managing temporary objects.
Definition: tmp.H:55
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
A class for handling words, derived from string.
Definition: word.H:62
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:346
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.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
tmp< VolField< Type > > Su(const VolField< Type > &su, const VolField< Type > &vf)
Definition: fvcSup.C:44
addToRunTimeSelectionTable(laminarFlameSpeed, constant, dictionary)
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
const dimensionSet dimless
dimensionedScalar log(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
IOerror FatalIOError
const dimensionSet dimVelocity
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dictionary dict
volScalarField & p