InjectionModel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "InjectionModel.H"
27 #include "mathematicalConstants.H"
28 #include "meshTools.H"
29 #include "volFields.H"
30 
31 using namespace Foam::constant::mathematical;
32 
33 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
34 
35 template<class CloudType>
37 (
38  const scalar time,
39  label& newParcels,
40  scalar& newVolumeFraction
41 )
42 {
43  // Initialise values
44  newParcels = 0;
45  newVolumeFraction = 0.0;
46  bool validInjection = false;
47 
48  // Return if not started injection event
49  if (time < SOI_)
50  {
51  timeStep0_ = time;
52  return validInjection;
53  }
54 
55  // Make times relative to SOI
56  scalar t0 = timeStep0_ - SOI_;
57  scalar t1 = time - SOI_;
58 
59  // Number of parcels to inject
60  newParcels = this->parcelsToInject(t0, t1);
61 
62  // Volume of parcels to inject
63  newVolumeFraction =
64  this->volumeToInject(t0, t1)
65  /(volumeTotal_ + ROOTVSMALL);
66 
67  if (newVolumeFraction > 0)
68  {
69  if (newParcels > 0)
70  {
71  timeStep0_ = time;
72  validInjection = true;
73  }
74  else
75  {
76  // injection should have started, but not sufficient volume to
77  // produce (at least) 1 parcel - hold value of timeStep0_
78  validInjection = false;
79  }
80  }
81  else
82  {
83  timeStep0_ = time;
84  validInjection = false;
85  }
86 
87  return validInjection;
88 }
89 
90 
91 template<class CloudType>
93 (
94  label& cellI,
95  label& tetFaceI,
96  label& tetPtI,
97  vector& position,
98  bool errorOnNotFound
99 )
100 {
101  const volVectorField& cellCentres = this->owner().mesh().C();
102 
103  const vector p0 = position;
104 
105  this->owner().mesh().findCellFacePt
106  (
107  position,
108  cellI,
109  tetFaceI,
110  tetPtI
111  );
112 
113  label procI = -1;
114 
115  if (cellI >= 0)
116  {
117  procI = Pstream::myProcNo();
118  }
119 
120  reduce(procI, maxOp<label>());
121 
122  // Ensure that only one processor attempts to insert this Parcel
123 
124  if (procI != Pstream::myProcNo())
125  {
126  cellI = -1;
127  tetFaceI = -1;
128  tetPtI = -1;
129  }
130 
131  // Last chance - find nearest cell and try that one - the point is
132  // probably on an edge
133  if (procI == -1)
134  {
135  cellI = this->owner().mesh().findNearestCell(position);
136 
137  if (cellI >= 0)
138  {
139  position += SMALL*(cellCentres[cellI] - position);
140 
141  if (this->owner().mesh().pointInCell(position, cellI))
142  {
143  procI = Pstream::myProcNo();
144  }
145  }
146 
147  reduce(procI, maxOp<label>());
148 
149  if (procI != Pstream::myProcNo())
150  {
151  cellI = -1;
152  tetFaceI = -1;
153  tetPtI = -1;
154  }
155  }
156 
157  if (procI == -1)
158  {
159  if (errorOnNotFound)
160  {
162  (
163  "Foam::InjectionModel<CloudType>::findCellAtPosition"
164  "("
165  "label&, "
166  "label&, "
167  "label&, "
168  "vector&, "
169  "bool"
170  ")"
171  ) << "Cannot find parcel injection cell. "
172  << "Parcel position = " << p0 << nl
173  << abort(FatalError);
174  }
175  else
176  {
177  return false;
178  }
179  }
180 
181  return true;
182 }
183 
184 
185 template<class CloudType>
187 (
188  const label parcels,
189  const scalar volumeFraction,
190  const scalar diameter,
191  const scalar rho
192 )
193 {
194  scalar nP = 0.0;
195  switch (parcelBasis_)
196  {
197  case pbMass:
198  {
199  scalar volumep = pi/6.0*pow3(diameter);
200  scalar volumeTot = massTotal_/rho;
201 
202  nP = (volumeFraction*volumeTot + delayedVolume_)/(parcels*volumep);
203  break;
204  }
205  case pbNumber:
206  {
207  nP = massTotal_/(rho*volumeTotal_);
208  break;
209  }
210  case pbFixed:
211  {
212  nP = nParticleFixed_;
213  break;
214  }
215  default:
216  {
217  nP = 0.0;
219  (
220  "Foam::scalar "
221  "Foam::InjectionModel<CloudType>::setNumberOfParticles"
222  "("
223  "const label, "
224  "const scalar, "
225  "const scalar, "
226  "const scalar"
227  ")"
228  )<< "Unknown parcelBasis type" << nl
229  << exit(FatalError);
230  }
231  }
232 
233  return nP;
234 }
235 
236 
237 template<class CloudType>
239 (
240  const label parcelsAdded,
241  const scalar massAdded
242 )
243 {
244  const label allParcelsAdded = returnReduce(parcelsAdded, sumOp<label>());
245 
246  if (allParcelsAdded > 0)
247  {
248  Info<< nl
249  << "Cloud: " << this->owner().name()
250  << " injector: " << this->modelName() << nl
251  << " Added " << allParcelsAdded << " new parcels" << nl << endl;
252  }
253 
254  // Increment total number of parcels added
255  parcelsAddedTotal_ += allParcelsAdded;
256 
257  // Increment total mass injected
258  massInjected_ += returnReduce(massAdded, sumOp<scalar>());
259 
260  // Update time for start of next injection
261  time0_ = this->owner().db().time().value();
262 
263  // Increment number of injections
264  nInjections_++;
265 }
266 
267 
268 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
269 
270 template<class CloudType>
272 :
274  SOI_(0.0),
275  volumeTotal_(0.0),
276  massTotal_(0.0),
277  massFlowRate_(owner.db().time(), "massFlowRate"),
278  massInjected_(this->template getModelProperty<scalar>("massInjected")),
279  nInjections_(this->template getModelProperty<label>("nInjections")),
280  parcelsAddedTotal_
281  (
282  this->template getModelProperty<scalar>("parcelsAddedTotal")
283  ),
284  parcelBasis_(pbNumber),
285  nParticleFixed_(0.0),
286  time0_(0.0),
287  timeStep0_(this->template getModelProperty<scalar>("timeStep0")),
288  delayedVolume_(0.0)
289 {}
290 
291 
292 template<class CloudType>
294 (
295  const dictionary& dict,
296  CloudType& owner,
297  const word& modelName,
298  const word& modelType
299 )
300 :
302  SOI_(0.0),
303  volumeTotal_(0.0),
304  massTotal_(0.0),
305  massFlowRate_(owner.db().time(), "massFlowRate"),
306  massInjected_(this->template getModelProperty<scalar>("massInjected")),
307  nInjections_(this->template getModelProperty<scalar>("nInjections")),
309  (
310  this->template getModelProperty<scalar>("parcelsAddedTotal")
311  ),
313  nParticleFixed_(0.0),
314  time0_(owner.db().time().value()),
315  timeStep0_(this->template getModelProperty<scalar>("timeStep0")),
316  delayedVolume_(0.0)
317 {
318  // Provide some info
319  // - also serves to initialise mesh dimensions - needed for parallel runs
320  // due to lazy evaluation of valid mesh dimensions
321  Info<< " Constructing " << owner.mesh().nGeometricD() << "-D injection"
322  << endl;
323 
324  if (owner.solution().transient())
325  {
326  this->coeffDict().lookup("massTotal") >> massTotal_;
327  this->coeffDict().lookup("SOI") >> SOI_;
328  SOI_ = owner.db().time().userTimeToTime(SOI_);
329  }
330  else
331  {
332  massFlowRate_.reset(this->coeffDict());
333  massTotal_ = massFlowRate_.value(owner.db().time().value());
334  }
335 
336  const word parcelBasisType = this->coeffDict().lookup("parcelBasisType");
337 
338  if (parcelBasisType == "mass")
339  {
341  }
342  else if (parcelBasisType == "number")
343  {
345  }
346  else if (parcelBasisType == "fixed")
347  {
349 
350  Info<< " Choosing nParticle to be a fixed value, massTotal "
351  << "variable now does not determine anything."
352  << endl;
353 
354  nParticleFixed_ = readScalar(this->coeffDict().lookup("nParticle"));
355  }
356  else
357  {
359  (
360  "Foam::InjectionModel<CloudType>::InjectionModel"
361  "("
362  "const dictionary&, "
363  "CloudType&, "
364  "const word&"
365  ")"
366  )<< "parcelBasisType must be either 'number', 'mass' or 'fixed'" << nl
367  << exit(FatalError);
368  }
369 }
370 
371 
372 template<class CloudType>
374 (
375  const InjectionModel<CloudType>& im
376 )
377 :
379  SOI_(im.SOI_),
388  time0_(im.time0_),
391 {}
392 
393 
394 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
395 
396 template<class CloudType>
398 {}
399 
400 
401 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
402 
403 template<class CloudType>
405 {
406  // do nothing
407 }
408 
409 
410 template<class CloudType>
412 {
413  label nTotal = 0.0;
414  if (this->owner().solution().transient())
415  {
416  nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
417  }
418  else
419  {
420  nTotal = parcelsToInject(0.0, 1.0);
421  }
422 
423  return massTotal_/nTotal;
424 }
425 
426 
427 template<class CloudType>
428 template<class TrackData>
430 {
431  if (!this->active())
432  {
433  return;
434  }
435 
436  const scalar time = this->owner().db().time().value();
437 
438  // Prepare for next time step
439  label parcelsAdded = 0;
440  scalar massAdded = 0.0;
441  label newParcels = 0;
442  scalar newVolumeFraction = 0.0;
443 
444  if (prepareForNextTimeStep(time, newParcels, newVolumeFraction))
445  {
446  scalar delayedVolume = 0;
447 
448  const scalar trackTime = this->owner().solution().trackTime();
449  const polyMesh& mesh = this->owner().mesh();
450  typename TrackData::cloudType& cloud = td.cloud();
451 
452  // Duration of injection period during this timestep
453  const scalar deltaT =
454  max(0.0, min(trackTime, min(time - SOI_, timeEnd() - time0_)));
455 
456  // Pad injection time if injection starts during this timestep
457  const scalar padTime = max(0.0, SOI_ - time0_);
458 
459  // Introduce new parcels linearly across carrier phase timestep
460  for (label parcelI = 0; parcelI < newParcels; parcelI++)
461  {
462  if (validInjection(parcelI))
463  {
464  // Calculate the pseudo time of injection for parcel 'parcelI'
465  scalar timeInj = time0_ + padTime + deltaT*parcelI/newParcels;
466 
467  // Determine the injection position and owner cell,
468  // tetFace and tetPt
469  label cellI = -1;
470  label tetFaceI = -1;
471  label tetPtI = -1;
472 
474 
476  (
477  parcelI,
478  newParcels,
479  timeInj,
480  pos,
481  cellI,
482  tetFaceI,
483  tetPtI
484  );
485 
486  if (cellI > -1)
487  {
488  // Lagrangian timestep
489  const scalar dt = time - timeInj;
490 
491  // Apply corrections to position for 2-D cases
493 
494  // Create a new parcel
495  parcelType* pPtr =
496  new parcelType(mesh, pos, cellI, tetFaceI, tetPtI);
497 
498  // Check/set new parcel thermo properties
499  cloud.setParcelThermoProperties(*pPtr, dt);
500 
501  // Assign new parcel properties in injection model
502  setProperties(parcelI, newParcels, timeInj, *pPtr);
503 
504  // Check/set new parcel injection properties
505  cloud.checkParcelProperties(*pPtr, dt, fullyDescribed());
506 
507  // Apply correction to velocity for 2-D cases
509  (
510  mesh,
511  mesh.solutionD(),
512  pPtr->U()
513  );
514 
515  // Number of particles per parcel
516  pPtr->nParticle() =
518  (
519  newParcels,
520  newVolumeFraction,
521  pPtr->d(),
522  pPtr->rho()
523  );
524 
525  if (pPtr->nParticle() >= 1.0)
526  {
527  parcelsAdded++;
528  massAdded += pPtr->nParticle()*pPtr->mass();
529 
530  if (pPtr->move(td, dt))
531  {
532  td.cloud().addParticle(pPtr);
533  }
534  else
535  {
536  delete pPtr;
537  }
538  }
539  else
540  {
541  delayedVolume += pPtr->nParticle()*pPtr->volume();
542  delete pPtr;
543  }
544  }
545  }
546  }
547 
548  delayedVolume_ = delayedVolume;
549  }
550 
551  postInjectCheck(parcelsAdded, massAdded);
552 }
553 
554 
555 template<class CloudType>
556 template<class TrackData>
558 (
559  TrackData& td,
560  const scalar trackTime
561 )
562 {
563  if (!this->active())
564  {
565  return;
566  }
567 
568  const polyMesh& mesh = this->owner().mesh();
569  typename TrackData::cloudType& cloud = td.cloud();
570 
572 
573  // Reset counters
574  time0_ = 0.0;
575  label parcelsAdded = 0;
576  scalar massAdded = 0.0;
577 
578  // Set number of new parcels to inject based on first second of injection
579  label newParcels = parcelsToInject(0.0, 1.0);
580 
581  // Inject new parcels
582  for (label parcelI = 0; parcelI < newParcels; parcelI++)
583  {
584  // Volume to inject is split equally amongst all parcel streams
585  scalar newVolumeFraction = 1.0/scalar(newParcels);
586 
587  // Determine the injection position and owner cell,
588  // tetFace and tetPt
589  label cellI = -1;
590  label tetFaceI = -1;
591  label tetPtI = -1;
592 
594 
596  (
597  parcelI,
598  newParcels,
599  0.0,
600  pos,
601  cellI,
602  tetFaceI,
603  tetPtI
604  );
605 
606  if (cellI > -1)
607  {
608  // Apply corrections to position for 2-D cases
610 
611  // Create a new parcel
612  parcelType* pPtr =
613  new parcelType(mesh, pos, cellI, tetFaceI, tetPtI);
614 
615  // Check/set new parcel thermo properties
616  cloud.setParcelThermoProperties(*pPtr, 0.0);
617 
618  // Assign new parcel properties in injection model
619  setProperties(parcelI, newParcels, 0.0, *pPtr);
620 
621  // Check/set new parcel injection properties
622  cloud.checkParcelProperties(*pPtr, 0.0, fullyDescribed());
623 
624  // Apply correction to velocity for 2-D cases
625  meshTools::constrainDirection(mesh, mesh.solutionD(), pPtr->U());
626 
627  // Number of particles per parcel
628  pPtr->nParticle() =
630  (
631  1,
632  newVolumeFraction,
633  pPtr->d(),
634  pPtr->rho()
635  );
636 
637  // Add the new parcel
638  td.cloud().addParticle(pPtr);
639 
640  massAdded += pPtr->nParticle()*pPtr->mass();
641  parcelsAdded++;
642  }
643  }
644 
645  postInjectCheck(parcelsAdded, massAdded);
646 }
647 
648 
649 template<class CloudType>
651 {
652  os << " " << this->modelName() << ":" << nl
653  << " number of parcels added = " << parcelsAddedTotal_ << nl
654  << " mass introduced = " << massInjected_ << nl;
655 
656  if (this->outputTime())
657  {
658  this->setModelProperty("massInjected", massInjected_);
659  this->setModelProperty("nInjections", nInjections_);
660  this->setModelProperty("parcelsAddedTotal", parcelsAddedTotal_);
661  this->setModelProperty("timeStep0", timeStep0_);
662  }
663 }
664 
665 
666 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
667 
668 #include "InjectionModelNew.C"
669 
670 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar pow3(const dimensionedScalar &ds)
virtual void postInjectCheck(const label parcelsAdded, const scalar massAdded)
Post injection checks.
scalar delayedVolume_
Volume that should have been injected, but would lead to.
const CloudType & owner() const
Return const access to the owner cloud.
scalar timeStart() const
Return the start-of-injection time.
virtual bool active() const
Return the model &#39;active&#39; status - default active = true.
Definition: subModelBase.C:154
parcelBasis parcelBasis_
Parcel basis enumeration.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const word & modelType() const
Return const access to the sub-model type.
Definition: subModelBase.C:122
label parcelsAddedTotal_
Running counter of total number of parcels added.
virtual bool fullyDescribed() const =0
Flag to identify whether model fully describes the parcel.
virtual Type value(const scalar x) const
Return value as a function of (scalar) independent variable.
virtual scalar setNumberOfParticles(const label parcels, const scalar volumeFraction, const scalar diameter, const scalar rho)
Set number of particles to inject given parcel properties.
virtual void setPositionAndCell(const label parcelI, const label nParcels, const scalar time, vector &position, label &cellOwner, label &tetFaceI, label &tetPtI)=0
Set the injection position and owner cell, tetFace and tetPt.
void injectSteadyState(TrackData &td, const scalar trackTime)
Main injection loop - steady-state.
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
Definition: polyMesh.C:847
virtual void info(Ostream &os)
Write injection info to stream.
virtual ~InjectionModel()
Destructor.
virtual void updateMesh()
Update mesh.
A class for handling words, derived from string.
Definition: word.H:59
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
virtual bool validInjection(const label parcelI)=0
Additional flag to identify whether or not injection of parcelI is.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
messageStream Info
dynamicFvMesh & mesh
const word & modelName() const
Return const access to the name of the sub-model.
Definition: subModelBase.C:104
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
TimeDataEntry< scalar > massFlowRate_
Mass flow rate profile for steady calculations.
void reset(const dictionary &dict)
Reset entry by re-reading from dictionary.
Definition: TimeDataEntry.C:82
Selector class for relaxation factors, solver type and solution.
Definition: solution.H:48
virtual bool outputTime() const
Flag to indicate when to write a property.
scalar time0_
Continuous phase time at start of injection time step [s].
virtual bool findCellAtPosition(label &cellI, label &tetFaceI, label &tetPtI, vector &position, bool errorOnNotFound=true)
Find the cell that contains the supplied position.
label nInjections_
Number of injections counter.
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
void constrainDirection(const polyMesh &mesh, const Vector< label > &dirs, vector &d)
Set the constrained components of directions/velocity to zero.
Definition: meshTools.C:693
void inject(TrackData &td)
Main injection loop.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
void setModelProperty(const word &entryName, const Type &value)
Add generic property to the sub-model.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
const dictionary & dict() const
Return const access to the cloud dictionary.
Definition: subModelBase.C:110
Templated injection model class.
scalar timeStep0_
Time at start of injection time step [s].
A cloud is a collection of lagrangian particles.
Definition: cloud.H:51
InjectionModel(CloudType &owner)
Construct null from owner.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:404
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
CloudType::parcelType parcelType
Convenience typedef for parcelType.
virtual label parcelsToInject(const scalar time0, const scalar time1)=0
Number of parcels to introduce relative to SOI.
scalar massTotal_
Total mass to inject [kg].
scalar massInjected_
Total mass injected to date [kg].
virtual bool prepareForNextTimeStep(const scalar time, label &newParcels, scalar &newVolumeFraction)
Determine properties for next time step/injection interval.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
error FatalError
Base class for cloud sub-models.
const dictionary & coeffDict() const
Return const access to the coefficients dictionary.
Definition: subModelBase.C:128
static const Vector zero
Definition: Vector.H:80
scalar nParticleFixed_
nParticle to assign to parcels when the &#39;fixed&#39; basis
void constrainToMeshCentre(const polyMesh &mesh, point &pt)
Set the constrained components of position to mesh centre.
Definition: meshTools.C:634
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, parcelType &parcel)=0
Set the parcel properties.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
scalar volumeTotal_
Total volume of particles introduced by this injector [m^3].
virtual scalar timeEnd() const =0
Return the end-of-injection time.
const fvMesh & mesh() const
Return refernce to the mesh.
Definition: DSMCCloudI.H:41
const Time & time() const
Return time.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
dimensionedScalar pos(const dimensionedScalar &ds)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
mathematical constants.
scalar SOI_
Start of injection [s].
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:245
virtual scalar averageParcelMass()
Return the average parcel mass over the injection period.
const Type & value() const
Return const reference to value.