Commit 0e2b526c authored by Andrew Heather's avatar Andrew Heather
Browse files

more updates

parent 7677c87f
......@@ -216,6 +216,11 @@ public:
//- Return particle properties dictionary
inline const IOdictionary& particleProperties() const;
//- Return the constant properties
inline const typename ParcelType::constantProperties&
constProps() const;
//- Return coupled flag
inline const Switch coupled() const;
......
......@@ -55,6 +55,14 @@ Foam::KinematicCloud<ParcelType>::particleProperties() const
}
template<class ParcelType>
inline const typename ParcelType::constantProperties&
Foam::KinematicCloud<ParcelType>::constProps() const
{
return constProps_;
}
template<class ParcelType>
inline const Foam::Switch Foam::KinematicCloud<ParcelType>::coupled() const
{
......
......@@ -143,6 +143,10 @@ public:
// Access
//- Return the constant properties
inline const typename ParcelType::constantProperties&
constProps() const;
//- Return const access to carrier phase thermo package
inline const hCombustionThermo& carrierThermo() const;
......
......@@ -26,6 +26,14 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ParcelType>
inline const typename ParcelType::constantProperties&
Foam::ReactingCloud<ParcelType>::constProps() const
{
return constProps_;
}
template<class ParcelType>
inline const Foam::hCombustionThermo&
Foam::ReactingCloud<ParcelType>::carrierThermo() const
......
......@@ -137,6 +137,11 @@ public:
// Access
//- Return the constant properties
inline const typename ParcelType::constantProperties&
constProps() const;
// Sub-models
//- Return reference to devolatilisation model
......
......@@ -26,6 +26,14 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ParcelType>
inline const typename ParcelType::constantProperties&
Foam::ReactingMultiphaseCloud<ParcelType>::constProps() const
{
return constProps_;
}
template<class ParcelType>
inline const Foam::DevolatilisationModel
<
......
......@@ -142,6 +142,10 @@ public:
// Access
//- Return the constant properties
inline const typename ParcelType::constantProperties&
constProps() const;
//- Return const access to thermo package
inline const basicThermo& carrierThermo() const;
......
......@@ -28,6 +28,14 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ParcelType>
inline const typename ParcelType::constantProperties&
Foam::ThermoCloud<ParcelType>::constProps() const
{
return constProps_;
}
template<class ParcelType>
inline const Foam::basicThermo&
Foam::ThermoCloud<ParcelType>::carrierThermo() const
......
......@@ -64,20 +64,19 @@ void Foam::KinematicParcel<ParcelType>::calc
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U, dUTrans
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 1. Calculate velocity - update U, dUTrans
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
vector dUTrans = vector::zero;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 2. Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().coupled())
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Update momentum transfer
td.cloud().UTrans()[cellI] += nParticle_*dUTrans;
......@@ -86,9 +85,9 @@ void Foam::KinematicParcel<ParcelType>::calc
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
this->U() = U1;
}
......
......@@ -26,6 +26,40 @@ License
#include "ReactingMultiphaseParcel.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::cpEff
(
TrackData& td,
const scalar p,
const scalar T
) const
{
return
this->Y_[0]*td.cloud().composition().cp(0, YGas_, p, T)
+ this->Y_[1]*td.cloud().composition().cp(0, YLiquid_, p, T)
+ this->Y_[2]*td.cloud().composition().cp(0, YSolid_, p, T);
}
template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
(
TrackData& td,
const scalar p,
const scalar T
) const
{
return
this->Y_[0]*td.cloud().composition().H(0, YGas_, p, T)
+ this->Y_[1]*td.cloud().composition().H(0, YLiquid_, p, T)
+ this->Y_[2]*td.cloud().composition().H(0, YSolid_, p, T);
}
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class ParcelType>
......@@ -57,8 +91,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const scalar mass0 = this->mass();
const scalar cp0 = this->cp_;
const scalar np0 = this->nParticle_;
const scalar T0 = this->T_;
scalarField& YMix = this->Y_;
scalar T0 = this->T_;
const scalar pc = this->pc_;
scalarField& Y = this->Y_;
scalar mass1 = mass0;
label idGas = td.cloud().composition().idGas();
label idLiquid = td.cloud().composition().idLiquid();
......@@ -97,16 +133,60 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~
// Calculate phase change
// ~~~~~~~~~~~~~~~~~~~~~~
scalar dMassPC = calcPhaseChange(td, dt, cellI, T1, dMassMT);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate phase change - update mass, Y, cp, T, dhTrans
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar dMassPC = calcPhaseChange(td, dt, cellI, T1, Y[idLiquid], dMassMT);
// Update particle mass
mass1 -= dMassPC;
// Update particle liquid component mass fractions
this->updateMassFraction(mass0, dMassMT, YLiquid_);
// New specific heat capacity of mixture
scalar cp1 = cpEff(td, pc, T1);
// Correct temperature due to evaporated components
// TODO: use hl function in liquidMixture???
// scalar dhPC = -dMassPCTot*td.cloud().composition().L(0, Y_, pc, T0);
scalar Lvap = td.cloud().composition().L(idLiquid, this->YLiquid_, pc, T0);
T1 -= Lvap*dMassPC/(0.5*(mass0 + mass1)*cp1);
// Correct dhTrans to account for the change in enthalpy due to the
// liquid phase change
dhTrans +=
dMassPC
*td.cloud().composition().H(idLiquid, YLiquid_, pc, 0.5*(T0 + T1));
//????????????????????????????????????????????????????????????????????????????????
// Store temperature for the start of the next process
T0 = T1;
//????????????????????????????????????????????????????????????????????????????????
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate Devolatilisation
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
calcDevolatilisation(td, dt, T0, T1, dMassMT);
const scalar dMassD = calcDevolatilisation(td, dt, T0, dMassMT);
// Update particle mass
mass1 -= dMassD;
// New specific heat capacity of mixture
cp1 = cpEff(td, pc, T1);
// Update gas and solid component mass fractions
// updateMassFraction(mass0, mass1, dMassMT, YGas_);
// updateMassFraction(mass0, mass1, dMassMT, YSolid_);
// Correct particle temperature to account for latent heat of
// devolatilisation
T1 -=
td.constProps().Ldevol()
*sum(dMassMT)
/(0.5*(mass0 + mass1)*cp1);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate surface reactions
......@@ -118,14 +198,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
calcSurfaceReactions(td, dt, cellI, T0, T1, dMassMT, dMassSR, dhRet);
// New total mass
const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
mass1 -= sum(dMassSR);
// Correct particle temperature to account for latent heat of
// devolatilisation
T1 -=
td.constProps().Ldevol()
*sum(dMassMT)
/(0.5*(mass0 + mass1)*cp0);
// Enthalpy retention divided between particle and carrier phase by the
// fraction fh and (1 - fh)
......@@ -135,36 +209,39 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Correct dhTrans to account for enthalpy of evolved volatiles
dhTrans +=
sum(dMassMT)
*td.cloud().composition().H(idGas, YGas_, this->pc_, 0.5*(T0 + T1));
*td.cloud().composition().H(idGas, YGas_, pc, 0.5*(T0 + T1));
// Correct dhTrans to account for enthalpy of consumed solids
dhTrans +=
sum(dMassSR)
*td.cloud().composition().H(idSolid, YSolid_, this->pc_, 0.5*(T0 + T1));
*td.cloud().composition().H(idSolid, YSolid_, pc, 0.5*(T0 + T1));
// ~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate source terms
// ~~~~~~~~~~~~~~~~~~~~~~~
// Transfer mass lost from particle to carrier mass source
forAll(dMassMT, i)
if (td.cloud().coupled())
{
td.cloud().rhoTrans(i)[cellI] += np0*(dMassMT[i] + dMassSR[i]);
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Transfer mass lost from particle to carrier mass source
forAll(dMassMT, i)
{
td.cloud().rhoTrans(i)[cellI] += np0*(dMassMT[i] + dMassSR[i]);
}
// Accumulate coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Update enthalpy transfer
// - enthalpy of lost solids already accounted for
td.cloud().hTrans()[cellI] += np0*dhTrans;
// Coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
// Accumulate coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
// Update enthalpy transfer
// - enthalpy of lost solids already accounted for
td.cloud().hTrans()[cellI] += np0*dhTrans;
// Coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
......@@ -179,14 +256,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
{
td.cloud().rhoTrans(i)[cellI] += np0*dMassMT[i];
}
td.cloud().hTrans()[cellI] +=
np0*mass1
*(
YMix[0]
*td.cloud().composition().H(idGas, YGas_, this->pc_, T1)
+ YMix[2]
*td.cloud().composition().H(idSolid, YSolid_, this->pc_, T1)
);
td.cloud().hTrans()[cellI] += np0*mass1*HEff(td, pc, T1);
td.cloud().UTrans()[cellI] += np0*mass1*U1;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
......@@ -196,10 +266,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
{
this->U_ = U1;
this->T_ = T1;
this->cp_ =
YMix[0]*td.cloud().composition().cp(idGas, YGas_, this->pc_, T1)
+ YMix[1]*td.cloud().composition().cp(idLiquid, YLiquid_, this->pc_, T1)
+ YMix[2]*td.cloud().composition().cp(idSolid, YSolid_, this->pc_, T1);
this->cp_ = cpEff(td, pc, T1);
// Update particle density or diameter
if (td.constProps().constantVolume())
......@@ -216,33 +283,14 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
template<class ParcelType>
template<class TrackData>
void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
(
TrackData& td,
const scalar dt,
const scalar T0,
const scalar T1,
const scalar T,
scalarList& dMassMT
)
{
if (td.cloud().composition().YMixture0()[1]>SMALL)
{
notImplemented
(
"void Foam::ReactingMultiphaseParcel<ParcelType>::"
"calcDevolatilisation \n"
"(\n"
" TrackData&,\n"
" const scalar,\n"
" const scalar,\n"
" const scalar,\n"
" scalarList&\n"
")\n"
"no treatment currently available for particles containing "
"liquid species"
);
}
// Check that model is active, and that the parcel temperature is
// within necessary limits for devolatilisation to occur
if
......@@ -252,27 +300,27 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
|| this->T_<td.constProps().Tbp()
)
{
return;
return 0.0;
}
// Determine mass to add to carrier phase
const scalar mass = this->mass();
scalarField& YMix = this->Y_;
scalarField& Y = this->Y_;
const scalar dMassTot = td.cloud().devolatilisation().calculate
(
dt,
this->mass0_,
mass,
td.cloud().composition().YMixture0(),
YMix,
T0,
Y,
T,
canCombust_
);
// Update (total) mass fractions
YMix[0] = (YMix[0]*mass - dMassTot)/(mass - dMassTot);
YMix[1] = YMix[1]*mass/(mass - dMassTot);
YMix[2] = 1.0 - YMix[0] - YMix[1];
Y[0] = (Y[0]*mass - dMassTot)/(mass - dMassTot);
Y[1] = Y[1]*mass/(mass - dMassTot);
Y[2] = 1.0 - Y[0] - Y[1];
// Add to cummulative mass transfer
label idGas = td.cloud().composition().idGas();
......@@ -286,6 +334,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
}
td.cloud().addToMassDevolatilisation(dMassTot);
return dMassTot;
}
......
......@@ -143,6 +143,19 @@ public:
};
private:
// Private member functions
//- Return the mixture effective specific heat capacity
template<class TrackData>
scalar cpEff(TrackData& td, const scalar p, const scalar T) const;
//- Return the mixture effective enthalpy
template<class TrackData>
scalar HEff(TrackData& td, const scalar p, const scalar T) const;
protected:
// Protected data
......@@ -167,12 +180,11 @@ protected:
//- Calculate Devolatilisation
template<class TrackData>
void calcDevolatilisation
scalar calcDevolatilisation
(
TrackData& td,
const scalar dt,
const scalar T0,
const scalar T1,
const scalar T,
scalarList& dMassMT
);
......
......@@ -77,11 +77,11 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
// Check state of Istream
is.check
(
"ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel\n"
"(\n"
" const Cloud<ParcelType>&,\n"
" Istream&,\n"
" bool\n"
"ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel"
"("
"const Cloud<ParcelType>&, "
"Istream&, "
"bool"
")"
);
}
......@@ -102,11 +102,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
// Get names and sizes for each Y...
const label idGas = c.composition().idGas();
const wordList gasNames = c.composition().componentNames(idGas);
const wordList& gasNames = c.composition().componentNames(idGas);
const label idLiquid = c.composition().idLiquid();
const wordList liquidNames = c.composition().componentNames(idLiquid);
const wordList& liquidNames = c.composition().componentNames(idLiquid);
const label idSolid = c.composition().idSolid();
const wordList solidNames = c.composition().componentNames(idSolid);
const wordList& solidNames = c.composition().componentNames(idSolid);
const wordList& stateLabels = c.composition().stateLabels();
// Set storage for each Y... for each parcel
forAllIter(typename Cloud<ParcelType>, c, iter)
......@@ -122,7 +123,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
{
IOField<scalar> YGas
(
c.fieldIOobject("Y" + gasNames[j], IOobject::MUST_READ)
c.fieldIOobject
(
"Y" + gasNames[j] + stateLabels[idGas],
IOobject::MUST_READ
)
);
label i = 0;
......@@ -137,7 +142,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
{
IOField<scalar> YLiquid
(
c.fieldIOobject("Y" + liquidNames[j], IOobject::MUST_READ)
c.fieldIOobject
(
"Y" + liquidNames[j] + stateLabels[idLiquid],
IOobject::MUST_READ
)
);
label i = 0;
......@@ -152,7 +161,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
{
IOField<scalar> YSolid
(
c.fieldIOobject("Y" + solidNames[j], IOobject::MUST_READ)
c.fieldIOobject
(
"Y" + solidNames[j] + stateLabels[idSolid],
IOobject::MUST_READ
)
);
label i = 0;
......@@ -178,13 +191,19 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
// Write the composition fractions
if (np > 0)
{
const wordList& stateLabels = c.composition().stateLabels();
const label idGas = c.composition().idGas();
const wordList gasNames = c.composition().componentNames(idGas);
const wordList& gasNames = c.composition().componentNames(idGas);
forAll(gasNames, j)
{
IOField<scalar> YGas