src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-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 "TAB.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 template <class CloudType>
  31 Foam::TAB<CloudType>::TAB
  32 (
  33     const dictionary& dict,
  34     CloudType& owner
  35 )
  36 :
  37     BreakupModel<CloudType>(dict, owner, typeName),
  38     Cmu_(BreakupModel<CloudType>::TABCmu_),
  39     Comega_(BreakupModel<CloudType>::TABComega_),
  40     WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
  41     SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod"))
  42 {
  43     // calculate the inverse function of the Rossin-Rammler Distribution
  44     const scalar xx0 = 12.0;
  45     const scalar rrd100 =
  46         1.0/(1.0-exp(-xx0)*(1.0 + xx0 + sqr(xx0)/2.0 + pow3(xx0)/6.0));
  47
  48     forAll(rrd_, n)
  49     {
  50         scalar xx = 0.12*(n + 1);
  51         rrd_[n] =
  52             (1.0 - exp(-xx)*(1.0 + xx + sqr(xx)/2.0 + pow3(xx)/6.0))*rrd100;
  53     }
  54
  55     if (!BreakupModel<CloudType>::solveOscillationEq_)
  56     {
  57         Info<< "Warning: solveOscillationEq is set to "
  58             << BreakupModel<CloudType>::solveOscillationEq_ << nl
  59             << " Setting it to true in order for the TAB model to work."
  60             << endl;
  61         BreakupModel<CloudType>::solveOscillationEq_ = true;
  62     }
  63
  64     if (SMDCalcMethod_ == "method1")
  65     {
  66         SMDMethod_ = method1;
  67     }
  68     else if (SMDCalcMethod_ == "method2")
  69     {
  70         SMDMethod_ = method2;
  71     }
  72     else
  73     {
  74         SMDMethod_ = method2;
  75         Info<< "Warning: SMDCalculationMethod unknown. Options are "
  76             "(method1 | method2). Using method2" << endl;
  77     }
  78 }
  79
  80
  81 template <class CloudType>
  82 Foam::TAB<CloudType>::TAB(const TAB<CloudType>& bum)
  83 :
  84     BreakupModel<CloudType>(bum),
  85     Cmu_(bum.Cmu_),
  86     Comega_(bum.Comega_),
  87     WeCrit_(bum.WeCrit_),
  88     SMDCalcMethod_(bum.SMDCalcMethod_)
  89 {}
  90
  91
  92 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  93
  94 template <class CloudType>
  95 Foam::TAB<CloudType>::~TAB()
  96 {}
  97
  98
  99 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 100
 101 template<class CloudType>
 102 bool Foam::TAB<CloudType>::update
 103 (
 104     const scalar dt,
 105     const vector& g,
 106     scalar& d,
 107     scalar& tc,
 108     scalar& ms,
 109     scalar& nParticle,
 110     scalar& KHindex,
 111     scalar& y,
 112     scalar& yDot,
 113     const scalar d0,
 114     const scalar rho,
 115     const scalar mu,
 116     const scalar sigma,
 117     const vector& U,
 118     const scalar rhoc,
 119     const scalar muc,
 120     const vector& Urel,
 121     const scalar Urmag,
 122     const scalar tMom,
 123     const scalar averageParcelMass,
 124     scalar& dChild,
 125     scalar& massChild,
 126     cachedRandom& rndGen
 127 ) const
 128 {
 129     scalar r = 0.5*d;
 130     scalar r2 = r*r;
 131     scalar r3 = r*r2;
 132
 133     scalar semiMass = nParticle*pow3(d);
 134
 135     // inverse of characteristic viscous damping time
 136     scalar rtd = 0.5*Cmu_*mu/(rho*r2);
 137
 138     // oscillation frequency (squared)
 139     scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
 140
 141     if (omega2 > 0)
 142     {
 143         scalar omega = sqrt(omega2);
 144         scalar We = rhoc*sqr(Urmag)*r/sigma;
 145         scalar Wetmp = We/WeCrit_;
 146
 147         scalar y1 = y - Wetmp;
 148         scalar y2 = yDot/omega;
 149
 150         scalar a = sqrt(y1*y1 + y2*y2);
 151
 152         // scotty we may have break-up
 153         if (a+Wetmp > 1.0)
 154         {
 155             scalar phic = y1/a;
 156
 157             // constrain phic within -1 to 1
 158             phic = max(min(phic, 1), -1);
 159
 160             scalar phit = acos(phic);
 161             scalar phi = phit;
 162             scalar quad = -y2/a;
 163             if (quad < 0)
 164             {
 165                 phi = 2.0*constant::mathematical::pi - phit;
 166             }
 167
 168             scalar tb = 0;
 169
 170             if (mag(y) < 1.0)
 171             {
 172                 scalar coste = 1.0;
 173                 if ((Wetmp - a < -1) && (yDot < 0))
 174                 {
 175                     coste = -1.0;
 176                 }
 177
 178                 scalar theta = acos((coste-Wetmp)/a);
 179
 180                 if (theta < phi)
 181                 {
 182                     if (2*constant::mathematical::pi-theta >= phi)
 183                     {
 184                         theta = -theta;
 185                     }
 186                     theta += 2*constant::mathematical::pi;
 187                 }
 188                 tb = (theta-phi)/omega;
 189
 190                 // breakup occurs
 191                 if (dt > tb)
 192                 {
 193                     y = 1.0;
 194                     yDot = -a*omega*sin(omega*tb + phi);
 195                 }
 196
 197             }
 198
 199             // update droplet size
 200             if (dt > tb)
 201             {
 202                 scalar rs =
 203                     r/(1.0 + (4.0/3.0)*sqr(y) + rho*r3/(8*sigma)*sqr(yDot));
 204
 205                 label n = 0;
 206                 scalar rNew = 0.0;
 207                 switch (SMDMethod_)
 208                 {
 209                     case method1:
 210                     {
 211                         #include "TABSMDCalcMethod1.H"
 212                         break;
 213                     }
 214                     case method2:
 215                     {
 216                         #include "TABSMDCalcMethod2.H"
 217                         break;
 218                     }
 219                 }
 220
 221                 if (rNew < r)
 222                 {
 223                     d = 2*rNew;
 224                     y = 0;
 225                     yDot = 0;
 226                 }
 227             }
 228         }
 229     }
 230     else
 231     {
 232         // reset droplet distortion parameters
 233         y = 0;
 234         yDot = 0;
 235     }
 236
 237     // update the nParticle count to conserve mass
 238     nParticle = semiMass/pow3(d);
 239
 240     // Do not add child parcel
 241     return false;
 242 }
 243
 244
 245 // ************************************************************************* //