FreeFOAM The Cross-Platform CFD Toolkit
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) 2009-2010 OpenCFD Ltd.
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 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 template <class CloudType>
32 (
33  scalar ChiA,
34  scalar ChiB
35 )
36 {
37  CloudType& cloud(this->owner());
38 
39  Random& rndGen(cloud.rndGen());
40 
41  scalar ChiAMinusOne = ChiA - 1;
42 
43  scalar ChiBMinusOne = ChiB - 1;
44 
45  if (ChiAMinusOne < SMALL && ChiBMinusOne < SMALL)
46  {
47  return rndGen.scalar01();
48  }
49 
50  scalar energyRatio;
51 
52  scalar P;
53 
54  do
55  {
56  P = 0;
57 
58  energyRatio = rndGen.scalar01();
59 
60  if (ChiAMinusOne < SMALL)
61  {
62  P = 1.0 - pow(energyRatio, ChiB);
63  }
64  else if (ChiBMinusOne < SMALL)
65  {
66  P = 1.0 - pow(energyRatio, ChiA);
67  }
68  else
69  {
70  P =
71  pow
72  (
73  (ChiAMinusOne + ChiBMinusOne)*energyRatio/ChiAMinusOne,
74  ChiAMinusOne
75  )
76  *pow
77  (
78  (ChiAMinusOne + ChiBMinusOne)*(1 - energyRatio)
79  /ChiBMinusOne,
80  ChiBMinusOne
81  );
82  }
83  } while (P < rndGen.scalar01());
84 
85  return energyRatio;
86 }
87 
88 
89 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
90 
91 template <class CloudType>
93 (
94  const dictionary& dict,
95  CloudType& cloud
96 )
97 :
98  BinaryCollisionModel<CloudType>(dict, cloud, typeName),
99  Tref_(readScalar(this->coeffDict().lookup("Tref"))),
100  relaxationCollisionNumber_
101  (
102  readScalar(this->coeffDict().lookup("relaxationCollisionNumber"))
103  )
104 {}
105 
106 
107 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
108 
109 template <class CloudType>
112 {}
113 
114 
115 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
116 
117 
118 template <class CloudType>
120 (
121  label typeIdP,
122  label typeIdQ,
123  const vector& UP,
124  const vector& UQ
125 ) const
126 {
127  const CloudType& cloud(this->owner());
128 
129  scalar dPQ =
130  0.5
131  *(
132  cloud.constProps(typeIdP).d()
133  + cloud.constProps(typeIdQ).d()
134  );
135 
136  scalar omegaPQ =
137  0.5
138  *(
139  cloud.constProps(typeIdP).omega()
140  + cloud.constProps(typeIdQ).omega()
141  );
142 
143  scalar cR = mag(UP - UQ);
144 
145  if (cR < VSMALL)
146  {
147  return 0;
148  }
149 
150  scalar mP = cloud.constProps(typeIdP).mass();
151 
152  scalar mQ = cloud.constProps(typeIdQ).mass();
153 
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 =
159  *pow(2.0*CloudType::kb*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  label typeIdP,
170  label typeIdQ,
171  vector& UP,
172  vector& UQ,
173  scalar& EiP,
174  scalar& EiQ
175 )
176 {
177  CloudType& cloud(this->owner());
178 
179  Random& rndGen(cloud.rndGen());
180 
181  scalar inverseCollisionNumber = 1/relaxationCollisionNumber_;
182 
183  // Larsen Borgnakke internal energy redistribution part. Using the serial
184  // application of the LB method, as per the INELRS subroutine in Bird's
185  // DSMC0R.FOR
186 
187  scalar preCollisionEiP = EiP;
188 
189  scalar preCollisionEiQ = EiQ;
190 
191  scalar iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
192 
193  scalar iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
194 
195  scalar omegaPQ =
196  0.5
197  *(
198  cloud.constProps(typeIdP).omega()
199  + cloud.constProps(typeIdQ).omega()
200  );
201 
202  scalar mP = cloud.constProps(typeIdP).mass();
203 
204  scalar mQ = cloud.constProps(typeIdQ).mass();
205 
206  scalar mR = mP*mQ/(mP + mQ);
207 
208  vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
209 
210  scalar cRsqr = magSqr(UP - UQ);
211 
212  scalar availableEnergy = 0.5*mR*cRsqr;
213 
214  scalar ChiB = 2.5 - omegaPQ;
215 
216  if (iDofP > 0)
217  {
218  if (inverseCollisionNumber > rndGen.scalar01())
219  {
220  availableEnergy += preCollisionEiP;
221 
222  scalar ChiA = 0.5*iDofP;
223 
224  EiP = energyRatio(ChiA, ChiB)*availableEnergy;
225 
226  availableEnergy -= EiP;
227  }
228  }
229 
230  if (iDofQ > 0)
231  {
232  if (inverseCollisionNumber > rndGen.scalar01())
233  {
234  availableEnergy += preCollisionEiQ;
235 
236  // Change to general LB ratio calculation
237  scalar energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));
238 
239  EiQ = energyRatio*availableEnergy;
240 
241  availableEnergy -= EiQ;
242  }
243  }
244 
245  // Rescale the translational energy
246  scalar cR = sqrt(2.0*availableEnergy/mR);
247 
248  // Variable Hard Sphere collision part
249 
250  scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
251 
252  scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
253 
254  scalar phi = 2.0*mathematicalConstant::pi*rndGen.scalar01();
255 
256  vector postCollisionRelU =
257  cR
258  *vector
259  (
260  cosTheta,
261  sinTheta*cos(phi),
262  sinTheta*sin(phi)
263  );
264 
265  UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
266 
267  UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
268 }
269 
270 
271 // ************************ vim: set sw=4 sts=4 et: ************************ //