FreeFOAM The Cross-Platform CFD Toolkit
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) 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 "KinematicParcel.H"
27 
28 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
29 
30 template<class ParcelType>
31 template<class TrackData>
33 (
34  TrackData& td,
35  const scalar dt,
36  const label cellI
37 )
38 {
39  rhoc_ = td.rhoInterp().interpolate(this->position(), cellI);
40  if (rhoc_ < td.constProps().rhoMin())
41  {
42  WarningIn
43  (
44  "void Foam::KinematicParcel<ParcelType>::setCellValues"
45  "("
46  "TrackData&, "
47  "const scalar, "
48  "const label"
49  ")"
50  ) << "Limiting observed density in cell " << cellI << " to "
51  << td.constProps().rhoMin() << nl << endl;
52 
53  rhoc_ = td.constProps().rhoMin();
54  }
55 
56  Uc_ = td.UInterp().interpolate(this->position(), cellI);
57 
58  muc_ = td.muInterp().interpolate(this->position(), cellI);
59 
60  // Apply dispersion components to carrier phase velocity
61  Uc_ = td.cloud().dispersion().update
62  (
63  dt,
64  cellI,
65  U_,
66  Uc_,
67  UTurb_,
68  tTurb_
69  );
70 }
71 
72 
73 template<class ParcelType>
74 template<class TrackData>
76 (
77  TrackData& td,
78  const scalar dt,
79  const label cellI
80 )
81 {
82  Uc_ += td.cloud().UTrans()[cellI]/massCell(cellI);
83 }
84 
85 
86 template<class ParcelType>
87 template<class TrackData>
89 (
90  TrackData& td,
91  const scalar dt,
92  const label cellI
93 )
94 {
95  // Define local properties at beginning of time step
96  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
97  const scalar np0 = nParticle_;
98  const scalar d0 = d_;
99  const vector U0 = U_;
100  const scalar rho0 = rho_;
101  const scalar mass0 = mass();
102 
103  // Reynolds number
104  const scalar Re = this->Re(U0, d0, rhoc_, muc_);
105 
106 
107  // Sources
108  //~~~~~~~~
109 
110  // Explicit momentum source for particle
111  vector Su = vector::zero;
112 
113  // Momentum transfer from the particle to the carrier phase
114  vector dUTrans = vector::zero;
115 
116 
117  // Motion
118  // ~~~~~~
119 
120  // Calculate new particle velocity
121  vector U1 =
122  calcVelocity(td, dt, cellI, Re, muc_, d0, U0, rho0, mass0, Su, dUTrans);
123 
124 
125  // Accumulate carrier phase source terms
126  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
127  if (td.cloud().coupled())
128  {
129  // Update momentum transfer
130  td.cloud().UTrans()[cellI] += np0*dUTrans;
131  }
132 
133 
134  // Set new particle properties
135  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
136  U_ = U1;
137 }
138 
139 
140 template<class ParcelType>
141 template<class TrackData>
143 (
144  TrackData& td,
145  const scalar dt,
146  const label cellI,
147  const scalar Re,
148  const scalar mu,
149  const scalar d,
150  const vector& U,
151  const scalar rho,
152  const scalar mass,
153  const vector& Su,
154  vector& dUTrans
155 ) const
156 {
157  const polyMesh& mesh = this->cloud().pMesh();
158 
159  // Momentum transfer coefficient
160  const scalar utc = td.cloud().drag().utc(Re, d, mu) + ROOTVSMALL;
161 
162  // Momentum source due to particle forces
163  const vector FCoupled =
164  mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U);
165  const vector FNonCoupled =
166  mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U);
167 
168 
169  // New particle velocity
170  //~~~~~~~~~~~~~~~~~~~~~~
171 
172  // Update velocity - treat as 3-D
173  const scalar As = this->areaS(d);
174  const vector ap = Uc_ + (FCoupled + FNonCoupled + Su)/(utc*As);
175  const scalar bp = 6.0*utc/(rho*d);
176 
178  td.cloud().UIntegrator().integrate(U, dt, ap, bp);
179 
180  vector Unew = Ures.value();
181 
182  dUTrans += dt*(utc*As*(Ures.average() - Uc_) - FCoupled);
183 
184  // Apply correction to velocity and dUTrans for reduced-D cases
185  meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
186  meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);
187 
188  return Unew;
189 }
190 
191 
192 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
193 
194 template <class ParcelType>
196 (
198 )
199 :
200  Particle<ParcelType>(p),
201  typeId_(p.typeId_),
202  nParticle_(p.nParticle_),
203  d_(p.d_),
204  U_(p.U_),
205  rho_(p.rho_),
206  tTurb_(p.tTurb_),
207  UTurb_(p.UTurb_),
208  rhoc_(p.rhoc_),
209  Uc_(p.Uc_),
210  muc_(p.muc_)
211 {}
212 
213 
214 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
215 
216 template<class ParcelType>
217 template<class TrackData>
219 {
220  ParcelType& p = static_cast<ParcelType&>(*this);
221 
222  td.switchProcessor = false;
223  td.keepParticle = true;
224 
225  const polyMesh& mesh = td.cloud().pMesh();
226  const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
227 
228  const scalar deltaT = mesh.time().deltaTValue();
229  scalar tEnd = (1.0 - p.stepFraction())*deltaT;
230  const scalar dtMax = tEnd;
231 
232  while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
233  {
234  // Apply correction to position for reduced-D cases
235  meshTools::constrainToMeshCentre(mesh, p.position());
236 
237  // Set the Lagrangian time-step
238  scalar dt = min(dtMax, tEnd);
239 
240  // Remember which cell the Parcel is in since this will change if a
241  // face is hit
242  label cellI = p.cell();
243 
244  if (p.active())
245  {
246  dt *= p.trackToFace(p.position() + dt*U_, td);
247  }
248 
249  tEnd -= dt;
250  p.stepFraction() = 1.0 - tEnd/deltaT;
251 
252  // Avoid problems with extremely small timesteps
253  if (dt > ROOTVSMALL)
254  {
255  // Update cell based properties
256  p.setCellValues(td, dt, cellI);
257 
258  if (td.cloud().cellValueSourceCorrection())
259  {
260  p.cellValueSourceCorrection(td, dt, cellI);
261  }
262 
263  p.calc(td, dt, cellI);
264  }
265 
266  if (p.onBoundary() && td.keepParticle)
267  {
268  if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
269  {
270  td.switchProcessor = true;
271  }
272  }
273  }
274 
275  return td.keepParticle;
276 }
277 
278 
279 template<class ParcelType>
280 template<class TrackData>
282 (
283  const polyPatch& pp,
284  TrackData& td,
285  const label patchI
286 )
287 {
288  ParcelType& p = static_cast<ParcelType&>(*this);
289  td.cloud().postProcessing().postPatch(p, patchI);
290 
291  return td.cloud().patchInteraction().correct
292  (
293  pp,
294  this->face(),
295  td.keepParticle,
296  active_,
297  U_
298  );
299 }
300 
301 
302 template<class ParcelType>
304 (
305  const polyPatch& pp,
306  int& td,
307  const label patchI
308 )
309 {
310  return false;
311 }
312 
313 
314 template<class ParcelType>
315 template<class TrackData>
317 (
318  const processorPolyPatch&,
319  TrackData& td
320 )
321 {
322  td.switchProcessor = true;
323 }
324 
325 
326 template<class ParcelType>
328 (
329  const processorPolyPatch&,
330  int&
331 )
332 {}
333 
334 
335 template<class ParcelType>
336 template<class TrackData>
338 (
339  const wallPolyPatch& wpp,
340  TrackData& td
341 )
342 {
343  // Wall interactions handled by generic hitPatch function
344 }
345 
346 
347 template<class ParcelType>
349 (
350  const wallPolyPatch&,
351  int&
352 )
353 {}
354 
355 
356 template<class ParcelType>
357 template<class TrackData>
359 (
360  const polyPatch&,
361  TrackData& td
362 )
363 {
364  td.keepParticle = false;
365 }
366 
367 
368 template<class ParcelType>
370 {}
371 
372 
373 template<class ParcelType>
375 {
377  U_ = transform(T, U_);
378 }
379 
380 
381 template<class ParcelType>
383 (
384  const vector& separation
385 )
386 {
388 }
389 
390 
391 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
392 
394 
395 // ************************ vim: set sw=4 sts=4 et: ************************ //