updated coupling with carrier phase

This commit is contained in:
andy
2009-05-18 11:08:07 +01:00
parent 34254a29b0
commit f673d3e338
14 changed files with 333 additions and 81 deletions

View File

@ -96,6 +96,10 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
constProps_(particleProperties_),
parcelTypeId_(readLabel(particleProperties_.lookup("parcelTypeId"))),
coupled_(particleProperties_.lookup("coupled")),
cellValueSourceCorrection_
(
particleProperties_.lookup("cellValueSourceCorrection")
),
rndGen_(label(0)),
rho_(rho),
U_(U),

View File

@ -105,7 +105,11 @@ class KinematicCloud
//- Flag to indicate whether parcels are coupled to the carrier phase
// i.e. whether or not to generate source terms for carrier phase
Switch coupled_;
const Switch coupled_;
//- Flag to correct cell values with latest transfer information
// during the lagrangian timestep
const Switch cellValueSourceCorrection_;
//- Random number generator - used by some injection routines
Random rndGen_;
@ -224,6 +228,9 @@ public:
//- Return coupled flag
inline const Switch coupled() const;
//- Return cell value correction flag
inline const Switch cellValueSourceCorrection() const;
//- Return refernce to the random object
inline Random& rndGen();

View File

@ -72,6 +72,14 @@ inline const Foam::Switch Foam::KinematicCloud<ParcelType>::coupled() const
}
template <class ParcelType>
inline const Foam::Switch
Foam::KinematicCloud<ParcelType>::cellValueSourceCorrection() const
{
return cellValueSourceCorrection_;
}
template<class ParcelType>
inline const Foam::volScalarField&
Foam::KinematicCloud<ParcelType>::rho() const

View File

@ -102,11 +102,25 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
)
),
radiation_(this->particleProperties().lookup("radiation")),
hTrans_
hsTrans_
(
IOobject
(
this->name() + "hTrans",
this->name() + "hsTrans",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh(),
dimensionedScalar("zero", dimensionSet(1, 2, -2, 0, 0), 0.0)
),
hcTrans_
(
IOobject
(
this->name() + "hcTrans",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
@ -132,7 +146,8 @@ template<class ParcelType>
void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
{
KinematicCloud<ParcelType>::resetSourceTerms();
hTrans_.field() = 0.0;
hsTrans_.field() = 0.0;
hcTrans_.field() = 0.0;
}

View File

@ -99,8 +99,14 @@ class ThermoCloud
// Sources
//- Enthalpy transfer
DimensionedField<scalar, volMesh> hTrans_;
//- Sensible enthalpy transfer
DimensionedField<scalar, volMesh> hsTrans_;
//- Chemical enthalpy transfer
// - If solving for total enthalpy, the carrier phase enthalpy will
// receive the full enthalpy of reaction via creation of reaction
// products
DimensionedField<scalar, volMesh> hcTrans_;
// Private Member Functions
@ -173,10 +179,19 @@ public:
// Enthalpy
//- Return reference to enthalpy source
inline DimensionedField<scalar, volMesh>& hTrans();
//- Return reference to sensible enthalpy source
inline DimensionedField<scalar, volMesh>& hsTrans();
//- return tmp enthalpy source term - fully explicit
//- Return tmp total sensible enthalpy source term
inline tmp<DimensionedField<scalar, volMesh> > Shs() const;
//- Return reference to chemical enthalpy source
inline DimensionedField<scalar, volMesh>& hcTrans();
//- Return tmp chemical enthalpy source term
inline tmp<DimensionedField<scalar, volMesh> > Shc() const;
//- Return tmp total enthalpy source term
inline tmp<DimensionedField<scalar, volMesh> > Sh() const;

View File

@ -77,9 +77,85 @@ inline bool Foam::ThermoCloud<ParcelType>::radiation() const
template<class ParcelType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<ParcelType>::hTrans()
Foam::ThermoCloud<ParcelType>::hsTrans()
{
return hTrans_;
return hsTrans_;
}
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::ThermoCloud<ParcelType>::Shs() const
{
tmp<DimensionedField<scalar, volMesh> > tShs
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "Shs",
this->db().time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
this->mesh(),
dimensionedScalar
(
"zero",
dimMass/dimLength/pow3(dimTime),
0.0
)
)
);
scalarField& Shs = tShs().field();
Shs = hsTrans_/(this->mesh().V()*this->db().time().deltaT());
return tShs;
}
template<class ParcelType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<ParcelType>::hcTrans()
{
return hcTrans_;
}
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::ThermoCloud<ParcelType>::Shc() const
{
tmp<DimensionedField<scalar, volMesh> > tShc
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "Shc",
this->db().time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
this->mesh(),
dimensionedScalar
(
"zero",
dimMass/dimLength/pow3(dimTime),
0.0
)
)
);
scalarField& Shc = tShc().field();
Shc = hcTrans_/(this->mesh().V()*this->db().time().deltaT());
return tShc;
}
@ -111,7 +187,7 @@ Foam::ThermoCloud<ParcelType>::Sh() const
);
scalarField& Sh = tSh().field();
Sh = hTrans_/(this->mesh().V()*this->db().time().deltaT());
Sh = (hsTrans_ + hcTrans_)/(this->mesh().V()*this->db().time().deltaT());
return tSh;
}

View File

@ -31,7 +31,7 @@ License
template<class ParcelType>
template<class TrackData>
void Foam::KinematicParcel<ParcelType>::updateCellQuantities
void Foam::KinematicParcel<ParcelType>::setCellValues
(
TrackData& td,
const scalar dt,
@ -43,7 +43,7 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
{
WarningIn
(
"void Foam::KinematicParcel<ParcelType>::updateCellQuantities"
"void Foam::KinematicParcel<ParcelType>::setCellValues"
"("
"TrackData&, "
"const scalar, "
@ -57,9 +57,6 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
Uc_ = td.UInterp().interpolate(this->position(), cellI);
// Apply correction to cell velocity to account for momentum transfer
Uc_ += td.cloud().UTrans()[cellI]/(massCell(cellI));
muc_ = td.muInterp().interpolate(this->position(), cellI);
// Apply dispersion components to carrier phase velocity
@ -75,6 +72,19 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
}
template<class ParcelType>
template<class TrackData>
void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection
(
TrackData& td,
const scalar dt,
const label cellI
)
{
Uc_ += td.cloud().UTrans()[cellI]/massCell(cellI);
}
template<class ParcelType>
template<class TrackData>
void Foam::KinematicParcel<ParcelType>::calc
@ -183,8 +193,8 @@ bool Foam::KinematicParcel<ParcelType>::move(TrackData& td)
// Set the Lagrangian time-step
scalar dt = min(dtMax, tEnd);
// Remember which cell the Parcel is in
// since this will change if a face is hit
// Remember which cell the Parcel is in since this will change if a
// face is hit
label cellI = p.cell();
dt *= p.trackToFace(p.position() + dt*U_, td);
@ -192,12 +202,17 @@ bool Foam::KinematicParcel<ParcelType>::move(TrackData& td)
tEnd -= dt;
p.stepFraction() = 1.0 - tEnd/deltaT;
// Update cell based properties
p.updateCellQuantities(td, dt, cellI);
// Avoid problems with extremely small timesteps
if (dt > ROOTVSMALL)
{
// Update cell based properties
p.setCellValues(td, dt, cellI);
if (td.cloud().cellValueSourceCorrection())
{
p.cellValueSourceCorrection(td, dt, cellI);
}
p.calc(td, dt, cellI);
}

View File

@ -372,9 +372,18 @@ public:
// Main calculation loop
//- Update cell based quantities
//- Set cell values
template<class TrackData>
void updateCellQuantities
void setCellValues
(
TrackData& td,
const scalar dt,
const label cellI
);
//- Correct cell values using latest transfer information
template<class TrackData>
void cellValueSourceCorrection
(
TrackData& td,
const scalar dt,

View File

@ -129,14 +129,54 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
template<class ParcelType>
template<class TrackData>
void Foam::ReactingMultiphaseParcel<ParcelType>::updateCellQuantities
void Foam::ReactingMultiphaseParcel<ParcelType>::setCellValues
(
TrackData& td,
const scalar dt,
const label cellI
)
{
ReactingParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
ReactingParcel<ParcelType>::setCellValues(td, dt, cellI);
}
template<class ParcelType>
template<class TrackData>
void Foam::ReactingMultiphaseParcel<ParcelType>::cellValueSourceCorrection
(
TrackData& td,
const scalar dt,
const label cellI
)
{
scalar massCell = this->massCell(cellI);
scalar addedMass = 0.0;
forAll(td.cloud().rhoTrans(), i)
{
addedMass += td.cloud().rhoTrans(i)[cellI];
}
this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
scalar massCellNew = massCell + addedMass;
this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
scalar cpEff = 0;
forAll(td.cloud().rhoTrans(), i)
{
scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
cpEff += Y*td.cloud().gases()[i].Cp(this->Tc_);
}
const scalar cpc = td.cpInterp().psi()[cellI];
this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
const scalar fCarrier = -1.0/td.constProps().hRetentionCoeff();
const scalar dh =
td.cloud().hsTrans()[cellI] + fCarrier*td.cloud().hcTrans[cellI];
this->Tc_ += dh/(this->cpc_*massCellNew);
}
@ -229,7 +269,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
scalarField dMassSRCarrier(td.cloud().gases().size(), 0.0);
// Return enthalpy source and calc mass transfer(s) due to surface reaction
scalar HReaction =
scalar hReaction =
calcSurfaceReactions
(
td,
@ -250,11 +290,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
dMassSRCarrier
);
// Heat of reaction split between component retained by particle
const scalar ShSR = td.constProps().hRetentionCoeff()*HReaction;
// ...and component added to the carrier phase
const scalar ShSRc = (1.0 - td.constProps().hRetentionCoeff())*HReaction;
// Heat of reaction retained by particle
const scalar ShSR = td.constProps().hRetentionCoeff()*hReaction;
// Update component mass fractions
@ -310,12 +347,14 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
label id = td.cloud().composition().localToGlobalGasId(LIQ, i);
td.cloud().rhoTrans(id)[cellI] += np0*dMassLiquid[i];
}
// // No mapping between solid components and carrier phase
// forAll(YSolid_, i)
// {
// label id = td.cloud().composition().localToGlobalGasId(SLD, i);
// td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
// }
/*
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
label id = td.cloud().composition().localToGlobalGasId(SLD, i);
td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
}
*/
forAll(dMassSRCarrier, i)
{
td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
@ -324,8 +363,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
// Update enthalpy transfer
td.cloud().hTrans()[cellI] += np0*(mass0*H0 - (mass1*H1 + ShSRc));
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
// Update chemical enthalpy transfer
td.cloud().hcTrans()[cellI] -= np0*ShSR;
}
@ -346,21 +388,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
}
forAll(YLiquid_, i)
{
label id =
td.cloud().composition().localToGlobalGasId(LIQ, i);
label id = td.cloud().composition().localToGlobalGasId(LIQ, i);
td.cloud().rhoTrans(id)[cellI] +=
np0*mass1*YMix[LIQ]*YLiquid_[i];
}
// // No mapping between solid components and carrier phase
// forAll(YSolid_, i)
// {
// label id =
// td.cloud().composition().localToGlobalGasId(SLD, i);
// td.cloud().rhoTrans(id)[cellI] +=
// np0*mass1*YMix[SLD]*YSolid_[i];
// }
td.cloud().hTrans()[cellI] += np0*mass1*H1;
/*
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
label id = td.cloud().composition().localToGlobalGasId(SLD, i);
td.cloud().rhoTrans(id)[cellI] +=
np0*mass1*YMix[SLD]*YSolid_[i];
}
*/
td.cloud().hsTrans()[cellI] += np0*mass1*H1;
td.cloud().UTrans()[cellI] += np0*mass1*U1;
}
}
@ -466,7 +507,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
}
// Update surface reactions
const scalar HReaction = td.cloud().surfaceReaction().calculate
const scalar hReaction = td.cloud().surfaceReaction().calculate
(
dt,
cellI,
@ -493,7 +534,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
*(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))
);
return HReaction;
return hReaction;
}
@ -502,4 +543,3 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
#include "ReactingMultiphaseParcelIO.C"
// ************************************************************************* //

View File

@ -333,9 +333,18 @@ public:
// Main calculation loop
//- Update cell based quantities
//- Set cell values
template<class TrackData>
void updateCellQuantities
void setCellValues
(
TrackData& td,
const scalar dt,
const label cellI
);
//- Correct cell values using latest transfer information
template<class TrackData>
void cellValueSourceCorrection
(
TrackData& td,
const scalar dt,
@ -393,4 +402,3 @@ public:
#endif
// ************************************************************************* //

View File

@ -30,21 +30,21 @@ License
template<class ParcelType>
template<class TrackData>
void Foam::ReactingParcel<ParcelType>::updateCellQuantities
void Foam::ReactingParcel<ParcelType>::setCellValues
(
TrackData& td,
const scalar dt,
const label cellI
)
{
ThermoParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
ThermoParcel<ParcelType>::setCellValues(td, dt, cellI);
pc_ = td.pInterp().interpolate(this->position(), cellI);
if (pc_ < td.constProps().pMin())
{
WarningIn
(
"void Foam::ReactingParcel<ParcelType>::updateCellQuantities"
"void Foam::ReactingParcel<ParcelType>::setCellValues"
"("
"TrackData&, "
"const scalar, "
@ -55,14 +55,42 @@ void Foam::ReactingParcel<ParcelType>::updateCellQuantities
pc_ = td.constProps().pMin();
}
}
template<class ParcelType>
template<class TrackData>
void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
(
TrackData& td,
const scalar dt,
const label cellI
)
{
scalar massCell = this->massCell(cellI);
// Apply correction to cell density to account for mass transfer
scalar addedMass = 0.0;
forAll(td.cloud().rhoTrans(), i)
{
addedMass += td.cloud().rhoTrans(i)[cellI];
}
this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
scalar massCellNew = massCell + addedMass;
this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
scalar cpEff = 0;
forAll(td.cloud().rhoTrans(), i)
{
scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
cpEff += Y*td.cloud().gases()[i].Cp(this->Tc_);
}
const scalar cpc = td.cpInterp().psi()[cellI];
this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
this->Tc_ += td.cloud().hsTrans()[cellI]/(this->cpc_*massCellNew);
}
@ -164,8 +192,8 @@ void Foam::ReactingParcel<ParcelType>::calc
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
// Update enthalpy transfer
td.cloud().hTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
}
@ -184,7 +212,7 @@ void Foam::ReactingParcel<ParcelType>::calc
td.cloud().rhoTrans(id)[cellI] += np0*mass1*Y_[i];
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
td.cloud().hTrans()[cellI] += np0*mass1*H1;
td.cloud().hsTrans()[cellI] += np0*mass1*H1;
}
}

View File

@ -276,9 +276,18 @@ public:
// Main calculation loop
//- Update cell based quantities
//- Set cell values
template<class TrackData>
void updateCellQuantities
void setCellValues
(
TrackData& td,
const scalar dt,
const label cellI
);
//- Correct cell values using latest transfer information
template<class TrackData>
void cellValueSourceCorrection
(
TrackData& td,
const scalar dt,
@ -333,4 +342,3 @@ public:
#endif
// ************************************************************************* //

View File

@ -30,28 +30,24 @@ License
template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::updateCellQuantities
void Foam::ThermoParcel<ParcelType>::setCellValues
(
TrackData& td,
const scalar dt,
const label cellI
)
{
KinematicParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
KinematicParcel<ParcelType>::setCellValues(td, dt, cellI);
cpc_ = td.cpInterp().interpolate(this->position(), cellI);
Tc_ = td.TInterp().interpolate(this->position(), cellI);
// Apply correction to cell temperature to account for enthalpy transfer
scalar cpMean = td.cpInterp().psi()[cellI];
Tc_ += td.cloud().hTrans()[cellI]/(cpMean*this->massCell(cellI));
if (Tc_ < td.constProps().TMin())
{
WarningIn
(
"void Foam::ThermoParcel<ParcelType>::updateCellQuantities"
"void Foam::ThermoParcel<ParcelType>::setCellValues"
"("
"TrackData&, "
"const scalar, "
@ -65,6 +61,22 @@ void Foam::ThermoParcel<ParcelType>::updateCellQuantities
}
template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
(
TrackData& td,
const scalar dt,
const label cellI
)
{
this->Uc_ += td.cloud().UTrans()[cellI]/this->massCell(cellI);
scalar cpMean = td.cpInterp().psi()[cellI];
Tc_ += td.cloud().hsTrans()[cellI]/(cpMean*this->massCell(cellI));
}
template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::calc
@ -118,8 +130,8 @@ void Foam::ThermoParcel<ParcelType>::calc
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*mass0*(U0 - U1);
// Update enthalpy transfer
td.cloud().hTrans()[cellI] += np0*mass0*(H0 - H1);
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*mass0*(H0 - H1);
}
// Set new particle properties
@ -198,4 +210,3 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
#include "ThermoParcelIO.C"
// ************************************************************************* //

View File

@ -286,9 +286,18 @@ public:
// Main calculation loop
//- Update cell based quantities
//- Set cell values
template<class TrackData>
void updateCellQuantities
void setCellValues
(
TrackData& td,
const scalar dt,
const label cellI
);
//- Correct cell values using latest transfer information
template<class TrackData>
void cellValueSourceCorrection
(
TrackData& td,
const scalar dt,
@ -343,4 +352,3 @@ public:
#endif
// ************************************************************************* //