src/lagrangian/dieselSpray/spraySubModels/breakupModel/SHF/SHF.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 "SHF.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "mathematicalConstants.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 defineTypeNameAndDebug(SHF, 0);
  39
  40 addToRunTimeSelectionTable
  41 (
  42     breakupModel,
  43     SHF,
  44     dictionary
  45 );
  46
  47
  48 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  49
  50 // Construct from components
  51 SHF::SHF
  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     rndGen_(sm.rndGen()),
  61     weCorrCoeff_(readScalar(coeffsDict_.lookup("weCorrCoeff"))),
  62     weBuCrit_(readScalar(coeffsDict_.lookup("weBuCrit"))),
  63     weBuBag_(readScalar(coeffsDict_.lookup("weBuBag"))),
  64     weBuMM_(readScalar(coeffsDict_.lookup("weBuMM"))),
  65     ohnCoeffCrit_(readScalar(coeffsDict_.lookup("ohnCoeffCrit"))),
  66     ohnCoeffBag_(readScalar(coeffsDict_.lookup("ohnCoeffBag"))),
  67     ohnCoeffMM_(readScalar(coeffsDict_.lookup("ohnCoeffMM"))),
  68     ohnExpCrit_(readScalar(coeffsDict_.lookup("ohnExpCrit"))),
  69     ohnExpBag_(readScalar(coeffsDict_.lookup("ohnExpBag"))),
  70     ohnExpMM_(readScalar(coeffsDict_.lookup("ohnExpMM"))),
  71     cInit_(readScalar(coeffsDict_.lookup("Cinit"))),
  72     c1_(readScalar(coeffsDict_.lookup("C1"))),
  73     c2_(readScalar(coeffsDict_.lookup("C2"))),
  74     c3_(readScalar(coeffsDict_.lookup("C3"))),
  75     cExp1_(readScalar(coeffsDict_.lookup("Cexp1"))),
  76     cExp2_(readScalar(coeffsDict_.lookup("Cexp2"))),
  77     cExp3_(readScalar(coeffsDict_.lookup("Cexp3"))),
  78     weConst_(readScalar(coeffsDict_.lookup("Weconst"))),
  79     weCrit1_(readScalar(coeffsDict_.lookup("Wecrit1"))),
  80     weCrit2_(readScalar(coeffsDict_.lookup("Wecrit2"))),
  81     coeffD_(readScalar(coeffsDict_.lookup("CoeffD"))),
  82     onExpD_(readScalar(coeffsDict_.lookup("OnExpD"))),
  83     weExpD_(readScalar(coeffsDict_.lookup("WeExpD"))),
  84     mu_(readScalar(coeffsDict_.lookup("mu"))),
  85     sigma_(readScalar(coeffsDict_.lookup("sigma"))),
  86     d32Coeff_(readScalar(coeffsDict_.lookup("d32Coeff"))),
  87     cDmaxBM_(readScalar(coeffsDict_.lookup("cDmaxBM"))),
  88     cDmaxS_(readScalar(coeffsDict_.lookup("cDmaxS"))),
  89     corePerc_(readScalar(coeffsDict_.lookup("corePerc")))
  90 {}
  91
  92
  93 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  94
  95 SHF::~SHF()
  96 {}
  97
  98
  99 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 100
 101 void SHF::breakupParcel
 102 (
 103     parcel& p,
 104     const scalar deltaT,
 105     const vector& vel,
 106     const liquidMixture& fuels
 107 ) const
 108 {
 109
 110     label celli = p.cell();
 111     scalar T = p.T();
 112     scalar pc = spray_.p()[celli];
 113
 114     scalar sigma = fuels.sigma(pc, T, p.X());
 115     scalar rhoLiquid = fuels.rho(pc, T, p.X());
 116     scalar muLiquid = fuels.mu(pc, T, p.X());
 117     scalar rhoGas = spray_.rho()[celli];
 118
 119     scalar weGas      = p.We(vel, rhoGas, sigma);
 120     scalar weLiquid   = p.We(vel, rhoLiquid, sigma);
 121
 122     // correct the Reynolds number. Reitz is using radius instead of diameter
 123
 124     scalar reLiquid   = p.Re(rhoLiquid, vel, muLiquid);
 125     scalar ohnesorge  = sqrt(weLiquid)/(reLiquid + VSMALL);
 126
 127     vector acceleration = p.Urel(vel)/p.tMom();
 128     vector trajectory = p.U()/mag(p.U());
 129
 130     vector vRel = p.Urel(vel);
 131
 132     scalar weGasCorr = weGas/(1.0 + weCorrCoeff_ * ohnesorge);
 133
 134     // droplet deformation characteristic time
 135
 136     scalar tChar = p.d()/mag(vRel)*sqrt(rhoLiquid/rhoGas);
 137
 138     scalar tFirst = cInit_ * tChar;
 139
 140     scalar tSecond = 0;
 141     scalar tCharSecond = 0;
 142
 143     bool bag = false;
 144     bool multimode = false;
 145     bool shear = false;
 146     bool success = false;
 147
 148
 149     //  updating the droplet characteristic time
 150     p.ct() += deltaT;
 151
 152     if(weGas > weConst_)
 153     {
 154         if(weGas < weCrit1_)
 155         {
 156             tCharSecond = c1_*pow((weGas - weConst_),cExp1_);
 157         }
 158         else if(weGas >= weCrit1_ && weGas <= weCrit2_)
 159         {
 160             tCharSecond = c2_*pow((weGas - weConst_),cExp2_);
 161         }
 162         else
 163         {
 164             tCharSecond = c3_*pow((weGas - weConst_),cExp3_);
 165         }
 166     }
 167
 168     scalar weC = weBuCrit_*(1.0+ohnCoeffCrit_*pow(ohnesorge,ohnExpCrit_));
 169     scalar weB = weBuBag_*(1.0+ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
 170     scalar weMM = weBuMM_*(1.0+ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));
 171
 172     if(weGas > weC && weGas < weB)
 173     {
 174         bag = true;
 175     }
 176
 177     if(weGas >= weB && weGas <= weMM)
 178     {
 179         multimode = true;
 180     }
 181
 182     if(weGas > weMM)
 183     {
 184         shear = true;
 185     }
 186
 187     tSecond = tCharSecond * tChar;
 188
 189     scalar tBreakUP = tFirst + tSecond;
 190     if(p.ct() > tBreakUP)
 191     {
 192
 193         scalar d32 = coeffD_*p.d()*pow(ohnesorge,onExpD_)*pow(weGasCorr,weExpD_);
 194
 195         if(bag || multimode)
 196         {
 197
 198             scalar d05 = d32Coeff_ * d32;
 199
 200             scalar x = 0.0;
 201             scalar y = 0.0;
 202             scalar d = 0.0;
 203
 204             while(!success)
 205             {
 206                 x = cDmaxBM_*rndGen_.scalar01();
 207                 d = sqr(x)*d05;
 208                 y = rndGen_.scalar01();
 209
 210                 scalar p = x/(2.0*sqrt(2.0*mathematicalConstant::pi)*sigma_)*exp(-0.5*sqr((x-mu_)/sigma_));
 211
 212                 if (y<p)
 213                 {
 214                     success = true;
 215                 }
 216             }
 217
 218             p.d() = d;
 219             p.ct() = 0.0;
 220         }
 221
 222         if(shear)
 223         {
 224             scalar dC = weConst_*sigma/(rhoGas*sqr(mag(vRel)));
 225             scalar d32Red = 4.0*(d32 * dC)/(5.0 * dC - d32);
 226             scalar initMass = p.m();
 227
 228             scalar d05 = d32Coeff_ * d32Red;
 229
 230             scalar x = 0.0;
 231             scalar y = 0.0;
 232             scalar d = 0.0;
 233
 234             while(!success)
 235             {
 236
 237                 x = cDmaxS_*rndGen_.scalar01();
 238                 d = sqr(x)*d05;
 239                 y = rndGen_.scalar01();
 240
 241                 scalar p = x/(2.0*sqrt(2.0*mathematicalConstant::pi)*sigma_)*exp(-0.5*sqr((x-mu_)/sigma_));
 242
 243                 if (y<p)
 244                 {
 245                     success = true;
 246                 }
 247             }
 248
 249             p.d() = dC;
 250             p.m() = corePerc_ * initMass;
 251
 252             spray_.addParticle
 253             (
 254                 new parcel
 255                 (
 256                     spray_,
 257                     p.position(),
 258                     p.cell(),
 259                     p.n(),
 260                     d,
 261                     p.T(),
 262                     (1.0 - corePerc_)* initMass,
 263                     0.0,
 264                     0.0,
 265                     0.0,
 266                     -GREAT,
 267                     p.tTurb(),
 268                     0.0,
 269                     scalar(p.injector()),
 270                     p.U(),
 271                     p.Uturb(),
 272                     p.X(),
 273                     p.fuelNames()
 274                 )
 275             );
 276
 277             p.ct() = 0.0;
 278         }
 279     }
 280 }
 281
 282
 283 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 284
 285 } // End namespace Foam
 286
 287 // ************************************************************************* //