src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.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 "reactingOneDim.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "directMappedPatchBase.H"
  29 #include "mapDistribute.H"
  30 #include "zeroGradientFvPatchFields.H"
  31 #include "surfaceInterpolate.H"
  32 #include "fvm.H"
  33 #include "fvcDiv.H"
  34 #include "fvcVolumeIntegrate.H"
  35 #include "fvMatrices.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41 namespace regionModels
  42 {
  43 namespace pyrolysisModels
  44 {
  45
  46 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  47
  48 defineTypeNameAndDebug(reactingOneDim, 0);
  49
  50 addToRunTimeSelectionTable(pyrolysisModel, reactingOneDim, mesh);
  51 addToRunTimeSelectionTable(pyrolysisModel, reactingOneDim, dictionary);
  52
  53 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
  54
  55 void reactingOneDim::readReactingOneDimControls()
  56 {
  57     const dictionary& solution = this->solution().subDict("SIMPLE");
  58     solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
  59     time_.controlDict().lookup("maxDi") >> maxDiff_;
  60
  61     coeffs().lookup("radFluxName") >> primaryRadFluxName_;
  62     coeffs().lookup("minimumDelta") >> minimumDelta_;
  63 }
  64
  65
  66 bool reactingOneDim::read()
  67 {
  68     if (pyrolysisModel::read())
  69     {
  70         readReactingOneDimControls();
  71         return true;
  72     }
  73     else
  74     {
  75         return false;
  76     }
  77 }
  78
  79
  80 bool reactingOneDim::read(const dictionary& dict)
  81 {
  82     if (pyrolysisModel::read(dict))
  83     {
  84         readReactingOneDimControls();
  85         return true;
  86     }
  87     else
  88     {
  89         return false;
  90     }
  91 }
  92
  93
  94 void reactingOneDim::updateQr()
  95 {
  96     // Retrieve field from coupled region using mapped boundary conditions
  97     QrCoupled_.correctBoundaryConditions();
  98
  99     // Update local Qr from coupled Qr field
 100     Qr_ == dimensionedScalar("zero", Qr_.dimensions(), 0.0);
 101     forAll(intCoupledPatchIDs_, i)
 102     {
 103         const label patchI = intCoupledPatchIDs_[i];
 104
 105         scalarField& Qrp = Qr_.boundaryField()[patchI];
 106
 107         // Qr is negative going out the solid
 108         // If the surface is emitting the radiative flux is set to zero
 109         Qrp = max(Qrp, scalar(0.0));
 110     }
 111
 112     // Propagate Qr through 1-D regions
 113     forAll(intCoupledPatchIDs_, i)
 114     {
 115         const label patchI = intCoupledPatchIDs_[i];
 116
 117         const scalarField& Qrp = Qr_.boundaryField()[patchI];
 118         const vectorField& Cf = regionMesh().Cf().boundaryField()[patchI];
 119
 120         forAll(Qrp, faceI)
 121         {
 122             const scalar Qr0 = Qrp[faceI];
 123             point Cf0 = Cf[faceI];
 124             const labelList& cells = boundaryFaceCells_[faceI];
 125             scalar kappaInt = 0.0;
 126             forAll(cells, k)
 127             {
 128                 const label cellI = cells[k];
 129                 const point& Cf1 = regionMesh().cellCentres()[cellI];
 130                 const scalar delta = mag(Cf1 - Cf0);
 131                 kappaInt += kappa_[cellI]*delta;
 132                 Qr_[cellI] = Qr0*exp(-kappaInt);
 133                 Cf0 = Cf1;
 134             }
 135         }
 136     }
 137
 138     Qr_.correctBoundaryConditions();
 139 }
 140
 141
 142 void reactingOneDim::updatePhiGas()
 143 {
 144     phiHsGas_ ==  dimensionedScalar("zero", phiHsGas_.dimensions(), 0.0);
 145     phiGas_ == dimensionedScalar("zero", phiGas_.dimensions(), 0.0);
 146
 147     const speciesTable& gasTable = solidChemistry_->gasTable();
 148
 149     forAll(gasTable, gasI)
 150     {
 151         tmp<volScalarField> tHsiGas = solidChemistry_->gasHs(T_, gasI);
 152         tmp<volScalarField> tRRiGas = solidChemistry_->RRg(gasI);
 153
 154         const volScalarField& HsiGas = tHsiGas();
 155         const volScalarField& RRiGas = tRRiGas();
 156
 157         const surfaceScalarField HsiGasf(fvc::interpolate(HsiGas));
 158         const surfaceScalarField RRiGasf(fvc::interpolate(RRiGas));
 159
 160         forAll(intCoupledPatchIDs_, i)
 161         {
 162             const label patchI = intCoupledPatchIDs_[i];
 163             const scalarField& phiGasp = phiHsGas_.boundaryField()[patchI];
 164
 165             forAll(phiGasp, faceI)
 166             {
 167                 const labelList& cells = boundaryFaceCells_[faceI];
 168                 scalar massInt = 0.0;
 169                 forAllReverse(cells, k)
 170                 {
 171                     const label cellI = cells[k];
 172                     massInt += RRiGas[cellI]*regionMesh().V()[cellI];
 173                     phiHsGas_[cellI] += massInt*HsiGas[cellI];
 174                 }
 175
 176                 phiGas_.boundaryField()[patchI][faceI] += massInt;
 177
 178                 if (debug)
 179                 {
 180                     Info<< " Gas : " << gasTable[gasI]
 181                         << " on patch : " << patchI
 182                         << " mass produced at face(local) : "
 183                         <<  faceI
 184                         << " is : " << massInt
 185                         << " [kg/s] " << endl;
 186                 }
 187             }
 188         }
 189         tHsiGas().clear();
 190     }
 191 }
 192
 193
 194 void reactingOneDim::updateFields()
 195 {
 196     updateQr();
 197
 198     updatePhiGas();
 199 }
 200
 201
 202 void reactingOneDim::updateMesh(const scalarField& mass0)
 203 {
 204     if (!moveMesh_)
 205     {
 206         return;
 207     }
 208
 209     const scalarField newV(mass0/rho_);
 210
 211     Info<< "Initial/final volumes = " << gSum(regionMesh().V()) << ", "
 212         << gSum(newV) << " [m3]" << endl;
 213
 214     // move the mesh
 215     const labelList moveMap = moveMesh(regionMesh().V() - newV, minimumDelta_);
 216
 217     // flag any cells that have not moved as non-reacting
 218     forAll(moveMap, i)
 219     {
 220         if (moveMap[i] == 0)
 221         {
 222             solidChemistry_->setCellReacting(i, false);
 223         }
 224     }
 225 }
 226
 227
 228 void reactingOneDim::solveContinuity()
 229 {
 230     if (debug)
 231     {
 232         Info<< "reactingOneDim::solveContinuity()" << endl;
 233     }
 234
 235     solve
 236     (
 237         fvm::ddt(rho_)
 238      ==
 239       - solidChemistry_->RRg()
 240     );
 241 }
 242
 243
 244 void reactingOneDim::solveSpeciesMass()
 245 {
 246     if (debug)
 247     {
 248         Info<< "reactingOneDim::solveSpeciesMass()" << endl;
 249     }
 250
 251     volScalarField Yt(0.0*Ys_[0]);
 252
 253     for (label i=0; i<Ys_.size()-1; i++)
 254     {
 255         volScalarField& Yi = Ys_[i];
 256
 257         fvScalarMatrix YiEqn
 258         (
 259             fvm::ddt(rho_, Yi)
 260          ==
 261             solidChemistry_->RRs(i)
 262         );
 263
 264         if (moveMesh_)
 265         {
 266             surfaceScalarField phiRhoMesh
 267             (
 268                 fvc::interpolate(Yi*rho_)*regionMesh().phi()
 269             );
 270
 271             YiEqn -= fvc::div(phiRhoMesh);
 272         }
 273
 274         YiEqn.solve(regionMesh().solver("Yi"));
 275         Yi.max(0.0);
 276         Yt += Yi;
 277     }
 278
 279     Ys_[Ys_.size() - 1] = 1.0 - Yt;
 280 }
 281
 282
 283 void reactingOneDim::solveEnergy()
 284 {
 285     if (debug)
 286     {
 287         Info<< "reactingOneDim::solveEnergy()" << endl;
 288     }
 289
 290     const volScalarField rhoCp(rho_*solidThermo_.Cp());
 291
 292     const surfaceScalarField phiQr(fvc::interpolate(Qr_)*nMagSf());
 293
 294     const surfaceScalarField phiGas(fvc::interpolate(phiHsGas_));
 295
 296     fvScalarMatrix TEqn
 297     (
 298         fvm::ddt(rhoCp, T_)
 299       - fvm::laplacian(K_, T_)
 300      ==
 301         chemistrySh_
 302       + fvc::div(phiQr)
 303       + fvc::div(phiGas)
 304     );
 305
 306     if (moveMesh_)
 307     {
 308         surfaceScalarField phiMesh
 309         (
 310             fvc::interpolate(rhoCp*T_)*regionMesh().phi()
 311         );
 312
 313         TEqn -= fvc::div(phiMesh);
 314     }
 315
 316     TEqn.relax();
 317     TEqn.solve();
 318
 319     Info<< "pyrolysis min/max(T) = " << min(T_).value() << ", "
 320         << max(T_).value() << endl;
 321 }
 322
 323
 324 void reactingOneDim::calculateMassTransfer()
 325 {
 326     totalGasMassFlux_ = 0;
 327     forAll(intCoupledPatchIDs_, i)
 328     {
 329         const label patchI = intCoupledPatchIDs_[i];
 330         totalGasMassFlux_ += gSum(phiGas_.boundaryField()[patchI]);
 331     }
 332
 333     if (infoOutput_)
 334     {
 335         totalHeatRR_ = fvc::domainIntegrate(chemistrySh_);
 336
 337         addedGasMass_ +=
 338             fvc::domainIntegrate(solidChemistry_->RRg())*time_.deltaT();
 339         lostSolidMass_ +=
 340             fvc::domainIntegrate(solidChemistry_->RRs())*time_.deltaT();
 341     }
 342 }
 343
 344
 345 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 346
 347 reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
 348 :
 349     pyrolysisModel(modelType, mesh),
 350     solidChemistry_(solidChemistryModel::New(regionMesh())),
 351     solidThermo_(solidChemistry_->solidThermo()),
 352     kappa_(solidThermo_.kappa()),
 353     K_(solidThermo_.K()),
 354     rho_(solidThermo_.rho()),
 355     Ys_(solidThermo_.composition().Y()),
 356     T_(solidThermo_.T()),
 357     primaryRadFluxName_(coeffs().lookupOrDefault<word>("radFluxName", "Qr")),
 358     nNonOrthCorr_(-1),
 359     maxDiff_(10),
 360     minimumDelta_(1e-4),
 361
 362     phiGas_
 363     (
 364         IOobject
 365         (
 366             "phiGas",
 367             time().timeName(),
 368             regionMesh(),
 369             IOobject::NO_READ,
 370             IOobject::AUTO_WRITE
 371         ),
 372         regionMesh(),
 373         dimensionedScalar("zero", dimMass/dimTime, 0.0)
 374     ),
 375
 376     phiHsGas_
 377     (
 378         IOobject
 379         (
 380             "phiHsGas",
 381             time().timeName(),
 382             regionMesh(),
 383             IOobject::NO_READ,
 384             IOobject::NO_WRITE
 385         ),
 386         regionMesh(),
 387         dimensionedScalar("zero", dimEnergy/dimTime, 0.0)
 388     ),
 389
 390     chemistrySh_
 391     (
 392         IOobject
 393         (
 394             "chemistrySh",
 395             time().timeName(),
 396             regionMesh(),
 397             IOobject::NO_READ,
 398             IOobject::AUTO_WRITE
 399         ),
 400         regionMesh(),
 401         dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
 402     ),
 403
 404     QrCoupled_
 405     (
 406         IOobject
 407         (
 408             primaryRadFluxName_,
 409             time().timeName(),
 410             regionMesh(),
 411             IOobject::MUST_READ,
 412             IOobject::AUTO_WRITE
 413         ),
 414         regionMesh()
 415     ),
 416
 417     Qr_
 418     (
 419         IOobject
 420         (
 421             "QrPyr",
 422             time().timeName(),
 423             regionMesh(),
 424             IOobject::NO_READ,
 425             IOobject::AUTO_WRITE
 426         ),
 427         regionMesh(),
 428         dimensionedScalar("zero", dimEnergy/dimArea/dimTime, 0.0),
 429         zeroGradientFvPatchVectorField::typeName
 430     ),
 431
 432     lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
 433     addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
 434     totalGasMassFlux_(0.0),
 435     totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0))
 436 {
 437     if (active_)
 438     {
 439         read();
 440     }
 441 }
 442
 443
 444 reactingOneDim::reactingOneDim
 445 (
 446     const word& modelType,
 447     const fvMesh& mesh,
 448     const dictionary& dict
 449 )
 450 :
 451     pyrolysisModel(modelType, mesh, dict),
 452     solidChemistry_(solidChemistryModel::New(regionMesh())),
 453     solidThermo_(solidChemistry_->solidThermo()),
 454     kappa_(solidThermo_.kappa()),
 455     K_(solidThermo_.K()),
 456     rho_(solidThermo_.rho()),
 457     Ys_(solidThermo_.composition().Y()),
 458     T_(solidThermo_.T()),
 459     primaryRadFluxName_(dict.lookupOrDefault<word>("radFluxName", "Qr")),
 460     nNonOrthCorr_(-1),
 461     maxDiff_(10),
 462     minimumDelta_(1e-4),
 463
 464     phiGas_
 465     (
 466         IOobject
 467         (
 468             "phiGas",
 469             time().timeName(),
 470             regionMesh(),
 471             IOobject::NO_READ,
 472             IOobject::AUTO_WRITE
 473         ),
 474         regionMesh(),
 475         dimensionedScalar("zero", dimMass/dimTime, 0.0)
 476     ),
 477
 478     phiHsGas_
 479     (
 480         IOobject
 481         (
 482             "phiHsGas",
 483             time().timeName(),
 484             regionMesh(),
 485             IOobject::NO_READ,
 486             IOobject::NO_WRITE
 487         ),
 488         regionMesh(),
 489         dimensionedScalar("zero", dimEnergy/dimTime, 0.0)
 490     ),
 491
 492     chemistrySh_
 493     (
 494         IOobject
 495         (
 496             "chemistrySh",
 497             time().timeName(),
 498             regionMesh(),
 499             IOobject::NO_READ,
 500             IOobject::AUTO_WRITE
 501         ),
 502         regionMesh(),
 503         dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
 504     ),
 505
 506     QrCoupled_
 507     (
 508         IOobject
 509         (
 510             primaryRadFluxName_,
 511             time().timeName(),
 512             regionMesh(),
 513             IOobject::MUST_READ,
 514             IOobject::AUTO_WRITE
 515         ),
 516         regionMesh()
 517     ),
 518
 519     Qr_
 520     (
 521         IOobject
 522         (
 523             "QrPyr",
 524             time().timeName(),
 525             regionMesh(),
 526             IOobject::NO_READ,
 527             IOobject::AUTO_WRITE
 528         ),
 529         regionMesh(),
 530         dimensionedScalar("zero", dimEnergy/dimArea/dimTime, 0.0),
 531         zeroGradientFvPatchVectorField::typeName
 532     ),
 533
 534     lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
 535     addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
 536     totalGasMassFlux_(0.0),
 537     totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0))
 538 {
 539     if (active_)
 540     {
 541         read(dict);
 542     }
 543 }
 544
 545
 546 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 547
 548 reactingOneDim::~reactingOneDim()
 549 {}
 550
 551
 552 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 553
 554 scalar reactingOneDim::addMassSources(const label patchI, const label faceI)
 555 {
 556     label index = 0;
 557     forAll(primaryPatchIDs_, i)
 558     {
 559         if (primaryPatchIDs_[i] == patchI)
 560         {
 561             index = i;
 562             break;
 563         }
 564     }
 565
 566     const label localPatchId =  intCoupledPatchIDs_[index];
 567
 568     const scalar massAdded = phiGas_.boundaryField()[localPatchId][faceI];
 569
 570     if (debug)
 571     {
 572         Info<< "\nPyrolysis region: " << type() << "added mass : "
 573             << massAdded << endl;
 574     }
 575
 576     return massAdded;
 577 }
 578
 579
 580 scalar reactingOneDim::solidRegionDiffNo() const
 581 {
 582     scalar DiNum = 0.0;
 583
 584     if (regionMesh().nInternalFaces() > 0)
 585     {
 586         surfaceScalarField KrhoCpbyDelta
 587         (
 588             regionMesh().surfaceInterpolation::deltaCoeffs()
 589           * fvc::interpolate(K_)
 590           / fvc::interpolate(Cp()*rho_)
 591         );
 592
 593         DiNum = max(KrhoCpbyDelta.internalField())*time_.deltaTValue();
 594     }
 595
 596     return DiNum;
 597 }
 598
 599
 600 scalar reactingOneDim::maxDiff() const
 601 {
 602     return maxDiff_;
 603 }
 604
 605
 606 const volScalarField& reactingOneDim::rho() const
 607 {
 608     return rho_;
 609 }
 610
 611
 612 const volScalarField& reactingOneDim::T() const
 613 {
 614     return T_;
 615 }
 616
 617
 618 const tmp<volScalarField> reactingOneDim::Cp() const
 619 {
 620     return solidThermo_.Cp();
 621 }
 622
 623
 624 const volScalarField& reactingOneDim::kappa() const
 625 {
 626     return kappa_;
 627 }
 628
 629
 630 const volScalarField& reactingOneDim::K() const
 631 {
 632     return K_;
 633 }
 634
 635
 636 const surfaceScalarField& reactingOneDim::phiGas() const
 637 {
 638     return phiGas_;
 639 }
 640
 641
 642 void reactingOneDim::preEvolveRegion()
 643 {
 644     pyrolysisModel::preEvolveRegion();
 645
 646     // Initialise all cells as able to react
 647     forAll(T_, cellI)
 648     {
 649         solidChemistry_->setCellReacting(cellI, true);
 650     }
 651
 652     // De-activate reactions if pyrolysis region coupled to (valid) film
 653     if (filmCoupled_)
 654     {
 655         const volScalarField& filmDelta = filmDeltaPtr_();
 656
 657         forAll(intCoupledPatchIDs_, i)
 658         {
 659             const label patchI = intCoupledPatchIDs_[i];
 660             const scalarField& filmDeltap = filmDelta.boundaryField()[patchI];
 661
 662             forAll(filmDeltap, faceI)
 663             {
 664                 const scalar filmDelta0 = filmDeltap[faceI];
 665                 if (filmDelta0 > reactionDeltaMin_)
 666                 {
 667                     const labelList& cells = boundaryFaceCells_[faceI];
 668
 669                     // TODO: only limit cell adjacent to film?
 670                     //solidChemistry_->setCellNoReacting(cells[0])
 671
 672                     // Propagate flag through 1-D region
 673                     forAll(cells, k)
 674                     {
 675                         solidChemistry_->setCellReacting(cells[k], false);
 676                     }
 677                 }
 678             }
 679         }
 680     }
 681 }
 682
 683
 684 void reactingOneDim::evolveRegion()
 685 {
 686
 687     Info<< "\nEvolving pyrolysis in region: " << regionMesh().name() << endl;
 688     const scalarField mass0 = rho_*regionMesh().V();
 689
 690     solidChemistry_->solve
 691     (
 692         time().value() - time().deltaTValue(),
 693         time().deltaTValue()
 694     );
 695
 696     solveContinuity();
 697
 698     updateMesh(mass0);
 699
 700     chemistrySh_ = solidChemistry_->Sh()();
 701
 702     updateFields();
 703
 704     solveSpeciesMass();
 705
 706     for (int nonOrth=0; nonOrth<=nNonOrthCorr_; nonOrth++)
 707     {
 708         solveEnergy();
 709     }
 710
 711     calculateMassTransfer();
 712
 713     solidThermo_.correct();
 714 }
 715
 716
 717 void reactingOneDim::info() const
 718 {
 719     Info<< "\nPyrolysis in region: " << regionMesh().name() << endl;
 720
 721     Info<< indent << "Total gas mass produced  [kg] = "
 722         << addedGasMass_.value() << nl
 723         << indent << "Total solid mass lost    [kg] = "
 724         << lostSolidMass_.value() << nl
 725         << indent << "Total pyrolysis gases  [kg/s] = "
 726         << totalGasMassFlux_ << nl
 727         << indent << "Total heat release rate [J/s] = "
 728         << totalHeatRR_.value() << nl;
 729 }
 730
 731
 732 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 733
 734 } // End namespace Foam
 735 } // End namespace regionModels
 736 } // End namespace pyrolysisModels
 737
 738 // ************************************************************************* //