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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "ORourkeCollision.H"
30 #include "SLGThermo.H"
31 #include "CompactListList.H"
32 #include "mathematicalConstants.H"
33 
34 using namespace Foam::constant::mathematical;
35 
36 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
37 
38 template<class CloudType>
40 (
41  typename CloudType::parcelType::trackingData& td,
42  const scalar dt
43 )
44 {
45  // Create the occupancy list for the cells
46  labelList occupancy(this->owner().mesh().nCells(), Zero);
47  for (const parcelType& p : this->owner())
48  {
49  occupancy[p.cell()]++;
50  }
51 
52  // Initialize the sizes of the lists of parcels in each cell
53  CompactListList<parcelType*> pInCell(occupancy);
54 
55  // Reset the occupancy to use as a counter
56  occupancy = 0;
57 
58  // Set the parcel pointer lists for each cell
59  for (parcelType& p : this->owner())
60  {
61  pInCell(p.cell(), occupancy[p.cell()]++) = &p;
62  }
63 
64  for (label celli=0; celli<this->owner().mesh().nCells(); celli++)
65  {
66  UList<parcelType*> pInCelli(pInCell[celli]);
67 
68  if (pInCelli.size() >= 2)
69  {
70  forAll(pInCelli, i)
71  {
72  for (label j=i+1; j<pInCelli.size(); j++)
73  {
74  parcelType& p1 = *pInCelli[i];
75  parcelType& p2 = *pInCelli[j];
76 
77  scalar m1 = p1.nParticle()*p1.mass();
78  scalar m2 = p2.nParticle()*p2.mass();
79 
80  bool massChanged = collideParcels(dt, p1, p2, m1, m2);
81 
82  if (massChanged)
83  {
84  if (m1 > ROOTVSMALL)
85  {
86  const scalarField X(liquids_.X(p1.Y()));
87  p1.setCellValues(this->owner(), td);
88  p1.rho() = liquids_.rho(td.pc(), p1.T(), X);
89  p1.Cp() = liquids_.Cp(td.pc(), p1.T(), X);
90  p1.sigma() = liquids_.sigma(td.pc(), p1.T(), X);
91  p1.mu() = liquids_.mu(td.pc(), p1.T(), X);
92  p1.d() = cbrt(6.0*m1/(p1.nParticle()*p1.rho()*pi));
93  }
94 
95  if (m2 > ROOTVSMALL)
96  {
97  const scalarField X(liquids_.X(p2.Y()));
98  p2.setCellValues(this->owner(), td);
99  p2.rho() = liquids_.rho(td.pc(), p2.T(), X);
100  p2.Cp() = liquids_.Cp(td.pc(), p2.T(), X);
101  p2.sigma() = liquids_.sigma(td.pc(), p2.T(), X);
102  p2.mu() = liquids_.mu(td.pc(), p2.T(), X);
103  p2.d() = cbrt(6.0*m2/(p2.nParticle()*p2.rho()*pi));
104  }
105  }
106  }
107  }
108  }
109  }
110 
111  // Remove coalesced parcels that fall below minimum mass threshold
112  for (parcelType& p : this->owner())
113  {
114  scalar mass = p.nParticle()*p.mass();
115 
116  if (mass < this->owner().constProps().minParcelMass())
117  {
118  this->owner().deleteParticle(p);
119  }
120  }
121 }
122 
123 
124 template<class CloudType>
126 (
127  const scalar dt,
128  parcelType& p1,
129  parcelType& p2,
130  scalar& m1,
131  scalar& m2
132 )
133 {
134  // Return if parcel masses are ~0
135  if ((m1 < ROOTVSMALL) || (m2 < ROOTVSMALL))
136  {
137  return false;
138  }
139 
140  const scalar Vc = this->owner().mesh().V()[p1.cell()];
141  const scalar d1 = p1.d();
142  const scalar d2 = p2.d();
143 
144  scalar magUrel = mag(p1.U() - p2.U());
145  scalar sumD = d1 + d2;
146  scalar nu0 = 0.25*constant::mathematical::pi*sqr(sumD)*magUrel*dt/Vc;
147  scalar nMin = min(p1.nParticle(), p2.nParticle());
148  scalar nu = nMin*nu0;
149  scalar collProb = exp(-nu);
150  scalar xx = this->owner().rndGen().template sample01<scalar>();
151 
152  // Collision occurs
153  if (xx > collProb)
154  {
155  if (d1 > d2)
156  {
157  return collideSorted(dt, p1, p2, m1, m2);
158  }
159  else
160  {
161  return collideSorted(dt, p2, p1, m2, m1);
162  }
163  }
164 
165  return false;
166 }
167 
168 
169 template<class CloudType>
171 (
172  const scalar dt,
173  parcelType& p1,
174  parcelType& p2,
175  scalar& m1,
176  scalar& m2
177 )
178 {
179  const scalar nP1 = p1.nParticle();
180  const scalar nP2 = p2.nParticle();
181 
182  const scalar sigma1 = p1.sigma();
183  const scalar sigma2 = p2.sigma();
184 
185  const scalar d1 = p1.d();
186  const scalar d2 = p2.d();
187 
188  const scalar T1 = p1.T();
189  const scalar T2 = p2.T();
190 
191  const scalar rho1 = p1.rho();
192  const scalar rho2 = p2.rho();
193 
194  const vector& U1 = p1.U();
195  const vector& U2 = p2.U();
196 
197 
198  vector URel = U1 - U2;
199  scalar magURel = mag(URel);
200 
201  scalar mTot = m1 + m2;
202 
203  scalar gamma = d1/max(ROOTVSMALL, d2);
204  scalar f = pow3(gamma) + 2.7*gamma - 2.4*sqr(gamma);
205 
206  // Mass-averaged temperature
207  scalar Tave = (T1*m1 + T2*m2)/mTot;
208 
209  // Interpolate to find average surface tension
210  scalar sigmaAve = sigma1;
211  if (mag(T2 - T1) > SMALL)
212  {
213  sigmaAve += (sigma2 - sigma1)*(Tave - T1)/(T2 - T1);
214  }
215 
216  scalar Vtot = m1/rho1 + m2/rho2;
217  scalar rhoAve = mTot/Vtot;
218 
219  scalar dAve = sqrt(d1*d2);
220  scalar WeColl = 0.5*rhoAve*sqr(magURel)*dAve/max(ROOTVSMALL, sigmaAve);
221 
222  scalar coalesceProb = min(1.0, 2.4*f/max(ROOTVSMALL, WeColl));
223 
224  scalar prob = this->owner().rndGen().template sample01<scalar>();
225 
226  // Coalescence
227  if (coalescence_ && prob < coalesceProb)
228  {
229  // Number of the droplets that coalesce
230  scalar nProb = prob*nP2/nP1;
231 
232  // Conservation of mass, momentum and energy
233  scalar m1Org = m1;
234  scalar m2Org = m2;
235 
236  scalar dm = nP1*nProb*m2/scalar(nP2);
237 
238  m1 += dm;
239  m2 -= dm;
240 
241  p1.T() = (Tave*mTot - m2*T2)/m1;
242 
243  p1.U() = (m1*U1 + (1.0 - m2/m2Org)*m2*U2)/m1;
244 
245  p1.Y() = (m1Org*p1.Y() + dm*p2.Y())/m1;
246 
247  p2.nParticle() = m2/(rho2*p2.volume());
248 
249  return true;
250  }
251  // Grazing collision (no coalescence)
252  else
253  {
254  scalar gf = sqrt(prob) - sqrt(coalesceProb);
255  scalar denom = 1.0 - sqrt(coalesceProb);
256  if (denom < 1.0e-5)
257  {
258  denom = 1.0;
259  }
260  gf /= denom;
261 
262  // If gf negative, this means that coalescence is turned off
263  // and these parcels should have coalesced
264  gf = max(0.0, gf);
265 
266  // gf -> 1 => v1p -> U1 ...
267  // gf -> 0 => v1p -> momentum/mTot
268 
269  vector mr = m1*U1 + m2*U2;
270  vector v1p = (mr + m2*gf*URel)/mTot;
271  vector v2p = (mr - m1*gf*URel)/mTot;
272 
273  if (nP1 < nP2)
274  {
275  p1.U() = v1p;
276  p2.U() = (nP1*v2p + (nP2 - nP1)*U2)/nP2;
277  }
278  else
279  {
280  p1.U() = (nP2*v1p + (nP1 - nP2)*U1)/nP1;
281  p2.U() = v2p;
282  }
283 
284  return false;
285  }
286 }
287 
288 
289 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
290 
291 template<class CloudType>
293 (
294  const dictionary& dict,
295  CloudType& owner,
296  const word& modelName
297 )
298 :
299  StochasticCollisionModel<CloudType>(dict, owner, modelName),
300  liquids_
301  (
302  owner.db().template lookupObject<SLGThermo>("SLGThermo").liquids()
303  ),
304  coalescence_
305  (
306  this->coeffDict().template get<bool>("coalescence")
307  )
308 {}
309 
310 
311 template<class CloudType>
313 (
315 )
316 :
318  liquids_(cm.liquids_),
319  coalescence_(cm.coalescence_)
320 {}
321 
322 
323 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
324 
325 template<class CloudType>
327 {}
328 
329 
330 // ************************************************************************* //
U1
volVectorField & U1
Definition: setRegionFluidFields.H:11
mathematicalConstants.H
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::CompactListList
A packed storage unstructured matrix of objects of type <T> using an offset table for access.
Definition: CompactListList.H:63
Foam::constant::mathematical::e
constexpr scalar e(M_E)
Foam::exp
dimensionedScalar exp(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:261
Foam::ORourkeCollision::~ORourkeCollision
virtual ~ORourkeCollision()
Destructor.
Definition: ORourkeCollision.C:326
SLGThermo.H
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::ORourkeCollision
Collision model by P.J. O'Rourke.
Definition: ORourkeCollision.H:49
Foam::DSMCCloud::rndGen
Random & rndGen()
Return reference to the random object.
Definition: DSMCCloudI.H:124
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
nu
volScalarField & nu
Definition: readMechanicalProperties.H:176
Foam::primitiveMesh::nCells
label nCells() const noexcept
Number of mesh cells.
Definition: primitiveMeshI.H:96
Foam::StochasticCollisionModel
Templated stochastic collision model class.
Definition: KinematicCloud.H:95
Foam::ORourkeCollision::ORourkeCollision
ORourkeCollision(const dictionary &dict, CloudType &cloud, const word &modelName=typeName)
Construct from dictionary.
Definition: ORourkeCollision.C:293
Foam::Field< scalar >
Foam::pow3
dimensionedScalar pow3(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:89
Foam::DSMCCloud::mesh
const fvMesh & mesh() const
Return reference to the mesh.
Definition: DSMCCloudI.H:44
Foam::Cloud::deleteParticle
void deleteParticle(ParticleType &p)
Remove particle from cloud and delete.
Definition: Cloud.C:111
rho2
volScalarField & rho2
Definition: setRegionFluidFields.H:30
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
rho1
volScalarField & rho1
Definition: setRegionFluidFields.H:27
Foam::DSMCCloud
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:71
U2
volVectorField & U2
Definition: setRegionFluidFields.H:15
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::ORourkeCollision::parcelType
CloudType::parcelType parcelType
Convenience typedef to the cloud's parcel type.
Definition: ORourkeCollision.H:58
Foam::ORourkeCollision::collideSorted
virtual bool collideSorted(const scalar dt, parcelType &p1, parcelType &p2, scalar &m1, scalar &m2)
Definition: ORourkeCollision.C:171
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical
Mathematical constants.
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
f
labelList f(nPoints)
Foam::Vector< scalar >
Foam::List< label >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
gamma
const scalar gamma
Definition: EEqn.H:9
Foam::UList
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:103
CompactListList.H
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Foam::cbrt
dimensionedScalar cbrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:155
mTot
const scalar mTot
Definition: readInitialConditions.H:64
ORourkeCollision.H
Foam::ORourkeCollision::collideParcels
virtual bool collideParcels(const scalar dt, parcelType &p1, parcelType &p2, scalar &m1, scalar &m2)
Collide parcels and return true if mass has changed.
Definition: ORourkeCollision.C:126
Foam::ORourkeCollision::collide
virtual void collide(typename CloudType::parcelType::trackingData &td, const scalar dt)
Main collision routine.
Definition: ORourkeCollision.C:40
Foam::ORourkeCollision::liquids_
const liquidMixtureProperties & liquids_
Definition: ORourkeCollision.H:60
Foam::fvMesh::V
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
Definition: fvMeshGeometry.C:179
Foam::ORourkeCollision::coalescence_
bool coalescence_
Coalescence activated?
Definition: ORourkeCollision.H:63