src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.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 "ETAB.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 template<class CloudType>
  31 Foam::ETAB<CloudType>::ETAB
  32 (
  33     const dictionary& dict,
  34     CloudType& owner
  35 )
  36 :
  37     BreakupModel<CloudType>(dict, owner, typeName),
  38     Cmu_(10.0),
  39     Comega_(8.0),
  40     k1_(0.2),
  41     k2_(0.2),
  42     WeCrit_(12.0),
  43     WeTransition_(100.0),
  44     AWe_(0.0)
  45 {
  46     if (!this->defaultCoeffs(true))
  47     {
  48         this->coeffDict().lookup("Cmu") >> Cmu_;
  49         this->coeffDict().lookup("Comega") >> Comega_;
  50         this->coeffDict().lookup("k1") >> k1_;
  51         this->coeffDict().lookup("k2") >> k2_;
  52         this->coeffDict().lookup("WeCrit") >> WeCrit_;
  53         this->coeffDict().lookup("WeTransition") >> WeTransition_;
  54     }
  55
  56     scalar k21 = k2_/k1_;
  57     AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow4(WeTransition_);
  58
  59     if (!BreakupModel<CloudType>::solveOscillationEq_)
  60     {
  61         Info<< "Warning: solveOscillationEq is set to "
  62             << BreakupModel<CloudType>::solveOscillationEq_ << nl
  63             << " Setting it to true in order for the ETAB model to work."
  64             << endl;
  65         BreakupModel<CloudType>::solveOscillationEq_ = true;
  66     }
  67 }
  68
  69
  70 template<class CloudType>
  71 Foam::ETAB<CloudType>::ETAB(const ETAB<CloudType>& bum)
  72 :
  73     BreakupModel<CloudType>(bum),
  74     Cmu_(bum.Cmu_),
  75     Comega_(bum.Comega_),
  76     k1_(bum.k1_),
  77     k2_(bum.k2_),
  78     WeCrit_(bum.WeCrit_),
  79     WeTransition_(bum.WeTransition_),
  80     AWe_(bum.AWe_)
  81 {}
  82
  83
  84 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  85
  86 template<class CloudType>
  87 Foam::ETAB<CloudType>::~ETAB()
  88 {}
  89
  90
  91 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  92
  93 template<class CloudType>
  94 bool Foam::ETAB<CloudType>::update
  95 (
  96     const scalar dt,
  97     const vector& g,
  98     scalar& d,
  99     scalar& tc,
 100     scalar& ms,
 101     scalar& nParticle,
 102     scalar& KHindex,
 103     scalar& y,
 104     scalar& yDot,
 105     const scalar d0,
 106     const scalar rho,
 107     const scalar mu,
 108     const scalar sigma,
 109     const vector& U,
 110     const scalar rhoc,
 111     const scalar muc,
 112     const vector& Urel,
 113     const scalar Urmag,
 114     const scalar tMom,
 115     const scalar averageParcelMass,
 116     scalar& dChild,
 117     scalar& massChild,
 118     cachedRandom& rndGen
 119 ) const
 120 {
 121     scalar r  = 0.5*d;
 122     scalar r2 = r*r;
 123     scalar r3 = r*r2;
 124
 125     scalar semiMass = nParticle*pow3(d);
 126
 127     // inverse of characteristic viscous damping time
 128     scalar rtd = 0.5*Cmu_*mu/(rho*r2);
 129
 130     // oscillation frequency (squared)
 131     scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
 132
 133     if (omega2 > 0)
 134     {
 135         scalar omega = sqrt(omega2);
 136         scalar romega = 1.0/omega;
 137
 138         scalar We = rhoc*sqr(Urmag)*r/sigma;
 139         scalar Wetmp = We/WeCrit_;
 140
 141         scalar y1 = y - Wetmp;
 142         scalar y2 = yDot*romega;
 143
 144         scalar a = sqrt(y1*y1 + y2*y2);
 145
 146         // scotty we may have break-up
 147         if (a + Wetmp > 1.0)
 148         {
 149             scalar phic = y1/a;
 150
 151             // constrain phic within -1 to 1
 152             phic = max(min(phic, 1), -1);
 153
 154             scalar phit = acos(phic);
 155             scalar phi = phit;
 156             scalar quad = -y2/a;
 157             if (quad < 0)
 158             {
 159                 phi = 2*constant::mathematical::pi - phit;
 160             }
 161
 162             scalar tb = 0;
 163
 164             if (mag(y) < 1.0)
 165             {
 166                 scalar theta = acos((1.0 - Wetmp)/a);
 167
 168                 if (theta < phi)
 169                 {
 170                     if (2*constant::mathematical::pi - theta >= phi)
 171                     {
 172                         theta = -theta;
 173                     }
 174                     theta += 2*constant::mathematical::pi;
 175                 }
 176                 tb = (theta - phi)*romega;
 177
 178                 // breakup occurs
 179                 if (dt > tb)
 180                 {
 181                     y = 1.0;
 182                     yDot = -a*omega*sin(omega*tb + phi);
 183                 }
 184             }
 185
 186             // update droplet size
 187             if (dt > tb)
 188             {
 189                 scalar sqrtWe = AWe_*pow4(We) + 1.0;
 190                 scalar Kbr = k1_*omega*sqrtWe;
 191
 192                 if (We > WeTransition_)
 193                 {
 194                     sqrtWe = sqrt(We);
 195                     Kbr =k2_*omega*sqrtWe;
 196                 }
 197
 198                 scalar rWetmp = 1.0/Wetmp;
 199                 scalar cosdtbu = max(-1.0, min(1.0, 1.0 - rWetmp));
 200                 scalar dtbu = romega*acos(cosdtbu);
 201                 scalar decay = exp(-Kbr*dtbu);
 202
 203                 scalar rNew = decay*r;
 204                 if (rNew < r)
 205                 {
 206                     d = 2.0*rNew;
 207                     y = 0.0;
 208                     yDot = 0.0;
 209                 }
 210             }
 211         }
 212     }
 213     else
 214     {
 215         // reset droplet distortion parameters
 216         y = 0;
 217         yDot = 0;
 218     }
 219
 220     // update the nParticle count to conserve mass
 221     nParticle = semiMass/pow3(d);
 222
 223     // Do not add child parcel
 224     return false;
 225 }
 226
 227
 228 // ************************************************************************* //