ReactingHeterogeneousParcel.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2018-2020 OpenCFD Ltd.
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
29 #include "mathematicalConstants.H"
30 
31 using namespace Foam::constant::mathematical;
32 
33 
34 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
35 
36 template<class ParcelType>
37 template<class TrackCloudType>
39 (
40  TrackCloudType& cloud,
41  trackingData& td,
42  const scalar p,
43  const scalar T,
44  const label idS
45 ) const
46 {
47  return cloud.composition().Cp(idS, this->Y_, p, T);
48 }
49 
50 
51 template<class ParcelType>
52 template<class TrackCloudType>
54 (
55  TrackCloudType& cloud,
56  trackingData& td,
57  const scalar p,
58  const scalar T,
59  const label idS
60 ) const
61 {
62  return cloud.composition().Hs(idS, this->Y_, p, T);
63 }
64 
65 
66 template<class ParcelType>
67 template<class TrackCloudType>
69 (
70  TrackCloudType& cloud,
71  trackingData& td,
72  const scalar p,
73  const scalar T,
74  const label idS
75 ) const
76 {
77  return cloud.composition().L(idS, this->Y_, p, T);
78 }
79 
80 
81 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
82 
83 
84 template<class ParcelType>
85 template<class TrackCloudType>
87 (
88  TrackCloudType& cloud,
89  const scalarField& dMass,
90  const scalar p,
91  const scalar T
92 )
93 {
94  const auto& composition = cloud.composition();
95 
96  scalarField dVolSolid(dMass.size(), Zero);
97  forAll(dVolSolid, i)
98  {
99  dVolSolid[i] =
100  dMass[i]/composition.solids().properties()[i].rho();
101  }
102 
103  return sum(dVolSolid);
104 }
105 
106 
107 template<class ParcelType>
108 template<class TrackCloudType>
110 (
111  TrackCloudType& cloud,
112  trackingData& td,
113  const scalar dt
114 )
115 {
116  typedef typename TrackCloudType::reactingCloudType reactingCloudType;
117 
119  cloud.composition();
120 
121  // Define local properties at beginning of timestep
122  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
123 
124  const scalar np0 = this->nParticle_;
125  const scalar d0 = this->d_;
126  const vector& U0 = this->U_;
127  const scalar T0 = this->T_;
128  const scalar mass0 = this->mass();
129 
130  const scalar pc = td.pc();
131 
132  const label idS = composition.idSolid();
133 
134  // Calc surface values
135  scalar Ts, rhos, mus, Prs, kappas;
136  this->calcSurfaceValues(cloud, td, T0, Ts, rhos, mus, Prs, kappas);
137  scalar Res = this->Re(rhos, U0, td.Uc(), d0, mus);
138 
139  // Sources
140  //~~~~~~~~
141 
142  // Explicit momentum source for particle
143  vector Su = Zero;
144 
145  // Linearised momentum source coefficient
146  scalar Spu = 0.0;
147 
148  // Momentum transfer from the particle to the carrier phase
149  vector dUTrans = Zero;
150 
151  // Explicit enthalpy source for particle
152  scalar Sh = 0.0;
153 
154  // Linearised enthalpy source coefficient
155  scalar Sph = 0.0;
156 
157  // Sensible enthalpy transfer from the particle to the carrier phase
158  scalar dhsTrans = 0.0;
159 
160  // Molar flux of species emitted from the particle (kmol/m^2/s)
161  scalar Ne = 0.0;
162 
163  // Surface concentrations of emitted species
164  scalarField Cs(composition.carrier().species().size(), 0.0);
165 
166  // Sum Ni*Cpi*Wi of emission species
167  scalar NCpW = 0.0;
168 
169 
170  // Heterogeneous reactions
171  // ~~~~~~~~~~~~~~~~~
172 
173  // Change in carrier phase composition due to surface reactions
174  scalarField dMassSRSolid(this->Y_.size(), 0.0);
175  scalarField dMassSRCarrier(composition.carrier().species().size(), 0.0);
176 
177  // Calc mass and enthalpy transfer due to reactions
178  calcHeterogeneousReactions
179  (
180  cloud,
181  td,
182  dt,
183  Res,
184  mus/rhos,
185  d0,
186  T0,
187  mass0,
188  canCombust_,
189  Ne,
190  NCpW,
191  this->Y_,
192  F_,
193  dMassSRSolid,
194  dMassSRCarrier,
195  Sh,
196  dhsTrans
197  );
198 
199  // 2. Update the parcel properties due to change in mass
200  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201 
202  scalarField dMassSolid(dMassSRSolid);
203 
204  scalar mass1 = mass0 - sum(dMassSolid);
205 
206  // Remove the particle when mass falls below minimum threshold
207  if (np0*mass1 < cloud.constProps().minParcelMass())
208  {
209  td.keepParticle = false;
210 
211  return;
212  }
213 
214  // Only solid is used to update mass fractions
215  (void)this->updateMassFraction(mass0, dMassSolid, this->Y_);
216 
217 
218  if
219  (
220  cloud.constProps().volUpdateType()
221  == constantProperties::volumeUpdateType::mUndefined
222  )
223  {
224  if (cloud.constProps().constantVolume())
225  {
226  this->rho_ = mass1/this->volume();
227  }
228  else
229  {
230  this->d_ = cbrt(mass1/this->rho_*6/pi);
231  }
232  }
233  else
234  {
235  switch (cloud.constProps().volUpdateType())
236  {
237  case constantProperties::volumeUpdateType::mConstRho :
238  {
239  this->d_ = cbrt(mass1/this->rho_*6/pi);
240  break;
241  }
242  case constantProperties::volumeUpdateType::mConstVol :
243  {
244  this->rho_ = mass1/this->volume();
245  break;
246  }
247  case constantProperties::volumeUpdateType::mUpdateRhoAndVol :
248  {
249  scalar deltaVol =
250  updatedDeltaVolume
251  (
252  cloud,
253  dMassSolid,
254  pc,
255  T0
256  );
257 
258  this->rho_ = mass1/(this->volume() + deltaVol);
259  this->d_ = cbrt(mass1/this->rho_*6.0/pi);
260  break;
261  }
262  }
263  }
264  // Correct surface values due to emitted species
265  this->correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
266  Res = this->Re(rhos, U0, td.Uc(), this->d_, mus);
267 
268 
269  // 3. Compute heat- and momentum transfers
270  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
271 
272  // Heat transfer
273  // ~~~~~~~~~~~~~
274 
275  // Calculate new particle temperature
276  this->T_ =
277  this->calcHeatTransfer
278  (
279  cloud,
280  td,
281  dt,
282  Res,
283  Prs,
284  kappas,
285  NCpW,
286  Sh,
287  dhsTrans,
288  Sph
289  );
290 
291  //DebugVar(np0);
292 
293  this->Cp_ = CpEff(cloud, td, pc, this->T_, idS);
294 
295  // Motion
296  // ~~~~~~
297 
298  // Calculate new particle velocity
299  this->U_ =
300  this->calcVelocity(cloud, td, dt, Res, mus, mass1, Su, dUTrans, Spu);
301 
302 
303  // 4. Accumulate carrier phase source terms
304  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
305 
306  if (cloud.solution().coupled())
307  {
308  // No mapping between solid components and carrier phase
309  /*
310  forAll(this->Y_, i)
311  {
312  scalar dm = np0*dMassSolid[i];
313  label gid = composition.localToCarrierId(SLD, i);
314  scalar hs = composition.carrier().Hs(gid, pc, T0);
315  cloud.rhoTrans(gid)[this->cell()] += dm;
316  cloud.UTrans()[this->cell()] += dm*U0;
317  cloud.hsTrans()[this->cell()] += dm*hs;
318  }
319  */
320 
321  forAll(dMassSRCarrier, i)
322  {
323  scalar dm = np0*dMassSRCarrier[i];
324  scalar hs = composition.carrier().Hs(i, pc, T0);
325  cloud.rhoTrans(i)[this->cell()] += dm;
326  cloud.UTrans()[this->cell()] += dm*U0;
327  cloud.hsTrans()[this->cell()] += dm*hs;
328  }
329 
330  // Update momentum transfer
331  cloud.UTrans()[this->cell()] += np0*dUTrans;
332  cloud.UCoeff()[this->cell()] += np0*Spu;
333 
334  // Update sensible enthalpy transfer
335  cloud.hsTrans()[this->cell()] += np0*dhsTrans;
336  cloud.hsCoeff()[this->cell()] += np0*Sph;
337 
338  // Update radiation fields
339  if (cloud.radiation())
340  {
341  const scalar ap = this->areaP();
342  const scalar T4 = pow4(T0);
343  cloud.radAreaP()[this->cell()] += dt*np0*ap;
344  cloud.radT4()[this->cell()] += dt*np0*T4;
345  cloud.radAreaPT4()[this->cell()] += dt*np0*ap*T4;
346  }
347  }
348 }
349 
350 
351 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
352 
353 template<class ParcelType>
354 template<class TrackCloudType>
356 (
357  TrackCloudType& cloud,
358  trackingData& td,
359  const scalar dt,
360  const scalar Re,
361  const scalar nu,
362  const scalar d,
363  const scalar T,
364  const scalar mass,
365  const label canCombust,
366  const scalar Ne,
367  scalar& NCpW,
368  const scalarField& YSolid,
369  scalarField& F,
370  scalarField& dMassSRSolid,
371  scalarField& dMassSRCarrier,
372  scalar& Sh,
373  scalar& dhsTrans
374 ) const
375 {
376  // Check that model is active
377  if (!cloud.heterogeneousReaction().active())
378  {
379  return;
380  }
381 
382  // Initialise demand-driven constants
383  (void)cloud.constProps().hRetentionCoeff();
384  (void)cloud.constProps().TMax();
385 
386  // Check that model is active
387  if (canCombust != 1)
388  {
389  return;
390  }
391 
392  // Update reactions
393  const scalar hReaction = cloud.heterogeneousReaction().calculate
394  (
395  dt,
396  Re,
397  nu,
398  this->cell(),
399  d,
400  T,
401  td.Tc(),
402  td.pc(),
403  td.rhoc(),
404  mass,
405  YSolid,
406  F,
407  Ne,
408  NCpW,
409  dMassSRSolid,
410  dMassSRCarrier
411  );
412 
413  cloud.heterogeneousReaction().addToSurfaceReactionMass
414  (
415  this->nParticle_*sum(dMassSRSolid)
416  );
417 
418  const scalar xsi = min(T/cloud.constProps().TMax(), 1.0);
419  const scalar coeff =
420  (1.0 - xsi*xsi)*cloud.constProps().hRetentionCoeff();
421 
422  Sh += coeff*hReaction/dt;
423 
424  dhsTrans += (1.0 - coeff)*hReaction;
425 }
426 
427 
428 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
429 
430 template<class ParcelType>
432 (
434 )
435 :
436  ParcelType(p),
437  F_(p.F_),
438  canCombust_(p.canCombust_)
439 {}
440 
441 
442 template<class ParcelType>
444 (
445  const ReactingHeterogeneousParcel<ParcelType>& p,
446  const polyMesh& mesh
447 )
448 :
449  ParcelType(p, mesh),
450  F_(p.F_),
451  canCombust_(p.canCombust_)
452 {}
453 
454 
455 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
456 
458 
459 // ************************************************************************* //
mathematicalConstants.H
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
F
volVectorField F(fluid.F())
Foam::fieldTypes::volume
const wordList volume
Standard volume field types (scalar, vector, tensor, etc)
Foam::ReactingHeterogeneousParcel::trackingData
ParcelType::trackingData trackingData
Use base tracking data.
Definition: ReactingHeterogeneousParcel.H:109
composition
basicSpecieMixture & composition
Definition: createFieldRefs.H:6
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
nu
volScalarField & nu
Definition: readMechanicalProperties.H:176
Su
zeroField Su
Definition: alphaSuSp.H:1
Foam::ReactingHeterogeneousParcel
Reacting heterogeneous Parcel.
Definition: ReactingHeterogeneousParcel.H:52
Foam::CompositionModel
Templated reacting parcel composition model class Consists of carrier species (via thermo package),...
Definition: ReactingCloud.H:58
Foam::pow4
dimensionedScalar pow4(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:100
Foam::Field< scalar >
ReactingHeterogeneousParcelIO.C
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:58
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::ReactingHeterogeneousParcel::ReactingHeterogeneousParcel
ReactingHeterogeneousParcel(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from mesh, position and topology.
Definition: ReactingHeterogeneousParcelI.H:65
Foam::cloud
A cloud is a registry collection of lagrangian particles.
Definition: cloud.H:57
Foam::constant::mathematical
Mathematical constants.
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::Vector< scalar >
Foam::ReactingHeterogeneousParcel::calcHeterogeneousReactions
void calcHeterogeneousReactions(TrackCloudType &cloud, trackingData &td, const scalar dt, const scalar Res, const scalar nu, const scalar d, const scalar T, const scalar mass, const label canCombust, const scalar N, scalar &NCpW, const scalarField &YSolid, scalarField &F, scalarField &dMassSRSolid, scalarField &dMassSRCarrier, scalar &Sh, scalar &dhsTrans) const
Calculate surface reactions.
Definition: ReactingHeterogeneousParcel.C:356
Foam::Re
scalarField Re(const UList< complex > &cf)
Extract real component.
Definition: complexField.C:159
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
Cs
const scalarField & Cs
Definition: solveBulkSurfactant.H:10
Foam::ReactingHeterogeneousParcel::updatedDeltaVolume
scalar updatedDeltaVolume(TrackCloudType &cloud, const scalarField &dMass, const scalar p, const scalar T)
Return change of volume due to mass exchange.
Foam::cbrt
dimensionedScalar cbrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:155
T0
scalar T0
Definition: createFields.H:22
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:54
Foam::ReactingHeterogeneousParcel::calc
void calc(TrackCloudType &cloud, trackingData &td, const scalar dt)
Update parcel properties over the time interval.
Definition: ReactingHeterogeneousParcel.C:110
ReactingHeterogeneousParcel.H