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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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
  26 #include "LarsenBorgnakkeVariableHardSphere.H"
  27 #include "constants.H"
  28
  29 using namespace Foam::constant::mathematical;
  30
  31 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  32
  33 template <class CloudType>
  34 Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
  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
  46     scalar ChiBMinusOne = ChiB - 1;
  47
  48     if (ChiAMinusOne < SMALL && ChiBMinusOne < SMALL)
  49     {
  50         return rndGen.scalar01();
  51     }
  52
  53     scalar energyRatio;
  54
  55     scalar P;
  56
  57     do
  58     {
  59         P = 0;
  60
  61         energyRatio = rndGen.scalar01();
  62
  63         if (ChiAMinusOne < SMALL)
  64         {
  65             P = 1.0 - pow(energyRatio, ChiB);
  66         }
  67         else if (ChiBMinusOne < SMALL)
  68         {
  69             P = 1.0 - pow(energyRatio, ChiA);
  70         }
  71         else
  72         {
  73             P =
  74                 pow
  75                 (
  76                     (ChiAMinusOne + ChiBMinusOne)*energyRatio/ChiAMinusOne,
  77                     ChiAMinusOne
  78                 )
  79                *pow
  80                 (
  81                     (ChiAMinusOne + ChiBMinusOne)*(1 - energyRatio)
  82                     /ChiBMinusOne,
  83                     ChiBMinusOne
  84                 );
  85         }
  86     } while (P < rndGen.scalar01());
  87
  88     return energyRatio;
  89 }
  90
  91
  92 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  93
  94 template <class CloudType>
  95 Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::
  96 LarsenBorgnakkeVariableHardSphere
  97 (
  98     const dictionary& dict,
  99     CloudType& cloud
 100 )
 101 :
 102     BinaryCollisionModel<CloudType>(dict, cloud, typeName),
 103     Tref_(readScalar(this->coeffDict().lookup("Tref"))),
 104     relaxationCollisionNumber_
 105     (
 106         readScalar(this->coeffDict().lookup("relaxationCollisionNumber"))
 107     )
 108 {}
 109
 110
 111 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 112
 113 template <class CloudType>
 114 Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::
 115 ~LarsenBorgnakkeVariableHardSphere()
 116 {}
 117
 118
 119 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 120
 121 template<class CloudType>
 122 bool Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::active() const
 123 {
 124     return true;
 125 }
 126
 127
 128 template <class CloudType>
 129 Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::sigmaTcR
 130 (
 131     const typename CloudType::parcelType& pP,
 132     const typename CloudType::parcelType& pQ
 133 ) const
 134 {
 135     const CloudType& cloud(this->owner());
 136
 137     label typeIdP = pP.typeId();
 138     label typeIdQ = pQ.typeId();
 139
 140     scalar dPQ =
 141         0.5
 142        *(
 143             cloud.constProps(typeIdP).d()
 144           + cloud.constProps(typeIdQ).d()
 145         );
 146
 147     scalar omegaPQ =
 148         0.5
 149        *(
 150             cloud.constProps(typeIdP).omega()
 151           + cloud.constProps(typeIdQ).omega()
 152         );
 153
 154     scalar cR = mag(pP.U() - pQ.U());
 155
 156     if (cR < VSMALL)
 157     {
 158         return 0;
 159     }
 160
 161     scalar mP = cloud.constProps(typeIdP).mass();
 162
 163     scalar mQ = cloud.constProps(typeIdQ).mass();
 164
 165     scalar mR = mP*mQ/(mP + mQ);
 166
 167     // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
 168     scalar sigmaTPQ =
 169         pi*dPQ*dPQ
 170        *pow(2.0*physicoChemical::k.value()*Tref_/(mR*cR*cR), omegaPQ - 0.5)
 171        /exp(Foam::lgamma(2.5 - omegaPQ));
 172
 173     return sigmaTPQ*cR;
 174 }
 175
 176
 177 template <class CloudType>
 178 void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
 179 (
 180     typename CloudType::parcelType& pP,
 181     typename CloudType::parcelType& pQ
 182 )
 183 {
 184     CloudType& cloud(this->owner());
 185
 186     label typeIdP = pP.typeId();
 187     label typeIdQ = pQ.typeId();
 188     vector& UP = pP.U();
 189     vector& UQ = pQ.U();
 190     scalar& EiP = pP.Ei();
 191     scalar& EiQ = pQ.Ei();
 192
 193     Random& rndGen(cloud.rndGen());
 194
 195     scalar inverseCollisionNumber = 1/relaxationCollisionNumber_;
 196
 197     // Larsen Borgnakke internal energy redistribution part.  Using the serial
 198     // application of the LB method, as per the INELRS subroutine in Bird's
 199     // DSMC0R.FOR
 200
 201     scalar preCollisionEiP = EiP;
 202
 203     scalar preCollisionEiQ = EiQ;
 204
 205     scalar iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
 206
 207     scalar iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
 208
 209     scalar omegaPQ =
 210         0.5
 211        *(
 212             cloud.constProps(typeIdP).omega()
 213           + cloud.constProps(typeIdQ).omega()
 214         );
 215
 216     scalar mP = cloud.constProps(typeIdP).mass();
 217
 218     scalar mQ = cloud.constProps(typeIdQ).mass();
 219
 220     scalar mR = mP*mQ/(mP + mQ);
 221
 222     vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
 223
 224     scalar cRsqr = magSqr(UP - UQ);
 225
 226     scalar availableEnergy = 0.5*mR*cRsqr;
 227
 228     scalar ChiB = 2.5 - omegaPQ;
 229
 230     if (iDofP > 0)
 231     {
 232         if (inverseCollisionNumber > rndGen.scalar01())
 233         {
 234             availableEnergy += preCollisionEiP;
 235
 236             scalar ChiA = 0.5*iDofP;
 237
 238             EiP = energyRatio(ChiA, ChiB)*availableEnergy;
 239
 240             availableEnergy -= EiP;
 241         }
 242     }
 243
 244     if (iDofQ > 0)
 245     {
 246         if (inverseCollisionNumber > rndGen.scalar01())
 247         {
 248             availableEnergy += preCollisionEiQ;
 249
 250             // Change to general LB ratio calculation
 251             scalar energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));
 252
 253             EiQ = energyRatio*availableEnergy;
 254
 255             availableEnergy -= EiQ;
 256         }
 257     }
 258
 259     // Rescale the translational energy
 260     scalar cR = sqrt(2.0*availableEnergy/mR);
 261
 262     // Variable Hard Sphere collision part
 263
 264     scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
 265
 266     scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
 267
 268     scalar phi = twoPi*rndGen.scalar01();
 269
 270     vector postCollisionRelU =
 271         cR
 272        *vector
 273         (
 274             cosTheta,
 275             sinTheta*cos(phi),
 276             sinTheta*sin(phi)
 277         );
 278
 279     UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
 280
 281     UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
 282 }
 283
 284
 285 // ************************************************************************* //