src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "ReactingMultiphaseParcel.H"
  27 #include "mathematicalConstants.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 template<class ParcelType>
  32 const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::GAS(0);
  33
  34 template<class ParcelType>
  35 const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::LIQ(1);
  36
  37 template<class ParcelType>
  38 const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::SLD(2);
  39
  40
  41 // * * * * * * * * * * * *  Private Member Functions * * * * * * * * * * * * //
  42
  43 template<class ParcelType>
  44 template<class TrackData>
  45 Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::cpEff
  46 (
  47     TrackData& td,
  48     const scalar p,
  49     const scalar T,
  50     const label idG,
  51     const label idL,
  52     const label idS
  53 ) const
  54 {
  55     return
  56         this->Y_[GAS]*td.cloud().composition().cp(idG, YGas_, p, T)
  57       + this->Y_[LIQ]*td.cloud().composition().cp(idL, YLiquid_, p, T)
  58       + this->Y_[SLD]*td.cloud().composition().cp(idS, YSolid_, p, T);
  59 }
  60
  61
  62 template<class ParcelType>
  63 template<class TrackData>
  64 Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
  65 (
  66     TrackData& td,
  67     const scalar p,
  68     const scalar T,
  69     const label idG,
  70     const label idL,
  71     const label idS
  72 ) const
  73 {
  74     return
  75         this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
  76       + this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
  77       + this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
  78 }
  79
  80
  81 template<class ParcelType>
  82 template<class TrackData>
  83 Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::LEff
  84 (
  85     TrackData& td,
  86     const scalar p,
  87     const scalar T,
  88     const label idG,
  89     const label idL,
  90     const label idS
  91 ) const
  92 {
  93     return
  94         this->Y_[GAS]*td.cloud().composition().L(idG, YGas_, p, T)
  95       + this->Y_[LIQ]*td.cloud().composition().L(idL, YLiquid_, p, T)
  96       + this->Y_[SLD]*td.cloud().composition().L(idS, YSolid_, p, T);
  97 }
  98
  99
 100 template<class ParcelType>
 101 Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
 102 (
 103     const scalar mass0,
 104     const scalarField& dMassGas,
 105     const scalarField& dMassLiquid,
 106     const scalarField& dMassSolid
 107 )
 108 {
 109     scalarField& YMix = this->Y_;
 110
 111     scalar massGas =
 112         this->updateMassFraction(mass0*YMix[GAS], dMassGas, YGas_);
 113     scalar massLiquid =
 114         this->updateMassFraction(mass0*YMix[LIQ], dMassLiquid, YLiquid_);
 115     scalar massSolid =
 116         this->updateMassFraction(mass0*YMix[SLD], dMassSolid, YSolid_);
 117
 118     scalar massNew = max(massGas + massLiquid + massSolid, ROOTVSMALL);
 119
 120     YMix[GAS] = massGas/massNew;
 121     YMix[LIQ] = massLiquid/massNew;
 122     YMix[SLD] = 1.0 - YMix[GAS] - YMix[LIQ];
 123
 124     return massNew;
 125 }
 126
 127
 128 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 129
 130 template<class ParcelType>
 131 template<class TrackData>
 132 void Foam::ReactingMultiphaseParcel<ParcelType>::setCellValues
 133 (
 134     TrackData& td,
 135     const scalar dt,
 136     const label cellI
 137 )
 138 {
 139     ReactingParcel<ParcelType>::setCellValues(td, dt, cellI);
 140 }
 141
 142
 143 template<class ParcelType>
 144 template<class TrackData>
 145 void Foam::ReactingMultiphaseParcel<ParcelType>::cellValueSourceCorrection
 146 (
 147     TrackData& td,
 148     const scalar dt,
 149     const label cellI
 150 )
 151 {
 152     scalar massCell = this->massCell(cellI);
 153
 154     scalar addedMass = 0.0;
 155     forAll(td.cloud().rhoTrans(), i)
 156     {
 157         addedMass += td.cloud().rhoTrans(i)[cellI];
 158     }
 159
 160     this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
 161
 162     scalar massCellNew = massCell + addedMass;
 163     this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
 164
 165     scalar cpEff = 0;
 166     if (addedMass > ROOTVSMALL)
 167     {
 168         forAll(td.cloud().rhoTrans(), i)
 169         {
 170             scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
 171             cpEff +=
 172                 Y*td.cloud().mcCarrierThermo().speciesData()[i].Cp(this->Tc_);
 173         }
 174     }
 175     const scalar cpc = td.cpInterp().psi()[cellI];
 176     this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
 177
 178     this->Tc_ += td.cloud().hsTrans()[cellI]/(this->cpc_*massCellNew);
 179 }
 180
 181
 182 template<class ParcelType>
 183 template<class TrackData>
 184 void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 185 (
 186     TrackData& td,
 187     const scalar dt,
 188     const label cellI
 189 )
 190 {
 191     // Define local properties at beginning of timestep
 192     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 193     const scalar np0 = this->nParticle_;
 194     const scalar d0 = this->d_;
 195     const vector& U0 = this->U_;
 196     const scalar rho0 = this->rho_;
 197     const scalar T0 = this->T_;
 198     const scalar cp0 = this->cp_;
 199     const scalar mass0 = this->mass();
 200
 201     const scalar pc = this->pc_;
 202
 203     const scalarField& YMix = this->Y_;
 204     const label idG = td.cloud().composition().idGas();
 205     const label idL = td.cloud().composition().idLiquid();
 206     const label idS = td.cloud().composition().idSolid();
 207
 208
 209     // Calc surface values
 210     // ~~~~~~~~~~~~~~~~~~~
 211     scalar Ts, rhos, mus, Pr, kappa;
 212     ThermoParcel<ParcelType>::
 213         calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
 214
 215     // Reynolds number
 216     scalar Re = this->Re(U0, d0, rhos, mus);
 217
 218
 219     // Sources
 220     //~~~~~~~~
 221
 222     // Explicit momentum source for particle
 223     vector Su = vector::zero;
 224
 225     // Momentum transfer from the particle to the carrier phase
 226     vector dUTrans = vector::zero;
 227
 228     // Explicit enthalpy source for particle
 229     scalar Sh = 0.0;
 230
 231     // Sensible enthalpy transfer from the particle to the carrier phase
 232     scalar dhsTrans = 0.0;
 233
 234
 235     // Phase change in liquid phase
 236     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 237
 238     // Mass transfer due to phase change
 239     scalarField dMassPC(YLiquid_.size(), 0.0);
 240
 241     // Molar flux of species emitted from the particle (kmol/m^2/s)
 242     scalar Ne = 0.0;
 243
 244     // Sum Ni*Cpi*Wi of emission species
 245     scalar NCpW = 0.0;
 246
 247     // Surface concentrations of emitted species
 248     scalarField Cs(td.cloud().mcCarrierThermo().species().size(), 0.0);
 249
 250     // Calc mass and enthalpy transfer due to phase change
 251     this->calcPhaseChange
 252     (
 253         td,
 254         dt,
 255         cellI,
 256         Re,
 257         Ts,
 258         mus/rhos,
 259         d0,
 260         T0,
 261         mass0,
 262         idL,
 263         YMix[LIQ],
 264         YLiquid_,
 265         dMassPC,
 266         Sh,
 267         Ne,
 268         NCpW,
 269         Cs
 270     );
 271
 272
 273     // Devolatilisation
 274     // ~~~~~~~~~~~~~~~~
 275
 276     // Mass transfer due to devolatilisation
 277     scalarField dMassDV(YGas_.size(), 0.0);
 278
 279     // Calc mass and enthalpy transfer due to devolatilisation
 280     this->calcDevolatilisation
 281     (
 282         td,
 283         dt,
 284         Ts,
 285         d0,
 286         T0,
 287         mass0,
 288         this->mass0_,
 289         idG,
 290         YMix[GAS],
 291         YGas_,
 292         canCombust_,
 293         dMassDV,
 294         Sh,
 295         Ne,
 296         NCpW,
 297         Cs
 298     );
 299
 300     // Correct surface values due to emitted species
 301     this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Pr, kappa);
 302
 303
 304     // Surface reactions
 305     // ~~~~~~~~~~~~~~~~~
 306
 307     // Change in carrier phase composition due to surface reactions
 308     scalarField dMassSRGas(YGas_.size(), 0.0);
 309     scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
 310     scalarField dMassSRSolid(YSolid_.size(), 0.0);
 311     scalarField
 312         dMassSRCarrier
 313         (
 314             td.cloud().mcCarrierThermo().species().size(),
 315             0.0
 316         );
 317
 318     // Clac mass and enthalpy transfer due to surface reactions
 319     this->calcSurfaceReactions
 320     (
 321         td,
 322         dt,
 323         cellI,
 324         d0,
 325         T0,
 326         mass0,
 327         canCombust_,
 328         Ne,
 329         YMix,
 330         YGas_,
 331         YLiquid_,
 332         YSolid_,
 333         dMassSRGas,
 334         dMassSRLiquid,
 335         dMassSRSolid,
 336         dMassSRCarrier,
 337         Sh,
 338         dhsTrans
 339     );
 340
 341
 342     // Update component mass fractions
 343     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 344
 345     scalarField dMassGas = dMassDV + dMassSRGas;
 346     scalarField dMassLiquid = dMassPC + dMassSRLiquid;
 347     scalarField dMassSolid = dMassSRSolid;
 348
 349     scalar mass1 =
 350         updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
 351
 352
 353     // Heat transfer
 354     // ~~~~~~~~~~~~~
 355
 356     // Calculate new particle temperature
 357     scalar T1 =
 358         this->calcHeatTransfer
 359         (
 360             td,
 361             dt,
 362             cellI,
 363             Re,
 364             Pr,
 365             kappa,
 366             d0,
 367             rho0,
 368             T0,
 369             cp0,
 370             NCpW,
 371             Sh,
 372             dhsTrans
 373         );
 374
 375
 376     // Motion
 377     // ~~~~~~
 378
 379     // Calculate new particle velocity
 380     vector U1 =
 381         this->calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
 382
 383     dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
 384
 385
 386     // Accumulate carrier phase source terms
 387     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 388
 389     if (td.cloud().coupled())
 390     {
 391         // Transfer mass lost from particle to carrier mass source
 392         forAll(YGas_, i)
 393         {
 394             label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
 395             td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
 396         }
 397         forAll(YLiquid_, i)
 398         {
 399             label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
 400             td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
 401         }
 402 /*
 403         // No mapping between solid components and carrier phase
 404         forAll(YSolid_, i)
 405         {
 406             label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
 407             td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
 408         }
 409 */
 410         forAll(dMassSRCarrier, i)
 411         {
 412             td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
 413         }
 414
 415         // Update momentum transfer
 416         td.cloud().UTrans()[cellI] += np0*dUTrans;
 417
 418         // Update sensible enthalpy transfer
 419         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 420     }
 421
 422
 423     // Remove the particle when mass falls below minimum threshold
 424     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 425
 426     if (mass1 < td.constProps().minParticleMass())
 427     {
 428         td.keepParticle = false;
 429
 430         if (td.cloud().coupled())
 431         {
 432             // Absorb parcel into carrier phase
 433             forAll(YGas_, i)
 434             {
 435                 label gid =
 436                     td.cloud().composition().localToGlobalCarrierId(GAS, i);
 437                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
 438             }
 439             forAll(YLiquid_, i)
 440             {
 441                 label gid =
 442                     td.cloud().composition().localToGlobalCarrierId(LIQ, i);
 443                 td.cloud().rhoTrans(gid)[cellI] +=
 444                     np0*mass1*YMix[LIQ]*YLiquid_[i];
 445             }
 446 /*
 447             // No mapping between solid components and carrier phase
 448             forAll(YSolid_, i)
 449             {
 450                 label gid =
 451                     td.cloud().composition().localToGlobalCarrierId(SLD, i);
 452                 td.cloud().rhoTrans(gid)[cellI] +=
 453                     np0*mass1*YMix[SLD]*YSolid_[i];
 454             }
 455 */
 456             td.cloud().UTrans()[cellI] += np0*mass1*U1;
 457             td.cloud().hsTrans()[cellI] +=
 458                 np0*mass1*HEff(td, pc, T1, idG, idL, idS); // using total h
 459         }
 460     }
 461
 462
 463     // Set new particle properties
 464     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 465
 466     else
 467     {
 468         this->cp_ = cpEff(td, pc, T1, idG, idL, idS);
 469         this->T_ = T1;
 470         this->U_ = U1;
 471
 472         // Update particle density or diameter
 473         if (td.constProps().constantVolume())
 474         {
 475             this->rho_ = mass1/this->volume();
 476         }
 477         else
 478         {
 479             this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
 480         }
 481     }
 482 }
 483
 484
 485 template<class ParcelType>
 486 template<class TrackData>
 487 void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 488 (
 489     TrackData& td,
 490     const scalar dt,
 491     const scalar Ts,
 492     const scalar d,
 493     const scalar T,
 494     const scalar mass,
 495     const scalar mass0,
 496     const label idVolatile,
 497     const scalar YVolatileTot,
 498     const scalarField& YVolatile,
 499     bool& canCombust,
 500     scalarField& dMassDV,
 501     scalar& Sh,
 502     scalar& N,
 503     scalar& NCpW,
 504     scalarField& Cs
 505 ) const
 506 {
 507     // Check that model is active, and that the parcel temperature is
 508     // within necessary limits for devolatilisation to occur
 509     if
 510     (
 511         !td.cloud().devolatilisation().active()
 512      || T < td.constProps().Tvap()
 513     )
 514     {
 515         return;
 516     }
 517
 518     // Total mass of volatiles evolved
 519     const scalar dMassTot = td.cloud().devolatilisation().calculate
 520     (
 521         dt,
 522         mass0,
 523         mass,
 524         T,
 525         td.cloud().composition().YMixture0()[idVolatile],
 526         YVolatileTot,
 527         canCombust
 528     );
 529
 530     // Volatile mass transfer - equal components of each volatile specie
 531     dMassDV = YVolatile*dMassTot;
 532
 533     td.cloud().addToMassDevolatilisation(this->nParticle_*dMassTot);
 534
 535     Sh -= dMassTot*td.constProps().LDevol()/dt;
 536
 537     // Molar average molecular weight of carrier mix
 538     const scalar Wc = this->rhoc_*specie::RR()*this->Tc_/this->pc_;
 539
 540     // Update molar emissions
 541     forAll(dMassDV, i)
 542     {
 543         // Note: hardcoded gaseous diffusivities for now
 544         // TODO: add to carrier thermo
 545         const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
 546         const label id =
 547             td.cloud().composition().localToGlobalCarrierId(GAS, i);
 548         const scalar Cp = td.cloud().mcCarrierThermo().speciesData()[id].Cp(Ts);
 549         const scalar W = td.cloud().mcCarrierThermo().speciesData()[id].W();
 550         const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 551
 552         // Dab calc'd using API vapour mass diffusivity function
 553         const scalar Dab =
 554             3.6059e-3*(pow(1.8*Ts, 1.75))*sqrt(1.0/W + 1.0/Wc)/(this->pc_*beta);
 555
 556         N += Ni;
 557         NCpW += Ni*Cp*W;
 558         Cs[id] += Ni*d/(2.0*Dab);
 559      }
 560 }
 561
 562
 563 template<class ParcelType>
 564 template<class TrackData>
 565 void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
 566 (
 567     TrackData& td,
 568     const scalar dt,
 569     const label cellI,
 570     const scalar d,
 571     const scalar T,
 572     const scalar mass,
 573     const bool canCombust,
 574     const scalar N,
 575     const scalarField& YMix,
 576     const scalarField& YGas,
 577     const scalarField& YLiquid,
 578     const scalarField& YSolid,
 579     scalarField& dMassSRGas,
 580     scalarField& dMassSRLiquid,
 581     scalarField& dMassSRSolid,
 582     scalarField& dMassSRCarrier,
 583     scalar& Sh,
 584     scalar& dhsTrans
 585 ) const
 586 {
 587     // Check that model is active
 588     if (!td.cloud().surfaceReaction().active() || !canCombust)
 589     {
 590         return;
 591     }
 592
 593     // Update surface reactions
 594     const scalar hReaction = td.cloud().surfaceReaction().calculate
 595     (
 596         dt,
 597         cellI,
 598         d,
 599         T,
 600         this->Tc_,
 601         this->pc_,
 602         this->rhoc_,
 603         mass,
 604         YGas,
 605         YLiquid,
 606         YSolid,
 607         YMix,
 608         N,
 609         dMassSRGas,
 610         dMassSRLiquid,
 611         dMassSRSolid,
 612         dMassSRCarrier
 613     );
 614
 615     td.cloud().addToMassSurfaceReaction
 616     (
 617         this->nParticle_
 618        *(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))
 619     );
 620
 621     const scalar xsi = min(T/5000.0, 1.0);
 622     const scalar hRetentionCoeffMod =
 623         (1.0 - xsi*xsi)*td.constProps().hRetentionCoeff();
 624
 625     Sh += hRetentionCoeffMod*hReaction/dt;
 626
 627     dhsTrans += (1.0 - hRetentionCoeffMod)*hReaction;
 628 }
 629
 630
 631 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 632
 633 template <class ParcelType>
 634 Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 635 (
 636     const ReactingMultiphaseParcel<ParcelType>& p
 637 )
 638 :
 639     ReactingParcel<ParcelType>(p),
 640     YGas_(p.YGas_),
 641     YLiquid_(p.YLiquid_),
 642     YSolid_(p.YSolid_)
 643 {}
 644
 645
 646 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 647
 648 #include "ReactingMultiphaseParcelIO.C"
 649
 650 // ************************************************************************* //