zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

updated coupling with carrier phase

Browse Source
This commit is contained in:
andy
2009-05-14 13:26:41 +01:00
parent d89df4142e
commit 30c7e8aa6c
16 changed files with 170 additions and 404 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
View File

@ -155,21 +155,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
false
),
mesh_,
dimensionedVector("zero", dimensionSet(1, 1, -1, 0, 0), vector::zero)
),
UCoeff_
(
IOobject
(
this->name() + "UCoeff",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimensionedScalar("zero", dimensionSet(1, 0, -1, 0, 0), 0.0)
dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
)
{}
@ -187,7 +173,6 @@ template<class ParcelType>
void Foam::KinematicCloud<ParcelType>::resetSourceTerms()
{
UTrans_.field() = vector::zero;
UCoeff_.field() = 0.0;
}

11
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
View File

@ -165,9 +165,6 @@ class KinematicCloud
//- Momentum
DimensionedField<vector, volMesh> UTrans_;
//- Coefficient for carrier phase U equation
DimensionedField<scalar, volMesh> UCoeff_;
// Private Member Functions
@ -309,14 +306,8 @@ public:
//- Return reference to momentum source
inline DimensionedField<vector, volMesh>& UTrans();
//- Coefficient for carrier phase U equation
inline DimensionedField<scalar, volMesh>& UCoeff();
//- Return tmp momentum source term - fully explicit
inline tmp<DimensionedField<vector, volMesh> > SU1() const;
//- Return tmp momentum source term - semi-implicit
inline tmp<fvVectorMatrix> SU2(volVectorField& U) const;
inline tmp<DimensionedField<vector, volMesh> > SU() const;
// Check

39
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
View File

@ -203,25 +203,17 @@ Foam::KinematicCloud<ParcelType>::UTrans()
}
template<class ParcelType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::KinematicCloud<ParcelType>::UCoeff()
{
return UCoeff_;
}
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
Foam::KinematicCloud<ParcelType>::SU1() const
Foam::KinematicCloud<ParcelType>::SU() const
{
tmp<DimensionedField<vector, volMesh> > tSU1
tmp<DimensionedField<vector, volMesh> > tSU
(
new DimensionedField<vector, volMesh>
(
IOobject
(
this->name() + "SU1",
this->name() + "SU",
this->db().time().timeName(),
this->mesh(),
IOobject::NO_READ,
@ -237,29 +229,10 @@ Foam::KinematicCloud<ParcelType>::SU1() const
)
);
vectorField& SU1 = tSU1().field();
SU1 = UTrans_/(mesh_.V()*this->db().time().deltaT());
vectorField& SU = tSU().field();
SU = UTrans_/(mesh_.V()*this->db().time().deltaT());
return tSU1;
}
template<class ParcelType>
inline Foam::tmp<Foam::fvVectorMatrix>
Foam::KinematicCloud<ParcelType>::SU2(volVectorField& U) const
{
if (debug)
{
Info<< "UTrans min/max = "
<< min(UTrans_) << ", " << max(UTrans_) << endl;
Info<< "UCoeff min/max = "
<< min(UCoeff_) << ", " << max(UCoeff_) << endl;
}
return
UTrans_/(mesh_.V()*this->db().time().deltaT())
- fvm::Sp(UCoeff_/mesh_.V(), U)
+ UCoeff_/mesh_.V()*U;
return tSU;
}

10
src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
View File

@ -180,15 +180,13 @@ public:
inline PtrList<DimensionedField<scalar, volMesh> >&
rhoTrans();
//- Return tmp mass source for field i
// Fully explicit
//- Return tmp mass source for field i - fully explicit
inline tmp<DimensionedField<scalar, volMesh> >
Srho1(const label i) const;
Srho(const label i) const;
//- Return tmp total mass source for carrier phase
// Fully explicit
inline tmp<DimensionedField<scalar, volMesh> >
Srho1() const;
// - fully explicit
inline tmp<DimensionedField<scalar, volMesh> > Srho() const;
// Check

4
src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
View File

@ -92,7 +92,7 @@ Foam::ReactingCloud<ParcelType>::rhoTrans()
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::ReactingCloud<ParcelType>::Srho1(const label i) const
Foam::ReactingCloud<ParcelType>::Srho(const label i) const
{
return rhoTrans_[i]/(this->db().time().deltaT()*this->mesh().V());
}
@ -100,7 +100,7 @@ Foam::ReactingCloud<ParcelType>::Srho1(const label i) const
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::ReactingCloud<ParcelType>::Srho1() const
Foam::ReactingCloud<ParcelType>::Srho() const
{
tmp<DimensionedField<scalar, volMesh> > trhoTrans
(

15
src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
View File

@ -115,20 +115,6 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
),
this->mesh(),
dimensionedScalar("zero", dimensionSet(1, 2, -2, 0, 0), 0.0)
),
hCoeff_
(
IOobject
(
this->name() + "hCoeff",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh(),
dimensionedScalar("zero", dimensionSet(1, 2, -3, -1, 0), 0.0)
)
{}
@ -147,7 +133,6 @@ void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
{
KinematicCloud<ParcelType>::resetSourceTerms();
hTrans_.field() = 0.0;
hCoeff_.field() = 0.0;
}

11
src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
View File

@ -102,9 +102,6 @@ class ThermoCloud
//- Enthalpy transfer
DimensionedField<scalar, volMesh> hTrans_;
//- Coefficient for carrier phase h equation
DimensionedField<scalar, volMesh> hCoeff_;
// Private Member Functions
@ -179,14 +176,8 @@ public:
//- Return reference to enthalpy source
inline DimensionedField<scalar, volMesh>& hTrans();
//- Coefficient for carrier phase h equation
inline DimensionedField<scalar, volMesh>& hCoeff();
//- return tmp enthalpy source term - fully explicit
inline tmp<DimensionedField<scalar, volMesh> > Sh1() const;
//- Return tmp enthalpy source term - semi-implicit
inline tmp<fvScalarMatrix> Sh2(volScalarField& h) const;
inline tmp<DimensionedField<scalar, volMesh> > Sh() const;
// Radiation - overrides thermoCloud virtual abstract members

40
src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
View File

@ -83,25 +83,17 @@ Foam::ThermoCloud<ParcelType>::hTrans()
}
template<class ParcelType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<ParcelType>::hCoeff()
{
return hCoeff_;
}
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::ThermoCloud<ParcelType>::Sh1() const
Foam::ThermoCloud<ParcelType>::Sh() const
{
tmp<DimensionedField<scalar, volMesh> > tSh1
tmp<DimensionedField<scalar, volMesh> > tSh
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "Sh1",
this->name() + "Sh",
this->db().time().timeName(),
this->mesh(),
IOobject::NO_READ,
@ -118,30 +110,10 @@ Foam::ThermoCloud<ParcelType>::Sh1() const
)
);
scalarField& Sh1 = tSh1().field();
Sh1 = hTrans_/(this->mesh().V()*this->db().time().deltaT());
scalarField& Sh = tSh().field();
Sh = hTrans_/(this->mesh().V()*this->db().time().deltaT());
return tSh1;
}
template<class ParcelType>
inline Foam::tmp<Foam::fvScalarMatrix>
Foam::ThermoCloud<ParcelType>::Sh2(volScalarField& h) const
{
const volScalarField cp = carrierThermo_.Cp();
if (debug)
{
Info<< "hTrans min/max = "
<< min(hTrans_) << ", " << max(hTrans_) << endl;
Info<< "hCoeff min/max = "
<< min(hCoeff_) << ", " << max(hCoeff_) << endl;
}
return hTrans_/(this->mesh().V()*this->db().time().deltaT())
- fvm::Sp(hCoeff_/(cp*this->mesh().V()), h)
+ hCoeff_/(cp*this->mesh().V())*h;
return tSh;
}

28
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
View File

@ -92,15 +92,6 @@ void Foam::KinematicParcel<ParcelType>::calc
const scalar rho0 = rho_;
const scalar mass0 = mass();
const polyMesh& mesh = this->cloud().pMesh();
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum
vector dUTrans = vector::zero;
// Motion
// ~~~~~~
@ -109,9 +100,7 @@ void Foam::KinematicParcel<ParcelType>::calc
vector Fx = vector::zero;
// Calculate new particle velocity
scalar Cud = 0.0;
vector U1 =
calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx, Cud, dUTrans);
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
// Accumulate carrier phase source terms
@ -119,10 +108,7 @@ void Foam::KinematicParcel<ParcelType>::calc
if (td.cloud().coupled())
{
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
td.cloud().UTrans()[cellI] += np0*mass0*(U0 - U1);
}
@ -143,15 +129,13 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
const vector& U,
const scalar rho,
const scalar mass,
const vector& Fx,
scalar& Cud,
vector& dUTrans
const vector& Fx
) const
{
const polyMesh& mesh = this->cloud().pMesh();
// Return linearised term from drag model
Cud = td.cloud().drag().Cu(U - Uc_, d, rhoc_, rho, muc_);
scalar Cud = td.cloud().drag().Cu(U - Uc_, d, rhoc_, rho, muc_);
// Calculate particle forces
vector Ftot = td.cloud().forces().calc(cellI, dt, rhoc_, rho, Uc_, U);
@ -169,10 +153,6 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
// Apply correction to velocity for reduced-D cases
meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
// Calculate the momentum transfer to the continuous phase
// - do not include gravity impulse
dUTrans = -mass*(Unew - U - dt*td.g());
return Unew;
}

4
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
View File

@ -243,9 +243,7 @@ protected:
const vector& U, // velocity
const scalar rho, // density
const scalar mass, // mass
const vector& Fx, // additional forces
scalar& Cud, // linearised drag term coeff
vector& dUTrans // momentum transfer to carrier phase
const vector& Fx // additional forces
) const;

218
src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
View File

@ -32,10 +32,10 @@ template<class ParcelType>
const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::GAS(0);
template<class ParcelType>
const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::LIQUID(1);
const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::LIQ(1);
template<class ParcelType>
const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::SOLID(2);
const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::SLD(2);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@ -54,8 +54,8 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::cpEff
{
return
this->Y_[GAS]*td.cloud().composition().cp(idG, YGas_, p, T)
+ this->Y_[LIQUID]*td.cloud().composition().cp(idL, YLiquid_, p, T)
+ this->Y_[SOLID]*td.cloud().composition().cp(idS, YSolid_, p, T);
+ this->Y_[LIQ]*td.cloud().composition().cp(idL, YLiquid_, p, T)
+ this->Y_[SLD]*td.cloud().composition().cp(idS, YSolid_, p, T);
}
@ -73,8 +73,27 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
{
return
this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
+ this->Y_[LIQUID]*td.cloud().composition().H(idL, YLiquid_, p, T)
+ this->Y_[SOLID]*td.cloud().composition().H(idS, YSolid_, p, T);
+ this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
+ this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
}
template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::LEff
(
TrackData& td,
const scalar p,
const scalar T,
const label idG,
const label idL,
const label idS
) const
{
return
this->Y_[GAS]*td.cloud().composition().L(idG, YGas_, p, T)
+ this->Y_[LIQ]*td.cloud().composition().L(idL, YLiquid_, p, T)
+ this->Y_[SLD]*td.cloud().composition().L(idS, YSolid_, p, T);
}
@ -94,14 +113,14 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
scalar dMassSolidTot = sum(dMassSolid);
this->updateMassFraction(mass0*YMix[GAS], dMassGas, YGas_);
this->updateMassFraction(mass0*YMix[LIQUID], dMassLiquid, YLiquid_);
this->updateMassFraction(mass0*YMix[SOLID], dMassSolid, YSolid_);
this->updateMassFraction(mass0*YMix[LIQ], dMassLiquid, YLiquid_);
this->updateMassFraction(mass0*YMix[SLD], dMassSolid, YSolid_);
scalar massNew = mass0 - (dMassGasTot + dMassLiquidTot + dMassSolidTot);
YMix[GAS] = (mass0*YMix[GAS] - dMassGasTot)/massNew;
YMix[LIQUID] = (mass0*YMix[LIQUID] - dMassLiquidTot)/massNew;
YMix[SOLID] = 1.0 - YMix[GAS] - YMix[LIQUID];
YMix[LIQ] = (mass0*YMix[LIQ] - dMassLiquidTot)/massNew;
YMix[SLD] = 1.0 - YMix[GAS] - YMix[LIQ];
return massNew;
}
@ -148,32 +167,18 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const label idL = td.cloud().composition().idLiquid();
const label idS = td.cloud().composition().idSolid();
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum
vector dUTrans = vector::zero;
// Enthalpy
scalar dhTrans = 0.0;
// Mass transfer due to phase change
scalarField dMassPC(YLiquid_.size(), 0.0);
// Mass transfer due to devolatilisation
scalarField dMassDV(YGas_.size(), 0.0);
// Change in carrier phase composition due to surface reactions
scalarField dMassSRGas(YGas_.size(), 0.0);
scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
scalarField dMassSRSolid(YSolid_.size(), 0.0);
scalarField dMassSRCarrier(td.cloud().gases().size(), 0.0);
// Intial ethalpy state
scalar H0H = HEff(td, pc, T0, idG, idL, idS);
scalar H0L = LEff(td, pc, T0, idG, idL, idS);
scalar H0 = H0H - H0L;
// Phase change in liquid phase
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Mass transfer due to phase change
scalarField dMassPC(YLiquid_.size(), 0.0);
// Return enthalpy source and calc mass transfer due to phase change
scalar ShPC =
calcPhaseChange
@ -185,7 +190,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
T0,
U0,
idL,
YMix[LIQUID],
YMix[LIQ],
YLiquid_,
dMassPC
);
@ -194,6 +199,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Devolatilisation
// ~~~~~~~~~~~~~~~~
// Mass transfer due to devolatilisation
scalarField dMassDV(YGas_.size(), 0.0);
// Return enthalpy source and calc mass transfer due to devolatilisation
scalar ShDV =
calcDevolatilisation
@ -214,8 +222,14 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Surface reactions
// ~~~~~~~~~~~~~~~~~
// Change in carrier phase composition due to surface reactions
scalarField dMassSRGas(YGas_.size(), 0.0);
scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
scalarField dMassSRSolid(YSolid_.size(), 0.0);
scalarField dMassSRCarrier(td.cloud().gases().size(), 0.0);
// Return enthalpy source and calc mass transfer(s) due to surface reaction
scalar ShSR =
scalar HReaction =
calcSurfaceReactions
(
td,
@ -233,34 +247,14 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
dMassSRGas,
dMassSRLiquid,
dMassSRSolid,
dMassSRCarrier,
dhTrans
dMassSRCarrier
);
// Heat of reaction split between component retained by particle
const scalar ShSR = td.constProps().hRetentionCoeff()*HReaction;
// Heat transfer
// ~~~~~~~~~~~~~
// Total enthalpy source
scalar Sh = ShPC + ShDV + ShSR;
// Calculate new particle temperature
scalar htc = 0.0;
scalar T1 =
calcHeatTransfer
(
td,
dt,
cellI,
d0,
U0,
rho0,
T0,
cp0,
Sh,
htc,
dhTrans
);
// ...and component added to the carrier phase
const scalar ShSRc = (1.0 - td.constProps().hRetentionCoeff())*HReaction;
// Update component mass fractions
@ -274,6 +268,22 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
// Heat transfer
// ~~~~~~~~~~~~~
// Total enthalpy source
scalar Sh = ShPC + ShDV + ShSR;
// Calculate new particle temperature
scalar T1 = calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh);
// Calculate new enthalpy state based on updated composition at new
// temperature
scalar H1H = HEff(td, pc, T1, idG, idL, idS);
scalar H1L = LEff(td, pc, T1, idG, idL, idS);
scalar H1 = H1H - H1L;
// Motion
// ~~~~~~
@ -281,21 +291,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
vector Fx = vector::zero;
// Calculate new particle velocity
scalar Cud = 0.0;
vector U1 =
calcVelocity
(
td,
dt,
cellI,
d0,
U0,
rho0,
0.5*(mass0 + mass1),
Fx,
Cud,
dUTrans
);
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
// Accumulate carrier phase source terms
@ -311,13 +307,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
}
forAll(YLiquid_, i)
{
label id = td.cloud().composition().localToGlobalGasId(LIQUID, i);
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(SOLID, i);
// label id = td.cloud().composition().localToGlobalGasId(SLD, i);
// td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
// }
forAll(dMassSRCarrier, i)
@ -326,16 +322,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
}
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
// Update enthalpy transfer
td.cloud().hTrans()[cellI] += np0*dhTrans;
// Coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
td.cloud().hTrans()[cellI] += np0*(mass0*H0 - (mass1*H1 + ShSRc));
}
@ -357,29 +347,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
forAll(YLiquid_, i)
{
label id =
td.cloud().composition().localToGlobalGasId(LIQUID, i);
td.cloud().composition().localToGlobalGasId(LIQ, i);
td.cloud().rhoTrans(id)[cellI] +=
np0
*mass1
*YMix[LIQUID]
*YLiquid_[i];
np0*mass1*YMix[LIQ]*YLiquid_[i];
}
// // No mapping between solid components and carrier phase
// forAll(YSolid_, i)
// {
// label id =
// td.cloud().composition().localToGlobalGasId(SOLID, i);
// td.cloud().composition().localToGlobalGasId(SLD, i);
// td.cloud().rhoTrans(id)[cellI] +=
// np0
// *mass1
// *YMix[SOLID]
// *YSolid_[i];
// np0*mass1*YMix[SLD]*YSolid_[i];
// }
td.cloud().hTrans()[cellI] +=
np0
*mass1
*HEff(td, pc, T1, idG, idL, idS);
td.cloud().hTrans()[cellI] += np0*mass1*H1;
td.cloud().UTrans()[cellI] += np0*mass1*U1;
}
}
@ -421,7 +402,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
const scalarField& YVolatile,
bool& canCombust,
scalarField& dMassDV
)
) const
{
// Check that model is active, and that the parcel temperature is
// within necessary limits for devolatilisation to occur
@ -448,10 +429,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
);
// Volatile mass transfer - equal components of each volatile specie
forAll(YVolatile, i)
{
dMassDV[i] = YVolatile[i]*dMassTot;
}
dMassDV = YVolatile*dMassTot;
td.cloud().addToMassDevolatilisation(this->nParticle_*dMassTot);
@ -478,9 +456,8 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
scalarField& dMassSRGas,
scalarField& dMassSRLiquid,
scalarField& dMassSRSolid,
scalarField& dMassSRCarrier,
scalar& dhTrans
)
scalarField& dMassSRCarrier
) const
{
// Check that model is active
if (!td.cloud().surfaceReaction().active() || !canCombust)
@ -489,18 +466,15 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
}
// Update surface reactions
const scalar pc = this->pc_;
const scalar Tc = this->Tc_;
const scalar rhoc = this->rhoc_;
const scalar HReaction = td.cloud().surfaceReaction().calculate
(
dt,
cellI,
d,
T,
Tc,
pc,
rhoc,
this->Tc_,
this->pc_,
this->rhoc_,
mass,
YGas,
YLiquid,
@ -519,27 +493,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
*(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))
);
// Add enthalpy of consumed components to the carrier phase enthalpy
// transfer
const label idG = td.cloud().composition().idGas();
const label idL = td.cloud().composition().idLiquid();
const label idS = td.cloud().composition().idSolid();
scalar dhGas =
sum(dMassSRGas)*td.cloud().composition().H(idG, YGas, pc, T);
scalar dhLiquid =
sum(dMassSRLiquid)*td.cloud().composition().H(idL, YLiquid, pc, T);
scalar dhSolid =
sum(dMassSRSolid)*td.cloud().composition().H(idS, YSolid, pc, T);
dhTrans += dhGas + dhLiquid + dhSolid;
// Heat of reaction divided between particle and carrier phase by the
// fraction fh and (1 - fh)
dhTrans += (1.0 - td.constProps().hRetentionCoeff())*HReaction;
return td.constProps().hRetentionCoeff()*HReaction;
return HReaction;
}

27
src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
View File

@ -71,11 +71,11 @@ public:
// IDs of phases in ReacingParcel phase list (Y)
static const label GAS;
static const label LIQUID;
static const label SOLID;
static const label LIQ;
static const label SLD;
//- Class to hold reacting particle constant properties
//- Class to hold reacting multiphase particle constant properties
class constantProperties
:
public ReactingParcel<ParcelType>::constantProperties
@ -177,6 +177,18 @@ private:
const label idS
) const;
//- Return the mixture effective latent heat
template<class TrackData>
scalar LEff
(
TrackData& td,
const scalar p,
const scalar T,
const label idG,
const label idL,
const label idS
) const;
//- Update the mass fractions (Y, YGas, YLiquid, YSolid)
scalar updateMassFractions
(
@ -223,7 +235,7 @@ protected:
const scalarField& YVolatile, // volatile component mass fractions
bool& canCombust, // 'can combust' flag
scalarField& dMassDV // mass transfer - local to particle
);
) const;
//- Calculate surface reactions
template<class TrackData>
@ -237,16 +249,15 @@ protected:
const scalar mass, // mass
const bool canCombust, // 'can combust' flag
const scalarField& dMassVolatile, // mass transfer of volatiles
const scalarField& YMix, // mixture mass fractions
const scalarField& YMix, // mixture mass fractions
const scalarField& YGas, // gas-phase mass fractions
const scalarField& YLiquid,// liquid-phase mass fractions
const scalarField& YSolid, // solid-phase mass fractions
scalarField& dMassSRGas, // gas-phase mass transfer - local
scalarField& dMassSRLiquid,// liquid-phase mass transfer - local
scalarField& dMassSRSolid, // solid-phase mass transfer - local
scalarField& dMassSRCarrier,// carrier phase mass transfer
scalar& dhTrans // enthalpy transfer to carrier phase
);
scalarField& dMassSRCarrier// carrier phase mass transfer
) const;
public:

12
src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
View File

@ -70,8 +70,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
// scale the mass fractions
const scalarField& YMix = this->Y_;
YGas_ /= YMix[GAS] + ROOTVSMALL;
YLiquid_ /= YMix[LIQUID] + ROOTVSMALL;
YSolid_ /= YMix[SOLID] + ROOTVSMALL;
YLiquid_ /= YMix[LIQ] + ROOTVSMALL;
YSolid_ /= YMix[SLD] + ROOTVSMALL;
}
// Check state of Istream
@ -153,7 +153,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingMultiphaseParcel<ParcelType>& p = iter();
p.YLiquid_[j] = YLiquid[i++]/(p.Y()[LIQUID] + ROOTVSMALL);
p.YLiquid_[j] = YLiquid[i++]/(p.Y()[LIQ] + ROOTVSMALL);
}
}
// Populate YSolid for each parcel
@ -172,7 +172,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingMultiphaseParcel<ParcelType>& p = iter();
p.YSolid_[j] = YSolid[i++]/(p.Y()[SOLID] + ROOTVSMALL);
p.YSolid_[j] = YSolid[i++]/(p.Y()[SLD] + ROOTVSMALL);
}
}
}
@ -235,7 +235,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
forAllConstIter(typename Cloud<ParcelType>, c, iter)
{
const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[LIQUID];
YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[LIQ];
}
YLiquid.write();
@ -259,7 +259,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
forAllConstIter(typename Cloud<ParcelType>, c, iter)
{
const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
YSolid[i++] = p0.YSolid()[j]*p0.Y()[SOLID];
YSolid[i++] = p0.YSolid()[j]*p0.Y()[SLD];
}
YSolid.write();

76
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
View File

@ -102,23 +102,18 @@ void Foam::ReactingParcel<ParcelType>::calc
const scalar cp0 = this->cp_;
const scalar mass0 = this->mass();
// Intialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum
vector dUTrans = vector::zero;
// Enthalpy
scalar dhTrans = 0.0;
// Mass transfer due to phase change
scalarField dMassPC(Y_.size(), 0.0);
// Intial ethalpy state
scalar H0H = td.cloud().composition().H(0, Y_, pc_, T0);
scalar H0L = td.cloud().composition().L(0, Y_, pc_, T0);
scalar H0 = H0H - H0L;
// Phase change
// ~~~~~~~~~~~~
// Mass transfer due to phase change
scalarField dMassPC(Y_.size(), 0.0);
// Return enthalpy source and calc mass transfer due to phase change
scalar ShPC =
calcPhaseChange(td, dt, cellI, d0, T0, U0, 0, 1.0, Y_, dMassPC);
@ -126,54 +121,30 @@ void Foam::ReactingParcel<ParcelType>::calc
// Update particle component mass fractions
updateMassFraction(mass0, dMassPC, Y_);
// Update mass
scalar mass1 = mass0 - sum(dMassPC);
// Heat transfer
// ~~~~~~~~~~~~~
// Calculate new particle temperature
scalar htc = 0.0;
scalar T1 =
calcHeatTransfer
(
td,
dt,
cellI,
d0,
U0,
rho0,
T0,
cp0,
ShPC,
htc,
dhTrans
);
scalar T1 = calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, ShPC);
// Calculate new enthalpy state
scalar H1H = td.cloud().composition().H(0, Y_, pc_, T1);
scalar H1L = td.cloud().composition().L(0, Y_, pc_, T1);
scalar H1 = H1H - H1L;
// Motion
// ~~~~~~
// Update mass
scalar mass1 = mass0 - sum(dMassPC);
// No additional forces
vector Fx = vector::zero;
// Calculate new particle velocity
scalar Cud = 0.0;
vector U1 =
calcVelocity
(
td,
dt,
cellI,
d0,
U0,
rho0,
0.5*(mass0 + mass1),
Fx,
Cud,
dUTrans
);
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
// Accumulate carrier phase source terms
@ -188,16 +159,10 @@ void Foam::ReactingParcel<ParcelType>::calc
}
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
// Update enthalpy transfer
td.cloud().hTrans()[cellI] += np0*dhTrans;
// Coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
td.cloud().hTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
}
@ -216,8 +181,7 @@ void Foam::ReactingParcel<ParcelType>::calc
td.cloud().rhoTrans(id)[cellI] += np0*mass1*Y_[i];
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
scalar HEff = td.cloud().composition().H(0, Y_, pc_, T1);
td.cloud().hTrans()[cellI] += np0*mass1*HEff;
td.cloud().hTrans()[cellI] += np0*mass1*H1;
}
}

58
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
View File

@ -84,12 +84,8 @@ void Foam::ThermoParcel<ParcelType>::calc
const scalar cp0 = this->cp_;
const scalar mass0 = this->mass();
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
vector dUTrans = vector::zero;
scalar dhTrans = 0.0;
// Initial enthalpy state
scalar H0 = cp0*T0;
// Heat transfer
@ -99,22 +95,10 @@ void Foam::ThermoParcel<ParcelType>::calc
scalar Sh = 0.0;
// Calculate new particle velocity
scalar htc = 0.0;
scalar T1 =
calcHeatTransfer
(
td,
dt,
cellI,
d0,
U0,
rho0,
T0,
cp0,
Sh,
htc,
dhTrans
);
scalar T1 = calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh);
// Calculate new enthalpy state
scalar H1 = cp0 - T1;
// Motion
@ -124,9 +108,7 @@ void Foam::ThermoParcel<ParcelType>::calc
vector Fx = vector::zero;
// Calculate new particle velocity
scalar Cud = 0.0;
vector U1 =
calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx, Cud, dUTrans);
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
// Accumulate carrier phase source terms
@ -134,16 +116,10 @@ void Foam::ThermoParcel<ParcelType>::calc
if (td.cloud().coupled())
{
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
td.cloud().UTrans()[cellI] += np0*mass0*(U0 - U1);
// Update enthalpy transfer
td.cloud().hTrans()[cellI] += np0*dhTrans;
// Coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
td.cloud().hTrans()[cellI] += np0*mass0*(H0 - H1);
}
// Set new particle properties
@ -165,19 +141,16 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
const scalar rho,
const scalar T,
const scalar cp,
const scalar Sh,
scalar& htc,
scalar& dhTrans
const scalar Sh
)
{
if (!td.cloud().heatTransfer().active())
{
htc = 0.0;
return T;
}
// Calc heat transfer coefficient
htc = td.cloud().heatTransfer().h
scalar htc = td.cloud().heatTransfer().h
(
d,
U - this->Uc_,
@ -192,11 +165,8 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
// Determine new particle temperature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Particle area
scalar Ap = this->areaS(d);
// Determine ap and bp coefficients
scalar ap = Tc_ + Sh/(htc*Ap + ROOTVSMALL);
scalar ap = Tc_ + Sh/(htc*this->areaS(d) + ROOTVSMALL);
scalar bp = 6.0*htc/(rho*d*cp);
if (td.cloud().radiation())
{
@ -219,10 +189,6 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
IntegrationScheme<scalar>::integrationResult Tres =
td.cloud().TIntegrator().integrate(T, dt, ap, bp);
// Enthalpy transfer
// - Using average particle temperature
dhTrans = dt*Ap*htc*(Tres.average() - Tc_);
return Tres.value();
}

4
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
View File

@ -222,9 +222,7 @@ protected:
const scalar rho, // density
const scalar T, // temperature
const scalar cp, // specific heat capacity
const scalar Sh, // additional enthalpy sources
scalar& htc, // heat transfer coeff
scalar& dhTrans // enthalpy transfer to carrier phase
const scalar Sh // additional enthalpy sources
);