src/lagrangian/dieselSpray/spraySubModels/breakupModel/ETAB/ETAB.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 "ETAB.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "mathematicalConstants.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(ETAB, 0);
  38
  39 addToRunTimeSelectionTable
  40 (
  41     breakupModel,
  42     ETAB,
  43     dictionary
  44 );
  45
  46 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  47
  48 // Construct from components
  49 ETAB::ETAB
  50 (
  51     const dictionary& dict,
  52     spray& sm
  53 )
  54 :
  55     breakupModel(dict, sm),
  56     coeffsDict_(dict.subDict(typeName + "Coeffs")),
  57     Cmu_(readScalar(coeffsDict_.lookup("Cmu"))),
  58     Comega_(readScalar(coeffsDict_.lookup("Comega"))),
  59     k1_(readScalar(coeffsDict_.lookup("k1"))),
  60     k2_(readScalar(coeffsDict_.lookup("k2"))),
  61     WeCrit_(readScalar(coeffsDict_.lookup("WeCrit"))),
  62     WeTransition_(readScalar(coeffsDict_.lookup("WeTransition"))),
  63     AWe_(0.0)
  64 {
  65     scalar k21 = k2_/k1_;
  66     AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow(WeTransition_, 4.0);
  67 }
  68
  69
  70 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  71
  72 ETAB::~ETAB()
  73 {}
  74
  75
  76 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  77
  78 void ETAB::breakupParcel
  79 (
  80     parcel& p,
  81     const scalar deltaT,
  82     const vector& Ug,
  83     const liquidMixture& fuels
  84 ) const
  85 {
  86
  87     scalar T  = p.T();
  88     scalar pc  = spray_.p()[p.cell()];
  89     scalar r  = 0.5*p.d();
  90     scalar r2 = r*r;
  91     scalar r3 = r*r2;
  92
  93     scalar rho   = fuels.rho(pc, T, p.X());
  94     scalar sigma = fuels.sigma(pc, T, p.X());
  95     scalar mu    = fuels.mu(pc, T, p.X());
  96
  97     // inverse of characteristic viscous damping time
  98     scalar rtd = 0.5*Cmu_*mu/(rho*r2);
  99
 100     // oscillation frequency (squared)
 101     scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
 102
 103     if (omega2 > 0)
 104     {
 105         scalar omega = sqrt(omega2);
 106         scalar romega = 1.0/omega;
 107
 108         scalar rhog = spray_.rho()[p.cell()];
 109         scalar We = p.We(Ug, rhog, sigma);
 110         scalar Wetmp = We/WeCrit_;
 111
 112         scalar y1 = p.dev() - Wetmp;
 113         scalar y2 = p.ddev()*romega;
 114
 115         scalar a = sqrt(y1*y1 + y2*y2);
 116
 117         // scotty we may have break-up
 118         if (a+Wetmp > 1.0)
 119         {
 120             scalar phic = y1/a;
 121
 122             // constrain phic within -1 to 1
 123             phic = max(min(phic, 1), -1);
 124
 125             scalar phit = acos(phic);
 126             scalar phi = phit;
 127             scalar quad = -y2/a;
 128             if (quad < 0)
 129             {
 130                 phi = 2*mathematicalConstant::pi - phit;
 131             }
 132
 133             scalar tb = 0;
 134
 135             if (mag(p.dev()) < 1.0)
 136             {
 137                 scalar theta = acos((1.0 - Wetmp)/a);
 138
 139                 if (theta < phi)
 140                 {
 141                     if (2*mathematicalConstant::pi-theta >= phi)
 142                     {
 143                         theta = -theta;
 144                     }
 145                     theta += 2*mathematicalConstant::pi;
 146                 }
 147                 tb = (theta-phi)*romega;
 148
 149                 // breakup occurs
 150                 if (deltaT > tb)
 151                 {
 152                     p.dev() = 1.0;
 153                     p.ddev() = -a*omega*sin(omega*tb + phi);
 154                 }
 155             }
 156
 157             // update droplet size
 158             if (deltaT > tb)
 159             {
 160                 scalar sqrtWe = AWe_*pow(We, 4.0) + 1.0;
 161                 scalar Kbr = k1_*omega*sqrtWe;
 162
 163                 if (We > WeTransition_)
 164                 {
 165                     sqrtWe = sqrt(We);
 166                     Kbr =k2_*omega*sqrtWe;
 167                 }
 168
 169                 scalar rWetmp = 1.0/Wetmp;
 170                 scalar cosdtbu = max(-1.0, min(1.0, 1.0-rWetmp));
 171                 scalar dtbu = romega*acos(cosdtbu);
 172                 scalar decay = exp(-Kbr*dtbu);
 173
 174                 scalar rNew = decay*r;
 175                 if (rNew < r)
 176                 {
 177                     p.d() = 2*rNew;
 178                     p.dev() = 0;
 179                     p.ddev() = 0;
 180                 }
 181             }
 182         }
 183     }
 184     else
 185     {
 186         // reset droplet distortion parameters
 187         p.dev() = 0;
 188         p.ddev() = 0;
 189     }
 190 }
 191
 192 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 193
 194 } // End namespace Foam
 195
 196 // ************************************************************************* //