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