src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "GradientDispersionRAS.H"
  28
  29 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  30
  31 template<class CloudType>
  32 Foam::GradientDispersionRAS<CloudType>::GradientDispersionRAS
  33 (
  34     const dictionary& dict,
  35     CloudType& owner
  36 )
  37 :
  38     DispersionRASModel<CloudType>(dict, owner),
  39     gradkPtr_(NULL)
  40 {}
  41
  42
  43 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  44
  45 template<class CloudType>
  46 Foam::GradientDispersionRAS<CloudType>::~GradientDispersionRAS()
  47 {
  48     cacheFields(false);
  49 }
  50
  51
  52 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  53
  54 template<class CloudType>
  55 bool Foam::GradientDispersionRAS<CloudType>::active() const
  56 {
  57     return true;
  58 }
  59
  60
  61 template<class CloudType>
  62 void Foam::GradientDispersionRAS<CloudType>::cacheFields(const bool store)
  63 {
  64     DispersionRASModel<CloudType>::cacheFields(store);
  65
  66     if (store)
  67     {
  68         gradkPtr_ = fvc::grad(*this->kPtr_).ptr();
  69     }
  70     else
  71     {
  72         if (gradkPtr_)
  73         {
  74             delete gradkPtr_;
  75             gradkPtr_ = NULL;
  76         }
  77     }
  78 }
  79
  80
  81 template<class CloudType>
  82 Foam::vector Foam::GradientDispersionRAS<CloudType>::update
  83 (
  84     const scalar dt,
  85     const label celli,
  86     const vector& U,
  87     const vector& Uc,
  88     vector& UTurb,
  89     scalar& tTurb
  90 )
  91 {
  92     const scalar cps = 0.16432;
  93
  94     const volScalarField& k = *this->kPtr_;
  95     const volScalarField& epsilon = *this->epsilonPtr_;
  96     const volVectorField& gradk = *this->gradkPtr_;
  97
  98     const scalar UrelMag = mag(U - Uc - UTurb);
  99
 100     const scalar tTurbLoc = min
 101     (
 102         k[celli]/epsilon[celli],
 103         cps*pow(k[celli], 1.5)/epsilon[celli]/(UrelMag + SMALL)
 104     );
 105
 106     // Parcel is perturbed by the turbulence
 107     if (dt < tTurbLoc)
 108     {
 109         tTurb += dt;
 110
 111         if (tTurb > tTurbLoc)
 112         {
 113             tTurb = 0.0;
 114
 115             scalar sigma = sqrt(2.0*k[celli]/3.0);
 116             vector dir = -gradk[celli]/(mag(gradk[celli]) + SMALL);
 117
 118             // Numerical Recipes... Ch. 7. Random Numbers...
 119             scalar x1 = 0.0;
 120             scalar x2 = 0.0;
 121             scalar rsq = 10.0;
 122             while ((rsq > 1.0) || (rsq == 0.0))
 123             {
 124                 x1 = 2.0*this->owner().rndGen().scalar01() - 1.0;
 125                 x2 = 2.0*this->owner().rndGen().scalar01() - 1.0;
 126                 rsq = x1*x1 + x2*x2;
 127             }
 128
 129             scalar fac = sqrt(-2.0*log(rsq)/rsq);
 130
 131             // In 2D calculations the -grad(k) is always
 132             // away from the axis of symmetry
 133             // This creates a 'hole' in the spray and to
 134             // prevent this we let x1 be both negative/positive
 135             if (this->owner().mesh().nSolutionD() == 2)
 136             {
 137                 fac *= x1;
 138             }
 139             else
 140             {
 141                 fac *= mag(x1);
 142             }
 143
 144             UTurb = sigma*fac*dir;
 145         }
 146     }
 147     else
 148     {
 149         tTurb = GREAT;
 150         UTurb = vector::zero;
 151     }
 152
 153     return Uc + UTurb;
 154 }
 155
 156
 157 // ************************************************************************* //