ReactingParcel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
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"
27 #include "specie.H"
28 #include "CompositionModel.H"
29 #include "PhaseChangeModel.H"
30 #include "mathematicalConstants.H"
31 
32 using namespace Foam::constant::mathematical;
33 
34 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
35 
36 template<class ParcelType>
37 template<class TrackCloudType>
39 (
40  TrackCloudType& cloud,
41  trackingData& td,
42  const scalar dt,
43  const scalar Re,
44  const scalar Pr,
45  const scalar Ts,
46  const scalar nus,
47  const scalar d,
48  const scalar T,
49  const scalar mass,
50  const label idPhase,
51  const scalar YPhase,
52  const scalarField& Y,
53  scalarField& dMassPC,
54  scalar& Sh,
55  scalar& N,
56  scalar& NCpW,
57  scalarField& Cs
58 )
59 {
60  typedef typename TrackCloudType::reactingCloudType reactingCloudType;
62  cloud.composition();
63  PhaseChangeModel<reactingCloudType>& phaseChange = cloud.phaseChange();
64 
65  if (!phaseChange.active() || (YPhase < small))
66  {
67  return;
68  }
69 
70  scalarField X(composition.liquids().X(Y));
71 
72  scalar Tvap = phaseChange.Tvap(X);
73 
74  if (T < Tvap)
75  {
76  return;
77  }
78 
79  const scalar TMax = phaseChange.TMax(td.pc(), X);
80  const scalar Tdash = min(T, TMax);
81  const scalar Tsdash = min(Ts, TMax);
82 
83  scalarField hmm(dMassPC);
84 
85  // Calculate mass transfer due to phase change
86  phaseChange.calculate
87  (
88  dt,
89  this->cell(),
90  Re,
91  Pr,
92  d,
93  nus,
94  Tdash,
95  Tsdash,
96  td.pc(),
97  td.Tc(),
98  X,
99  dMassPC
100  );
101 
102  // Limit phase change mass by availability of each specie
103  dMassPC = min(mass*YPhase*Y, dMassPC);
104 
105  const scalar dMassTot = sum(dMassPC);
106 
107  // Add to cumulative phase change mass
108  phaseChange.addToPhaseChangeMass(this->nParticle_*dMassTot);
109 
110  forAll(dMassPC, i)
111  {
112  const label cid = composition.localToCarrierId(idPhase, i);
113 
114  const scalar dh = phaseChange.dh(cid, i, td.pc(), Tdash);
115  Sh -= dMassPC[i]*dh/dt;
116  }
117 
118 
119  // Update molar emissions
120  if (cloud.heatTransfer().BirdCorrection())
121  {
122  // Average molecular weight of carrier mix - assumes perfect gas
123  const scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
124 
125  forAll(dMassPC, i)
126  {
127  const label cid = composition.localToCarrierId(idPhase, i);
128 
129  const scalar Cp = composition.carrier().Cp(cid, td.pc(), Tsdash);
130  const scalar W = composition.carrier().Wi(cid);
131  const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
132 
133  const scalar Dab =
134  composition.liquids().properties()[i].D(td.pc(), Tsdash, Wc);
135 
136  // Molar flux of species coming from the particle (kmol/m^2/s)
137  N += Ni;
138 
139  // Sum of Ni*Cpi*Wi of emission species
140  NCpW += Ni*Cp*W;
141 
142  // Concentrations of emission species
143  Cs[cid] += Ni*d/(2.0*Dab);
144  }
145  }
146 }
147 
148 
149 template<class ParcelType>
151 (
152  const scalar mass0,
153  const scalarField& dMass,
154  scalarField& Y
155 ) const
156 {
157  scalar mass1 = mass0 - sum(dMass);
158 
159  // only update the mass fractions if the new particle mass is finite
160  if (mass1 > rootVSmall)
161  {
162  forAll(Y, i)
163  {
164  Y[i] = (Y[i]*mass0 - dMass[i])/mass1;
165  }
166  }
167 
168  return mass1;
169 }
170 
171 
172 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
173 
174 template<class ParcelType>
176 (
178 )
179 :
180  ParcelType(p),
181  mass0_(p.mass0_),
182  Y_(p.Y_)
183 {}
184 
185 
186 template<class ParcelType>
188 (
190  const polyMesh& mesh
191 )
192 :
193  ParcelType(p, mesh),
194  mass0_(p.mass0_),
195  Y_(p.Y_)
196 {}
197 
198 
199 // * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * * //
200 
201 template<class ParcelType>
202 template<class TrackCloudType>
204 (
205  TrackCloudType& cloud,
206  trackingData& td
207 )
208 {
209  ParcelType::setCellValues(cloud, td);
210 
211  td.pc() = td.pInterp().interpolate
212  (
213  this->coordinates(),
214  this->currentTetIndices()
215  );
216 
217  if (td.pc() < cloud.constProps().pMin())
218  {
219  if (debug)
220  {
222  << "Limiting observed pressure in cell " << this->cell()
223  << " to " << cloud.constProps().pMin() << nl << endl;
224  }
225 
226  td.pc() = cloud.constProps().pMin();
227  }
228 }
229 
230 
231 template<class ParcelType>
232 template<class TrackCloudType>
234 (
235  TrackCloudType& cloud,
236  trackingData& td,
237  const scalar dt
238 )
239 {
240  scalar addedMass = 0.0;
241  scalar maxMassI = 0.0;
242  forAll(cloud.rhoTrans(), i)
243  {
244  scalar dm = cloud.rhoTrans(i)[this->cell()];
245  maxMassI = max(maxMassI, mag(dm));
246  addedMass += dm;
247  }
248 
249  if (maxMassI < rootVSmall)
250  {
251  return;
252  }
253 
254  const scalar massCell = this->massCell(td);
255 
256  td.rhoc() += addedMass/cloud.pMesh().cellVolumes()[this->cell()];
257 
258  const scalar massCellNew = massCell + addedMass;
259  td.Uc() = (td.Uc()*massCell + cloud.UTrans()[this->cell()])/massCellNew;
260 
261  scalar CpEff = 0.0;
262  forAll(cloud.rhoTrans(), i)
263  {
264  scalar Y = cloud.rhoTrans(i)[this->cell()]/addedMass;
265  CpEff += Y*cloud.composition().carrier().Cp(i, td.pc(), td.Tc());
266  }
267 
268  const scalar Cpc = td.CpInterp().psi()[this->cell()];
269  td.Cpc() = (massCell*Cpc + addedMass*CpEff)/massCellNew;
270 
271  td.Tc() += cloud.hsTrans()[this->cell()]/(td.Cpc()*massCellNew);
272 
273  if (td.Tc() < cloud.constProps().TMin())
274  {
275  if (debug)
276  {
278  << "Limiting observed temperature in cell " << this->cell()
279  << " to " << cloud.constProps().TMin() << nl << endl;
280  }
281 
282  td.Tc() = cloud.constProps().TMin();
283  }
284 }
285 
286 
287 template<class ParcelType>
288 template<class TrackCloudType>
290 (
291  TrackCloudType& cloud,
292  trackingData& td,
293  const scalar T,
294  const scalarField& Cs,
295  scalar& rhos,
296  scalar& mus,
297  scalar& Prs,
298  scalar& kappas
299 )
300 {
301  // No correction if total concentration of emitted species is small
302  if (!cloud.heatTransfer().BirdCorrection() || (sum(Cs) < small))
303  {
304  return;
305  }
306 
307  const SLGThermo& thermo = cloud.thermo();
308 
309  // Far field carrier molar fractions
310  scalarField Xinf(thermo.carrier().species().size());
311 
312  forAll(Xinf, i)
313  {
314  Xinf[i] = thermo.carrier().Y(i)[this->cell()]/thermo.carrier().Wi(i);
315  }
316  Xinf /= sum(Xinf);
317 
318  // Molar fraction of far field species at particle surface
319  const scalar Xsff = 1.0 - min(sum(Cs)*RR*this->T_/td.pc(), 1.0);
320 
321  // Surface carrier total molar concentration
322  const scalar CsTot = td.pc()/(RR*this->T_);
323 
324  // Surface carrier composition (molar fraction)
325  scalarField Xs(Xinf.size());
326 
327  // Surface carrier composition (mass fraction)
328  scalarField Ys(Xinf.size());
329 
330  forAll(Xs, i)
331  {
332  // Molar concentration of species at particle surface
333  const scalar Csi = Cs[i] + Xsff*Xinf[i]*CsTot;
334 
335  Xs[i] = (2.0*Csi + Xinf[i]*CsTot)/3.0;
336  Ys[i] = Xs[i]*thermo.carrier().Wi(i);
337  }
338  Xs /= sum(Xs);
339  Ys /= sum(Ys);
340 
341 
342  rhos = 0;
343  mus = 0;
344  kappas = 0;
345  scalar Cps = 0;
346  scalar sumYiSqrtW = 0;
347  scalar sumYiCbrtW = 0;
348 
349  forAll(Ys, i)
350  {
351  const scalar W = thermo.carrier().Wi(i);
352  const scalar sqrtW = sqrt(W);
353  const scalar cbrtW = cbrt(W);
354 
355  rhos += Xs[i]*W;
356  mus += Ys[i]*sqrtW*thermo.carrier().mu(i, td.pc(), T);
357  kappas += Ys[i]*cbrtW*thermo.carrier().kappa(i, td.pc(), T);
358  Cps += Xs[i]*thermo.carrier().Cp(i, td.pc(), T);
359 
360  sumYiSqrtW += Ys[i]*sqrtW;
361  sumYiCbrtW += Ys[i]*cbrtW;
362  }
363 
364  Cps = max(Cps, rootVSmall);
365 
366  rhos *= td.pc()/(RR*T);
367  rhos = max(rhos, rootVSmall);
368 
369  mus /= sumYiSqrtW;
370  mus = max(mus, rootVSmall);
371 
372  kappas /= sumYiCbrtW;
373  kappas = max(kappas, rootVSmall);
374 
375  Prs = Cps*mus/kappas;
376 }
377 
378 
379 template<class ParcelType>
380 template<class TrackCloudType>
382 (
383  TrackCloudType& cloud,
384  trackingData& td,
385  const scalar dt
386 )
387 {
388  typedef typename TrackCloudType::reactingCloudType reactingCloudType;
390  cloud.composition();
391 
392 
393  // Define local properties at beginning of time step
394  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
395 
396  const scalar np0 = this->nParticle_;
397  const scalar d0 = this->d_;
398  const vector& U0 = this->U_;
399  const scalar T0 = this->T_;
400  const scalar mass0 = this->mass();
401 
402 
403  // Calc surface values
404  scalar Ts, rhos, mus, Prs, kappas;
405  this->calcSurfaceValues(cloud, td, T0, Ts, rhos, mus, Prs, kappas);
406  scalar Res = this->Re(rhos, U0, td.Uc(), d0, mus);
407 
408 
409  // Sources
410  // ~~~~~~~
411 
412  // Explicit momentum source for particle
413  vector Su = Zero;
414 
415  // Linearised momentum source coefficient
416  scalar Spu = 0.0;
417 
418  // Momentum transfer from the particle to the carrier phase
419  vector dUTrans = Zero;
420 
421  // Explicit enthalpy source for particle
422  scalar Sh = 0.0;
423 
424  // Linearised enthalpy source coefficient
425  scalar Sph = 0.0;
426 
427  // Sensible enthalpy transfer from the particle to the carrier phase
428  scalar dhsTrans = 0.0;
429 
430 
431  // 1. Compute models that contribute to mass transfer - U, T held constant
432  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
433 
434  // Phase change
435  // ~~~~~~~~~~~~
436 
437  // Mass transfer due to phase change
438  scalarField dMassPC(Y_.size(), 0.0);
439 
440  // Molar flux of species emitted from the particle (kmol/m^2/s)
441  scalar Ne = 0.0;
442 
443  // Sum Ni*Cpi*Wi of emission species
444  scalar NCpW = 0.0;
445 
446  // Surface concentrations of emitted species
447  scalarField Cs(composition.carrier().species().size(), 0.0);
448 
449  // Calc mass and enthalpy transfer due to phase change
450  calcPhaseChange
451  (
452  cloud,
453  td,
454  dt,
455  Res,
456  Prs,
457  Ts,
458  mus/rhos,
459  d0,
460  T0,
461  mass0,
462  0,
463  1.0,
464  Y_,
465  dMassPC,
466  Sh,
467  Ne,
468  NCpW,
469  Cs
470  );
471 
472 
473  // 2. Update the parcel properties due to change in mass
474  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
475 
476  scalarField dMass(dMassPC);
477  scalar mass1 = updateMassFraction(mass0, dMass, Y_);
478 
479  this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
480 
481  // Update particle density or diameter
482  if (cloud.constProps().constantVolume())
483  {
484  this->rho_ = mass1/this->volume();
485  }
486  else
487  {
488  this->d_ = cbrt(mass1/this->rho_*6.0/pi);
489  }
490 
491  // Remove the particle when mass falls below minimum threshold
492  if (np0*mass1 < cloud.constProps().minParcelMass())
493  {
494  td.keepParticle = false;
495 
496  if (cloud.solution().coupled())
497  {
498  scalar dm = np0*mass0;
499 
500  // Absorb parcel into carrier phase
501  forAll(Y_, i)
502  {
503  scalar dmi = dm*Y_[i];
504  label gid = composition.localToCarrierId(0, i);
505  scalar hs = composition.carrier().Hs(gid, td.pc(), T0);
506 
507  cloud.rhoTrans(gid)[this->cell()] += dmi;
508  cloud.hsTrans()[this->cell()] += dmi*hs;
509  }
510  cloud.UTrans()[this->cell()] += dm*U0;
511 
512  cloud.phaseChange().addToPhaseChangeMass(np0*mass1);
513  }
514 
515  return;
516  }
517 
518  // Correct surface values due to emitted species
519  correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
520  Res = this->Re(rhos, U0, td.Uc(), this->d(), mus);
521 
522 
523  // 3. Compute heat- and momentum transfers
524  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
525 
526  // Heat transfer
527  // ~~~~~~~~~~~~~
528 
529  // Calculate new particle temperature
530  this->T_ =
531  this->calcHeatTransfer
532  (
533  cloud,
534  td,
535  dt,
536  Res,
537  Prs,
538  kappas,
539  NCpW,
540  Sh,
541  dhsTrans,
542  Sph
543  );
544 
545  this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
546 
547 
548  // Motion
549  // ~~~~~~
550 
551  // Calculate new particle velocity
552  this->U_ =
553  this->calcVelocity(cloud, td, dt, Res, mus, mass1, Su, dUTrans, Spu);
554 
555 
556  // 4. Accumulate carrier phase source terms
557  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
558 
559  if (cloud.solution().coupled())
560  {
561  // Transfer mass lost to carrier mass, momentum and enthalpy sources
562  forAll(dMass, i)
563  {
564  scalar dm = np0*dMass[i];
565  label gid = composition.localToCarrierId(0, i);
566  scalar hs = composition.carrier().Hs(gid, td.pc(), T0);
567 
568  cloud.rhoTrans(gid)[this->cell()] += dm;
569  cloud.UTrans()[this->cell()] += dm*U0;
570  cloud.hsTrans()[this->cell()] += dm*hs;
571  }
572 
573  // Update momentum transfer
574  cloud.UTrans()[this->cell()] += np0*dUTrans;
575  cloud.UCoeff()[this->cell()] += np0*Spu;
576 
577  // Update sensible enthalpy transfer
578  cloud.hsTrans()[this->cell()] += np0*dhsTrans;
579  cloud.hsCoeff()[this->cell()] += np0*Sph;
580 
581  // Update radiation fields
582  if (cloud.radiation())
583  {
584  const scalar ap = this->areaP();
585  const scalar T4 = pow4(T0);
586  cloud.radAreaP()[this->cell()] += dt*np0*ap;
587  cloud.radT4()[this->cell()] += dt*np0*T4;
588  cloud.radAreaPT4()[this->cell()] += dt*np0*ap*T4;
589  }
590  }
591 }
592 
593 
594 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
595 
596 #include "ReactingParcelIO.C"
597 
598 // ************************************************************************* //
virtual scalar TMax(const scalar p, const scalarField &X) const
Return maximum/limiting temperature.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Templated phase change model class.
Definition: ReactingCloud.H:55
zeroField Su
Definition: alphaSuSp.H:1
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual scalar kappa(const label speciei, const scalar p, const scalar T) const =0
Thermal conductivity [W/m/K].
void calcPhaseChange(TrackCloudType &cloud, trackingData &td, const scalar dt, const scalar Re, const scalar Pr, const scalar Ts, const scalar nus, const scalar d, const scalar T, const scalar mass, const label idPhase, const scalar YPhase, const scalarField &YComponents, scalarField &dMassPC, scalar &Sh, scalar &N, scalar &NCpW, scalarField &Cs)
Calculate Phase change.
basicSpecieMixture & composition
scalar pc() const
Return the continuous phase pressure.
scalar updateMassFraction(const scalar mass0, const scalarField &dMass, scalarField &Y) const
Update mass fraction.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
scalarField Y_
Mass fractions of mixture [].
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
const speciesTable & species() const
Return the table of species.
const interpolation< scalar > & pInterp() const
Return const access to the interpolator for continuous phase.
ReactingParcel(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from mesh, coordinates and topology.
virtual scalar Wi(const label speciei) const =0
Molecular weight of the given specie [kg/kmol].
label localToCarrierId(const label phaseI, const label id, const bool allowNotFound=false) const
Return carrier id of component given local id.
void addToPhaseChangeMass(const scalar dMass)
Add to phase change mass.
rhoReactionThermo & thermo
Definition: createFields.H:28
virtual scalar Cp(const label phaseI, const scalarField &Y, const scalar p, const scalar T) const
Return specific heat caoacity for the phase phaseI.
PtrList< volScalarField > & Y()
Return the mass-fraction fields.
void cellValueSourceCorrection(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct cell values using latest transfer information.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
virtual Type interpolate(const vector &position, const label celli, const label facei=-1) const =0
Interpolate field to the given point in the given cell.
virtual scalar Cp(const label speciei, const scalar p, const scalar T) const =0
Heat capacity at constant pressure [J/kg/K].
virtual scalar mu(const label speciei, const scalar p, const scalar T) const =0
Dynamic viscosity [kg/m/s].
mathematical constants.
dimensionedScalar cbrt(const dimensionedScalar &ds)
void correctSurfaceValues(TrackCloudType &cloud, trackingData &td, const scalar T, const scalarField &Cs, scalar &rhos, scalar &mus, scalar &Prs, scalar &kappas)
Correct surface values due to emitted species.
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package, and provides:
Definition: SLGThermo.H:62
static const zero Zero
Definition: zero.H:97
static const char nl
Definition: Ostream.H:265
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const basicSpecieMixture & carrier() const
Return reference to the gaseous components.
Definition: SLGThermo.C:108
virtual scalar dh(const label idc, const label idl, const scalar p, const scalar T) const
Return the enthalpy per unit mass.
PtrList< coordinateSystem > coordinates(solidRegions.size())
#define WarningInFunction
Report a warning using Foam::Warning.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
const scalarList W(::W(thermo))
dimensionedScalar pow4(const dimensionedScalar &ds)
scalar Cp(const scalar p, const scalar T) const
Definition: EtoHthermo.H:2
void setCellValues(TrackCloudType &cloud, trackingData &td)
Set cell values.
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual void calculate(const scalar dt, const label celli, const scalar Re, const scalar Pr, const scalar d, const scalar nu, const scalar T, const scalar Ts, const scalar pc, const scalar Tc, const scalarField &X, scalarField &dMassPC) const =0
Update model.
const scalar RR
Universal gas constant (default in [J/kmol/K])
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Reacting parcel class with one/two-way coupling with the continuous phase.
scalar mass0_
Initial mass [kg].
Templated reacting parcel composition model class Consists of carrier species (via thermo package)...
Definition: ReactingCloud.H:52
virtual scalar Tvap(const scalarField &X) const
Return vapourisation temperature.
virtual bool active() const
Return the model &#39;active&#39; status - default active = true.
Definition: subModelBase.C:154
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97
const basicSpecieMixture & carrier() const
Return the carrier components (wrapper function)
void calc(TrackCloudType &cloud, trackingData &td, const scalar dt)
Update parcel properties over the time interval.
const liquidMixtureProperties & liquids() const
Return the global (additional) liquids.