From 531e22bc24dc2014b9ed905dd080153031dba1fb Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 19 Aug 2008 16:22:27 +0100 Subject: [PATCH] updates to the reacting parcel/cloud --- .../Templates/ReactingParcel/ReactingParcel.C | 146 +++++++++--------- .../Templates/ReactingParcel/ReactingParcel.H | 11 +- .../ReactingParcel/ReactingParcelI.H | 11 +- .../CompositionModel/CompositionModel.H | 9 +- .../SingleMixtureFraction.C | 23 +++ .../SingleMixtureFraction.H | 5 +- .../NoSurfaceReaction/NoSurfaceReaction.C | 10 +- .../NoSurfaceReaction/NoSurfaceReaction.H | 6 +- .../SurfaceReactionModel.H | 6 +- 9 files changed, 139 insertions(+), 88 deletions(-) diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C index d9d47894a5..4d1702817d 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C @@ -59,7 +59,6 @@ void Foam::ReactingParcel::calcCoupled const scalar mass0 = this->mass(); const scalar np0 = this->nParticle_; const scalar T0 = this->T_; - const scalar cp0 = this->cp_; // ~~~~~~~~~~~~~~~~~~~~~~~~~ // Initialise transfer terms @@ -79,9 +78,12 @@ void Foam::ReactingParcel::calcCoupled // - components do not necessarily exist in particle already scalarList dMassSR(td.cloud().gases().size(), 0.0); - // Total mass lost from particle due to surface reactions - // - sub-model will adjust component mass fractions - scalar dMassMTSR = 0.0; + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Calculate velocity - update U + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + scalar Cud = 0.0; + const vector U1 = calcVelocity(td, dt, Cud, dUTrans); // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -91,13 +93,6 @@ void Foam::ReactingParcel::calcCoupled scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans); - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // Calculate velocity - update U - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - scalar Cud = 0.0; - const vector U1 = calcVelocity(td, dt, Cud, dUTrans); - - // ~~~~~~~~~~~~~~~~~~~~~~~ // Calculate mass transfer // ~~~~~~~~~~~~~~~~~~~~~~~ @@ -107,27 +102,36 @@ void Foam::ReactingParcel::calcCoupled // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Calculate surface reactions // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ - calcSurfaceReactions(td, dt, celli, T0, T1, dMassMTSR, dMassSR); + + // Initialise enthalpy retention to zero + scalar dhRet = 0.0; + + calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet); // New total mass - const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR; + const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR); - // Ratio of mass devolatilised to the total volatile mass of the particle - const scalar fVol = 1 - - (YMixture_[0]*mass1) - /(td.cloud().composition().YMixture0()[0]*mass0_); + // Correct particle temperature to account for latent heat of + // devolatilisation + T1 -= + td.constProps().Ldevol() + *sum(dMassMT) + /(0.5*(mass0 + mass1)*this->cp_); - // Specific heat capacity of non-volatile components - const scalar cpNonVolatile = - ( - YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_) - + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_) - )/(YMixture_[1] + YMixture_[2]); + // Add retained enthalpy from surface reaction to particle and remove + // from gas + T1 += dhRet/(0.5*(mass0 + mass1)*this->cp_); + dhTrans -= dhRet; - // New specific heat capacity - linear variation until volatiles - // have evolved - const scalar cp1 = (cpNonVolatile - td.constProps().cp0())*fVol - + td.constProps().cp0(); + // Correct dhTrans to account for enthalpy of evolved volatiles + dhTrans += + sum(dMassMT) + *td.cloud().composition().HGas(YGas_, 0.5*(T0 + T1)); + + // Correct dhTrans to account for enthalpy of consumed solids + dhTrans += + sum(dMassSR) + *td.cloud().composition().HSolid(YSolid_, 0.5*(T0 + T1)); // ~~~~~~~~~~~~~~~~~~~~~~~ @@ -147,8 +151,8 @@ void Foam::ReactingParcel::calcCoupled td.cloud().UCoeff()[celli] += np0*mass0*Cud; // Update enthalpy transfer -// td.cloud().hTrans()[celli] += np0*(mass0*cp0*T0 - mass1*cp1*T1); - td.cloud().hTrans()[celli] += np0*((mass0*cp0 - mass1*cp1)*T0 + dhTrans); + // - enthalpy of lost solids already accounted for + td.cloud().hTrans()[celli] += np0*dhTrans; // Accumulate coefficient to be applied in carrier phase enthalpy coupling td.cloud().hCoeff()[celli] += np0*htc*this->areaS(); @@ -166,7 +170,12 @@ void Foam::ReactingParcel::calcCoupled { td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i]; } - td.cloud().hTrans()[celli] += np0*mass1*cp1*T1; + td.cloud().hTrans()[celli] += + np0*mass1 + *( + YMixture_[0]*td.cloud().composition().HGas(YGas_, T1) + + YMixture_[2]*td.cloud().composition().HSolid(YSolid_, T1) + ); td.cloud().UTrans()[celli] += np0*mass1*U1; } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -176,7 +185,10 @@ void Foam::ReactingParcel::calcCoupled { this->U_ = U1; this->T_ = T1; - this->cp_ = cp1; + this->cp_ = + YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1) + + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1) + + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_); // Update particle density or diameter if (td.cloud().massTransfer().changesVolume()) @@ -205,7 +217,7 @@ void Foam::ReactingParcel::calcUncoupled // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ const scalar T0 = this->T_; const scalar mass0 = this->mass(); -// const scalar cp0 = this->cp(); + const scalar cp0 = this->cp_; // ~~~~~~~~~~~~~~~~~~~~~~~~~ // Initialise transfer terms @@ -225,9 +237,12 @@ void Foam::ReactingParcel::calcUncoupled // - components do not necessarily exist in particle already scalarList dMassSR(td.cloud().gases().size(), 0.0); - // Total mass lost from particle due to surface reactions - // - sub-model will adjust component mass fractions - scalar dMassMTSR = 0.0; + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Calculate velocity - update U + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + scalar Cud = 0.0; + const vector U1 = calcVelocity(td, dt, Cud, dUTrans); // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -240,13 +255,6 @@ void Foam::ReactingParcel::calcUncoupled T1 = min(td.constProps().Tvap(), T1); - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // Calculate velocity - update U - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - scalar Cud = 0.0; - const vector U1 = calcVelocity(td, dt, Cud, dUTrans); - - // ~~~~~~~~~~~~~~~~~~~~~~~ // Calculate mass transfer // ~~~~~~~~~~~~~~~~~~~~~~~ @@ -256,37 +264,23 @@ void Foam::ReactingParcel::calcUncoupled // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Calculate surface reactions // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ - calcSurfaceReactions - ( - td, - dt, - celli, - T0, - T1, - dMassMTSR, - dMassSR - ); + + // Initialise enthalpy retention to zero + scalar dhRet = 0.0; + + calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet); // New total mass - const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR; + const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR); - // Ratio of mass devolatilised to the total volatile mass of the particle - const scalar fVol = 1 - - (YMixture_[0]*mass1) - /(td.cloud().composition().YMixture0()[0]*mass0_); - - // Specific heat capacity of non-volatile components - const scalar cpNonVolatile = - ( - YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_) - + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_) - )/(YMixture_[1] + YMixture_[2]); - - // New specific heat capacity - linear variation until volatiles - // have evolved - const scalar cp1 = (cpNonVolatile - td.constProps().cp0())*fVol - + td.constProps().cp0(); + // New specific heat capacity + const scalar cp1 = + YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1) + + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1) + + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_); + // Add retained enthalpy to particle + T1 += dhRet/(mass0*0.5*(cp0 + cp1)); // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Remove the particle when mass falls below minimum threshold @@ -389,8 +383,9 @@ void Foam::ReactingParcel::calcSurfaceReactions const label celli, const scalar T0, const scalar T1, - scalar& dMassMTSR, - scalarList& dMassMT + const scalarList& dMassMT, + scalarList& dMassSR, + scalar& dhRet ) { // Check that model is active @@ -409,21 +404,20 @@ void Foam::ReactingParcel::calcSurfaceReactions T0, T1, this->Tc_, + pc_, this->rhoc_, this->mass(), + dMassMT, YGas_, YLiquid_, YSolid_, YMixture_, - dMassMTSR, - dMassMT + dMassSR, + dhRet ); } -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - - // * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * // #include "ReactingParcelIO.C" diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H index 43e125f13b..a6b8eb5ad5 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H @@ -89,6 +89,9 @@ public: //- Boiling point [K] const scalar Tbp_; + //- Latent heat of devolatilisation [J/kg] + const scalar Ldevol_; + public: @@ -102,6 +105,9 @@ public: //- Return const access to the boiling point inline scalar Tbp() const; + + //- Return const access to the latent heat of devolatilisation + inline scalar Ldevol() const; }; @@ -210,8 +216,9 @@ protected: const label celli, const scalar T0, const scalar T1, - scalar& dMassMTSR, - scalarList& dMassMT + const scalarList& dMassMT, + scalarList& dMassSR, + scalar& dhRet ); diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H index 7138666b2e..0c50b0a853 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H @@ -34,7 +34,8 @@ inline Foam::ReactingParcel::constantProperties::constantProperties : ThermoParcel::constantProperties(dict), Tvap_(dimensionedScalar(dict.lookup("Tvap")).value()), - Tbp_(dimensionedScalar(dict.lookup("Tbp")).value()) + Tbp_(dimensionedScalar(dict.lookup("Tbp")).value()), + Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value()) {} @@ -127,6 +128,14 @@ Foam::ReactingParcel::constantProperties::Tbp() const } +template +inline Foam::scalar +Foam::ReactingParcel::constantProperties::Ldevol() const +{ + return Ldevol_; +} + + // * * * * * * * * * * * trackData Member Functions * * * * * * * * * * * * // template diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H index ce85b19a5c..9caea94a34 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H +++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H @@ -209,13 +209,20 @@ public: //- Return the gas constant for the gas mixture virtual scalar RGas(const scalarField& YGas) const = 0; - //- Return enthalpy for the gas mixture + //- Return enthalpy for the gas mixture [energy per unit mass] virtual scalar HGas ( const scalarField& YGas, const scalar T ) const = 0; + //- Return enthalpy for the solid mixture [energy per unit mass] + virtual scalar HSolid + ( + const scalarField& YSolid, + const scalar T + ) const = 0; + //- Return specific heat caparcity for the gas mixture virtual scalar cpGas ( diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C index 350546b5fe..74fb92338c 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C +++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C @@ -519,6 +519,29 @@ const } +template +Foam::scalar Foam::SingleMixtureFraction::HSolid +( + const scalarField& YSolid, + const scalar T +) +const +{ + scalar HMixture = 0.0; + forAll(YSolid, i) + { + label id = solidGlobalIds_[i]; + HMixture += + YSolid[i] + *( + this->solids().properties()[id].Hf() + + this->solids().properties()[id].cp()*T + ); + } + return HMixture; +} + + template Foam::scalar Foam::SingleMixtureFraction::cpGas ( diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H index bd03d69456..de1fb8daf7 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H +++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H @@ -188,9 +188,12 @@ public: //- Return the gas constant for the gas mixture scalar RGas(const scalarField& YGas) const; - //- Return enthalpy for the gas mixture + //- Return enthalpy for the gas mixture [energy per unit mass] scalar HGas(const scalarField& YGas, const scalar T) const; + //- Return enthalpy for the solid mixture [energy per unit mass] + scalar HSolid(const scalarField& YSolid, const scalar T) const; + //- Return specific heat caparcity for the gas mixture scalar cpGas(const scalarField& YGas, const scalar T) const; diff --git a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C index 8ea8816f04..d164a500c2 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C +++ b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C @@ -64,16 +64,20 @@ void Foam::NoSurfaceReaction::calculate const scalar T0, const scalar T1, const scalar Tc, + const scalar pc, const scalar rhoc, const scalar massp, + const scalarList& dMassMT, scalarField& YGas, scalarField& YLiquid, scalarField& YSolid, scalarField& YMixture, - scalar& dMassMTSR, - scalarList& dMassSR + scalarList& dMassSR, + scalar& dhRet ) const -{} +{ + // do nothing +} // ************************************************************************* // diff --git a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H index f3072df598..8a58352b5d 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H +++ b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H @@ -84,14 +84,16 @@ public: const scalar T0, const scalar T1, const scalar Tc, + const scalar pc, const scalar rhoc, const scalar massp, + const scalarList& dMassMT, scalarField& YGas, scalarField& YLiquid, scalarField& YSolid, scalarField& YMixture, - scalar& dMassMTSR, - scalarList& dMassSR + scalarList& dMassSR, + scalar& dhRet ) const; }; diff --git a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H index b50b5a90b4..a3b3e4e00f 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H +++ b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H @@ -141,14 +141,16 @@ public: const scalar T0, const scalar T1, const scalar Tc, + const scalar pc, const scalar rhoc, const scalar massp, + const scalarList& dMassMT, scalarField& YGas, scalarField& YLiquid, scalarField& YSolid, scalarField& YMixture, - scalar& dMassMTSR, - scalarList& dMassSR + scalarList& dMassSR, + scalar& dhRet ) const = 0; };