ConeInjection.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2015-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "ConeInjection.H"
30 #include "Function1.H"
31 #include "mathematicalConstants.H"
32 #include "unitConversion.H"
33 
34 using namespace Foam::constant::mathematical;
35 
36 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
37 
38 template<class CloudType>
40 (
41  const dictionary& dict,
42  CloudType& owner,
43  const word& modelName
44 )
45 :
46  InjectionModel<CloudType>(dict, owner, modelName, typeName),
47  positionAxis_(this->coeffDict().lookup("positionAxis")),
48  injectorCells_(positionAxis_.size()),
49  injectorTetFaces_(positionAxis_.size()),
50  injectorTetPts_(positionAxis_.size()),
51  duration_(this->coeffDict().getScalar("duration")),
52  parcelsPerInjector_
53  (
54  this->coeffDict().getScalar("parcelsPerInjector")
55  ),
56  flowRateProfile_
57  (
59  (
60  "flowRateProfile",
61  this->coeffDict(),
62  &owner.mesh()
63  )
64  ),
65  Umag_
66  (
68  (
69  "Umag",
70  this->coeffDict(),
71  &owner.mesh()
72  )
73  ),
74  thetaInner_
75  (
77  (
78  "thetaInner",
79  this->coeffDict(),
80  &owner.mesh()
81  )
82  ),
83  thetaOuter_
84  (
86  (
87  "thetaOuter",
88  this->coeffDict(),
89  &owner.mesh()
90  )
91  ),
92  sizeDistribution_
93  (
95  (
96  this->coeffDict().subDict("sizeDistribution"), owner.rndGen()
97  )
98  ),
99  nInjected_(this->parcelsAddedTotal()),
100  injectorOrder_(identity(positionAxis_.size())),
101  tanVec1_(),
102  tanVec2_()
103 {
104  updateMesh();
105 
106  tanVec1_.setSize(positionAxis_.size());
107  tanVec2_.setSize(positionAxis_.size());
108 
109  // Convert from user time to reduce the number of time conversion calls
110  const Time& time = owner.db().time();
111  duration_ = time.userTimeToTime(duration_);
112  flowRateProfile_->userTimeToTime(time);
113  Umag_->userTimeToTime(time);
114  thetaInner_->userTimeToTime(time);
115  thetaOuter_->userTimeToTime(time);
116 
117  // Normalise direction vector and determine direction vectors
118  // tangential to injector axis direction
119  forAll(positionAxis_, i)
120  {
121  vector& axis = positionAxis_[i].second();
122  axis.normalise();
123 
124  vector tangent = Zero;
125  scalar magTangent = 0.0;
126 
127  Random& rnd = this->owner().rndGen();
128  while (magTangent < SMALL)
129  {
130  vector v = rnd.sample01<vector>();
131 
132  tangent = v - (v & axis)*axis;
133  magTangent = mag(tangent);
134  }
135 
136  tanVec1_[i] = tangent/magTangent;
137  tanVec2_[i] = axis^tanVec1_[i];
138  }
139 
140  // Set total volume to inject
141  this->volumeTotal_ = flowRateProfile_->integrate(0.0, duration_);
142 }
143 
144 
145 template<class CloudType>
147 (
148  const ConeInjection<CloudType>& im
149 )
150 :
152  positionAxis_(im.positionAxis_),
153  injectorCells_(im.injectorCells_),
154  injectorTetFaces_(im.injectorTetFaces_),
155  injectorTetPts_(im.injectorTetPts_),
156  duration_(im.duration_),
157  parcelsPerInjector_(im.parcelsPerInjector_),
158  flowRateProfile_(im.flowRateProfile_.clone()),
159  Umag_(im.Umag_.clone()),
160  thetaInner_(im.thetaInner_.clone()),
161  thetaOuter_(im.thetaOuter_.clone()),
162  sizeDistribution_(im.sizeDistribution_.clone()),
163  nInjected_(im.nInjected_),
164  injectorOrder_(im.injectorOrder_),
165  tanVec1_(im.tanVec1_),
166  tanVec2_(im.tanVec2_)
167 {}
168 
169 
170 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
171 
172 template<class CloudType>
174 {
175  bitSet reject(positionAxis_.size());
176 
177  // Set/cache the injector cells
178  forAll(positionAxis_, i)
179  {
180  if
181  (
182  !this->findCellAtPosition
183  (
184  injectorCells_[i],
185  injectorTetFaces_[i],
186  injectorTetPts_[i],
187  positionAxis_[i].first(),
188  !this->ignoreOutOfBounds_
189 
190  )
191  )
192  {
193  reject.set(i);
194  }
195  }
196 
197  const label nRejected = reject.count();
198 
199  if (nRejected)
200  {
201  reject.flip();
202  inplaceSubset(reject, injectorCells_);
203  inplaceSubset(reject, injectorTetFaces_);
204  inplaceSubset(reject, injectorTetPts_);
205  inplaceSubset(reject, positionAxis_);
206 
207  Info<< " " << nRejected
208  << " positions rejected, out of bounds" << endl;
209  }
210 }
211 
212 
213 template<class CloudType>
215 {
216  return this->SOI_ + duration_;
217 }
218 
219 
220 template<class CloudType>
222 (
223  const scalar time0,
224  const scalar time1
225 )
226 {
227  if ((time0 >= 0.0) && (time0 < duration_))
228  {
229  const scalar targetVolume = flowRateProfile_->integrate(0, time1);
230 
231  const scalar volumeFraction = targetVolume/this->volumeTotal_;
232 
233  const label targetParcels =
234  ceil(positionAxis_.size()*parcelsPerInjector_*volumeFraction);
235 
236  return targetParcels - nInjected_;
237  }
238 
239  return 0;
240 }
241 
242 
243 template<class CloudType>
245 (
246  const scalar time0,
247  const scalar time1
248 )
249 {
250  if ((time0 >= 0.0) && (time0 < duration_))
251  {
252  return flowRateProfile_->integrate(time0, time1);
253  }
254 
255  return 0.0;
256 }
257 
258 
259 template<class CloudType>
261 (
262  const label parcelI,
263  const label,
264  const scalar,
265  vector& position,
266  label& cellOwner,
267  label& tetFacei,
268  label& tetPti
269 )
270 {
271  Random& rnd = this->owner().rndGen();
272  rnd.shuffle(injectorOrder_);
273 
274  const label i = injectorOrder_[parcelI % positionAxis_.size()];
275  position = positionAxis_[i].first();
276  cellOwner = injectorCells_[i];
277  tetFacei = injectorTetFaces_[i];
278  tetPti = injectorTetPts_[i];
279 }
280 
281 
282 template<class CloudType>
284 (
285  const label parcelI,
286  const label,
287  const scalar time,
288  typename CloudType::parcelType& parcel
289 )
290 {
291  Random& rnd = this->owner().rndGen();
292 
293  const label i = injectorOrder_[parcelI % positionAxis_.size()];
294 
295  // Set direction vectors for position i
296  scalar t = time - this->SOI_;
297  scalar ti = thetaInner_->value(t);
298  scalar to = thetaOuter_->value(t);
299  scalar coneAngle = degToRad(rnd.position<scalar>(ti, to));
300 
301  scalar alpha = sin(coneAngle);
302  scalar dcorr = cos(coneAngle);
303  scalar beta = twoPi*rnd.sample01<scalar>();
304 
305  vector normal = alpha*(tanVec1_[i]*cos(beta) + tanVec2_[i]*sin(beta));
306  vector dirVec = dcorr*positionAxis_[i].second();
307  dirVec += normal;
308  dirVec.normalise();
309 
310  // Set particle velocity
311  parcel.U() = Umag_->value(t)*dirVec;
312 
313  // Set particle diameter
314  parcel.d() = sizeDistribution_().sample();
315 
316  // Increment number of particles injected
317  nInjected_++;
318 }
319 
320 
321 template<class CloudType>
323 {
324  return false;
325 }
326 
327 
328 template<class CloudType>
330 {
331  return true;
332 }
333 
334 
335 // ************************************************************************* //
Foam::Random
Random number generator.
Definition: Random.H:59
mathematicalConstants.H
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:73
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::bitSet
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:63
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:212
Foam::sin
dimensionedScalar sin(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:264
Function1.H
Foam::Random::sample01
Type sample01()
Return a sample whose components lie in the range [0,1].
Definition: RandomTemplates.C:36
Foam::TimeState::userTimeToTime
virtual scalar userTimeToTime(const scalar theta) const
Convert the user-time (e.g. CA deg) to real-time (s).
Definition: TimeState.C:49
Foam::ConeInjection::validInjection
virtual bool validInjection(const label parcelI)
Return flag to identify whether or not injection of parcelI is.
Definition: ConeInjection.C:329
Foam::Random::shuffle
void shuffle(UList< Type > &values)
Shuffle the values in the list.
Definition: RandomTemplates.C:83
unitConversion.H
Unit conversion functions.
Foam::InjectionModel
Templated injection model class.
Definition: InjectionModel.H:73
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::Vector::normalise
Vector< Cmpt > & normalise()
Normalise the vector by its magnitude.
Definition: VectorI.H:123
Foam::Function1
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: propellerInfo.H:291
Foam::DSMCCloud::rndGen
Random & rndGen()
Return reference to the random object.
Definition: DSMCCloudI.H:124
Foam::ConeInjection::setPositionAndCell
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.
Definition: ConeInjection.C:261
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
ConeInjection.H
Foam::ConeInjection::volumeToInject
virtual scalar volumeToInject(const scalar time0, const scalar time1)
Volume of parcels to introduce relative to SOI.
Definition: ConeInjection.C:245
Foam::DSMCCloud::mesh
const fvMesh & mesh() const
Return reference to the mesh.
Definition: DSMCCloudI.H:44
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
beta
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
Foam::ConeInjection::fullyDescribed
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
Definition: ConeInjection.C:322
Foam::ConeInjection::ConeInjection
ConeInjection(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
Definition: ConeInjection.C:40
Foam::radiation::lookup
Lookup type of boundary radiation properties.
Definition: lookup.H:63
Foam::DSMCCloud
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:71
Foam::constant::mathematical::twoPi
constexpr scalar twoPi(2 *M_PI)
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::inplaceSubset
void inplaceSubset(const BoolListType &select, ListType &input, const bool invert=false)
Inplace extract elements of the input list when select is true.
Definition: ListOpsTemplates.C:583
Foam::degToRad
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Definition: unitConversion.H:48
Foam::distributionModel::New
static autoPtr< distributionModel > New(const dictionary &dict, Random &rndGen)
Selector.
Definition: distributionModelNew.C:34
Foam::constant::mathematical
Mathematical constants.
Foam::ConeInjection::parcelsToInject
virtual label parcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
Definition: ConeInjection.C:222
Foam::ConeInjection
Multi-point cone injection model.
Definition: ConeInjection.H:68
Foam::Random::position
Type position(const Type &start, const Type &end)
Return a sample on the interval [start,end].
Definition: RandomTemplates.C:62
Foam::Vector< scalar >
Foam::ConeInjection::setProperties
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, typename CloudType::parcelType &parcel)
Set the parcel properties.
Definition: ConeInjection.C:284
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::identity
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
Definition: labelList.C:38
Foam::ConeInjection::updateMesh
virtual void updateMesh()
Set injector locations when mesh is updated.
Definition: ConeInjection.C:173
Foam::DSMCCloud::parcelType
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:220
Foam::ConeInjection::timeEnd
scalar timeEnd() const
Return the end-of-injection time.
Definition: ConeInjection.C:214
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:265