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>
34 Foam::label Foam::MomentumParcel<ParcelType>::maxTrackAttempts = 1;
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();
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 template<class ParcelType>
255 (
257  const polyMesh& mesh
258 )
259 :
260  ParcelType(p, mesh),
261  moving_(p.moving_),
262  typeId_(p.typeId_),
263  nParticle_(p.nParticle_),
264  d_(p.d_),
265  dTarget_(p.dTarget_),
266  U_(p.U_),
267  rho_(p.rho_),
268  age_(p.age_),
269  tTurb_(p.tTurb_),
270  UTurb_(p.UTurb_)
271 {}
272 
273 
274 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
275 
276 template<class ParcelType>
277 template<class TrackCloudType>
279 (
280  TrackCloudType& cloud,
281  trackingData& td,
282  const scalar trackTime
283 )
284 {
285  typename TrackCloudType::parcelType& p =
286  static_cast<typename TrackCloudType::parcelType&>(*this);
287  typename TrackCloudType::parcelType::trackingData& ttd =
288  static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
289 
290  ttd.keepParticle = true;
291  td.sendToProc = -1;
292 
293  const scalarField& cellLengthScale = cloud.cellLengthScale();
294  const scalar maxCo = cloud.solution().maxCo();
295 
296  while (ttd.keepParticle && ttd.sendToProc == -1 && p.stepFraction() < 1)
297  {
298  // Cache the current position, cell and step-fraction
299  const point start = p.position();
300  const scalar sfrac = p.stepFraction();
301 
302  // Total displacement over the time-step
303  const vector s = trackTime*U_;
304 
305  // Cell length scale
306  const scalar l = cellLengthScale[p.cell()];
307 
308  // Deviation from the mesh centre for reduced-D cases
309  const vector d = p.deviationFromMeshCentre();
310 
311  // Fraction of the displacement to track in this loop. This is limited
312  // to ensure that the both the time and distance tracked is less than
313  // maxCo times the total value.
314  scalar f = 1 - p.stepFraction();
315  f = min(f, maxCo);
316  f = min(f, maxCo/min(max(mag(s)/l, rootSmall), rootGreat));
317  if (p.moving())
318  {
319  // Track to the next face
320  p.trackToFace(f*s - d, f);
321  }
322  else
323  {
324  // At present the only thing that sets moving_ to false is a stick
325  // wall interaction. We want the position of the particle to remain
326  // the same relative to the face that it is on. The local
327  // coordinates therefore do not change. We still advance in time and
328  // perform the relevant interactions with the fixed particle.
329  p.stepFraction() += f;
330  }
331 
332  const scalar dt = (p.stepFraction() - sfrac)*trackTime;
333 
334  // Avoid problems with extremely small timesteps
335  if (dt > rootVSmall)
336  {
337  // Update cell based properties
338  p.setCellValues(cloud, ttd);
339 
340  p.calcDispersion(cloud, ttd, dt);
341 
342  if (cloud.solution().cellValueSourceCorrection())
343  {
344  p.cellValueSourceCorrection(cloud, ttd, dt);
345  }
346 
347  p.calc(cloud, ttd, dt);
348  }
349 
350  p.age() += dt;
351 
352  if (p.moving() && p.onFace())
353  {
354  cloud.functions().postFace(p, ttd.keepParticle);
355  }
356 
357  cloud.functions().postMove(p, dt, start, ttd.keepParticle);
358 
359  if (p.moving() && p.onFace() && ttd.keepParticle)
360  {
361  p.hitFace(f*s - d, f, cloud, ttd);
362  }
363  }
364 
365  return ttd.keepParticle;
366 }
367 
368 
369 template<class ParcelType>
370 template<class TrackCloudType>
372 (
373  TrackCloudType& cloud,
374  trackingData& td
375 )
376 {
377  typename TrackCloudType::parcelType& p =
378  static_cast<typename TrackCloudType::parcelType&>(*this);
379 
380  const polyPatch& pp = p.mesh().boundaryMesh()[p.patch()];
381 
382  // Invoke post-processing model
383  cloud.functions().postPatch(p, pp, td.keepParticle);
384 
385  // Invoke surface film model
386  if (cloud.surfaceFilm().transferParcel(p, pp, td.keepParticle))
387  {
388  // All interactions done
389  return true;
390  }
391  else if (pp.coupled())
392  {
393  // Don't apply the patchInteraction models to coupled boundaries
394  return false;
395  }
396  else
397  {
398  // Invoke patch interaction model
399  return cloud.patchInteraction().correct(p, pp, td.keepParticle);
400  }
401 }
402 
403 
404 template<class ParcelType>
405 template<class TrackCloudType>
407 (
408  TrackCloudType&,
409  trackingData&
410 )
411 {
412  const polyPatch& pp = this->mesh().boundaryMesh()[this->patch()];
413 
415  << "Particle " << this->origId() << " hit " << pp.type() << " patch "
416  << pp.name() << " at " << this->position()
417  << " but no interaction model was specified for this patch"
418  << exit(FatalError);
419 }
420 
421 
422 template<class ParcelType>
424 (
425  const transformer& transform
426 )
427 {
428  ParcelType::transformProperties(transform);
429  U_ = transform.transform(U_);
430 }
431 
432 
433 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
434 
435 #include "MomentumParcelIO.C"
436 
437 // ************************************************************************* //
bool moving_
Moving flag - tracking stopped when moving = false.
scalar d_
Diameter [m].
void hitWallPatch(TrackCloudType &cloud, trackingData &td)
Overridable function to handle the particle hitting a wallPatch.
const interpolation< vector > & UInterp() const
Return const access to the interpolator for continuous.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
const word & name() const
Return name.
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space...
Definition: transformer.H:83
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
scalar dTarget_
Target diameter [m].
void setCellValues(TrackCloudType &cloud, trackingData &td)
Set cell values.
const vector & Su() const
Return const access to the explicit contribution [kg m/s^2].
Definition: forceSuSpI.H:56
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
bool hitPatch(TrackCloudType &cloud, trackingData &td)
Overridable function to handle the particle hitting a patch.
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:297
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
Definition: polyMesh.C:911
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label typeId_
Parcel type id.
fvMesh & mesh
void cellValueSourceCorrection(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct cell values using latest transfer information.
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.
Helper container for force Su and Sp terms.
Definition: forceSuSp.H:61
scalar rho_
Density [kg/m^3].
scalar nParticle_
Number of particles in Parcel.
void calc(TrackCloudType &cloud, trackingData &td, const scalar dt)
Update parcel properties over the time interval.
virtual Type interpolate(const vector &position, const label celli, const label facei=-1) const =0
Interpolate field to the given point in the given cell.
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:316
const interpolation< scalar > & rhoInterp() const
Return const access to the interpolator for continuous.
scalar age_
Age [s].
scalar tTurb_
Time spent in turbulent eddy [s].
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
static const zero Zero
Definition: zero.H:97
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
scalar muc() const
Return the continuous phase viscosity.
vector U_
Velocity of Parcel [m/s].
Momentum parcel class with rotational motion (as spherical particles only) and one/two-way coupling w...
scalar Sp() const
Return const access to the implicit coefficient [kg/s].
Definition: forceSuSpI.H:62
scalar maxCo
static const char nl
Definition: Ostream.H:260
labelList f(nPoints)
virtual void transformProperties(const transformer &)
Transform the physical properties of the particle.
void calcDispersion(TrackCloudType &cloud, trackingData &td, const scalar dt)
Apply dispersion to the carrier phase velocity and update.
MomentumParcel(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from mesh, coordinates and topology.
scalar rhoc() const
Return the continuous phase density.
#define WarningInFunction
Report a warning using Foam::Warning.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
const interpolation< scalar > & muInterp() const
Return const access to the interpolator for continuous.
dimensioned< scalar > mag(const dimensioned< Type > &)
Type transform(const Type &) const
Transform the given type.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
const tmp< volScalarField::Internal > & Su
Definition: alphaSuSp.H:6
bool move(TrackCloudType &cloud, trackingData &td, const scalar trackTime)
Move the parcel.
void constrainDirection(const polyMesh &mesh, const Vector< label > &dirs, vector &d)
Set the constrained components of directions/velocity to zero.
Definition: meshTools.C:671
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
const vector & Uc() const
Return the continuous phase velocity.
vector UTurb_
Turbulent velocity fluctuation [m/s].
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97