src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "reitzKHRT.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "mathematicalConstants.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(reitzKHRT, 0);
  38
  39 addToRunTimeSelectionTable
  40 (
  41     breakupModel,
  42     reitzKHRT,
  43     dictionary
  44 );
  45
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 // Construct from components
  50 reitzKHRT::reitzKHRT
  51 (
  52     const dictionary& dict,
  53     spray& sm
  54 )
  55 :
  56     breakupModel(dict, sm),
  57     coeffsDict_(dict.subDict(typeName + "Coeffs")),
  58     g_(sm.g()),
  59     b0_(readScalar(coeffsDict_.lookup("B0"))),
  60     b1_(readScalar(coeffsDict_.lookup("B1"))),
  61     cTau_(readScalar(coeffsDict_.lookup("Ctau"))),
  62     cRT_(readScalar(coeffsDict_.lookup("CRT"))),
  63     msLimit_(readScalar(coeffsDict_.lookup("msLimit"))),
  64     weberLimit_(readScalar(coeffsDict_.lookup("WeberLimit")))
  65 {}
  66
  67
  68 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  69
  70 reitzKHRT::~reitzKHRT()
  71 {}
  72
  73
  74 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  75
  76 void reitzKHRT::breakupParcel
  77 (
  78     parcel& p,
  79     const scalar deltaT,
  80     const vector& vel,
  81     const liquidMixture& fuels
  82 ) const
  83 {
  84
  85     label celli = p.cell();
  86     scalar T = p.T();
  87     scalar r = 0.5*p.d();
  88     scalar pc = spray_.p()[celli];
  89
  90     scalar sigma = fuels.sigma(pc, T, p.X());
  91     scalar rhoLiquid = fuels.rho(pc, T, p.X());
  92     scalar muLiquid = fuels.mu(pc, T, p.X());
  93     scalar rhoGas = spray_.rho()[celli];
  94     scalar Np = p.N(rhoLiquid);
  95     scalar semiMass = Np*pow(p.d(), 3.0);
  96
  97     scalar weGas      = p.We(vel, rhoGas, sigma);
  98     scalar weLiquid   = p.We(vel, rhoLiquid, sigma);
  99     // correct the Reynolds number. Reitz is using radius instead of diameter
 100     scalar reLiquid   = 0.5*p.Re(rhoLiquid, vel, muLiquid);
 101     scalar ohnesorge  = sqrt(weLiquid)/(reLiquid + VSMALL);
 102     scalar taylor     = ohnesorge*sqrt(weGas);
 103
 104     vector acceleration = p.Urel(vel)/p.tMom();
 105     vector trajectory = p.U()/mag(p.U());
 106     scalar gt = (g_ + acceleration) & trajectory;
 107
 108     // frequency of the fastest growing KH-wave
 109     scalar omegaKH =
 110     (
 111         0.34 + 0.38*pow(weGas, 1.5)
 112     )/
 113     (
 114         (1 + ohnesorge)*(1 + 1.4*pow(taylor, 0.6))
 115     )*sqrt(sigma/(rhoLiquid*pow(r, 3)) );
 116
 117     // corresponding KH wave-length.
 118     scalar lambdaKH =
 119     9.02*r*
 120     (
 121         1.0 + 0.45*sqrt(ohnesorge)
 122     )*
 123     (
 124         1.0 + 0.4*pow(taylor, 0.7)
 125     )/
 126     pow
 127     (
 128         (
 129             1.0 + 0.865*pow(weGas, 1.67)
 130         ),
 131         0.6
 132     );
 133
 134     // characteristic Kelvin-Helmholtz breakup time
 135     scalar tauKH = 3.726*b1_*r/(omegaKH*lambdaKH);
 136
 137     // stable KH diameter
 138     scalar dc = 2.0*b0_*lambdaKH;
 139
 140     // the frequency of the fastest growing RT wavelength.
 141     scalar helpVariable = mag(gt*(rhoLiquid - rhoGas));
 142     scalar omegaRT = sqrt
 143     (
 144         2.0*pow(helpVariable, 1.5)
 145        /(3.0*sqrt(3.0*sigma)*(rhoGas + rhoLiquid))
 146     );
 147
 148     // RT wave number
 149     scalar KRT = sqrt(helpVariable/(3.0*sigma + VSMALL));
 150
 151     // Wavelength of the fastest growing Raleigh-Taylor frequency
 152     scalar lambdaRT = 2.0*mathematicalConstant::pi*cRT_/(KRT + VSMALL);
 153
 154     // if lambdaRT < diameter, then RT waves are growing on the surface
 155     // and we start to keep track of how long they have been growing
 156     if ((p.ct() > 0) || (lambdaRT < p.d()))
 157     {
 158         p.ct() += deltaT;
 159     }
 160
 161     // characteristic RT breakup time
 162     scalar tauRT = cTau_/(omegaRT + VSMALL);
 163
 164     // check if we have RT breakup
 165     if ((p.ct() > tauRT) && (lambdaRT < p.d()))
 166     {
 167         // the RT breakup creates diameter/lambdaRT new droplets
 168         p.ct() = -GREAT;
 169         scalar multiplier = p.d()/lambdaRT;
 170         scalar nDrops = multiplier*Np;
 171         p.d() = cbrt(semiMass/nDrops);
 172     }
 173     // otherwise check for KH breakup
 174     else if (dc < p.d())
 175     {
 176         // no breakup below Weber = 12
 177         if (weGas > weberLimit_)
 178         {
 179
 180             label injector = label(p.injector());
 181             scalar fraction = deltaT/tauKH;
 182
 183             // reduce the diameter according to the rate-equation
 184             p.d() = (fraction*dc + p.d())/(1.0 + fraction);
 185
 186             scalar dc3 = pow(dc, 3.0);
 187             scalar ms = rhoLiquid*Np*dc3*mathematicalConstant::pi/6.0;
 188             p.ms() += ms;
 189
 190             label nParcels = spray_.injectors()[injector].properties()->nParcelsToInject
 191             (
 192                 spray_.injectors()[injector].properties()->tsoi(),
 193                 spray_.injectors()[injector].properties()->teoi()
 194             );
 195
 196             scalar averageParcelMass = spray_.injectors()[injector].properties()->mass()/nParcels;
 197
 198             if
 199             (
 200                 (p.ms()/averageParcelMass > msLimit_)
 201             )
 202             {
 203                 // set the initial ms value to -GREAT. This prevents
 204                 // new droplets from being formed from the child droplet
 205                 // from the KH instability
 206
 207                 // mass of stripped child parcel
 208                 scalar mc = p.ms();
 209                 // Prevent child parcel from taking too much mass
 210                 if (mc > 0.5*p.m())
 211                 {
 212                     mc = 0.5*p.m();
 213                 }
 214
 215                 spray_.addParticle
 216                 (
 217                     new parcel
 218                     (
 219                         spray_,
 220                         p.position(),
 221                         p.cell(),
 222                         p.n(),
 223                         dc,
 224                         p.T(),
 225                         mc,
 226                         0.0,
 227                         0.0,
 228                         0.0,
 229                         -GREAT,
 230                         p.tTurb(),
 231                         0.0,
 232                         p.injector(),
 233                         p.U(),
 234                         p.Uturb(),
 235                         p.X(),
 236                         p.fuelNames()
 237                     )
 238                 );
 239
 240                 p.m() -= mc;
 241                 p.ms() = 0.0;
 242             }
 243         }
 244     }
 245 }
 246
 247
 248 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 249
 250 } // End namespace Foam
 251
 252 // ************************************************************************* //