src/lagrangian/intermediate/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011-2011 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
  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 "BrownianMotionForce.H"
  27 #include "mathematicalConstants.H"
  28 #include "demandDrivenData.H"
  29
  30 using namespace Foam::constant;
  31
  32 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  33
  34 template<class CloudType>
  35 Foam::scalar Foam::BrownianMotionForce<CloudType>::erfInv(const scalar y) const
  36 {
  37     const scalar a = 0.147;
  38     scalar k = 2.0/(mathematical::pi*a) +  0.5*log(1.0 - y*y);
  39     scalar h = log(1.0 - y*y)/a;
  40     scalar x = sqrt(-k + sqrt(k*k - h));
  41
  42     if (y < 0.0)
  43     {
  44         return -x;
  45     }
  46     else
  47     {
  48         return x;
  49     }
  50 }
  51
  52
  53 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  54
  55 template<class CloudType>
  56 Foam::BrownianMotionForce<CloudType>::BrownianMotionForce
  57 (
  58     CloudType& owner,
  59     const fvMesh& mesh,
  60     const dictionary& dict
  61 )
  62 :
  63     ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
  64     rndGen_(owner.rndGen()),
  65     lambda_(readScalar(this->coeffs().lookup("lambda"))),
  66     turbulence_(readBool(this->coeffs().lookup("turbulence"))),
  67     turbulenceModelPtr_(NULL),
  68     kPtr_(NULL),
  69     ownK_(false)
  70 {
  71     if (turbulence_)
  72     {
  73         HashTable<const compressible::turbulenceModel*> models =
  74             this->mesh().objectRegistry::template lookupClass
  75             <
  76                 compressible::turbulenceModel
  77             >();
  78
  79         if (models.size() == 1)
  80         {
  81             turbulenceModelPtr_ = models.begin()();
  82         }
  83         else
  84         {
  85             FatalErrorIn
  86             (
  87                 "Foam::BrownianMotionForce<CloudType>::BrownianMotionForce"
  88                 "("
  89                     "CloudType&, "
  90                     "const fvMesh&, "
  91                     "const dictionary&"
  92                 ")"
  93             )   << "Unable to find a valid turbulence model in mesh database"
  94                 << exit(FatalError);
  95         }
  96     }
  97 }
  98
  99
 100 template<class CloudType>
 101 Foam::BrownianMotionForce<CloudType>::BrownianMotionForce
 102 (
 103     const BrownianMotionForce& bmf
 104 )
 105 :
 106     ParticleForce<CloudType>(bmf),
 107     rndGen_(bmf.rndGen_),
 108     lambda_(bmf.lambda_),
 109     turbulence_(bmf.turbulence_),
 110     turbulenceModelPtr_(NULL),
 111     kPtr_(NULL),
 112     ownK_(false)
 113 {}
 114
 115
 116 // * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
 117
 118 template<class CloudType>
 119 Foam::BrownianMotionForce<CloudType>::~BrownianMotionForce()
 120 {
 121     turbulenceModelPtr_ = NULL;
 122 }
 123
 124
 125 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 126
 127 template<class CloudType>
 128 void Foam::BrownianMotionForce<CloudType>::cacheFields(const bool store)
 129 {
 130     if (turbulence_)
 131     {
 132         if (store)
 133         {
 134             tmp<volScalarField> tk = turbulenceModelPtr_->k();
 135             if (tk.isTmp())
 136             {
 137                 kPtr_ = tk.ptr();
 138                 ownK_ = true;
 139             }
 140             else
 141             {
 142                 kPtr_ = tk.operator->();
 143                 ownK_ = false;
 144             }
 145         }
 146         else
 147         {
 148             if (ownK_ && kPtr_)
 149             {
 150                 deleteDemandDrivenData(kPtr_);
 151                 ownK_ = false;
 152             }
 153         }
 154     }
 155 }
 156
 157
 158 template<class CloudType>
 159 Foam::forceSuSp Foam::BrownianMotionForce<CloudType>::calcCoupled
 160 (
 161     const typename CloudType::parcelType& p,
 162     const scalar dt,
 163     const scalar mass,
 164     const scalar Re,
 165     const scalar muc
 166 ) const
 167 {
 168     forceSuSp value(vector::zero, 0.0);
 169
 170     const scalar dp = p.d();
 171     const scalar Tc = p.Tc();
 172
 173     const scalar eta = rndGen_.sample01<scalar>();
 174     const scalar alpha = 2.0*lambda_/dp;
 175     const scalar cc = 1.0 + alpha*(1.257 + 0.4*exp(-1.1/alpha));
 176
 177     const scalar sigma = physicoChemical::sigma.value();
 178
 179     scalar f = 0.0;
 180     if (turbulence_)
 181     {
 182         const label cellI = p.cell();
 183         const volScalarField& k = *kPtr_;
 184         const scalar kc = k[cellI];
 185         const scalar Dp = sigma*Tc*cc/(3*mathematical::pi*muc*dp);
 186         f = eta/mass*sqrt(2.0*sqr(kc)*sqr(Tc)/(Dp*dt));
 187     }
 188     else
 189     {
 190         const scalar rhoRatio = p.rho()/p.rhoc();
 191         const scalar s0 =
 192             216*muc*sigma*Tc/(sqr(mathematical::pi)*pow5(dp)*(rhoRatio)*cc);
 193         f = eta*sqrt(mathematical::pi*s0/dt);
 194     }
 195
 196     const scalar sqrt2 = sqrt(2.0);
 197     for (label i = 0; i < 3; i++)
 198     {
 199         const scalar x = rndGen_.sample01<scalar>();
 200         const scalar eta = sqrt2*erfInv(2*x - 1.0);
 201         value.Su()[i] = mass*f*eta;
 202     }
 203
 204     return value;
 205 }
 206
 207
 208 // ************************************************************************* //