src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-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 "ODEChemistryModel.H"
  28 #include "chemistrySolver.H"
  29 #include "reactingMixture.H"
  30
  31 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  32
  33 template<class CompType, class ThermoType>
  34 Foam::ODEChemistryModel<CompType, ThermoType>::ODEChemistryModel
  35 (
  36     const fvMesh& mesh,
  37     const word& compTypeName,
  38     const word& thermoTypeName
  39 )
  40 :
  41     CompType(mesh, thermoTypeName),
  42
  43     ODE(),
  44
  45     Y_(this->thermo().composition().Y()),
  46
  47     reactions_
  48     (
  49         dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo())
  50     ),
  51     specieThermo_
  52     (
  53         dynamic_cast<const reactingMixture<ThermoType>&>
  54             (this->thermo()).speciesData()
  55     ),
  56
  57     nSpecie_(Y_.size()),
  58     nReaction_(reactions_.size()),
  59
  60     solver_
  61     (
  62         chemistrySolver<CompType, ThermoType>::New
  63         (
  64             *this,
  65             compTypeName,
  66             thermoTypeName
  67         )
  68     ),
  69
  70     RR_(nSpecie_)
  71 {
  72     // create the fields for the chemistry sources
  73     forAll(RR_, fieldI)
  74     {
  75         RR_.set
  76         (
  77             fieldI,
  78             new scalarField(mesh.nCells(), 0.0)
  79         );
  80     }
  81
  82     Info<< "ODEChemistryModel: Number of species = " << nSpecie_
  83         << " and reactions = " << nReaction_ << endl;
  84 }
  85
  86
  87 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  88
  89 template<class CompType, class ThermoType>
  90 Foam::ODEChemistryModel<CompType, ThermoType>::~ODEChemistryModel()
  91 {}
  92
  93
  94 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  95
  96 template<class CompType, class ThermoType>
  97 Foam::scalarField Foam::ODEChemistryModel<CompType, ThermoType>::omega
  98 (
  99     const scalarField& c,
 100     const scalar T,
 101     const scalar p
 102 ) const
 103 {
 104     scalar pf, cf, pr, cr;
 105     label lRef, rRef;
 106
 107     scalarField om(nEqns(), 0.0);
 108
 109     forAll(reactions_, i)
 110     {
 111         const Reaction<ThermoType>& R = reactions_[i];
 112
 113         scalar omegai = omega
 114         (
 115             R, c, T, p, pf, cf, lRef, pr, cr, rRef
 116         );
 117
 118         forAll(R.lhs(), s)
 119         {
 120             label si = R.lhs()[s].index;
 121             scalar sl = R.lhs()[s].stoichCoeff;
 122             om[si] -= sl*omegai;
 123         }
 124
 125         forAll(R.rhs(), s)
 126         {
 127             label si = R.rhs()[s].index;
 128             scalar sr = R.rhs()[s].stoichCoeff;
 129             om[si] += sr*omegai;
 130         }
 131     }
 132
 133     return om;
 134 }
 135
 136
 137 template<class CompType, class ThermoType>
 138 Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
 139 (
 140     const Reaction<ThermoType>& R,
 141     const scalarField& c,
 142     const scalar T,
 143     const scalar p,
 144     scalar& pf,
 145     scalar& cf,
 146     label& lRef,
 147     scalar& pr,
 148     scalar& cr,
 149     label& rRef
 150 ) const
 151 {
 152     scalarField c2(nSpecie_, 0.0);
 153     for (label i=0; i<nSpecie_; i++)
 154     {
 155         c2[i] = max(0.0, c[i]);
 156     }
 157
 158     scalar kf = R.kf(T, p, c2);
 159     scalar kr = R.kr(kf, T, p, c2);
 160
 161     pf = 1.0;
 162     pr = 1.0;
 163
 164     label Nl = R.lhs().size();
 165     label Nr = R.rhs().size();
 166
 167     label slRef = 0;
 168     lRef = R.lhs()[slRef].index;
 169
 170     pf = kf;
 171     for (label s=1; s<Nl; s++)
 172     {
 173         label si = R.lhs()[s].index;
 174
 175         if (c[si] < c[lRef])
 176         {
 177             scalar exp = R.lhs()[slRef].exponent;
 178             pf *= pow(max(0.0, c[lRef]), exp);
 179             lRef = si;
 180             slRef = s;
 181         }
 182         else
 183         {
 184             scalar exp = R.lhs()[s].exponent;
 185             pf *= pow(max(0.0, c[si]), exp);
 186         }
 187     }
 188     cf = max(0.0, c[lRef]);
 189
 190     {
 191         scalar exp = R.lhs()[slRef].exponent;
 192         if (exp<1.0)
 193         {
 194             if (cf > SMALL)
 195             {
 196                 pf *= pow(cf, exp - 1.0);
 197             }
 198             else
 199             {
 200                 pf = 0.0;
 201             }
 202         }
 203         else
 204         {
 205             pf *= pow(cf, exp - 1.0);
 206         }
 207     }
 208
 209     label srRef = 0;
 210     rRef = R.rhs()[srRef].index;
 211
 212     // find the matrix element and element position for the rhs
 213     pr = kr;
 214     for (label s=1; s<Nr; s++)
 215     {
 216         label si = R.rhs()[s].index;
 217         if (c[si] < c[rRef])
 218         {
 219             scalar exp = R.rhs()[srRef].exponent;
 220             pr *= pow(max(0.0, c[rRef]), exp);
 221             rRef = si;
 222             srRef = s;
 223         }
 224         else
 225         {
 226             scalar exp = R.rhs()[s].exponent;
 227             pr *= pow(max(0.0, c[si]), exp);
 228         }
 229     }
 230     cr = max(0.0, c[rRef]);
 231
 232     {
 233         scalar exp = R.rhs()[srRef].exponent;
 234         if (exp<1.0)
 235         {
 236             if (cr>SMALL)
 237             {
 238                 pr *= pow(cr, exp - 1.0);
 239             }
 240             else
 241             {
 242                 pr = 0.0;
 243             }
 244         }
 245         else
 246         {
 247             pr *= pow(cr, exp - 1.0);
 248         }
 249     }
 250
 251     return pf*cf - pr*cr;
 252 }
 253
 254
 255 template<class CompType, class ThermoType>
 256 void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
 257 (
 258     const scalar time,
 259     const scalarField &c,
 260     scalarField& dcdt
 261 ) const
 262 {
 263     scalar T = c[nSpecie_];
 264     scalar p = c[nSpecie_ + 1];
 265
 266     dcdt = omega(c, T, p);
 267
 268     // constant pressure
 269     // dT/dt = ...
 270     scalar rho = 0.0;
 271     scalar cSum = 0.0;
 272     for (label i=0; i<nSpecie_; i++)
 273     {
 274         scalar W = specieThermo_[i].W();
 275         cSum += c[i];
 276         rho += W*c[i];
 277     }
 278     scalar mw = rho/cSum;
 279     scalar cp = 0.0;
 280     for (label i=0; i<nSpecie_; i++)
 281     {
 282         scalar cpi = specieThermo_[i].cp(T);
 283         scalar Xi = c[i]/rho;
 284         cp += Xi*cpi;
 285     }
 286     cp /= mw;
 287
 288     scalar dT = 0.0;
 289     for (label i=0; i<nSpecie_; i++)
 290     {
 291         scalar hi = specieThermo_[i].h(T);
 292         dT += hi*dcdt[i];
 293     }
 294     dT /= rho*cp;
 295
 296     // limit the time-derivative, this is more stable for the ODE
 297     // solver when calculating the allowed time step
 298     scalar dtMag = min(500.0, mag(dT));
 299     dcdt[nSpecie_] = -dT*dtMag/(mag(dT) + 1.0e-10);
 300
 301     // dp/dt = ...
 302     dcdt[nSpecie_+1] = 0.0;
 303 }
 304
 305
 306 template<class CompType, class ThermoType>
 307 void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
 308 (
 309     const scalar t,
 310     const scalarField& c,
 311     scalarField& dcdt,
 312     scalarSquareMatrix& dfdc
 313 ) const
 314 {
 315     scalar T = c[nSpecie_];
 316     scalar p = c[nSpecie_ + 1];
 317
 318     scalarField c2(nSpecie_, 0.0);
 319     for (label i=0; i<nSpecie_; i++)
 320     {
 321         c2[i] = max(c[i], 0.0);
 322     }
 323
 324     for (label i=0; i<nEqns(); i++)
 325     {
 326         for (label j=0; j<nEqns(); j++)
 327         {
 328             dfdc[i][j] = 0.0;
 329         }
 330     }
 331
 332     // length of the first argument must be nSpecie()
 333     dcdt = omega(c2, T, p);
 334
 335     for (label ri=0; ri<reactions_.size(); ri++)
 336     {
 337         const Reaction<ThermoType>& R = reactions_[ri];
 338
 339         scalar kf0 = R.kf(T, p, c2);
 340         scalar kr0 = R.kr(T, p, c2);
 341
 342         forAll(R.lhs(), j)
 343         {
 344             label sj = R.lhs()[j].index;
 345             scalar kf = kf0;
 346             forAll(R.lhs(), i)
 347             {
 348                 label si = R.lhs()[i].index;
 349                 scalar el = R.lhs()[i].exponent;
 350                 if (i == j)
 351                 {
 352                     if (el < 1.0)
 353                     {
 354                         if (c2[si]>SMALL)
 355                         {
 356                             kf *= el*pow(c2[si] + VSMALL, el - 1.0);
 357                         }
 358                         else
 359                         {
 360                             kf = 0.0;
 361                         }
 362                     }
 363                     else
 364                     {
 365                         kf *= el*pow(c2[si], el - 1.0);
 366                     }
 367                 }
 368                 else
 369                 {
 370                     kf *= pow(c2[si], el);
 371                 }
 372             }
 373
 374             forAll(R.lhs(), i)
 375             {
 376                 label si = R.lhs()[i].index;
 377                 scalar sl = R.lhs()[i].stoichCoeff;
 378                 dfdc[si][sj] -= sl*kf;
 379             }
 380             forAll(R.rhs(), i)
 381             {
 382                 label si = R.rhs()[i].index;
 383                 scalar sr = R.rhs()[i].stoichCoeff;
 384                 dfdc[si][sj] += sr*kf;
 385             }
 386         }
 387
 388         forAll(R.rhs(), j)
 389         {
 390             label sj = R.rhs()[j].index;
 391             scalar kr = kr0;
 392             forAll(R.rhs(), i)
 393             {
 394                 label si = R.rhs()[i].index;
 395                 scalar er = R.rhs()[i].exponent;
 396                 if (i==j)
 397                 {
 398                     if (er<1.0)
 399                     {
 400                         if (c2[si]>SMALL)
 401                         {
 402                             kr *= er*pow(c2[si] + VSMALL, er - 1.0);
 403                         }
 404                         else
 405                         {
 406                             kr = 0.0;
 407                         }
 408                     }
 409                     else
 410                     {
 411                         kr *= er*pow(c2[si], er - 1.0);
 412                     }
 413                 }
 414                 else
 415                 {
 416                     kr *= pow(c2[si], er);
 417                 }
 418             }
 419
 420             forAll(R.lhs(), i)
 421             {
 422                 label si = R.lhs()[i].index;
 423                 scalar sl = R.lhs()[i].stoichCoeff;
 424                 dfdc[si][sj] += sl*kr;
 425             }
 426             forAll(R.rhs(), i)
 427             {
 428                 label si = R.rhs()[i].index;
 429                 scalar sr = R.rhs()[i].stoichCoeff;
 430                 dfdc[si][sj] -= sr*kr;
 431             }
 432         }
 433     }
 434
 435     // calculate the dcdT elements numerically
 436     scalar delta = 1.0e-8;
 437     scalarField dcdT0 = omega(c2, T - delta, p);
 438     scalarField dcdT1 = omega(c2, T + delta, p);
 439
 440     for (label i=0; i<nEqns(); i++)
 441     {
 442         dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
 443     }
 444
 445 }
 446
 447
 448 template<class CompType, class ThermoType>
 449 Foam::tmp<Foam::volScalarField>
 450 Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
 451 {
 452     scalar pf,cf,pr,cr;
 453     label lRef, rRef;
 454
 455     const volScalarField rho
 456     (
 457         IOobject
 458         (
 459             "rho",
 460             this->time().timeName(),
 461             this->mesh(),
 462             IOobject::NO_READ,
 463             IOobject::NO_WRITE,
 464             false
 465         ),
 466         this->thermo().rho()
 467     );
 468
 469     tmp<volScalarField> tsource
 470     (
 471         new volScalarField
 472         (
 473             IOobject
 474             (
 475                 "tc",
 476                 this->time().timeName(),
 477                 this->mesh(),
 478                 IOobject::NO_READ,
 479                 IOobject::NO_WRITE
 480             ),
 481             this->mesh(),
 482             dimensionedScalar("zero", dimTime, SMALL),
 483             zeroGradientFvPatchScalarField::typeName
 484         )
 485     );
 486
 487     scalarField& t = tsource();
 488
 489     label nReaction = reactions_.size();
 490
 491     if (this->chemistry_)
 492     {
 493         forAll(rho, celli)
 494         {
 495             scalar rhoi = rho[celli];
 496             scalar Ti = this->thermo().T()[celli];
 497             scalar pi = this->thermo().p()[celli];
 498             scalarField c(nSpecie_);
 499             scalar cSum = 0.0;
 500
 501             for (label i=0; i<nSpecie_; i++)
 502             {
 503                 scalar Yi = Y_[i][celli];
 504                 c[i] = rhoi*Yi/specieThermo_[i].W();
 505                 cSum += c[i];
 506             }
 507
 508             forAll(reactions_, i)
 509             {
 510                 const Reaction<ThermoType>& R = reactions_[i];
 511
 512                 omega
 513                 (
 514                     R, c, Ti, pi, pf, cf, lRef, pr, cr, rRef
 515                 );
 516
 517                 forAll(R.rhs(), s)
 518                 {
 519                     scalar sr = R.rhs()[s].stoichCoeff;
 520                     t[celli] += sr*pf*cf;
 521                 }
 522             }
 523             t[celli] = nReaction*cSum/t[celli];
 524         }
 525     }
 526
 527
 528     tsource().correctBoundaryConditions();
 529
 530     return tsource;
 531 }
 532
 533
 534 template<class CompType, class ThermoType>
 535 Foam::tmp<Foam::volScalarField>
 536 Foam::ODEChemistryModel<CompType, ThermoType>::dQ() const
 537 {
 538     tmp<volScalarField> tdQ
 539     (
 540         new volScalarField
 541         (
 542             IOobject
 543             (
 544                 "dQ",
 545                 this->mesh_.time().timeName(),
 546                 this->mesh_,
 547                 IOobject::NO_READ,
 548                 IOobject::NO_WRITE
 549             ),
 550             this->mesh_,
 551             dimensionedScalar
 552             (
 553                 "zero",
 554                 dimensionSet(1, -3, -1 , 0, 0, 0, 0),
 555                 0.0
 556             )
 557         )
 558     );
 559
 560     if (this->chemistry_)
 561     {
 562         scalarField& dQ = tdQ();
 563
 564         scalarField cp(dQ.size(), 0.0);
 565
 566         forAll(Y_, i)
 567         {
 568             forAll(dQ, cellI)
 569             {
 570                 scalar Ti = this->thermo().T()[cellI];
 571                 cp[cellI] += Y_[i][cellI]*specieThermo_[i].Cp(Ti);
 572                 scalar hi = specieThermo_[i].h(Ti);
 573                 dQ[cellI] -= hi*RR_[i][cellI];
 574             }
 575         }
 576
 577         dQ /= cp;
 578     }
 579
 580     return tdQ;
 581 }
 582
 583
 584 template<class CompType, class ThermoType>
 585 Foam::label Foam::ODEChemistryModel<CompType, ThermoType>::nEqns() const
 586 {
 587     // nEqns = number of species + temperature + pressure
 588     return nSpecie_ + 2;
 589 }
 590
 591
 592 template<class CompType, class ThermoType>
 593 void Foam::ODEChemistryModel<CompType, ThermoType>::calculate()
 594 {
 595     const volScalarField rho
 596     (
 597         IOobject
 598         (
 599             "rho",
 600             this->time().timeName(),
 601             this->mesh(),
 602             IOobject::NO_READ,
 603             IOobject::NO_WRITE,
 604             false
 605         ),
 606         this->thermo().rho()
 607     );
 608
 609     for (label i=0; i<nSpecie_; i++)
 610     {
 611         RR_[i].setSize(rho.size());
 612     }
 613
 614     if (this->chemistry_)
 615     {
 616         forAll(rho, celli)
 617         {
 618             for (label i=0; i<nSpecie_; i++)
 619             {
 620                 RR_[i][celli] = 0.0;
 621             }
 622
 623             scalar rhoi = rho[celli];
 624             scalar Ti = this->thermo().T()[celli];
 625             scalar pi = this->thermo().p()[celli];
 626
 627             scalarField c(nSpecie_);
 628             scalarField dcdt(nEqns(), 0.0);
 629
 630             for (label i=0; i<nSpecie_; i++)
 631             {
 632                 scalar Yi = Y_[i][celli];
 633                 c[i] = rhoi*Yi/specieThermo_[i].W();
 634             }
 635
 636             dcdt = omega(c, Ti, pi);
 637
 638             for (label i=0; i<nSpecie_; i++)
 639             {
 640                 RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
 641             }
 642         }
 643     }
 644 }
 645
 646
 647 template<class CompType, class ThermoType>
 648 Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 649 (
 650     const scalar t0,
 651     const scalar deltaT
 652 )
 653 {
 654     const volScalarField rho
 655     (
 656         IOobject
 657         (
 658             "rho",
 659             this->time().timeName(),
 660             this->mesh(),
 661             IOobject::NO_READ,
 662             IOobject::NO_WRITE,
 663             false
 664         ),
 665         this->thermo().rho()
 666     );
 667
 668     for (label i=0; i<nSpecie_; i++)
 669     {
 670         RR_[i].setSize(rho.size());
 671     }
 672
 673     if (!this->chemistry_)
 674     {
 675         return GREAT;
 676     }
 677
 678     scalar deltaTMin = GREAT;
 679
 680     forAll(rho, celli)
 681     {
 682         for (label i=0; i<nSpecie_; i++)
 683         {
 684             RR_[i][celli] = 0.0;
 685         }
 686
 687         scalar rhoi = rho[celli];
 688         scalar Ti = this->thermo().T()[celli];
 689         scalar hi = this->thermo().h()[celli];
 690         scalar pi = this->thermo().p()[celli];
 691
 692         scalarField c(nSpecie_);
 693         scalarField c0(nSpecie_);
 694         scalarField dc(nSpecie_, 0.0);
 695
 696         for (label i=0; i<nSpecie_; i++)
 697         {
 698             c[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
 699         }
 700         c0 = c;
 701
 702         scalar t = t0;
 703         scalar tauC = this->deltaTChem_[celli];
 704         scalar dt = min(deltaT, tauC);
 705         scalar timeLeft = deltaT;
 706
 707         // calculate the chemical source terms
 708         scalar cTot = 0.0;
 709
 710         while(timeLeft > SMALL)
 711         {
 712             tauC = solver().solve(c, Ti, pi, t, dt);
 713             t += dt;
 714
 715             // update the temperature
 716             cTot = sum(c);
 717             ThermoType mixture(0.0*specieThermo_[0]);
 718             for (label i=0; i<nSpecie_; i++)
 719             {
 720                 mixture += (c[i]/cTot)*specieThermo_[i];
 721             }
 722             Ti = mixture.TH(hi, Ti);
 723
 724             timeLeft -= dt;
 725             this->deltaTChem_[celli] = tauC;
 726             dt = min(timeLeft, tauC);
 727             dt = max(dt, SMALL);
 728         }
 729         deltaTMin = min(tauC, deltaTMin);
 730
 731         dc = c - c0;
 732         scalar WTot = 0.0;
 733         for (label i=0; i<nSpecie_; i++)
 734         {
 735             WTot += c[i]*specieThermo_[i].W();
 736         }
 737         WTot /= cTot;
 738
 739         for (label i=0; i<nSpecie_; i++)
 740         {
 741             RR_[i][celli] = dc[i]*specieThermo_[i].W()/deltaT;
 742         }
 743     }
 744
 745     // Don't allow the time-step to change more than a factor of 2
 746     deltaTMin = min(deltaTMin, 2*deltaT);
 747
 748     return deltaTMin;
 749 }
 750
 751
 752 // ************************************************************************* //