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

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