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