FreeFOAM The Cross-Platform CFD Toolkit
ReactingMultiphaseParcel.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) 2008-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 
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 template<class ParcelType>
33 
34 template<class ParcelType>
36 
37 template<class ParcelType>
39 
40 
41 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
42 
43 template<class ParcelType>
44 template<class TrackData>
46 (
47  TrackData& td,
48  const scalar p,
49  const scalar T,
50  const label idG,
51  const label idL,
52  const label idS
53 ) const
54 {
55  return
56  this->Y_[GAS]*td.cloud().composition().cp(idG, YGas_, p, T)
57  + this->Y_[LIQ]*td.cloud().composition().cp(idL, YLiquid_, p, T)
58  + this->Y_[SLD]*td.cloud().composition().cp(idS, YSolid_, p, T);
59 }
60 
61 
62 template<class ParcelType>
63 template<class TrackData>
65 (
66  TrackData& td,
67  const scalar p,
68  const scalar T,
69  const label idG,
70  const label idL,
71  const label idS
72 ) const
73 {
74  return
75  this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
76  + this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
77  + this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
78 }
79 
80 
81 template<class ParcelType>
82 template<class TrackData>
84 (
85  TrackData& td,
86  const scalar p,
87  const scalar T,
88  const label idG,
89  const label idL,
90  const label idS
91 ) const
92 {
93  return
94  this->Y_[GAS]*td.cloud().composition().L(idG, YGas_, p, T)
95  + this->Y_[LIQ]*td.cloud().composition().L(idL, YLiquid_, p, T)
96  + this->Y_[SLD]*td.cloud().composition().L(idS, YSolid_, p, T);
97 }
98 
99 
100 template<class ParcelType>
102 (
103  const scalar mass0,
104  const scalarField& dMassGas,
105  const scalarField& dMassLiquid,
106  const scalarField& dMassSolid
107 )
108 {
109  scalarField& YMix = this->Y_;
110 
111  scalar massGas =
112  this->updateMassFraction(mass0*YMix[GAS], dMassGas, YGas_);
113  scalar massLiquid =
114  this->updateMassFraction(mass0*YMix[LIQ], dMassLiquid, YLiquid_);
115  scalar massSolid =
116  this->updateMassFraction(mass0*YMix[SLD], dMassSolid, YSolid_);
117 
118  scalar massNew = max(massGas + massLiquid + massSolid, ROOTVSMALL);
119 
120  YMix[GAS] = massGas/massNew;
121  YMix[LIQ] = massLiquid/massNew;
122  YMix[SLD] = 1.0 - YMix[GAS] - YMix[LIQ];
123 
124  return massNew;
125 }
126 
127 
128 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
129 
130 template<class ParcelType>
131 template<class TrackData>
133 (
134  TrackData& td,
135  const scalar dt,
136  const label cellI
137 )
138 {
140 }
141 
142 
143 template<class ParcelType>
144 template<class TrackData>
146 (
147  TrackData& td,
148  const scalar dt,
149  const label cellI
150 )
151 {
152  scalar massCell = this->massCell(cellI);
153 
154  scalar addedMass = 0.0;
155  forAll(td.cloud().rhoTrans(), i)
156  {
157  addedMass += td.cloud().rhoTrans(i)[cellI];
158  }
159 
160  this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
161 
162  scalar massCellNew = massCell + addedMass;
163  this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
164 
165  scalar cpEff = 0;
166  if (addedMass > ROOTVSMALL)
167  {
168  forAll(td.cloud().rhoTrans(), i)
169  {
170  scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
171  cpEff +=
172  Y*td.cloud().mcCarrierThermo().speciesData()[i].Cp(this->Tc_);
173  }
174  }
175  const scalar cpc = td.cpInterp().psi()[cellI];
176  this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
177 
178  this->Tc_ += td.cloud().hsTrans()[cellI]/(this->cpc_*massCellNew);
179 }
180 
181 
182 template<class ParcelType>
183 template<class TrackData>
185 (
186  TrackData& td,
187  const scalar dt,
188  const label cellI
189 )
190 {
191  // Define local properties at beginning of timestep
192  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
193  const scalar np0 = this->nParticle_;
194  const scalar d0 = this->d_;
195  const vector& U0 = this->U_;
196  const scalar rho0 = this->rho_;
197  const scalar T0 = this->T_;
198  const scalar cp0 = this->cp_;
199  const scalar mass0 = this->mass();
200 
201  const scalar pc = this->pc_;
202 
203  const scalarField& YMix = this->Y_;
204  const label idG = td.cloud().composition().idGas();
205  const label idL = td.cloud().composition().idLiquid();
206  const label idS = td.cloud().composition().idSolid();
207 
208 
209  // Calc surface values
210  // ~~~~~~~~~~~~~~~~~~~
211  scalar Ts, rhos, mus, Pr, kappa;
213  calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
214 
215  // Reynolds number
216  scalar Re = this->Re(U0, d0, rhos, mus);
217 
218 
219  // Sources
220  //~~~~~~~~
221 
222  // Explicit momentum source for particle
223  vector Su = vector::zero;
224 
225  // Momentum transfer from the particle to the carrier phase
226  vector dUTrans = vector::zero;
227 
228  // Explicit enthalpy source for particle
229  scalar Sh = 0.0;
230 
231  // Sensible enthalpy transfer from the particle to the carrier phase
232  scalar dhsTrans = 0.0;
233 
234 
235  // Phase change in liquid phase
236  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
237 
238  // Mass transfer due to phase change
239  scalarField dMassPC(YLiquid_.size(), 0.0);
240 
241  // Molar flux of species emitted from the particle (kmol/m^2/s)
242  scalar Ne = 0.0;
243 
244  // Sum Ni*Cpi*Wi of emission species
245  scalar NCpW = 0.0;
246 
247  // Surface concentrations of emitted species
248  scalarField Cs(td.cloud().mcCarrierThermo().species().size(), 0.0);
249 
250  // Calc mass and enthalpy transfer due to phase change
251  calcPhaseChange
252  (
253  td,
254  dt,
255  cellI,
256  Re,
257  Ts,
258  mus/rhos,
259  d0,
260  T0,
261  mass0,
262  idL,
263  YMix[LIQ],
264  YLiquid_,
265  dMassPC,
266  Sh,
267  Ne,
268  NCpW,
269  Cs
270  );
271 
272 
273  // Devolatilisation
274  // ~~~~~~~~~~~~~~~~
275 
276  // Mass transfer due to devolatilisation
277  scalarField dMassDV(YGas_.size(), 0.0);
278 
279  // Calc mass and enthalpy transfer due to devolatilisation
280  calcDevolatilisation
281  (
282  td,
283  dt,
284  Ts,
285  d0,
286  T0,
287  mass0,
288  this->mass0_,
289  idG,
290  YMix[GAS],
291  YGas_,
292  canCombust_,
293  dMassDV,
294  Sh,
295  Ne,
296  NCpW,
297  Cs
298  );
299 
300  // Correct surface values due to emitted species
301  correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Pr, kappa);
302 
303 
304  // Surface reactions
305  // ~~~~~~~~~~~~~~~~~
306 
307  // Change in carrier phase composition due to surface reactions
308  scalarField dMassSRGas(YGas_.size(), 0.0);
309  scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
310  scalarField dMassSRSolid(YSolid_.size(), 0.0);
312  dMassSRCarrier
313  (
314  td.cloud().mcCarrierThermo().species().size(),
315  0.0
316  );
317 
318  // Clac mass and enthalpy transfer due to surface reactions
319  calcSurfaceReactions
320  (
321  td,
322  dt,
323  cellI,
324  d0,
325  T0,
326  mass0,
327  canCombust_,
328  Ne,
329  YMix,
330  YGas_,
331  YLiquid_,
332  YSolid_,
333  dMassSRGas,
334  dMassSRLiquid,
335  dMassSRSolid,
336  dMassSRCarrier,
337  Sh,
338  dhsTrans
339  );
340 
341 
342  // Update component mass fractions
343  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
344 
345  scalarField dMassGas = dMassDV + dMassSRGas;
346  scalarField dMassLiquid = dMassPC + dMassSRLiquid;
347  scalarField dMassSolid = dMassSRSolid;
348 
349  scalar mass1 =
350  updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
351 
352 
353  // Heat transfer
354  // ~~~~~~~~~~~~~
355 
356  // Calculate new particle temperature
357  scalar T1 =
358  calcHeatTransfer
359  (
360  td,
361  dt,
362  cellI,
363  Re,
364  Pr,
365  kappa,
366  d0,
367  rho0,
368  T0,
369  cp0,
370  NCpW,
371  Sh,
372  dhsTrans
373  );
374 
375 
376  // Motion
377  // ~~~~~~
378 
379  // Calculate new particle velocity
380  vector U1 =
381  calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
382 
383  dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
384 
385 
386  // Accumulate carrier phase source terms
387  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
388 
389  if (td.cloud().coupled())
390  {
391  // Transfer mass lost from particle to carrier mass source
392  forAll(YGas_, i)
393  {
394  label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
395  td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
396  }
397  forAll(YLiquid_, i)
398  {
399  label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
400  td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
401  }
402 /*
403  // No mapping between solid components and carrier phase
404  forAll(YSolid_, i)
405  {
406  label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
407  td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
408  }
409 */
410  forAll(dMassSRCarrier, i)
411  {
412  td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
413  }
414 
415  // Update momentum transfer
416  td.cloud().UTrans()[cellI] += np0*dUTrans;
417 
418  // Update sensible enthalpy transfer
419  td.cloud().hsTrans()[cellI] += np0*dhsTrans;
420  }
421 
422 
423  // Remove the particle when mass falls below minimum threshold
424  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
425 
426  if (mass1 < td.constProps().minParticleMass())
427  {
428  td.keepParticle = false;
429 
430  if (td.cloud().coupled())
431  {
432  // Absorb parcel into carrier phase
433  forAll(YGas_, i)
434  {
435  label gid =
436  td.cloud().composition().localToGlobalCarrierId(GAS, i);
437  td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
438  }
439  forAll(YLiquid_, i)
440  {
441  label gid =
442  td.cloud().composition().localToGlobalCarrierId(LIQ, i);
443  td.cloud().rhoTrans(gid)[cellI] +=
444  np0*mass1*YMix[LIQ]*YLiquid_[i];
445  }
446 /*
447  // No mapping between solid components and carrier phase
448  forAll(YSolid_, i)
449  {
450  label gid =
451  td.cloud().composition().localToGlobalCarrierId(SLD, i);
452  td.cloud().rhoTrans(gid)[cellI] +=
453  np0*mass1*YMix[SLD]*YSolid_[i];
454  }
455 */
456  td.cloud().UTrans()[cellI] += np0*mass1*U1;
457  td.cloud().hsTrans()[cellI] +=
458  np0*mass1*HEff(td, pc, T1, idG, idL, idS); // using total h
459  }
460  }
461 
462 
463  // Set new particle properties
464  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
465 
466  else
467  {
468  this->cp_ = cpEff(td, pc, T1, idG, idL, idS);
469  this->T_ = T1;
470  this->U_ = U1;
471 
472  // Update particle density or diameter
473  if (td.constProps().constantVolume())
474  {
475  this->rho_ = mass1/this->volume();
476  }
477  else
478  {
479  this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
480  }
481  }
482 }
483 
484 
485 template<class ParcelType>
486 template<class TrackData>
488 (
489  TrackData& td,
490  const scalar dt,
491  const scalar Ts,
492  const scalar d,
493  const scalar T,
494  const scalar mass,
495  const scalar mass0,
496  const label idVolatile,
497  const scalar YVolatileTot,
498  const scalarField& YVolatile,
499  bool& canCombust,
500  scalarField& dMassDV,
501  scalar& Sh,
502  scalar& N,
503  scalar& NCpW,
504  scalarField& Cs
505 ) const
506 {
507  // Check that model is active, and that the parcel temperature is
508  // within necessary limits for devolatilisation to occur
509  if
510  (
511  !td.cloud().devolatilisation().active()
512  || T < td.constProps().Tvap()
513  )
514  {
515  return;
516  }
517 
518  // Total mass of volatiles evolved
519  const scalar dMassTot = td.cloud().devolatilisation().calculate
520  (
521  dt,
522  mass0,
523  mass,
524  T,
525  td.cloud().composition().YMixture0()[idVolatile],
526  YVolatileTot,
527  canCombust
528  );
529 
530  // Volatile mass transfer - equal components of each volatile specie
531  dMassDV = YVolatile*dMassTot;
532 
533  td.cloud().addToMassDevolatilisation(this->nParticle_*dMassTot);
534 
535  Sh -= dMassTot*td.constProps().LDevol()/dt;
536 
537  // Molar average molecular weight of carrier mix
538  const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
539 
540  // Update molar emissions
541  forAll(dMassDV, i)
542  {
543  // Note: hardcoded gaseous diffusivities for now
544  // TODO: add to carrier thermo
545  const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
546  const label id =
547  td.cloud().composition().localToGlobalCarrierId(GAS, i);
548  const scalar Cp = td.cloud().mcCarrierThermo().speciesData()[id].Cp(Ts);
549  const scalar W = td.cloud().mcCarrierThermo().speciesData()[id].W();
550  const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
551 
552  // Dab calc'd using API vapour mass diffusivity function
553  const scalar Dab =
554  3.6059e-3*(pow(1.8*Ts, 1.75))*sqrt(1.0/W + 1.0/Wc)/(this->pc_*beta);
555 
556  N += Ni;
557  NCpW += Ni*Cp*W;
558  Cs[id] += Ni*d/(2.0*Dab);
559  }
560 }
561 
562 
563 template<class ParcelType>
564 template<class TrackData>
566 (
567  TrackData& td,
568  const scalar dt,
569  const label cellI,
570  const scalar d,
571  const scalar T,
572  const scalar mass,
573  const bool canCombust,
574  const scalar N,
575  const scalarField& YMix,
576  const scalarField& YGas,
577  const scalarField& YLiquid,
578  const scalarField& YSolid,
579  scalarField& dMassSRGas,
580  scalarField& dMassSRLiquid,
581  scalarField& dMassSRSolid,
582  scalarField& dMassSRCarrier,
583  scalar& Sh,
584  scalar& dhsTrans
585 ) const
586 {
587  // Check that model is active
588  if (!td.cloud().surfaceReaction().active() || !canCombust)
589  {
590  return;
591  }
592 
593  // Update surface reactions
594  const scalar hReaction = td.cloud().surfaceReaction().calculate
595  (
596  dt,
597  cellI,
598  d,
599  T,
600  this->Tc_,
601  this->pc_,
602  this->rhoc_,
603  mass,
604  YGas,
605  YLiquid,
606  YSolid,
607  YMix,
608  N,
609  dMassSRGas,
610  dMassSRLiquid,
611  dMassSRSolid,
612  dMassSRCarrier
613  );
614 
615  td.cloud().addToMassSurfaceReaction
616  (
617  this->nParticle_
618  *(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))
619  );
620 
621  const scalar xsi = min(T/5000.0, 1.0);
622  const scalar hRetentionCoeffMod =
623  (1.0 - xsi*xsi)*td.constProps().hRetentionCoeff();
624 
625  Sh += hRetentionCoeffMod*hReaction/dt;
626 
627  dhsTrans += (1.0 - hRetentionCoeffMod)*hReaction;
628 }
629 
630 
631 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
632 
633 template <class ParcelType>
635 (
637 )
638 :
639  ReactingParcel<ParcelType>(p),
640  YGas_(p.YGas_),
641  YLiquid_(p.YLiquid_),
642  YSolid_(p.YSolid_)
643 {}
644 
645 
646 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
647 
649 
650 // ************************ vim: set sw=4 sts=4 et: ************************ //