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