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-2008 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>::updateCellQuantities
  34 (
  35     TrackData& td,
  36     const scalar dt,
  37     const label celli
  38 )
  39 {
  40     ThermoParcel<ParcelType>::updateCellQuantities(td, dt, celli);
  41
  42     pc_ = td.pInterp().interpolate(this->position(), celli);
  43 }
  44
  45
  46 template<class ParcelType>
  47 template<class TrackData>
  48 void Foam::ReactingParcel<ParcelType>::calcCoupled
  49 (
  50     TrackData& td,
  51     const scalar dt,
  52     const label celli
  53 )
  54 {
  55     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  56     // Define local properties at beginning of timestep
  57     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  58     const vector U0 = this->U_;
  59     const scalar mass0 = this->mass();
  60     const scalar np0 = this->nParticle_;
  61     const scalar T0 = this->T_;
  62
  63     // ~~~~~~~~~~~~~~~~~~~~~~~~~
  64     // Initialise transfer terms
  65     // ~~~~~~~~~~~~~~~~~~~~~~~~~
  66
  67     // Momentum transfer from the particle to the carrier phase
  68     vector dUTrans = vector::zero;
  69
  70     // Enthalpy transfer from the particle to the carrier phase
  71     scalar dhTrans = 0.0;
  72
  73     // Mass transfer from particle to carrier phase
  74     // - components exist in particle already
  75     scalarList dMassMT(td.cloud().gases().size(), 0.0);
  76
  77     // Mass transfer due to surface reactions from particle to carrier phase
  78     // - components do not necessarily exist in particle already
  79     scalarList dMassSR(td.cloud().gases().size(), 0.0);
  80
  81
  82     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  83     // Calculate velocity - update U
  84     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  85     scalar Cud = 0.0;
  86     const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
  87
  88
  89     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  90     // Calculate heat transfer - update T
  91     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  92     scalar htc = 0.0;
  93     scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
  94
  95
  96     // ~~~~~~~~~~~~~~~~~~~~~~~
  97     // Calculate mass transfer
  98     // ~~~~~~~~~~~~~~~~~~~~~~~
  99     calcMassTransfer(td, dt, T0, T1, dMassMT);
 100
 101
 102     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 103     // Calculate surface reactions
 104     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 105
 106     // Initialise enthalpy retention to zero
 107     scalar dhRet = 0.0;
 108
 109     calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
 110
 111     // New total mass
 112     const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
 113
 114     // Correct particle temperature to account for latent heat of
 115     // devolatilisation
 116     T1 -=
 117         td.constProps().Ldevol()
 118        *sum(dMassMT)
 119        /(0.5*(mass0 + mass1)*this->cp_);
 120
 121     // Add retained enthalpy from surface reaction to particle and remove
 122     // from gas
 123     T1 += dhRet/(0.5*(mass0 + mass1)*this->cp_);
 124     dhTrans -= dhRet;
 125
 126     // Correct dhTrans to account for enthalpy of evolved volatiles
 127     dhTrans +=
 128         sum(dMassMT)
 129        *td.cloud().composition().HGas(YGas_, 0.5*(T0 + T1));
 130
 131     // Correct dhTrans to account for enthalpy of consumed solids
 132     dhTrans +=
 133         sum(dMassSR)
 134        *td.cloud().composition().HSolid(YSolid_, 0.5*(T0 + T1));
 135
 136
 137     // ~~~~~~~~~~~~~~~~~~~~~~~
 138     // Accumulate source terms
 139     // ~~~~~~~~~~~~~~~~~~~~~~~
 140
 141     // Transfer mass lost from particle to carrier mass source
 142     forAll(dMassMT, i)
 143     {
 144         td.cloud().rhoTrans(i)[celli] += np0*(dMassMT[i] + dMassSR[i]);
 145     }
 146
 147     // Update momentum transfer
 148     td.cloud().UTrans()[celli] += np0*dUTrans;
 149
 150     // Accumulate coefficient to be applied in carrier phase momentum coupling
 151     td.cloud().UCoeff()[celli] += np0*mass0*Cud;
 152
 153     // Update enthalpy transfer
 154     // - enthalpy of lost solids already accounted for
 155     td.cloud().hTrans()[celli] += np0*dhTrans;
 156
 157     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
 158     td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
 159
 160
 161     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 162     // Remove the particle when mass falls below minimum threshold
 163     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 164     if (mass1 < td.constProps().minParticleMass())
 165     {
 166         td.keepParticle = false;
 167
 168         // Absorb particle(s) into carrier phase
 169         forAll(dMassMT, i)
 170         {
 171             td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
 172         }
 173         td.cloud().hTrans()[celli] +=
 174             np0*mass1
 175            *(
 176                 YMixture_[0]*td.cloud().composition().HGas(YGas_, T1)
 177               + YMixture_[2]*td.cloud().composition().HSolid(YSolid_, T1)
 178             );
 179         td.cloud().UTrans()[celli] += np0*mass1*U1;
 180     }
 181     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 182     // Set new particle properties
 183     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 184     else
 185     {
 186         this->U_ = U1;
 187         this->T_ = T1;
 188         this->cp_ =
 189             YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
 190           + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
 191           + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
 192
 193         // Update particle density or diameter
 194         if (td.cloud().massTransfer().changesVolume())
 195         {
 196             this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
 197         }
 198         else
 199         {
 200             this->rho_ = mass1/this->volume();
 201         }
 202     }
 203 }
 204
 205
 206 template<class ParcelType>
 207 template<class TrackData>
 208 void Foam::ReactingParcel<ParcelType>::calcUncoupled
 209 (
 210     TrackData& td,
 211     const scalar dt,
 212     const label celli
 213 )
 214 {
 215     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 216     // Define local properties at beginning of timestep
 217     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 218     const scalar T0 = this->T_;
 219     const scalar mass0 = this->mass();
 220     const scalar cp0 = this->cp_;
 221
 222     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 223     // Initialise transfer terms
 224     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 225
 226     // Momentum transfer from the particle to the carrier phase
 227     vector dUTrans = vector::zero;
 228
 229     // Enthalpy transfer from the particle to the carrier phase
 230     scalar dhTrans = 0.0;
 231
 232     // Mass transfer from particle to carrier phase
 233     // - components exist in particle already
 234     scalarList dMassMT(td.cloud().gases().size(), 0.0);
 235
 236     // Mass transfer due to surface reactions from particle to carrier phase
 237     // - components do not necessarily exist in particle already
 238     scalarList dMassSR(td.cloud().gases().size(), 0.0);
 239
 240
 241     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 242     // Calculate velocity - update U
 243     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 244     scalar Cud = 0.0;
 245     const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 246
 247
 248     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 249     // Calculate heat transfer - update T
 250     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 251     scalar htc = 0.0;
 252     scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 253
 254
 255     // ~~~~~~~~~~~~~~~~~~~~~~~
 256     // Calculate mass transfer
 257     // ~~~~~~~~~~~~~~~~~~~~~~~
 258     calcMassTransfer(td, dt, T0, T1, dMassMT);
 259
 260
 261     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 262     // Calculate surface reactions
 263     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 264
 265     // Initialise enthalpy retention to zero
 266     scalar dhRet = 0.0;
 267
 268     calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
 269
 270     // New total mass
 271     const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
 272
 273     // New specific heat capacity
 274     const scalar cp1 =
 275         YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
 276       + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
 277       + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
 278
 279     // Add retained enthalpy to particle
 280     T1 += dhRet/(mass0*0.5*(cp0 + cp1));
 281
 282     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 283     // Remove the particle when mass falls below minimum threshold
 284     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 285     if (mass1 < td.constProps().minParticleMass())
 286     {
 287         td.keepParticle = false;
 288     }
 289     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 290     // Set new particle properties
 291     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 292     else
 293     {
 294         this->U_ = U1;
 295         this->T_ = T1;
 296         this->cp_ = cp1;
 297
 298         // Update particle density or diameter
 299         if (td.cloud().massTransfer().changesVolume())
 300         {
 301             this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
 302         }
 303         else
 304         {
 305             this->rho_ = mass1/this->volume();
 306         }
 307     }
 308 }
 309
 310
 311 template<class ParcelType>
 312 template<class TrackData>
 313 void Foam::ReactingParcel<ParcelType>::calcMassTransfer
 314 (
 315     TrackData& td,
 316     const scalar dt,
 317     const scalar T0,
 318     const scalar T1,
 319     scalarList& dMassMT
 320 )
 321 {
 322     if (td.cloud().composition().YMixture0()[1]>SMALL)
 323     {
 324         notImplemented
 325         (
 326             "void Foam::ReactingParcel<ParcelType>::"
 327             "calcMassTransfer(...): no treatment currently "
 328             "available for particles containing liquid species"
 329         )
 330     }
 331
 332     // Check that model is active, and that the parcel temperature is
 333     // within necessary limits for mass transfer to occur
 334     if
 335     (
 336         !td.cloud().massTransfer().active()
 337      || this->T_<td.constProps().Tvap()
 338      || this->T_<td.constProps().Tbp()
 339     )
 340     {
 341         return;
 342     }
 343
 344     // Determine mass to add to carrier phase
 345     const scalar mass = this->mass();
 346     const scalar dMassTot = td.cloud().massTransfer().calculate
 347     (
 348         dt,
 349         mass0_,
 350         mass,
 351         td.cloud().composition().YMixture0(),
 352         YMixture_,
 353         T0,
 354         canCombust_
 355     );
 356
 357     // Update (total) mass fractions
 358     YMixture_[0] = (YMixture_[0]*mass - dMassTot)/(mass - dMassTot);
 359     YMixture_[1] = YMixture_[1]*mass/(mass - dMassTot);
 360     YMixture_[2] = 1.0 - YMixture_[0] - YMixture_[1];
 361
 362     // Add to cummulative mass transfer
 363     forAll (YGas_, i)
 364     {
 365         label id = td.cloud().composition().gasGlobalIds()[i];
 366
 367         // Mass transfer
 368         scalar volatileMass = YGas_[i]*dMassTot;
 369         dMassMT[id] += volatileMass;
 370     }
 371 }
 372
 373
 374 template<class ParcelType>
 375 template<class TrackData>
 376 void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
 377 (
 378     TrackData& td,
 379     const scalar dt,
 380     const label celli,
 381     const scalar T0,
 382     const scalar T1,
 383     const scalarList& dMassMT,
 384     scalarList& dMassSR,
 385     scalar& dhRet
 386 )
 387 {
 388     // Check that model is active
 389     if (!td.cloud().surfaceReaction().active() || !canCombust_)
 390     {
 391         return;
 392     }
 393
 394     // Update mass transfer(s)
 395     // - Also updates Y()'s
 396     td.cloud().surfaceReaction().calculate
 397     (
 398         dt,
 399         celli,
 400         this->d_,
 401         T0,
 402         T1,
 403         this->Tc_,
 404         pc_,
 405         this->rhoc_,
 406         this->mass(),
 407         dMassMT,
 408         YGas_,
 409         YLiquid_,
 410         YSolid_,
 411         YMixture_,
 412         dMassSR,
 413         dhRet
 414     );
 415 }
 416
 417
 418 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 419
 420 #include "ReactingParcelIO.C"
 421
 422 // ************************************************************************* //
 423