KinematicParcel.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) 2011-2014 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 "KinematicParcel.H"
27 #include "forceSuSp.H"
28 #include "IntegrationScheme.H"
29 #include "meshTools.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 template<class ParcelType>
35 
36 
37 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
38 
39 template<class ParcelType>
40 template<class TrackData>
42 (
43  TrackData& td,
44  const scalar dt,
45  const label cellI
46 )
47 {
48  tetIndices tetIs = this->currentTetIndices();
49 
50  rhoc_ = td.rhoInterp().interpolate(this->position(), tetIs);
51 
52  if (rhoc_ < td.cloud().constProps().rhoMin())
53  {
54  if (debug)
55  {
56  WarningIn
57  (
58  "void Foam::KinematicParcel<ParcelType>::setCellValues"
59  "("
60  "TrackData&, "
61  "const scalar, "
62  "const label"
63  ")"
64  ) << "Limiting observed density in cell " << cellI << " to "
65  << td.cloud().constProps().rhoMin() << nl << endl;
66  }
67 
68  rhoc_ = td.cloud().constProps().rhoMin();
69  }
70 
71  Uc_ = td.UInterp().interpolate(this->position(), tetIs);
72 
73  muc_ = td.muInterp().interpolate(this->position(), tetIs);
74 
75  // Apply dispersion components to carrier phase velocity
76  Uc_ = td.cloud().dispersion().update
77  (
78  dt,
79  cellI,
80  U_,
81  Uc_,
82  UTurb_,
83  tTurb_
84  );
85 }
86 
87 
88 template<class ParcelType>
89 template<class TrackData>
91 (
92  TrackData& td,
93  const scalar dt,
94  const label cellI
95 )
96 {
97  Uc_ += td.cloud().UTrans()[cellI]/massCell(cellI);
98 }
99 
100 
101 template<class ParcelType>
102 template<class TrackData>
104 (
105  TrackData& td,
106  const scalar dt,
107  const label cellI
108 )
109 {
110  // Define local properties at beginning of time step
111  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
112  const scalar np0 = nParticle_;
113  const scalar mass0 = mass();
114 
115  // Reynolds number
116  const scalar Re = this->Re(U_, d_, rhoc_, muc_);
117 
118 
119  // Sources
120  //~~~~~~~~
121 
122  // Explicit momentum source for particle
123  vector Su = vector::zero;
124 
125  // Linearised momentum source coefficient
126  scalar Spu = 0.0;
127 
128  // Momentum transfer from the particle to the carrier phase
129  vector dUTrans = vector::zero;
130 
131 
132  // Motion
133  // ~~~~~~
134 
135  // Calculate new particle velocity
136  this->U_ = calcVelocity(td, dt, cellI, Re, muc_, mass0, Su, dUTrans, Spu);
137 
138 
139  // Accumulate carrier phase source terms
140  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
141  if (td.cloud().solution().coupled())
142  {
143  // Update momentum transfer
144  td.cloud().UTrans()[cellI] += np0*dUTrans;
145 
146  // Update momentum transfer coefficient
147  td.cloud().UCoeff()[cellI] += np0*Spu;
148  }
149 }
150 
151 
152 template<class ParcelType>
153 template<class TrackData>
155 (
156  TrackData& td,
157  const scalar dt,
158  const label cellI,
159  const scalar Re,
160  const scalar mu,
161  const scalar mass,
162  const vector& Su,
163  vector& dUTrans,
164  scalar& Spu
165 ) const
166 {
167  typedef typename TrackData::cloudType cloudType;
168  typedef typename cloudType::parcelType parcelType;
169  typedef typename cloudType::forceType forceType;
170 
171  const forceType& forces = td.cloud().forces();
172 
173  // Momentum source due to particle forces
174  const parcelType& p = static_cast<const parcelType&>(*this);
175  const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
176  const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
177  const forceSuSp Feff = Fcp + Fncp;
178  const scalar massEff = forces.massEff(p, mass);
179 
180 
181  // New particle velocity
182  //~~~~~~~~~~~~~~~~~~~~~~
183 
184  // Update velocity - treat as 3-D
185  const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff;
186  const scalar bp = Feff.Sp()/massEff;
187 
188  Spu = dt*Feff.Sp();
189 
191  td.cloud().UIntegrator().integrate(U_, dt, abp, bp);
192 
193  vector Unew = Ures.value();
194 
195  // note: Feff.Sp() and Fc.Sp() must be the same
196  dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su());
197 
198  // Apply correction to velocity and dUTrans for reduced-D cases
199  const polyMesh& mesh = td.cloud().pMesh();
200  meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
201  meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);
202 
203  return Unew;
204 }
205 
206 
207 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
208 
209 template<class ParcelType>
211 (
213 )
214 :
215  ParcelType(p),
216  active_(p.active_),
217  typeId_(p.typeId_),
218  nParticle_(p.nParticle_),
219  d_(p.d_),
220  dTarget_(p.dTarget_),
221  U_(p.U_),
222  rho_(p.rho_),
223  age_(p.age_),
224  tTurb_(p.tTurb_),
225  UTurb_(p.UTurb_),
226  rhoc_(p.rhoc_),
227  Uc_(p.Uc_),
228  muc_(p.muc_)
229 {}
230 
231 
232 template<class ParcelType>
234 (
236  const polyMesh& mesh
237 )
238 :
239  ParcelType(p, mesh),
240  active_(p.active_),
241  typeId_(p.typeId_),
242  nParticle_(p.nParticle_),
243  d_(p.d_),
244  dTarget_(p.dTarget_),
245  U_(p.U_),
246  rho_(p.rho_),
247  age_(p.age_),
248  tTurb_(p.tTurb_),
249  UTurb_(p.UTurb_),
250  rhoc_(p.rhoc_),
251  Uc_(p.Uc_),
252  muc_(p.muc_)
253 {}
254 
255 
256 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
257 
258 template<class ParcelType>
259 template<class TrackData>
261 (
262  TrackData& td,
263  const scalar trackTime
264 )
265 {
266  typename TrackData::cloudType::parcelType& p =
267  static_cast<typename TrackData::cloudType::parcelType&>(*this);
268 
269  td.switchProcessor = false;
270  td.keepParticle = true;
271 
272  const polyMesh& mesh = td.cloud().pMesh();
273  const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
274  const scalarField& cellLengthScale = td.cloud().cellLengthScale();
275  const scalar maxCo = td.cloud().solution().maxCo();
276 
277  scalar tEnd = (1.0 - p.stepFraction())*trackTime;
278  scalar dtMax = trackTime;
279  if (td.cloud().solution().transient())
280  {
281  dtMax *= maxCo;
282  }
283 
284  bool tracking = true;
285  label nTrackingStalled = 0;
286 
287  while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
288  {
289  // Apply correction to position for reduced-D cases
290  meshTools::constrainToMeshCentre(mesh, p.position());
291 
292  const point start(p.position());
293 
294  // Set the Lagrangian time-step
295  scalar dt = min(dtMax, tEnd);
296 
297  // Cache the parcel current cell as this will change if a face is hit
298  const label cellI = p.cell();
299 
300  const scalar magU = mag(U_);
301  if (p.active() && tracking && (magU > ROOTVSMALL))
302  {
303  const scalar d = dt*magU;
304  const scalar dCorr = min(d, maxCo*cellLengthScale[cellI]);
305  dt *=
306  dCorr/d
307  *p.trackToFace(p.position() + dCorr*U_/magU, td);
308  }
309 
310  tEnd -= dt;
311 
312  scalar newStepFraction = 1.0 - tEnd/trackTime;
313 
314  if (tracking)
315  {
316  if
317  (
318  mag(p.stepFraction() - newStepFraction)
319  < particle::minStepFractionTol
320  )
321  {
322  nTrackingStalled++;
323 
324  if (nTrackingStalled > maxTrackAttempts)
325  {
326  tracking = false;
327  }
328  }
329  else
330  {
331  nTrackingStalled = 0;
332  }
333  }
334 
335  p.stepFraction() = newStepFraction;
336 
337  bool calcParcel = true;
338  if (!tracking && td.cloud().solution().steadyState())
339  {
340  calcParcel = false;
341  }
342 
343  // Avoid problems with extremely small timesteps
344  if ((dt > ROOTVSMALL) && calcParcel)
345  {
346  // Update cell based properties
347  p.setCellValues(td, dt, cellI);
348 
349  if (td.cloud().solution().cellValueSourceCorrection())
350  {
351  p.cellValueSourceCorrection(td, dt, cellI);
352  }
353 
354  p.calc(td, dt, cellI);
355  }
356 
357  if (p.onBoundary() && td.keepParticle)
358  {
359  if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
360  {
361  td.switchProcessor = true;
362  }
363  }
364 
365  p.age() += dt;
366 
367  td.cloud().functions().postMove(p, cellI, dt, start, td.keepParticle);
368  }
369 
370  return td.keepParticle;
371 }
372 
373 
374 template<class ParcelType>
375 template<class TrackData>
377 {
378  typename TrackData::cloudType::parcelType& p =
379  static_cast<typename TrackData::cloudType::parcelType&>(*this);
380 
381  td.cloud().functions().postFace(p, p.face(), td.keepParticle);
382 }
383 
384 
385 template<class ParcelType>
387 {}
388 
389 
390 template<class ParcelType>
391 template<class TrackData>
393 (
394  const polyPatch& pp,
395  TrackData& td,
396  const label patchI,
397  const scalar trackFraction,
398  const tetIndices& tetIs
399 )
400 {
401  typename TrackData::cloudType::parcelType& p =
402  static_cast<typename TrackData::cloudType::parcelType&>(*this);
403 
404  // Invoke post-processing model
405  td.cloud().functions().postPatch
406  (
407  p,
408  pp,
409  trackFraction,
410  tetIs,
411  td.keepParticle
412  );
413 
414  // Invoke surface film model
415  if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))
416  {
417  // All interactions done
418  return true;
419  }
420  else
421  {
422  // Invoke patch interaction model
423  return td.cloud().patchInteraction().correct
424  (
425  p,
426  pp,
427  td.keepParticle,
428  trackFraction,
429  tetIs
430  );
431  }
432 }
433 
434 
435 template<class ParcelType>
436 template<class TrackData>
438 (
439  const processorPolyPatch&,
440  TrackData& td
441 )
442 {
443  td.switchProcessor = true;
444 }
445 
446 
447 template<class ParcelType>
448 template<class TrackData>
450 (
451  const wallPolyPatch& wpp,
452  TrackData& td,
453  const tetIndices&
454 )
455 {
456  // Wall interactions handled by generic hitPatch function
457 }
458 
459 
460 template<class ParcelType>
461 template<class TrackData>
463 (
464  const polyPatch&,
465  TrackData& td
466 )
467 {
468  td.keepParticle = false;
469 }
470 
471 
472 template<class ParcelType>
474 {
475  ParcelType::transformProperties(T);
476 
477  U_ = transform(T, U_);
478 }
479 
480 
481 template<class ParcelType>
483 (
484  const vector& separation
485 )
486 {
487  ParcelType::transformProperties(separation);
488 }
489 
490 
491 template<class ParcelType>
493 (
494  const vector&
495 ) const
496 {
497  return 0.5*d_;
498 }
499 
500 
501 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
502 
503 #include "KinematicParcelIO.C"
504 
505 // ************************************************************************* //
void hitProcessorPatch(const processorPolyPatch &, TrackData &td)
Overridable function to handle the particle hitting a.
scalar d_
Diameter [m].
vector Uc_
Velocity [m/s].
void hitFace(int &td)
Overridable function to handle the particle hitting a face.
bool active_
Active flag - tracking inactive when active = false.
vector UTurb_
Turbulent velocity fluctuation [m/s].
vector U_
Velocity of Parcel [m/s].
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
scalar age_
Age [s].
This function object calculates the forces and moments by integrating the pressure and skin-friction ...
Definition: forces.H:206
Helper class to supply results of integration.
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
scalar Sp() const
Return const access to the implicit coefficient [kg/s].
Definition: forceSuSpI.H:62
scalar dTarget_
Target diameter [m].
forces(const forces &)
Disallow default bitwise copy construct.
const volScalarField & T
Definition: createFields.H:25
scalar rho_
Density [kg/m3].
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
scalar tTurb_
Time spent in turbulent eddy [s].
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
const vector & Su() const
Return const access to the explicit contribution [kg.m/s2].
Definition: forceSuSpI.H:56
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Neighbour processor patch.
#define WarningIn(functionName)
Report a warning using Foam::Warning.
void constrainDirection(const polyMesh &mesh, const Vector< label > &dirs, vector &d)
Set the constrained components of directions/velocity to zero.
Definition: meshTools.C:693
Foam::polyBoundaryMesh.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
label typeId_
Parcel type id.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
void calc(TrackData &td, const scalar dt, const label cellI)
Update parcel properties over the time interval.
Type average() const
Return const access to the average.
void setCellValues(TrackData &td, const scalar dt, const label cellI)
Set cell values.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Foam::wallPolyPatch.
Definition: wallPolyPatch.H:48
Type value() const
Return const access to the value.
bool hitPatch(const polyPatch &p, TrackData &td, const label patchI, const scalar trackFraction, const tetIndices &tetIs)
Overridable function to handle the particle hitting a patch.
scalar nParticle_
Number of particles in Parcel.
void constrainToMeshCentre(const polyMesh &mesh, point &pt)
Set the constrained components of position to mesh centre.
Definition: meshTools.C:634
bool move(TrackData &td, const scalar trackTime)
Move the parcel.
Kinematic parcel class with rotational motion (as spherical particles only) and one/two-way coupling ...
const vector calcVelocity(TrackData &td, const scalar dt, const label cellI, const scalar Re, const scalar mu, const scalar mass, const vector &Su, vector &dUTrans, scalar &Spu) const
Calculate new particle velocity.
void hitWallPatch(const wallPolyPatch &, TrackData &td, const tetIndices &)
Overridable function to handle the particle hitting a wallPatch.
KinematicParcel(const polyMesh &mesh, const vector &position, const label cellI, const label tetFaceI, const label tetPtI)
Construct from owner, position, and cloud owner.
void cellValueSourceCorrection(TrackData &td, const scalar dt, const label cellI)
Correct cell values using latest transfer information.
scalar rhoc_
Density [kg/m3].
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual scalar wallImpactDistance(const vector &n) const
The nearest distance to a wall that the particle can be.
scalar maxCo
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:73
scalar muc_
Viscosity [Pa.s].
Helper container for force Su and Sp terms.
Definition: forceSuSp.H:58
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97