src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.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 "error.H"
  27
  28 #include "RutlandFlashBoil.H"
  29 #include "addToRunTimeSelectionTable.H"
  30 #include "mathematicalConstants.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36     defineTypeNameAndDebug(RutlandFlashBoil, 0);
  37
  38     addToRunTimeSelectionTable
  39     (
  40         evaporationModel,
  41         RutlandFlashBoil,
  42         dictionary
  43     );
  44 }
  45
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 Foam::RutlandFlashBoil::RutlandFlashBoil( const dictionary& dict)
  50 :
  51     evaporationModel(dict),
  52     evapDict_(dict.subDict(typeName + "Coeffs")),
  53     preReScFactor_(readScalar(evapDict_.lookup("preReScFactor"))),
  54     ReExponent_(readScalar(evapDict_.lookup("ReExponent"))),
  55     ScExponent_(readScalar(evapDict_.lookup("ScExponent"))),
  56     evaporationScheme_(evapDict_.lookup("evaporationScheme")),
  57     nEvapIter_(0)
  58 {
  59     if (evaporationScheme_ == "implicit")
  60     {
  61         nEvapIter_ = 2;
  62     }
  63     else if (evaporationScheme_ == "explicit")
  64     {
  65         nEvapIter_ = 1;
  66     }
  67     else
  68     {
  69         FatalError
  70             << "evaporationScheme type " << evaporationScheme_
  71             << " unknown. Use implicit or explicit." << nl
  72             << abort(FatalError);
  73     }
  74 }
  75
  76
  77 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  78
  79 Foam::RutlandFlashBoil::~RutlandFlashBoil()
  80 {}
  81
  82
  83 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  84
  85 bool Foam::RutlandFlashBoil::evaporation() const
  86 {
  87     return true;
  88 }
  89
  90
  91 // Correlation for the Sherwood Number
  92 Foam::scalar Foam::RutlandFlashBoil::Sh
  93 (
  94     const scalar ReynoldsNumber,
  95     const scalar SchmidtNumber
  96 ) const
  97 {
  98     return
  99         2.0
 100       + preReScFactor_
 101        *pow(ReynoldsNumber, ReExponent_)
 102        *pow(SchmidtNumber,ScExponent_);
 103 }
 104
 105
 106 Foam::scalar Foam::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!
 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     // proposed correction to sherwood number, implemented
 165
 166     Sherwood = 2.0 + (Sherwood - 2.0)/Fb;
 167
 168     scalar denominator =
 169         6.0*massDiffusionCoefficient*Sherwood*rhoFuelVapor*lgExpr;
 170
 171     if (denominator > SMALL)
 172     {
 173         time = max(VSMALL, liquidDensity*sqr(diameter)/denominator);
 174     }
 175
 176     return time;
 177 }
 178
 179
 180 Foam::scalar Foam::RutlandFlashBoil::boilingTime
 181 (
 182     const scalar liquidDensity,
 183     const scalar cpFuel,
 184     const scalar heatOfVapour,
 185     const scalar kappa,
 186     const scalar Nusselt,
 187     const scalar deltaTemp,
 188     const scalar diameter,
 189     const scalar liquidCore,
 190     const scalar ct,
 191     const scalar tDrop,
 192     const scalar tBoilingSurface,
 193     const scalar vapourSurfaceEnthalpy,
 194     const scalar vapourFarEnthalpy,
 195     const scalar cpGas,
 196     const scalar temperature,
 197     const scalar kLiq
 198 ) const
 199 {
 200     scalar time = GREAT;
 201
 202     // the deltaTemperature is limited to not go below .5 deg K
 203     // for numerical reasons.
 204     // This is probably not important, since it results in an upper
 205     // limit for the boiling time... which we have anyway.
 206
 207     // k set to the k value at the droplet temperature, not as in the Rutland
 208     // Paper
 209
 210     if (liquidCore > 0.5)
 211     {
 212         if (tDrop > tBoilingSurface)
 213         {
 214             //  Evaporation of the liquid sheet
 215
 216             scalar psi = 2.72;
 217             scalar kIncreased = psi*kLiq;
 218             scalar alfa = psi*kIncreased/(liquidDensity*cpFuel);
 219             scalar F = alfa*ct/sqr(0.5*diameter);
 220
 221             scalar expSum = 0.0;
 222             scalar expSumOld = expSum;
 223
 224             label Niter = 200;
 225
 226             for (label k=0; k < Niter; k++)
 227             {
 228                 expSum += exp(sqr(-k*constant::mathematical::pi*sqrt(F)/2.0));
 229                 if (mag(expSum-expSumOld)/expSum < 1.0e-3)
 230                 {
 231                     break;
 232                 }
 233                 expSumOld = expSum;
 234             }
 235         }
 236     }
 237     else
 238     {
 239         scalar dTLB = min(0.5, tDrop -  tBoilingSurface);
 240         scalar alfaS = 0.0;
 241
 242         if (dTLB >= 0.0 && dTLB < 5.0)
 243         {
 244             alfaS = 0.76*pow(dTLB, 0.26);
 245         }
 246         if (dTLB >= 5.0 &&  dTLB < 25.0)
 247         {
 248             alfaS = 0.027*pow(dTLB, 2.33);
 249         }
 250         if (dTLB >= 25.0)
 251         {
 252             alfaS = 13.8*pow(dTLB, 0.39);
 253         }
 254
 255         scalar Gf =
 256             4.0*alfaS*dTLB
 257            *constant::mathematical::pi*sqr(diameter/2.0)
 258            /heatOfVapour;
 259
 260         // calculation of the heat transfer vapourization at superheated
 261         // conditions (temperature>tBoilingSurface)
 262         scalar G = 0.0;
 263         if (temperature > tBoilingSurface)
 264         {
 265             scalar NusseltCorr = Nusselt ;
 266             scalar A =
 267                 mag((vapourFarEnthalpy-vapourSurfaceEnthalpy)/heatOfVapour);
 268
 269             // 2.0? or 1.0? try 1!
 270             scalar B =
 271                 1.0*constant::mathematical::pi*kappa/cpGas*diameter*NusseltCorr;
 272             scalar nPos = B*log(1.0 + A)/Gf + 1.0;
 273             scalar nNeg = (1.0/A)*(exp(Gf/B) - 1.0 - A) + 1.0;
 274
 275             scalar Gpos = Gf*nPos;
 276             scalar Gneg = Gf/nNeg;
 277
 278             //scalar FgPos = Gpos + Gf - B * log( 1.0 + ( 1.0 + Gf/Gpos )*A);
 279             scalar FgNeg = Gneg + Gf - B*log(1.0 + (1.0 + Gf/Gneg )*A);
 280
 281             if (FgNeg > 0.0)
 282             {
 283                 for (label j = 0; j < 20; j++)
 284                 {
 285                     Gneg = Gneg/10.0;
 286                     Gneg = max(Gneg, VSMALL);
 287                     FgNeg = Gneg + Gf - B*log(1.0 + (1.0 + Gf/Gneg)*A);
 288                     if (FgNeg < 0.0)
 289                     {
 290                         break;
 291                     }
 292                 }
 293             }
 294
 295             FgNeg = Gneg + Gf - B*log( 1.0 + ( 1.0 + Gf/Gneg)*A);
 296
 297             G = 0.5*(Gpos + Gneg);
 298             scalar Gold = -100;
 299
 300             label Niter = 200;
 301             label k=0;
 302
 303             if (FgNeg > 0.0)
 304             {
 305                 Info<< "no convergence" << endl;
 306             }
 307
 308
 309             if (FgNeg < 0.0)
 310             {
 311                 for (k=0; k<Niter; k++)
 312                 {
 313                     scalar Fg = G + Gf - B*log( 1.0 + ( 1.0 + Gf/G )*A);
 314
 315                     if (Fg > 0)
 316                     {
 317                         Gpos = G;
 318                         G = 0.5*(Gpos + Gneg);
 319                     }
 320                     else
 321                     {
 322                         Gneg = G;
 323                         G = 0.5*(Gpos + Gneg);
 324                     }
 325
 326                     Gold = G;
 327                     if (mag(G - Gold)/Gold < 1.0e-3)
 328                     {
 329                         break;
 330                     }
 331                 }
 332
 333                 if (k >= Niter - 1)
 334                 {
 335                     Info<< " No convergence for G " << endl;
 336                 }
 337             }
 338             else
 339             {
 340                 G = 0.0;
 341             }
 342         }
 343
 344         time =
 345             (constant::mathematical::pi*pow3(diameter)/6.0)
 346            *liquidDensity
 347            /(G + Gf);
 348
 349         time = max(VSMALL, time);
 350     }
 351
 352     return time;
 353 }
 354
 355
 356 // ************************************************************************* //