src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "parcel.H"
  28
  29 #include "spray.H"
  30 #include "dragModel.H"
  31 #include "evaporationModel.H"
  32 #include "heatTransferModel.H"
  33 #include "basicMultiComponentMixture.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  41
  42 void parcel::setRelaxationTimes
  43 (
  44     label celli,
  45     scalar& tauMomentum,
  46     scalarField& tauEvaporation,
  47     scalar& tauHeatTransfer,
  48     scalarField& tauBoiling,
  49     const spray& sDB,
  50     const scalar rho,
  51     const vector& Up,
  52     const scalar temperature,
  53     const scalar pressure,
  54     const scalarField& Yfg,
  55     const scalarField& m0,
  56     const scalar dt
  57 )
  58 {
  59
  60     const liquidMixture& fuels = sDB.fuels();
  61
  62     scalar mCell = rho*sDB.mesh().V()[cell()];
  63     scalarField mfg(Yfg*mCell);
  64
  65     label Ns = sDB.composition().Y().size();
  66     label Nf = fuels.components().size();
  67
  68     // Tf is based on the 1/3 rule
  69     scalar Tf  = T() + (temperature - T())/3.0;
  70
  71     // calculate mixture properties
  72     scalar W = 0.0;
  73     scalar kMixture = 0.0;
  74     scalar cpMixture = 0.0;
  75     scalar muf = 0.0;
  76
  77     for(label i=0; i<Ns; i++)
  78     {
  79         scalar Y = sDB.composition().Y()[i][celli];
  80         W += Y/sDB.gasProperties()[i].W();
  81         // Using mass-fractions to average...
  82         kMixture += Y*sDB.gasProperties()[i].kappa(Tf);
  83         cpMixture += Y*sDB.gasProperties()[i].Cp(Tf);
  84         muf += Y*sDB.gasProperties()[i].mu(Tf);
  85     }
  86     W = 1.0/W;
  87
  88     scalarField Xf(Nf, 0.0);
  89     scalarField Yf(Nf, 0.0);
  90     scalarField psat(Nf, 0.0);
  91     scalarField msat(Nf, 0.0);
  92
  93     for(label i=0; i<Nf; i++)
  94     {
  95         label j = sDB.liquidToGasIndex()[i];
  96         scalar Y = sDB.composition().Y()[j][celli];
  97         scalar Wi = sDB.gasProperties()[j].W();
  98         Yf[i] = Y;
  99         Xf[i] = Y*W/Wi;
 100         psat[i] = fuels.properties()[i].pv(pressure, temperature);
 101         msat[i] = min(1.0, psat[i]/pressure)*Wi/W;
 102     }
 103
 104     scalar nuf = muf/rho;
 105
 106     scalar liquidDensity = fuels.rho(pressure, T(), X());
 107     scalar liquidcL = fuels.cp(pressure, T(), X());
 108     scalar heatOfVapour = fuels.hl(pressure, T(), X());
 109
 110     // calculate the partial rho of the fuel vapour
 111     // alternative is to use the mass fraction
 112     // however, if rhoFuelVap is small (zero)
 113     // d(mass)/dt = 0 => no evaporation... hmmm... is that good? NO!
 114
 115     // Assume equilibrium at drop-surface => pressure @ surface
 116     // = vapour pressure to calculate fuel-vapour density @ surface
 117     scalar pressureAtSurface = fuels.pv(pressure, T(), X());
 118     scalar rhoFuelVap = pressureAtSurface*fuels.W(X())/(specie::RR*Tf);
 119
 120     scalarField Xs(sDB.fuels().Xs(pressure, temperature, T(), Xf, X()));
 121     scalarField Ys(Nf, 0.0);
 122     scalar Wliq = 0.0;
 123
 124     for(label i=0; i<Nf; i++)
 125     {
 126         label j = sDB.liquidToGasIndex()[i];
 127         scalar Wi = sDB.gasProperties()[j].W();
 128         Wliq += Xs[i]*Wi;
 129     }
 130
 131     for(label i=0; i<Nf; i++)
 132     {
 133         label j = sDB.liquidToGasIndex()[i];
 134         scalar Wi = sDB.gasProperties()[j].W();
 135         Ys[i] = Xs[i]*Wi/Wliq;
 136     }
 137
 138     scalar Reynolds = Re(Up, nuf);
 139     scalar Prandtl = Pr(cpMixture, muf, kMixture);
 140
 141     // calculate the characteritic times
 142
 143     if(liquidCore_> 0.5)
 144     {
 145 //      no drag for parcels in the liquid core..
 146         tauMomentum = GREAT;
 147     }
 148     else
 149     {
 150         tauMomentum = sDB.drag().relaxationTime
 151         (
 152             Urel(Up),
 153             d(),
 154             rho,
 155             liquidDensity,
 156             nuf,
 157             dev()
 158         );
 159     }
 160
 161     // store the relaxationTime since it is needed in some breakup models.
 162     tMom_ = tauMomentum;
 163
 164     tauHeatTransfer = sDB.heatTransfer().relaxationTime
 165     (
 166         liquidDensity,
 167         d(),
 168         liquidcL,
 169         kMixture,
 170         Reynolds,
 171         Prandtl
 172     );
 173
 174     // evaporation-properties are evaluated at averaged temperature
 175     // set the boiling conditions true if pressure @ surface is 99.9%
 176     // of the pressure
 177     // this is mainly to put a limit on the evaporation time,
 178     // since tauEvaporation is very very small close to the boiling point.
 179
 180     for(label i=0; i<Nf; i++)
 181     {
 182         scalar Td = min(T(), 0.999*fuels.properties()[i].Tc());
 183         bool boiling = fuels.properties()[i].pv(pressure, Td)
 184             >= 0.999*pressure;
 185
 186         scalar Di = fuels.properties()[i].D(pressure, Td);
 187         scalar Schmidt = Sc(nuf, Di);
 188
 189         scalar partialPressure = Xf[i]*pressure;
 190
 191 //      saturated vapour
 192         if(partialPressure > psat[i])
 193         {
 194             tauEvaporation[i] = GREAT;
 195         }
 196 //      not saturated vapour
 197         else
 198         {
 199             if (!boiling)
 200             {
 201                 // For saturation evaporation, only use 99.99% for
 202                 // numerical robustness
 203                 scalar dm = max(SMALL, 0.9999*msat[i] - mfg[i]);
 204
 205                 tauEvaporation[i] = sDB.evaporation().relaxationTime
 206                 (
 207                     d(),
 208                     fuels.properties()[i].rho(pressure, Td),
 209                     rhoFuelVap,
 210                     Di,
 211                     Reynolds,
 212                     Schmidt,
 213                     Xs[i],
 214                     Xf[i],
 215                     m0[i],
 216                     dm,
 217                     dt
 218                 );
 219             }
 220             else
 221             {
 222                 scalar Nusselt =
 223                     sDB.heatTransfer().Nu(Reynolds, Prandtl);
 224
 225 //              calculating the boiling temperature of the liquid at ambient pressure
 226                 scalar tBoilingSurface = Td;
 227
 228                 label Niter = 0;
 229                 scalar deltaT = 10.0;
 230                 scalar dp0 = fuels.properties()[i].pv(pressure, tBoilingSurface) - pressure;
 231                 while ((Niter < 200) && (mag(deltaT) > 1.0e-3))
 232                 {
 233                     Niter++;
 234                     scalar pBoil = fuels.properties()[i].pv(pressure, tBoilingSurface);
 235                     scalar dp = pBoil - pressure;
 236                     if ( (dp > 0.0) && (dp0 > 0.0) )
 237                     {
 238                         tBoilingSurface -= deltaT;
 239                     }
 240                     else
 241                     {
 242                         if ( (dp < 0.0) && (dp0 < 0.0) )
 243                         {
 244                             tBoilingSurface += deltaT;
 245                         }
 246                         else
 247                         {
 248                             deltaT *= 0.5;
 249                             if ( (dp > 0.0) && (dp0 < 0.0) )
 250                             {
 251                                 tBoilingSurface -= deltaT;
 252                             }
 253                             else
 254                             {
 255                                 tBoilingSurface += deltaT;
 256                             }
 257                         }
 258                     }
 259                     dp0 = dp;
 260                 }
 261
 262                 scalar vapourSurfaceEnthalpy = 0.0;
 263                 scalar vapourFarEnthalpy = 0.0;
 264
 265                 for(label k = 0; k < sDB.gasProperties().size(); k++)
 266                 {
 267                     vapourSurfaceEnthalpy += sDB.composition().Y()[k][celli]*sDB.gasProperties()[k].H(tBoilingSurface);
 268                     vapourFarEnthalpy += sDB.composition().Y()[k][celli]*sDB.gasProperties()[k].H(temperature);
 269                 }
 270
 271                 scalar kLiquid = fuels.properties()[i].K(pressure, 0.5*(tBoilingSurface+T()));
 272
 273                 tauBoiling[i] = sDB.evaporation().boilingTime
 274                 (
 275                     fuels.properties()[i].rho(pressure, Td),
 276                     fuels.properties()[i].cp(pressure, Td),
 277                     heatOfVapour,
 278                     kMixture,
 279                     Nusselt,
 280                     temperature - T(),
 281                     d(),
 282                     liquidCore(),
 283                     sDB.runTime().value() - ct(),
 284                     Td,
 285                     tBoilingSurface,
 286                     vapourSurfaceEnthalpy,
 287                     vapourFarEnthalpy,
 288                     cpMixture,
 289                     temperature,
 290                     kLiquid
 291                 );
 292             }
 293
 294         }
 295     }
 296 }
 297
 298
 299 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 300
 301 } // End namespace Foam
 302
 303 // ************************************************************************* //