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) 1991-2009 OpenCFD Ltd.
   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
  29 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
  30
  31 template<class ParcelType>
  32 template<class TrackData>
  33 void Foam::ReactingParcel<ParcelType>::setCellValues
  34 (
  35     TrackData& td,
  36     const scalar dt,
  37     const label cellI
  38 )
  39 {
  40     ThermoParcel<ParcelType>::setCellValues(td, dt, cellI);
  41
  42     pc_ = td.pInterp().interpolate(this->position(), cellI);
  43     if (pc_ < td.constProps().pMin())
  44     {
  45         WarningIn
  46         (
  47             "void Foam::ReactingParcel<ParcelType>::setCellValues"
  48             "("
  49                 "TrackData&, "
  50                 "const scalar, "
  51                 "const label"
  52             ")"
  53         )   << "Limiting observed pressure in cell " << cellI << " to "
  54             << td.constProps().pMin() <<  nl << endl;
  55
  56         pc_ = td.constProps().pMin();
  57     }
  58 }
  59
  60
  61 template<class ParcelType>
  62 template<class TrackData>
  63 void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
  64 (
  65     TrackData& td,
  66     const scalar dt,
  67     const label cellI
  68 )
  69 {
  70     scalar massCell = this->massCell(cellI);
  71
  72     scalar addedMass = 0.0;
  73     forAll(td.cloud().rhoTrans(), i)
  74     {
  75         addedMass += td.cloud().rhoTrans(i)[cellI];
  76     }
  77
  78     this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
  79
  80     scalar massCellNew = massCell + addedMass;
  81     this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
  82
  83     scalar cpEff = 0;
  84     if (addedMass > ROOTVSMALL)
  85     {
  86         forAll(td.cloud().rhoTrans(), i)
  87         {
  88             scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
  89             cpEff +=
  90                 Y*td.cloud().mcCarrierThermo().speciesData()[i].Cp(this->Tc_);
  91         }
  92     }
  93     const scalar cpc = td.cpInterp().psi()[cellI];
  94     this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
  95
  96     this->Tc_ += td.cloud().hsTrans()[cellI]/(this->cpc_*massCellNew);
  97 }
  98
  99
 100 template<class ParcelType>
 101 Foam::scalar Foam::ReactingParcel<ParcelType>::updateMassFraction
 102 (
 103     const scalar mass0,
 104     const scalarField& dMass,
 105     scalarField& Y
 106 ) const
 107 {
 108     scalar mass1 = mass0 - sum(dMass);
 109
 110     // only update the mass fractions if the new particle mass is finite
 111     if (mass1 > ROOTVSMALL)
 112     {
 113         forAll(Y, i)
 114         {
 115             Y[i] = (Y[i]*mass0 - dMass[i])/mass1;
 116         }
 117     }
 118
 119     return mass1;
 120 }
 121
 122
 123 template<class ParcelType>
 124 template<class TrackData>
 125 void Foam::ReactingParcel<ParcelType>::calc
 126 (
 127     TrackData& td,
 128     const scalar dt,
 129     const label cellI
 130 )
 131 {
 132     // Define local properties at beginning of time step
 133     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 134     const scalar np0 = this->nParticle_;
 135     const scalar d0 = this->d_;
 136     const vector& U0 = this->U_;
 137     const scalar rho0 = this->rho_;
 138     const scalar T0 = this->T_;
 139     const scalar cp0 = this->cp_;
 140     const scalar mass0 = this->mass();
 141
 142     // Explicit momentum source for particle
 143     vector Su = vector::zero;
 144
 145     // Momentum transfer from the particle to the carrier phase
 146     vector dUTrans = vector::zero;
 147
 148     // Explicit enthalpy source for particle
 149     scalar Sh = 0.0;
 150
 151     // Sensible enthalpy transfer from the particle to the carrier phase
 152     scalar dhsTrans = 0.0;
 153
 154
 155     // Phase change
 156     // ~~~~~~~~~~~~
 157
 158     // Mass transfer due to phase change
 159     scalarField dMassPC(Y_.size(), 0.0);
 160
 161     // Calc mass and enthalpy transfer due to phase change
 162     calcPhaseChange
 163     (
 164         td,
 165         dt,
 166         cellI,
 167         d0,
 168         T0,
 169         U0,
 170         mass0,
 171         0,
 172         1.0,
 173         Y_,
 174         dMassPC,
 175         Sh,
 176         dhsTrans
 177     );
 178
 179     // Update particle component mass and mass fractions
 180     scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
 181
 182
 183     // Heat transfer
 184     // ~~~~~~~~~~~~~
 185
 186     // Calculate new particle temperature
 187     scalar T1 =
 188         calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh, dhsTrans);
 189
 190
 191     // Motion
 192     // ~~~~~~
 193
 194     // Calculate new particle velocity
 195     vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
 196
 197     dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
 198
 199     // Accumulate carrier phase source terms
 200     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 201     if (td.cloud().coupled())
 202     {
 203         // Transfer mass lost from particle to carrier mass source
 204         forAll(dMassPC, i)
 205         {
 206             label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
 207             td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
 208             td.cloud().hcTrans()[cellI] +=
 209                 np0
 210                *dMassPC[i]
 211                *td.cloud().mcCarrierThermo().speciesData()[gid].H(T0);
 212         }
 213
 214         // Update momentum transfer
 215         td.cloud().UTrans()[cellI] += np0*dUTrans;
 216
 217         // Update sensible enthalpy transfer
 218         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 219     }
 220
 221
 222     // Remove the particle when mass falls below minimum threshold
 223     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 224     if (mass1 < td.constProps().minParticleMass())
 225     {
 226         td.keepParticle = false;
 227
 228         if (td.cloud().coupled())
 229         {
 230             // Absorb parcel into carrier phase
 231             forAll(Y_, i)
 232             {
 233                 label gid =
 234                     td.cloud().composition().localToGlobalCarrierId(0, i);
 235                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
 236             }
 237             td.cloud().UTrans()[cellI] += np0*mass1*U1;
 238             td.cloud().hcTrans()[cellI] +=
 239                 np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
 240         }
 241     }
 242
 243
 244     // Set new particle properties
 245     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 246
 247     else
 248     {
 249         this->cp_ = td.cloud().composition().cp(0, Y_, pc_, T1);
 250         this->T_ = T1;
 251         this->U_ = U1;
 252
 253         // Update particle density or diameter
 254         if (td.constProps().constantVolume())
 255         {
 256             this->rho_ = mass1/this->volume();
 257         }
 258         else
 259         {
 260             this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
 261         }
 262     }
 263 }
 264
 265
 266 template<class ParcelType>
 267 template<class TrackData>
 268 void Foam::ReactingParcel<ParcelType>::calcPhaseChange
 269 (
 270     TrackData& td,
 271     const scalar dt,
 272     const label cellI,
 273     const scalar d,
 274     const scalar T,
 275     const vector& U,
 276     const scalar mass,
 277     const label idPhase,
 278     const scalar YPhase,
 279     const scalarField& YComponents,
 280     scalarField& dMassPC,
 281     scalar& Sh,
 282     scalar& dhsTrans
 283 )
 284 {
 285     if
 286     (
 287         !td.cloud().phaseChange().active()
 288      || T < td.constProps().Tvap()
 289      || YPhase < SMALL
 290     )
 291     {
 292         return;
 293     }
 294
 295     // Calculate mass transfer due to phase change
 296     td.cloud().phaseChange().calculate
 297     (
 298         dt,
 299         cellI,
 300         d,
 301         min(T, td.constProps().Tbp()), // Limit to boiling temperature
 302         pc_,
 303         this->Tc_,
 304         this->muc_/(this->rhoc_ + ROOTVSMALL),
 305         U - this->Uc_,
 306         dMassPC
 307     );
 308
 309     // Limit phase change mass by availability of each specie
 310     dMassPC = min(mass*YPhase*YComponents, dMassPC);
 311
 312     scalar dMassTot = sum(dMassPC);
 313
 314     // Add to cumulative phase change mass
 315     td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
 316
 317     // Enthalphy transfer to carrier phase
 318     label id;
 319     forAll(YComponents, i)
 320     {
 321         id = td.cloud().composition().localToGlobalCarrierId(idPhase, i);
 322         const scalar hv = td.cloud().mcCarrierThermo().speciesData()[id].H(T);
 323
 324         id = td.cloud().composition().globalIds(idPhase)[i];
 325         const scalar hl =
 326             td.cloud().composition().liquids().properties()[id].h(pc_, T);
 327
 328         Sh += dMassPC[i]*(hl - hv)/dt;
 329     }
 330 }
 331
 332
 333 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 334
 335 template <class ParcelType>
 336 Foam::ReactingParcel<ParcelType>::ReactingParcel
 337 (
 338     const ReactingParcel<ParcelType>& p
 339 )
 340 :
 341     ThermoParcel<ParcelType>(p),
 342     mass0_(p.mass0_),
 343     Y_(p.Y_),
 344     pc_(p.pc_)
 345 {}
 346
 347
 348 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 349
 350 #include "ReactingParcelIO.C"
 351
 352 // ************************************************************************* //
 353