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