FreeFOAM The Cross-Platform CFD Toolkit
parcel.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 "parcel.H"
27 
28 #include <dieselSpray/spray.H>
29 #include <dieselSpray/dragModel.H>
32 #include <dieselSpray/wallModel.H>
33 #include <OpenFOAM/wallPolyPatch.H>
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
44 }
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
49 (
50  const Cloud<parcel>& cloud,
51  const vector& position,
52  const label cellI,
53  const vector& n,
54  const scalar d,
55  const scalar T,
56  const scalar m,
57  const scalar y,
58  const scalar yDot,
59  const scalar ct,
60  const scalar ms,
61  const scalar tTurb,
62  const scalar liquidCore,
63  const scalar injector,
64  const vector& U,
65  const vector& Uturb,
66  const scalarField& X,
67  const List<word>& liquidNames
68 )
69 :
70  Particle<parcel>(cloud, position, cellI),
71  liquidComponents_
72  (
73  liquidNames
74  ),
75  d_(d),
76  T_(T),
77  m_(m),
78  y_(y),
79  yDot_(yDot),
80  ct_(ct),
81  ms_(ms),
82  tTurb_(tTurb),
83  liquidCore_(liquidCore),
84  injector_(injector),
85  U_(U),
86  Uturb_(Uturb),
87  n_(n),
88  X_(X),
89  tMom_(GREAT)
90 {}
91 
92 
93 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
94 
96 {
97  const polyMesh& mesh = cloud().pMesh();
98  const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
99 
100  const liquidMixture& fuels = sDB.fuels();
101 
102  scalar deltaT = sDB.runTime().deltaT().value();
103  label Nf = fuels.components().size();
104  label Ns = sDB.composition().Y().size();
105 
106  // Calculate the interpolated gas properties at the position of the parcel
107  vector Up = sDB.UInterpolator().interpolate(position(), cell())
108  + Uturb();
109  scalar rhog = sDB.rhoInterpolator().interpolate(position(), cell());
110  scalar pg = sDB.pInterpolator().interpolate(position(), cell());
111  scalar Tg = sDB.TInterpolator().interpolate(position(), cell());
112 
113  scalarField Yfg(Nf, 0.0);
114 
115  scalar cpMixture = 0.0;
116  for(label i=0; i<Ns; i++)
117  {
118  const volScalarField& Yi = sDB.composition().Y()[i];
119  if (sDB.isLiquidFuel()[i])
120  {
121  label j = sDB.gasToLiquidIndex()[i];
122  scalar Yicelli = Yi[cell()];
123  Yfg[j] = Yicelli;
124  }
125  cpMixture += Yi[cell()]*sDB.gasProperties()[i].Cp(Tg);
126  }
127 
128  // correct the gaseous temperature for evaporated fuel
129 
130  scalar cellV = sDB.mesh().V()[cell()];
131  scalar cellMass = rhog*cellV;
132  Tg += sDB.shs()[cell()]/(cpMixture*cellMass);
133  Tg = max(200.0, Tg);
134 
135  scalar tauMomentum = GREAT;
136  scalar tauHeatTransfer = GREAT;
137  scalarField tauEvaporation(Nf, GREAT);
138  scalarField tauBoiling(Nf, GREAT);
139 
140  bool keepParcel = true;
141 
142  setRelaxationTimes
143  (
144  cell(),
145  tauMomentum,
146  tauEvaporation,
147  tauHeatTransfer,
148  tauBoiling,
149  sDB,
150  rhog,
151  Up,
152  Tg,
153  pg,
154  Yfg,
155  m()*fuels.Y(X()),
156  deltaT
157  );
158 
159 
160  // set the end-time for the track
161  scalar tEnd = (1.0 - stepFraction())*deltaT;
162 
163  // set the maximum time step for this parcel
164  scalar dtMax = min
165  (
166  tEnd,
167  min
168  (
169  tauMomentum,
170  min
171  (
172  1.0e+10*mag(min(tauEvaporation)), // evaporation is not an issue
173  min
174  (
175  mag(tauHeatTransfer),
176  mag(min(tauBoiling))
177  )
178  )
179  )
180  )/sDB.subCycles();
181 
182  // prevent the number of subcycles from being too many
183  // (10 000 seems high enough)
184  dtMax = max(dtMax, 1.0e-4*tEnd);
185 
186  bool switchProcessor = false;
187  vector planeNormal = vector::zero;
188  if (sDB.twoD())
189  {
190  planeNormal = n() ^ sDB.axisOfSymmetry();
191  planeNormal /= mag(planeNormal);
192  }
193 
194  // move the parcel until there is no 'timeLeft'
195  while (keepParcel && tEnd > SMALL && !switchProcessor)
196  {
197  // set the lagrangian time-step
198  scalar dt = min(dtMax, tEnd);
199 
200  // remember which cell the parcel is in
201  // since this will change if a face is hit
202  label celli = cell();
203  scalar p = sDB.p()[celli];
204 
205  // track parcel to face, or end of trajectory
206  if (keepParcel)
207  {
208  // Track and adjust the time step if the trajectory is not completed
209  dt *= trackToFace(position() + dt*U_, sDB);
210 
211  // Decrement the end-time acording to how much time the track took
212  tEnd -= dt;
213 
214  // Set the current time-step fraction.
215  stepFraction() = 1.0 - tEnd/deltaT;
216 
217  if (onBoundary()) // hit face
218  {
220  }
221  }
222 
223  if (keepParcel && sDB.twoD())
224  {
225  scalar z = position() & sDB.axisOfSymmetry();
226  vector r = position() - z*sDB.axisOfSymmetry();
227  if (mag(r) > SMALL)
228  {
230  }
231  }
232 
233  // **** calculate the lagrangian source terms ****
234  // First we get the 'old' properties.
235  // and then 'update' them to get the 'new'
236  // properties.
237  // The difference is then added to the source terms.
238 
239  scalar oRho = fuels.rho(p, T(), X());
240  scalarField oMass(Nf, 0.0);
241  scalar oHg = 0.0;
242  scalar oTotMass = m();
243  scalarField oYf(fuels.Y(X()));
244 
245  forAll(oMass, i)
246  {
247  oMass[i] = m()*oYf[i];
248  label j = sDB.liquidToGasIndex()[i];
249  oHg += oYf[i]*sDB.gasProperties()[j].Hs(T());
250  }
251 
252  vector oMom = m()*U();
253  scalar oHv = fuels.hl(p, T(), X());
254  scalar oH = oHg - oHv;
255  scalar oPE = (p - fuels.pv(p, T(), X()))/oRho;
256 
257  // update the parcel properties (U, T, D)
258  updateParcelProperties
259  (
260  dt,
261  sDB,
262  celli,
263  face()
264  );
265 
266  scalar nRho = fuels.rho(p, T(), X());
267  scalar nHg = 0.0;
268  scalarField nMass(Nf, 0.0);
269  scalarField nYf(fuels.Y(X()));
270 
271  forAll(nMass, i)
272  {
273  nMass[i] = m()*nYf[i];
274  label j = sDB.liquidToGasIndex()[i];
275  nHg += nYf[i]*sDB.gasProperties()[j].Hs(T());
276  }
277 
278  vector nMom = m()*U();
279  scalar nHv = fuels.hl(p, T(), X());
280  scalar nH = nHg - nHv;
281  scalar nPE = (p - fuels.pv(p, T(), X()))/nRho;
282 
283  // Update the Spray Source Terms
284  forAll(nMass, i)
285  {
286  sDB.srhos()[i][celli] += oMass[i] - nMass[i];
287  }
288  sDB.sms()[celli] += oMom - nMom;
289 
290  sDB.shs()[celli] +=
291  oTotMass*(oH + oPE)
292  - m()*(nH + nPE);
293 
294  // Remove evaporated mass from stripped mass
295  ms() -= ms()*(oTotMass-m())/oTotMass;
296 
297  // remove parcel if it is 'small'
298  if (m() < 1.0e-12)
299  {
300  keepParcel = false;
301 
302  // ... and add the removed 'stuff' to the gas
303  forAll(nMass, i)
304  {
305  sDB.srhos()[i][celli] += nMass[i];
306  }
307 
308  sDB.sms()[celli] += nMom;
309  sDB.shs()[celli] += m()*(nH + nPE);
310  }
311 
312  if (onBoundary() && keepParcel)
313  {
314  if (face() > -1)
315  {
316  if (isA<processorPolyPatch>(pbMesh[patch(face())]))
317  {
318  switchProcessor = true;
319  }
320  }
321  }
322  }
323 
324  return keepParcel;
325 }
326 
327 
328 void Foam::parcel::updateParcelProperties
329 (
330  const scalar dt,
331  spray& sDB,
332  const label celli,
333  const label facei
334 )
335 {
336  const liquidMixture& fuels = sDB.fuels();
337 
338  label Nf = sDB.fuels().components().size();
339  label Ns = sDB.composition().Y().size();
340 
341  // calculate mean molecular weight
342  scalar W = 0.0;
343  for(label i=0; i<Ns; i++)
344  {
345  W += sDB.composition().Y()[i][celli]/sDB.gasProperties()[i].W();
346 
347  }
348  W = 1.0/W;
349 
350  // Calculate the interpolated gas properties at the position of the parcel
351  vector Up = sDB.UInterpolator().interpolate(position(), celli, facei)
352  + Uturb();
353  scalar rhog = sDB.rhoInterpolator().interpolate(position(), celli, facei);
354  scalar pg = sDB.pInterpolator().interpolate(position(), celli, facei);
355  scalar Tg0 = sDB.TInterpolator().interpolate(position(), celli, facei);
356 
357  // correct the gaseous temperature for evaporated fuel
358  scalar cpMix = 0.0;
359  for(label i=0; i<Ns; i++)
360  {
361  cpMix += sDB.composition().Y()[i][celli]
362  *sDB.gasProperties()[i].Cp(T());
363  }
364  scalar cellV = sDB.mesh().V()[celli];
365  scalar rho = sDB.rho()[celli];
366  scalar cellMass = rho*cellV;
367  scalar dh = sDB.shs()[celli];
368  scalarField addedMass(Nf, 0.0);
369 
370  forAll(addedMass, i)
371  {
372  addedMass[i] += sDB.srhos()[i][celli]*cellV;
373  }
374 
375  scalar Tg = Tg0 + dh/(cpMix*cellMass);
376  Tg = max(200.0, Tg);
377 
378  scalarField Yfg(Nf, 0.0);
379  forAll(Yfg, i)
380  {
381  label j = sDB.liquidToGasIndex()[i];
382  const volScalarField& Yj = sDB.composition().Y()[j];
383  scalar Yfg0 = Yj[celli];
384  Yfg[i] = (Yfg0*cellMass + addedMass[i])/(addedMass[i] + cellMass);
385  }
386 
387  scalar tauMomentum = GREAT;
388  scalar tauHeatTransfer = GREAT;
389  scalarField tauEvaporation(Nf, GREAT);
390  scalarField tauBoiling(Nf, GREAT);
391 
392  setRelaxationTimes
393  (
394  celli,
395  tauMomentum,
396  tauEvaporation,
397  tauHeatTransfer,
398  tauBoiling,
399  sDB,
400  rhog,
401  Up,
402  Tg,
403  pg,
404  Yfg,
405  m()*fuels.Y(X()),
406  dt
407  );
408 
409  scalar timeRatio = dt/tauMomentum;
410 
411  vector Ucorr = Up;
412  vector gcorr = sDB.g();
413 
414  if (sDB.twoD())
415  {
416  // remove the tangential velocity component
417  scalar v1 = Up & sDB.axisOfSymmetry();
418  scalar v2 = Up & n();
419  Ucorr = v1*sDB.axisOfSymmetry() + v2*n();
420 
421  // Remove the tangential gravity component
422  scalar g1 = gcorr & sDB.axisOfSymmetry();
423  scalar g2 = gcorr & n();
424  gcorr = g1*sDB.axisOfSymmetry() + g2*n();
425  }
426 
427  U() = (U() + timeRatio*Ucorr + gcorr*dt)/(1.0 + timeRatio);
428 
429  if (sDB.twoD())
430  {
431  vector normal = n() ^ sDB.axisOfSymmetry();
432  normal /= mag(normal);
433  scalar dU = U() & normal;
434  U() -= dU*normal;
435  }
436 
437  scalar TDroplet = T();
438  scalar oldDensity = fuels.rho(pg, T(), X());
439  scalar oldMass = m();
440  scalarField Yf0(fuels.Y(X()));
441  scalarField mi0(m()*Yf0);
442  scalarField mi(mi0);
443 
444  scalar oldhg = 0.0;
445  for (label i=0; i<Nf; i++)
446  {
447  label j = sDB.liquidToGasIndex()[i];
448  oldhg += Yf0[i]*sDB.gasProperties()[j].Hs(T());
449  }
450 
451  scalar oldhv = fuels.hl(pg, T(), X());
452  scalar Np = N(oldDensity);
453 
454  scalar newDensity = oldDensity;
455  scalar newMass = oldMass;
456  scalar newhg = oldhg;
457  scalar newhv = oldhv;
458 
459  scalar Tnew = T();
460 
461  // NN.
462  // first calculate the new temperature and droplet mass,
463  // then calculate the energy source and correct the
464  // gaseous temperature, Tg, and mass fraction, Yfg,
465  // to calculate the new properties for the parcel
466  // This procedure seems to be more stable
467  label n = 0;
468  while ((n < sDB.evaporation().nEvapIter()) && (m() > VSMALL))
469  {
470  n++;
471  // new characteristic times does not need to be calculated the first time
472  if (n > 1)
473  {
474  newDensity = fuels.rho(pg, Tnew, X());
475  newMass = m();
476  newhg = 0.0;
477  scalarField Ynew(fuels.Y(X()));
478  for(label i=0; i<Nf; i++)
479  {
480  label j = sDB.liquidToGasIndex()[i];
481  newhg += Ynew[i]*sDB.gasProperties()[j].Hs(Tnew);
482  }
483 
484  newhv = fuels.hl(pg, Tnew, X());
485 
486  scalar dm = oldMass - newMass;
487  scalar dhNew = oldMass*(oldhg-oldhv) - newMass*(newhg-newhv);
488 
489  // Prediction of new gaseous temperature and fuel mass fraction
490  Tg = Tg0 + (dh+dhNew)/(cpMix*cellMass);
491  Tg = max(200.0, Tg);
492 
493  forAll(Yfg, i)
494  {
495  label j = sDB.liquidToGasIndex()[i];
496  const volScalarField& Yj = sDB.composition().Y()[j];
497  scalar Yfg0 = Yj[celli];
498  Yfg[i] = (Yfg0*cellMass + addedMass[i] + dm)
499  /(addedMass[i] + cellMass + dm);
500  }
501 
502  setRelaxationTimes
503  (
504  celli,
505  tauMomentum,
506  tauEvaporation,
507  tauHeatTransfer,
508  tauBoiling,
509  sDB,
510  rhog,
511  Up,
512  Tg,
513  pg,
514  Yfg,
515  m()*fuels.Y(X()),
516  dt
517  );
518 
519  }
520 
521  scalar Taverage = TDroplet + (Tg - TDroplet)/3.0;
522  // for a liquid Cl \approx Cp
523  scalar liquidcL = sDB.fuels().cp(pg, TDroplet, X());
524 
525  cpMix = 0.0;
526  for(label i=0; i<Ns; i++)
527  {
528  if (sDB.isLiquidFuel()[i])
529  {
530  label j = sDB.gasToLiquidIndex()[i];
531  cpMix += Yfg[j]*sDB.gasProperties()[i].Cp(Taverage);
532  }
533  else
534  {
535  scalar Y = sDB.composition().Y()[i][celli];
536  cpMix += Y*sDB.gasProperties()[i].Cp(Taverage);
537  }
538  }
539 
540  scalar evaporationSource = 0.0;
541  scalar z = 0.0;
542  scalar tauEvap = 0.0;
543 
544  for(label i=0; i<Nf; i++)
545  {
546  tauEvap += X()[i]*fuels.properties()[i].W()/tauEvaporation[i];
547  }
548  tauEvap = fuels.W(X())/tauEvap;
549 
550 
551  if (sDB.evaporation().evaporation())
552  {
553  scalar hv = fuels.hl(pg, TDroplet, X());
554  evaporationSource =
555  hv/liquidcL/tauEvap;
556 
557  z = cpMix*tauHeatTransfer/liquidcL/tauEvap;
558  }
559 
560  if (sDB.heatTransfer().heatTransfer())
561  {
562  scalar fCorrect =
563  sDB.heatTransfer().fCorrection(z)/tauHeatTransfer;
564 
565  Tnew =
566  (TDroplet + dt*(fCorrect * Tg - evaporationSource))
567  /(1.0 + dt*fCorrect);
568 
569  // Prevent droplet temperature to go above critial value
570  Tnew = min(Tnew, fuels.Tc(X()));
571 
572  // Prevent droplet temperature to go too low
573  // Mainly a numerical stability issue
574  Tnew = max(200.0, Tnew);
575  scalar Td = Tnew;
576 
577  scalar pAtSurface = fuels.pv(pg, Td, X());
578  scalar pCompare = 0.999*pg;
579  scalar boiling = pAtSurface >= pCompare;
580  if (boiling)
581  {
582  // can not go above boiling temperature
583  scalar Terr = 1.0e-3;
584  label n=0;
585  scalar dT = 1.0;
586  scalar pOld = pAtSurface;
587  while (dT > Terr)
588  {
589  n++;
590  pAtSurface = fuels.pv(pg, Td, X());
591  if ((pAtSurface < pCompare) && (pOld < pCompare))
592  {
593  Td += dT;
594  }
595  else
596  {
597  if ((pAtSurface > pCompare) && (pOld > pCompare))
598  {
599  Td -= dT;
600  }
601  else
602  {
603  dT *= 0.5;
604  if ((pAtSurface > pCompare) && (pOld < pCompare))
605  {
606  Td -= dT;
607  }
608  else
609  {
610  Td += dT;
611  }
612  }
613  }
614  pOld = pAtSurface;
615  if (debug)
616  {
617  if (n>100)
618  {
619  Info << "n = " << n << ", T = " << Td << ", pv = " << pAtSurface << endl;
620  }
621  }
622  }
623  Tnew = Td;
624  }
625  }
626 
627  // Evaporate droplet!
628  // if the droplet is NOT boiling use implicit scheme.
629  if (sDB.evaporation().evaporation())
630  {
631  for(label i=0; i<Nf; i++)
632  {
633  // immediately evaporate mass that has reached critical
634  // condition
635  if (mag(Tnew - fuels.Tc(X())) < SMALL)
636  {
637  mi[i] = 0.0;
638  }
639  else
640  {
641  scalar Td = min(Tnew, 0.999*fuels.properties()[i].Tc());
642 
643  scalar pAtSurface = fuels.properties()[i].pv(pg, Td);
644  scalar boiling = pAtSurface >= 0.999*pg;
645 
646  if (!boiling)
647  {
648  scalar fr = dt/tauEvaporation[i];
649  mi[i] = mi0[i]/(1.0 + fr);
650  }
651  else
652  {
653  scalar fr = dt/tauBoiling[i];
654  mi[i] = mi0[i]/(1.0 + fr);
655  }
656  }
657  }
658 
659  scalar mTot = sum(mi);
660  if (mTot > VSMALL)
661  {
662  scalarField Ynew(mi/mTot);
663  scalarField Xnew(sDB.fuels().X(Ynew));
664  forAll(Xnew, i)
665  {
666  X()[i] = Xnew[i];
667  }
668  m() = mTot;
669  }
670  else
671  {
672  m() = 0.0;
673  }
674  }
675  T() = Tnew;
676  scalar rhod = fuels.rho(pg, T(), X());
677  d() = cbrt(6.0*m()/(Np*rhod*M_PI));
678  }
679 
680  T() = Tnew;
681  scalar rhod = fuels.rho(pg, T(), X());
682  m() = sum(mi);
683  d() = cbrt(6.0*m()/(Np*rhod*M_PI));
684 }
685 
686 
688 {
689  U_ = transform(T, U_);
690 }
691 
692 
694 {}
695 
696 
697 // ************************ vim: set sw=4 sts=4 et: ************************ //