src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud_.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 <lagrangianIntermediate/KinematicCloud_.H>
  28 #include <lagrangianIntermediate/DispersionModel.H>
  29 #include <lagrangianIntermediate/DragModel.H>
  30 #include <lagrangianIntermediate/InjectionModel.H>
  31 #include <lagrangianIntermediate/WallInteractionModel.H>
  32
  33 #include <lagrangianIntermediate/IntegrationScheme.H>
  34
  35 #include <finiteVolume/interpolationCellPoint.H>
  36
  37 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
  38
  39 template<class ParcelType>
  40 void Foam::KinematicCloud<ParcelType>::setInjectorCellAndPosition
  41 (
  42     label& pCell,
  43     vector& pPosition
  44 )
  45 {
  46     const vector originalPosition = pPosition;
  47
  48     bool foundCell = false;
  49
  50     pCell = mesh_.findCell(pPosition);
  51
  52     if (pCell >= 0)
  53     {
  54         const vector& C = mesh_.C()[pCell];
  55         pPosition += 1.0e-6*(C - pPosition);
  56
  57         foundCell = mesh_.pointInCell
  58         (
  59             pPosition,
  60             pCell
  61         );
  62     }
  63     reduce(foundCell, orOp<bool>());
  64
  65     // Last chance - find nearest cell and try that one
  66     // - the point is probably on an edge
  67     if (!foundCell)
  68     {
  69         pCell =  mesh_.findNearestCell(pPosition);
  70
  71         if (pCell >= 0)
  72         {
  73             const vector& C = mesh_.C()[pCell];
  74             pPosition += 1.0e-6*(C - pPosition);
  75
  76             foundCell = mesh_.pointInCell
  77             (
  78                 pPosition,
  79                 pCell
  80             );
  81         }
  82         reduce(foundCell, orOp<bool>());
  83     }
  84
  85     if (!foundCell)
  86     {
  87         FatalErrorIn
  88         (
  89             "void KinematicCloud<ParcelType>::findInjectorCell"
  90             "(label&, vector&)"
  91         )<< "Cannot find parcel injection cell. "
  92          << "Parcel position = " << originalPosition << nl
  93          << abort(FatalError);
  94     }
  95 }
  96
  97
  98 template<class ParcelType>
  99 Foam::scalar Foam::KinematicCloud<ParcelType>::setNumberOfParticles
 100 (
 101     const label nParcels,
 102     const scalar pDiameter,
 103     const scalar pVolumeFraction,
 104     const scalar pRho,
 105     const scalar pVolume
 106 )
 107 {
 108     scalar nP = 0.0;
 109     switch (parcelBasis_)
 110     {
 111         case pbMass:
 112         {
 113             nP = pVolumeFraction*massTotal_/nParcels
 114                /(pRho*mathematicalConstant::pi/6.0*pow(pDiameter, 3));
 115             break;
 116         }
 117         case pbNumber:
 118         {
 119             nP = pVolumeFraction*massTotal_/(pRho*pVolume);
 120             break;
 121         }
 122         default:
 123         {
 124             nP = 0.0;
 125             FatalErrorIn
 126             (
 127                 "Foam::KinematicCloud<ParcelType>::setNumberOfParticles"
 128                 "(const label, const scalar, const scalar, const scalar, "
 129                 "const scalar)"
 130             )<< "Unknown parcelBasis type" << nl
 131              << exit(FatalError);
 132         }
 133     }
 134
 135     return nP;
 136 }
 137
 138
 139 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 140
 141 template<class ParcelType>
 142 Foam::KinematicCloud<ParcelType>::KinematicCloud
 143 (
 144     const word& cloudType,
 145     const volPointInterpolation& vpi,
 146     const volScalarField& rho,
 147     const volVectorField& U,
 148     const volScalarField& mu,
 149     const dimensionedVector& g
 150 )
 151 :
 152     Cloud<ParcelType>(rho.mesh(), cloudType, false),
 153     kinematicCloud(),
 154     cloudType_(cloudType),
 155     mesh_(rho.mesh()),
 156     vpi_(vpi),
 157     particleProperties_
 158     (
 159         IOobject
 160         (
 161             cloudType + "Properties",
 162             rho.mesh().time().constant(),
 163             rho.mesh(),
 164             IOobject::MUST_READ,
 165             IOobject::NO_WRITE
 166         )
 167     ),
 168     constProps_(particleProperties_),
 169     parcelTypeId_(readLabel(particleProperties_.lookup("parcelTypeId"))),
 170     coupled_(particleProperties_.lookup("coupled")),
 171     rndGen_(label(0)),
 172     time0_(this->db().time().value()),
 173     parcelBasisType_(particleProperties_.lookup("parcelBasisType")),
 174     parcelBasis_(pbNumber),
 175     massTotal_
 176     (
 177         dimensionedScalar(particleProperties_.lookup("massTotal")).value()
 178     ),
 179     massInjected_(0.0),
 180     rho_(rho),
 181     U_(U),
 182     mu_(mu),
 183     g_(g),
 184     interpolationSchemes_(particleProperties_.subDict("interpolationSchemes")),
 185     dispersionModel_
 186     (
 187         DispersionModel<KinematicCloud<ParcelType> >::New
 188         (
 189             particleProperties_,
 190             *this
 191         )
 192     ),
 193     dragModel_
 194     (
 195         DragModel<KinematicCloud<ParcelType> >::New
 196         (
 197             particleProperties_,
 198             *this
 199         )
 200     ),
 201     injectionModel_
 202     (
 203         InjectionModel<KinematicCloud<ParcelType> >::New
 204         (
 205             particleProperties_,
 206             *this
 207         )
 208     ),
 209     wallInteractionModel_
 210     (
 211         WallInteractionModel<KinematicCloud<ParcelType> >::New
 212         (
 213             particleProperties_,
 214             *this
 215         )
 216     ),
 217     UIntegrator_
 218     (
 219         vectorIntegrationScheme::New
 220         (
 221             "U",
 222             particleProperties_.subDict("integrationSchemes")
 223         )
 224     ),
 225     nInjections_(0),
 226     nParcelsAdded_(0),
 227     nParcelsAddedTotal_(0),
 228     UTrans_
 229     (
 230         IOobject
 231         (
 232             this->name() + "UTrans",
 233             this->db().time().timeName(),
 234             this->db(),
 235             IOobject::NO_READ,
 236             IOobject::NO_WRITE,
 237             false
 238         ),
 239         mesh_,
 240         dimensionedVector("zero", dimensionSet(1, 1, -1, 0, 0), vector::zero)
 241     ),
 242     UCoeff_
 243     (
 244         IOobject
 245         (
 246             this->name() + "UCoeff",
 247             this->db().time().timeName(),
 248             this->db(),
 249             IOobject::NO_READ,
 250             IOobject::NO_WRITE,
 251             false
 252         ),
 253         mesh_,
 254         dimensionedScalar("zero",  dimensionSet(1, 0, -1, 0, 0), 0.0)
 255     )
 256 {
 257     if (parcelBasisType_ == "mass")
 258     {
 259         parcelBasis_ = pbMass;
 260     }
 261     else if (parcelBasisType_ == "number")
 262     {
 263         parcelBasis_ = pbNumber;
 264     }
 265     else
 266     {
 267         FatalErrorIn
 268         (
 269             "Foam::KinematicCloud<ParcelType>::KinematicCloud"
 270             "(const word&, const volPointInterpolation&, const volScalarField&"
 271             ", const volVectorField&, const volScalarField&, const "
 272             "dimensionedVector&)"
 273         )<< "parcelBasisType must be either 'number' or 'mass'" << nl
 274          << exit(FatalError);
 275     }
 276 }
 277
 278
 279 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 280
 281 template<class ParcelType>
 282 Foam::KinematicCloud<ParcelType>::~KinematicCloud()
 283 {}
 284
 285
 286 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 287
 288 template<class ParcelType>
 289 void Foam::KinematicCloud<ParcelType>::resetSourceTerms()
 290 {
 291     UTrans_.field() = vector::zero;
 292     UCoeff_.field() = 0.0;
 293 }
 294
 295
 296 template<class ParcelType>
 297 void Foam::KinematicCloud<ParcelType>::evolve()
 298 {
 299     autoPtr<interpolation<scalar> > rhoInterpolator =
 300         interpolation<scalar>::New
 301         (
 302             interpolationSchemes_,
 303             vpi_,
 304             rho_
 305         );
 306
 307     autoPtr<interpolation<vector> > UInterpolator =
 308         interpolation<vector>::New
 309         (
 310             interpolationSchemes_,
 311             vpi_,
 312             U_
 313         );
 314
 315     autoPtr<interpolation<scalar> > muInterpolator =
 316         interpolation<scalar>::New
 317         (
 318             interpolationSchemes_,
 319             vpi_,
 320             mu_
 321         );
 322
 323     typename ParcelType::trackData td
 324     (
 325         *this,
 326         constProps_,
 327         rhoInterpolator(),
 328         UInterpolator(),
 329         muInterpolator(),
 330         g_.value()
 331     );
 332
 333     inject(td);
 334
 335     if (coupled_)
 336     {
 337         resetSourceTerms();
 338     }
 339
 340     Cloud<ParcelType>::move(td);
 341 }
 342
 343
 344 template<class ParcelType>
 345 template<class TrackingData>
 346 void Foam::KinematicCloud<ParcelType>::inject
 347 (
 348     TrackingData& td
 349 )
 350 {
 351     scalar time = this->db().time().value();
 352
 353     scalar pRho = td.constProps().rho0();
 354
 355     this->injection().prepareForNextTimeStep(time0_, time);
 356
 357     // Number of parcels to introduce during this timestep
 358     const label nParcels = this->injection().nParcels();
 359
 360     // Return if no parcels are required
 361     if (!nParcels)
 362     {
 363         this->postInjectCheck();
 364         return;
 365     }
 366
 367     // Volume of particles to introduce during this timestep
 368     scalar pVolume = this->injection().volume();
 369
 370     // Volume fraction to introduce during this timestep
 371     scalar pVolumeFraction = this->injection().volumeFraction();
 372
 373     // Duration of injection period during this timestep
 374     scalar deltaT = min
 375     (
 376         this->db().time().deltaT().value(),
 377         min
 378         (
 379             time - this->injection().timeStart(),
 380             this->injection().timeEnd() - time0_
 381         )
 382     );
 383
 384     // Pad injection time if injection starts during this timestep
 385     scalar padTime = max
 386     (
 387         0.0,
 388         this->injection().timeStart() - time0_
 389     );
 390
 391     // Introduce new parcels linearly with time
 392     for (label iParcel=0; iParcel<nParcels; iParcel++)
 393     {
 394         // Calculate the pseudo time of injection for parcel 'iParcel'
 395         scalar timeInj = time0_ + padTime + deltaT*iParcel/nParcels;
 396
 397         // Determine injected parcel properties
 398         vector pPosition = this->injection().position
 399         (
 400             iParcel,
 401             timeInj,
 402             this->meshInfo()
 403         );
 404
 405         // Diameter of parcels
 406         scalar pDiameter = this->injection().d0(iParcel, timeInj);
 407
 408         // Number of particles per parcel
 409         scalar pNumberOfParticles = setNumberOfParticles
 410         (
 411             nParcels,
 412             pDiameter,
 413             pVolumeFraction,
 414             pRho,
 415             pVolume
 416         );
 417
 418         // Velocity of parcels
 419         vector pU = this->injection().velocity
 420         (
 421             iParcel,
 422             timeInj,
 423             this->meshInfo()
 424         );
 425
 426         // Determine the injection cell
 427         label pCell = -1;
 428         setInjectorCellAndPosition(pCell, pPosition);
 429
 430         if (pCell >= 0)
 431         {
 432             // construct the parcel that is to be injected
 433             ParcelType* pPtr = new ParcelType
 434             (
 435                 td.cloud(),
 436                 parcelTypeId_,
 437                 pPosition,
 438                 pCell,
 439                 pDiameter,
 440                 pU,
 441                 pNumberOfParticles,
 442                 td.constProps()
 443             );
 444
 445             scalar dt = time - timeInj;
 446
 447             pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
 448                 /this->time().deltaT().value();
 449
 450             this->injectParcel(td, pPtr);
 451          }
 452     }
 453
 454     this->postInjectCheck();
 455
 456     if (debug)
 457     {
 458         this->dumpParticlePositions();
 459     }
 460 }
 461
 462
 463 template<class ParcelType>
 464 template<class TrackingData>
 465 void Foam::KinematicCloud<ParcelType>::injectParcel
 466 (
 467     TrackingData& td,
 468     ParcelType* p
 469 )
 470 {
 471     addParticle(p);
 472     nParcelsAdded_++;
 473     nParcelsAddedTotal_++;
 474     massInjected_ += p->mass()*p->nParticle();
 475 }
 476
 477
 478 template<class ParcelType>
 479 void Foam::KinematicCloud<ParcelType>::postInjectCheck()
 480 {
 481     if (nParcelsAdded_)
 482     {
 483         Pout<< "\n--> Cloud: " << this->name() << nl
 484             << "    Added " << nParcelsAdded_
 485             <<  " new parcels" << nl << endl;
 486     }
 487
 488     // Reset parcel counters
 489     nParcelsAdded_ = 0;
 490
 491     // Set time for start of next injection
 492     time0_ = this->db().time().value();
 493
 494     // Increment number of injections
 495     nInjections_++;
 496 }
 497
 498
 499 template<class ParcelType>
 500 void Foam::KinematicCloud<ParcelType>::info() const
 501 {
 502     Info<< "Cloud name: " << this->name() << nl
 503         << "    Parcels added during this run   = "
 504         << returnReduce(nParcelsAddedTotal_, sumOp<label>()) << nl
 505         << "    Mass introduced during this run = "
 506         << returnReduce(massInjected_, sumOp<scalar>()) << nl
 507         << "    Current number of parcels       = "
 508         << returnReduce(this->size(), sumOp<label>()) << nl
 509         << "    Current mass in system          = "
 510         << returnReduce(massInSystem(), sumOp<scalar>()) << nl
 511         << endl;
 512 }
 513
 514
 515 template<class ParcelType>
 516 void Foam::KinematicCloud<ParcelType>::dumpParticlePositions() const
 517 {
 518     OFstream pObj
 519     (
 520         this->db().time().path()/"parcelPositions_"
 521       + this->name() + "_"
 522       + ::Foam::name(this->nInjections_) + ".obj"
 523     );
 524
 525     forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)
 526     {
 527         const ParcelType& p = iter();
 528         pObj<< "v " << p.position().x() << " " << p.position().y() << " "
 529             << p.position().z() << nl;
 530     }
 531
 532     pObj.flush();
 533 }
 534
 535
 536 // ************************ vim: set sw=4 sts=4 et: ************************ //