InjectionModel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 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  this->owner().mesh().findCellFacePt
142  (
143  position,
144  celli,
145  tetFacei,
146  tetPti
147  );
148 
149  if (celli > 0)
150  {
151  proci = Pstream::myProcNo();
152  }
153  }
154 
155  reduce(proci, maxOp<label>());
156 
157  if (proci != Pstream::myProcNo())
158  {
159  celli = -1;
160  tetFacei = -1;
161  tetPti = -1;
162  }
163  }
164 
165  if (proci == -1)
166  {
167  if (errorOnNotFound)
168  {
170  << "Cannot find parcel injection cell. "
171  << "Parcel position = " << p0 << nl
172  << exit(FatalError);
173  }
174  else
175  {
176  return false;
177  }
178  }
179 
180  return true;
181 }
182 
183 
184 template<class CloudType>
186 (
187  const label parcels,
188  const scalar volumeFraction,
189  const scalar diameter,
190  const scalar rho
191 )
192 {
193  scalar nP = 0.0;
194  switch (parcelBasis_)
195  {
196  case pbMass:
197  {
198  scalar volumep = pi/6.0*pow3(diameter);
199  scalar volumeTot = massTotal_/rho;
200 
201  nP = volumeFraction*volumeTot/(parcels*volumep);
202  break;
203  }
204  case pbNumber:
205  {
206  nP = massTotal_/(rho*volumeTotal_);
207  break;
208  }
209  case pbFixed:
210  {
211  nP = nParticleFixed_;
212  break;
213  }
214  default:
215  {
216  nP = 0.0;
218  << "Unknown parcelBasis type" << nl
219  << exit(FatalError);
220  }
221  }
222 
223  return nP;
224 }
225 
226 
227 template<class CloudType>
229 (
230  const label parcelsAdded,
231  const scalar massAdded
232 )
233 {
234  const label allParcelsAdded = returnReduce(parcelsAdded, sumOp<label>());
235 
236  if (allParcelsAdded > 0)
237  {
238  Info<< nl
239  << "Cloud: " << this->owner().name()
240  << " injector: " << this->modelName() << nl
241  << " Added " << allParcelsAdded << " new parcels" << nl << endl;
242  }
243 
244  // Increment total number of parcels added
245  parcelsAddedTotal_ += allParcelsAdded;
246 
247  // Increment total mass injected
248  massInjected_ += returnReduce(massAdded, sumOp<scalar>());
249 
250  // Update time for start of next injection
251  time0_ = this->owner().db().time().value();
252 
253  // Increment number of injections
254  nInjections_++;
255 }
256 
257 
258 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
259 
260 template<class CloudType>
262 :
264  SOI_(0.0),
265  volumeTotal_(0.0),
266  massTotal_(0.0),
267  massFlowRate_(owner.db().time(), "massFlowRate"),
268  massInjected_(this->template getModelProperty<scalar>("massInjected")),
269  nInjections_(this->template getModelProperty<label>("nInjections")),
270  parcelsAddedTotal_
271  (
272  this->template getModelProperty<scalar>("parcelsAddedTotal")
273  ),
274  parcelBasis_(pbNumber),
275  nParticleFixed_(0.0),
276  time0_(0.0),
277  timeStep0_(this->template getModelProperty<scalar>("timeStep0"))
278 {}
279 
280 
281 template<class CloudType>
283 (
284  const dictionary& dict,
285  CloudType& owner,
286  const word& modelName,
287  const word& modelType
288 )
289 :
291  SOI_(0.0),
292  volumeTotal_(0.0),
293  massTotal_(0.0),
294  massFlowRate_(owner.db().time(), "massFlowRate"),
295  massInjected_(this->template getModelProperty<scalar>("massInjected")),
296  nInjections_(this->template getModelProperty<scalar>("nInjections")),
298  (
299  this->template getModelProperty<scalar>("parcelsAddedTotal")
300  ),
302  nParticleFixed_(0.0),
303  time0_(owner.db().time().value()),
304  timeStep0_(this->template getModelProperty<scalar>("timeStep0"))
305 {
306  // Provide some info
307  // - also serves to initialise mesh dimensions - needed for parallel runs
308  // due to lazy evaluation of valid mesh dimensions
309  Info<< " Constructing " << owner.mesh().nGeometricD() << "-D injection"
310  << endl;
311 
312  if (owner.solution().transient())
313  {
314  this->coeffDict().lookup("massTotal") >> massTotal_;
315  this->coeffDict().lookup("SOI") >> SOI_;
316  SOI_ = owner.db().time().userTimeToTime(SOI_);
317  }
318  else
319  {
320  massFlowRate_.reset(this->coeffDict());
321  massTotal_ = massFlowRate_.value(owner.db().time().value());
322  }
323 
324  const word parcelBasisType = this->coeffDict().lookup("parcelBasisType");
325 
326  if (parcelBasisType == "mass")
327  {
329  }
330  else if (parcelBasisType == "number")
331  {
333  }
334  else if (parcelBasisType == "fixed")
335  {
337 
338  Info<< " Choosing nParticle to be a fixed value, massTotal "
339  << "variable now does not determine anything."
340  << endl;
341 
343  this->coeffDict().template lookup<scalar>("nParticle");
344  }
345  else
346  {
348  << "parcelBasisType must be either 'number', 'mass' or 'fixed'"
349  << nl << exit(FatalError);
350  }
351 }
352 
353 
354 template<class CloudType>
356 (
357  const InjectionModel<CloudType>& im
358 )
359 :
361  SOI_(im.SOI_),
370  time0_(im.time0_),
372 {}
373 
374 
375 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
376 
377 template<class CloudType>
379 {}
380 
381 
382 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
383 
384 template<class CloudType>
386 {}
387 
388 
389 template<class CloudType>
391 {
392  label nTotal = 0.0;
393  if (this->owner().solution().transient())
394  {
395  nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
396  }
397  else
398  {
399  nTotal = parcelsToInject(0.0, 1.0);
400  }
401 
402  return massTotal_/nTotal;
403 }
404 
405 
406 template<class CloudType>
407 template<class TrackCloudType>
409 (
410  TrackCloudType& cloud,
411  typename CloudType::parcelType::trackingData& td
412 )
413 {
414  const scalar time = this->owner().db().time().value();
415 
416  // Prepare for next time step
417  label parcelsAdded = 0;
418  scalar massAdded = 0.0;
419  label newParcels = 0;
420  scalar newVolumeFraction = 0.0;
421 
422  if (prepareForNextTimeStep(time, newParcels, newVolumeFraction))
423  {
424  const scalar trackTime = this->owner().solution().trackTime();
425  const polyMesh& mesh = this->owner().mesh();
426 
427  // Duration of injection period during this timestep
428  const scalar deltaT =
429  max(0.0, min(trackTime, min(time - SOI_, timeEnd() - time0_)));
430 
431  // Pad injection time if injection starts during this timestep
432  const scalar padTime = max(0.0, SOI_ - time0_);
433 
434  // Introduce new parcels linearly across carrier phase timestep
435  for (label parcelI = 0; parcelI < newParcels; parcelI++)
436  {
437  if (validInjection(parcelI))
438  {
439  // Calculate the pseudo time of injection for parcel 'parcelI'
440  scalar timeInj = time0_ + padTime + deltaT*parcelI/newParcels;
441 
442  // Determine the injection position and owner cell,
443  // tetFace and tetPt
444  label celli = -1;
445  label tetFacei = -1;
446  label tetPti = -1;
447 
448  vector pos = Zero;
449 
451  (
452  parcelI,
453  newParcels,
454  timeInj,
455  pos,
456  celli,
457  tetFacei,
458  tetPti
459  );
460 
461  if (celli > -1)
462  {
463  // Lagrangian timestep
464  const scalar dt = time - timeInj;
465 
466  // Apply corrections to position for 2-D cases
468 
469  // Create a new parcel
470  parcelType* pPtr = new parcelType(mesh, pos, celli);
471 
472  // Check/set new parcel thermo properties
473  cloud.setParcelThermoProperties(*pPtr, dt);
474 
475  // Assign new parcel properties in injection model
476  setProperties(parcelI, newParcels, timeInj, *pPtr);
477 
478  // Check/set new parcel injection properties
479  cloud.checkParcelProperties(*pPtr, dt, fullyDescribed());
480 
481  // Apply correction to velocity for 2-D cases
483  (
484  mesh,
485  mesh.solutionD(),
486  pPtr->U()
487  );
488 
489  // Number of particles per parcel
490  pPtr->nParticle() =
492  (
493  newParcels,
494  newVolumeFraction,
495  pPtr->d(),
496  pPtr->rho()
497  );
498 
499  parcelsAdded ++;
500 
501  massAdded += pPtr->nParticle()*pPtr->mass();
502 
503  if (pPtr->move(cloud, td, dt))
504  {
505  cloud.addParticle(pPtr);
506  }
507  else
508  {
509  delete pPtr;
510  }
511  }
512  }
513  }
514  }
515 
516  postInjectCheck(parcelsAdded, massAdded);
517 }
518 
519 
520 template<class CloudType>
521 template<class TrackCloudType>
523 (
524  TrackCloudType& cloud,
525  typename CloudType::parcelType::trackingData& td,
526  const scalar trackTime
527 )
528 {
529  const polyMesh& mesh = this->owner().mesh();
530 
532 
533  // Reset counters
534  time0_ = 0.0;
535  label parcelsAdded = 0;
536  scalar massAdded = 0.0;
537 
538  // Set number of new parcels to inject based on first second of injection
539  label newParcels = parcelsToInject(0.0, 1.0);
540 
541  // Inject new parcels
542  for (label parcelI = 0; parcelI < newParcels; parcelI++)
543  {
544  // Volume to inject is split equally amongst all parcel streams
545  scalar newVolumeFraction = 1.0/scalar(newParcels);
546 
547  // Determine the injection position and owner cell,
548  // tetFace and tetPt
549  label celli = -1;
550  label tetFacei = -1;
551  label tetPti = -1;
552 
553  vector pos = Zero;
554 
556  (
557  parcelI,
558  newParcels,
559  0.0,
560  pos,
561  celli,
562  tetFacei,
563  tetPti
564  );
565 
566  if (celli > -1)
567  {
568  // Apply corrections to position for 2-D cases
570 
571  // Create a new parcel
572  parcelType* pPtr = new parcelType(mesh, pos, celli);
573 
574  // Check/set new parcel thermo properties
575  cloud.setParcelThermoProperties(*pPtr, 0.0);
576 
577  // Assign new parcel properties in injection model
578  setProperties(parcelI, newParcels, 0.0, *pPtr);
579 
580  // Check/set new parcel injection properties
581  cloud.checkParcelProperties(*pPtr, 0.0, fullyDescribed());
582 
583  // Apply correction to velocity for 2-D cases
584  meshTools::constrainDirection(mesh, mesh.solutionD(), pPtr->U());
585 
586  // Number of particles per parcel
587  pPtr->nParticle() =
589  (
590  1,
591  newVolumeFraction,
592  pPtr->d(),
593  pPtr->rho()
594  );
595 
596  // Add the new parcel
597  cloud.addParticle(pPtr);
598 
599  massAdded += pPtr->nParticle()*pPtr->mass();
600  parcelsAdded++;
601  }
602  }
603 
604  postInjectCheck(parcelsAdded, massAdded);
605 }
606 
607 
608 template<class CloudType>
610 {
611  os << " " << this->modelName() << ":" << nl
612  << " number of parcels added = " << parcelsAddedTotal_ << nl
613  << " mass introduced = " << massInjected_ << nl;
614 
615  if (this->writeTime())
616  {
617  this->setModelProperty("massInjected", massInjected_);
618  this->setModelProperty("nInjections", nInjections_);
619  this->setModelProperty("parcelsAddedTotal", parcelsAddedTotal_);
620  this->setModelProperty("timeStep0", timeStep0_);
621  }
622 }
623 
624 
625 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
626 
627 #include "InjectionModelNew.C"
628 
629 // ************************************************************************* //
scalar massTotal_
Total mass to inject [kg].
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
void injectSteadyState(TrackCloudType &cloud, typename CloudType::parcelType::trackingData &td, const scalar trackTime)
Main injection loop - steady-state.
scalar time0_
Continuous phase time at start of injection time step [s].
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
virtual void postInjectCheck(const label parcelsAdded, const scalar massAdded)
Post injection checks.
virtual scalar averageParcelMass()
Return the average parcel mass over the injection period.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
Templated injection model class.
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
Definition: polyMesh.C:831
CloudType::parcelType parcelType
Convenience typedef for parcelType.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
label parcelsAddedTotal_
Running counter of total number of parcels added.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
scalar SOI_
Start of injection [s].
Base class for cloud sub-models.
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, parcelType &parcel)=0
Set the parcel properties.
void inject(TrackCloudType &cloud, typename CloudType::parcelType::trackingData &td)
Main injection loop.
virtual scalar timeEnd() const =0
Return the end-of-injection time.
virtual bool writeTime() const
Flag to indicate when to write a property.
const word & modelName() const
Return const access to the name of the sub-model.
Definition: subModelBase.C:104
const dictionary & dict() const
Return const access to the cloud dictionary.
Definition: subModelBase.C:110
virtual ~InjectionModel()
Destructor.
parcelBasis parcelBasis_
Parcel basis enumeration.
const CloudType & owner() const
Return const access to the owner cloud.
dimensionedScalar pos(const dimensionedScalar &ds)
scalar nParticleFixed_
nParticle to assign to parcels when the &#39;fixed&#39; basis
label nInjections_
Number of injections counter.
dynamicFvMesh & mesh
void setModelProperty(const word &entryName, const Type &value)
Add generic property to the sub-model.
mathematical constants.
A class for handling words, derived from string.
Definition: word.H:59
virtual void info(Ostream &os)
Write injection info to stream.
virtual bool validInjection(const label parcelI)=0
Additional flag to identify whether or not injection of parcelI is.
const Type & value() const
Return const reference to value.
static const zero Zero
Definition: zero.H:97
void reset(const dictionary &dict)
Reset entry by re-reading from dictionary.
Definition: TimeFunction1.C:79
const dictionary & coeffDict() const
Return const access to the coefficients dictionary.
Definition: subModelBase.C:128
virtual bool findCellAtPosition(label &celli, label &tetFacei, label &tetPti, vector &position, bool errorOnNotFound=true)
Find the cell that contains the supplied position.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
InjectionModel(CloudType &owner)
Construct null from owner.
static const char nl
Definition: Ostream.H:260
virtual void updateMesh()
Update mesh.
const word & modelType() const
Return const access to the sub-model type.
Definition: subModelBase.C:122
virtual Type value(const scalar x) const
Return value as a function of scalar x.
Definition: TimeFunction1.C:86
const Time & time() const
Return time.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
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.
virtual bool prepareForNextTimeStep(const scalar time, label &newParcels, scalar &newVolumeFraction)
Determine properties for next time step/injection interval.
const fvMesh & mesh() const
Return references to the mesh.
Definition: DSMCCloudI.H:41
dimensionedScalar pow3(const dimensionedScalar &ds)
virtual scalar setNumberOfParticles(const label parcels, const scalar volumeFraction, const scalar diameter, const scalar rho)
Set number of particles to inject given parcel properties.
scalar massInjected_
Total mass injected to date [kg].
scalar timeStep0_
Time at start of injection time step [s].
TimeFunction1< scalar > massFlowRate_
Mass flow rate profile for steady calculations.
messageStream Info
virtual bool fullyDescribed() const =0
Flag to identify whether model fully describes the parcel.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
virtual label parcelsToInject(const scalar time0, const scalar time1)=0
Number of parcels to introduce relative to SOI.
Selector class for relaxation factors, solver type and solution.
Definition: solution.H:48
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
scalar volumeTotal_
Total volume of particles introduced by this injector [m^3].
void constrainDirection(const polyMesh &mesh, const Vector< label > &dirs, vector &d)
Set the constrained components of directions/velocity to zero.
Definition: meshTools.C:671
scalar timeStart() const
Return the start-of-injection time.
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:312
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
void constrainToMeshCentre(const polyMesh &mesh, point &pt)
Set the constrained components of position to mesh centre.
Definition: meshTools.C:613
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844