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

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