FreeFOAM The Cross-Platform CFD Toolkit
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) 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 \*---------------------------------------------------------------------------*/
25 
26 #include "SCOPELaminarFlameSpeed.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 namespace laminarFlameSpeedModels
34 {
35  defineTypeNameAndDebug(SCOPE, 0);
36 
38  (
39  laminarFlameSpeed,
40  SCOPE,
41  dictionary
42  );
43 }
44 }
45 
46 
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 
49 Foam::laminarFlameSpeedModels::SCOPE::polynomial::polynomial
50 (
51  const dictionary& polyDict
52 )
53 :
54  FixedList<scalar, 7>(polyDict.lookup("coefficients")),
55  ll(readScalar(polyDict.lookup("lowerLimit"))),
56  ul(readScalar(polyDict.lookup("upperLimit"))),
57  llv(polyPhi(ll, *this)),
58  ulv(polyPhi(ul, *this)),
59  lu(0)
60 {}
61 
62 
63 Foam::laminarFlameSpeedModels::SCOPE::SCOPE
64 (
65  const dictionary& dict,
66  const hhuCombustionThermo& ct
67 )
68 :
69  laminarFlameSpeed(dict, ct),
70 
71  coeffsDict_(dict.subDict(typeName + "Coeffs").subDict(fuel_)),
72  LFL_(readScalar(coeffsDict_.lookup("lowerFlamabilityLimit"))),
73  UFL_(readScalar(coeffsDict_.lookup("upperFlamabilityLimit"))),
74  SuPolyL_(coeffsDict_.subDict("lowerSuPolynomial")),
75  SuPolyU_(coeffsDict_.subDict("upperSuPolynomial")),
76  Texp_(readScalar(coeffsDict_.lookup("Texp"))),
77  pexp_(readScalar(coeffsDict_.lookup("pexp"))),
78  MaPolyL_(coeffsDict_.subDict("lowerMaPolynomial")),
79  MaPolyU_(coeffsDict_.subDict("upperMaPolynomial"))
80 {
81  SuPolyL_.ll = max(SuPolyL_.ll, LFL_) + SMALL;
82  SuPolyU_.ul = min(SuPolyU_.ul, UFL_) - SMALL;
83 
84  SuPolyL_.lu = 0.5*(SuPolyL_.ul + SuPolyU_.ll);
85  SuPolyU_.lu = SuPolyL_.lu - SMALL;
86 
87  MaPolyL_.lu = 0.5*(MaPolyL_.ul + MaPolyU_.ll);
88  MaPolyU_.lu = MaPolyL_.lu - SMALL;
89 
90  if (debug)
91  {
92  Info<< "phi Su (T = Tref, p = pref)" << endl;
93  label n = 200;
94  for (int i=0; i<n; i++)
95  {
96  scalar phi = (2.0*i)/n;
97  Info<< phi << token::TAB << SuRef(phi) << endl;
98  }
99  }
100 }
101 
102 
103 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
104 
106 {}
107 
108 
109 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
110 
111 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::polyPhi
112 (
113  scalar phi,
114  const polynomial& a
115 )
116 {
117  scalar x = phi - 1.0;
118 
119  return
120  a[0]
121  *(
122  scalar(1)
123  + x*(a[1] + x*(a[2] + x*(a[3] + x*(a[4] + x*(a[5] + x*a[6])))))
124  );
125 }
126 
127 
128 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::SuRef
129 (
130  scalar phi
131 ) const
132 {
133  if (phi < LFL_ || phi > UFL_)
134  {
135  // Return 0 beyond the flamibility limits
136  return scalar(0);
137  }
138  else if (phi < SuPolyL_.ll)
139  {
140  // Use linear interpolation between the low end of the
141  // lower polynomial and the lower flamibility limit
142  return SuPolyL_.llv*(phi - LFL_)/(SuPolyL_.ll - LFL_);
143  }
144  else if (phi > SuPolyU_.ul)
145  {
146  // Use linear interpolation between the upper end of the
147  // upper polynomial and the upper flamibility limit
148  return SuPolyU_.ulv*(UFL_ - phi)/(UFL_ - SuPolyU_.ul);
149  }
150  else if (phi < SuPolyL_.lu)
151  {
152  // Evaluate the lower polynomial
153  return polyPhi(phi, SuPolyL_);
154  }
155  else if (phi > SuPolyU_.lu)
156  {
157  // Evaluate the upper polynomial
158  return polyPhi(phi, SuPolyU_);
159  }
160  else
161  {
162  FatalErrorIn("laminarFlameSpeedModels::SCOPE::SuRef(scalar phi)")
163  << "phi = " << phi
164  << " cannot be handled by SCOPE function with the "
165  "given coefficients"
166  << exit(FatalError);
167 
168  return scalar(0);
169  }
170 }
171 
172 
174 (
175  scalar phi
176 ) const
177 {
178  if (phi < MaPolyL_.ll)
179  {
180  // Beyond the lower limit assume Ma is constant
181  return MaPolyL_.llv;
182  }
183  else if (phi > MaPolyU_.ul)
184  {
185  // Beyond the upper limit assume Ma is constant
186  return MaPolyU_.ulv;
187  }
188  else if (phi < SuPolyL_.lu)
189  {
190  // Evaluate the lower polynomial
191  return polyPhi(phi, MaPolyL_);
192  }
193  else if (phi > SuPolyU_.lu)
194  {
195  // Evaluate the upper polynomial
196  return polyPhi(phi, MaPolyU_);
197  }
198  else
199  {
200  FatalErrorIn("laminarFlameSpeedModels::SCOPE::Ma(scalar phi)")
201  << "phi = " << phi
202  << " cannot be handled by SCOPE function with the "
203  "given coefficients"
204  << exit(FatalError);
205 
206  return scalar(0);
207  }
208 }
209 
210 
211 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
212 (
213  scalar p,
214  scalar Tu,
215  scalar phi
216 ) const
217 {
218  static const scalar Tref = 300.0;
219  static const scalar pRef = 1.013e5;
220 
221  return SuRef(phi)*pow((Tu/Tref), Texp_)*pow((p/pRef), pexp_);
222 }
223 
224 
225 Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
226 (
227  const volScalarField& p,
228  const volScalarField& Tu,
229  scalar phi
230 ) const
231 {
232  tmp<volScalarField> tSu0
233  (
234  new volScalarField
235  (
236  IOobject
237  (
238  "Su0",
239  p.time().timeName(),
240  p.db(),
243  ),
244  p.mesh(),
245  dimensionedScalar("Su0", dimVelocity, 0.0)
246  )
247  );
248 
249  volScalarField& Su0 = tSu0();
250 
251  forAll(Su0, celli)
252  {
253  Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi);
254  }
255 
256  forAll(Su0.boundaryField(), patchi)
257  {
258  scalarField& Su0p = Su0.boundaryField()[patchi];
259  const scalarField& pp = p.boundaryField()[patchi];
260  const scalarField& Tup = Tu.boundaryField()[patchi];
261 
262  forAll(Su0p, facei)
263  {
264  Su0p[facei] = Su0pTphi(pp[facei], Tup[facei], phi);
265  }
266  }
267 
268  return tSu0;
269 }
270 
271 
272 Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
273 (
274  const volScalarField& p,
275  const volScalarField& Tu,
276  const volScalarField& phi
277 ) const
278 {
279  tmp<volScalarField> tSu0
280  (
281  new volScalarField
282  (
283  IOobject
284  (
285  "Su0",
286  p.time().timeName(),
287  p.db(),
290  ),
291  p.mesh(),
292  dimensionedScalar("Su0", dimVelocity, 0.0)
293  )
294  );
295 
296  volScalarField& Su0 = tSu0();
297 
298  forAll(Su0, celli)
299  {
300  Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi[celli]);
301  }
302 
303  forAll(Su0.boundaryField(), patchi)
304  {
305  scalarField& Su0p = Su0.boundaryField()[patchi];
306  const scalarField& pp = p.boundaryField()[patchi];
307  const scalarField& Tup = Tu.boundaryField()[patchi];
308  const scalarField& phip = phi.boundaryField()[patchi];
309 
310  forAll(Su0p, facei)
311  {
312  Su0p[facei] =
313  Su0pTphi
314  (
315  pp[facei],
316  Tup[facei],
317  phip[facei]
318  );
319  }
320  }
321 
322  return tSu0;
323 }
324 
325 
327 (
328  const volScalarField& phi
329 ) const
330 {
331  tmp<volScalarField> tMa
332  (
333  new volScalarField
334  (
335  IOobject
336  (
337  "Ma",
338  phi.time().timeName(),
339  phi.db(),
342  ),
343  phi.mesh(),
344  dimensionedScalar("Ma", dimless, 0.0)
345  )
346  );
347 
348  volScalarField& ma = tMa();
349 
350  forAll(ma, celli)
351  {
352  ma[celli] = Ma(phi[celli]);
353  }
354 
355  forAll(ma.boundaryField(), patchi)
356  {
357  scalarField& map = ma.boundaryField()[patchi];
358  const scalarField& phip = phi.boundaryField()[patchi];
359 
360  forAll(map, facei)
361  {
362  map[facei] = Ma(phip[facei]);
363  }
364  }
365 
366  return tMa;
367 }
368 
369 
372 {
373  if (hhuCombustionThermo_.composition().contains("ft"))
374  {
375  const volScalarField& ft = hhuCombustionThermo_.composition().Y("ft");
376 
377  return Ma
378  (
380  (
381  hhuCombustionThermo_.lookup("stoichiometricAirFuelMassRatio")
382  )*ft/(scalar(1) - ft)
383  );
384  }
385  else
386  {
387  const fvMesh& mesh = hhuCombustionThermo_.p().mesh();
388 
389  return tmp<volScalarField>
390  (
391  new volScalarField
392  (
393  IOobject
394  (
395  "Ma",
396  mesh.time().timeName(),
397  mesh,
400  ),
401  mesh,
402  dimensionedScalar("Ma", dimless, Ma(equivalenceRatio_))
403  )
404  );
405  }
406 }
407 
408 
411 {
412  if (hhuCombustionThermo_.composition().contains("ft"))
413  {
414  const volScalarField& ft = hhuCombustionThermo_.composition().Y("ft");
415 
416  return Su0pTphi
417  (
418  hhuCombustionThermo_.p(),
419  hhuCombustionThermo_.Tu(),
421  (
422  hhuCombustionThermo_.lookup("stoichiometricAirFuelMassRatio")
423  )*ft/(scalar(1) - ft)
424  );
425  }
426  else
427  {
428  return Su0pTphi
429  (
430  hhuCombustionThermo_.p(),
431  hhuCombustionThermo_.Tu(),
432  equivalenceRatio_
433  );
434  }
435 }
436 
437 
438 // ************************ vim: set sw=4 sts=4 et: ************************ //