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-2018 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 
342  nParticleFixed_ = readScalar(this->coeffDict().lookup("nParticle"));
343  }
344  else
345  {
347  << "parcelBasisType must be either 'number', 'mass' or 'fixed'"
348  << nl << exit(FatalError);
349  }
350 }
351 
352 
353 template<class CloudType>
355 (
356  const InjectionModel<CloudType>& im
357 )
358 :
360  SOI_(im.SOI_),
369  time0_(im.time0_),
371 {}
372 
373 
374 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
375 
376 template<class CloudType>
378 {}
379 
380 
381 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
382 
383 template<class CloudType>
385 {}
386 
387 
388 template<class CloudType>
390 {
391  label nTotal = 0.0;
392  if (this->owner().solution().transient())
393  {
394  nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
395  }
396  else
397  {
398  nTotal = parcelsToInject(0.0, 1.0);
399  }
400 
401  return massTotal_/nTotal;
402 }
403 
404 
405 template<class CloudType>
406 template<class TrackCloudType>
408 (
409  TrackCloudType& cloud,
410  typename CloudType::parcelType::trackingData& td
411 )
412 {
413  if (!this->active())
414  {
415  return;
416  }
417 
418  const scalar time = this->owner().db().time().value();
419 
420  // Prepare for next time step
421  label parcelsAdded = 0;
422  scalar massAdded = 0.0;
423  label newParcels = 0;
424  scalar newVolumeFraction = 0.0;
425 
426  if (prepareForNextTimeStep(time, newParcels, newVolumeFraction))
427  {
428  const scalar trackTime = this->owner().solution().trackTime();
429  const polyMesh& mesh = this->owner().mesh();
430 
431  // Duration of injection period during this timestep
432  const scalar deltaT =
433  max(0.0, min(trackTime, min(time - SOI_, timeEnd() - time0_)));
434 
435  // Pad injection time if injection starts during this timestep
436  const scalar padTime = max(0.0, SOI_ - time0_);
437 
438  // Introduce new parcels linearly across carrier phase timestep
439  for (label parcelI = 0; parcelI < newParcels; parcelI++)
440  {
441  if (validInjection(parcelI))
442  {
443  // Calculate the pseudo time of injection for parcel 'parcelI'
444  scalar timeInj = time0_ + padTime + deltaT*parcelI/newParcels;
445 
446  // Determine the injection position and owner cell,
447  // tetFace and tetPt
448  label celli = -1;
449  label tetFacei = -1;
450  label tetPti = -1;
451 
452  vector pos = Zero;
453 
455  (
456  parcelI,
457  newParcels,
458  timeInj,
459  pos,
460  celli,
461  tetFacei,
462  tetPti
463  );
464 
465  if (celli > -1)
466  {
467  // Lagrangian timestep
468  const scalar dt = time - timeInj;
469 
470  // Apply corrections to position for 2-D cases
472 
473  // Create a new parcel
474  parcelType* pPtr = new parcelType(mesh, pos, celli);
475 
476  // Check/set new parcel thermo properties
477  cloud.setParcelThermoProperties(*pPtr, dt);
478 
479  // Assign new parcel properties in injection model
480  setProperties(parcelI, newParcels, timeInj, *pPtr);
481 
482  // Check/set new parcel injection properties
483  cloud.checkParcelProperties(*pPtr, dt, fullyDescribed());
484 
485  // Apply correction to velocity for 2-D cases
487  (
488  mesh,
489  mesh.solutionD(),
490  pPtr->U()
491  );
492 
493  // Number of particles per parcel
494  pPtr->nParticle() =
496  (
497  newParcels,
498  newVolumeFraction,
499  pPtr->d(),
500  pPtr->rho()
501  );
502 
503  parcelsAdded ++;
504 
505  massAdded += pPtr->nParticle()*pPtr->mass();
506 
507  if (pPtr->move(cloud, td, dt))
508  {
509  cloud.addParticle(pPtr);
510  }
511  else
512  {
513  delete pPtr;
514  }
515  }
516  }
517  }
518  }
519 
520  postInjectCheck(parcelsAdded, massAdded);
521 }
522 
523 
524 template<class CloudType>
525 template<class TrackCloudType>
527 (
528  TrackCloudType& cloud,
529  typename CloudType::parcelType::trackingData& td,
530  const scalar trackTime
531 )
532 {
533  if (!this->active())
534  {
535  return;
536  }
537 
538  const polyMesh& mesh = this->owner().mesh();
539 
541 
542  // Reset counters
543  time0_ = 0.0;
544  label parcelsAdded = 0;
545  scalar massAdded = 0.0;
546 
547  // Set number of new parcels to inject based on first second of injection
548  label newParcels = parcelsToInject(0.0, 1.0);
549 
550  // Inject new parcels
551  for (label parcelI = 0; parcelI < newParcels; parcelI++)
552  {
553  // Volume to inject is split equally amongst all parcel streams
554  scalar newVolumeFraction = 1.0/scalar(newParcels);
555 
556  // Determine the injection position and owner cell,
557  // tetFace and tetPt
558  label celli = -1;
559  label tetFacei = -1;
560  label tetPti = -1;
561 
562  vector pos = Zero;
563 
565  (
566  parcelI,
567  newParcels,
568  0.0,
569  pos,
570  celli,
571  tetFacei,
572  tetPti
573  );
574 
575  if (celli > -1)
576  {
577  // Apply corrections to position for 2-D cases
579 
580  // Create a new parcel
581  parcelType* pPtr = new parcelType(mesh, pos, celli);
582 
583  // Check/set new parcel thermo properties
584  cloud.setParcelThermoProperties(*pPtr, 0.0);
585 
586  // Assign new parcel properties in injection model
587  setProperties(parcelI, newParcels, 0.0, *pPtr);
588 
589  // Check/set new parcel injection properties
590  cloud.checkParcelProperties(*pPtr, 0.0, fullyDescribed());
591 
592  // Apply correction to velocity for 2-D cases
593  meshTools::constrainDirection(mesh, mesh.solutionD(), pPtr->U());
594 
595  // Number of particles per parcel
596  pPtr->nParticle() =
598  (
599  1,
600  newVolumeFraction,
601  pPtr->d(),
602  pPtr->rho()
603  );
604 
605  // Add the new parcel
606  cloud.addParticle(pPtr);
607 
608  massAdded += pPtr->nParticle()*pPtr->mass();
609  parcelsAdded++;
610  }
611  }
612 
613  postInjectCheck(parcelsAdded, massAdded);
614 }
615 
616 
617 template<class CloudType>
619 {
620  os << " " << this->modelName() << ":" << nl
621  << " number of parcels added = " << parcelsAddedTotal_ << nl
622  << " mass introduced = " << massInjected_ << nl;
623 
624  if (this->writeTime())
625  {
626  this->setModelProperty("massInjected", massInjected_);
627  this->setModelProperty("nInjections", nInjections_);
628  this->setModelProperty("parcelsAddedTotal", parcelsAddedTotal_);
629  this->setModelProperty("timeStep0", timeStep0_);
630  }
631 }
632 
633 
634 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
635 
636 #include "InjectionModelNew.C"
637 
638 // ************************************************************************* //
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:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Templated injection model class.
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
Definition: polyMesh.C:824
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:427
label parcelsAddedTotal_
Running counter of total number of parcels added.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
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
stressControl lookup("compactNormalStress") >> compactNormalStress
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
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
void reset(const dictionary &dict)
Reset entry by re-reading from dictionary.
Definition: TimeFunction1.C:82
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:53
InjectionModel(CloudType &owner)
Construct null from owner.
static const char nl
Definition: Ostream.H:265
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) independent variable.
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:361
virtual bool active() const
Return the model &#39;active&#39; status - default active = true.
Definition: subModelBase.C:154
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
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:576