FreeFOAM The Cross-Platform CFD Toolkit
ETAB.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "ETAB.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 defineTypeNameAndDebug(ETAB, 0);
38 
40 (
41  breakupModel,
42  ETAB,
43  dictionary
44 );
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
48 // Construct from components
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 
73 {}
74 
75 
76 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
77 
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 // ************************ vim: set sw=4 sts=4 et: ************************ //