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