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