30 template <
class CloudType>
37 CloudType& cloud(this->owner());
39 Random& rndGen(cloud.rndGen());
41 scalar ChiAMinusOne = ChiA - 1;
43 scalar ChiBMinusOne = ChiB - 1;
45 if (ChiAMinusOne < SMALL && ChiBMinusOne < SMALL)
47 return rndGen.scalar01();
58 energyRatio = rndGen.scalar01();
60 if (ChiAMinusOne < SMALL)
62 P = 1.0 -
pow(energyRatio, ChiB);
64 else if (ChiBMinusOne < SMALL)
66 P = 1.0 -
pow(energyRatio, ChiA);
73 (ChiAMinusOne + ChiBMinusOne)*energyRatio/ChiAMinusOne,
78 (ChiAMinusOne + ChiBMinusOne)*(1 - energyRatio)
83 }
while (P < rndGen.scalar01());
91 template <
class CloudType>
99 Tref_(
readScalar(this->coeffDict().lookup(
"Tref"))),
100 relaxationCollisionNumber_
102 readScalar(this->coeffDict().lookup(
"relaxationCollisionNumber"))
109 template <
class CloudType>
118 template <
class CloudType>
127 const CloudType& cloud(this->owner());
132 cloud.constProps(typeIdP).d()
133 + cloud.constProps(typeIdQ).d()
139 cloud.constProps(typeIdP).omega()
140 + cloud.constProps(typeIdQ).omega()
143 scalar cR =
mag(UP - UQ);
150 scalar mP = cloud.constProps(typeIdP).mass();
152 scalar mQ = cloud.constProps(typeIdQ).mass();
154 scalar mR = mP*mQ/(mP + mQ);
159 *
pow(2.0*CloudType::kb*Tref_/(mR*cR*cR), omegaPQ - 0.5)
166 template <
class CloudType>
177 CloudType& cloud(this->owner());
179 Random& rndGen(cloud.rndGen());
181 scalar inverseCollisionNumber = 1/relaxationCollisionNumber_;
187 scalar preCollisionEiP = EiP;
189 scalar preCollisionEiQ = EiQ;
191 scalar iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
193 scalar iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
198 cloud.constProps(typeIdP).omega()
199 + cloud.constProps(typeIdQ).omega()
202 scalar mP = cloud.constProps(typeIdP).mass();
204 scalar mQ = cloud.constProps(typeIdQ).mass();
206 scalar mR = mP*mQ/(mP + mQ);
208 vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
210 scalar cRsqr =
magSqr(UP - UQ);
212 scalar availableEnergy = 0.5*mR*cRsqr;
214 scalar ChiB = 2.5 - omegaPQ;
218 if (inverseCollisionNumber > rndGen.scalar01())
220 availableEnergy += preCollisionEiP;
222 scalar ChiA = 0.5*iDofP;
224 EiP = energyRatio(ChiA, ChiB)*availableEnergy;
226 availableEnergy -= EiP;
232 if (inverseCollisionNumber > rndGen.scalar01())
234 availableEnergy += preCollisionEiQ;
237 scalar energyRatio = 1.0 -
pow(rndGen.scalar01(),(1.0/ChiB));
239 EiQ = energyRatio*availableEnergy;
241 availableEnergy -= EiQ;
246 scalar cR =
sqrt(2.0*availableEnergy/mR);
250 scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
252 scalar sinTheta =
sqrt(1.0 - cosTheta*cosTheta);
256 vector postCollisionRelU =
265 UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
267 UQ = Ucm - postCollisionRelU*mP/(mP + mQ);