FreeFOAM The Cross-Platform CFD Toolkit
LISA.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 <OpenFOAM/error.H>
27 
28 #include "LISA.H"
31 
32 #include <pdf/RosinRammler.H>
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 defineTypeNameAndDebug(LISA, 0);
42 
44 (
45  atomizationModel,
46  LISA,
47  dictionary
48 );
49 
50 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
51 
52 // Construct from components
54 (
55  const dictionary& dict,
56  spray& sm
57 )
58 :
59  atomizationModel(dict, sm),
60  coeffsDict_(dict.subDict(typeName+"Coeffs")),
61  rndGen_(sm.rndGen()),
62  Cl_(readScalar(coeffsDict_.lookup("Cl"))),
63  cTau_(readScalar(coeffsDict_.lookup("cTau"))),
64  Q_(readScalar(coeffsDict_.lookup("Q"))),
65  J_(readScalar(coeffsDict_.lookup("J")))
66 {}
67 
68 
69 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
70 
72 {}
73 
74 
75 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
76 
78 (
79  parcel& p,
80  const scalar deltaT,
81  const vector& vel,
82  const liquidMixture& fuels
83 ) const
84 {
85 
86 
87  const PtrList<volScalarField>& Y = spray_.composition().Y();
88 
89  label Ns = Y.size();
90  label cellI = p.cell();
91  scalar pressure = spray_.p()[cellI];
92  scalar temperature = spray_.T()[cellI];
93  scalar Taverage = p.T() + (temperature - p.T())/3.0;
94  scalar Winv = 0.0;
95 
96  for(label i=0; i<Ns; i++)
97  {
98  Winv += Y[i][cellI]/spray_.gasProperties()[i].W();
99  }
100 
101  scalar R = specie::RR*Winv;
102 
103  // ideal gas law to evaluate density
104  scalar rhoAverage = pressure/R/Taverage;
105  //scalar nuAverage = muAverage/rhoAverage;
106  scalar sigma = fuels.sigma(pressure, p.T(), p.X());
107 
108  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
109  // The We and Re numbers are to be evaluated using the 1/3 rule.
110  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
111 
112  scalar WeberNumber = p.We(vel, rhoAverage, sigma);
113 
114  scalar tau = 0.0;
115  scalar dL = 0.0;
116  scalar k = 0.0;
117  scalar muFuel = fuels.mu(pressure, p.T(), p.X());
118  scalar rhoFuel = fuels.rho(1.0e+5, p.T(), p.X());
119  scalar nuFuel = muFuel/rhoFuel;
120 
121  vector uDir = p.U()/mag(p.U());
122 
123  scalar uGas = mag(vel & uDir);
124  vector Ug = uGas*uDir;
125 
126 /*
127  TL
128  It might be the relative velocity between Liquid and Gas, but I use the
129  absolute velocity of the parcel as suggested by the authors
130 */
131 
132 // scalar U = mag(p.Urel(vel));
133  scalar U = mag(p.U());
134 
135  p.ct() += deltaT;
136 
137  scalar Q = rhoAverage/rhoFuel;
138 
139  const injectorType& it =
140  spray_.injectors()[label(p.injector())].properties();
141 
142  if (it.nHoles() > 1)
143  {
144  Info << "Warning: This atomization model is not suitable for multihole injector." << endl
145  << "Only the first hole will be used." << endl;
146  }
147 
148  const vector itPosition = it.position(0);
149  scalar pWalk = mag(p.position() - itPosition);
150 
151 // Updating liquid sheet tickness... that is the droplet diameter
152 
153  const vector direction = it.direction(0, spray_.runTime().value());
154 
155  scalar h = (p.position() - itPosition) & direction;
156 
157  scalar d = sqrt(sqr(pWalk)-sqr(h));
158 
159  scalar time = pWalk/mag(p.U());
160 
161  scalar elapsedTime = spray_.runTime().value();
162 
163  scalar massFlow = it.massFlowRate(max(0.0,elapsedTime-time));
164 
165  scalar hSheet = massFlow/(mathematicalConstant::pi*d*rhoFuel*mag(p.U()));
166 
167  p.d() = min(hSheet,p.d());
168 
169  if(WeberNumber > 27.0/16.0)
170  {
171 
172  scalar kPos = 0.0;
173  scalar kNeg = Q*pow(U, 2.0)*rhoFuel/sigma;
174 
175  scalar derivativePos = sqrt
176  (
177  Q*pow(U,2.0)
178  );
179 
180  scalar derivativeNeg =
181  (
182  8.0*pow(nuFuel, 2.0)*pow(kNeg, 3.0)
183  + Q*pow(U, 2.0)*kNeg
184  - 3.0*sigma/2.0/rhoFuel*pow(kNeg, 2.0)
185  )
186  /
187  sqrt
188  (
189  4.0*pow(nuFuel, 2.0)*pow(kNeg, 4.0)
190  + Q*pow(U, 2.0)*pow(kNeg, 2.0)
191  - sigma*pow(kNeg, 3.0)/rhoFuel
192  )
193  -
194  4.0*nuFuel*kNeg;
195 
196  scalar kOld = 0.0;
197 
198 
199  for(label i=0; i<40; i++)
200  {
201 
202  k = kPos - (derivativePos/((derivativeNeg-derivativePos)/(kNeg-kPos)));
203 
204  scalar derivativek =
205  (
206  8.0*pow(nuFuel, 2.0)*pow(k, 3.0)
207  + Q*pow(U, 2.0)*k
208  - 3.0*sigma/2.0/rhoFuel*pow(k, 2.0)
209  )
210  /
211  sqrt
212  (
213  4.0*pow(nuFuel, 2.0)*pow(k, 4.0)
214  + Q*pow(U, 2.0)*pow(k, 2.0)
215  - sigma*pow(k, 3.0)/rhoFuel
216  )
217  -
218  4.0*nuFuel*k;
219 
220  if(derivativek > 0)
221  {
222  derivativePos = derivativek;
223  kPos = k;
224  }
225  else
226  {
227  derivativeNeg = derivativek;
228  kNeg = k;
229  }
230 
231  if(mag(k-kOld)/k < 1e-4)
232  {
233  break;
234  }
235 
236  kOld = k;
237 
238  }
239 
240  scalar omegaS =
241  - 2.0 * nuFuel * pow(k, 2.0)
242  + sqrt
243  (
244  4.0*pow(nuFuel, 2.0)*pow(k, 4.0)
245  + Q*pow(U, 2.0)*pow(k, 2.0)
246  - sigma*pow(k, 3.0)/rhoFuel
247  );
248 
249  tau = cTau_/omegaS;
250 
251  dL = sqrt(8.0*p.d()/k);
252 
253  }
254  else
255  {
256 
257  k =
258  rhoAverage*pow(U, 2.0)
259  /
260  2.0*sigma;
261 
262  scalar J = pWalk*p.d()/2.0;
263 
264  tau = pow(3.0*cTau_,2.0/3.0)*cbrt(J*sigma/(sqr(Q)*pow(U,4.0)*rhoFuel));
265 
266  dL = sqrt(4.0*p.d()/k);
267  }
268 
269 
270 
271  scalar kL =
272  1.0
273  /
274  (
275  dL *
276  pow(0.5 + 1.5 * muFuel/pow((rhoFuel*sigma*dL), 0.5), 0.5)
277  );
278 
279  scalar dD = cbrt(3.0*mathematicalConstant::pi*pow(dL, 2.0)/kL);
280 
281  scalar lisaExp = 0.27;
282  scalar ambientPressure = 1.0e+5;
283 
284  scalar pRatio = spray_.ambientPressure()/ambientPressure;
285 
286  dD = dD*pow(pRatio,lisaExp);
287 
288 // modifications to take account of the flash boiling on primary breakUp
289 
290  scalar pExp = 0.135;
291 
292  scalar chi = 0.0;
293 
294  label Nf = fuels.components().size();
295 
296  scalar Td = p.T();
297 
298  for(label i = 0; i < Nf ; i++)
299  {
300 
301  if(fuels.properties()[i].pv(spray_.ambientPressure(), Td) >= 0.999*spray_.ambientPressure())
302  {
303 
304 // The fuel is boiling.....
305 // Calculation of the boiling temperature
306 
307  scalar tBoilingSurface = Td;
308 
309  label Niter = 200;
310 
311  for(label k=0; k< Niter ; k++)
312  {
313  scalar pBoil = fuels.properties()[i].pv(pressure, tBoilingSurface);
314 
315  if(pBoil > pressure)
316  {
317  tBoilingSurface = tBoilingSurface - (Td-temperature)/Niter;
318  }
319  else
320  {
321  break;
322  }
323  }
324 
325  scalar hl = fuels.properties()[i].hl(spray_.ambientPressure(), tBoilingSurface);
326  scalar iTp = fuels.properties()[i].h(spray_.ambientPressure(), Td) - spray_.ambientPressure()/fuels.properties()[i].rho(spray_.ambientPressure(), Td);
327  scalar iTb = fuels.properties()[i].h(spray_.ambientPressure(), tBoilingSurface) - spray_.ambientPressure()/fuels.properties()[i].rho(spray_.ambientPressure(), tBoilingSurface);
328 
329  chi += p.X()[i]*(iTp-iTb)/hl;
330 
331  }
332  }
333 
334  // bounding chi
335 
336  chi = max(chi, 0.0);
337  chi = min(chi, 1.0);
338 
339  // modifing dD to take account of flash boiling
340 
341  dD = dD*(1.0 - chi*pow(pRatio, -pExp));
342 
343  scalar lBU = Cl_ * mag(p.U())*tau;
344 
345  if(pWalk > lBU)
346  {
347 
348  p.liquidCore() = 0.0;
349 
350 // calculate the new diameter with a Rosin Rammler distribution
351 
352  scalar minValue = min(p.d(),dD/10.0);
353 
354  scalar maxValue = dD;
355 
356  if(maxValue - minValue < SMALL)
357  {
358  minValue = p.d()/10.0;
359  }
360 
361  scalar range = maxValue - minValue;
362 
363  scalar y = 0;
364  scalar x = 0;
365 
366  bool success = false;
367 
368  while(!success)
369  {
370 
371  x = minValue + range*rndGen_.scalar01();
372  y = rndGen_.scalar01();
373 
374  scalar p = 0.0;
375 
376  scalar nExp = 1;
377 
378  scalar xx = pow(x/dD, nExp);
379 
380  p = xx*exp(-xx);
381 
382  if (y<p)
383  {
384  success = true;
385  }
386 
387  }
388 
389 // New droplet diameter
390 
391  p.d() = x;
392  p.ct() = 0.0;
393 
394  }
395 
396 
397 }
398 
399 
400 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
401 
402 } // End namespace Foam
403 
404 // ************************ vim: set sw=4 sts=4 et: ************************ //