src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.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 "InjectionModel.H"
  28 #include "mathematicalConstants.H"
  29 #include "meshTools.H"
  30
  31 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
  32
  33 template<class CloudType>
  34 void Foam::InjectionModel<CloudType>::readProps()
  35 {
  36     IOobject propsDictHeader
  37     (
  38         "injectionProperties",
  39         owner_.db().time().timeName(),
  40         "uniform"/cloud::prefix/owner_.name(),
  41         owner_.db(),
  42         IOobject::MUST_READ,
  43         IOobject::NO_WRITE,
  44         false
  45     );
  46
  47     if (propsDictHeader.headerOk())
  48     {
  49         const IOdictionary propsDict(propsDictHeader);
  50
  51         propsDict.readIfPresent("massInjected", massInjected_);
  52         propsDict.readIfPresent("nInjections", nInjections_);
  53         propsDict.readIfPresent("parcelsAddedTotal", parcelsAddedTotal_);
  54         propsDict.readIfPresent("timeStep0", timeStep0_);
  55     }
  56 }
  57
  58
  59 template<class CloudType>
  60 void Foam::InjectionModel<CloudType>::writeProps()
  61 {
  62     if (owner_.db().time().outputTime())
  63     {
  64         IOdictionary propsDict
  65         (
  66             IOobject
  67             (
  68                 "injectionProperties",
  69                 owner_.db().time().timeName(),
  70                 "uniform"/cloud::prefix/owner_.name(),
  71                 owner_.db(),
  72                 IOobject::NO_READ,
  73                 IOobject::NO_WRITE,
  74                 false
  75             )
  76         );
  77
  78         propsDict.add("massInjected", massInjected_);
  79         propsDict.add("nInjections", nInjections_);
  80         propsDict.add("parcelsAddedTotal", parcelsAddedTotal_);
  81         propsDict.add("timeStep0", timeStep0_);
  82
  83         propsDict.regIOobject::write();
  84     }
  85 }
  86
  87
  88 template<class CloudType>
  89 void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
  90 (
  91     const scalar time,
  92     label& newParcels,
  93     scalar& newVolume
  94 )
  95 {
  96     // Initialise values
  97     newParcels = 0;
  98     newVolume = 0.0;
  99
 100     // Return if not started injection event
 101     if (time < SOI_)
 102     {
 103         timeStep0_ = time;
 104         return;
 105     }
 106
 107     // Make times relative to SOI
 108     scalar t0 = timeStep0_ - SOI_;
 109     scalar t1 = time - SOI_;
 110
 111     // Number of parcels to inject
 112     newParcels = parcelsToInject(t0, t1);
 113
 114     // Volume of parcels to inject
 115     newVolume = volumeToInject(t0, t1);
 116
 117     // Hold previous time if no parcels, but non-zero volume fraction
 118     if ((newParcels == 0) && (newVolume > 0.0))
 119     {
 120         // hold value of timeStep0_
 121     }
 122     else
 123     {
 124         // advance value of timeStep0_
 125         timeStep0_ = time;
 126     }
 127 }
 128
 129
 130 template<class CloudType>
 131 void Foam::InjectionModel<CloudType>::findCellAtPosition
 132 (
 133     label& cellI,
 134     vector& position
 135 )
 136 {
 137     const vector p0 = position;
 138
 139     bool foundCell = false;
 140
 141     cellI = owner_.mesh().findCell(position);
 142
 143     if (cellI >= 0)
 144     {
 145         const vector& C = owner_.mesh().C()[cellI];
 146         position += SMALL*(C - position);
 147
 148         foundCell = owner_.mesh().pointInCell(position, cellI);
 149     }
 150     reduce(foundCell, orOp<bool>());
 151
 152     // Last chance - find nearest cell and try that one
 153     // - the point is probably on an edge
 154     if (!foundCell)
 155     {
 156         cellI = owner_.mesh().findNearestCell(position);
 157
 158         if (cellI >= 0)
 159         {
 160             const vector& C = owner_.mesh().C()[cellI];
 161             position += SMALL*(C - position);
 162
 163             foundCell = owner_.mesh().pointInCell(position, cellI);
 164         }
 165         reduce(foundCell, orOp<bool>());
 166     }
 167
 168     if (!foundCell)
 169     {
 170         FatalErrorIn
 171         (
 172             "Foam::InjectionModel<CloudType>::findCellAtPosition"
 173             "("
 174                 "label&, "
 175                 "vector&"
 176             ")"
 177         )<< "Cannot find parcel injection cell. "
 178          << "Parcel position = " << p0 << nl
 179          << abort(FatalError);
 180     }
 181 }
 182
 183
 184 template<class CloudType>
 185 Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
 186 (
 187     const label parcels,
 188     const scalar volume,
 189     const scalar diameter,
 190     const scalar rho
 191 )
 192 {
 193     scalar nP = 0.0;
 194     switch (parcelBasis_)
 195     {
 196         case pbMass:
 197         {
 198             nP = volume/volumeTotal_
 199                 *massTotal_/rho
 200                /(parcels*mathematicalConstant::pi/6.0*pow3(diameter));
 201             break;
 202         }
 203         case pbNumber:
 204         {
 205             nP = massTotal_/(rho*volumeTotal_*parcels);
 206             break;
 207         }
 208         default:
 209         {
 210             nP = 0.0;
 211             FatalErrorIn
 212             (
 213                 "Foam::scalar "
 214                 "Foam::InjectionModel<CloudType>::setNumberOfParticles"
 215                 "("
 216                 "    const label, "
 217                 "    const scalar, "
 218                 "    const scalar, "
 219                 "    const scalar, "
 220                 "    const scalar"
 221                 ")"
 222             )<< "Unknown parcelBasis type" << nl
 223              << exit(FatalError);
 224         }
 225     }
 226
 227     return nP;
 228 }
 229
 230
 231 template<class CloudType>
 232 void Foam::InjectionModel<CloudType>::postInjectCheck(const label parcelsAdded)
 233 {
 234     if (parcelsAdded > 0)
 235     {
 236         Pout<< nl
 237             << "--> Cloud: " << owner_.name() << nl
 238             << "    Added " << parcelsAdded
 239             << " new parcels" << nl << endl;
 240     }
 241
 242     // Increment total number of parcels added
 243     parcelsAddedTotal_ += parcelsAdded;
 244
 245     // Update time for start of next injection
 246     time0_ = owner_.db().time().value();
 247
 248     // Increment number of injections
 249     nInjections_++;
 250
 251     // Write current state to properties file
 252     writeProps();
 253 }
 254
 255
 256 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 257
 258 template<class CloudType>
 259 Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
 260 :
 261     dict_(dictionary::null),
 262     owner_(owner),
 263     coeffDict_(dictionary::null),
 264     SOI_(0.0),
 265     volumeTotal_(0.0),
 266     massTotal_(0.0),
 267     massInjected_(0.0),
 268     nInjections_(0),
 269     parcelsAddedTotal_(0),
 270     parcelBasis_(pbNumber),
 271     time0_(0.0),
 272     timeStep0_(0.0)
 273 {
 274     readProps();
 275 }
 276
 277
 278 template<class CloudType>
 279 Foam::InjectionModel<CloudType>::InjectionModel
 280 (
 281     const dictionary& dict,
 282     CloudType& owner,
 283     const word& type
 284 )
 285 :
 286     dict_(dict),
 287     owner_(owner),
 288     coeffDict_(dict.subDict(type + "Coeffs")),
 289     SOI_(readScalar(coeffDict_.lookup("SOI"))),
 290     volumeTotal_(0.0),
 291     massTotal_(dimensionedScalar(coeffDict_.lookup("massTotal")).value()),
 292     massInjected_(0.0),
 293     nInjections_(0),
 294     parcelsAddedTotal_(0),
 295     parcelBasis_(pbNumber),
 296     time0_(owner.db().time().value()),
 297     timeStep0_(0.0)
 298 {
 299     // Provide some info
 300     // - also serves to initialise mesh dimensions - needed for parallel runs
 301     //   due to lazy evaluation of valid mesh dimensions
 302     Info<< "    Constructing " << owner.mesh().nGeometricD() << "-D injection"
 303         << endl;
 304
 305     word parcelBasisType = coeffDict_.lookup("parcelBasisType");
 306     if (parcelBasisType == "mass")
 307     {
 308         parcelBasis_ = pbMass;
 309     }
 310     else if (parcelBasisType == "number")
 311     {
 312         parcelBasis_ = pbNumber;
 313     }
 314     else
 315     {
 316         FatalErrorIn
 317         (
 318             "Foam::InjectionModel<CloudType>::InjectionModel"
 319             "("
 320                 "const dictionary&, "
 321                 "CloudType&, "
 322                 "const word&"
 323             ")"
 324         )<< "parcelBasisType must be either 'number' or 'mass'" << nl
 325          << exit(FatalError);
 326     }
 327
 328     readProps();
 329 }
 330
 331
 332 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 333
 334 template<class CloudType>
 335 Foam::InjectionModel<CloudType>::~InjectionModel()
 336 {}
 337
 338
 339 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 340
 341 template<class CloudType>
 342 template<class TrackData>
 343 void Foam::InjectionModel<CloudType>::inject(TrackData& td)
 344 {
 345     if (!active())
 346     {
 347         return;
 348     }
 349
 350     const scalar time = owner_.db().time().value();
 351     const scalar carrierDt = owner_.db().time().deltaT().value();
 352     const polyMesh& mesh = owner_.mesh();
 353
 354     // Prepare for next time step
 355     label newParcels = 0;
 356     scalar newVolume = 0.0;
 357     prepareForNextTimeStep(time, newParcels, newVolume);
 358
 359     // Return if no parcels are required
 360     if (newParcels == 0)
 361     {
 362         postInjectCheck(0);
 363         return;
 364     }
 365
 366     // Duration of injection period during this timestep
 367     const scalar deltaT =
 368         max(0.0, min(carrierDt, min(time - SOI_, timeEnd() - time0_)));
 369
 370     // Pad injection time if injection starts during this timestep
 371     const scalar padTime = max(0.0, SOI_ - time0_);
 372
 373     // Introduce new parcels linearly across carrier phase timestep
 374     label parcelsAdded = 0;
 375     for (label parcelI=0; parcelI<newParcels; parcelI++)
 376     {
 377         if (validInjection(parcelI))
 378         {
 379             // Calculate the pseudo time of injection for parcel 'parcelI'
 380             scalar timeInj = time0_ + padTime + deltaT*parcelI/newParcels;
 381
 382             // Determine the injection position and owner cell
 383             label cellI = -1;
 384             vector pos = vector::zero;
 385             setPositionAndCell(parcelI, newParcels, timeInj, pos, cellI);
 386
 387             if (cellI > -1)
 388             {
 389                 // Lagrangian timestep
 390                 scalar dt = time - timeInj;
 391
 392                 // Apply corrections to position for 2-D cases
 393                 meshTools::constrainToMeshCentre(mesh, pos);
 394
 395                 // Create a new parcel
 396                 parcelType* pPtr = new parcelType(td.cloud(), pos, cellI);
 397
 398                 // Assign new parcel properties in injection model
 399                 setProperties(parcelI, newParcels, timeInj, *pPtr);
 400
 401                 // Check new parcel properties
 402                 td.cloud().checkParcelProperties(*pPtr, dt, fullyDescribed());
 403
 404                 // Apply correction to velocity for 2-D cases
 405                 meshTools::constrainDirection
 406                 (
 407                     mesh,
 408                     mesh.solutionD(),
 409                     pPtr->U()
 410                 );
 411
 412                 // Number of particles per parcel
 413                 pPtr->nParticle() =
 414                     setNumberOfParticles
 415                     (
 416                         newParcels,
 417                         newVolume,
 418                         pPtr->d(),
 419                         pPtr->rho()
 420                     );
 421
 422                 // Add the new parcel
 423                 td.cloud().addParticle(pPtr);
 424
 425                 massInjected_ += pPtr->nParticle()*pPtr->mass();
 426                 parcelsAdded++;
 427             }
 428         }
 429     }
 430
 431     postInjectCheck(parcelsAdded);
 432 }
 433
 434
 435 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 436
 437 #include "NewInjectionModel.C"
 438
 439 // ************************************************************************* //