src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.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 "VariableHardSphere.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 template <class CloudType>
  31 Foam::VariableHardSphere<CloudType>::VariableHardSphere
  32 (
  33     const dictionary& dict,
  34     CloudType& cloud
  35 )
  36 :
  37     BinaryCollisionModel<CloudType>(dict, cloud, typeName),
  38     Tref_(readScalar(this->coeffDict().lookup("Tref")))
  39 {}
  40
  41
  42 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  43
  44 template <class CloudType>
  45 Foam::VariableHardSphere<CloudType>::~VariableHardSphere()
  46 {}
  47
  48
  49 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  50
  51
  52 template <class CloudType>
  53 Foam::scalar Foam::VariableHardSphere<CloudType>::sigmaTcR
  54 (
  55     label typeIdP,
  56     label typeIdQ,
  57     const vector& UP,
  58     const vector& UQ
  59 ) const
  60 {
  61     const CloudType& cloud(this->owner());
  62
  63     scalar dPQ =
  64         0.5
  65        *(
  66             cloud.constProps(typeIdP).d()
  67           + cloud.constProps(typeIdQ).d()
  68         );
  69
  70     scalar omegaPQ =
  71         0.5
  72        *(
  73             cloud.constProps(typeIdP).omega()
  74           + cloud.constProps(typeIdQ).omega()
  75         );
  76
  77     scalar cR = mag(UP - UQ);
  78
  79     if (cR < VSMALL)
  80     {
  81         return 0;
  82     }
  83
  84     scalar mP = cloud.constProps(typeIdP).mass();
  85
  86     scalar mQ = cloud.constProps(typeIdQ).mass();
  87
  88     scalar mR = mP*mQ/(mP + mQ);
  89
  90     // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
  91     scalar sigmaTPQ =
  92         mathematicalConstant::pi*dPQ*dPQ
  93        *pow(2.0*CloudType::kb*Tref_/(mR*cR*cR), omegaPQ - 0.5)
  94        /exp(Foam::lgamma(2.5 - omegaPQ));
  95
  96     return sigmaTPQ*cR;
  97 }
  98
  99
 100 template <class CloudType>
 101 void Foam::VariableHardSphere<CloudType>::collide
 102 (
 103     label typeIdP,
 104     label typeIdQ,
 105     vector& UP,
 106     vector& UQ,
 107     scalar& EiP,
 108     scalar& EiQ
 109 )
 110 {
 111     CloudType& cloud(this->owner());
 112
 113     Random& rndGen(cloud.rndGen());
 114
 115     scalar mP = cloud.constProps(typeIdP).mass();
 116
 117     scalar mQ = cloud.constProps(typeIdQ).mass();
 118
 119     vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
 120
 121     scalar cR = mag(UP - UQ);
 122
 123     scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
 124
 125     scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
 126
 127     scalar phi = 2.0*mathematicalConstant::pi*rndGen.scalar01();
 128
 129     vector postCollisionRelU =
 130         cR
 131        *vector
 132         (
 133             cosTheta,
 134             sinTheta*cos(phi),
 135             sinTheta*sin(phi)
 136         );
 137
 138     UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
 139
 140     UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
 141 }
 142
 143
 144 // ************************************************************************* //