MomentumParcel.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-2022 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 "MomentumParcel.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 TrackCloudType>
42 (
43  TrackCloudType& cloud,
44  trackingData& td
45 )
46 {
47  tetIndices tetIs = this->currentTetIndices(td.mesh);
48 
49  td.rhoc() = td.rhoInterp().interpolate(this->coordinates(), tetIs);
50 
51  if (td.rhoc() < cloud.constProps().rhoMin())
52  {
53  if (debug)
54  {
56  << "Limiting observed density in cell " << this->cell()
57  << " to " << cloud.constProps().rhoMin() << nl << endl;
58  }
59 
60  td.rhoc() = cloud.constProps().rhoMin();
61  }
62 
63  td.Uc() = td.UInterp().interpolate(this->coordinates(), tetIs);
64 
65  td.muc() = td.muInterp().interpolate(this->coordinates(), tetIs);
66 }
67 
68 
69 template<class ParcelType>
70 template<class TrackCloudType>
72 (
73  TrackCloudType& cloud,
74  trackingData& td,
75  const scalar dt
76 )
77 {
78  td.Uc() = cloud.dispersion().update
79  (
80  dt,
81  this->cell(),
82  U_,
83  td.Uc(),
84  UTurb_,
85  tTurb_
86  );
87 }
88 
89 
90 template<class ParcelType>
91 template<class TrackCloudType>
93 (
94  TrackCloudType& cloud,
95  trackingData& td,
96  const scalar dt
97 )
98 {
99  td.Uc() += cloud.UTransRef()[this->cell()]/massCell(td);
100 }
101 
102 
103 template<class ParcelType>
104 template<class TrackCloudType>
106 (
107  TrackCloudType& cloud,
108  trackingData& td,
109  const scalar dt
110 )
111 {
112  // Define local properties at beginning of time step
113  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
114  const scalar np0 = nParticle_;
115  const scalar mass0 = mass();
116 
117  // Reynolds number
118  const scalar Re = this->Re(td);
119 
120 
121  // Sources
122  //~~~~~~~~
123 
124  // Explicit momentum source for particle
125  vector Su = Zero;
126 
127  // Linearised momentum source coefficient
128  scalar Spu = 0.0;
129 
130  // Momentum transfer from the particle to the carrier phase
131  vector dUTrans = Zero;
132 
133 
134  // Motion
135  // ~~~~~~
136 
137  // Calculate new particle velocity
138  this->U_ =
139  calcVelocity(cloud, td, dt, Re, td.muc(), mass0, Su, dUTrans, Spu);
140 
141 
142  // Accumulate carrier phase source terms
143  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
144  if (cloud.solution().coupled())
145  {
146  // Update momentum transfer
147  cloud.UTransRef()[this->cell()] += np0*dUTrans;
148 
149  // Update momentum transfer coefficient
150  cloud.UCoeffRef()[this->cell()] += np0*Spu;
151  }
152 }
153 
154 
155 template<class ParcelType>
156 template<class TrackCloudType>
158 (
159  TrackCloudType& cloud,
160  trackingData& td,
161  const scalar dt,
162  const scalar Re,
163  const scalar mu,
164  const scalar mass,
165  const vector& Su,
166  vector& dUTrans,
167  scalar& Spu
168 ) const
169 {
170  const typename TrackCloudType::parcelType& p =
171  static_cast<const typename TrackCloudType::parcelType&>(*this);
172  typename TrackCloudType::parcelType::trackingData& ttd =
173  static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
174 
175  const typename TrackCloudType::forceType& forces = cloud.forces();
176 
177  // Momentum source due to particle forces
178  const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, mu);
179  const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, mu);
180  const scalar massEff = forces.massEff(p, ttd, mass);
181 
182  /*
183  // Proper splitting ...
184  // Calculate the integration coefficients
185  const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff;
186  const vector ancp = (Fncp.Sp()*td.Uc() + Fncp.Su() + Su)/massEff;
187  const scalar bcp = Fcp.Sp()/massEff;
188  const scalar bncp = Fncp.Sp()/massEff;
189 
190  // Integrate to find the new parcel velocity
191  const vector deltaUcp =
192  cloud.UIntegrator().partialDelta
193  (
194  U_, dt, acp + ancp, bcp + bncp, acp, bcp
195  );
196  const vector deltaUncp =
197  cloud.UIntegrator().partialDelta
198  (
199  U_, dt, acp + ancp, bcp + bncp, ancp, bncp
200  );
201  const vector deltaT = deltaUcp + deltaUncp;
202  */
203 
204  // Shortcut splitting assuming no implicit non-coupled force ...
205  // Calculate the integration coefficients
206  const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff;
207  const vector ancp = (Fncp.Su() + Su)/massEff;
208  const scalar bcp = Fcp.Sp()/massEff;
209 
210  // Integrate to find the new parcel velocity
211  const vector deltaU = cloud.UIntegrator().delta(U_, dt, acp + ancp, bcp);
212  const vector deltaUncp = ancp*dt;
213  const vector deltaUcp = deltaU - deltaUncp;
214 
215  // Calculate the new velocity and the momentum transfer terms
216  vector Unew = U_ + deltaU;
217 
218  dUTrans -= massEff*deltaUcp;
219 
220  Spu = dt*Fcp.Sp();
221 
222  // Apply correction to velocity and dUTrans for reduced-D cases
223  const polyMesh& mesh = cloud.pMesh();
224  meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
225  meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);
226 
227  return Unew;
228 }
229 
230 
231 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
232 
233 template<class ParcelType>
235 (
237 )
238 :
239  ParcelType(p),
240  moving_(p.moving_),
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 {}
251 
252 
253 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
254 
255 template<class ParcelType>
256 template<class TrackCloudType>
258 (
259  TrackCloudType& cloud,
260  trackingData& td
261 )
262 {
263  typename TrackCloudType::parcelType& p =
264  static_cast<typename TrackCloudType::parcelType&>(*this);
265  typename TrackCloudType::parcelType::trackingData& ttd =
266  static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
267 
268  ttd.keepParticle = true;
269  td.sendToProc = -1;
270 
271  const scalarField& cellLengthScale = cloud.cellLengthScale();
272  const scalar maxCo = cloud.solution().maxCo();
273 
274  while
275  (
276  ttd.keepParticle
277  && ttd.sendToProc == -1
278  && p.stepFraction() < ttd.stepFractionRange().second()
279  )
280  {
281  if (p.moving() && p.onFace())
282  {
283  cloud.functions().postFace(p);
284  }
285 
286  // Cache the current position, cell and step-fraction
287  const point start = p.position(td.mesh);
288  const scalar sfrac = p.stepFraction();
289 
290  // Total displacement over the time-step
291  const vector s = ttd.trackTime()*U_;
292 
293  // Cell length scale
294  const scalar l = cellLengthScale[p.cell()];
295 
296  // Deviation from the mesh centre for reduced-D cases
297  const vector d = p.deviationFromMeshCentre(td.mesh);
298 
299  // Fraction of the displacement to track in this loop. This is limited
300  // to ensure that the both the time and distance tracked is less than
301  // maxCo times the total value.
302  scalar f = ttd.stepFractionRange().second() - p.stepFraction();
303  f = min(f, maxCo);
304  f = min(f, maxCo/min(max(mag(s)/l, rootSmall), rootGreat));
305  if (p.moving())
306  {
307  // Track to the next face
308  p.trackToFace(td.mesh, f*s - d, f);
309  }
310  else
311  {
312  // At present the only thing that sets moving_ to false is a stick
313  // wall interaction. We want the position of the particle to remain
314  // the same relative to the face that it is on. The local
315  // coordinates therefore do not change. We still advance in time and
316  // perform the relevant interactions with the fixed particle.
317  p.stepFraction() += f;
318  }
319 
320  const scalar dt = (p.stepFraction() - sfrac)*ttd.trackTime();
321 
322  // Avoid problems with extremely small timesteps
323  if (dt > rootVSmall)
324  {
325  // Update cell based properties
326  p.setCellValues(cloud, ttd);
327 
328  p.calcDispersion(cloud, ttd, dt);
329 
330  if (cloud.solution().cellValueSourceCorrection())
331  {
332  p.cellValueSourceCorrection(cloud, ttd, dt);
333  }
334 
335  p.calc(cloud, ttd, dt);
336  }
337 
338  p.age() += dt;
339 
340  cloud.functions().postMove(p, dt, start, ttd.keepParticle);
341 
342  if (p.moving() && p.onFace() && ttd.keepParticle)
343  {
344  cloud.functions().preFace(p);
345 
346  p.hitFace(f*s - d, f, cloud, ttd);
347  }
348  }
349 
350  return ttd.keepParticle;
351 }
352 
353 
354 template<class ParcelType>
356 (
357  const transformer& transform
358 )
359 {
360  ParcelType::transformProperties(transform);
361 
362  U_ = transform.transform(U_);
363 }
364 
365 
366 template<class ParcelType>
367 template<class TrackCloudType>
369 (
370  TrackCloudType& cloud,
371  trackingData& td
372 )
373 {
374  ParcelType::correctAfterParallelTransfer(cloud, td);
375 
376  typename TrackCloudType::parcelType& p =
377  static_cast<typename TrackCloudType::parcelType&>(*this);
378 
379  const polyPatch& pp = td.mesh.boundaryMesh()[td.sendToPatch];
380 
381  cloud.functions().postPatch(p, pp);
382 }
383 
384 
385 template<class ParcelType>
386 template<class TrackCloudType>
388 (
389  TrackCloudType& cloud,
390  trackingData& td
391 )
392 {
393  typename TrackCloudType::parcelType& p =
394  static_cast<typename TrackCloudType::parcelType&>(*this);
395 
396  const polyPatch& pp = td.mesh.boundaryMesh()[p.patch(td.mesh)];
397 
398  // Allow a surface film model to consume the parcel
399  if (cloud.surfaceFilm().transferParcel(p, pp, td.keepParticle))
400  {
401  cloud.functions().postPatch(p, pp);
402  return true;
403  }
404 
405  // Pass to the patch interaction model
406  if (cloud.patchInteraction().correct(p, pp, td.keepParticle))
407  {
408  cloud.functions().postPatch(p, pp);
409  return true;
410  }
411 
412  return false;
413 }
414 
415 
416 template<class ParcelType>
417 template<class TrackCloudType>
419 (
420  TrackCloudType& cloud,
421  trackingData& td
422 )
423 {
424  ParcelType::hitWedgePatch(cloud, td);
425 
426  typename TrackCloudType::parcelType& p =
427  static_cast<typename TrackCloudType::parcelType&>(*this);
428 
429  const polyPatch& pp = td.mesh.boundaryMesh()[p.patch(td.mesh)];
430 
431  cloud.functions().postPatch(p, pp);
432 }
433 
434 
435 template<class ParcelType>
436 template<class TrackCloudType>
438 (
439  TrackCloudType& cloud,
440  trackingData& td
441 )
442 {
443  ParcelType::hitSymmetryPlanePatch(cloud, td);
444 
445  typename TrackCloudType::parcelType& p =
446  static_cast<typename TrackCloudType::parcelType&>(*this);
447 
448  const polyPatch& pp = td.mesh.boundaryMesh()[p.patch(td.mesh)];
449 
450  cloud.functions().postPatch(p, pp);
451 }
452 
453 
454 template<class ParcelType>
455 template<class TrackCloudType>
457 (
458  TrackCloudType& cloud,
459  trackingData& td
460 )
461 {
462  ParcelType::hitSymmetryPatch(cloud, td);
463 
464  typename TrackCloudType::parcelType& p =
465  static_cast<typename TrackCloudType::parcelType&>(*this);
466 
467  const polyPatch& pp = td.mesh.boundaryMesh()[p.patch(td.mesh)];
468 
469  cloud.functions().postPatch(p, pp);
470 }
471 
472 
473 template<class ParcelType>
474 template<class TrackCloudType>
476 (
477  TrackCloudType& cloud,
478  trackingData& td
479 )
480 {
481  ParcelType::hitCyclicPatch(cloud, td);
482 
483  typename TrackCloudType::parcelType& p =
484  static_cast<typename TrackCloudType::parcelType&>(*this);
485 
486  const polyPatch& pp = td.mesh.boundaryMesh()[p.patch(td.mesh)];
487 
488  cloud.functions().postPatch(p, pp);
489 }
490 
491 
492 template<class ParcelType>
493 template<class TrackCloudType>
495 (
496  const vector& displacement,
497  const scalar fraction,
498  const label patchi,
499  TrackCloudType& cloud,
500  trackingData& td
501 )
502 {
503  const bool result =
504  ParcelType::hitNonConformalCyclicPatch
505  (
506  displacement,
507  fraction,
508  patchi,
509  cloud,
510  td
511  );
512 
513  if (td.sendToProc == -1 && result)
514  {
515  typename TrackCloudType::parcelType& p =
516  static_cast<typename TrackCloudType::parcelType&>(*this);
517 
518  const polyPatch& pp = td.mesh.boundaryMesh()[patchi];
519 
520  cloud.functions().postPatch(p, pp);
521  }
522 
523  return result;
524 }
525 
526 
527 template<class ParcelType>
528 template<class TrackCloudType>
530 (
531  TrackCloudType&,
532  trackingData& td
533 )
534 {
535  const polyPatch& pp = td.mesh.boundaryMesh()[this->patch(td.mesh)];
536 
538  << "Particle " << this->origId() << " hit " << pp.type() << " patch "
539  << pp.name() << " at " << this->position(td.mesh)
540  << " but no interaction model was specified for this patch"
541  << exit(FatalError);
542 }
543 
544 
545 template<class ParcelType>
546 template<class TrackCloudType>
548 (
549  TrackCloudType& cloud,
550  trackingData& td
551 )
552 {
553  ParcelType::hitBasicPatch(cloud, td);
554 
555  typename TrackCloudType::parcelType& p =
556  static_cast<typename TrackCloudType::parcelType&>(*this);
557 
558  const polyPatch& pp = td.mesh.boundaryMesh()[p.patch(td.mesh)];
559 
560  cloud.functions().postPatch(p, pp);
561 }
562 
563 
564 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
565 
566 #include "MomentumParcelIO.C"
567 
568 // ************************************************************************* //
const vector & Uc() const
Return the continuous phase velocity.
scalar muc() const
Return the continuous phase viscosity.
const interpolation< vector > & UInterp() const
Return const access to the interpolator for continuous.
scalar rhoc() const
Return the continuous phase density.
const interpolation< scalar > & muInterp() const
Return const access to the interpolator for continuous.
const interpolation< scalar > & rhoInterp() const
Return const access to the interpolator for continuous.
Momentum parcel class with rotational motion (as spherical particles only) and one/two-way coupling w...
void calcDispersion(TrackCloudType &cloud, trackingData &td, const scalar dt)
Apply dispersion to the carrier phase velocity and update.
bool move(TrackCloudType &cloud, trackingData &td)
Move the parcel.
void cellValueSourceCorrection(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct cell values using latest transfer information.
void hitCyclicPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a.
const vector calcVelocity(TrackCloudType &cloud, trackingData &td, const scalar dt, const scalar Re, const scalar mu, const scalar mass, const vector &Su, vector &dUTrans, scalar &Spu) const
Calculate new particle velocity.
void hitSymmetryPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a.
void hitSymmetryPlanePatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a.
MomentumParcel(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti, const label facei)
Construct from mesh, coordinates and topology.
bool hitPatch(TrackCloudType &cloud, trackingData &td)
Overridable function to handle the particle hitting a patch.
void hitWedgePatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a wedgePatch.
void setCellValues(TrackCloudType &cloud, trackingData &td)
Set cell values.
void correctAfterParallelTransfer(TrackCloudType &, trackingData &)
Make changes following a parallel transfer.
void hitWallPatch(TrackCloudType &cloud, trackingData &td)
Overridable function to handle the particle hitting a wallPatch.
bool hitNonConformalCyclicPatch(const vector &displacement, const scalar fraction, const label patchi, TrackCloudType &cloud, trackingData &td)
Overridable function to handle the particle hitting an.
virtual void transformProperties(const transformer &)
Transform the physical properties of the particle.
void hitBasicPatch(TrackCloudType &cloud, trackingData &td)
Overridable function to handle the particle hitting a basic.
void calc(TrackCloudType &cloud, trackingData &td, const scalar dt)
Update parcel properties over the time interval.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:60
A cloud is a collection of lagrangian particles.
Definition: cloud.H:55
friend dimensionSet transform(const dimensionSet &)
Helper container for force Su and Sp terms.
Definition: forceSuSp.H:64
const vector & Su() const
Return const access to the explicit contribution [kg m/s^2].
Definition: forceSuSpI.H:56
scalar Sp() const
Return const access to the implicit coefficient [kg/s].
Definition: forceSuSpI.H:62
const word & name() const
Return name.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
Definition: polyMesh.C:1006
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:70
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
Definition: tetIndices.H:82
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space.
Definition: transformer.H:84
scalar maxCo
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
label patchi
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar mu
Atomic mass unit.
tmp< VolField< Type > > Su(const VolField< Type > &su, const VolField< Type > &vf)
Definition: fvcSup.C:44
void constrainDirection(const polyMesh &mesh, const Vector< label > &dirs, vector &d)
Set the constrained components of directions/velocity to zero.
Definition: meshTools.C:671
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static const zero Zero
Definition: zero.H:97
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:504
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
error FatalError
static const char nl
Definition: Ostream.H:266
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97
labelList f(nPoints)
volScalarField & p