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

multiple updates

Browse Source
This commit is contained in:
andy
2009-03-16 17:43:12 +00:00
parent 0ccaf7aae7
commit c5c16150b7
23 changed files with 263 additions and 252 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -232,7 +232,7 @@ void Foam::KinematicCloud<ParcelType>::evolve()
if (debug)
{
this->dumpParticlePositions();
this->writePositions();
}
if (coupled_)
@ -261,25 +261,4 @@ void Foam::KinematicCloud<ParcelType>::info() const
}
template<class ParcelType>
void Foam::KinematicCloud<ParcelType>::dumpParticlePositions() const
{
OFstream pObj
(
this->db().time().path()/"parcelPositions_"
+ this->name() + "_"
+ name(this->injection().nInjections()) + ".obj"
);
forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)
{
const ParcelType& p = iter();
pObj<< "v " << p.position().x() << " " << p.position().y() << " "
<< p.position().z() << nl;
}
pObj.flush();
}
// ************************************************************************* //

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

@ -338,9 +338,6 @@ public:
//- Print cloud information
void info() const;
//- Dump particle positions to .obj file
void dumpParticlePositions() const;
// Fields

2
src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
View File

@ -204,7 +204,7 @@ void Foam::ReactingCloud<ParcelType>::evolve()
if (debug)
{
this->dumpParticlePositions();
this->writePositions();
}
if (this->coupled())

2
src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
View File

@ -181,7 +181,7 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::evolve()
if (debug)
{
this->dumpParticlePositions();
this->writePositions();
}
if (this->coupled())

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

@ -203,7 +203,7 @@ void Foam::ThermoCloud<ParcelType>::evolve()
if (debug)
{
this->dumpParticlePositions();
this->writePositions();
}
if (this->coupled())

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

@ -27,8 +27,6 @@ License
#include "KinematicParcel.H"
#include "dimensionedConstants.H"
#include "fvcGrad.H"
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class ParcelType>
@ -66,30 +64,43 @@ void Foam::KinematicParcel<ParcelType>::calc
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 1. Calculate velocity - update U, dUTrans
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
// Define local properties at beginning of time step
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar mass0 = mass();
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum
vector dUTrans = vector::zero;
const vector U1 = calcVelocity(td, dt, vector::zero, mass(), Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 2. Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Motion
// ~~~~~~
// No additional forces
vector Fx = vector::zero;
// Calculate new particle velocity
scalar Cud = 0.0;
vector U1 = calcVelocity(td, dt, cellI, Fx, mass0, Cud, dUTrans);
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().coupled())
{
// Update momentum transfer
td.cloud().UTrans()[cellI] += nParticle_*dUTrans;
// Coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += nParticle_*mass()*Cud;
td.cloud().UCoeff()[cellI] += nParticle_*mass0*Cud;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
U_ = U1;
}
@ -100,10 +111,11 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
TrackData& td,
const scalar dt,
const label cellI,
const vector& Fx,
const scalar mass,
scalar& Cud,
vector& dUTrans,
vector& dUTrans
) const
{
// Return linearised term from drag model
@ -127,12 +139,13 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
// Pressure gradient force
if (td.cloud().forcePressureGradient())
{
Ftot += rhoc_/rho_*(U_ & fvc::grad(Uc_))
const vector& d = this->mesh().deltaCoeffs()[cellI];
Ftot += rhoc_/rho_*(U_ & (d^Uc_));
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle velocity
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// New particle velocity
//~~~~~~~~~~~~~~~~~~~~~~
// Update velocity - treat as 3-D
const vector ap = Uc_ + (Ftot + Fx)/(Cud + VSMALL);
@ -142,8 +155,6 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
// Calculate the momentum transfer to the continuous phase
// - do not include gravity impulse
// TODO: USE AVERAGE PARTICLE MASS
dUTrans = -mass*(Unew - U_ - dt*td.g());
return Unew;

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

@ -84,7 +84,6 @@ public:
//- Class to hold kinematic particle constant properties
class constantProperties
{
// Private data
//- Particle density [kg/m3] (constant)
@ -120,7 +119,6 @@ public:
:
public Particle<ParcelType>::trackData
{
// Private data
//- Reference to the cloud containing this particle
@ -234,6 +232,7 @@ protected:
(
TrackData& td,
const scalar dt,
const label cellI,
const vector& Fx,
const scalar mass,
scalar& Cud,

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

@ -84,7 +84,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Define local properties at beginning of timestep
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar mass0 = this->mass();
@ -92,97 +91,94 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const scalar T0 = this->T_;
const scalar pc = this->pc_;
scalarField& YMix = this->Y_;
const label idG = td.cloud().composition().idGas();
const label idL = td.cloud().composition().idLiquid();
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate phase change in liquid phase
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Mass transfer from particle to carrier phase
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum
vector dUTrans = vector::zero;
// Enthalpy
scalar dhTrans = 0.0;
// Mass transfer due to phase change
scalarList dMassPC(td.cloud().gases().size(), 0.0);
scalar shPC =
// Mass transfer due to devolatilisation
scalarList dMassDV(td.cloud().gases().size(), 0.0);
// Change in carrier phase composition due to surface reactions
scalarList dMassSRc(td.cloud().gases().size(), 0.0);
// Phase change in liquid phase
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Return enthalpy source and calc mass transfer due to phase change
scalar ShPC =
calcPhaseChange(td, dt, cellI, T0, idL, YMix[idL], YLiquid_, dMassPC);
// Update particle component mass fractions
updateMassFraction(mass0, dMassPC, YLiquid_);
this->updateMassFraction(mass0, dMassPC, YLiquid_);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate Devolatilisation
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
// Mass transfer from particle to carrier phase
// - components exist in particle already
scalarList dMassDV(td.cloud().gases().size(), 0.0);
scalar shDV = calcDevolatilisation(td, dt, T0, mass0, idGas, YMix, dMassDV);
// Devolatilisation
// ~~~~~~~~~~~~~~~~
// Return enthalpy source and calc mass transfer due to devolatilisation
scalar ShDV = calcDevolatilisation(td, dt, T0, mass0, idG, YMix, dMassDV);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate surface reactions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Surface reactions
// ~~~~~~~~~~~~~~~~~
// Mass transfer of volatile components from particle to carrier phase
const scalarList dMassMT = dMassPC + dMassDV;
// Mass transfer due to surface reactions from particle to carrier phase
// - components do not necessarily exist in particle already
scalarList dMassSR(td.cloud().gases().size(), 0.0);
// Initialise enthalpy retention to zero
scalar dhRet = 0.0;
calcSurfaceReactions(td, dt, cellI, T0, dMassMT, dMassSR, dhRet);
// Enthalpy retention divided between particle and carrier phase by the
// fraction fh and (1 - fh)
scalar ShSR = td.constProps().fh()*dhRet;
dhTrans -= (1.0 - td.constProps().fh())*dhRet;
// Correct dhTrans to account for enthalpy of consumed solids
dhTrans +=
sum(dMassSR)*td.cloud().composition().H(idSolid, YSolid_, pc, T0);
// Correct dhTrans to account for enthalpy of evolved volatiles
dhTrans +=
sum(dMassMT)*td.cloud().composition().H(idGas, YGas_, pc, T0);
// Return enthalpy source and calc mass transfer(s) due to surface reaction
scalar ShSR =
calcSurfaceReactions(td, dt, cellI, T0, dMassMT, dMassSRc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~
// Calculate heat transfer
// ~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
// Heat transfer
// ~~~~~~~~~~~~~
// Total enthalpy source
scalar Sh = ShPC + ShDV + ShSR;
scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, shHT);
// Calculate new particle temperature
scalar htc = 0.0;
scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~
// Calculate velocity
// ~~~~~~~~~~~~~~~~~~
// Update mass
scalar mass1 = mass0 - massPC - massD - massSR;
scalar Cud = 0.0;
vector dUTrans = vector::zero;
// Motion
// ~~~~~~
// Update mass (not to include cMassSR)
scalar mass1 = mass0 - sum(dMassPC) - sum(dMassDV);
// No additional forces
vector Fx = vector::zero;
vector U1 = calcVelocity(td, dt, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
// Calculate new particle velocity
scalar Cud = 0.0;
vector U1 =
calcVelocity(td, dt, cellI, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Collect contributions to determine new particle thermo properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Update specific heat capacity
cp1 = cpEff(td, pc, T1);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().coupled())
{
// Transfer mass lost from particle to carrier mass source
forAll(dMassMT, i)
{
td.cloud().rhoTrans(i)[cellI] +=
np0*(dMassPC[i] + dMassDV[i] + dMassSR[i]);
np0*(dMassPC[i] + dMassDV[i] + dMassSRc[i]);
}
// Update momentum transfer
@ -200,9 +196,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mass1 < td.constProps().minParticleMass())
{
td.keepParticle = false;
@ -215,14 +211,16 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
td.cloud().hTrans()[cellI] += np0*mass1*HEff(td, pc, T1);
td.cloud().UTrans()[cellI] += np0*mass1*U1;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
else
{
this->U_ = U1;
this->T_ = T1;
this->cp_ = cp1;
this->cp_ = cpEff(td, pc, T1);
// Update particle density or diameter
if (td.constProps().constantVolume())
@ -246,8 +244,8 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
const scalar T,
const scalar mass,
const label idVolatile,
scalarField_ YMixture,
scalarList& dMassMT
scalarField& YMixture,
scalarList& dMassDV
)
{
// Check that model is active, and that the parcel temperature is
@ -268,16 +266,16 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
dt,
this->mass0_,
mass,
td.cloud().composition().YMixture0()[idVolatile],
YMix[0],
T,
td.cloud().composition().YMixture0()[idVolatile],
YMixture[0],
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];
YMixture[0] = (YMixture[0]*mass - dMassTot)/(mass - dMassTot);
YMixture[1] = YMixture[1]*mass/(mass - dMassTot);
YMixture[2] = 1.0 - YMixture[0] - YMixture[1];
// Add to cummulative mass transfer
forAll (YGas_, i)
@ -286,26 +284,26 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
// Volatiles mass transfer
scalar volatileMass = YGas_[i]*dMassTot;
dMassMT[id] += volatileMass;
dMassDV[id] += volatileMass;
}
td.cloud().addToMassDevolatilisation(dMassTot);
return = td.constProps().Ldevol()*dMassTot;
return td.constProps().Ldevol()*dMassTot;
}
template<class ParcelType>
template<class TrackData>
void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
(
TrackData& td,
const scalar dt,
const label cellI,
const scalar T,
const scalarList& dMassMT,
scalarList& dMassSR,
scalar& dhRet
scalarList& dMassSRc,
scalar& dhTrans
)
{
// Check that model is active
@ -316,7 +314,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
// Update mass transfer(s)
// - Also updates Y()'s
td.cloud().surfaceReaction().calculate
scalar HReaction = td.cloud().surfaceReaction().calculate
(
dt,
cellI,
@ -331,11 +329,27 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
YLiquid_,
YSolid_,
this->Y_,
dMassSR,
dhRet
dMassSRc,
HReaction
);
// Heat of reaction divided between particle and carrier phase by the
// fraction fh and (1 - fh)
dhTrans -= (1.0 - td.constProps().fh())*HReaction;
/*
// Correct dhTrans to account for enthalpy of consumed solids
dhTrans +=
sum(dMassSR)*td.cloud().composition().H(idSolid, YSolid_, pc, T0);
// Correct dhTrans to account for enthalpy of evolved volatiles
dhTrans +=
sum(dMassMT)*td.cloud().composition().H(idGas, YGas_, pc, T0);
// TODO: td.cloud().addToMassSurfaceReaction(sum(dMassSR));
*/
return td.constProps().fh()*HReaction;
}

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

@ -185,23 +185,23 @@ protected:
TrackData& td,
const scalar dt,
const scalar T,
const scalar mass,
const label idVolatile,
scalarList& dMassMT
const scalar mass,
const label idVolatile,
scalarField& YMixture,
scalarList& dMassDV
);
//- Calculate surface reactions
template<class TrackData>
void calcSurfaceReactions
scalar calcSurfaceReactions
(
TrackData& td,
const scalar dt,
const label cellI,
const scalar T0,
const scalar T1,
const scalar T,
const scalarList& dMassMT,
scalarList& dMassSR,
scalar& dhRet
scalarList& dMassSRc,
scalar& dhTrans
);

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

@ -69,7 +69,6 @@ void Foam::ReactingParcel<ParcelType>::calc
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Define local properties at beginning of time step
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar mass0 = this->mass();
@ -77,51 +76,54 @@ void Foam::ReactingParcel<ParcelType>::calc
const scalar T0 = this->T_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// 1. Calculate phase change
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Mass transfer from particle to carrier phase
// Intialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum
vector dUTrans = vector::zero;
// Enthalpy
scalar dhTrans = 0.0;
// Mass transfer due to phase change
scalarList dMassPC(td.cloud().gases().size(), 0.0);
// Phase change
// ~~~~~~~~~~~~
// Return enthalpy source and calc mass transfer due to phase change
scalar ShPC = calcPhaseChange(td, dt, cellI, T0, 0, 1.0, Y_, dMassPC);
// Update particle component mass fractions
updateMassFraction(mass0, dMassPC, Y_);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
// 2. Calculate heat transfer
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
// Heat transfer
// ~~~~~~~~~~~~~
// Calculate new particle temperature
scalar htc = 0.0;
scalar ShHT = 0.0;
scalar T1 = calcHeatTransfer(td, dt, cellI, ShPC, htc, ShHT);
scalar T1 = calcHeatTransfer(td, dt, cellI, ShPC, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~
// 3. Calculate velocity
// ~~~~~~~~~~~~~~~~~~~~~
// Motion
// ~~~~~~
// Update mass
scalar mass1 = mass0 - sum(dMassPC);
scalar Cud = 0.0;
vector dUTrans = vector::zero;
// No additional forces
vector Fx = vector::zero;
vector U1 = calcVelocity(td, dt, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
// Calculate new particle velocity
scalar Cud = 0.0;
vector U1 =
calcVelocity(td, dt, cellI, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 4. Collect contributions to determine new particle thermo properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Total enthalpy transfer from the particle to the carrier phase
scalar dhTrans = ShHT + ShPC;
// New specific heat capacity of mixture - using new temperature
scalar cp1 = td.cloud().composition().cp(0, Y_, pc_, T1);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 5. Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().coupled())
{
// Transfer mass lost from particle to carrier mass source
@ -137,16 +139,15 @@ void Foam::ReactingParcel<ParcelType>::calc
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
// Update enthalpy transfer
td.cloud().hTrans()[cellI] += np0*dhTrans;
td.cloud().hTrans()[cellI] += np0*(dhTrans + ShPC);
// Coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 6a. Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mass1 < td.constProps().minParticleMass())
{
td.keepParticle = false;
@ -159,18 +160,20 @@ void Foam::ReactingParcel<ParcelType>::calc
td.cloud().rhoTrans(i)[cellI] += np0*mass1*Y_[i];
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
scalar HEff = td.cloud().composition().H(0, YComponents, pc_, T1);
scalar HEff = td.cloud().composition().H(0, Y_, pc_, T1);
td.cloud().hTrans()[cellI] += np0*mass1*HEff;
}
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 6b. Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
else
{
this->U_ = U1;
this->T_ = T1;
this->cp_ = cp1;
this->cp_ = td.cloud().composition().cp(0, Y_, pc_, T1);
// Update particle density or diameter
if (td.constProps().constantVolume())
@ -195,7 +198,7 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
const scalar T,
const label idPhase,
const scalar YPhase,
scalarField& YComponents,
const scalarField& YComponents,
scalarList& dMass
)
{

2
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
View File

@ -190,7 +190,7 @@ protected:
const scalar T,
const label idPhase,
const scalar YPhase,
const scalarField& YComponents,
const scalarField& YComponents,
scalarList& dMass
);

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

@ -53,7 +53,6 @@ void Foam::ThermoParcel<ParcelType>::calc
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Define local properties at beginning of time step
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const vector U0 = this->U_;
@ -61,25 +60,37 @@ void Foam::ThermoParcel<ParcelType>::calc
const scalar np0 = this->nParticle_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 1. Calculate heat transfer - update T
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
scalar ShHT = 0.0;
const scalar T1 = calcHeatTransfer(td, dt, cellI, 0.0, htc, ShHT);
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 2. Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
vector dUTrans = vector::zero;
const vector U1 = calcVelocity(td, dt, vector::zero, Cud, mass0, dUTrans);
scalar dhTrans = 0.0;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Heat transfer
// ~~~~~~~~~~~~~
// No additional enthalpy sources
vector Sh = 0.0;
// Calculate new particle velocity
scalar htc = 0.0;
scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, dhTrans);
// Motion
// ~~~~~~
// No additional forces
vector Fx = vector::zero;
// Calculate new particle velocity
scalar Cud = 0.0;
vector U1 = calcVelocity(td, dt, cellI, Fx, Cud, mass0, dUTrans);
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().coupled())
{
// Update momentum transfer
@ -89,16 +100,15 @@ void Foam::ThermoParcel<ParcelType>::calc
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
// Update enthalpy transfer
td.cloud().hTrans()[cellI] += np0*ShHT;
td.cloud().hTrans()[cellI] += np0*dhTrans;
// Coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 5. Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
this->U() = U1;
// Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
this->U_ = U1;
T_ = T1;
}
@ -118,7 +128,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
if (!td.cloud().heatTransfer().active())
{
htc = 0.0;
dhTrans = 0.0;
ShHT = 0.0;
return T_;
}
@ -134,7 +144,11 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
this->muc_
);
//- Assuming diameter = diameter at start of time step
// Set new particle temperature
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Assuming diameter = diameter at start of time step
scalar Ap = this->areasS();
// Determine ap and bp coefficients
@ -161,7 +175,8 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
IntegrationScheme<scalar>::integrationResult Tres =
td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
// Using average parcel temperature for enthalpy transfer calculation
// Enthalpy transfer
// - Using average particle temperature
ShHT = dt*Ap*htc*(Tres.average() - Tc_);
return Tres.value();

3
src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
View File

@ -199,9 +199,6 @@ void Foam::LiquidEvaporation<CloudType>::calculate
scalarList& dMass
) const
{
// initialise total mass transferred from the particle to carrier phase
scalar dMassTot = 0.0;
// construct carrier phase species volume fractions for cell, cellI
scalarField Xc = calcXc(cellI);

4
src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C
View File

@ -66,9 +66,9 @@ Foam::scalar Foam::ConstantRateDevolatilisation<CloudType>::calculate
const scalar dt,
const scalar mass0,
const scalar mass,
const scalar YVolatile0,
const scalarField& YVolatile,
const scalar T,
const scalar YVolatile0,
const scalar YVolatile,
bool& canCombust
) const
{

4
src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H
View File

@ -93,9 +93,9 @@ public:
const scalar dt,
const scalar mass0,
const scalar mass,
const scalar YVolatile0,
const scalarField& YVolatile,
const scalar T,
const scalar YVolatile0,
const scalar YVolatile,
bool& canCombust
) const;
};

4
src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H
View File

@ -135,9 +135,9 @@ public:
const scalar dt,
const scalar mass0,
const scalar mass,
const scalarField& YMixture0,
const scalarField& YMixture,
const scalar T,
const scalar YVolatile0,
const scalar YVolatile,
bool& canCombust
) const = 0;
};

4
src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C
View File

@ -61,8 +61,8 @@ Foam::scalar Foam::NoDevolatilisation<CloudType>::calculate
const scalar,
const scalar,
const scalar,
const scalarField&,
const scalarField&,
const scalar,
const scalar,
const scalar,
bool& canCombust
) const

4
src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H
View File

@ -76,8 +76,8 @@ public:
const scalar,
const scalar,
const scalar,
const scalarField&,
const scalarField&,
const scalar,
const scalar,
const scalar,
bool&
) const;

4
src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C
View File

@ -69,9 +69,9 @@ Foam::scalar Foam::SingleKineticRateDevolatilisation<CloudType>::calculate
const scalar dt,
const scalar mass0,
const scalar mass,
const scalar YVolatile0,
const scalarField& YVolatile,
const scalar T,
const scalar YVolatile0,
const scalar YVolatile,
bool& canCombust
) const
{

4
src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H
View File

@ -95,9 +95,9 @@ public:
const scalar dt,
const scalar mass0,
const scalar mass,
const scalar YVolatile0,
const scalarField& YVolatile,
const scalar T,
const scalar YVolatile0,
const scalar YVolatile,
bool& canCombust
) const;
};

33
src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
View File

@ -56,27 +56,26 @@ bool Foam::NoSurfaceReaction<CloudType>::active() const
template<class CloudType>
void Foam::NoSurfaceReaction<CloudType>::calculate
Foam::scalar Foam::NoSurfaceReaction<CloudType>::calculate
(
const scalar dt,
const label celli,
const scalar dp,
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,
scalarList& dMassSR,
scalar& dhRet
const scalar,
const label,
const scalar,
const scalar,
const scalar,
const scalar,
const scalar,
const scalar,
const scalarList&,
scalarField&,
scalarField&,
scalarField&,
scalarField&,
scalarList&
) const
{
// do nothing
return 0.0;
}

14
src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
View File

@ -74,24 +74,22 @@ public:
virtual bool active() const;
//- Update surface reactions
virtual void calculate
virtual scalar calculate
(
const scalar dt,
const label celli,
const scalar dp,
const scalar T0,
const scalar T1,
const label cellI,
const scalar d,
const scalar T,
const scalar Tc,
const scalar pc,
const scalar rhoc,
const scalar massp,
const scalar mass,
const scalarList& dMassMT,
scalarField& YGas,
scalarField& YLiquid,
scalarField& YSolid,
scalarField& YMixture,
scalarList& dMassSR,
scalar& dhRet
scalarList& dMassSRc
) const;
};

15
src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
View File

@ -132,24 +132,23 @@ public:
virtual bool active() const = 0;
//- Update surface reactions
virtual void calculate
// Returns the heat of reaction
virtual scalar calculate
(
const scalar dt,
const label celli,
const scalar dp,
const scalar T0,
const scalar T1,
const label cellI,
const scalar d,
const scalar T,
const scalar Tc,
const scalar pc,
const scalar rhoc,
const scalar massp,
const scalar mass,
const scalarList& dMassMT,
scalarField& YGas,
scalarField& YLiquid,
scalarField& YSolid,
scalarField& YMixture,
scalarList& dMassSR,
scalar& dhRet
scalarList& dMassSRc
) const = 0;
};