src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.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 "error.H"
  27
  28 #include "RutlandFlashBoil.H"
  29 #include "addToRunTimeSelectionTable.H"
  30 #include "mathematicalConstants.H"
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 defineTypeNameAndDebug(RutlandFlashBoil, 0);
  39
  40 addToRunTimeSelectionTable
  41 (
  42     evaporationModel,
  43     RutlandFlashBoil,
  44     dictionary
  45 );
  46
  47
  48 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  49
  50 // Construct from dictionary
  51 RutlandFlashBoil::RutlandFlashBoil
  52 (
  53     const dictionary& dict
  54 )
  55 :
  56     evaporationModel(dict),
  57     evapDict_(dict.subDict(typeName + "Coeffs")),
  58     preReScFactor_(readScalar(evapDict_.lookup("preReScFactor"))),
  59     ReExponent_(readScalar(evapDict_.lookup("ReExponent"))),
  60     ScExponent_(readScalar(evapDict_.lookup("ScExponent"))),
  61     evaporationScheme_(evapDict_.lookup("evaporationScheme")),
  62     nEvapIter_(0)
  63 {
  64     if (evaporationScheme_ == "implicit")
  65     {
  66         nEvapIter_ = 2;
  67     }
  68     else if (evaporationScheme_ == "explicit")
  69     {
  70         nEvapIter_ = 1;
  71     }
  72     else
  73     {
  74         FatalError
  75             << "evaporationScheme type " << evaporationScheme_
  76             << " unknown.\n"
  77             << "Use implicit or explicit."
  78             << abort(FatalError);
  79     }
  80 }
  81
  82
  83 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  84
  85 RutlandFlashBoil::~RutlandFlashBoil()
  86 {}
  87
  88
  89 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  90
  91 bool RutlandFlashBoil::evaporation() const
  92 {
  93     return true;
  94 }
  95
  96 // Correlation for the Sherwood Number
  97 scalar RutlandFlashBoil::Sh
  98 (
  99     const scalar ReynoldsNumber,
 100     const scalar SchmidtNumber
 101 ) const
 102 {
 103     return 2.0 + preReScFactor_*pow(ReynoldsNumber,ReExponent_)*pow(SchmidtNumber,ScExponent_);
 104 }
 105
 106 scalar RutlandFlashBoil::relaxationTime
 107 (
 108     const scalar diameter,
 109     const scalar liquidDensity,
 110     const scalar rhoFuelVapor,
 111     const scalar massDiffusionCoefficient,
 112     const scalar ReynoldsNumber,
 113     const scalar SchmidtNumber,
 114     const scalar Xs,
 115     const scalar Xf,
 116     const scalar m0,
 117     const scalar dm,
 118     const scalar dt
 119 ) const
 120 {
 121     scalar time = GREAT;
 122     scalar lgExpr = 0.0;
 123
 124     /*
 125         (pressure - partialFuelVaporPressure)/
 126         (pressure - pressureAtSurface)
 127       = 1 + Xratio
 128
 129         if the pressure @ Surface > pressure
 130         this lead to boiling
 131         and Xratio -> infinity (as it should)
 132         ... this is numerically nasty
 133
 134     NB! by N. Nordin
 135         X_v,s = (p_v,s/p) X_v,d
 136         where X_v,d = 1 for single component fuel
 137         according to eq (3.136)
 138         in D. Clerides Thesis
 139     */
 140
 141     scalar Xratio = (Xs - Xf)/max(SMALL, 1.0 - Xs);
 142
 143     if (Xratio > 0.0)
 144     {
 145         lgExpr = log(1.0 + Xratio);
 146     }
 147
 148     // From Equation (3.79) in C. Kralj's Thesis:
 149     // Note that the 2.0 (instead of 6.0) below is correct, since evaporation
 150     // is d(diameter)/dt and not d(mass)/dt
 151
 152     scalar Sherwood = Sh(ReynoldsNumber, SchmidtNumber);
 153
 154     scalar FbExp = 0.7;
 155
 156     scalar logXratio = log(1.0+Xratio);
 157     scalar Fb = 1.0;
 158
 159     if(logXratio > SMALL)
 160     {
 161         Fb = pow((1.0 + Xratio),FbExp) * log(1.0+Xratio)/Xratio;
 162     }
 163
 164 //  TL: proposed correction to sherwood number, implemented
 165
 166     Sherwood = 2.0 + (Sherwood - 2.0)/Fb;
 167
 168     scalar denominator =
 169         6.0 * massDiffusionCoefficient
 170       * Sherwood
 171       * rhoFuelVapor * lgExpr;
 172
 173     if (denominator > SMALL)
 174     {
 175         time = max(VSMALL, liquidDensity * pow(diameter, 2.0)/denominator);
 176     }
 177
 178     return time;
 179 }
 180
 181
 182 scalar RutlandFlashBoil::boilingTime
 183 (
 184     const scalar liquidDensity,
 185     const scalar cpFuel,
 186     const scalar heatOfVapour,
 187     const scalar kappa,
 188     const scalar Nusselt,
 189     const scalar deltaTemp,
 190     const scalar diameter,
 191     const scalar liquidCore,
 192     const scalar ct,
 193     const scalar tDrop,
 194     const scalar tBoilingSurface,
 195     const scalar vapourSurfaceEnthalpy,
 196     const scalar vapourFarEnthalpy,
 197     const scalar cpGas,
 198     const scalar temperature,
 199     const scalar kLiq
 200 ) const
 201 {
 202
 203     scalar time = GREAT;
 204
 205     // the deltaTemperature is limited to not go below .5 deg K
 206     // for numerical reasons.
 207     // This is probably not important, since it results in an upper
 208     // limit for the boiling time... which we have anyway.
 209
 210     //  TL kSet to the k value at the droplet temperature, not as in the Rutland Paper
 211
 212     if(liquidCore > 0.5)
 213     {
 214         if(tDrop > tBoilingSurface)
 215         {
 216             //  Evaporation of the liquid sheet
 217
 218             scalar psi = 2.72;
 219             scalar kIncreased = psi * kLiq;
 220             scalar alfa = psi * kIncreased/(liquidDensity * cpFuel);
 221             scalar F = alfa * ct/sqr(0.5 * diameter);
 222
 223             scalar expSum = 0.0;
 224             scalar expSumOld = expSum;
 225
 226             label Niter = 200;
 227
 228             for(label k=0; k < Niter ; k++)
 229             {
 230                 expSum += exp(sqr(-k*mathematicalConstant::pi*sqrt(F)/2.0));
 231                 if(mag(expSum-expSumOld)/expSum < 1.0e-3)
 232                 {
 233                     break;
 234                 }
 235                 expSumOld = expSum;
 236             }
 237         }
 238     }
 239     else
 240     {
 241         scalar dTLB = min(0.5, tDrop -  tBoilingSurface);
 242         scalar alfaS = 0.0;
 243
 244         if(dTLB >= 0.0 &&  dTLB < 5.0)
 245         {
 246             alfaS = 0.76 * pow(dTLB, 0.26);
 247         }
 248         if(dTLB >= 5.0 &&  dTLB < 25.0)
 249         {
 250             alfaS = 0.027 * pow(dTLB, 2.33);
 251         }
 252         if(dTLB >= 25.0)
 253         {
 254             alfaS = 13.8 * pow(dTLB, 0.39);
 255         }
 256
 257         scalar Gf =
 258         (
 259             4.0 * alfaS * dTLB * mathematicalConstant::pi * sqr(diameter/2.0)
 260         )
 261         /
 262         heatOfVapour;
 263
 264         //  calculation of the heat transfer vapourization at superheated conditions (temperature>tBoilingSurface)
 265         scalar G = 0.0;
 266         if(temperature > tBoilingSurface)
 267         {
 268             scalar NusseltCorr = Nusselt ;
 269             scalar A = mag((vapourFarEnthalpy-vapourSurfaceEnthalpy)/heatOfVapour);
 270
 271             // TL : 2.0? or 1.0? try 1!
 272             scalar B = 1.0*mathematicalConstant::pi*kappa/cpGas*diameter*NusseltCorr;
 273             scalar nPos = B * log(1.0 + A)/Gf + 1.0;
 274             scalar nNeg = (1.0/A)*(exp(Gf/B) - 1.0 - A) + 1.0;
 275
 276             scalar Gpos = Gf*nPos;
 277             scalar Gneg = Gf/nNeg;
 278
 279             //scalar FgPos = Gpos + Gf - B * log( 1.0 + ( 1.0 + Gf/Gpos ) * A);
 280             scalar FgNeg = Gneg + Gf - B * log( 1.0 + ( 1.0 + Gf/Gneg ) * A);
 281
 282             if(FgNeg > 0.0)
 283             {
 284                 for(label j = 0; j < 20; j++)
 285                 {
 286                     Gneg = Gneg/10.0;
 287                     Gneg = max(Gneg, VSMALL);
 288                     FgNeg = Gneg + Gf - B * log( 1.0 + ( 1.0 + Gf/Gneg ) * A);
 289                     if(FgNeg < 0.0)
 290                     {
 291                         break;
 292                     }
 293                 }
 294             }
 295
 296             FgNeg = Gneg + Gf - B * log( 1.0 + ( 1.0 + Gf/Gneg ) * A);
 297
 298             G = 0.5*(Gpos+Gneg);
 299             scalar Gold = -100;
 300
 301             label Niter = 200;
 302             label k=0;
 303
 304             if(FgNeg > 0.0)
 305             {
 306                 Info << "no convergence" << endl;
 307             }
 308
 309
 310             if(FgNeg < 0.0)
 311             {
 312
 313                 for(k=0; k<Niter ; k++)
 314                 {
 315
 316                     scalar Fg = G + Gf - B * log( 1.0 + ( 1.0 + Gf/G ) * A);
 317
 318                     if(Fg > 0)
 319                     {
 320                         Gpos = G;
 321                         G = 0.5*(Gpos+Gneg);
 322                     }
 323                     else
 324                     {
 325                         Gneg = G;
 326                         G = 0.5*(Gpos+Gneg);
 327                     }
 328
 329                     Gold = G;
 330                     if(mag(G-Gold)/Gold < 1.0e-3)
 331                     {
 332                         break;
 333                     }
 334                 }
 335
 336                 if(k >= Niter - 1)
 337                 {
 338                     Info << " No convergence for G " << endl;
 339                 }
 340             }
 341             else
 342             {
 343                 G = 0.0;
 344             }
 345         }
 346
 347         time = ((4.0/3.0)*mathematicalConstant::pi*pow(diameter/2.0,3.0))*liquidDensity/(G+Gf);
 348         time = max(VSMALL, time);
 349     }
 350
 351     return time;
 352 }
 353
 354 inline label RutlandFlashBoil::nEvapIter() const
 355 {
 356     return nEvapIter_;
 357 }
 358 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 359
 360 } // End namespace Foam
 361
 362 // ************************************************************************* //