src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.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 "definedPressureSwirl.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "mathematicalConstants.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(definedPressureSwirlInjector, 0);
  38
  39 addToRunTimeSelectionTable
  40 (
  41     injectorModel,
  42     definedPressureSwirlInjector,
  43     dictionary
  44 );
  45
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 // Construct from components
  50 definedPressureSwirlInjector::definedPressureSwirlInjector
  51 (
  52     const dictionary& dict,
  53     spray& sm
  54 )
  55 :
  56     injectorModel(dict, sm),
  57     definedPressureSwirlInjectorDict_(dict.subDict(typeName + "Coeffs")),
  58
  59     coneAngle_(definedPressureSwirlInjectorDict_.lookup("ConeAngle")),
  60     coneInterval_(definedPressureSwirlInjectorDict_.lookup("ConeInterval")),
  61     maxKv_(definedPressureSwirlInjectorDict_.lookup("maxKv")),
  62
  63     angle_(0.0)
  64 {
  65
  66     scalar referencePressure = sm.p().average().value();
  67
  68     // correct velocityProfile
  69     forAll(sm.injectors(), i)
  70     {
  71         sm.injectors()[i].properties()->correctProfiles(sm.fuels(), referencePressure);
  72     }
  73
  74 }
  75
  76
  77 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  78
  79 definedPressureSwirlInjector::~definedPressureSwirlInjector()
  80 {}
  81
  82
  83 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  84
  85 scalar definedPressureSwirlInjector::d0
  86 (
  87     const label n,
  88     const scalar t
  89 ) const
  90 {
  91     const injectorType& it = injectors_[n].properties();
  92
  93     scalar c = rndGen_.scalar01();
  94     scalar coneAngle = it.getTableValue(coneAngle_, t);
  95     scalar coneInterval = it.getTableValue(coneInterval_, t);
  96     angle_ = coneAngle ;
  97
  98 //  modifications to take account of flash boiling....
  99
 100     const liquidMixture& fuels = sm_.fuels();
 101     scalar chi = 0.0;
 102     scalar Tinj = it.T(t);
 103     label Nf = fuels.components().size();
 104     scalar temperature = sm_.ambientTemperature();
 105     scalar pressure = sm_.ambientPressure();
 106
 107
 108     for(label i = 0; i < Nf ; i++)
 109     {
 110
 111         if(fuels.properties()[i].pv(sm_.ambientPressure(), Tinj) >= 0.999*sm_.ambientPressure())
 112         {
 113
 114 //          The fuel is boiling.....
 115 //          Calculation of the boiling temperature
 116
 117             scalar tBoilingSurface = Tinj ;
 118
 119             label Niter = 200;
 120
 121             for(label k=0; k< Niter ; k++)
 122             {
 123
 124                 scalar pBoil = fuels.properties()[i].pv(pressure, tBoilingSurface);
 125
 126                 if(pBoil > pressure)
 127                 {
 128                     tBoilingSurface = tBoilingSurface - (Tinj-temperature)/Niter;
 129                 }
 130                 else
 131                 {
 132                     break;
 133                 }
 134
 135             }
 136
 137             scalar hl = fuels.properties()[i].hl(sm_.ambientPressure(), tBoilingSurface);
 138             scalar iTp = fuels.properties()[i].h(sm_.ambientPressure(), Tinj) - sm_.ambientPressure()/fuels.properties()[i].rho(sm_.ambientPressure(), Tinj);
 139             scalar iTb = fuels.properties()[i].h(sm_.ambientPressure(), tBoilingSurface) - sm_.ambientPressure()/fuels.properties()[i].rho(sm_.ambientPressure(), tBoilingSurface);
 140
 141             chi += it.X()[i]*(iTp-iTb)/hl;
 142
 143         }
 144     }
 145
 146     //  bounding chi
 147
 148     chi = max(chi, 0.0);
 149     chi = min(chi, 1.0);
 150
 151     angle_ = angle_ + (144.0 - angle_) * sqr(chi) + 2.0 * coneInterval * (0.5 - c);
 152
 153 //  end modifications
 154
 155     angle_ *= mathematicalConstant::pi/360.0;
 156
 157     scalar injectedMassFlow = it.massFlowRate(t);
 158
 159     scalar cosAngle = cos(angle_);
 160
 161     scalar rhoFuel = sm_.fuels().rho(sm_.ambientPressure(), it.T(t), it.X());
 162     scalar injectorDiameter = it.d();
 163
 164     scalar deltaPressure = deltaPressureInj(t,n);
 165
 166     scalar kV = kv(n, injectedMassFlow, deltaPressure, t);
 167
 168     scalar v = kV * sqrt(2.0*deltaPressure/rhoFuel);
 169
 170     u_ = v * cosAngle;
 171
 172     scalar A = injectedMassFlow/(mathematicalConstant::pi*rhoFuel*u_);
 173
 174 /*
 175
 176     TL
 177     The formula for the sheet tickness proposed by the authors is,
 178     in my opinion, "strange".....
 179     I modified it multiplying the sheet tickness for the cone angle cosinus.
 180
 181 */
 182
 183     scalar angleT = angle_;
 184     return (injectorDiameter-sqrt(pow(injectorDiameter,2.0)-4.0*A))*cos(angleT)/2.0;
 185
 186 //  original implementation
 187
 188 /*
 189     return (injectorDiameter-sqrt(pow(injectorDiameter,2)-4.0*A))/2.0;
 190 */
 191
 192
 193 }
 194
 195 vector definedPressureSwirlInjector::direction
 196 (
 197     const label n,
 198     const label hole,
 199     const scalar time,
 200     const scalar d
 201 ) const
 202 {
 203
 204     scalar alpha = sin(angle_);
 205     scalar dcorr = cos(angle_);
 206     scalar beta = 2.0*mathematicalConstant::pi*rndGen_.scalar01();
 207
 208     // randomly distributed vector normal to the injection vector
 209     vector normal = vector::zero;
 210
 211     if (sm_.twoD())
 212     {
 213         scalar reduce = 0.01;
 214         // correct beta if this is a 2D run
 215         // map it onto the 'angleOfWedge'
 216
 217         beta *= (1.0-2.0*reduce)*sm_.angleOfWedge()/(2.0*mathematicalConstant::pi);
 218         beta += reduce*sm_.angleOfWedge();
 219         normal = alpha*
 220         (
 221             sm_.axisOfWedge()*cos(beta) +
 222             sm_.axisOfWedgeNormal()*sin(beta)
 223         );
 224     }
 225     else
 226     {
 227         normal = alpha*
 228         (
 229             injectors_[n].properties()->tan1(hole)*cos(beta) +
 230             injectors_[n].properties()->tan2(hole)*sin(beta)
 231         );
 232     }
 233
 234     // set the direction of injection by adding the normal vector
 235     vector dir = dcorr*injectors_[n].properties()->direction(hole, time) + normal;
 236     dir /= mag(dir);
 237
 238     return dir;
 239 }
 240
 241
 242 scalar definedPressureSwirlInjector::velocity
 243 (
 244     const label i,
 245     const scalar time
 246 ) const
 247 {
 248     return u_*sqrt(1.0 + pow(tan(angle_),2.0));
 249 }
 250
 251 scalar definedPressureSwirlInjector::averageVelocity
 252 (
 253     const label i
 254 ) const
 255 {
 256
 257     const injectorType& it = sm_.injectors()[i].properties();
 258
 259     scalar dt = it.teoi() - it.tsoi();
 260
 261     scalar injectedMassFlow = it.mass()/(it.teoi()-it.tsoi());
 262
 263     scalar injectionPressure = averagePressure(i);
 264
 265     scalar Tav = it.integrateTable(it.T())/dt;
 266     scalar rhoFuel = sm_.fuels().rho(sm_.ambientPressure(), Tav, it.X());
 267
 268     scalar kV = kv(i, injectedMassFlow, injectionPressure, 0);
 269
 270     return  kV*sqrt(2.0*(injectionPressure-sm_.ambientPressure())/rhoFuel);
 271
 272 }
 273
 274
 275 scalar definedPressureSwirlInjector::kv
 276 (
 277     const label inj,
 278     const scalar massFlow,
 279     const scalar dPressure,
 280     const scalar t
 281 ) const
 282 {
 283
 284     const injectorType& it = injectors_[inj].properties();
 285
 286     scalar coneAngle = it.getTableValue(coneAngle_, t);
 287
 288     coneAngle *= mathematicalConstant::pi/360.0;
 289
 290     scalar cosAngle = cos(coneAngle);
 291     scalar Tav = it.integrateTable(it.T())/(it.teoi()-it.tsoi());
 292
 293     scalar rhoFuel = sm_.fuels().rho(sm_.ambientPressure(), Tav, it.X());
 294     scalar injectorDiameter = it.d();
 295
 296     scalar kv = max
 297     (
 298         it.getTableValue(maxKv_, t),
 299         4.0*massFlow
 300         *
 301         sqrt(rhoFuel/2.0/dPressure)
 302         /
 303         (mathematicalConstant::pi*pow(injectorDiameter, 2.0)*rhoFuel*cosAngle)
 304     );
 305
 306     return min(1.0,kv);
 307 }
 308
 309
 310
 311
 312 scalar definedPressureSwirlInjector::deltaPressureInj(const scalar time, const label inj) const
 313 {
 314     return injectors_[inj].properties()->injectionPressure(time) - sm_.ambientPressure();
 315 }
 316
 317 scalar definedPressureSwirlInjector::averagePressure(const label inj) const
 318 {
 319
 320     const injectorType& it = sm_.injectors()[inj].properties();
 321
 322     scalar dt = it.teoi() - it.tsoi();
 323     return it.integrateTable(it.injectionPressureProfile())/dt;
 324 }
 325
 326 } // End namespace Foam
 327
 328 // ************************************************************************* //