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-------------------------------------------------------------------------------
10License
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
30
31using namespace Foam::constant::mathematical;
32
33
34// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
35
36template<class ParcelType>
37template<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
51template<class ParcelType>
52template<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
66template<class ParcelType>
67template<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
84template<class ParcelType>
85template<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
107template<class ParcelType>
108template<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
353template<class ParcelType>
354template<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
430template<class ParcelType>
432(
434)
435:
436 ParcelType(p),
437 F_(p.F_),
438 canCombust_(p.canCombust_)
439{}
440
441
442template<class ParcelType>
444(
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// ************************************************************************* //
Templated reacting parcel composition model class Consists of carrier species (via thermo package),...
scalar updatedDeltaVolume(TrackCloudType &cloud, const scalarField &dMass, const scalar p, const scalar T)
Return change of volume due to mass exchange.
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.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
A cloud is a registry collection of lagrangian particles.
Definition: cloud.H:60
Class used to pass data into container.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
basicSpecieMixture & composition
volScalarField & p
const volScalarField & T
dynamicFvMesh & mesh
zeroField Su
Definition: alphaSuSp.H:1
volVectorField F(fluid.F())
Mathematical constants.
constexpr scalar pi(M_PI)
scalarField Re(const UList< complex > &cf)
Extract real component.
Definition: complexField.C:159
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
dimensionedScalar pow4(const dimensionedScalar &ds)
dimensionedScalar cbrt(const dimensionedScalar &ds)
volScalarField & nu
const scalarField & Cs
scalar T0
Definition: createFields.H:22
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333