All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ConeInjection.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-2024 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 "ConeInjection.H"
27 #include "Constant.H"
28 #include "mathematicalConstants.H"
29 
30 using namespace Foam::constant::mathematical;
31 
32 // * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
33 
34 template<class CloudType>
36 {
37  const word injectionMethod =
38  this->coeffDict().template lookupOrDefault<word>
39  (
40  "injectionMethod",
42  );
43 
44  if (injectionMethod == "point" || injectionMethod == word::null)
45  {
46  injectionMethod_ = imPoint;
47 
48  topoChange();
49  }
50  else if (injectionMethod == "disc")
51  {
52  injectionMethod_ = imDisc;
53 
54  dInner_ =
55  this->coeffDict().template lookup<scalar>("dInner", dimLength);
56  dOuter_ =
57  this->coeffDict().template lookup<scalar>("dOuter", dimLength);
58  }
59  else
60  {
62  << "injectionMethod must be either 'point' or 'disc'"
63  << exit(FatalError);
64  }
65 }
66 
67 
68 template<class CloudType>
70 {
71  const word flowType =
72  this->coeffDict().template lookupOrDefault<word>
73  (
74  "flowType",
76  );
77 
78  if (flowType == "constantVelocity" || flowType == word::null)
79  {
80  flowType_ = ftConstantVelocity;
81 
82  Umag_.reset
83  (
85  (
86  "Umag",
87  this->owner().db().time().userUnits(),
89  this->coeffDict()
90  ).ptr()
91  );
92  }
93  else if (flowType == "pressureDrivenVelocity")
94  {
95  flowType_ = ftPressureDrivenVelocity;
96 
97  Pinj_.reset
98  (
100  (
101  "Pinj",
102  this->owner().db().time().userUnits(),
103  dimPressure,
104  this->coeffDict()
105  ).ptr()
106  );
107  }
108  else if (flowType == "flowRateAndDischarge")
109  {
110  flowType_ = ftFlowRateAndDischarge;
111 
112  dInner_ =
113  this->coeffDict().template lookup<scalar>("dInner", dimLength);
114  dOuter_ =
115  this->coeffDict().template lookup<scalar>("dOuter", dimLength);
116 
117  Cd_.reset
118  (
120  (
121  "Cd",
122  this->owner().db().time().userUnits(),
123  dimless,
124  this->coeffDict()
125  ).ptr()
126  );
127  }
128  else
129  {
131  << "flowType must be either 'constantVelocity', "
132  << "'pressureDrivenVelocity' or 'flowRateAndDischarge'"
133  << exit(FatalError);
134  }
135 }
136 
137 
138 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
139 
140 template<class CloudType>
142 (
143  const dictionary& dict,
144  CloudType& owner,
145  const word& modelName
146 )
147 :
148  InjectionModel<CloudType>(dict, owner, modelName, typeName),
149  injectionMethod_(imPoint),
150  flowType_(ftConstantVelocity),
151  position_
152  (
154  (
155  "position",
156  this->owner().db().time().userUnits(),
157  dimLength,
158  this->coeffDict()
159  )
160  ),
161  direction_
162  (
164  (
165  "direction",
166  this->owner().db().time().userUnits(),
167  dimless,
168  this->coeffDict()
169  )
170  ),
171  injectorCoordinates_(barycentric::uniform(NaN)),
172  injectorCell_(-1),
173  injectorTetFace_(-1),
174  injectorTetPt_(-1),
175  duration_(this->readDuration(dict, owner)),
176  massFlowRate_(this->readMassFlowRate(dict, owner, duration_)),
177  parcelsPerSecond_(this->readParcelsPerSecond(dict, owner)),
178  thetaInner_
179  (
180  Function1<scalar>::New
181  (
182  "thetaInner",
183  this->owner().db().time().userUnits(),
184  unitDegrees,
185  this->coeffDict()
186  )
187  ),
188  thetaOuter_
189  (
190  Function1<scalar>::New
191  (
192  "thetaOuter",
193  this->owner().db().time().userUnits(),
194  unitDegrees,
195  this->coeffDict()
196  )
197  ),
198  sizeDistribution_
199  (
201  (
202  dimLength,
203  this->coeffDict().subDict("sizeDistribution"),
204  this->sizeSampleQ(),
205  owner.rndGen().generator()
206  )
207  ),
208  dInner_(vGreat),
209  dOuter_(vGreat),
210  Umag_(nullptr),
211  Cd_(nullptr),
212  Pinj_(nullptr)
213 {
214  setInjectionMethod();
215 
216  setFlowType();
217 
218  topoChange();
219 }
220 
221 
222 template<class CloudType>
224 (
225  const ConeInjection<CloudType>& im
226 )
227 :
229  injectionMethod_(im.injectionMethod_),
230  flowType_(im.flowType_),
231  position_(im.position_, false),
232  direction_(im.direction_, false),
233  injectorCoordinates_(im.injectorCoordinates_),
234  injectorCell_(im.injectorCell_),
235  injectorTetFace_(im.injectorTetFace_),
236  injectorTetPt_(im.injectorTetPt_),
237  duration_(im.duration_),
238  massFlowRate_(im.massFlowRate_, false),
239  parcelsPerSecond_(im.parcelsPerSecond_, false),
240  thetaInner_(im.thetaInner_, false),
241  thetaOuter_(im.thetaOuter_, false),
242  sizeDistribution_(im.sizeDistribution_, false),
243  dInner_(im.dInner_),
244  dOuter_(im.dOuter_),
245  Umag_(im.Umag_, false),
246  Cd_(im.Cd_, false),
247  Pinj_(im.Pinj_, false)
248 {}
249 
250 
251 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
252 
253 template<class CloudType>
255 {}
256 
257 
258 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
259 
260 template<class CloudType>
262 {
263  if (injectionMethod_ == imPoint && position_->constant())
264  {
265  vector position = position_->value(0);
266  this->findCellAtPosition
267  (
268  position,
269  injectorCoordinates_,
270  injectorCell_,
271  injectorTetFace_,
272  injectorTetPt_
273  );
274  }
275 }
276 
277 
278 template<class CloudType>
280 {
281  return this->SOI_ + duration_;
282 }
283 
284 
285 template<class CloudType>
287 (
288  const scalar time0,
289  const scalar time1
290 )
291 {
292  if (time0 >= 0 && time0 < duration_)
293  {
295  //return floor(parcelsPerSecond_->integral(time0, time1));
296 
297  // Modified calculation to make numbers exact
298  return
299  floor
300  (
301  parcelsPerSecond_->integral(0, time1)
302  - this->parcelsAddedTotal()
303  );
304  }
305  else
306  {
307  return 0;
308  }
309 }
310 
311 
312 template<class CloudType>
314 (
315  const scalar time0,
316  const scalar time1
317 )
318 {
319  if (time0 >= 0 && time0 < duration_)
320  {
321  return massFlowRate_->integral(time0, time1);
322  }
323  else
324  {
325  return 0;
326  }
327 }
328 
329 
330 template<class CloudType>
332 (
333  const label parcelI,
334  const label,
335  const scalar time,
336  barycentric& coordinates,
337  label& celli,
338  label& tetFacei,
339  label& tetPti,
340  label& facei
341 )
342 {
343  randomGenerator& rndGen = this->owner().rndGen();
344 
345  const scalar t = time - this->SOI_;
346 
347  switch (injectionMethod_)
348  {
349  case imPoint:
350  {
351  const point pos = position_->value(t);
352  if (position_->constant())
353  {
354  coordinates = injectorCoordinates_;
355  celli = injectorCell_;
356  tetFacei = injectorTetFace_;
357  tetPti = injectorTetPt_;
358  }
359  else
360  {
361  this->findCellAtPosition
362  (
363  pos,
364  coordinates,
365  celli,
366  tetFacei,
367  tetPti,
368  false
369  );
370  }
371  break;
372  }
373  case imDisc:
374  {
375  const scalar beta = twoPi*this->globalScalar01(rndGen);
376  const scalar frac = this->globalScalar01(rndGen);
377  const vector n = normalised(direction_->value(t));
378  const vector t1 = normalised(perpendicular(n));
379  const vector t2 = normalised(n ^ t1);
380  const vector tanVec = t1*cos(beta) + t2*sin(beta);
381  const scalar d = sqrt((1 - frac)*sqr(dInner_) + frac*sqr(dOuter_));
382  const point pos = position_->value(t) + d/2*tanVec;
383  this->findCellAtPosition
384  (
385  pos,
386  coordinates,
387  celli,
388  tetFacei,
389  tetPti,
390  false
391  );
392  break;
393  }
394  default:
395  {
396  break;
397  }
398  }
399 }
400 
401 
402 template<class CloudType>
404 (
405  const label parcelI,
406  const label,
407  const scalar time,
408  typename CloudType::parcelType::trackingData& td,
409  typename CloudType::parcelType& parcel
410 )
411 {
412  const polyMesh& mesh = this->owner().mesh();
413 
414  randomGenerator& rndGen = this->owner().rndGen();
415 
416  const scalar t = time - this->SOI_;
417 
418  // Get the angle from the axis and the vector perpendicular from the axis.
419  // If injecting at a point, then these are calculated from two new random
420  // numbers. If a disc, then these calculations have already been done in
421  // setPositionAndCell, so the angle and vector can be reverse engineered
422  // from the position.
423  scalar theta = vGreat;
424  vector tanVec = vector::max;
425  switch (injectionMethod_)
426  {
427  case imPoint:
428  {
429  const scalar beta = twoPi*rndGen.scalar01();
430  const scalar frac = rndGen.scalar01();
431  const vector n = normalised(direction_->value(t));
432  const vector t1 = normalised(perpendicular(n));
433  const vector t2 = normalised(n ^ t1);
434  tanVec = t1*cos(beta) + t2*sin(beta);
435  theta =
436  sqrt
437  (
438  (1 - frac)*sqr(thetaInner_->value(t))
439  + frac*sqr(thetaOuter_->value(t))
440  );
441  break;
442  }
443  case imDisc:
444  {
445  const scalar r = mag(parcel.position(mesh) - position_->value(t));
446  const scalar frac = (2*r - dInner_)/(dOuter_ - dInner_);
447  tanVec = normalised(parcel.position(mesh) - position_->value(t));
448  theta =
449  (1 - frac)*thetaInner_->value(t)
450  + frac*thetaOuter_->value(t);
451  break;
452  }
453  default:
454  {
455  break;
456  }
457  }
458 
459  // The direction of injection
460  const vector dirVec =
461  normalised
462  (
463  cos(theta)*normalised(direction_->value(t))
464  + sin(theta)*tanVec
465  );
466 
467  // Set the velocity
468  switch (flowType_)
469  {
470  case ftConstantVelocity:
471  {
472  parcel.U() = Umag_->value(t)*dirVec;
473  break;
474  }
475  case ftPressureDrivenVelocity:
476  {
477  const scalar pAmbient = this->owner().pAmbient();
478  const scalar rho = parcel.rho();
479  const scalar Umag = ::sqrt(2*(Pinj_->value(t) - pAmbient)/rho);
480  parcel.U() = Umag*dirVec;
481  break;
482  }
483  case ftFlowRateAndDischarge:
484  {
485  const scalar A = 0.25*pi*(sqr(dOuter_) - sqr(dInner_));
486  const scalar Umag =
487  massFlowRate_->value(t)/(parcel.rho()*Cd_->value(t)*A);
488  parcel.U() = Umag*dirVec;
489  break;
490  }
491  default:
492  {
493  break;
494  }
495  }
496 
497  // Set the particle diameter
498  parcel.d() = sizeDistribution_->sample();
499 }
500 
501 
502 template<class CloudType>
504 {
505  return false;
506 }
507 
508 
509 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
label n
This injector injects particles in a number of cones. The user specifies a position and a direction t...
virtual ~ConeInjection()
Destructor.
virtual void topoChange()
Set injector locations when mesh is updated.
ConeInjection(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, typename CloudType::parcelType::trackingData &td, typename CloudType::parcelType &parcel)
Set the parcel properties.
virtual void setPositionAndCell(const label parcelI, const label nParcels, const scalar time, barycentric &coordinates, label &celli, label &tetFacei, label &tetPti, label &facei)
Set the injection position and owner cell, tetFace and tetPt.
virtual label nParcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
scalar timeEnd() const
Return the end-of-injection time.
virtual scalar massToInject(const scalar time0, const scalar time1)
Parcel mass to introduce relative to SOI.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:80
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:225
Run-time selectable general function of one variable.
Definition: Function1.H:125
static autoPtr< Function1< scalar > > New(const word &name, const Function1s::unitConversions &units, const dictionary &dict)
Select from dictionary.
Definition: Function1New.C:32
Templated injection model class.
static const Form max
Definition: VectorSpace.H:120
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
const Type & value() const
Return const reference to value.
Base class for statistical distributions.
Definition: distribution.H:76
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
Random number generator.
scalar scalar01()
Return a scalar uniformly distributed between zero and one.
A class for handling words, derived from string.
Definition: word.H:62
static const word null
An empty word.
Definition: word.H:77
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
mathematical constants.
const scalar twoPi(2 *pi)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
dimensionedScalar pos(const dimensionedScalar &ds)
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
const dimensionSet dimPressure
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimless
dimensionedScalar sin(const dimensionedScalar &ds)
const dimensionSet dimLength
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensionSet normalised(const dimensionSet &)
Definition: dimensionSet.C:510
dimensioned< scalar > mag(const dimensioned< Type > &)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
const dimensionSet dimVelocity
error FatalError
dimensionSet perpendicular(const dimensionSet &)
Definition: dimensionSet.C:516
dimensionedScalar cos(const dimensionedScalar &ds)
const unitConversion unitDegrees
dictionary dict
randomGenerator rndGen(653213)