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