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