src/lagrangian/dieselSpray/spraySubModels/breakupModel/TAB/TAB.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 "TAB.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "mathematicalConstants.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 defineTypeNameAndDebug(TAB, 0);
  39
  40 addToRunTimeSelectionTable
  41 (
  42     breakupModel,
  43     TAB,
  44     dictionary
  45 );
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 // Construct from components
  50 TAB::TAB
  51 (
  52     const dictionary& dict,
  53     spray& sm
  54 )
  55 :
  56     breakupModel(dict, sm),
  57     coeffsDict_(dict.subDict(typeName + "Coeffs")),
  58     Cmu_(readScalar(coeffsDict_.lookup("Cmu"))),
  59     Comega_(readScalar(coeffsDict_.lookup("Comega"))),
  60     WeCrit_(readScalar(coeffsDict_.lookup("WeCrit")))
  61 {
  62
  63     // calculate the inverse function of the Rossin-Rammler Distribution
  64     const scalar xx0 = 12.0;
  65     const scalar rrd100 = 1.0/(1.0-exp(-xx0)*(1+xx0+pow(xx0, 2)/2+pow(xx0, 3)/6));
  66
  67     for(label n=0; n<100; n++)
  68     {
  69         scalar xx = 0.12*(n+1);
  70         rrd_[n] = (1-exp(-xx)*(1 + xx + pow(xx, 2)/2 + pow(xx, 3)/6))*rrd100;
  71     }
  72 }
  73
  74
  75 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  76
  77 TAB::~TAB()
  78 {}
  79
  80
  81 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  82
  83 void TAB::breakupParcel
  84 (
  85     parcel& p,
  86     const scalar deltaT,
  87     const vector& Ug,
  88     const liquidMixture& fuels
  89 ) const
  90 {
  91
  92     scalar T = p.T();
  93     scalar pc = spray_.p()[p.cell()];
  94     scalar r = 0.5*p.d();
  95     scalar r2 = r*r;
  96     scalar r3 = r*r2;
  97
  98     scalar rho = fuels.rho(pc, T, p.X());
  99     scalar sigma = fuels.sigma(pc, T, p.X());
 100     scalar mu = fuels.mu(pc, T, p.X());
 101
 102     // inverse of characteristic viscous damping time
 103     scalar rtd = 0.5*Cmu_*mu/(rho*r2);
 104
 105     // oscillation frequency (squared)
 106     scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
 107
 108     if (omega2 > 0)
 109     {
 110         scalar omega = sqrt(omega2);
 111         scalar rhog = spray_.rho()[p.cell()];
 112         scalar We = p.We(Ug, rhog, sigma);
 113         scalar Wetmp = We/WeCrit_;
 114
 115         scalar y1 = p.dev() - Wetmp;
 116         scalar y2 = p.ddev()/omega;
 117
 118         scalar a = sqrt(y1*y1 + y2*y2);
 119
 120         // scotty we may have break-up
 121         if (a+Wetmp > 1.0)
 122         {
 123             scalar phic = y1/a;
 124
 125             // constrain phic within -1 to 1
 126             phic = max(min(phic, 1), -1);
 127
 128             scalar phit = acos(phic);
 129             scalar phi = phit;
 130             scalar quad = -y2/a;
 131             if (quad < 0)
 132             {
 133                 phi = 2*mathematicalConstant::pi - phit;
 134             }
 135
 136             scalar tb = 0;
 137
 138             if (mag(p.dev()) < 1.0)
 139             {
 140                 scalar coste = 1.0;
 141                 if
 142                 (
 143                     (Wetmp - a < -1)
 144                  && (p.ddev() < 0)
 145                 )
 146                 {
 147                     coste = -1.0;
 148                 }
 149
 150                 scalar theta = acos((coste-Wetmp)/a);
 151
 152                 if (theta < phi)
 153                 {
 154                     if (2*mathematicalConstant::pi-theta >= phi)
 155                     {
 156                         theta = -theta;
 157                     }
 158                     theta += 2*mathematicalConstant::pi;
 159                 }
 160                 tb = (theta-phi)/omega;
 161
 162                 // breakup occurs
 163                 if (deltaT > tb)
 164                 {
 165                     p.dev() = 1.0;
 166                     p.ddev() = -a*omega*sin(omega*tb + phi);
 167                 }
 168
 169             }
 170
 171             // update droplet size
 172             if (deltaT > tb)
 173             {
 174                 scalar rs = r/
 175                 (
 176                     1
 177                   + (4.0/3.0)*pow(p.dev(), 2)
 178                   + rho*r3/(8*sigma)*pow(p.ddev(), 2)
 179                 );
 180
 181                 label n = 0;
 182                 bool found = false;
 183                 scalar random = rndGen_.scalar01();
 184                 while (!found && (n<99))
 185                 {
 186                     if (rrd_[n]>random)
 187                     {
 188                         found = true;
 189                     }
 190                     n++;
 191
 192                 }
 193                 scalar rNew = 0.04*n*rs;
 194                 if (rNew < r)
 195                 {
 196                     p.d() = 2*rNew;
 197                     p.dev() = 0;
 198                     p.ddev() = 0;
 199                 }
 200
 201             }
 202
 203         }
 204
 205     }
 206     else
 207     {
 208         // reset droplet distortion parameters
 209         p.dev() = 0;
 210         p.ddev() = 0;
 211     }
 212
 213 }
 214
 215 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 216
 217 } // End namespace Foam
 218
 219 // ************************************************************************* //