LarsenBorgnakkeVariableHardSphere.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-2020 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_(this->coeffDict().template lookup<scalar>("Tref")),
102  relaxationCollisionNumber_
103  (
104  this->coeffDict().template lookup<scalar>("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  const typename CloudType::parcelType& pP,
123  const typename CloudType::parcelType& pQ
124 ) const
125 {
126  const CloudType& cloud(this->owner());
127 
128  label typeIdP = pP.typeId();
129  label typeIdQ = pQ.typeId();
130 
131  scalar dPQ =
132  0.5
133  *(
134  cloud.constProps(typeIdP).d()
135  + cloud.constProps(typeIdQ).d()
136  );
137 
138  scalar omegaPQ =
139  0.5
140  *(
141  cloud.constProps(typeIdP).omega()
142  + cloud.constProps(typeIdQ).omega()
143  );
144 
145  scalar cR = mag(pP.U() - pQ.U());
146 
147  if (cR < vSmall)
148  {
149  return 0;
150  }
151 
152  scalar mP = cloud.constProps(typeIdP).mass();
153  scalar mQ = cloud.constProps(typeIdQ).mass();
154  scalar mR = mP*mQ/(mP + mQ);
155 
156  // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
157  scalar sigmaTPQ =
158  pi*dPQ*dPQ
159  *pow(2.0*physicoChemical::k.value()*Tref_/(mR*cR*cR), omegaPQ - 0.5)
160  /exp(Foam::lgamma(2.5 - omegaPQ));
161 
162  return sigmaTPQ*cR;
163 }
164 
165 
166 template<class CloudType>
168 (
169  typename CloudType::parcelType& pP,
170  typename CloudType::parcelType& pQ
171 )
172 {
173  CloudType& cloud(this->owner());
174 
175  label typeIdP = pP.typeId();
176  label typeIdQ = pQ.typeId();
177  vector& UP = pP.U();
178  vector& UQ = pQ.U();
179  scalar& EiP = pP.Ei();
180  scalar& EiQ = pQ.Ei();
181 
182  Random& rndGen(cloud.rndGen());
183 
184  scalar inverseCollisionNumber = 1/relaxationCollisionNumber_;
185 
186  // Larsen Borgnakke internal energy redistribution part. Using the serial
187  // application of the LB method, as per the INELRS subroutine in Bird's
188  // DSMC0R.FOR
189 
190  scalar preCollisionEiP = EiP;
191  scalar preCollisionEiQ = EiQ;
192 
193  direction iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
194  direction iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
195 
196  scalar omegaPQ =
197  0.5
198  *(
199  cloud.constProps(typeIdP).omega()
200  + cloud.constProps(typeIdQ).omega()
201  );
202 
203  scalar mP = cloud.constProps(typeIdP).mass();
204  scalar mQ = cloud.constProps(typeIdQ).mass();
205  scalar mR = mP*mQ/(mP + mQ);
206  vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
207  scalar cRsqr = magSqr(UP - UQ);
208  scalar availableEnergy = 0.5*mR*cRsqr;
209  scalar ChiB = 2.5 - omegaPQ;
210 
211  if (iDofP > 0)
212  {
213  if (inverseCollisionNumber > rndGen.scalar01())
214  {
215  availableEnergy += preCollisionEiP;
216 
217  if (iDofP == 2)
218  {
219  scalar energyRatio = 1.0 - pow(rndGen.scalar01(), (1.0/ChiB));
220  EiP = energyRatio*availableEnergy;
221  }
222  else
223  {
224  scalar ChiA = 0.5*iDofP;
225  EiP = energyRatio(ChiA, ChiB)*availableEnergy;
226  }
227 
228  availableEnergy -= EiP;
229  }
230  }
231 
232  if (iDofQ > 0)
233  {
234  if (inverseCollisionNumber > rndGen.scalar01())
235  {
236  availableEnergy += preCollisionEiQ;
237 
238  if (iDofQ == 2)
239  {
240  scalar energyRatio = 1.0 - pow(rndGen.scalar01(), (1.0/ChiB));
241  EiQ = energyRatio*availableEnergy;
242  }
243  else
244  {
245  scalar ChiA = 0.5*iDofQ;
246  EiQ = energyRatio(ChiA, ChiB)*availableEnergy;
247  }
248 
249  availableEnergy -= EiQ;
250  }
251  }
252 
253  // Rescale the translational energy
254  scalar cR = sqrt(2.0*availableEnergy/mR);
255 
256  // Variable Hard Sphere collision part
257  scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
258  scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
259  scalar phi = twoPi*rndGen.scalar01();
260 
261  vector postCollisionRelU =
262  cR
263  *vector
264  (
265  cosTheta,
266  sinTheta*cos(phi),
267  sinTheta*sin(phi)
268  );
269 
270  UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
271  UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
272 }
273 
274 
275 // ************************************************************************* //
virtual void collide(typename CloudType::parcelType &pP, typename CloudType::parcelType &pQ)
Apply collision.
dictionary dict
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
DSMCCloud< dsmcParcel > CloudType
virtual scalar sigmaTcR(const typename CloudType::parcelType &pP, const typename CloudType::parcelType &pQ) const
Return the collision cross section * relative velocity product.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
uint8_t direction
Definition: direction.H:45
Templated DSMC particle collision class.
Definition: DSMCCloud.H:55
dimensionedScalar sqrt(const dimensionedScalar &ds)
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Random rndGen(label(0))
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar exp(const dimensionedScalar &ds)
mathematical constants.
Random & rndGen()
Return references to the random object.
Definition: DSMCCloudI.H:121
phi
Definition: correctPhi.H:3
Random number generator.
Definition: Random.H:57
const scalar twoPi(2 *pi)
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:221
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar sin(const dimensionedScalar &ds)
dimensionedScalar lgamma(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
LarsenBorgnakkeVariableHardSphere(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
const dimensionedScalar k
Boltzmann constant.
dimensioned< scalar > mag(const dimensioned< Type > &)
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
const List< typename ParcelType::constantProperties > & constProps() const
Return all of the constant properties.
Definition: DSMCCloudI.H:95
Variable Hard Sphere BinaryCollision Model with Larsen Borgnakke internal energy redistribution. Based on the INELRS subroutine in Bird&#39;s DSMC0R.FOR.