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