FreeFOAM The Cross-Platform CFD Toolkit
ReactingParcel_.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 "ReactingParcel_.H"
28 #include <specie/specie.H>
29 
30 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
31 
32 template<class ParcelType>
33 template<class TrackData>
35 (
36  TrackData& td,
37  const scalar dt,
38  const label cellI
39 )
40 {
42 
43  pc_ = td.pInterp().interpolate(this->position(), cellI);
44  if (pc_ < td.constProps().pMin())
45  {
46  WarningIn
47  (
48  "void Foam::ReactingParcel<ParcelType>::setCellValues"
49  "("
50  "TrackData&, "
51  "const scalar, "
52  "const label"
53  ")"
54  ) << "Limiting observed pressure in cell " << cellI << " to "
55  << td.constProps().pMin() << nl << endl;
56 
57  pc_ = td.constProps().pMin();
58  }
59 }
60 
61 
62 template<class ParcelType>
63 template<class TrackData>
65 (
66  TrackData& td,
67  const scalar dt,
68  const label cellI
69 )
70 {
71  scalar massCell = this->massCell(cellI);
72 
73  scalar addedMass = 0.0;
74  forAll(td.cloud().rhoTrans(), i)
75  {
76  addedMass += td.cloud().rhoTrans(i)[cellI];
77  }
78 
79  this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
80 
81  scalar massCellNew = massCell + addedMass;
82  this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
83 
84  scalar cpEff = 0;
85  if (addedMass > ROOTVSMALL)
86  {
87  forAll(td.cloud().rhoTrans(), i)
88  {
89  scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
90  cpEff +=
91  Y*td.cloud().mcCarrierThermo().speciesData()[i].Cp(this->Tc_);
92  }
93  }
94  const scalar cpc = td.cpInterp().psi()[cellI];
95  this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
96 
97  this->Tc_ += td.cloud().hsTrans()[cellI]/(this->cpc_*massCellNew);
98 }
99 
100 
101 template<class ParcelType>
102 template<class TrackData>
104 (
105  TrackData& td,
106  const label cellI,
107  const scalar T,
108  const scalarField& Cs,
109  scalar& rhos,
110  scalar& mus,
111  scalar& Pr,
112  scalar& kappa
113 )
114 {
115  // No correction if total concentration of emitted species is small
116  if (sum(Cs) < SMALL)
117  {
118  return;
119  }
120 
121  // Far field carrier molar fractions
122  scalarField Xinf(td.cloud().mcCarrierThermo().speciesData().size());
123 
124  forAll(Xinf, i)
125  {
126  Xinf[i] =
127  td.cloud().mcCarrierThermo().Y(i)[cellI]
128  /td.cloud().mcCarrierThermo().speciesData()[i].W();
129  }
130  Xinf /= sum(Xinf);
131 
132  // Molar fraction of far field species at particle surface
133  const scalar Xsff = 1.0 - min(sum(Cs)*specie::RR*this->T_/pc_, 1.0);
134 
135  // Surface carrier total molar concentration
136  const scalar CsTot = pc_/(specie::RR*this->T_);
137 
138  // Surface carrier composition (molar fraction)
139  scalarField Xs(Xinf.size());
140 
141  // Surface carrier composition (mass fraction)
142  scalarField Ys(Xinf.size());
143 
144  forAll(Xs, i)
145  {
146  // Molar concentration of species at particle surface
147  const scalar Csi = Cs[i] + Xsff*Xinf[i]*CsTot;
148 
149  Xs[i] = (2.0*Csi + Xinf[i]*CsTot)/3.0;
150  Ys[i] = Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].W();
151  }
152  Xs /= sum(Xs);
153  Ys /= sum(Ys);
154 
155 
156  rhos = 0;
157  mus = 0;
158  kappa = 0;
159  scalar cps = 0;
160  scalar sumYiSqrtW = 0;
161  scalar sumYiCbrtW = 0;
162 
163  forAll(Ys, i)
164  {
165  const scalar sqrtW =
166  sqrt(td.cloud().mcCarrierThermo().speciesData()[i].W());
167  const scalar cbrtW =
168  cbrt(td.cloud().mcCarrierThermo().speciesData()[i].W());
169 
170  rhos += Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].W();
171  mus += Ys[i]*sqrtW*td.cloud().mcCarrierThermo().speciesData()[i].mu(T);
172  kappa +=
173  Ys[i]*cbrtW*td.cloud().mcCarrierThermo().speciesData()[i].kappa(T);
174  cps += Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].Cp(T);
175 
176  sumYiSqrtW += Ys[i]*sqrtW;
177  sumYiCbrtW += Ys[i]*cbrtW;
178  }
179 
180  rhos *= pc_/(specie::RR*T);
181  mus /= sumYiSqrtW;
182  kappa /= sumYiCbrtW;
183  Pr = cps*mus/kappa;
184 }
185 
186 
187 template<class ParcelType>
189 (
190  const scalar mass0,
191  const scalarField& dMass,
192  scalarField& Y
193 ) const
194 {
195  scalar mass1 = mass0 - sum(dMass);
196 
197  // only update the mass fractions if the new particle mass is finite
198  if (mass1 > ROOTVSMALL)
199  {
200  forAll(Y, i)
201  {
202  Y[i] = (Y[i]*mass0 - dMass[i])/mass1;
203  }
204  }
205 
206  return mass1;
207 }
208 
209 
210 template<class ParcelType>
211 template<class TrackData>
213 (
214  TrackData& td,
215  const scalar dt,
216  const label cellI
217 )
218 {
219  // Define local properties at beginning of time step
220  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
221  const scalar np0 = this->nParticle_;
222  const scalar d0 = this->d_;
223  const vector& U0 = this->U_;
224  const scalar rho0 = this->rho_;
225  const scalar T0 = this->T_;
226  const scalar cp0 = this->cp_;
227  const scalar mass0 = this->mass();
228 
229 
230  // Calc surface values
231  // ~~~~~~~~~~~~~~~~~~~
232  scalar Ts, rhos, mus, Pr, kappa;
233  this->calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
234 
235  // Reynolds number
236  scalar Re = this->Re(U0, d0, rhos, mus);
237 
238 
239  // Sources
240  //~~~~~~~~
241 
242  // Explicit momentum source for particle
243  vector Su = vector::zero;
244 
245  // Momentum transfer from the particle to the carrier phase
246  vector dUTrans = vector::zero;
247 
248  // Explicit enthalpy source for particle
249  scalar Sh = 0.0;
250 
251  // Sensible enthalpy transfer from the particle to the carrier phase
252  scalar dhsTrans = 0.0;
253 
254 
255  // Phase change
256  // ~~~~~~~~~~~~
257 
258  // Mass transfer due to phase change
259  scalarField dMassPC(Y_.size(), 0.0);
260 
261  // Molar flux of species emitted from the particle (kmol/m^2/s)
262  scalar Ne = 0.0;
263 
264  // Sum Ni*Cpi*Wi of emission species
265  scalar NCpW = 0.0;
266 
267  // Surface concentrations of emitted species
268  scalarField Cs(td.cloud().mcCarrierThermo().species().size(), 0.0);
269 
270  // Calc mass and enthalpy transfer due to phase change
271  calcPhaseChange
272  (
273  td,
274  dt,
275  cellI,
276  Re,
277  Ts,
278  mus/rhos,
279  d0,
280  T0,
281  mass0,
282  0,
283  1.0,
284  Y_,
285  dMassPC,
286  Sh,
287  Ne,
288  NCpW,
289  Cs
290  );
291 
292  // Correct surface values due to emitted species
293  correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Pr, kappa);
294 
295  // Update particle component mass and mass fractions
296  scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
297 
298 
299  // Heat transfer
300  // ~~~~~~~~~~~~~
301 
302  // Calculate new particle temperature
303  scalar T1 =
304  calcHeatTransfer
305  (
306  td,
307  dt,
308  cellI,
309  Re,
310  Pr,
311  kappa,
312  d0,
313  rho0,
314  T0,
315  cp0,
316  NCpW,
317  Sh,
318  dhsTrans
319  );
320 
321 
322  // Motion
323  // ~~~~~~
324 
325  // Calculate new particle velocity
326  vector U1 =
327  calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
328 
329  dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
330 
331  // Accumulate carrier phase source terms
332  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
333  if (td.cloud().coupled())
334  {
335  // Transfer mass lost from particle to carrier mass source
336  forAll(dMassPC, i)
337  {
338  label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
339  td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
340  }
341 
342  // Update momentum transfer
343  td.cloud().UTrans()[cellI] += np0*dUTrans;
344 
345  // Update sensible enthalpy transfer
346  td.cloud().hsTrans()[cellI] += np0*dhsTrans;
347  }
348 
349 
350  // Remove the particle when mass falls below minimum threshold
351  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
352  if (mass1 < td.constProps().minParticleMass())
353  {
354  td.keepParticle = false;
355 
356  if (td.cloud().coupled())
357  {
358  // Absorb parcel into carrier phase
359  forAll(Y_, i)
360  {
361  label gid =
362  td.cloud().composition().localToGlobalCarrierId(0, i);
363  td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
364  }
365  td.cloud().UTrans()[cellI] += np0*mass1*U1;
366  td.cloud().hsTrans()[cellI] +=
367  np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
368  }
369  }
370 
371 
372  // Set new particle properties
373  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
374 
375  else
376  {
377  this->cp_ = td.cloud().composition().cp(0, Y_, pc_, T1);
378  this->T_ = T1;
379  this->U_ = U1;
380 
381  // Update particle density or diameter
382  if (td.constProps().constantVolume())
383  {
384  this->rho_ = mass1/this->volume();
385  }
386  else
387  {
388  this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
389  }
390  }
391 }
392 
393 
394 template<class ParcelType>
395 template<class TrackData>
397 (
398  TrackData& td,
399  const scalar dt,
400  const label cellI,
401  const scalar Re,
402  const scalar Ts,
403  const scalar nus,
404  const scalar d,
405  const scalar T,
406  const scalar mass,
407  const label idPhase,
408  const scalar YPhase,
409  const scalarField& YComponents,
410  scalarField& dMassPC,
411  scalar& Sh,
412  scalar& N,
413  scalar& NCpW,
414  scalarField& Cs
415 )
416 {
417  typedef PhaseChangeModel
418  <
420  > phaseChangeModelType;
421 
422  if
423  (
424  !td.cloud().phaseChange().active()
425  || T < td.constProps().Tvap()
426  || YPhase < SMALL
427  )
428  {
429  return;
430  }
431 
432  // Calculate mass transfer due to phase change
433  td.cloud().phaseChange().calculate
434  (
435  dt,
436  cellI,
437  Re,
438  d,
439  nus,
440  T,
441  Ts,
442  pc_,
443  dMassPC
444  );
445 
446  // Limit phase change mass by availability of each specie
447  dMassPC = min(mass*YPhase*YComponents, dMassPC);
448 
449  scalar dMassTot = sum(dMassPC);
450 
451  // Add to cumulative phase change mass
452  td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
453 
454  // Average molecular weight of carrier mix - assumes perfect gas
455  scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
456 
457  forAll(YComponents, i)
458  {
459  const label idc =
460  td.cloud().composition().localToGlobalCarrierId(idPhase, i);
461  const label idl = td.cloud().composition().globalIds(idPhase)[i];
462 
463  // Calculate enthalpy transfer
464  if
465  (
466  td.cloud().phaseChange().enthalpyTransfer()
467  == phaseChangeModelType::etLatentHeat
468  )
469  {
470  scalar hlp =
471  td.cloud().composition().liquids().properties()[idl].hl(pc_, T);
472 
473  Sh -= dMassPC[i]*hlp/dt;
474  }
475  else
476  {
477  // Note: enthalpies of both phases must use the same reference
478  scalar hc = td.cloud().mcCarrierThermo().speciesData()[idc].H(T);
479  scalar hp =
480  td.cloud().composition().liquids().properties()[idl].h(pc_, T);
481 
482  Sh -= dMassPC[i]*(hc - hp)/dt;
483  }
484 
485  // Update particle surface thermo properties
486  const scalar Dab =
487  td.cloud().composition().liquids().properties()[idl].D(pc_, Ts, Wc);
488 
489  const scalar Cp =
490  td.cloud().mcCarrierThermo().speciesData()[idc].Cp(Ts);
491  const scalar W = td.cloud().mcCarrierThermo().speciesData()[idc].W();
492  const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
493 
494  // Molar flux of species coming from the particle (kmol/m^2/s)
495  N += Ni;
496 
497  // Sum of Ni*Cpi*Wi of emission species
498  NCpW += Ni*Cp*W;
499 
500  // Concentrations of emission species
501  Cs[idc] += Ni*d/(2.0*Dab);
502  }
503 }
504 
505 
506 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
507 
508 template <class ParcelType>
510 (
512 )
513 :
514  ThermoParcel<ParcelType>(p),
515  mass0_(p.mass0_),
516  Y_(p.Y_),
517  pc_(p.pc_)
518 {}
519 
520 
521 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
522 
523 #include "ReactingParcelIO.C"
524 
525 // ************************ vim: set sw=4 sts=4 et: ************************ //
526