src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.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 "SHF.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 template <class CloudType>
  31 Foam::SHF<CloudType>::SHF
  32 (
  33     const dictionary& dict,
  34     CloudType& owner
  35 )
  36 :
  37     BreakupModel<CloudType>(dict, owner, typeName),
  38     weCorrCoeff_(readScalar(this->coeffDict().lookup("weCorrCoeff"))),
  39     weBuCrit_(readScalar(this->coeffDict().lookup("weBuCrit"))),
  40     weBuBag_(readScalar(this->coeffDict().lookup("weBuBag"))),
  41     weBuMM_(readScalar(this->coeffDict().lookup("weBuMM"))),
  42     ohnCoeffCrit_(readScalar(this->coeffDict().lookup("ohnCoeffCrit"))),
  43     ohnCoeffBag_(readScalar(this->coeffDict().lookup("ohnCoeffBag"))),
  44     ohnCoeffMM_(readScalar(this->coeffDict().lookup("ohnCoeffMM"))),
  45     ohnExpCrit_(readScalar(this->coeffDict().lookup("ohnExpCrit"))),
  46     ohnExpBag_(readScalar(this->coeffDict().lookup("ohnExpBag"))),
  47     ohnExpMM_(readScalar(this->coeffDict().lookup("ohnExpMM"))),
  48     cInit_(readScalar(this->coeffDict().lookup("Cinit"))),
  49     c1_(readScalar(this->coeffDict().lookup("C1"))),
  50     c2_(readScalar(this->coeffDict().lookup("C2"))),
  51     c3_(readScalar(this->coeffDict().lookup("C3"))),
  52     cExp1_(readScalar(this->coeffDict().lookup("Cexp1"))),
  53     cExp2_(readScalar(this->coeffDict().lookup("Cexp2"))),
  54     cExp3_(readScalar(this->coeffDict().lookup("Cexp3"))),
  55     weConst_(readScalar(this->coeffDict().lookup("Weconst"))),
  56     weCrit1_(readScalar(this->coeffDict().lookup("Wecrit1"))),
  57     weCrit2_(readScalar(this->coeffDict().lookup("Wecrit2"))),
  58     coeffD_(readScalar(this->coeffDict().lookup("CoeffD"))),
  59     onExpD_(readScalar(this->coeffDict().lookup("OnExpD"))),
  60     weExpD_(readScalar(this->coeffDict().lookup("WeExpD"))),
  61     mu_(readScalar(this->coeffDict().lookup("mu"))),
  62     sigma_(readScalar(this->coeffDict().lookup("sigma"))),
  63     d32Coeff_(readScalar(this->coeffDict().lookup("d32Coeff"))),
  64     cDmaxBM_(readScalar(this->coeffDict().lookup("cDmaxBM"))),
  65     cDmaxS_(readScalar(this->coeffDict().lookup("cDmaxS"))),
  66     corePerc_(readScalar(this->coeffDict().lookup("corePerc")))
  67 {}
  68
  69
  70 template <class CloudType>
  71 Foam::SHF<CloudType>::SHF(const SHF<CloudType>& bum)
  72 :
  73     BreakupModel<CloudType>(bum),
  74     weCorrCoeff_(bum.weCorrCoeff_),
  75     weBuCrit_(bum.weBuCrit_),
  76     weBuBag_(bum.weBuBag_),
  77     weBuMM_(bum.weBuMM_),
  78     ohnCoeffCrit_(bum.ohnCoeffCrit_),
  79     ohnCoeffBag_(bum.ohnCoeffBag_),
  80     ohnCoeffMM_(bum.ohnCoeffMM_),
  81     ohnExpCrit_(bum.ohnExpCrit_),
  82     ohnExpBag_(bum.ohnExpBag_),
  83     ohnExpMM_(bum.ohnExpMM_),
  84     cInit_(bum.cInit_),
  85     c1_(bum.c1_),
  86     c2_(bum.c2_),
  87     c3_(bum.c3_),
  88     cExp1_(bum.cExp1_),
  89     cExp2_(bum.cExp2_),
  90     cExp3_(bum.cExp3_),
  91     weConst_(bum.weConst_),
  92     weCrit1_(bum.weCrit1_),
  93     weCrit2_(bum.weCrit2_),
  94     coeffD_(bum.coeffD_),
  95     onExpD_(bum.onExpD_),
  96     weExpD_(bum.weExpD_),
  97     mu_(bum.mu_),
  98     sigma_(bum.sigma_),
  99     d32Coeff_(bum.d32Coeff_),
 100     cDmaxBM_(bum.cDmaxBM_),
 101     cDmaxS_(bum.cDmaxS_),
 102     corePerc_(bum.corePerc_)
 103 {}
 104
 105
 106 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 107
 108 template <class CloudType>
 109 Foam::SHF<CloudType>::~SHF()
 110 {}
 111
 112
 113 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 114
 115 template<class CloudType>
 116 bool Foam::SHF<CloudType>::update
 117 (
 118     const scalar dt,
 119     const vector& g,
 120     scalar& d,
 121     scalar& tc,
 122     scalar& ms,
 123     scalar& nParticle,
 124     scalar& KHindex,
 125     scalar& y,
 126     scalar& yDot,
 127     const scalar d0,
 128     const scalar rho,
 129     const scalar mu,
 130     const scalar sigma,
 131     const vector& U,
 132     const scalar rhoc,
 133     const scalar muc,
 134     const vector& Urel,
 135     const scalar Urmag,
 136     const scalar tMom,
 137     const scalar averageParcelMass,
 138     scalar& dChild,
 139     scalar& massChild,
 140     cachedRandom& rndGen
 141 ) const
 142 {
 143     bool addChild = false;
 144
 145     scalar d03 = pow3(d);
 146     scalar rhopi6 = rho*constant::mathematical::pi/6.0;
 147     scalar mass0 = nParticle*rhopi6*d03;
 148     scalar mass = mass0;
 149
 150     scalar weGas = 0.5*rhoc*sqr(Urmag)*d/sigma;
 151     scalar weLiquid = 0.5*rho*sqr(Urmag)*d/sigma;
 152
 153     // correct the Reynolds number. Reitz is using radius instead of diameter
 154     scalar reLiquid   = 0.5*Urmag*d/mu;
 155     scalar ohnesorge  = sqrt(weLiquid)/(reLiquid + VSMALL);
 156
 157     scalar weGasCorr = weGas/(1.0 + weCorrCoeff_*ohnesorge);
 158
 159     // droplet deformation characteristic time
 160
 161     scalar tChar = d/Urmag*sqrt(rho/rhoc);
 162
 163     scalar tFirst = cInit_*tChar;
 164
 165     scalar tSecond = 0;
 166     scalar tCharSecond = 0;
 167
 168     bool bag = false;
 169     bool multimode = false;
 170     bool shear = false;
 171     bool success = false;
 172
 173
 174     // update the droplet characteristic time
 175     tc += dt;
 176
 177     if (weGas > weConst_)
 178     {
 179         if(weGas < weCrit1_)
 180         {
 181             tCharSecond = c1_*pow((weGas - weConst_), cExp1_);
 182         }
 183         else if(weGas >= weCrit1_ && weGas <= weCrit2_)
 184         {
 185             tCharSecond = c2_*pow((weGas - weConst_), cExp2_);
 186         }
 187         else
 188         {
 189             tCharSecond = c3_*pow((weGas - weConst_), cExp3_);
 190         }
 191     }
 192
 193     scalar weC = weBuCrit_*(1.0 + ohnCoeffCrit_*pow(ohnesorge, ohnExpCrit_));
 194     scalar weB = weBuBag_*(1.0 + ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
 195     scalar weMM = weBuMM_*(1.0 + ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));
 196
 197     if (weGas > weC && weGas < weB)
 198     {
 199         bag = true;
 200     }
 201
 202     if (weGas >= weB && weGas <= weMM)
 203     {
 204         multimode = true;
 205     }
 206
 207     if (weGas > weMM)
 208     {
 209         shear = true;
 210     }
 211
 212     tSecond = tCharSecond*tChar;
 213
 214     scalar tBreakUP = tFirst + tSecond;
 215     if (tc > tBreakUP)
 216     {
 217         scalar d32 = coeffD_*d*pow(ohnesorge, onExpD_)*pow(weGasCorr, weExpD_);
 218
 219         if (bag || multimode)
 220         {
 221             scalar d05 = d32Coeff_ * d32;
 222
 223             scalar x = 0.0;
 224             scalar yGuess = 0.0;
 225             scalar dGuess = 0.0;
 226
 227             while(!success)
 228             {
 229                 x = cDmaxBM_*rndGen.sample01<scalar>();
 230                 dGuess = sqr(x)*d05;
 231                 yGuess = rndGen.sample01<scalar>();
 232
 233                 scalar p =
 234                     x
 235                    /(2.0*sqrt(2.0*constant::mathematical::pi)*sigma_)
 236                    *exp(-0.5*sqr((x - mu_)/sigma_));
 237
 238                 if (yGuess < p)
 239                 {
 240                     success = true;
 241                 }
 242             }
 243
 244             d = dGuess;
 245             tc = 0.0;
 246         }
 247
 248         if (shear)
 249         {
 250             scalar dC = weConst_*sigma/(rhoc*sqr(Urmag));
 251             scalar d32Red = 4.0*(d32*dC)/(5.0*dC - d32);
 252
 253             scalar d05 = d32Coeff_ * d32Red;
 254
 255             scalar x = 0.0;
 256             scalar yGuess = 0.0;
 257             scalar dGuess = 0.0;
 258
 259             while(!success)
 260             {
 261                 x = cDmaxS_*rndGen.sample01<scalar>();
 262                 dGuess = sqr(x)*d05;
 263                 yGuess = rndGen.sample01<scalar>();
 264
 265                 scalar p =
 266                     x
 267                    /(2.0*sqrt(2.0*constant::mathematical::pi)*sigma_)
 268                    *exp(-0.5*sqr((x - mu_)/sigma_));
 269
 270                 if (yGuess<p)
 271                 {
 272                     success = true;
 273                 }
 274             }
 275
 276             d = dC;
 277             dChild = dGuess;
 278             massChild = corePerc_*mass0;
 279             mass -= massChild;
 280
 281             addChild = true;
 282             // reset timer
 283             tc = 0.0;
 284         }
 285
 286         // correct nParticle to conserve mass
 287         nParticle = mass/(rhopi6*pow3(d));
 288     }
 289
 290     return addChild;
 291 }
 292
 293
 294 // ************************************************************************* //