src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "LarsenBorgnakkeVariableHardSphere.H"
  27
  28 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  29
  30 template <class CloudType>
  31 Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
  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>
  92 Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::LarsenBorgnakkeVariableHardSphere
  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>
 110 Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::
 111 ~LarsenBorgnakkeVariableHardSphere()
 112 {}
 113
 114
 115 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 116
 117
 118 template <class CloudType>
 119 Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::sigmaTcR
 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 =
 158         mathematicalConstant::pi*dPQ*dPQ
 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>
 167 void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
 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 // ************************************************************************* //