src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.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 "ReactingParcel.H"
  27 #include "specie.H"
  28 #include "CompositionModel.H"
  29 #include "mathematicalConstants.H"
  30
  31 using namespace Foam::constant::mathematical;
  32
  33 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
  34
  35 template<class ParcelType>
  36 template<class TrackData>
  37 void Foam::ReactingParcel<ParcelType>::setCellValues
  38 (
  39     TrackData& td,
  40     const scalar dt,
  41     const label cellI
  42 )
  43 {
  44     ParcelType::setCellValues(td, dt, cellI);
  45
  46     pc_ = td.pInterp().interpolate
  47     (
  48         this->position(),
  49         this->currentTetIndices()
  50     );
  51
  52     if (pc_ < td.cloud().constProps().pMin())
  53     {
  54         if (debug)
  55         {
  56             WarningIn
  57             (
  58                 "void Foam::ReactingParcel<ParcelType>::setCellValues"
  59                 "("
  60                     "TrackData&, "
  61                     "const scalar, "
  62                     "const label"
  63                 ")"
  64             )   << "Limiting observed pressure in cell " << cellI << " to "
  65                 << td.cloud().constProps().pMin() <<  nl << endl;
  66         }
  67
  68         pc_ = td.cloud().constProps().pMin();
  69     }
  70 }
  71
  72
  73 template<class ParcelType>
  74 template<class TrackData>
  75 void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
  76 (
  77     TrackData& td,
  78     const scalar dt,
  79     const label cellI
  80 )
  81 {
  82     scalar addedMass = 0.0;
  83     forAll(td.cloud().rhoTrans(), i)
  84     {
  85         addedMass += td.cloud().rhoTrans(i)[cellI];
  86     }
  87
  88     if (addedMass < ROOTVSMALL)
  89     {
  90         return;
  91     }
  92
  93     const scalar massCell = this->massCell(cellI);
  94
  95     this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
  96
  97     const scalar massCellNew = massCell + addedMass;
  98     this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
  99
 100     scalar CpEff = 0.0;
 101     if (addedMass > ROOTVSMALL)
 102     {
 103         forAll(td.cloud().rhoTrans(), i)
 104         {
 105             scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
 106             CpEff += Y*td.cloud().composition().carrier().Cp(i, this->Tc_);
 107         }
 108     }
 109
 110     const scalar Cpc = td.CpInterp().psi()[cellI];
 111     this->Cpc_ = (massCell*Cpc + addedMass*CpEff)/massCellNew;
 112
 113     this->Tc_ += td.cloud().hsTrans()[cellI]/(this->Cpc_*massCellNew);
 114
 115     if (debug && (this->Tc_ < td.cloud().constProps().TMin()))
 116     {
 117         if (debug)
 118         {
 119             WarningIn
 120             (
 121                 "void Foam::ReactingParcel<ParcelType>::"
 122                 "cellValueSourceCorrection"
 123                 "("
 124                     "TrackData&, "
 125                     "const scalar, "
 126                     "const label"
 127                 ")"
 128             )   << "Limiting observed temperature in cell " << cellI << " to "
 129                 << td.cloud().constProps().TMin() <<  nl << endl;
 130         }
 131
 132         this->Tc_ = td.cloud().constProps().TMin();
 133     }
 134
 135 //  constant pressure
 136 //  this->pc_ = this->pc_;
 137 }
 138
 139
 140 template<class ParcelType>
 141 template<class TrackData>
 142 void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
 143 (
 144     TrackData& td,
 145     const label cellI,
 146     const scalar T,
 147     const scalarField& Cs,
 148     scalar& rhos,
 149     scalar& mus,
 150     scalar& Prs,
 151     scalar& kappas
 152 )
 153 {
 154     // No correction if total concentration of emitted species is small
 155     if (sum(Cs) < SMALL)
 156     {
 157         return;
 158     }
 159
 160     const SLGThermo& thermo = td.cloud().thermo();
 161
 162     // Far field carrier  molar fractions
 163     scalarField Xinf(td.cloud().thermo().carrier().species().size());
 164
 165     forAll(Xinf, i)
 166     {
 167         Xinf[i] = thermo.carrier().Y(i)[cellI]/thermo.carrier().W(i);
 168     }
 169     Xinf /= sum(Xinf);
 170
 171     // Molar fraction of far field species at particle surface
 172     const scalar Xsff = 1.0 - min(sum(Cs)*specie::RR*this->T_/pc_, 1.0);
 173
 174     // Surface carrier total molar concentration
 175     const scalar CsTot = pc_/(specie::RR*this->T_);
 176
 177     // Surface carrier composition (molar fraction)
 178     scalarField Xs(Xinf.size());
 179
 180     // Surface carrier composition (mass fraction)
 181     scalarField Ys(Xinf.size());
 182
 183     forAll(Xs, i)
 184     {
 185         // Molar concentration of species at particle surface
 186         const scalar Csi = Cs[i] + Xsff*Xinf[i]*CsTot;
 187
 188         Xs[i] = (2.0*Csi + Xinf[i]*CsTot)/3.0;
 189         Ys[i] = Xs[i]*thermo.carrier().W(i);
 190     }
 191     Xs /= sum(Xs);
 192     Ys /= sum(Ys);
 193
 194
 195     rhos = 0;
 196     mus = 0;
 197     kappas = 0;
 198     scalar Cps = 0;
 199     scalar sumYiSqrtW = 0;
 200     scalar sumYiCbrtW = 0;
 201
 202     forAll(Ys, i)
 203     {
 204         const scalar W = thermo.carrier().W(i);
 205         const scalar sqrtW = sqrt(W);
 206         const scalar cbrtW = cbrt(W);
 207
 208         rhos += Xs[i]*W;
 209         mus += Ys[i]*sqrtW*thermo.carrier().mu(i, T);
 210         kappas += Ys[i]*cbrtW*thermo.carrier().kappa(i, T);
 211         Cps += Xs[i]*thermo.carrier().Cp(i, T);
 212
 213         sumYiSqrtW += Ys[i]*sqrtW;
 214         sumYiCbrtW += Ys[i]*cbrtW;
 215     }
 216
 217     Cps = max(Cps, ROOTVSMALL);
 218
 219     rhos *= pc_/(specie::RR*T);
 220     rhos = max(rhos, ROOTVSMALL);
 221
 222     mus /= sumYiSqrtW;
 223     mus = max(mus, ROOTVSMALL);
 224
 225     kappas /= sumYiCbrtW;
 226     kappas = max(kappas, ROOTVSMALL);
 227
 228     Prs = Cps*mus/kappas;
 229 }
 230
 231
 232 template<class ParcelType>
 233 Foam::scalar Foam::ReactingParcel<ParcelType>::updateMassFraction
 234 (
 235     const scalar mass0,
 236     const scalarField& dMass,
 237     scalarField& Y
 238 ) const
 239 {
 240     scalar mass1 = mass0 - sum(dMass);
 241
 242     // only update the mass fractions if the new particle mass is finite
 243     if (mass1 > ROOTVSMALL)
 244     {
 245         forAll(Y, i)
 246         {
 247             Y[i] = (Y[i]*mass0 - dMass[i])/mass1;
 248         }
 249     }
 250
 251     return mass1;
 252 }
 253
 254
 255 template<class ParcelType>
 256 template<class TrackData>
 257 void Foam::ReactingParcel<ParcelType>::calc
 258 (
 259     TrackData& td,
 260     const scalar dt,
 261     const label cellI
 262 )
 263 {
 264     typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
 265     const CompositionModel<reactingCloudType>& composition =
 266         td.cloud().composition();
 267
 268
 269     // Define local properties at beginning of time step
 270     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 271
 272     const scalar np0 = this->nParticle_;
 273     const scalar d0 = this->d_;
 274     const vector& U0 = this->U_;
 275     const scalar T0 = this->T_;
 276     const scalar mass0 = this->mass();
 277
 278
 279     // Calc surface values
 280     scalar Ts, rhos, mus, Prs, kappas;
 281     this->calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Prs, kappas);
 282     scalar Res = this->Re(U0, d0, rhos, mus);
 283
 284
 285     // Sources
 286     //~~~~~~~~
 287
 288     // Explicit momentum source for particle
 289     vector Su = vector::zero;
 290
 291     // Linearised momentum source coefficient
 292     scalar Spu = 0.0;
 293
 294     // Momentum transfer from the particle to the carrier phase
 295     vector dUTrans = vector::zero;
 296
 297     // Explicit enthalpy source for particle
 298     scalar Sh = 0.0;
 299
 300     // Linearised enthalpy source coefficient
 301     scalar Sph = 0.0;
 302
 303     // Sensible enthalpy transfer from the particle to the carrier phase
 304     scalar dhsTrans = 0.0;
 305
 306
 307     // 1. Compute models that contribute to mass transfer - U, T held constant
 308     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 309
 310     // Phase change
 311     // ~~~~~~~~~~~~
 312
 313     // Mass transfer due to phase change
 314     scalarField dMassPC(Y_.size(), 0.0);
 315
 316     // Molar flux of species emitted from the particle (kmol/m^2/s)
 317     scalar Ne = 0.0;
 318
 319     // Sum Ni*Cpi*Wi of emission species
 320     scalar NCpW = 0.0;
 321
 322     // Surface concentrations of emitted species
 323     scalarField Cs(composition.carrier().species().size(), 0.0);
 324
 325     // Calc mass and enthalpy transfer due to phase change
 326     calcPhaseChange
 327     (
 328         td,
 329         dt,
 330         cellI,
 331         Res,
 332         Ts,
 333         mus/rhos,
 334         d0,
 335         T0,
 336         mass0,
 337         0,
 338         1.0,
 339         Y_,
 340         dMassPC,
 341         Sh,
 342         Ne,
 343         NCpW,
 344         Cs
 345     );
 346
 347
 348     // 2. Update the parcel properties due to change in mass
 349     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 350
 351     scalarField dMass(dMassPC);
 352     scalar mass1 = updateMassFraction(mass0, dMass, Y_);
 353
 354     this->Cp_ = composition.Cp(0, Y_, pc_, T0);
 355
 356     // Update particle density or diameter
 357     if (td.cloud().constProps().constantVolume())
 358     {
 359         this->rho_ = mass1/this->volume();
 360     }
 361     else
 362     {
 363         this->d_ = cbrt(mass1/this->rho_*6.0/pi);
 364     }
 365
 366     // Remove the particle when mass falls below minimum threshold
 367     if (np0*mass1 < td.cloud().constProps().minParticleMass())
 368     {
 369         td.keepParticle = false;
 370
 371         if (td.cloud().solution().coupled())
 372         {
 373             scalar dm = np0*mass1;
 374
 375             // Absorb parcel into carrier phase
 376             forAll(Y_, i)
 377             {
 378                 scalar dmi = dm*Y_[i];
 379                 label gid = composition.localToGlobalCarrierId(0, i);
 380                 scalar hs = composition.carrier().Hs(gid, T0);
 381
 382                 td.cloud().rhoTrans(gid)[cellI] += dmi;
 383                 td.cloud().hsTrans()[cellI] += dmi*hs;
 384             }
 385             td.cloud().UTrans()[cellI] += dm*U0;
 386
 387             td.cloud().addToMassPhaseChange(dm);
 388         }
 389
 390         return;
 391     }
 392
 393     // Correct surface values due to emitted species
 394     correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
 395     Res = this->Re(U0, this->d_, rhos, mus);
 396
 397
 398     // 3. Compute heat- and momentum transfers
 399     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 400
 401     // Heat transfer
 402     // ~~~~~~~~~~~~~
 403
 404     // Calculate new particle temperature
 405     this->T_ =
 406         this->calcHeatTransfer
 407         (
 408             td,
 409             dt,
 410             cellI,
 411             Res,
 412             Prs,
 413             kappas,
 414             NCpW,
 415             Sh,
 416             dhsTrans,
 417             Sph
 418         );
 419
 420     this->Cp_ = composition.Cp(0, Y_, pc_, T0);
 421
 422
 423     // Motion
 424     // ~~~~~~
 425
 426     // Calculate new particle velocity
 427     this->U_ =
 428         this->calcVelocity(td, dt, cellI, Res, mus, mass1, Su, dUTrans, Spu);
 429
 430
 431     // 4. Accumulate carrier phase source terms
 432     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 433
 434     if (td.cloud().solution().coupled())
 435     {
 436         // Transfer mass lost to carrier mass, momentum and enthalpy sources
 437         forAll(dMass, i)
 438         {
 439             scalar dm = np0*dMass[i];
 440             label gid = composition.localToGlobalCarrierId(0, i);
 441             scalar hs = composition.carrier().Hs(gid, T0);
 442
 443             td.cloud().rhoTrans(gid)[cellI] += dm;
 444             td.cloud().UTrans()[cellI] += dm*U0;
 445             td.cloud().hsTrans()[cellI] += dm*hs;
 446         }
 447
 448         // Update momentum transfer
 449         td.cloud().UTrans()[cellI] += np0*dUTrans;
 450         td.cloud().UCoeff()[cellI] += np0*Spu;
 451
 452         // Update sensible enthalpy transfer
 453         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 454         td.cloud().hsCoeff()[cellI] += np0*Sph;
 455     }
 456 }
 457
 458
 459 template<class ParcelType>
 460 template<class TrackData>
 461 void Foam::ReactingParcel<ParcelType>::calcPhaseChange
 462 (
 463     TrackData& td,
 464     const scalar dt,
 465     const label cellI,
 466     const scalar Re,
 467     const scalar Ts,
 468     const scalar nus,
 469     const scalar d,
 470     const scalar T,
 471     const scalar mass,
 472     const label idPhase,
 473     const scalar YPhase,
 474     const scalarField& YComponents,
 475     scalarField& dMassPC,
 476     scalar& Sh,
 477     scalar& N,
 478     scalar& NCpW,
 479     scalarField& Cs
 480 )
 481 {
 482     if
 483     (
 484         !td.cloud().phaseChange().active()
 485      || T < td.cloud().constProps().Tvap()
 486      || YPhase < SMALL
 487     )
 488     {
 489         return;
 490     }
 491
 492     typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
 493     const CompositionModel<reactingCloudType>& composition =
 494         td.cloud().composition();
 495
 496
 497     // Calculate mass transfer due to phase change
 498     td.cloud().phaseChange().calculate
 499     (
 500         dt,
 501         cellI,
 502         Re,
 503         d,
 504         nus,
 505         T,
 506         Ts,
 507         pc_,
 508         dMassPC
 509     );
 510
 511     // Limit phase change mass by availability of each specie
 512     dMassPC = min(mass*YPhase*YComponents, dMassPC);
 513
 514     const scalar dMassTot = sum(dMassPC);
 515
 516     // Add to cumulative phase change mass
 517     td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
 518
 519     forAll(dMassPC, i)
 520     {
 521         const label idc = composition.localToGlobalCarrierId(idPhase, i);
 522         const label idl = composition.globalIds(idPhase)[i];
 523
 524         const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
 525         Sh -= dMassPC[i]*dh/dt;
 526     }
 527
 528
 529     // Update molar emissions
 530     if (td.cloud().heatTransfer().BirdCorrection())
 531     {
 532         // Average molecular weight of carrier mix - assumes perfect gas
 533         const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
 534
 535
 536         forAll(dMassPC, i)
 537         {
 538             const label idc = composition.localToGlobalCarrierId(idPhase, i);
 539             const label idl = composition.globalIds(idPhase)[i];
 540
 541             const scalar Cp = composition.carrier().Cp(idc, Ts);
 542             const scalar W = composition.carrier().W(idc);
 543             const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
 544
 545             const scalar Dab =
 546                 composition.liquids().properties()[idl].D(pc_, Ts, Wc);
 547
 548             // Molar flux of species coming from the particle (kmol/m^2/s)
 549             N += Ni;
 550
 551             // Sum of Ni*Cpi*Wi of emission species
 552             NCpW += Ni*Cp*W;
 553
 554             // Concentrations of emission species
 555             Cs[idc] += Ni*d/(2.0*Dab);
 556         }
 557     }
 558 }
 559
 560
 561 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 562
 563 template<class ParcelType>
 564 Foam::ReactingParcel<ParcelType>::ReactingParcel
 565 (
 566     const ReactingParcel<ParcelType>& p
 567 )
 568 :
 569     ParcelType(p),
 570     mass0_(p.mass0_),
 571     Y_(p.Y_),
 572     pc_(p.pc_)
 573 {}
 574
 575
 576 template<class ParcelType>
 577 Foam::ReactingParcel<ParcelType>::ReactingParcel
 578 (
 579     const ReactingParcel<ParcelType>& p,
 580     const polyMesh& mesh
 581 )
 582 :
 583     ParcelType(p, mesh),
 584     mass0_(p.mass0_),
 585     Y_(p.Y_),
 586     pc_(p.pc_)
 587 {}
 588
 589
 590 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 591
 592 #include "ReactingParcelIO.C"
 593
 594 // ************************************************************************* //
 595