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-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 "ConeInjection.H"
27 #include "TimeFunction1.H"
28 #include "Constant.H"
29 #include "mathematicalConstants.H"
30 #include "unitConversion.H"
31 
32 using namespace Foam::constant::mathematical;
33 
34 // * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
35 
36 template<class CloudType>
38 {
39  const word injectionMethod =
40  this->coeffDict().template lookupOrDefault<word>
41  (
42  "injectionMethod",
44  );
45 
46  if (injectionMethod == "point" || injectionMethod == word::null)
47  {
48  injectionMethod_ = imPoint;
49 
50  updateMesh();
51  }
52  else if (injectionMethod == "disc")
53  {
54  injectionMethod_ = imDisc;
55 
56  this->coeffDict().lookup("dInner") >> dInner_;
57  this->coeffDict().lookup("dOuter") >> dOuter_;
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(this->coeffDict());
83  }
84  else if (flowType == "pressureDrivenVelocity")
85  {
86  flowType_ = ftPressureDrivenVelocity;
87 
88  Pinj_.reset(this->coeffDict());
89  }
90  else if (flowType == "flowRateAndDischarge")
91  {
92  flowType_ = ftFlowRateAndDischarge;
93 
94  this->coeffDict().lookup("dInner") >> dInner_;
95  this->coeffDict().lookup("dOuter") >> dOuter_;
96 
97  Cd_.reset(this->coeffDict());
98  }
99  else
100  {
102  << "flowType must be either 'constantVelocity', "
103  << "'pressureDrivenVelocity' or 'flowRateAndDischarge'"
104  << exit(FatalError);
105  }
106 }
107 
108 
109 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
110 
111 template<class CloudType>
113 (
114  const dictionary& dict,
115  CloudType& owner,
116  const word& modelName
117 )
118 :
119  InjectionModel<CloudType>(dict, owner, modelName, typeName),
120  injectionMethod_(imPoint),
121  flowType_(ftConstantVelocity),
122  position_
123  (
125  (
126  owner.db().time(),
127  "position",
128  this->coeffDict()
129  )
130  ),
131  positionIsConstant_(isA<Function1s::Constant<vector>>(position_)),
132  direction_
133  (
135  (
136  owner.db().time(),
137  "direction",
138  this->coeffDict()
139  )
140  ),
141  injectorCell_(-1),
142  injectorTetFace_(-1),
143  injectorTetPt_(-1),
144  duration_(this->coeffDict().template lookup<scalar>("duration")),
145  parcelsPerSecond_
146  (
147  this->coeffDict().template lookup<scalar>("parcelsPerSecond")
148  ),
149  flowRateProfile_
150  (
152  (
153  owner.db().time(),
154  "flowRateProfile",
155  this->coeffDict()
156  )
157  ),
158  thetaInner_
159  (
161  (
162  owner.db().time(),
163  "thetaInner",
164  this->coeffDict()
165  )
166  ),
167  thetaOuter_
168  (
170  (
171  owner.db().time(),
172  "thetaOuter",
173  this->coeffDict()
174  )
175  ),
176  sizeDistribution_
177  (
179  (
180  this->coeffDict().subDict("sizeDistribution"), owner.rndGen()
181  )
182  ),
183  dInner_(vGreat),
184  dOuter_(vGreat),
185  Umag_(owner.db().time(), "Umag"),
186  Cd_(owner.db().time(), "Cd"),
187  Pinj_(owner.db().time(), "Pinj")
188 {
189  duration_ = owner.db().time().userTimeToTime(duration_);
190 
191  setInjectionMethod();
192 
193  setFlowType();
194 
195  // Set total volume to inject
196  this->volumeTotal_ = flowRateProfile_.integral(0, duration_);
197 
198  updateMesh();
199 }
200 
201 
202 template<class CloudType>
204 (
205  const ConeInjection<CloudType>& im
206 )
207 :
209  injectionMethod_(im.injectionMethod_),
210  flowType_(im.flowType_),
211  position_(im.position_),
212  positionIsConstant_(im.positionIsConstant_),
213  direction_(im.direction_),
214  injectorCell_(im.injectorCell_),
215  injectorTetFace_(im.injectorTetFace_),
216  injectorTetPt_(im.injectorTetPt_),
217  duration_(im.duration_),
218  parcelsPerSecond_(im.parcelsPerSecond_),
219  flowRateProfile_(im.flowRateProfile_),
220  thetaInner_(im.thetaInner_),
221  thetaOuter_(im.thetaOuter_),
222  sizeDistribution_(im.sizeDistribution_().clone().ptr()),
223  dInner_(im.dInner_),
224  dOuter_(im.dOuter_),
225  Umag_(im.Umag_),
226  Cd_(im.Cd_),
227  Pinj_(im.Pinj_)
228 {}
229 
230 
231 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
232 
233 template<class CloudType>
235 {}
236 
237 
238 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
239 
240 template<class CloudType>
242 {
243  if (injectionMethod_ == imPoint && positionIsConstant_)
244  {
245  vector position = position_.value(0);
246  this->findCellAtPosition
247  (
248  injectorCell_,
249  injectorTetFace_,
250  injectorTetPt_,
251  position
252  );
253  }
254 }
255 
256 
257 template<class CloudType>
259 {
260  return this->SOI_ + duration_;
261 }
262 
263 
264 template<class CloudType>
266 (
267  const scalar time0,
268  const scalar time1
269 )
270 {
271  if (time0 >= 0 && time0 < duration_)
272  {
274  //return floor(parcelsPerSecond_*(time1 - time0));
275 
276  // Modified calculation to make numbers exact
277  return floor(parcelsPerSecond_*time1 - this->parcelsAddedTotal());
278  }
279  else
280  {
281  return 0;
282  }
283 }
284 
285 
286 template<class CloudType>
288 (
289  const scalar time0,
290  const scalar time1
291 )
292 {
293  if (time0 >= 0 && time0 < duration_)
294  {
295  return flowRateProfile_.integral(time0, time1);
296  }
297  else
298  {
299  return 0;
300  }
301 }
302 
303 
304 template<class CloudType>
306 (
307  const label parcelI,
308  const label,
309  const scalar time,
310  vector& position,
311  label& cellOwner,
312  label& tetFacei,
313  label& tetPti
314 )
315 {
316  Random& rndGen = this->owner().rndGen();
317 
318  const scalar t = time - this->SOI_;
319 
320  switch (injectionMethod_)
321  {
322  case imPoint:
323  {
324  position = position_.value(t);
325  if (positionIsConstant_)
326  {
327  cellOwner = injectorCell_;
328  tetFacei = injectorTetFace_;
329  tetPti = injectorTetPt_;
330  }
331  else
332  {
333  this->findCellAtPosition
334  (
335  cellOwner,
336  tetFacei,
337  tetPti,
338  position,
339  false
340  );
341  }
342  break;
343  }
344  case imDisc:
345  {
346  const scalar beta = twoPi*rndGen.globalScalar01();
347  const scalar frac = rndGen.globalScalar01();
348  const vector n = normalised(direction_.value(t));
349  const vector t1 = normalised(perpendicular(n));
350  const vector t2 = normalised(n ^ t1);
351  const vector tanVec = t1*cos(beta) + t2*sin(beta);
352  const scalar d = sqrt((1 - frac)*sqr(dInner_) + frac*sqr(dOuter_));
353  position = position_.value(t) + d/2*tanVec;
354  this->findCellAtPosition
355  (
356  cellOwner,
357  tetFacei,
358  tetPti,
359  position,
360  false
361  );
362  break;
363  }
364  default:
365  {
366  break;
367  }
368  }
369 }
370 
371 
372 template<class CloudType>
374 (
375  const label parcelI,
376  const label,
377  const scalar time,
378  typename CloudType::parcelType& parcel
379 )
380 {
381  Random& rndGen = this->owner().rndGen();
382 
383  const scalar t = time - this->SOI_;
384 
385  // Get the angle from the axis and the vector perpendicular from the axis.
386  // If injecting at a point, then these are calculated from two new random
387  // numbers. If a disc, then these calculations have already been done in
388  // setPositionAndCell, so the angle and vector can be reverse engineered
389  // from the position.
390  scalar theta = vGreat;
391  vector tanVec = vector::max;
392  switch (injectionMethod_)
393  {
394  case imPoint:
395  {
396  const scalar beta = twoPi*rndGen.scalar01();
397  const scalar frac = rndGen.scalar01();
398  const vector n = normalised(direction_.value(t));
399  const vector t1 = normalised(perpendicular(n));
400  const vector t2 = normalised(n ^ t1);
401  tanVec = t1*cos(beta) + t2*sin(beta);
402  theta =
403  degToRad
404  (
405  sqrt
406  (
407  (1 - frac)*sqr(thetaInner_.value(t))
408  + frac*sqr(thetaOuter_.value(t))
409  )
410  );
411  break;
412  }
413  case imDisc:
414  {
415  const scalar r = mag(parcel.position() - position_.value(t));
416  const scalar frac = (2*r - dInner_)/(dOuter_ - dInner_);
417  tanVec = normalised(parcel.position() - position_.value(t));
418  theta =
419  degToRad
420  (
421  (1 - frac)*thetaInner_.value(t)
422  + frac*thetaOuter_.value(t)
423  );
424  break;
425  }
426  default:
427  {
428  break;
429  }
430  }
431 
432  // The direction of injection
433  const vector dirVec =
434  normalised
435  (
436  cos(theta)*normalised(direction_.value(t))
437  + sin(theta)*tanVec
438  );
439 
440  // Set the velocity
441  switch (flowType_)
442  {
443  case ftConstantVelocity:
444  {
445  parcel.U() = Umag_.value(t)*dirVec;
446  break;
447  }
448  case ftPressureDrivenVelocity:
449  {
450  const scalar pAmbient = this->owner().pAmbient();
451  const scalar rho = parcel.rho();
452  const scalar Umag = ::sqrt(2*(Pinj_.value(t) - pAmbient)/rho);
453  parcel.U() = Umag*dirVec;
454  break;
455  }
456  case ftFlowRateAndDischarge:
457  {
458  const scalar A = 0.25*pi*(sqr(dOuter_) - sqr(dInner_));
459  const scalar massFlowRate =
460  this->massTotal()*flowRateProfile_.value(t)/this->volumeTotal();
461  const scalar Umag =
462  massFlowRate/(parcel.rho()*Cd_.value(t)*A);
463  parcel.U() = Umag*dirVec;
464  break;
465  }
466  default:
467  {
468  break;
469  }
470  }
471 
472  // Set the particle diameter
473  parcel.d() = sizeDistribution_->sample();
474 }
475 
476 
477 template<class CloudType>
479 {
480  return false;
481 }
482 
483 
484 template<class CloudType>
486 {
487  return true;
488 }
489 
490 
491 // ************************************************************************* //
bool isA(const Type &t)
Check if a dynamic_cast to typeid is possible.
Definition: typeInfo.H:134
ConeInjection(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
dictionary dict
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
Templated injection model class.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
virtual bool validInjection(const label parcelI)
Return flag to identify whether or not injection of parcelI is.
virtual label parcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
dimensionedScalar sqrt(const dimensionedScalar &ds)
virtual void setPositionAndCell(const label parcelI, const label nParcels, const scalar time, vector &position, label &cellOwner, label &tetFacei, label &tetPti)
Set the injection position and owner cell, tetFace and tetPt.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
Random rndGen(label(0))
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
dimensionedScalar cos(const dimensionedScalar &ds)
mathematical constants.
A class for handling words, derived from string.
Definition: word.H:59
Light wrapper around Function1 to provide a mechanism to update time-based entries.
scalar timeEnd() const
Return the end-of-injection time.
static const word null
An empty word.
Definition: word.H:77
virtual scalar volumeToInject(const scalar time0, const scalar time1)
Volume of parcels to introduce relative to SOI.
scalar globalScalar01()
Advance the state and return a scalar sample from a uniform.
Definition: Random.C:94
Random number generator.
Definition: Random.H:57
const scalar twoPi(2 *pi)
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:221
dimensionedScalar sin(const dimensionedScalar &ds)
virtual ~ConeInjection()
Destructor.
Vector< Cmpt > perpendicular(const Vector< Cmpt > &v)
Definition: VectorI.H:166
static autoPtr< distributionModel > New(const dictionary &dict, Random &rndGen)
Selector.
Templated function that returns a constant value.
Definition: Constant.H:56
virtual void updateMesh()
Set injector locations when mesh is updated.
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, typename CloudType::parcelType &parcel)
Set the parcel properties.
scalar scalar01()
Advance the state and return a scalar sample from a uniform.
Definition: RandomI.H:57
This injector injects particles in a number of cones. The user specifies a position and a direction t...