applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C

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