FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

TAB.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "TAB.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/mathematicalConstants.H>
00029 
00030 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00031 
00032 namespace Foam
00033 {
00034 
00035 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00036 
00037 defineTypeNameAndDebug(TAB, 0);
00038 
00039 addToRunTimeSelectionTable
00040 (
00041     breakupModel,
00042     TAB,
00043     dictionary
00044 );
00045 
00046 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00047 
00048 // Construct from components
00049 TAB::TAB
00050 (
00051     const dictionary& dict,
00052     spray& sm
00053 )
00054 :
00055     breakupModel(dict, sm),
00056     coeffsDict_(dict.subDict(typeName + "Coeffs")),
00057     Cmu_(readScalar(coeffsDict_.lookup("Cmu"))),
00058     Comega_(readScalar(coeffsDict_.lookup("Comega"))),
00059     WeCrit_(readScalar(coeffsDict_.lookup("WeCrit")))
00060 {
00061 
00062     // calculate the inverse function of the Rossin-Rammler Distribution
00063     const scalar xx0 = 12.0;
00064     const scalar rrd100 = 1.0/(1.0-exp(-xx0)*(1+xx0+pow(xx0, 2)/2+pow(xx0, 3)/6));
00065 
00066     for(label n=0; n<100; n++)
00067     {
00068         scalar xx = 0.12*(n+1);
00069         rrd_[n] = (1-exp(-xx)*(1 + xx + pow(xx, 2)/2 + pow(xx, 3)/6))*rrd100;
00070     }
00071 }
00072 
00073 
00074 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00075 
00076 TAB::~TAB()
00077 {}
00078 
00079 
00080 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00081 
00082 void TAB::breakupParcel
00083 (
00084     parcel& p,
00085     const scalar deltaT,
00086     const vector& Ug,
00087     const liquidMixture& fuels
00088 ) const
00089 {
00090 
00091     scalar T = p.T();
00092     scalar pc = spray_.p()[p.cell()];
00093     scalar r = 0.5*p.d();
00094     scalar r2 = r*r;
00095     scalar r3 = r*r2;
00096 
00097     scalar rho = fuels.rho(pc, T, p.X());
00098     scalar sigma = fuels.sigma(pc, T, p.X());
00099     scalar mu = fuels.mu(pc, T, p.X());
00100 
00101     // inverse of characteristic viscous damping time
00102     scalar rtd = 0.5*Cmu_*mu/(rho*r2);
00103     
00104     // oscillation frequency (squared)
00105     scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
00106 
00107     if (omega2 > 0)
00108     {
00109         scalar omega = sqrt(omega2);
00110         scalar rhog = spray_.rho()[p.cell()];
00111         scalar We = p.We(Ug, rhog, sigma);
00112         scalar Wetmp = We/WeCrit_;
00113 
00114         scalar y1 = p.dev() - Wetmp;
00115         scalar y2 = p.ddev()/omega;
00116 
00117         scalar a = sqrt(y1*y1 + y2*y2);
00118 
00119         // scotty we may have break-up
00120         if (a+Wetmp > 1.0)
00121         {
00122             scalar phic = y1/a;
00123 
00124             // constrain phic within -1 to 1
00125             phic = max(min(phic, 1), -1);
00126 
00127             scalar phit = acos(phic);
00128             scalar phi = phit;
00129             scalar quad = -y2/a;
00130             if (quad < 0)
00131             {
00132                 phi = 2*mathematicalConstant::pi - phit;
00133             }
00134             
00135             scalar tb = 0;
00136             
00137             if (mag(p.dev()) < 1.0)
00138             {
00139                 scalar coste = 1.0;
00140                 if
00141                 (
00142                     (Wetmp - a < -1)
00143                  && (p.ddev() < 0)
00144                 )
00145                 {
00146                     coste = -1.0;
00147                 }
00148                 
00149                 scalar theta = acos((coste-Wetmp)/a);
00150                 
00151                 if (theta < phi)
00152                 {
00153                     if (2*mathematicalConstant::pi-theta >= phi)
00154                     {
00155                         theta = -theta;
00156                     }
00157                     theta += 2*mathematicalConstant::pi;
00158                 }
00159                 tb = (theta-phi)/omega;
00160 
00161                 // breakup occurs
00162                 if (deltaT > tb)
00163                 {
00164                     p.dev() = 1.0;
00165                     p.ddev() = -a*omega*sin(omega*tb + phi);
00166                 }
00167 
00168             }
00169 
00170             // update droplet size
00171             if (deltaT > tb)
00172             {
00173                 scalar rs = r/
00174                 (
00175                     1 
00176                   + (4.0/3.0)*pow(p.dev(), 2)
00177                   + rho*r3/(8*sigma)*pow(p.ddev(), 2)
00178                 );
00179                 
00180                 label n = 0;
00181                 bool found = false;
00182                 scalar random = rndGen_.scalar01();
00183                 while (!found && (n<99))
00184                 {
00185                     if (rrd_[n]>random)
00186                     {
00187                         found = true;
00188                     }
00189                     n++;
00190 
00191                 }
00192                 scalar rNew = 0.04*n*rs;
00193                 if (rNew < r)
00194                 {
00195                     p.d() = 2*rNew;
00196                     p.dev() = 0;
00197                     p.ddev() = 0;
00198                 }
00199 
00200             }
00201 
00202         }
00203        
00204     }
00205     else
00206     {
00207         // reset droplet distortion parameters
00208         p.dev() = 0;
00209         p.ddev() = 0;
00210     }
00211 
00212 }
00213 
00214 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00215 
00216 } // End namespace Foam
00217 
00218 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines