FreeFOAM The Cross-Platform CFD Toolkit
RutlandFlashBoil.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 "RutlandFlashBoil.H"
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 defineTypeNameAndDebug(RutlandFlashBoil, 0);
39 
41 (
42  evaporationModel,
43  RutlandFlashBoil,
44  dictionary
45 );
46 
47 
48 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
49 
50 // Construct from dictionary
52 (
53  const dictionary& dict
54 )
55 :
56  evaporationModel(dict),
57  evapDict_(dict.subDict(typeName + "Coeffs")),
58  preReScFactor_(readScalar(evapDict_.lookup("preReScFactor"))),
59  ReExponent_(readScalar(evapDict_.lookup("ReExponent"))),
60  ScExponent_(readScalar(evapDict_.lookup("ScExponent"))),
61  evaporationScheme_(evapDict_.lookup("evaporationScheme")),
62  nEvapIter_(0)
63 {
64  if (evaporationScheme_ == "implicit")
65  {
66  nEvapIter_ = 2;
67  }
68  else if (evaporationScheme_ == "explicit")
69  {
70  nEvapIter_ = 1;
71  }
72  else
73  {
75  << "evaporationScheme type " << evaporationScheme_
76  << " unknown.\n"
77  << "Use implicit or explicit."
78  << abort(FatalError);
79  }
80 }
81 
82 
83 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
84 
86 {}
87 
88 
89 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
90 
92 {
93  return true;
94 }
95 
96 // Correlation for the Sherwood Number
98 (
99  const scalar ReynoldsNumber,
100  const scalar SchmidtNumber
101 ) const
102 {
103  return 2.0 + preReScFactor_*pow(ReynoldsNumber,ReExponent_)*pow(SchmidtNumber,ScExponent_);
104 }
105 
107 (
108  const scalar diameter,
109  const scalar liquidDensity,
110  const scalar rhoFuelVapor,
111  const scalar massDiffusionCoefficient,
112  const scalar ReynoldsNumber,
113  const scalar SchmidtNumber,
114  const scalar Xs,
115  const scalar Xf,
116  const scalar m0,
117  const scalar dm,
118  const scalar dt
119 ) const
120 {
121  scalar time = GREAT;
122  scalar lgExpr = 0.0;
123 
124  /*
125  (pressure - partialFuelVaporPressure)/
126  (pressure - pressureAtSurface)
127  = 1 + Xratio
128 
129  if the pressure @ Surface > pressure
130  this lead to boiling
131  and Xratio -> infinity (as it should)
132  ... this is numerically nasty
133 
134  NB! by N. Nordin
135  X_v,s = (p_v,s/p) X_v,d
136  where X_v,d = 1 for single component fuel
137  according to eq (3.136)
138  in D. Clerides Thesis
139  */
140 
141  scalar Xratio = (Xs - Xf)/max(SMALL, 1.0 - Xs);
142 
143  if (Xratio > 0.0)
144  {
145  lgExpr = log(1.0 + Xratio);
146  }
147 
148  // From Equation (3.79) in C. Kralj's Thesis:
149  // Note that the 2.0 (instead of 6.0) below is correct, since evaporation
150  // is d(diameter)/dt and not d(mass)/dt
151 
152  scalar Sherwood = Sh(ReynoldsNumber, SchmidtNumber);
153 
154  scalar FbExp = 0.7;
155 
156  scalar logXratio = log(1.0+Xratio);
157  scalar Fb = 1.0;
158 
159  if(logXratio > SMALL)
160  {
161  Fb = pow((1.0 + Xratio),FbExp) * log(1.0+Xratio)/Xratio;
162  }
163 
164 // TL: proposed correction to sherwood number, implemented
165 
166  Sherwood = 2.0 + (Sherwood - 2.0)/Fb;
167 
168  scalar denominator =
169  6.0 * massDiffusionCoefficient
170  * Sherwood
171  * rhoFuelVapor * lgExpr;
172 
173  if (denominator > SMALL)
174  {
175  time = max(VSMALL, liquidDensity * pow(diameter, 2.0)/denominator);
176  }
177 
178  return time;
179 }
180 
181 
183 (
184  const scalar liquidDensity,
185  const scalar cpFuel,
186  const scalar heatOfVapour,
187  const scalar kappa,
188  const scalar Nusselt,
189  const scalar deltaTemp,
190  const scalar diameter,
191  const scalar liquidCore,
192  const scalar ct,
193  const scalar tDrop,
194  const scalar tBoilingSurface,
195  const scalar vapourSurfaceEnthalpy,
196  const scalar vapourFarEnthalpy,
197  const scalar cpGas,
198  const scalar temperature,
199  const scalar kLiq
200 ) const
201 {
202 
203  scalar time = GREAT;
204 
205  // the deltaTemperature is limited to not go below .5 deg K
206  // for numerical reasons.
207  // This is probably not important, since it results in an upper
208  // limit for the boiling time... which we have anyway.
209 
210  // TL kSet to the k value at the droplet temperature, not as in the Rutland Paper
211 
212  if(liquidCore > 0.5)
213  {
214  if(tDrop > tBoilingSurface)
215  {
216  // Evaporation of the liquid sheet
217 
218  scalar psi = 2.72;
219  scalar kIncreased = psi * kLiq;
220  scalar alfa = psi * kIncreased/(liquidDensity * cpFuel);
221  scalar F = alfa * ct/sqr(0.5 * diameter);
222 
223  scalar expSum = 0.0;
224  scalar expSumOld = expSum;
225 
226  label Niter = 200;
227 
228  for(label k=0; k < Niter ; k++)
229  {
230  expSum += exp(sqr(-k*mathematicalConstant::pi*sqrt(F)/2.0));
231  if(mag(expSum-expSumOld)/expSum < 1.0e-3)
232  {
233  break;
234  }
235  expSumOld = expSum;
236  }
237  }
238  }
239  else
240  {
241  scalar dTLB = min(0.5, tDrop - tBoilingSurface);
242  scalar alfaS = 0.0;
243 
244  if(dTLB >= 0.0 && dTLB < 5.0)
245  {
246  alfaS = 0.76 * pow(dTLB, 0.26);
247  }
248  if(dTLB >= 5.0 && dTLB < 25.0)
249  {
250  alfaS = 0.027 * pow(dTLB, 2.33);
251  }
252  if(dTLB >= 25.0)
253  {
254  alfaS = 13.8 * pow(dTLB, 0.39);
255  }
256 
257  scalar Gf =
258  (
259  4.0 * alfaS * dTLB * mathematicalConstant::pi * sqr(diameter/2.0)
260  )
261  /
262  heatOfVapour;
263 
264  // calculation of the heat transfer vapourization at superheated conditions (temperature>tBoilingSurface)
265  scalar G = 0.0;
266  if(temperature > tBoilingSurface)
267  {
268  scalar NusseltCorr = Nusselt ;
269  scalar A = mag((vapourFarEnthalpy-vapourSurfaceEnthalpy)/heatOfVapour);
270 
271  // TL : 2.0? or 1.0? try 1!
272  scalar B = 1.0*mathematicalConstant::pi*kappa/cpGas*diameter*NusseltCorr;
273  scalar nPos = B * log(1.0 + A)/Gf + 1.0;
274  scalar nNeg = (1.0/A)*(exp(Gf/B) - 1.0 - A) + 1.0;
275 
276  scalar Gpos = Gf*nPos;
277  scalar Gneg = Gf/nNeg;
278 
279  //scalar FgPos = Gpos + Gf - B * log( 1.0 + ( 1.0 + Gf/Gpos ) * A);
280  scalar FgNeg = Gneg + Gf - B * log( 1.0 + ( 1.0 + Gf/Gneg ) * A);
281 
282  if(FgNeg > 0.0)
283  {
284  for(label j = 0; j < 20; j++)
285  {
286  Gneg = Gneg/10.0;
287  Gneg = max(Gneg, VSMALL);
288  FgNeg = Gneg + Gf - B * log( 1.0 + ( 1.0 + Gf/Gneg ) * A);
289  if(FgNeg < 0.0)
290  {
291  break;
292  }
293  }
294  }
295 
296  FgNeg = Gneg + Gf - B * log( 1.0 + ( 1.0 + Gf/Gneg ) * A);
297 
298  G = 0.5*(Gpos+Gneg);
299  scalar Gold = -100;
300 
301  label Niter = 200;
302  label k=0;
303 
304  if(FgNeg > 0.0)
305  {
306  Info << "no convergence" << endl;
307  }
308 
309 
310  if(FgNeg < 0.0)
311  {
312 
313  for(k=0; k<Niter ; k++)
314  {
315 
316  scalar Fg = G + Gf - B * log( 1.0 + ( 1.0 + Gf/G ) * A);
317 
318  if(Fg > 0)
319  {
320  Gpos = G;
321  G = 0.5*(Gpos+Gneg);
322  }
323  else
324  {
325  Gneg = G;
326  G = 0.5*(Gpos+Gneg);
327  }
328 
329  Gold = G;
330  if(mag(G-Gold)/Gold < 1.0e-3)
331  {
332  break;
333  }
334  }
335 
336  if(k >= Niter - 1)
337  {
338  Info << " No convergence for G " << endl;
339  }
340  }
341  else
342  {
343  G = 0.0;
344  }
345  }
346 
347  time = ((4.0/3.0)*mathematicalConstant::pi*pow(diameter/2.0,3.0))*liquidDensity/(G+Gf);
348  time = max(VSMALL, time);
349  }
350 
351  return time;
352 }
353 
354 inline label RutlandFlashBoil::nEvapIter() const
355 {
356  return nEvapIter_;
357 }
358 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
359 
360 } // End namespace Foam
361 
362 // ************************ vim: set sw=4 sts=4 et: ************************ //