FreeFOAM The Cross-Platform CFD Toolkit
setRelaxationTimes.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 <dieselSpray/parcel.H>
27 
28 #include <dieselSpray/spray.H>
29 #include <dieselSpray/dragModel.H>
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
40 
41 void parcel::setRelaxationTimes
42 (
43  label celli,
44  scalar& tauMomentum,
45  scalarField& tauEvaporation,
46  scalar& tauHeatTransfer,
47  scalarField& tauBoiling,
48  const spray& sDB,
49  const scalar rho,
50  const vector& Up,
51  const scalar temperature,
52  const scalar pressure,
53  const scalarField& Yfg,
54  const scalarField& m0,
55  const scalar dt
56 )
57 {
58 
59  const liquidMixture& fuels = sDB.fuels();
60 
61  scalar mCell = rho*sDB.mesh().V()[cell()];
62  scalarField mfg(Yfg*mCell);
63 
64  label Ns = sDB.composition().Y().size();
65  label Nf = fuels.components().size();
66 
67  // Tf is based on the 1/3 rule
68  scalar Tf = T() + (temperature - T())/3.0;
69 
70  // calculate mixture properties
71  scalar W = 0.0;
72  scalar kMixture = 0.0;
73  scalar cpMixture = 0.0;
74  scalar muf = 0.0;
75 
76  for(label i=0; i<Ns; i++)
77  {
78  scalar Y = sDB.composition().Y()[i][celli];
79  W += Y/sDB.gasProperties()[i].W();
80  // Using mass-fractions to average...
81  kMixture += Y*sDB.gasProperties()[i].kappa(Tf);
82  cpMixture += Y*sDB.gasProperties()[i].Cp(Tf);
83  muf += Y*sDB.gasProperties()[i].mu(Tf);
84  }
85  W = 1.0/W;
86 
87  scalarField Xf(Nf, 0.0);
88  scalarField Yf(Nf, 0.0);
89  scalarField psat(Nf, 0.0);
90  scalarField msat(Nf, 0.0);
91 
92  for(label i=0; i<Nf; i++)
93  {
94  label j = sDB.liquidToGasIndex()[i];
95  scalar Y = sDB.composition().Y()[j][celli];
96  scalar Wi = sDB.gasProperties()[j].W();
97  Yf[i] = Y;
98  Xf[i] = Y*W/Wi;
99  psat[i] = fuels.properties()[i].pv(pressure, temperature);
100  msat[i] = min(1.0, psat[i]/pressure)*Wi/W;
101  }
102 
103  scalar nuf = muf/rho;
104 
105  scalar liquidDensity = fuels.rho(pressure, T(), X());
106  scalar liquidcL = fuels.cp(pressure, T(), X());
107  scalar heatOfVapour = fuels.hl(pressure, T(), X());
108 
109  // calculate the partial rho of the fuel vapour
110  // alternative is to use the mass fraction
111  // however, if rhoFuelVap is small (zero)
112  // d(mass)/dt = 0 => no evaporation... hmmm... is that good? NO!
113 
114  // Assume equilibrium at drop-surface => pressure @ surface
115  // = vapour pressure to calculate fuel-vapour density @ surface
116  scalar pressureAtSurface = fuels.pv(pressure, T(), X());
117  scalar rhoFuelVap = pressureAtSurface*fuels.W(X())/(specie::RR*Tf);
118 
119  scalarField Xs(sDB.fuels().Xs(pressure, temperature, T(), Xf, X()));
120  scalarField Ys(Nf, 0.0);
121  scalar Wliq = 0.0;
122 
123  for(label i=0; i<Nf; i++)
124  {
125  label j = sDB.liquidToGasIndex()[i];
126  scalar Wi = sDB.gasProperties()[j].W();
127  Wliq += Xs[i]*Wi;
128  }
129 
130  for(label i=0; i<Nf; i++)
131  {
132  label j = sDB.liquidToGasIndex()[i];
133  scalar Wi = sDB.gasProperties()[j].W();
134  Ys[i] = Xs[i]*Wi/Wliq;
135  }
136 
137  scalar Reynolds = Re(Up, nuf);
138  scalar Prandtl = Pr(cpMixture, muf, kMixture);
139 
140  // calculate the characteritic times
141 
142  if(liquidCore_> 0.5)
143  {
144 // no drag for parcels in the liquid core..
145  tauMomentum = GREAT;
146  }
147  else
148  {
149  tauMomentum = sDB.drag().relaxationTime
150  (
151  Urel(Up),
152  d(),
153  rho,
154  liquidDensity,
155  nuf,
156  dev()
157  );
158  }
159 
160  // store the relaxationTime since it is needed in some breakup models.
161  tMom_ = tauMomentum;
162 
163  tauHeatTransfer = sDB.heatTransfer().relaxationTime
164  (
165  liquidDensity,
166  d(),
167  liquidcL,
168  kMixture,
169  Reynolds,
170  Prandtl
171  );
172 
173  // evaporation-properties are evaluated at averaged temperature
174  // set the boiling conditions true if pressure @ surface is 99.9%
175  // of the pressure
176  // this is mainly to put a limit on the evaporation time,
177  // since tauEvaporation is very very small close to the boiling point.
178 
179  for(label i=0; i<Nf; i++)
180  {
181  scalar Td = min(T(), 0.999*fuels.properties()[i].Tc());
182  bool boiling = fuels.properties()[i].pv(pressure, Td) >= 0.999*pressure;
183  scalar Di = fuels.properties()[i].D(pressure, Td);
184  scalar Schmidt = Sc(nuf, Di);
185 
186  scalar partialPressure = Xf[i]*pressure;
187 
188 // saturated vapour
189  if(partialPressure > psat[i])
190  {
191  tauEvaporation[i] = GREAT;
192  }
193 // not saturated vapour
194  else
195  {
196  if (!boiling)
197  {
198  // for saturation evaporation, only use 99.99% for numerical robustness
199  scalar dm = max(SMALL, 0.9999*msat[i] - mfg[i]);
200 
201  tauEvaporation[i] = sDB.evaporation().relaxationTime
202  (
203  d(),
204  fuels.properties()[i].rho(pressure, Td),
205  rhoFuelVap,
206  Di,
207  Reynolds,
208  Schmidt,
209  Xs[i],
210  Xf[i],
211  m0[i],
212  dm,
213  dt
214  );
215  }
216  else
217  {
218  scalar Nusselt =
219  sDB.heatTransfer().Nu(Reynolds, Prandtl);
220 
221 // calculating the boiling temperature of the liquid at ambient pressure
222  scalar tBoilingSurface = Td;
223 
224  label Niter = 0;
225  scalar deltaT = 10.0;
226  scalar dp0 = fuels.properties()[i].pv(pressure, tBoilingSurface) - pressure;
227  while ((Niter < 200) && (mag(deltaT) > 1.0e-3))
228  {
229  Niter++;
230  scalar pBoil = fuels.properties()[i].pv(pressure, tBoilingSurface);
231  scalar dp = pBoil - pressure;
232  if ( (dp > 0.0) && (dp0 > 0.0) )
233  {
234  tBoilingSurface -= deltaT;
235  }
236  else
237  {
238  if ( (dp < 0.0) && (dp0 < 0.0) )
239  {
240  tBoilingSurface += deltaT;
241  }
242  else
243  {
244  deltaT *= 0.5;
245  if ( (dp > 0.0) && (dp0 < 0.0) )
246  {
247  tBoilingSurface -= deltaT;
248  }
249  else
250  {
251  tBoilingSurface += deltaT;
252  }
253  }
254  }
255  dp0 = dp;
256  }
257 
258  scalar vapourSurfaceEnthalpy = 0.0;
259  scalar vapourFarEnthalpy = 0.0;
260 
261  for(label k = 0; k < sDB.gasProperties().size(); k++)
262  {
263  vapourSurfaceEnthalpy += sDB.composition().Y()[k][celli]*sDB.gasProperties()[k].H(tBoilingSurface);
264  vapourFarEnthalpy += sDB.composition().Y()[k][celli]*sDB.gasProperties()[k].H(temperature);
265  }
266 
267  scalar kLiquid = fuels.properties()[i].K(pressure, 0.5*(tBoilingSurface+T()));
268 
269  tauBoiling[i] = sDB.evaporation().boilingTime
270  (
271  fuels.properties()[i].rho(pressure, Td),
272  fuels.properties()[i].cp(pressure, Td),
273  heatOfVapour,
274  kMixture,
275  Nusselt,
276  temperature - T(),
277  d(),
278  liquidCore(),
279  sDB.runTime().value() - ct(),
280  Td,
281  tBoilingSurface,
282  vapourSurfaceEnthalpy,
283  vapourFarEnthalpy,
284  cpMixture,
285  temperature,
286  kLiquid
287  );
288  }
289 
290  }
291  }
292 }
293 
294 
295 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
296 
297 } // End namespace Foam
298 
299 // ************************ vim: set sw=4 sts=4 et: ************************ //