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

ETAB.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 "ETAB.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(ETAB, 0);
00038 
00039 addToRunTimeSelectionTable
00040 (
00041     breakupModel,
00042     ETAB,
00043     dictionary
00044 );
00045 
00046 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00047 
00048 // Construct from components
00049 ETAB::ETAB
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     k1_(readScalar(coeffsDict_.lookup("k1"))),
00060     k2_(readScalar(coeffsDict_.lookup("k2"))),
00061     WeCrit_(readScalar(coeffsDict_.lookup("WeCrit"))),
00062     WeTransition_(readScalar(coeffsDict_.lookup("WeTransition"))),
00063     AWe_(0.0)
00064 {
00065     scalar k21 = k2_/k1_;
00066     AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow(WeTransition_, 4.0);
00067 }
00068 
00069 
00070 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00071 
00072 ETAB::~ETAB()
00073 {}
00074 
00075 
00076 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00077 
00078 void ETAB::breakupParcel
00079 (
00080     parcel& p,
00081     const scalar deltaT,
00082     const vector& Ug,
00083     const liquidMixture& fuels
00084 ) const
00085 {
00086 
00087     scalar T  = p.T();
00088     scalar pc  = spray_.p()[p.cell()];
00089     scalar r  = 0.5*p.d();
00090     scalar r2 = r*r;
00091     scalar r3 = r*r2;
00092 
00093     scalar rho   = fuels.rho(pc, T, p.X());
00094     scalar sigma = fuels.sigma(pc, T, p.X());
00095     scalar mu    = fuels.mu(pc, T, p.X());
00096 
00097     // inverse of characteristic viscous damping time
00098     scalar rtd = 0.5*Cmu_*mu/(rho*r2);
00099 
00100     // oscillation frequency (squared)
00101     scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
00102 
00103     if (omega2 > 0)
00104     {
00105         scalar omega = sqrt(omega2);
00106         scalar romega = 1.0/omega;
00107 
00108         scalar rhog = spray_.rho()[p.cell()];
00109         scalar We = p.We(Ug, rhog, sigma);
00110         scalar Wetmp = We/WeCrit_;
00111 
00112         scalar y1 = p.dev() - Wetmp;
00113         scalar y2 = p.ddev()*romega;
00114 
00115         scalar a = sqrt(y1*y1 + y2*y2);
00116 
00117         // scotty we may have break-up
00118         if (a+Wetmp > 1.0)
00119         {
00120             scalar phic = y1/a;
00121 
00122             // constrain phic within -1 to 1
00123             phic = max(min(phic, 1), -1);
00124 
00125             scalar phit = acos(phic);
00126             scalar phi = phit;
00127             scalar quad = -y2/a;
00128             if (quad < 0)
00129             {
00130                 phi = 2*mathematicalConstant::pi - phit;
00131             }
00132 
00133             scalar tb = 0;
00134 
00135             if (mag(p.dev()) < 1.0)
00136             {
00137                 scalar theta = acos((1.0 - Wetmp)/a);
00138 
00139                 if (theta < phi)
00140                 {
00141                     if (2*mathematicalConstant::pi-theta >= phi)
00142                     {
00143                         theta = -theta;
00144                     }
00145                     theta += 2*mathematicalConstant::pi;
00146                 }
00147                 tb = (theta-phi)*romega;
00148 
00149                 // breakup occurs
00150                 if (deltaT > tb)
00151                 {
00152                     p.dev() = 1.0;
00153                     p.ddev() = -a*omega*sin(omega*tb + phi);
00154                 }
00155             }
00156 
00157             // update droplet size
00158             if (deltaT > tb)
00159             {
00160                 scalar sqrtWe = AWe_*pow(We, 4.0) + 1.0;
00161                 scalar Kbr = k1_*omega*sqrtWe;
00162 
00163                 if (We > WeTransition_)
00164                 {
00165                     sqrtWe = sqrt(We);
00166                     Kbr =k2_*omega*sqrtWe;
00167                 }
00168 
00169                 scalar rWetmp = 1.0/Wetmp;
00170                 scalar cosdtbu = max(-1.0, min(1.0, 1.0-rWetmp));
00171                 scalar dtbu = romega*acos(cosdtbu);
00172                 scalar decay = exp(-Kbr*dtbu);
00173 
00174                 scalar rNew = decay*r;
00175                 if (rNew < r)
00176                 {
00177                     p.d() = 2*rNew;
00178                     p.dev() = 0;
00179                     p.ddev() = 0;
00180                 }
00181             }
00182         }
00183     }
00184     else
00185     {
00186         // reset droplet distortion parameters
00187         p.dev() = 0;
00188         p.ddev() = 0;
00189     }
00190 }
00191 
00192 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00193 
00194 } // End namespace Foam
00195 
00196 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines