ORourkeCollision.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-2016 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 "ORourkeCollision.H"
27 #include "SLGThermo.H"
28 #include "CompactListList.H"
29 #include "mathematicalConstants.H"
30 
31 using namespace Foam::constant::mathematical;
32 
33 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
34 
35 template<class CloudType>
37 {
38  // Create the occupancy list for the cells
39  labelList occupancy(this->owner().mesh().nCells(), 0);
40  forAllIter(typename CloudType, this->owner(), iter)
41  {
42  occupancy[iter().cell()]++;
43  }
44 
45  // Initialize the sizes of the lists of parcels in each cell
46  CompactListList<parcelType*> pInCell(occupancy);
47 
48  // Reset the occupancy to use as a counter
49  occupancy = 0;
50 
51  // Set the parcel pointer lists for each cell
52  forAllIter(typename CloudType, this->owner(), iter)
53  {
54  pInCell(iter().cell(), occupancy[iter().cell()]++) = &iter();
55  }
56 
57  for (label celli=0; celli<this->owner().mesh().nCells(); celli++)
58  {
59  UList<parcelType*> pInCelli(pInCell[celli]);
60 
61  if (pInCelli.size() >= 2)
62  {
63  forAll(pInCelli, i)
64  {
65  for (label j=i+1; j<pInCelli.size(); j++)
66  {
67  parcelType& p1 = *pInCelli[i];
68  parcelType& p2 = *pInCelli[j];
69 
70  scalar m1 = p1.nParticle()*p1.mass();
71  scalar m2 = p2.nParticle()*p2.mass();
72 
73  bool massChanged = collideParcels(dt, p1, p2, m1, m2);
74 
75  if (massChanged)
76  {
77  if (m1 > ROOTVSMALL)
78  {
79  const scalarField X(liquids_.X(p1.Y()));
80  p1.rho() = liquids_.rho(p1.pc(), p1.T(), X);
81  p1.Cp() = liquids_.Cp(p1.pc(), p1.T(), X);
82  p1.sigma() = liquids_.sigma(p1.pc(), p1.T(), X);
83  p1.mu() = liquids_.mu(p1.pc(), p1.T(), X);
84  p1.d() = cbrt(6.0*m1/(p1.nParticle()*p1.rho()*pi));
85  }
86 
87  if (m2 > ROOTVSMALL)
88  {
89  const scalarField X(liquids_.X(p2.Y()));
90  p2.rho() = liquids_.rho(p2.pc(), p2.T(), X);
91  p2.Cp() = liquids_.Cp(p2.pc(), p2.T(), X);
92  p2.sigma() = liquids_.sigma(p2.pc(), p2.T(), X);
93  p2.mu() = liquids_.mu(p2.pc(), p2.T(), X);
94  p2.d() = cbrt(6.0*m2/(p2.nParticle()*p2.rho()*pi));
95  }
96  }
97  }
98  }
99  }
100  }
101 
102  // Remove coalesced parcels that fall below minimum mass threshold
103  forAllIter(typename CloudType, this->owner(), iter)
104  {
105  parcelType& p = iter();
106  scalar mass = p.nParticle()*p.mass();
107 
108  if (mass < this->owner().constProps().minParcelMass())
109  {
110  this->owner().deleteParticle(p);
111  }
112  }
113 }
114 
115 
116 template<class CloudType>
118 (
119  const scalar dt,
120  parcelType& p1,
121  parcelType& p2,
122  scalar& m1,
123  scalar& m2
124 )
125 {
126  // Return if parcel masses are ~0
127  if ((m1 < ROOTVSMALL) || (m2 < ROOTVSMALL))
128  {
129  return false;
130  }
131 
132  const scalar Vc = this->owner().mesh().V()[p1.cell()];
133  const scalar d1 = p1.d();
134  const scalar d2 = p2.d();
135 
136  scalar magUrel = mag(p1.U() - p2.U());
137  scalar sumD = d1 + d2;
138  scalar nu0 = 0.25*constant::mathematical::pi*sqr(sumD)*magUrel*dt/Vc;
139  scalar nMin = min(p1.nParticle(), p2.nParticle());
140  scalar nu = nMin*nu0;
141  scalar collProb = exp(-nu);
142  scalar xx = this->owner().rndGen().template sample01<scalar>();
143 
144  // Collision occurs
145  if (xx > collProb)
146  {
147  if (d1 > d2)
148  {
149  return collideSorted(dt, p1, p2, m1, m2);
150  }
151  else
152  {
153  return collideSorted(dt, p2, p1, m2, m1);
154  }
155  }
156  else
157  {
158  return false;
159  }
160 }
161 
162 
163 template<class CloudType>
165 (
166  const scalar dt,
167  parcelType& p1,
168  parcelType& p2,
169  scalar& m1,
170  scalar& m2
171 )
172 {
173  const scalar nP1 = p1.nParticle();
174  const scalar nP2 = p2.nParticle();
175 
176  const scalar sigma1 = p1.sigma();
177  const scalar sigma2 = p2.sigma();
178 
179  const scalar d1 = p1.d();
180  const scalar d2 = p2.d();
181 
182  const scalar T1 = p1.T();
183  const scalar T2 = p2.T();
184 
185  const scalar rho1 = p1.rho();
186  const scalar rho2 = p2.rho();
187 
188  const vector& U1 = p1.U();
189  const vector& U2 = p2.U();
190 
191 
192  vector URel = U1 - U2;
193  scalar magURel = mag(URel);
194 
195  scalar mTot = m1 + m2;
196 
197  scalar gamma = d1/max(ROOTVSMALL, d2);
198  scalar f = pow3(gamma) + 2.7*gamma - 2.4*sqr(gamma);
199 
200  // Mass-averaged temperature
201  scalar Tave = (T1*m1 + T2*m2)/mTot;
202 
203  // Interpolate to find average surface tension
204  scalar sigmaAve = sigma1;
205  if (mag(T2 - T1) > SMALL)
206  {
207  sigmaAve += (sigma2 - sigma1)*(Tave - T1)/(T2 - T1);
208  }
209 
210  scalar Vtot = m1/rho1 + m2/rho2;
211  scalar rhoAve = mTot/Vtot;
212 
213  scalar dAve = sqrt(d1*d2);
214  scalar WeColl = 0.5*rhoAve*sqr(magURel)*dAve/max(ROOTVSMALL, sigmaAve);
215 
216  scalar coalesceProb = min(1.0, 2.4*f/max(ROOTVSMALL, WeColl));
217 
218  scalar prob = this->owner().rndGen().template sample01<scalar>();
219 
220  // Coalescence
221  if (coalescence_ && prob < coalesceProb)
222  {
223  // Number of the droplets that coalesce
224  scalar nProb = prob*nP2/nP1;
225 
226  // Conservation of mass, momentum and energy
227  scalar m1Org = m1;
228  scalar m2Org = m2;
229 
230  scalar dm = nP1*nProb*m2/scalar(nP2);
231 
232  m1 += dm;
233  m2 -= dm;
234 
235  p1.T() = (Tave*mTot - m2*T2)/m1;
236 
237  p1.U() = (m1*U1 + (1.0 - m2/m2Org)*m2*U2)/m1;
238 
239  p1.Y() = (m1Org*p1.Y() + dm*p2.Y())/m1;
240 
241  p2.nParticle() = m2/(rho2*p2.volume());
242 
243  return true;
244  }
245  // Grazing collision (no coalescence)
246  else
247  {
248  scalar gf = sqrt(prob) - sqrt(coalesceProb);
249  scalar denom = 1.0 - sqrt(coalesceProb);
250  if (denom < 1.0e-5)
251  {
252  denom = 1.0;
253  }
254  gf /= denom;
255 
256  // If gf negative, this means that coalescence is turned off
257  // and these parcels should have coalesced
258  gf = max(0.0, gf);
259 
260  // gf -> 1 => v1p -> U1 ...
261  // gf -> 0 => v1p -> momentum/mTot
262 
263  vector mr = m1*U1 + m2*U2;
264  vector v1p = (mr + m2*gf*URel)/mTot;
265  vector v2p = (mr - m1*gf*URel)/mTot;
266 
267  if (nP1 < nP2)
268  {
269  p1.U() = v1p;
270  p2.U() = (nP1*v2p + (nP2 - nP1)*U2)/nP2;
271  }
272  else
273  {
274  p1.U() = (nP2*v1p + (nP1 - nP2)*U1)/nP1;
275  p2.U() = v2p;
276  }
277 
278  return false;
279  }
280 }
281 
282 
283 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
284 
285 template<class CloudType>
287 (
288  const dictionary& dict,
289  CloudType& owner,
290  const word& modelName
291 )
292 :
293  StochasticCollisionModel<CloudType>(dict, owner, modelName),
294  liquids_
295  (
296  owner.db().template lookupObject<SLGThermo>("SLGThermo").liquids()
297  ),
298  coalescence_(this->coeffDict().lookup("coalescence"))
299 {}
300 
301 
302 template<class CloudType>
304 (
306 )
307 :
309  liquids_(cm.liquids_),
310  coalescence_(cm.coalescence_)
311 {}
312 
313 
314 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
315 
316 template<class CloudType>
318 {}
319 
320 
321 // ************************************************************************* //
ORourkeCollision(const dictionary &dict, CloudType &cloud, const word &modelName=typeName)
Construct from dictionary.
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
const scalar mTot
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:453
dimensionedScalar sqrt(const dimensionedScalar &ds)
virtual ~ORourkeCollision()
Destructor.
rho1
Definition: pEqn.H:114
dynamicFvMesh & mesh
dimensionedScalar exp(const dimensionedScalar &ds)
mathematical constants.
rho2
Definition: pEqn.H:115
A class for handling words, derived from string.
Definition: word.H:59
dimensionedScalar cbrt(const dimensionedScalar &ds)
virtual bool collideParcels(const scalar dt, parcelType &p1, parcelType &p2, scalar &m1, scalar &m2)
Collide parcels and return true if mass has changed.
volVectorField & U1
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
A packed storage unstructured matrix of objects of type <T> using an offset table for access...
virtual void collide(const scalar dt)
Main collision routine.
CloudType::parcelType parcelType
Convenience typedef to the cloud&#39;s parcel type.
virtual bool collideSorted(const scalar dt, parcelType &p1, parcelType &p2, scalar &m1, scalar &m2)
labelList f(nPoints)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const liquidMixtureProperties & liquids_
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
Collision model by P.J. O&#39;Rourke.
dimensionedScalar pow3(const dimensionedScalar &ds)
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
Templated stochastic collision model class.
volVectorField & U2
dimensioned< scalar > mag(const dimensioned< Type > &)
volScalarField & p
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
volScalarField & nu
Switch coalescence_
Coalescence activation switch.