Go to the documentation of this file.
39 template<
class ParcelType>
40 template<
class TrackCloudType>
43 TrackCloudType&
cloud,
58 return sum(dVolLiquid);
62 template<
class ParcelType>
63 template<
class TrackCloudType>
66 TrackCloudType&
cloud,
89 auto& phaseChange =
cloud.phaseChange();
92 if (!phaseChange.active())
99 scalar Tvap = phaseChange.Tvap(X);
105 const scalar TMax = phaseChange.TMax(td.
pc(), X);
106 const scalar Tdash =
min(
T, TMax);
107 const scalar Tsdash =
min(Ts, TMax);
110 phaseChange.calculate
135 dMassPC[i] =
min(mass*YPhase*
Y[i], dMassPC[i]);
140 const scalar dMassTot =
sum(dMassPC);
143 phaseChange.addToPhaseChangeMass(this->nParticle_*dMassTot);
147 const label cid =
composition.localToCarrierId(idPhase, i);
149 const scalar dh = phaseChange.dh(cid, i, td.
pc(), Tdash);
150 Sh -= dMassPC[i]*dh/dt;
155 if (
cloud.heatTransfer().BirdCorrection())
158 const scalar Wc = td.rhoc()*
RR*td.Tc()/td.
pc();
162 const label cid =
composition.localToCarrierId(idPhase, i);
166 const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
169 composition.liquids().properties()[i].D(td.
pc(), Tsdash, Wc);
178 Cs[cid] += Ni*d/(2.0*Dab);
184 template<
class ParcelType>
192 scalar mass1 = mass0 -
sum(dMass);
195 if (mass1 > ROOTVSMALL)
199 Y[i] = (
Y[i]*mass0 - dMass[i])/mass1;
209 template<
class ParcelType>
221 template<
class ParcelType>
224 const ReactingParcel<ParcelType>&
p,
236 template<
class ParcelType>
237 template<
class TrackCloudType>
240 TrackCloudType&
cloud,
244 ParcelType::setCellValues(
cloud, td);
249 this->currentTetIndices()
252 if (td.
pc() <
cloud.constProps().pMin())
257 <<
"Limiting observed pressure in cell " << this->
cell()
258 <<
" to " <<
cloud.constProps().pMin() <<
nl <<
endl;
261 td.
pc() =
cloud.constProps().pMin();
266 template<
class ParcelType>
267 template<
class TrackCloudType>
270 TrackCloudType&
cloud,
275 scalar addedMass = 0.0;
276 scalar maxMassI = 0.0;
279 scalar dm =
cloud.rhoTrans(i)[this->
cell()];
280 maxMassI =
max(maxMassI,
mag(dm));
284 if (maxMassI < ROOTVSMALL)
289 const scalar massCell = this->massCell(td);
291 td.rhoc() += addedMass/
cloud.pMesh().cellVolumes()[this->
cell()];
293 const scalar massCellNew = massCell + addedMass;
294 td.Uc() = (td.Uc()*massCell +
cloud.UTrans()[this->
cell()])/massCellNew;
299 scalar
Y =
cloud.rhoTrans(i)[this->
cell()]/addedMass;
300 CpEff +=
Y*
cloud.composition().carrier().Cp(i, td.
pc(), td.Tc());
303 const scalar Cpc = td.CpInterp().psi()[this->
cell()];
304 td.Cpc() = (massCell*Cpc + addedMass*CpEff)/massCellNew;
306 td.Tc() +=
cloud.hsTrans()[this->
cell()]/(td.Cpc()*massCellNew);
308 if (td.Tc() <
cloud.constProps().TMin())
313 <<
"Limiting observed temperature in cell " << this->
cell()
314 <<
" to " <<
cloud.constProps().TMin() <<
nl <<
endl;
317 td.Tc() =
cloud.constProps().TMin();
322 template<
class ParcelType>
323 template<
class TrackCloudType>
326 TrackCloudType&
cloud,
337 if (!
cloud.heatTransfer().BirdCorrection() || (
sum(
Cs) < SMALL))
354 const scalar Xsff = 1.0 -
min(
sum(
Cs)*
RR*this->T_/td.
pc(), 1.0);
357 const scalar CsTot = td.
pc()/(
RR*this->T_);
368 const scalar Csi =
Cs[i] + Xsff*Xinf[i]*CsTot;
370 Xs[i] = (2.0*Csi + Xinf[i]*CsTot)/3.0;
371 Ys[i] = Xs[i]*
thermo.carrier().W(i);
381 scalar sumYiSqrtW = 0;
382 scalar sumYiCbrtW = 0;
386 const scalar W =
thermo.carrier().W(i);
387 const scalar sqrtW =
sqrt(W);
388 const scalar cbrtW =
cbrt(W);
391 mus += Ys[i]*sqrtW*
thermo.carrier().mu(i, td.
pc(),
T);
392 kappas += Ys[i]*cbrtW*
thermo.carrier().kappa(i, td.
pc(),
T);
393 Cps += Xs[i]*
thermo.carrier().Cp(i, td.
pc(),
T);
395 sumYiSqrtW += Ys[i]*sqrtW;
396 sumYiCbrtW += Ys[i]*cbrtW;
399 Cps =
max(Cps, ROOTVSMALL);
401 rhos *= td.
pc()/(
RR*
T);
402 rhos =
max(rhos, ROOTVSMALL);
405 mus =
max(mus, ROOTVSMALL);
407 kappas /= sumYiCbrtW;
408 kappas =
max(kappas, ROOTVSMALL);
410 Prs = Cps*mus/kappas;
414 template<
class ParcelType>
415 template<
class TrackCloudType>
418 TrackCloudType&
cloud,
429 const scalar np0 = this->nParticle_;
430 const scalar d0 = this->d_;
431 const vector& U0 = this->U_;
432 const scalar
T0 = this->T_;
433 const scalar mass0 = this->mass();
434 const scalar
rho0 = this->rho_;
438 scalar Ts, rhos, mus, Prs, kappas;
439 this->calcSurfaceValues(
cloud, td,
T0, Ts, rhos, mus, Prs, kappas);
440 scalar Res = this->
Re(rhos, U0, td.Uc(), d0, mus);
462 scalar dhsTrans = 0.0;
513 scalar mass1 = updateMassFraction(mass0, dMass, Y_);
519 cloud.constProps().volUpdateType()
520 == constantProperties::volumeUpdateType::mUndefined
524 if (
cloud.constProps().constantVolume())
526 this->rho_ = mass1/this->
volume();
530 this->d_ =
cbrt(mass1/this->rho_*6/
pi);
535 switch (
cloud.constProps().volUpdateType())
537 case constantProperties::volumeUpdateType::mConstRho :
539 this->d_ =
cbrt(mass1/this->rho_*6/
pi);
542 case constantProperties::volumeUpdateType::mConstVol :
544 this->rho_ = mass1/this->
volume();
547 case constantProperties::volumeUpdateType::mUpdateRhoAndVol :
557 this->rho_ = mass1/(this->
volume() - deltaVol);
558 this->d_ =
cbrt(mass1/this->rho_*6/
pi);
566 if (np0*mass1 <
cloud.constProps().minParcelMass())
568 td.keepParticle =
false;
570 if (
cloud.solution().coupled())
572 scalar dm = np0*mass0;
577 scalar dmi = dm*Y_[i];
581 cloud.rhoTrans(gid)[this->
cell()] += dmi;
586 cloud.phaseChange().addToPhaseChangeMass(np0*mass1);
593 correctSurfaceValues(
cloud, td, Ts,
Cs, rhos, mus, Prs, kappas);
594 Res = this->
Re(rhos, U0, td.Uc(), this->d(), mus);
605 this->calcHeatTransfer
627 this->calcVelocity(
cloud, td, dt, Res, mus, mass1,
Su, dUTrans, Spu);
633 if (
cloud.solution().coupled())
638 scalar dm = np0*dMass[i];
648 cloud.UTrans()[this->
cell()] += np0*dUTrans;
652 cloud.hsTrans()[this->
cell()] += np0*dhsTrans;
653 cloud.hsCoeff()[this->
cell()] += np0*Sph;
656 if (
cloud.radiation())
658 const scalar ap = this->areaP();
659 const scalar T4 =
pow4(
T0);
660 cloud.radAreaP()[this->
cell()] += dt*np0*ap;
661 cloud.radT4()[this->
cell()] += dt*np0*T4;
662 cloud.radAreaPT4()[this->
cell()] += dt*np0*ap*T4;
int debug
Static debugging option.
const scalar RR
Universal gas constant: default in [J/(kmol K)].
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const interpolation< scalar > & pInterp() const
Return const access to the interpolator for continuous phase.
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package,...
static constexpr const zero Zero
Global zero (0)
const speciesTable & species() const
Return the table of species.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
void calc(TrackCloudType &cloud, trackingData &td, const scalar dt)
Update parcel properties over the time interval.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const wordList volume
Standard volume field types (scalar, vector, tensor, etc)
basicSpecieMixture & composition
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
scalar pc() const
Return the continuous phase pressure.
virtual scalar Hs(const label speciei, const scalar p, const scalar T) const =0
Sensible enthalpy [J/kg].
#define forAll(list, i)
Loop across all elements in list.
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 scalar rho, const label idPhase, const scalar YPhase, const scalarField &YLiq, const scalarField &YSol, scalarField &dMassPC, scalar &Sh, scalar &N, scalar &NCpW, scalarField &Cs)
Calculate Phase change.
dimensionedScalar pow4(const dimensionedScalar &ds)
scalar updatedDeltaVolume(TrackCloudType &cloud, const scalarField &dMass, const scalar p, const scalar T)
Return change of volume due to mass exchange.
void setCellValues(TrackCloudType &cloud, trackingData &td)
Set cell values.
PtrList< coordinateSystem > coordinates(solidRegions.size())
dimensionedScalar Pr("Pr", dimless, laminarTransport)
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
ReactingParcel(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from mesh, coordinates and topology.
void cellValueSourceCorrection(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct cell values using latest transfer information.
PtrList< volScalarField > & Y
A cloud is a registry collection of lagrangian particles.
constexpr scalar pi(M_PI)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalarField Re(const UList< complex > &cf)
Extract real component.
const volScalarField & Cp
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const Vector< label > N(dict.get< Vector< label >>("N"))
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.
dimensionedScalar cbrt(const dimensionedScalar &ds)
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.
Reacting parcel class with one/two-way coupling with the continuous phase.
#define WarningInFunction
Report a warning using Foam::Warning.
A cell is defined as a list of faces with extra functionality.
scalar updateMassFraction(const scalar mass0, const scalarField &dMass, scalarField &Y) const
Update mass fraction.
virtual scalar Cp(const label speciei, const scalar p, const scalar T) const =0
Heat capacity at constant pressure [J/(kg K)].