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