src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C

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