Commit 95018c68 authored by Will Bainbridge's avatar Will Bainbridge Committed by Andrew Heather
Browse files

ThermoParcel, ReactingParcel: Removed continuous phase data

parent 3e85a5ab
......@@ -212,7 +212,7 @@ void Foam::KinematicCloud<CloudType>::evolveCloud
// before it is required.
cloud.motion(cloud, td);
stochasticCollision().update(solution_.trackTime());
stochasticCollision().update(td, solution_.trackTime());
}
else
{
......
......@@ -221,7 +221,6 @@ void Foam::ReactingCloud<CloudType>::setParcelThermoProperties
{
CloudType::setParcelThermoProperties(parcel, lagrangianDt);
parcel.pc() = this->thermo().thermo().p()[parcel.cell()];
parcel.Y() = composition().YMixture0();
}
......
......@@ -42,8 +42,7 @@ template<class TrackCloudType>
void Foam::KinematicParcel<ParcelType>::setCellValues
(
TrackCloudType& cloud,
trackingData& td,
const scalar dt
trackingData& td
)
{
tetIndices tetIs = this->currentTetIndices();
......@@ -65,8 +64,18 @@ void Foam::KinematicParcel<ParcelType>::setCellValues
td.Uc() = td.UInterp().interpolate(this->coordinates(), tetIs);
td.muc() = td.muInterp().interpolate(this->coordinates(), tetIs);
}
// Apply dispersion components to carrier phase velocity
template<class ParcelType>
template<class TrackCloudType>
void Foam::KinematicParcel<ParcelType>::calcDispersion
(
TrackCloudType& cloud,
trackingData& td,
const scalar dt
)
{
td.Uc() = cloud.dispersion().update
(
dt,
......@@ -307,7 +316,9 @@ bool Foam::KinematicParcel<ParcelType>::move
if (dt > ROOTVSMALL)
{
// Update cell based properties
p.setCellValues(cloud, ttd, dt);
p.setCellValues(cloud, ttd);
p.calcDispersion(cloud, ttd, dt);
if (solution.cellValueSourceCorrection())
{
......
......@@ -593,7 +593,12 @@ public:
//- Set cell values
template<class TrackCloudType>
void setCellValues
void setCellValues(TrackCloudType& cloud, trackingData& td);
//- Apply dispersion to the carrier phase velocity and update
// parcel turbulence parameters
template<class TrackCloudType>
void calcDispersion
(
TrackCloudType& cloud,
trackingData& td,
......
......@@ -143,11 +143,10 @@ template<class TrackCloudType>
void Foam::ReactingMultiphaseParcel<ParcelType>::setCellValues
(
TrackCloudType& cloud,
trackingData& td,
const scalar dt
trackingData& td
)
{
ParcelType::setCellValues(cloud, td, dt);
ParcelType::setCellValues(cloud, td);
}
......@@ -187,7 +186,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const scalar T0 = this->T_;
const scalar mass0 = this->mass();
const scalar pc = this->pc_;
const scalar pc = td.pc();
const scalarField& YMix = this->Y_;
const label idG = composition.idGas();
......@@ -571,8 +570,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
if (cloud.heatTransfer().BirdCorrection())
{
// Molar average molecular weight of carrier mix
const scalar Wc =
max(SMALL, td.rhoc()*RR*this->Tc_/this->pc_);
const scalar Wc = max(SMALL, td.rhoc()*RR*td.Tc()/td.pc());
// Note: hardcoded gaseous diffusivities for now
// TODO: add to carrier thermo
......@@ -581,7 +579,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
forAll(dMassDV, i)
{
const label id = composition.localToCarrierId(GAS, i);
const scalar Cp = composition.carrier().Cp(id, this->pc_, Ts);
const scalar Cp = composition.carrier().Cp(id, td.pc(), Ts);
const scalar W = composition.carrier().W(id);
const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
......@@ -589,7 +587,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
const scalar Dab =
3.6059e-3*(pow(1.8*Ts, 1.75))
*sqrt(1.0/W + 1.0/Wc)
/(this->pc_*beta);
/(td.pc()*beta);
N += Ni;
NCpW += Ni*Cp*W;
......@@ -647,8 +645,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
this->cell(),
d,
T,
this->Tc_,
this->pc_,
td.Tc(),
td.pc(),
td.rhoc(),
mass,
YGas,
......
......@@ -421,12 +421,7 @@ public:
//- Set cell values
template<class TrackCloudType>
void setCellValues
(
TrackCloudType& cloud,
trackingData& td,
const scalar dt
);
void setCellValues(TrackCloudType& cloud, trackingData& td);
//- Correct cell values using latest transfer information
template<class TrackCloudType>
......
......@@ -74,7 +74,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
return;
}
const scalar TMax = phaseChange.TMax(pc_, X);
const scalar TMax = phaseChange.TMax(td.pc(), X);
const scalar Tdash = min(T, TMax);
const scalar Tsdash = min(Ts, TMax);
......@@ -89,8 +89,8 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
nus,
Tdash,
Tsdash,
pc_,
this->Tc_,
td.pc(),
td.Tc(),
X,
dMassPC
);
......@@ -107,7 +107,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
{
const label cid = composition.localToCarrierId(idPhase, i);
const scalar dh = phaseChange.dh(cid, i, pc_, Tdash);
const scalar dh = phaseChange.dh(cid, i, td.pc(), Tdash);
Sh -= dMassPC[i]*dh/dt;
}
......@@ -116,18 +116,18 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
if (cloud.heatTransfer().BirdCorrection())
{
// Average molecular weight of carrier mix - assumes perfect gas
const scalar Wc = td.rhoc()*RR*this->Tc_/this->pc_;
const scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
forAll(dMassPC, i)
{
const label cid = composition.localToCarrierId(idPhase, i);
const scalar Cp = composition.carrier().Cp(cid, pc_, Tsdash);
const scalar Cp = composition.carrier().Cp(cid, td.pc(), Tsdash);
const scalar W = composition.carrier().W(cid);
const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
const scalar Dab =
composition.liquids().properties()[i].D(pc_, Tsdash, Wc);
composition.liquids().properties()[i].D(td.pc(), Tsdash, Wc);
// Molar flux of species coming from the particle (kmol/m^2/s)
N += Ni;
......@@ -175,8 +175,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
:
ParcelType(p),
mass0_(p.mass0_),
Y_(p.Y_),
pc_(p.pc_)
Y_(p.Y_)
{}
......@@ -189,8 +188,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
:
ParcelType(p, mesh),
mass0_(p.mass0_),
Y_(p.Y_),
pc_(p.pc_)
Y_(p.Y_)
{}
......@@ -201,19 +199,18 @@ template<class TrackCloudType>
void Foam::ReactingParcel<ParcelType>::setCellValues
(
TrackCloudType& cloud,
trackingData& td,
const scalar dt
trackingData& td
)
{
ParcelType::setCellValues(cloud, td, dt);
ParcelType::setCellValues(cloud, td);
pc_ = td.pInterp().interpolate
td.pc() = td.pInterp().interpolate
(
this->coordinates(),
this->currentTetIndices()
);
if (pc_ < cloud.constProps().pMin())
if (td.pc() < cloud.constProps().pMin())
{
if (debug)
{
......@@ -222,7 +219,7 @@ void Foam::ReactingParcel<ParcelType>::setCellValues
<< " to " << cloud.constProps().pMin() << nl << endl;
}
pc_ = cloud.constProps().pMin();
td.pc() = cloud.constProps().pMin();
}
}
......@@ -261,20 +258,15 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
forAll(cloud.rhoTrans(), i)
{
scalar Y = cloud.rhoTrans(i)[this->cell()]/addedMass;
CpEff += Y*cloud.composition().carrier().Cp
(
i,
this->pc_,
this->Tc_
);
CpEff += Y*cloud.composition().carrier().Cp(i, td.pc(), td.Tc());
}
const scalar Cpc = td.CpInterp().psi()[this->cell()];
this->Cpc_ = (massCell*Cpc + addedMass*CpEff)/massCellNew;
td.Cpc() = (massCell*Cpc + addedMass*CpEff)/massCellNew;
this->Tc_ += cloud.hsTrans()[this->cell()]/(this->Cpc_*massCellNew);
td.Tc() += cloud.hsTrans()[this->cell()]/(td.Cpc()*massCellNew);
if (this->Tc_ < cloud.constProps().TMin())
if (td.Tc() < cloud.constProps().TMin())
{
if (debug)
{
......@@ -283,7 +275,7 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
<< " to " << cloud.constProps().TMin() << nl << endl;
}
this->Tc_ = cloud.constProps().TMin();
td.Tc() = cloud.constProps().TMin();
}
}
......@@ -320,10 +312,10 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
Xinf /= sum(Xinf);
// Molar fraction of far field species at particle surface
const scalar Xsff = 1.0 - min(sum(Cs)*RR*this->T_/pc_, 1.0);
const scalar Xsff = 1.0 - min(sum(Cs)*RR*this->T_/td.pc(), 1.0);
// Surface carrier total molar concentration
const scalar CsTot = pc_/(RR*this->T_);
const scalar CsTot = td.pc()/(RR*this->T_);
// Surface carrier composition (molar fraction)
scalarField Xs(Xinf.size());
......@@ -357,9 +349,9 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
const scalar cbrtW = cbrt(W);
rhos += Xs[i]*W;
mus += Ys[i]*sqrtW*thermo.carrier().mu(i, pc_, T);
kappas += Ys[i]*cbrtW*thermo.carrier().kappa(i, pc_, T);
Cps += Xs[i]*thermo.carrier().Cp(i, pc_, T);
mus += Ys[i]*sqrtW*thermo.carrier().mu(i, td.pc(), T);
kappas += Ys[i]*cbrtW*thermo.carrier().kappa(i, td.pc(), T);
Cps += Xs[i]*thermo.carrier().Cp(i, td.pc(), T);
sumYiSqrtW += Ys[i]*sqrtW;
sumYiCbrtW += Ys[i]*cbrtW;
......@@ -367,7 +359,7 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
Cps = max(Cps, ROOTVSMALL);
rhos *= pc_/(RR*T);
rhos *= td.pc()/(RR*T);
rhos = max(rhos, ROOTVSMALL);
mus /= sumYiSqrtW;
......@@ -478,7 +470,7 @@ void Foam::ReactingParcel<ParcelType>::calc
scalarField dMass(dMassPC);
scalar mass1 = updateMassFraction(mass0, dMass, Y_);
this->Cp_ = composition.Cp(0, Y_, pc_, T0);
this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
// Update particle density or diameter
if (cloud.constProps().constantVolume())
......@@ -504,7 +496,7 @@ void Foam::ReactingParcel<ParcelType>::calc
{
scalar dmi = dm*Y_[i];
label gid = composition.localToCarrierId(0, i);
scalar hs = composition.carrier().Hs(gid, pc_, T0);
scalar hs = composition.carrier().Hs(gid, td.pc(), T0);
cloud.rhoTrans(gid)[this->cell()] += dmi;
cloud.hsTrans()[this->cell()] += dmi*hs;
......@@ -544,7 +536,7 @@ void Foam::ReactingParcel<ParcelType>::calc
Sph
);
this->Cp_ = composition.Cp(0, Y_, pc_, T0);
this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
// Motion
......@@ -565,7 +557,7 @@ void Foam::ReactingParcel<ParcelType>::calc
{
scalar dm = np0*dMass[i];
label gid = composition.localToCarrierId(0, i);
scalar hs = composition.carrier().Hs(gid, pc_, T0);
scalar hs = composition.carrier().Hs(gid, td.pc(), T0);
cloud.rhoTrans(gid)[this->cell()] += dm;
cloud.UTrans()[this->cell()] += dm*U0;
......
......@@ -128,6 +128,12 @@ public:
autoPtr<interpolation<scalar>> pInterp_;
// Cached continuous phase properties
//- Pressure [Pa]
scalar pc_;
public:
typedef typename ParcelType::trackingData::trackPart trackPart;
......@@ -148,6 +154,12 @@ public:
//- Return const access to the interpolator for continuous phase
// pressure field
inline const interpolation<scalar>& pInterp() const;
//- Return the continuous phase pressure
inline scalar pc() const;
//- Access the continuous phase pressure
inline scalar& pc();
};
......@@ -164,12 +176,6 @@ protected:
scalarField Y_;
// Cell-based quantities
//- Pressure [Pa]
scalar pc_;
// Protected Member Functions
//- Calculate Phase change
......@@ -336,12 +342,6 @@ public:
//- Return const access to mass fractions of mixture []
inline const scalarField& Y() const;
//- Return the owner cell pressure [Pa]
inline scalar pc() const;
//- Return reference to the owner cell pressure [Pa]
inline scalar& pc();
// Edit
......@@ -356,12 +356,7 @@ public:
//- Set cell values
template<class TrackCloudType>
void setCellValues
(
TrackCloudType& cloud,
trackingData& td,
const scalar dt
);
void setCellValues(TrackCloudType& cloud, trackingData& td);
//- Correct cell values using latest transfer information
template<class TrackCloudType>
......
......@@ -71,8 +71,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
:
ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
mass0_(0.0),
Y_(0),
pc_(0.0)
Y_(0)
{}
......@@ -86,8 +85,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
:
ParcelType(mesh, position, celli),
mass0_(0.0),
Y_(0),
pc_(0.0)
Y_(0)
{}
......@@ -129,8 +127,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
constProps
),
mass0_(0.0),
Y_(Y0),
pc_(0.0)
Y_(Y0)
{
// Set initial parcel mass
mass0_ = this->mass();
......@@ -171,20 +168,6 @@ inline const Foam::scalarField& Foam::ReactingParcel<ParcelType>::Y() const
}
template<class ParcelType>
inline Foam::scalar Foam::ReactingParcel<ParcelType>::pc() const
{
return pc_;
}
template<class ParcelType>
inline Foam::scalar& Foam::ReactingParcel<ParcelType>::pc()
{
return pc_;
}
template<class ParcelType>
inline Foam::scalar& Foam::ReactingParcel<ParcelType>::mass0()
{
......
......@@ -55,8 +55,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
:
ParcelType(mesh, is, readFields),
mass0_(0.0),
Y_(0),
pc_(0.0)
Y_(0)
{
if (readFields)
{
......
......@@ -39,7 +39,8 @@ inline Foam::ReactingParcel<ParcelType>::trackingData::trackingData
cloud.solution().interpolationSchemes(),
cloud.p()
)
)
),
pc_(Zero)
{}
......@@ -51,4 +52,18 @@ Foam::ReactingParcel<ParcelType>::trackingData::pInterp() const
}
template<class ParcelType>
inline Foam::scalar Foam::ReactingParcel<ParcelType>::trackingData::pc() const
{
return pc_;
}
template<class ParcelType>
inline Foam::scalar& Foam::ReactingParcel<ParcelType>::trackingData::pc()
{
return pc_;
}
// ************************************************************************* //
......@@ -35,19 +35,18 @@ template<class TrackCloudType>
void Foam::ThermoParcel<ParcelType>::setCellValues
(
TrackCloudType& cloud,
trackingData& td,
const scalar dt
trackingData& td
)
{
ParcelType::setCellValues(cloud, td, dt);
ParcelType::setCellValues(cloud, td);
tetIndices tetIs = this->currentTetIndices();
Cpc_ = td.CpInterp().interpolate(this->coordinates(), tetIs);
td.Cpc() = td.CpInterp().interpolate(this->coordinates(), tetIs);
Tc_ = td.TInterp().interpolate(this->coordinates(), tetIs);
td.Tc() = td.TInterp().interpolate(this->coordinates(), tetIs);
if (Tc_ < cloud.constProps().TMin())
if (td.Tc() < cloud.constProps().TMin())
{
if (debug)
{
......@@ -56,7 +55,7 @@ void Foam::ThermoParcel<ParcelType>::setCellValues
<< " to " << cloud.constProps().TMin() << nl << endl;
}
Tc_ = cloud.constProps().TMin();
td.Tc() = cloud.constProps().TMin();
}
}
......@@ -71,7 +70,7 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
)
{
const label celli = this->cell();
const scalar massCell = this->massCell(celli);
const scalar massCell = this->massCell(td);
td.Uc() += cloud.UTrans()[celli]/massCell;
......@@ -79,9 +78,9 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
// Tc_ = td.TInterp().interpolate(this->coordinates(), tetIs);
const scalar CpMean = td.CpInterp().psi()[celli];
Tc_ += cloud.hsTrans()[celli]/(CpMean*massCell);
td.Tc() += cloud.hsTrans()[celli]/(CpMean*massCell);
if (Tc_ < cloud.constProps().TMin())
if (td.Tc() < cloud.constProps().TMin())
{
if (debug)
{
......@@ -90,7 +89,7 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
<< " to " << cloud.constProps().TMin() << nl << endl;
}
Tc_ = cloud.constProps().TMin();
td.Tc() = cloud.constProps().TMin();
}
}
......@@ -110,7 +109,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
) const
{
// Surface temperature using two thirds rule
Ts = (2.0*T + Tc_)/3.0;
Ts = (2.0*T + td.Tc())/3.0;
if (Ts < cloud.constProps().TMin())
{
......@@ -125,7 +124,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues