LarsenBorgnakkeVariableHardSphere.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-2015 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 
27 #include "constants.H"
28 
29 using namespace Foam::constant::mathematical;
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 template<class CloudType>
35 (
36  scalar ChiA,
37  scalar ChiB
38 )
39 {
40  CloudType& cloud(this->owner());
41 
42  Random& rndGen(cloud.rndGen());
43 
44  scalar ChiAMinusOne = ChiA - 1;
45  scalar ChiBMinusOne = ChiB - 1;
46 
47  if (ChiAMinusOne < SMALL && ChiBMinusOne < SMALL)
48  {
49  return rndGen.scalar01();
50  }
51 
52  scalar energyRatio;
53  scalar P;
54 
55  do
56  {
57  P = 0;
58 
59  energyRatio = rndGen.scalar01();
60 
61  if (ChiAMinusOne < SMALL)
62  {
63  P = 1.0 - pow(energyRatio, ChiB);
64  }
65  else if (ChiBMinusOne < SMALL)
66  {
67  P = 1.0 - pow(energyRatio, ChiA);
68  }
69  else
70  {
71  P =
72  pow
73  (
74  (ChiAMinusOne + ChiBMinusOne)*energyRatio/ChiAMinusOne,
75  ChiAMinusOne
76  )
77  *pow
78  (
79  (ChiAMinusOne + ChiBMinusOne)*(1 - energyRatio)
80  /ChiBMinusOne,
81  ChiBMinusOne
82  );
83  }
84  } while (P < rndGen.scalar01());
85 
86  return energyRatio;
87 }
88 
89 
90 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
91 
92 template<class CloudType>
95 (
96  const dictionary& dict,
97  CloudType& cloud
98 )
99 :
100  BinaryCollisionModel<CloudType>(dict, cloud, typeName),
101  Tref_(readScalar(this->coeffDict().lookup("Tref"))),
102  relaxationCollisionNumber_
103  (
104  readScalar(this->coeffDict().lookup("relaxationCollisionNumber"))
105  )
106 {}
107 
108 
109 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
110 
111 template<class CloudType>
114 {}
115 
116 
117 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
118 
119 template<class CloudType>
121 {
122  return true;
123 }
124 
125 
126 template<class CloudType>
128 (
129  const typename CloudType::parcelType& pP,
130  const typename CloudType::parcelType& pQ
131 ) const
132 {
133  const CloudType& cloud(this->owner());
134 
135  label typeIdP = pP.typeId();
136  label typeIdQ = pQ.typeId();
137 
138  scalar dPQ =
139  0.5
140  *(
141  cloud.constProps(typeIdP).d()
142  + cloud.constProps(typeIdQ).d()
143  );
144 
145  scalar omegaPQ =
146  0.5
147  *(
148  cloud.constProps(typeIdP).omega()
149  + cloud.constProps(typeIdQ).omega()
150  );
151 
152  scalar cR = mag(pP.U() - pQ.U());
153 
154  if (cR < VSMALL)
155  {
156  return 0;
157  }
158 
159  scalar mP = cloud.constProps(typeIdP).mass();
160  scalar mQ = cloud.constProps(typeIdQ).mass();
161  scalar mR = mP*mQ/(mP + mQ);
162 
163  // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
164  scalar sigmaTPQ =
165  pi*dPQ*dPQ
166  *pow(2.0*physicoChemical::k.value()*Tref_/(mR*cR*cR), omegaPQ - 0.5)
167  /exp(Foam::lgamma(2.5 - omegaPQ));
168 
169  return sigmaTPQ*cR;
170 }
171 
172 
173 template<class CloudType>
175 (
176  typename CloudType::parcelType& pP,
177  typename CloudType::parcelType& pQ
178 )
179 {
180  CloudType& cloud(this->owner());
181 
182  label typeIdP = pP.typeId();
183  label typeIdQ = pQ.typeId();
184  vector& UP = pP.U();
185  vector& UQ = pQ.U();
186  scalar& EiP = pP.Ei();
187  scalar& EiQ = pQ.Ei();
188 
189  Random& rndGen(cloud.rndGen());
190 
191  scalar inverseCollisionNumber = 1/relaxationCollisionNumber_;
192 
193  // Larsen Borgnakke internal energy redistribution part. Using the serial
194  // application of the LB method, as per the INELRS subroutine in Bird's
195  // DSMC0R.FOR
196 
197  scalar preCollisionEiP = EiP;
198  scalar preCollisionEiQ = EiQ;
199 
200  direction iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
201  direction iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
202 
203  scalar omegaPQ =
204  0.5
205  *(
206  cloud.constProps(typeIdP).omega()
207  + cloud.constProps(typeIdQ).omega()
208  );
209 
210  scalar mP = cloud.constProps(typeIdP).mass();
211  scalar mQ = cloud.constProps(typeIdQ).mass();
212  scalar mR = mP*mQ/(mP + mQ);
213  vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
214  scalar cRsqr = magSqr(UP - UQ);
215  scalar availableEnergy = 0.5*mR*cRsqr;
216  scalar ChiB = 2.5 - omegaPQ;
217 
218  if (iDofP > 0)
219  {
220  if (inverseCollisionNumber > rndGen.scalar01())
221  {
222  availableEnergy += preCollisionEiP;
223 
224  if (iDofP == 2)
225  {
226  scalar energyRatio = 1.0 - pow(rndGen.scalar01(), (1.0/ChiB));
227  EiP = energyRatio*availableEnergy;
228  }
229  else
230  {
231  scalar ChiA = 0.5*iDofP;
232  EiP = energyRatio(ChiA, ChiB)*availableEnergy;
233  }
234 
235  availableEnergy -= EiP;
236  }
237  }
238 
239  if (iDofQ > 0)
240  {
241  if (inverseCollisionNumber > rndGen.scalar01())
242  {
243  availableEnergy += preCollisionEiQ;
244 
245  if (iDofQ == 2)
246  {
247  scalar energyRatio = 1.0 - pow(rndGen.scalar01(), (1.0/ChiB));
248  EiQ = energyRatio*availableEnergy;
249  }
250  else
251  {
252  scalar ChiA = 0.5*iDofQ;
253  EiQ = energyRatio(ChiA, ChiB)*availableEnergy;
254  }
255 
256  availableEnergy -= EiQ;
257  }
258  }
259 
260  // Rescale the translational energy
261  scalar cR = sqrt(2.0*availableEnergy/mR);
262 
263  // Variable Hard Sphere collision part
264  scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
265  scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
266  scalar phi = twoPi*rndGen.scalar01();
267 
268  vector postCollisionRelU =
269  cR
270  *vector
271  (
272  cosTheta,
273  sinTheta*cos(phi),
274  sinTheta*sin(phi)
275  );
276 
277  UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
278  UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
279 }
280 
281 
282 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
Simple random number generator.
Definition: Random.H:49
unsigned char direction
Definition: direction.H:43
cachedRandom rndGen(label(0),-1)
virtual void collide(typename CloudType::parcelType &pP, typename CloudType::parcelType &pQ)
Apply collision.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual bool active() const
Flag to indicate whether model activates collision model.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual scalar sigmaTcR(const typename CloudType::parcelType &pP, const typename CloudType::parcelType &pQ) const
Return the collision cross section * relative velocity product.
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 exp(const dimensionedScalar &ds)
DSMCCloud< dsmcParcel > CloudType
Variable Hard Sphere BinaryCollision Model with Larsen Borgnakke internal energy redistribution. Based on the INELRS subroutine in Bird&#39;s DSMC0R.FOR.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:217
dictionary dict
stressControl lookup("compactNormalStress") >> compactNormalStress
LarsenBorgnakkeVariableHardSphere(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar lgamma(const dimensionedScalar &ds)
surfaceScalarField & phi
Random & rndGen()
Return refernce to the random object.
Definition: DSMCCloudI.H:121
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Templated DSMC particle collision class.
Definition: DSMCCloud.H:55
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
mathematical constants.
const List< typename ParcelType::constantProperties > & constProps() const
Return all of the constant properties.
Definition: DSMCCloudI.H:95
const scalar twoPi(2 *pi)
const dimensionedScalar k
Boltzmann constant.
dimensionedScalar sin(const dimensionedScalar &ds)