Commit 402b8637 authored by Will Bainbridge's avatar Will Bainbridge Committed by Andrew Heather
Browse files

KinematicParcel: Removed continuous phase data

Interpolated continuous phase data is only needed during a track and
therefore shouldn't be stored on the parcel. The continuous velocity,
density and viscosity have been moved from the kinematic parcel to the
kinematic parcel tracking data. This reduces the memory usage of the
kinematic layer by about one third. The thermo and reacting layers still
require the same treatment.
parent 30c2f8fa
......@@ -48,9 +48,9 @@ void Foam::KinematicParcel<ParcelType>::setCellValues
{
tetIndices tetIs = this->currentTetIndices();
rhoc_ = td.rhoInterp().interpolate(this->coordinates(), tetIs);
td.rhoc() = td.rhoInterp().interpolate(this->coordinates(), tetIs);
if (rhoc_ < cloud.constProps().rhoMin())
if (td.rhoc() < cloud.constProps().rhoMin())
{
if (debug)
{
......@@ -59,20 +59,20 @@ void Foam::KinematicParcel<ParcelType>::setCellValues
<< " to " << cloud.constProps().rhoMin() << nl << endl;
}
rhoc_ = cloud.constProps().rhoMin();
td.rhoc() = cloud.constProps().rhoMin();
}
Uc_ = td.UInterp().interpolate(this->coordinates(), tetIs);
td.Uc() = td.UInterp().interpolate(this->coordinates(), tetIs);
muc_ = td.muInterp().interpolate(this->coordinates(), tetIs);
td.muc() = td.muInterp().interpolate(this->coordinates(), tetIs);
// Apply dispersion components to carrier phase velocity
Uc_ = cloud.dispersion().update
td.Uc() = cloud.dispersion().update
(
dt,
this->cell(),
U_,
Uc_,
td.Uc(),
UTurb_,
tTurb_
);
......@@ -88,7 +88,7 @@ void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection
const scalar dt
)
{
Uc_ += cloud.UTrans()[this->cell()]/massCell(this->cell());
td.Uc() += cloud.UTrans()[this->cell()]/massCell(td);
}
......@@ -107,7 +107,7 @@ void Foam::KinematicParcel<ParcelType>::calc
const scalar mass0 = mass();
// Reynolds number
const scalar Re = this->Re(U_, d_, rhoc_, muc_);
const scalar Re = this->Re(td);
// Sources
......@@ -127,7 +127,8 @@ void Foam::KinematicParcel<ParcelType>::calc
// ~~~~~~
// Calculate new particle velocity
this->U_ = calcVelocity(cloud, td, dt, Re, muc_, mass0, Su, dUTrans, Spu);
this->U_ =
calcVelocity(cloud, td, dt, Re, td.muc(), mass0, Su, dUTrans, Spu);
// Accumulate carrier phase source terms
......@@ -158,24 +159,25 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
scalar& Spu
) const
{
typedef typename TrackCloudType::parcelType parcelType;
typedef typename TrackCloudType::forceType forceType;
const typename TrackCloudType::parcelType& p =
static_cast<const typename TrackCloudType::parcelType&>(*this);
typename TrackCloudType::parcelType::trackingData& ttd =
static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
const forceType& forces = cloud.forces();
const typename TrackCloudType::forceType& forces = cloud.forces();
// Momentum source due to particle forces
const parcelType& p = static_cast<const parcelType&>(*this);
const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, mu);
const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, mu);
const forceSuSp Feff = Fcp + Fncp;
const scalar massEff = forces.massEff(p, mass);
const scalar massEff = forces.massEff(p, ttd, mass);
// New particle velocity
//~~~~~~~~~~~~~~~~~~~~~~
// Update velocity - treat as 3-D
const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff;
const vector abp = (Feff.Sp()*td.Uc() + (Feff.Su() + Su))/massEff;
const scalar bp = Feff.Sp()/massEff;
Spu = dt*Feff.Sp();
......@@ -186,7 +188,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
vector Unew = Ures.value();
// note: Feff.Sp() and Fc.Sp() must be the same
dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su());
dUTrans += dt*(Feff.Sp()*(Ures.average() - td.Uc()) - Fcp.Su());
// Apply correction to velocity and dUTrans for reduced-D cases
const polyMesh& mesh = cloud.pMesh();
......@@ -215,10 +217,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
rho_(p.rho_),
age_(p.age_),
tTurb_(p.tTurb_),
UTurb_(p.UTurb_),
rhoc_(p.rhoc_),
Uc_(p.Uc_),
muc_(p.muc_)
UTurb_(p.UTurb_)
{}
......@@ -239,10 +238,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
rho_(p.rho_),
age_(p.age_),
tTurb_(p.tTurb_),
UTurb_(p.UTurb_),
rhoc_(p.rhoc_),
Uc_(p.Uc_),
muc_(p.muc_)
UTurb_(p.UTurb_)
{}
......@@ -259,14 +255,12 @@ bool Foam::KinematicParcel<ParcelType>::move
{
typename TrackCloudType::parcelType& p =
static_cast<typename TrackCloudType::parcelType&>(*this);
typename TrackCloudType::particleType::trackingData& ttd =
static_cast<typename TrackCloudType::particleType::trackingData&>(td);
typename TrackCloudType::parcelType::trackingData& ttd =
static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
ttd.switchProcessor = false;
ttd.keepParticle = true;
const polyMesh& mesh = cloud.pMesh();
const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
const cloudSolution& solution = cloud.solution();
const scalarField& cellLengthScale = cloud.cellLengthScale();
const scalar maxCo = solution.maxCo();
......
......@@ -184,6 +184,18 @@ public:
autoPtr<interpolation<scalar>> muInterp_;
// Cached continuous phase properties
//- Density [kg/m3]
scalar rhoc_;
//- Velocity [m/s]
vector Uc_;
//- Viscosity [Pa.s]
scalar muc_;
//- Local gravitational or other body-force acceleration
const vector& g_;
......@@ -219,6 +231,24 @@ public:
// phase dynamic viscosity field
inline const interpolation<scalar>& muInterp() const;
//- Return the continuous phase density
inline scalar rhoc() const;
//- Access the continuous phase density
inline scalar& rhoc();
//- Return the continuous phase velocity
inline const vector& Uc() const;
//- Access the continuous phase velocity
inline vector& Uc();
//- Return the continuous phase viscosity
inline scalar muc() const;
//- Access the continuous phase viscosity
inline scalar& muc();
// Return const access to the gravitational acceleration vector
inline const vector& g() const;
......@@ -268,18 +298,6 @@ protected:
vector UTurb_;
// Cell-based quantities
//- Density [kg/m3]
scalar rhoc_;
//- Velocity [m/s]
vector Uc_;
//- Viscosity [Pa.s]
scalar muc_;
// Protected Member Functions
//- Calculate new particle velocity
......@@ -459,15 +477,6 @@ public:
//- Return const access to turbulent velocity fluctuation
inline const vector& UTurb() const;
//- Return const access to carrier density [kg/m3]
inline scalar rhoc() const;
//- Return const access to carrier velocity [m/s]
inline const vector& Uc() const;
//- Return const access to carrier viscosity [Pa.s]
inline scalar muc() const;
// Edit
......@@ -505,7 +514,7 @@ public:
// Helper functions
//- Cell owner mass
inline scalar massCell(const label celli) const;
inline scalar massCell(const trackingData& td) const;
//- Particle mass
inline scalar mass() const;
......@@ -532,31 +541,53 @@ public:
inline static scalar areaS(const scalar d);
//- Reynolds number
inline scalar Re
inline scalar Re(const trackingData& td) const;
//- Reynolds number for given conditions
inline static scalar Re
(
const vector& U, // particle velocity
const scalar d, // particle diameter
const scalar rhoc, // carrier density
const scalar muc // carrier dynamic viscosity
) const;
const scalar rhoc,
const vector& U,
const vector& Uc,
const scalar d,
const scalar muc
);
//- Weber number
inline scalar We
(
const vector& U, // particle velocity
const scalar d, // particle diameter
const scalar rhoc, // carrier density
const scalar sigma // particle surface tension
const trackingData& td,
const scalar sigma
) const;
//- Weber number for given conditions
inline static scalar We
(
const scalar rhoc,
const vector& U,
const vector& Uc,
const scalar d,
const scalar sigma
);
//- Eotvos number
inline scalar Eo
(
const vector& a, // acceleration
const scalar d, // particle diameter
const scalar sigma // particle surface tension
const trackingData& td,
const scalar sigma
) const;
//- Eotvos number for given conditions
inline static scalar Eo
(
const vector& g,
const scalar rho,
const scalar rhoc,
const vector& U,
const scalar d,
const scalar sigma
);
// Main calculation loop
......
......@@ -89,10 +89,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
rho_(0.0),
age_(0.0),
tTurb_(0.0),
UTurb_(Zero),
rhoc_(0.0),
Uc_(Zero),
muc_(0.0)
UTurb_(Zero)
{}
......@@ -114,10 +111,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
rho_(0.0),
age_(0.0),
tTurb_(0.0),
UTurb_(Zero),
rhoc_(0.0),
Uc_(Zero),
muc_(0.0)
UTurb_(Zero)
{}
......@@ -147,10 +141,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
rho_(constProps.rho0()),
age_(0.0),
tTurb_(0.0),
UTurb_(Zero),
rhoc_(0.0),
Uc_(Zero),
muc_(0.0)
UTurb_(Zero)
{}
......@@ -269,28 +260,7 @@ inline const Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb() const
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::rhoc() const
{
return rhoc_;
}
template<class ParcelType>
inline const Foam::vector& Foam::KinematicParcel<ParcelType>::Uc() const
{
return Uc_;
}
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::muc() const
{
return muc_;
}
template<class ParcelType>
inline void Foam::KinematicParcel<ParcelType>::active(const bool state)
inline bool& Foam::KinematicParcel<ParcelType>::active()
{
active_ = state;
}
......@@ -362,10 +332,10 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::massCell
(
const label celli
const trackingData& td
) const
{
return rhoc_*this->mesh().cellVolumes()[celli];
return td.rhoc()*this->mesh().cellVolumes()[this->cell()];
}
......@@ -428,39 +398,76 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::areaS(const scalar d)
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
(
const trackingData& td
) const
{
return Re(td.rhoc(), U_, td.Uc(), d_, td.muc());
}
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
(
const scalar rhoc,
const vector& U,
const vector& Uc,
const scalar d,
const scalar rhoc,
const scalar muc
)
{
return rhoc*mag(U - Uc)*d/max(muc, ROOTVSMALL);
}
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::We
(
const trackingData& td,
const scalar sigma
) const
{
return rhoc*mag(U - Uc_)*d/(muc + ROOTVSMALL);
return We(td.rhoc(), U_, td.Uc(), d_, sigma);
}
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::We
(
const scalar rhoc,
const vector& U,
const vector& Uc,
const scalar d,
const scalar rhoc,
const scalar sigma
)
{
return rhoc*magSqr(U - Uc)*d/max(sigma, ROOTVSMALL);
}
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::Eo
(
const trackingData& td,
const scalar sigma
) const
{
return rhoc*magSqr(U - Uc_)*d/(sigma + ROOTVSMALL);
return Eo(td.g(), rho_, td.rhoc(), U_, d_, sigma);
}
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::Eo
(
const vector& a,
const vector& g,
const scalar rho,
const scalar rhoc,
const vector& U,
const scalar d,
const scalar sigma
) const
)
{
vector dir = U_/(mag(U_) + ROOTVSMALL);
return mag(a & dir)*(rho_ - rhoc_)*sqr(d)/(sigma + ROOTVSMALL);
const vector dir = U/max(mag(U), ROOTVSMALL);
return mag(g & dir)*(rho - rhoc)*sqr(d)/max(sigma, ROOTVSMALL);
}
......
......@@ -41,7 +41,7 @@ Foam::string Foam::KinematicParcel<ParcelType>::propertyTypes_ =
template<class ParcelType>
const std::size_t Foam::KinematicParcel<ParcelType>::sizeofFields
(
offsetof(KinematicParcel<ParcelType>, rhoc_)
sizeof(KinematicParcel<ParcelType>)
- offsetof(KinematicParcel<ParcelType>, active_)
);
......@@ -66,10 +66,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
rho_(0.0),
age_(0.0),
tTurb_(0.0),
UTurb_(Zero),
rhoc_(0.0),
Uc_(Zero),
muc_(0.0)
UTurb_(Zero)
{
if (readFields)
{
......
......@@ -56,6 +56,9 @@ inline Foam::KinematicParcel<ParcelType>::trackingData::trackingData
cloud.mu()
)
),
rhoc_(Zero),
Uc_(Zero),
muc_(Zero),
g_(cloud.g().value()),
part_(part)
{}
......@@ -93,6 +96,50 @@ Foam::KinematicParcel<ParcelType>::trackingData::g() const
}
template<class ParcelType>
inline Foam::scalar
Foam::KinematicParcel<ParcelType>::trackingData::rhoc() const
{
return rhoc_;
}
template<class ParcelType>
inline Foam::scalar& Foam::KinematicParcel<ParcelType>::trackingData::rhoc()
{
return rhoc_;
}
template<class ParcelType>
inline const Foam::vector&
Foam::KinematicParcel<ParcelType>::trackingData::Uc() const
{
return Uc_;
}
template<class ParcelType>
inline Foam::vector& Foam::KinematicParcel<ParcelType>::trackingData::Uc()
{
return Uc_;
}
template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::trackingData::muc() const
{
return muc_;
}
template<class ParcelType>
inline Foam::scalar& Foam::KinematicParcel<ParcelType>::trackingData::muc()
{
return muc_;
}
template<class ParcelType>
inline typename Foam::KinematicParcel<ParcelType>::trackingData::trackPart
Foam::KinematicParcel<ParcelType>::trackingData::part() const
......
......@@ -198,7 +198,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Calc surface values
scalar Ts, rhos, mus, Prs, kappas;
this->calcSurfaceValues(cloud, td, T0, Ts, rhos, mus, Prs, kappas);
scalar Res = this->Re(U0, d0, rhos, mus);
scalar Res = this->Re(rhos, U0, td.Uc(), d0, mus);
// Sources
......@@ -390,7 +390,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Correct surface values due to emitted species
this->correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
Res = this->Re(U0, this->d_, rhos, mus);
Res = this->Re(rhos, U0, td.Uc(), this->d_, mus);
// 3. Compute heat- and momentum transfers
......@@ -572,7 +572,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
{
// Molar average molecular weight of carrier mix
const scalar Wc =
max(SMALL, this->rhoc_*RR*this->Tc_/this->pc_);
max(SMALL, td.rhoc()*RR*this->Tc_/this->pc_);
// Note: hardcoded gaseous diffusivities for now
// TODO: add to carrier thermo
......@@ -649,7 +649,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
T,
this->Tc_,
this->pc_,
this->rhoc_,
td.rhoc(),