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 t0,
289  const scalar t1
290 )
291 {
292  if (t1 >= 0 && t0 < duration_)
293  {
294  return parcelsPerSecond_->integral(max(t0, 0), min(t1, duration_));
295  }
296  else
297  {
298  return 0;
299  }
300 }
301 
302 
303 template<class CloudType>
305 (
306  const scalar t0,
307  const scalar t1
308 )
309 {
310  if (t1 >= 0 && t0 < duration_)
311  {
312  return massFlowRate_->integral(max(t0, 0), min(t1, duration_));
313  }
314  else
315  {
316  return 0;
317  }
318 }
319 
320 
321 template<class CloudType>
323 (
324  const label parcelI,
325  const label,
326  const scalar time,
328  label& celli,
329  label& tetFacei,
330  label& tetPti,
331  label& facei
332 )
333 {
334  randomGenerator& rndGen = this->owner().rndGen();
335 
336  const scalar t = time - this->SOI_;
337 
338  switch (injectionMethod_)
339  {
340  case imPoint:
341  {
342  const point pos = position_->value(t);
343  if (position_->constant())
344  {
345  coordinates = injectorCoordinates_;
346  celli = injectorCell_;
347  tetFacei = injectorTetFace_;
348  tetPti = injectorTetPt_;
349  }
350  else
351  {
352  this->findCellAtPosition
353  (
354  pos,
355  coordinates,
356  celli,
357  tetFacei,
358  tetPti,
359  false
360  );
361  }
362  break;
363  }
364  case imDisc:
365  {
366  const scalar beta = twoPi*this->globalScalar01(rndGen);
367  const scalar frac = this->globalScalar01(rndGen);
368  const vector n = normalised(direction_->value(t));
369  const vector t1 = normalised(perpendicular(n));
370  const vector t2 = normalised(n ^ t1);
371  const vector tanVec = t1*cos(beta) + t2*sin(beta);
372  const scalar d = sqrt((1 - frac)*sqr(dInner_) + frac*sqr(dOuter_));
373  const point pos = position_->value(t) + d/2*tanVec;
374  this->findCellAtPosition
375  (
376  pos,
377  coordinates,
378  celli,
379  tetFacei,
380  tetPti,
381  false
382  );
383  break;
384  }
385  default:
386  {
387  break;
388  }
389  }
390 }
391 
392 
393 template<class CloudType>
395 (
396  const label parcelI,
397  const label,
398  const scalar time,
399  typename CloudType::parcelType::trackingData& td,
400  typename CloudType::parcelType& parcel
401 )
402 {
403  const polyMesh& mesh = this->owner().mesh();
404 
405  randomGenerator& rndGen = this->owner().rndGen();
406 
407  const scalar t = time - this->SOI_;
408 
409  // Get the angle from the axis and the vector perpendicular from the axis.
410  // If injecting at a point, then these are calculated from two new random
411  // numbers. If a disc, then these calculations have already been done in
412  // setPositionAndCell, so the angle and vector can be reverse engineered
413  // from the position.
414  scalar theta = vGreat;
415  vector tanVec = vector::max;
416  switch (injectionMethod_)
417  {
418  case imPoint:
419  {
420  const scalar beta = twoPi*rndGen.scalar01();
421  const scalar frac = rndGen.scalar01();
422  const vector n = normalised(direction_->value(t));
423  const vector t1 = normalised(perpendicular(n));
424  const vector t2 = normalised(n ^ t1);
425  tanVec = t1*cos(beta) + t2*sin(beta);
426  theta =
427  sqrt
428  (
429  (1 - frac)*sqr(thetaInner_->value(t))
430  + frac*sqr(thetaOuter_->value(t))
431  );
432  break;
433  }
434  case imDisc:
435  {
436  const scalar r = mag(parcel.position(mesh) - position_->value(t));
437  const scalar frac = (2*r - dInner_)/(dOuter_ - dInner_);
438  tanVec = normalised(parcel.position(mesh) - position_->value(t));
439  theta =
440  (1 - frac)*thetaInner_->value(t)
441  + frac*thetaOuter_->value(t);
442  break;
443  }
444  default:
445  {
446  break;
447  }
448  }
449 
450  // The direction of injection
451  const vector dirVec =
452  normalised
453  (
454  cos(theta)*normalised(direction_->value(t))
455  + sin(theta)*tanVec
456  );
457 
458  // Set the velocity
459  switch (flowType_)
460  {
461  case ftConstantVelocity:
462  {
463  parcel.U() = Umag_->value(t)*dirVec;
464  break;
465  }
466  case ftPressureDrivenVelocity:
467  {
468  const scalar pAmbient = this->owner().pAmbient();
469  const scalar rho = parcel.rho();
470  const scalar Umag = ::sqrt(2*(Pinj_->value(t) - pAmbient)/rho);
471  parcel.U() = Umag*dirVec;
472  break;
473  }
474  case ftFlowRateAndDischarge:
475  {
476  const scalar A = 0.25*pi*(sqr(dOuter_) - sqr(dInner_));
477  const scalar Umag =
478  massFlowRate_->value(t)/(parcel.rho()*Cd_->value(t)*A);
479  parcel.U() = Umag*dirVec;
480  break;
481  }
482  default:
483  {
484  break;
485  }
486  }
487 
488  // Set the particle diameter
489  parcel.d() = sizeDistribution_->sample();
490 }
491 
492 
493 template<class CloudType>
495 {
496  return false;
497 }
498 
499 
500 // ************************************************************************* //
label n
This injector injects particles in a number of cones. The user specifies a position and a direction t...
virtual scalar nParcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
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 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 keywords followed by any number of values (e.g. words and numbers) or sub-dictionaries.
Definition: dictionary.H:162
const Type & value() const
Return const reference to value.
Base class for statistical distributions.
Definition: distribution.H:76
const polyMesh & mesh() const
Return reference to polyMesh.
Definition: fvMesh.H:443
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
Foam::fvMesh mesh(Foam::IOobject(regionName, runTime.name(), runTime, Foam::IOobject::MUST_READ), false)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
mathematical constants.
const scalar twoPi(2 *pi)
static const coefficient A("A", dimPressure, 611.21)
barycentric coordinates(const polyMesh &mesh, const point &position, const label celli, const label facei, const label faceTrii, const scalar stepFraction)
Return the coordinates given the position and tet topology.
Definition: tracking.C:1259
point position(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label facei, const label faceTrii, const scalar stepFraction)
Return the position given the coordinates and tet topology.
Definition: trackingI.H:224
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
const dimensionSet dimless
dimensionedScalar sin(const dimensionedScalar &ds)
const dimensionSet dimLength
void mag(LagrangianPatchField< scalar > &f, const LagrangianPatchField< Type > &f1)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionSet normalised(const dimensionSet &)
Definition: dimensionSet.C:507
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
void sqr(LagrangianPatchField< typename outerProduct< Type, Type >::type > &f, const LagrangianPatchField< Type > &f1)
const dimensionSet dimVelocity
error FatalError
dimensionSet perpendicular(const dimensionSet &)
Definition: dimensionSet.C:513
void sqrt(LagrangianPatchField< scalar > &f, const LagrangianPatchField< scalar > &f1)
tmp< DimensionedField< TypeR, GeoMesh, Field > > New(const tmp< DimensionedField< TypeR, GeoMesh, Field >> &tdf1, const word &name, const dimensionSet &dimensions)
dimensionedScalar cos(const dimensionedScalar &ds)
const unitConversion unitDegrees
dictionary dict
randomGenerator rndGen(653213)