FreeFOAM The Cross-Platform CFD Toolkit
InjectionModel.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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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"
28 #include <meshTools/meshTools.H>
29 
30 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
31 
32 template<class CloudType>
34 {
35  IOobject propsDictHeader
36  (
37  "injectionProperties",
38  owner_.db().time().timeName(),
39  "uniform"/cloud::prefix/owner_.name(),
40  owner_.db(),
41  IOobject::MUST_READ,
42  IOobject::NO_WRITE,
43  false
44  );
45 
46  if (propsDictHeader.headerOk())
47  {
48  const IOdictionary propsDict(propsDictHeader);
49 
50  propsDict.readIfPresent("massInjected", massInjected_);
51  propsDict.readIfPresent("nInjections", nInjections_);
52  propsDict.readIfPresent("parcelsAddedTotal", parcelsAddedTotal_);
53  propsDict.readIfPresent("timeStep0", timeStep0_);
54  }
55 }
56 
57 
58 template<class CloudType>
60 {
61  if (owner_.db().time().outputTime())
62  {
63  IOdictionary propsDict
64  (
65  IOobject
66  (
67  "injectionProperties",
68  owner_.db().time().timeName(),
69  "uniform"/cloud::prefix/owner_.name(),
70  owner_.db(),
71  IOobject::NO_READ,
72  IOobject::NO_WRITE,
73  false
74  )
75  );
76 
77  propsDict.add("massInjected", massInjected_);
78  propsDict.add("nInjections", nInjections_);
79  propsDict.add("parcelsAddedTotal", parcelsAddedTotal_);
80  propsDict.add("timeStep0", timeStep0_);
81 
82  propsDict.regIOobject::write();
83  }
84 }
85 
86 
87 template<class CloudType>
89 (
90  const scalar time,
91  label& newParcels,
92  scalar& newVolume
93 )
94 {
95  // Initialise values
96  newParcels = 0;
97  newVolume = 0.0;
98 
99  // Return if not started injection event
100  if (time < SOI_)
101  {
102  timeStep0_ = time;
103  return;
104  }
105 
106  // Make times relative to SOI
107  scalar t0 = timeStep0_ - SOI_;
108  scalar t1 = time - SOI_;
109 
110  // Number of parcels to inject
111  newParcels = parcelsToInject(t0, t1);
112 
113  // Volume of parcels to inject
114  newVolume = volumeToInject(t0, t1);
115 
116  // Hold previous time if no parcels, but non-zero volume fraction
117  if ((newParcels == 0) && (newVolume > 0.0))
118  {
119  // hold value of timeStep0_
120  }
121  else
122  {
123  // advance value of timeStep0_
124  timeStep0_ = time;
125  }
126 }
127 
128 
129 template<class CloudType>
131 (
132  label& cellI,
133  vector& position
134 )
135 {
136  const vector p0 = position;
137 
138  bool foundCell = false;
139 
140  cellI = owner_.mesh().findCell(position);
141 
142  if (cellI >= 0)
143  {
144  const vector& C = owner_.mesh().C()[cellI];
145  position += SMALL*(C - position);
146 
147  foundCell = owner_.mesh().pointInCell(position, cellI);
148  }
149  reduce(foundCell, orOp<bool>());
150 
151  // Last chance - find nearest cell and try that one
152  // - the point is probably on an edge
153  if (!foundCell)
154  {
155  cellI = owner_.mesh().findNearestCell(position);
156 
157  if (cellI >= 0)
158  {
159  const vector& C = owner_.mesh().C()[cellI];
160  position += SMALL*(C - position);
161 
162  foundCell = owner_.mesh().pointInCell(position, cellI);
163  }
164  reduce(foundCell, orOp<bool>());
165  }
166 
167  if (!foundCell)
168  {
170  (
171  "Foam::InjectionModel<CloudType>::findCellAtPosition"
172  "("
173  "label&, "
174  "vector&"
175  ")"
176  )<< "Cannot find parcel injection cell. "
177  << "Parcel position = " << p0 << nl
178  << abort(FatalError);
179  }
180 }
181 
182 
183 template<class CloudType>
185 (
186  const label parcels,
187  const scalar volume,
188  const scalar diameter,
189  const scalar rho
190 )
191 {
192  scalar nP = 0.0;
193  switch (parcelBasis_)
194  {
195  case pbMass:
196  {
197  nP =
198  volume/volumeTotal_
199  *massTotal_/rho
200  /(parcels*mathematicalConstant::pi/6.0*pow3(diameter));
201  break;
202  }
203  case pbNumber:
204  {
205  nP = massTotal_/(rho*volumeTotal_);
206  break;
207  }
208  default:
209  {
210  nP = 0.0;
212  (
213  "Foam::scalar "
214  "Foam::InjectionModel<CloudType>::setNumberOfParticles"
215  "("
216  "const label, "
217  "const scalar, "
218  "const scalar, "
219  "const scalar"
220  ")"
221  )<< "Unknown parcelBasis type" << nl
222  << exit(FatalError);
223  }
224  }
225 
226  return nP;
227 }
228 
229 
230 template<class CloudType>
232 (
233  const label parcelsAdded,
234  const scalar massAdded
235 )
236 {
237  const label allParcelsAdded = returnReduce(parcelsAdded, sumOp<label>());
238 
239  if (allParcelsAdded > 0)
240  {
241  Info<< nl
242  << "--> Cloud: " << owner_.name() << nl
243  << " Added " << allParcelsAdded << " new parcels" << nl << endl;
244  }
245 
246  // Increment total number of parcels added
247  parcelsAddedTotal_ += allParcelsAdded;
248 
249  // Increment total mass injected
250  massInjected_ += returnReduce(massAdded, sumOp<scalar>());
251 
252  // Update time for start of next injection
253  time0_ = owner_.db().time().value();
254 
255  // Increment number of injections
256  nInjections_++;
257 
258  // Write current state to properties file
259  writeProps();
260 }
261 
262 
263 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
264 
265 template<class CloudType>
267 :
268  dict_(dictionary::null),
269  owner_(owner),
270  coeffDict_(dictionary::null),
271  SOI_(0.0),
272  volumeTotal_(0.0),
273  massTotal_(0.0),
274  massInjected_(0.0),
275  nInjections_(0),
276  parcelsAddedTotal_(0),
277  parcelBasis_(pbNumber),
278  time0_(0.0),
279  timeStep0_(0.0)
280 {
281  readProps();
282 }
283 
284 
285 template<class CloudType>
287 (
288  const dictionary& dict,
289  CloudType& owner,
290  const word& type
291 )
292 :
293  dict_(dict),
294  owner_(owner),
295  coeffDict_(dict.subDict(type + "Coeffs")),
296  SOI_(readScalar(coeffDict_.lookup("SOI"))),
297  volumeTotal_(0.0),
298  massTotal_(dimensionedScalar(coeffDict_.lookup("massTotal")).value()),
299  massInjected_(0.0),
300  nInjections_(0),
301  parcelsAddedTotal_(0),
302  parcelBasis_(pbNumber),
303  time0_(owner.db().time().value()),
304  timeStep0_(0.0)
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  word parcelBasisType = coeffDict_.lookup("parcelBasisType");
313  if (parcelBasisType == "mass")
314  {
315  parcelBasis_ = pbMass;
316  }
317  else if (parcelBasisType == "number")
318  {
319  parcelBasis_ = pbNumber;
320  }
321  else
322  {
324  (
325  "Foam::InjectionModel<CloudType>::InjectionModel"
326  "("
327  "const dictionary&, "
328  "CloudType&, "
329  "const word&"
330  ")"
331  )<< "parcelBasisType must be either 'number' or 'mass'" << nl
332  << exit(FatalError);
333  }
334 
335  readProps();
336 }
337 
338 
339 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
340 
341 template<class CloudType>
343 {}
344 
345 
346 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
347 
348 template<class CloudType>
349 template<class TrackData>
351 {
352  if (!active())
353  {
354  return;
355  }
356 
357  const scalar time = owner_.db().time().value();
358  const scalar carrierDt = owner_.db().time().deltaTValue();
359  const polyMesh& mesh = owner_.mesh();
360 
361  // Prepare for next time step
362  label parcelsAdded = 0;
363  scalar massAdded = 0.0;
364  label newParcels = 0;
365  scalar newVolume = 0.0;
366 
367  prepareForNextTimeStep(time, newParcels, newVolume);
368 
369  // Duration of injection period during this timestep
370  const scalar deltaT =
371  max(0.0, min(carrierDt, min(time - SOI_, timeEnd() - time0_)));
372 
373  // Pad injection time if injection starts during this timestep
374  const scalar padTime = max(0.0, SOI_ - time0_);
375 
376  // Introduce new parcels linearly across carrier phase timestep
377  for (label parcelI=0; parcelI<newParcels; parcelI++)
378  {
379  if (validInjection(parcelI))
380  {
381  // Calculate the pseudo time of injection for parcel 'parcelI'
382  scalar timeInj = time0_ + padTime + deltaT*parcelI/newParcels;
383 
384  // Determine the injection position and owner cell
385  label cellI = -1;
387  setPositionAndCell(parcelI, newParcels, timeInj, pos, cellI);
388 
389  if (cellI > -1)
390  {
391  // Lagrangian timestep
392  scalar dt = time - timeInj;
393 
394  // Apply corrections to position for 2-D cases
396 
397  // Create a new parcel
398  parcelType* pPtr = new parcelType(td.cloud(), pos, cellI);
399 
400  // Assign new parcel properties in injection model
401  setProperties(parcelI, newParcels, timeInj, *pPtr);
402 
403  // Check new parcel properties
404  td.cloud().checkParcelProperties(*pPtr, dt, fullyDescribed());
405 
406  // Apply correction to velocity for 2-D cases
408  (
409  mesh,
410  mesh.solutionD(),
411  pPtr->U()
412  );
413 
414  // Number of particles per parcel
415  pPtr->nParticle() =
416  setNumberOfParticles
417  (
418  newParcels,
419  newVolume,
420  pPtr->d(),
421  pPtr->rho()
422  );
423 
424  // Add the new parcel
425  td.cloud().addParticle(pPtr);
426 
427  massAdded += pPtr->nParticle()*pPtr->mass();
428  parcelsAdded++;
429  }
430  }
431  }
432 
433  postInjectCheck(parcelsAdded, massAdded);
434 }
435 
436 
437 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
438 
439 #include "NewInjectionModel.C"
440 
441 // ************************ vim: set sw=4 sts=4 et: ************************ //