src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.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 "ReactingCloud.H"
  28 #include "CompositionModel.H"
  29 #include "MassTransferModel.H"
  30 #include "SurfaceReactionModel.H"
  31
  32 #include "interpolationCellPoint.H"
  33
  34 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  35
  36 template<class ParcelType>
  37 Foam::ReactingCloud<ParcelType>::ReactingCloud
  38 (
  39     const word& cloudType,
  40     const volPointInterpolation& vpi,
  41     const volScalarField& rho,
  42     const volVectorField& U,
  43     const dimensionedVector& g,
  44     hCombustionThermo& thermo,
  45     PtrList<specieReactingProperties>& gases
  46 )
  47 :
  48     ThermoCloud<ParcelType>(cloudType, vpi, rho, U, g, thermo),
  49     reactingCloud(),
  50     constProps_(this->particleProperties()),
  51     carrierThermo_(thermo),
  52     gases_(gases),
  53     compositionModel_
  54     (
  55         CompositionModel<ReactingCloud<ParcelType> >::New
  56         (
  57             this->particleProperties(),
  58             *this
  59         )
  60     ),
  61     massTransferModel_
  62     (
  63         MassTransferModel<ReactingCloud<ParcelType> >::New
  64         (
  65             this->particleProperties(),
  66             *this
  67         )
  68     ),
  69     surfaceReactionModel_
  70     (
  71         SurfaceReactionModel<ReactingCloud<ParcelType> >::New
  72         (
  73             this->particleProperties(),
  74             *this
  75         )
  76     ),
  77     rhoTrans_(thermo.composition().Y().size())
  78 {
  79     // Set storage for mass source fields and initialise to zero
  80     forAll(rhoTrans_, i)
  81     {
  82         rhoTrans_.set
  83         (
  84             i,
  85             new DimensionedField<scalar, volMesh>
  86             (
  87                 IOobject
  88                 (
  89                      this->name() + "rhoTrans" + Foam::name(i),
  90                      this->db().time().timeName(),
  91                      this->db(),
  92                      IOobject::NO_READ,
  93                      IOobject::NO_WRITE,
  94                      false
  95                 ),
  96                 this->mesh(),
  97                 dimensionedScalar("zero", dimMass, 0.0)
  98             )
  99         );
 100     }
 101 }
 102
 103
 104 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 105
 106 template<class ParcelType>
 107 Foam::ReactingCloud<ParcelType>::~ReactingCloud()
 108 {}
 109
 110
 111 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 112
 113 template<class ParcelType>
 114 void Foam::ReactingCloud<ParcelType>::resetSourceTerms()
 115 {
 116     ThermoCloud<ParcelType>::resetSourceTerms();
 117     forAll(rhoTrans_, i)
 118     {
 119         rhoTrans_[i].field() = 0.0;
 120     }
 121 }
 122
 123
 124 template<class ParcelType>
 125 void Foam::ReactingCloud<ParcelType>::evolve()
 126 {
 127     const volScalarField& T = carrierThermo_.T();
 128     const volScalarField cp = carrierThermo_.Cp();
 129     const volScalarField& p = carrierThermo_.p();
 130
 131     autoPtr<interpolation<scalar> > rhoInterpolator = interpolation<scalar>::New
 132     (
 133         this->interpolationSchemes(),
 134         this->vpi(),
 135         this->rho()
 136     );
 137
 138     autoPtr<interpolation<vector> > UInterpolator = interpolation<vector>::New
 139     (
 140         this->interpolationSchemes(),
 141         this->vpi(),
 142         this->U()
 143     );
 144
 145     autoPtr<interpolation<scalar> > muInterpolator = interpolation<scalar>::New
 146     (
 147         this->interpolationSchemes(),
 148         this->vpi(),
 149         this->mu()
 150     );
 151
 152     autoPtr<interpolation<scalar> > TInterpolator = interpolation<scalar>::New
 153     (
 154         this->interpolationSchemes(),
 155         this->vpi(),
 156         T
 157     );
 158
 159     autoPtr<interpolation<scalar> > cpInterpolator = interpolation<scalar>::New
 160     (
 161         this->interpolationSchemes(),
 162         this->vpi(),
 163         cp
 164     );
 165
 166     autoPtr<interpolation<scalar> > pInterpolator = interpolation<scalar>::New
 167     (
 168         this->interpolationSchemes(),
 169         this->vpi(),
 170         p
 171     );
 172
 173     typename ParcelType::trackData td
 174     (
 175         *this,
 176         constProps_,
 177         rhoInterpolator(),
 178         UInterpolator(),
 179         muInterpolator(),
 180         TInterpolator(),
 181         cpInterpolator(),
 182         pInterpolator(),
 183         this->g().value()
 184     );
 185
 186     inject(td);
 187
 188     if (this->coupled())
 189     {
 190         resetSourceTerms();
 191     }
 192
 193     Cloud<ParcelType>::move(td);
 194 }
 195
 196
 197 template<class ParcelType>
 198 template<class TrackingData>
 199 void Foam::ReactingCloud<ParcelType>::inject
 200 (
 201     TrackingData& td
 202 )
 203 {
 204     scalar time = this->db().time().value();
 205
 206     scalar pRho = td.constProps().rho0();
 207
 208     this->injection().prepareForNextTimeStep(this->time0(), time);
 209
 210     // Number of parcels to introduce during this timestep
 211     const label nParcels = this->injection().nParcels();
 212
 213     // Return if no parcels are required
 214     if (!nParcels)
 215     {
 216         this->postInjectCheck();
 217         return;
 218     }
 219
 220     // Volume of particles to introduce during this timestep
 221     scalar pVolume = this->injection().volume();
 222
 223     // Volume fraction to introduce during this timestep
 224     scalar pVolumeFraction = this->injection().volumeFraction();
 225
 226     // Duration of injection period during this timestep
 227     scalar deltaT = min
 228     (
 229         this->db().time().deltaT().value(),
 230         min
 231         (
 232             time - this->injection().timeStart(),
 233             this->injection().timeEnd() - this->time0()
 234         )
 235     );
 236
 237     // Pad injection time if injection starts during this timestep
 238     scalar padTime = max
 239     (
 240         0.0,
 241         this->injection().timeStart() - this->time0()
 242     );
 243
 244     // Introduce new parcels linearly with time
 245     for (label iParcel=0; iParcel<nParcels; iParcel++)
 246     {
 247         // Calculate the pseudo time of injection for parcel 'iParcel'
 248         scalar timeInj = this->time0() + padTime + deltaT*iParcel/nParcels;
 249
 250         // Determine injected parcel properties
 251         vector pPosition = this->injection().position
 252         (
 253             iParcel,
 254             timeInj,
 255             this->meshInfo()
 256         );
 257
 258         // Diameter of parcels
 259         scalar pDiameter = this->injection().d0(iParcel, timeInj);
 260
 261         // Number of particles per parcel
 262         scalar pNumberOfParticles = this->setNumberOfParticles
 263         (
 264             nParcels,
 265             pDiameter,
 266             pVolumeFraction,
 267             pRho,
 268             pVolume
 269         );
 270
 271         // Velocity of parcels
 272         vector pU = this->injection().velocity
 273         (
 274             iParcel,
 275             timeInj,
 276             this->meshInfo()
 277         );
 278
 279         // Determine the injection cell
 280         label pCell = -1;
 281         this->setInjectorCellAndPosition(pCell, pPosition);
 282
 283         if (pCell >= 0)
 284         {
 285             // construct the parcel that is to be injected
 286             ParcelType* pPtr = new ParcelType
 287             (
 288                 td.cloud(),
 289                 this->parcelTypeId(),
 290                 pPosition,
 291                 pCell,
 292                 pDiameter,
 293                 pU,
 294                 pNumberOfParticles,
 295                 composition().YGas0(),
 296                 composition().YLiquid0(),
 297                 composition().YSolid0(),
 298                 composition().YMixture0(),
 299                 td.constProps()
 300             );
 301
 302             scalar dt = time - timeInj;
 303
 304             pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
 305                 /this->db().time().deltaT().value();
 306
 307             this->injectParcel(td, pPtr);
 308          }
 309     }
 310
 311     this->postInjectCheck();
 312
 313     if (debug)
 314     {
 315         this->dumpParticlePositions();
 316     }
 317 }
 318
 319
 320 // ************************************************************************* //