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-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"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33 namespace laminarFlameSpeedModels
  34 {
  35     defineTypeNameAndDebug(SCOPE, 0);
  36
  37     addToRunTimeSelectionTable
  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
 105 Foam::laminarFlameSpeedModels::SCOPE::~SCOPE()
 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
 173 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::Ma
 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(),
 241                 IOobject::NO_READ,
 242                 IOobject::NO_WRITE
 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(),
 288                 IOobject::NO_READ,
 289                 IOobject::NO_WRITE
 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
 326 Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
 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(),
 340                 IOobject::NO_READ,
 341                 IOobject::NO_WRITE
 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
 370 Foam::tmp<Foam::volScalarField>
 371 Foam::laminarFlameSpeedModels::SCOPE::Ma() const
 372 {
 373     if (hhuCombustionThermo_.composition().contains("ft"))
 374     {
 375         const volScalarField& ft = hhuCombustionThermo_.composition().Y("ft");
 376
 377         return Ma
 378         (
 379             dimensionedScalar
 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,
 398                     IOobject::NO_READ,
 399                     IOobject::NO_WRITE
 400                 ),
 401                 mesh,
 402                 dimensionedScalar("Ma", dimless, Ma(equivalenceRatio_))
 403             )
 404         );
 405     }
 406 }
 407
 408
 409 Foam::tmp<Foam::volScalarField>
 410 Foam::laminarFlameSpeedModels::SCOPE::operator()() const
 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(),
 420             dimensionedScalar
 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 // ************************************************************************* //