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) 2004-2011 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 "IFstream.H"
  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_
  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
 118 Foam::laminarFlameSpeedModels::SCOPE::~SCOPE()
 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
 186 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::Ma
 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(),
 254                 IOobject::NO_READ,
 255                 IOobject::NO_WRITE
 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(),
 301                 IOobject::NO_READ,
 302                 IOobject::NO_WRITE
 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
 339 Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
 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(),
 353                 IOobject::NO_READ,
 354                 IOobject::NO_WRITE
 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
 383 Foam::tmp<Foam::volScalarField>
 384 Foam::laminarFlameSpeedModels::SCOPE::Ma() const
 385 {
 386     if (hhuCombustionThermo_.composition().contains("ft"))
 387     {
 388         const volScalarField& ft = hhuCombustionThermo_.composition().Y("ft");
 389
 390         return Ma
 391         (
 392             dimensionedScalar
 393             (
 394                 hhuCombustionThermo_.lookup("stoichiometricAirFuelMassRatio")
 395             )*ft/(scalar(1) - ft)
 396         );
 397     }
 398     else
 399     {
 400         const fvMesh& mesh = hhuCombustionThermo_.p().mesh();
 401
 402         return tmp<volScalarField>
 403         (
 404             new volScalarField
 405             (
 406                 IOobject
 407                 (
 408                     "Ma",
 409                     mesh.time().timeName(),
 410                     mesh,
 411                     IOobject::NO_READ,
 412                     IOobject::NO_WRITE
 413                 ),
 414                 mesh,
 415                 dimensionedScalar("Ma", dimless, Ma(equivalenceRatio_))
 416             )
 417         );
 418     }
 419 }
 420
 421
 422 Foam::tmp<Foam::volScalarField>
 423 Foam::laminarFlameSpeedModels::SCOPE::operator()() const
 424 {
 425     if (hhuCombustionThermo_.composition().contains("ft"))
 426     {
 427         const volScalarField& ft = hhuCombustionThermo_.composition().Y("ft");
 428
 429         return Su0pTphi
 430         (
 431             hhuCombustionThermo_.p(),
 432             hhuCombustionThermo_.Tu(),
 433             dimensionedScalar
 434             (
 435                 hhuCombustionThermo_.lookup("stoichiometricAirFuelMassRatio")
 436             )*ft/(scalar(1) - ft)
 437         );
 438     }
 439     else
 440     {
 441         return Su0pTphi
 442         (
 443             hhuCombustionThermo_.p(),
 444             hhuCombustionThermo_.Tu(),
 445             equivalenceRatio_
 446         );
 447     }
 448 }
 449
 450
 451 // ************************************************************************* //