FreeFOAM The Cross-Platform CFD Toolkit
reitzKHRT.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 "reitzKHRT.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 defineTypeNameAndDebug(reitzKHRT, 0);
38 
40 (
41  breakupModel,
42  reitzKHRT,
43  dictionary
44 );
45 
46 
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 
49 // Construct from components
51 (
52  const dictionary& dict,
53  spray& sm
54 )
55 :
56  breakupModel(dict, sm),
57  coeffsDict_(dict.subDict(typeName + "Coeffs")),
58  g_(sm.g()),
59  b0_(readScalar(coeffsDict_.lookup("B0"))),
60  b1_(readScalar(coeffsDict_.lookup("B1"))),
61  cTau_(readScalar(coeffsDict_.lookup("Ctau"))),
62  cRT_(readScalar(coeffsDict_.lookup("CRT"))),
63  msLimit_(readScalar(coeffsDict_.lookup("msLimit"))),
64  weberLimit_(readScalar(coeffsDict_.lookup("WeberLimit")))
65 {}
66 
67 
68 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
69 
71 {}
72 
73 
74 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
75 
77 (
78  parcel& p,
79  const scalar deltaT,
80  const vector& vel,
81  const liquidMixture& fuels
82 ) const
83 {
84 
85  label celli = p.cell();
86  scalar T = p.T();
87  scalar r = 0.5*p.d();
88  scalar pc = spray_.p()[celli];
89 
90  scalar sigma = fuels.sigma(pc, T, p.X());
91  scalar rhoLiquid = fuels.rho(pc, T, p.X());
92  scalar muLiquid = fuels.mu(pc, T, p.X());
93  scalar rhoGas = spray_.rho()[celli];
94  scalar Np = p.N(rhoLiquid);
95  scalar semiMass = Np*pow(p.d(), 3.0);
96 
97  scalar weGas = p.We(vel, rhoGas, sigma);
98  scalar weLiquid = p.We(vel, rhoLiquid, sigma);
99  // correct the Reynolds number. Reitz is using radius instead of diameter
100  scalar reLiquid = 0.5*p.Re(rhoLiquid, vel, muLiquid);
101  scalar ohnesorge = sqrt(weLiquid)/(reLiquid + VSMALL);
102  scalar taylor = ohnesorge*sqrt(weGas);
103 
104  vector acceleration = p.Urel(vel)/p.tMom();
105  vector trajectory = p.U()/mag(p.U());
106  scalar gt = (g_ + acceleration) & trajectory;
107 
108  // frequency of the fastest growing KH-wave
109  scalar omegaKH =
110  (0.34 + 0.38*pow(weGas, 1.5))
111  /((1 + ohnesorge)*(1 + 1.4*pow(taylor, 0.6)))
112  *sqrt(sigma/(rhoLiquid*pow(r, 3)));
113 
114  // corresponding KH wave-length.
115  scalar lambdaKH =
116  9.02
117  *r
118  *(1.0 + 0.45*sqrt(ohnesorge))
119  *(1.0 + 0.4*pow(taylor, 0.7))
120  /pow(1.0 + 0.865*pow(weGas, 1.67), 0.6);
121 
122  // characteristic Kelvin-Helmholtz breakup time
123  scalar tauKH = 3.726*b1_*r/(omegaKH*lambdaKH);
124 
125  // stable KH diameter
126  scalar dc = 2.0*b0_*lambdaKH;
127 
128  // the frequency of the fastest growing RT wavelength.
129  scalar helpVariable = mag(gt*(rhoLiquid - rhoGas));
130  scalar omegaRT = sqrt
131  (
132  2.0*pow(helpVariable, 1.5)
133  /(3.0*sqrt(3.0*sigma)*(rhoGas + rhoLiquid))
134  );
135 
136  // RT wave number
137  scalar KRT = sqrt(helpVariable/(3.0*sigma + VSMALL));
138 
139  // wavelength of the fastest growing RT frequency
140  scalar lambdaRT = 2.0*mathematicalConstant::pi*cRT_/(KRT + VSMALL);
141 
142  // if lambdaRT < diameter, then RT waves are growing on the surface
143  // and we start to keep track of how long they have been growing
144  if ((p.ct() > 0) || (lambdaRT < p.d()))
145  {
146  p.ct() += deltaT;
147  }
148 
149  // characteristic RT breakup time
150  scalar tauRT = cTau_/(omegaRT + VSMALL);
151 
152  // check if we have RT breakup
153  if ((p.ct() > tauRT) && (lambdaRT < p.d()))
154  {
155  // the RT breakup creates diameter/lambdaRT new droplets
156  p.ct() = -GREAT;
157  scalar multiplier = p.d()/lambdaRT;
158  scalar nDrops = multiplier*Np;
159  p.d() = cbrt(semiMass/nDrops);
160  }
161  // otherwise check for KH breakup
162  else if (dc < p.d())
163  {
164  // no breakup below Weber = 12
165  if (weGas > weberLimit_)
166  {
167 
168  label injector = label(p.injector());
169  scalar fraction = deltaT/tauKH;
170 
171  // reduce the diameter according to the rate-equation
172  p.d() = (fraction*dc + p.d())/(1.0 + fraction);
173 
174  scalar ms = rhoLiquid*Np*pow3(dc)*mathematicalConstant::pi/6.0;
175  p.ms() += ms;
176 
177  // Total number of parcels for the whole injection event
178  label nParcels =
179  spray_.injectors()[injector].properties()->nParcelsToInject
180  (
181  spray_.injectors()[injector].properties()->tsoi(),
182  spray_.injectors()[injector].properties()->teoi()
183  );
184 
185  scalar averageParcelMass =
186  spray_.injectors()[injector].properties()->mass()/nParcels;
187 
188  if (p.ms()/averageParcelMass > msLimit_)
189  {
190  // set the initial ms value to -GREAT. This prevents
191  // new droplets from being formed from the child droplet
192  // from the KH instability
193 
194  // mass of stripped child parcel
195  scalar mc = p.ms();
196  // Prevent child parcel from taking too much mass
197  if (mc > 0.5*p.m())
198  {
199  mc = 0.5*p.m();
200  }
201 
202  spray_.addParticle
203  (
204  new parcel
205  (
206  spray_,
207  p.position(),
208  p.cell(),
209  p.n(),
210  dc,
211  p.T(),
212  mc,
213  0.0,
214  0.0,
215  0.0,
216  -GREAT,
217  p.tTurb(),
218  0.0,
219  p.injector(),
220  p.U(),
221  p.Uturb(),
222  p.X(),
223  p.fuelNames()
224  )
225  );
226 
227  p.m() -= mc;
228  p.ms() = 0.0;
229  }
230  }
231  }
232 }
233 
234 
235 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
236 
237 } // End namespace Foam
238 
239 // ************************ vim: set sw=4 sts=4 et: ************************ //