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