further updates en route to new evap model(s)

This commit is contained in:
andy
2009-03-05 15:48:51 +00:00
parent d25686db61
commit 7677c87f15
32 changed files with 367 additions and 929 deletions

View File

@ -7,18 +7,21 @@ $(DERIVEDPARCELS)/basicThermoParcel/basicThermoParcel.C
$(DERIVEDPARCELS)/basicReactingParcel/basicReactingParcel.C
$(DERIVEDPARCELS)/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C
/* Cloud base classes */
clouds/baseClasses/kinematicCloud/kinematicCloud.C
clouds/baseClasses/thermoCloud/thermoCloud.C
clouds/baseClasses/reactingCloud/reactingCloud.C
clouds/baseClasses/reactingMultiphaseCloud/reactingMultiphaseCloud.C
/* Cloud container/injection mechanisms */
clouds/derived/basicKinematicCloud/basicKinematicCloud.C
clouds/derived/basicThermoCloud/basicThermoCloud.C
clouds/derived/basicReactingCloud/basicReactingCloud.C
clouds/derived/basicReactingMultiphaseCloud/basicReactingMultiphaseCloud.C
/* kinematic parcel sub-models */
KINEMATICPARCEL=$(DERIVEDPARCELS)/basicKinematicParcel
$(KINEMATICPARCEL)/defineBasicKinematicParcel.C
@ -27,6 +30,7 @@ $(KINEMATICPARCEL)/submodels/makeBasicKinematicParcelDragModels.C
$(KINEMATICPARCEL)/submodels/makeBasicKinematicParcelInjectionModels.C
$(KINEMATICPARCEL)/submodels/makeBasicKinematicParcelWallInteractionModels.C
/* thermo parcel sub-models */
THERMOPARCEL=$(DERIVEDPARCELS)/basicThermoParcel
$(THERMOPARCEL)/defineBasicThermoParcel.C
@ -36,6 +40,7 @@ $(THERMOPARCEL)/submodels/makeBasicThermoParcelHeatTransferModels.C
$(THERMOPARCEL)/submodels/makeBasicThermoParcelInjectionModels.C
$(THERMOPARCEL)/submodels/makeBasicThermoParcelWallInteractionModels.C
/* reacting parcel sub-models */
REACTINGPARCEL=$(DERIVEDPARCELS)/basicReactingParcel
$(REACTINGPARCEL)/defineBasicReactingParcel.C
@ -47,6 +52,7 @@ $(REACTINGPARCEL)/submodels/makeBasicReactingParcelInjectionModels.C
$(REACTINGPARCEL)/submodels/makeBasicReactingParcelPhaseChangeModels.C
$(REACTINGPARCEL)/submodels/makeBasicReactingParcelWallInteractionModels.C
/* reacting multiphase parcel sub-models */
REACTINGMPPARCEL=$(DERIVEDPARCELS)/basicReactingMultiphaseParcel
$(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
@ -60,25 +66,25 @@ $(REACTINGMPPARCEL)/submodels/makeBasicReactingMultiphaseParcelPhaseChangeModels
$(REACTINGMPPARCEL)/submodels/makeBasicReactingMultiphaseParcelSurfaceReactionModels.C
$(REACTINGMPPARCEL)/submodels/makeBasicReactingMultiphaseParcelWallInteractionModels.C
/* bolt-on models */
submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
/* integration schemes */
IntegrationScheme/makeIntegrationSchemes.C
/* phase properties */
phaseProperties/phaseProperties/phaseProperties.C
phaseProperties/phaseProperties/phasePropertiesIO.C
phaseProperties/phasePropertiesList/phasePropertiesList.C
/* evaporation properties */
evaporationProperties/evaporationProperties/evaporationProperties.C
evaporationProperties/evaporationProperties/evaporationPropertiesIO.C
/* data entries */
submodels/IO/DataEntry/makeDataEntries.C
submodels/IO/DataEntry/polynomial/polynomial.C
LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate

View File

@ -1,80 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "evaporationProperties.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::evaporationProperties::evaporationProperties()
:
name_("unknownSpecie"),
Dab_(0.0),
TvsPSat_()
{}
Foam::evaporationProperties::evaporationProperties
(
const evaporationProperties& pp
)
:
name_(pp.name_),
Dab_(pp.Dab_),
TvsPSat_(pp.TvsPSat_)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::evaporationProperties::~evaporationProperties()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::word& Foam::evaporationProperties::name() const
{
return name_;
}
Foam::scalar Foam::evaporationProperties::Dab() const
{
return Dab_;
}
const Foam::DataEntry<Foam::scalar>&
Foam::evaporationProperties::TvsPSat() const
{
return TvsPSat_();
}
// ************************************************************************* //

View File

@ -1,117 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::evaporationProperties
Description
Helper class to manage evaporation properties
SourceFiles
evaporationProperties.C
\*---------------------------------------------------------------------------*/
#ifndef evaporationProperties_H
#define evaporationProperties_H
#include "DataEntry.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class evaporationProperties Declaration
\*---------------------------------------------------------------------------*/
class evaporationProperties
{
// Private data
//- Name of specie
word name_;
//- Diffusion coefficient
scalar Dab_;
//- Data entry for saturation pressure as a function of temperature
autoPtr<DataEntry<scalar> > TvsPSat_;
// Private member functions
public:
// Constructors
//- Null constructor
evaporationProperties();
//- Construct from Istream
evaporationProperties(Istream&);
//- Construct as copy
evaporationProperties(const evaporationProperties&);
//- Destructor
~evaporationProperties();
// Public member functions
// Access
//- Return const access to the specie name
const word& name() const;
//- Return const access to the diffusion coefficient
scalar Dab() const;
//- Return const access to the saturation pressure as a function
// of temperature
const DataEntry<scalar>& TvsPSat() const;
// IOstream Operators
friend Istream& operator>>(Istream&, evaporationProperties&);
friend Ostream& operator<<(Ostream&, const evaporationProperties&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,130 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "evaporationProperties.H"
#include "dictionaryEntry.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::evaporationProperties::evaporationProperties(Istream& is)
:
name_("unknown"),
Dab_(0.0),
TvsPSat_(NULL)
{
is.check
(
"Foam::evaporationProperties::evaporationProperties(Istream& is)"
);
dictionaryEntry evapInfo(dictionary::null, is);
if (!evapInfo.isDict())
{
FatalErrorIn
(
"Foam::evaporationProperties::evaporationProperties(Istream& is)"
) << "Evaporation properties should be given in dictionary entries"
<< nl << exit(FatalError);
}
name_ = evapInfo.keyword();
evapInfo.lookup("Dab") >> Dab_;
TvsPSat_.reset(DataEntry<scalar>::New("TvsPSat", evapInfo.dict()).ptr());
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Istream& Foam::operator>>(Istream& is, evaporationProperties& ep)
{
is.check
(
"Foam::Istream& Foam::operator>>"
"("
"Foam::Istream&,"
"Foam::evaporationProperties&"
")"
);
dictionaryEntry evapInfo(dictionary::null, is);
if (!evapInfo.isDict())
{
FatalErrorIn
(
"Foam::Istream& Foam::operator>>"
"("
"Istream& is,"
"evaporationProperties& pp"
")"
) << "Evaporation properties should be given in dictionary entries"
<< nl << exit(FatalError);
}
ep.name_ = evapInfo.keyword();
evapInfo.lookup("Dab") >> ep.Dab_;
ep.TvsPSat_.reset(DataEntry<scalar>::New("TvsPSat", evapInfo.dict()).ptr());
return is;
}
Foam::Ostream& Foam::operator<<(Ostream& os, const evaporationProperties& ep)
{
os.check
(
"Foam::Ostream& Foam::operator<<"
"("
"Foam::Ostream&,"
"const Foam::evaporationProperties&"
")"
);
os << ep.name_ << nl << token::BEGIN_BLOCK << nl
<< incrIndent;
os.writeKeyword("Dab") << ep.Dab_ << token::END_STATEMENT << nl;
// os << ep.TvsPSat_() << nl;
os << decrIndent << token::END_BLOCK << nl;
os.check
(
"Foam::Ostream& Foam::operator<<"
"("
"Foam::Ostream&,"
"const Foam::evaporationProperties&"
")"
);
return os;
}
// ************************************************************************* //

View File

@ -35,18 +35,18 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
)
{
rhoc_ = td.rhoInterp().interpolate(this->position(), celli);
Uc_ = td.UInterp().interpolate(this->position(), celli);
muc_ = td.muInterp().interpolate(this->position(), celli);
rhoc_ = td.rhoInterp().interpolate(this->position(), cellI);
Uc_ = td.UInterp().interpolate(this->position(), cellI);
muc_ = td.muInterp().interpolate(this->position(), cellI);
// Apply dispersion components to carrier phase velocity
Uc_ = td.cloud().dispersion().update
(
dt,
celli,
cellI,
U_,
Uc_,
UTurb_,
@ -57,44 +57,33 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
template<class ParcelType>
template<class TrackData>
void Foam::KinematicParcel<ParcelType>::calcCoupled
void Foam::KinematicParcel<ParcelType>::calc
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Define local properties at beginning of timestep
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// const scalar mass0 = mass();
// const vector U0 = U_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U, dUTrans
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
vector dUTrans = vector::zero;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate source terms
// ~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().coupled())
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Update momentum transfer
td.cloud().UTrans()[celli] += nParticle_*dUTrans;
// Update momentum transfer
td.cloud().UTrans()[cellI] += nParticle_*dUTrans;
// Accumulate coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[celli] += nParticle_*mass()*Cud;
// Coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += nParticle_*mass()*Cud;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -106,42 +95,18 @@ void Foam::KinematicParcel<ParcelType>::calcCoupled
template<class ParcelType>
template<class TrackData>
void Foam::KinematicParcel<ParcelType>::calcUncoupled
(
TrackData& td,
const scalar dt,
const label
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
this->U() = calcVelocity(td, dt, Cud, dUTrans);
}
template<class ParcelType>
template<class TrackData>
Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
TrackData& td,
const scalar dt,
scalar& Cud,
vector& dUTrans
)
) const
{
// Return linearised term from drag model
Cud = td.cloud().drag().Cu(U_ - Uc_, d_, rhoc_, rho_, muc_);
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle velocity
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -188,7 +153,7 @@ bool Foam::KinematicParcel<ParcelType>::move
// Remember which cell the Parcel is in
// since this will change if a face is hit
label celli = p.cell();
label cellI = p.cell();
dt *= p.trackToFace(p.position() + dt*U_, td);
@ -196,19 +161,12 @@ bool Foam::KinematicParcel<ParcelType>::move
p.stepFraction() = 1.0 - tEnd/deltaT;
// Update cell based properties
p.updateCellQuantities(td, dt, celli);
p.updateCellQuantities(td, dt, cellI);
// Avoid problems with extremely small timesteps
if (dt > ROOTVSMALL)
{
if (td.cloud().coupled())
{
p.calcCoupled(td, dt, celli);
}
else
{
p.calcUncoupled(td, dt, celli);
}
p.calc(td, dt, cellI);
}
if (p.onBoundary() && td.keepParticle)

View File

@ -224,13 +224,13 @@ protected:
//- Calculate new particle velocity
template<class TrackData>
vector calcVelocity
const vector calcVelocity
(
TrackData& td,
const scalar dt,
scalar& Cud,
vector& dUTrans
);
) const;
public:
@ -249,7 +249,7 @@ public:
KinematicCloud<ParcelType>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
@ -355,25 +355,16 @@ public:
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
);
//- Coupled calculation with the continuous phase
//- Update parcel properties over the time interval
template<class TrackData>
void calcCoupled
void calc
(
TrackData& td,
const scalar dt,
const label celli
);
//- Uncoupled calculation with the continuous phase
template<class TrackData>
void calcUncoupled
(
TrackData& td,
const scalar dt,
const label
const label cellI
);

View File

@ -64,14 +64,14 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
KinematicCloud<ParcelType>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const constantProperties& constProps
)
:
Particle<ParcelType>(owner, position, celli),
Particle<ParcelType>(owner, position, cellI),
typeId_(typeId),
d_(d0),
U_(U0),

View File

@ -34,20 +34,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::updateCellQuantities
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
)
{
ReactingParcel<ParcelType>::updateCellQuantities(td, dt, celli);
ReactingParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
}
template<class ParcelType>
template<class TrackData>
void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
void Foam::ReactingMultiphaseParcel<ParcelType>::calc
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -58,7 +58,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
const scalar cp0 = this->cp_;
const scalar np0 = this->nParticle_;
const scalar T0 = this->T_;
scalarField& YMix = this->YMixture_;
scalarField& YMix = this->Y_;
label idGas = td.cloud().composition().idGas();
label idLiquid = td.cloud().composition().idLiquid();
@ -94,14 +94,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
// Calculate heat transfer - update T
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~
// Calculate phase change
// ~~~~~~~~~~~~~~~~~~~~~~
scalarField X = td.cloud().composition().X(idLiquid, YLiquid_);
calcPhaseChange(td, dt, T0, X, dMassMT);
scalar dMassPC = calcPhaseChange(td, dt, cellI, T1, dMassMT);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -116,7 +115,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
// Initialise enthalpy retention to zero
scalar dhRet = 0.0;
calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
calcSurfaceReactions(td, dt, cellI, T0, T1, dMassMT, dMassSR, dhRet);
// New total mass
const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
@ -128,10 +127,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
*sum(dMassMT)
/(0.5*(mass0 + mass1)*cp0);
// Add retained enthalpy from surface reaction to particle and remove
// from gas
T1 += dhRet/(0.5*(mass0 + mass1)*cp0);
dhTrans -= dhRet;
// Enthalpy retention divided between particle and carrier phase by the
// fraction fh and (1 - fh)
T1 += td.constProps().fh()*dhRet/(0.5*(mass0 + mass1)*cp0);
dhTrans -= (1.0 - td.constProps().fh())*dhRet;
// Correct dhTrans to account for enthalpy of evolved volatiles
dhTrans +=
@ -151,21 +150,21 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
// Transfer mass lost from particle to carrier mass source
forAll(dMassMT, i)
{
td.cloud().rhoTrans(i)[celli] += np0*(dMassMT[i] + dMassSR[i]);
td.cloud().rhoTrans(i)[cellI] += np0*(dMassMT[i] + dMassSR[i]);
}
// Update momentum transfer
td.cloud().UTrans()[celli] += np0*dUTrans;
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Accumulate coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[celli] += np0*mass0*Cud;
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
// Update enthalpy transfer
// - enthalpy of lost solids already accounted for
td.cloud().hTrans()[celli] += np0*dhTrans;
td.cloud().hTrans()[cellI] += np0*dhTrans;
// Accumulate coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -178,9 +177,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
// Absorb particle(s) into carrier phase
forAll(dMassMT, i)
{
td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
td.cloud().rhoTrans(i)[cellI] += np0*dMassMT[i];
}
td.cloud().hTrans()[celli] +=
td.cloud().hTrans()[cellI] +=
np0*mass1
*(
YMix[0]
@ -188,7 +187,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
+ YMix[2]
*td.cloud().composition().H(idSolid, YSolid_, this->pc_, T1)
);
td.cloud().UTrans()[celli] += np0*mass1*U1;
td.cloud().UTrans()[cellI] += np0*mass1*U1;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
@ -215,123 +214,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
}
template<class ParcelType>
template<class TrackData>
void Foam::ReactingMultiphaseParcel<ParcelType>::calcUncoupled
(
TrackData& td,
const scalar dt,
const label celli
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Define local properties at beginning of timestep
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar T0 = this->T_;
const scalar mass0 = this->mass();
const scalar cp0 = this->cp_;
scalarField& YMix = this->YMixture_;
label idGas = td.cloud().composition().idGas();
label idLiquid = td.cloud().composition().idLiquid();
label idSolid = td.cloud().composition().idSolid();
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// Enthalpy transfer from the particle to the carrier phase
scalar dhTrans = 0.0;
// Mass transfer from particle to carrier phase
// - components exist in particle already
scalarList dMassMT(td.cloud().gases().size(), 0.0);
// 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);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate heat transfer - update T
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~
// Calculate phase change
// ~~~~~~~~~~~~~~~~~~~~~~
scalarField X = td.cloud().composition().X(idLiquid, YLiquid_);
calcPhaseChange(td, dt, T0, YLiquid_, dMassMT);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate Devolatilisation
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
calcDevolatilisation(td, dt, T0, T1, dMassMT);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate surface reactions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 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) - sum(dMassSR);
// New specific heat capacity
const scalar cp1 =
YMix[0]*td.cloud().composition().cp(idGas, YGas_, this->pc_, T1)
+ YMix[1]*td.cloud().composition().cp(idLiquid, YLiquid_, this->pc_, T1)
+ YMix[2]*td.cloud().composition().cp(idSolid, YSolid_, this->pc_, T1);
// Add retained enthalpy to particle
T1 += dhRet/(mass0*0.5*(cp0 + cp1));
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mass1 < td.constProps().minParticleMass())
{
td.keepParticle = false;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
else
{
this->U_ = U1;
this->T_ = T1;
this->cp_ = cp1;
// Update particle density or diameter
if (td.constProps().constantVolume())
{
this->rho_ = mass1/this->volume();
}
else
{
this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
}
}
}
template<class ParcelType>
template<class TrackData>
void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
@ -358,7 +240,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
")\n"
"no treatment currently available for particles containing "
"liquid species"
)
);
}
// Check that model is active, and that the parcel temperature is
@ -375,7 +257,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
// Determine mass to add to carrier phase
const scalar mass = this->mass();
scalarField& YMix = this->YMixture_;
scalarField& YMix = this->Y_;
const scalar dMassTot = td.cloud().devolatilisation().calculate
(
dt,
@ -413,7 +295,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
(
TrackData& td,
const scalar dt,
const label celli,
const label cellI,
const scalar T0,
const scalar T1,
const scalarList& dMassMT,
@ -432,7 +314,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
td.cloud().surfaceReaction().calculate
(
dt,
celli,
cellI,
this->d_,
T0,
T1,
@ -444,12 +326,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
YGas_,
YLiquid_,
YSolid_,
this->YMixture_,
this->Y_,
dMassSR,
dhRet
);
// TODO: td.cloud().addToMassDevolatilisation(sum(dMassSR));
// TODO: td.cloud().addToMassSurfaceReaction(sum(dMassSR));
}

View File

@ -79,6 +79,10 @@ public:
//- Latent heat of devolatilisation [J/kg]
const scalar Ldevol_;
//- Fraction of enthalpy retained by parcel due to surface
// reactions
const scalar fh_;
public:
@ -89,6 +93,10 @@ public:
//- Return const access to the latent heat of devolatilisation
inline scalar Ldevol() const;
//- Return const access to the fraction of enthalpy retained by
// parcel due to surface reactions
inline scalar fh() const;
};
@ -174,7 +182,7 @@ protected:
(
TrackData& td,
const scalar dt,
const label celli,
const label cellI,
const scalar T0,
const scalar T1,
const scalarList& dMassMT,
@ -199,14 +207,14 @@ public:
ReactingMultiphaseCloud<ParcelType>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const scalarField& YGas0,
const scalarField& YLiquid0,
const scalarField& YSolid0,
const scalarField& YMixture0,
const scalarField& Y0,
const constantProperties& constProps
);
@ -262,25 +270,16 @@ public:
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
);
//- Coupled calculation with the continuous phase
//- Update parcel properties over the time interval
template<class TrackData>
void calcCoupled
void calc
(
TrackData& td,
const scalar dt,
const label celli
);
//- Uncoupled calculation with the continuous phase
template<class TrackData>
void calcUncoupled
(
TrackData& td,
const scalar dt,
const label
const label cellI
);

View File

@ -34,7 +34,8 @@ constantProperties
)
:
ReactingParcel<ParcelType>::constantProperties(dict),
Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value())
Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value()),
fh_(dimensionedScalar(dict.lookup("fh")).value())
{}
@ -75,14 +76,14 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
ReactingMultiphaseCloud<ParcelType>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const scalarField& YGas0,
const scalarField& YLiquid0,
const scalarField& YSolid0,
const scalarField& YMixture0,
const scalarField& Y0,
const constantProperties& constProps
)
:
@ -91,11 +92,11 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
owner,
typeId,
position,
celli,
cellI,
d0,
U0,
nParticle0,
YMixture0,
Y0,
constProps
),
YGas_(YGas0),
@ -114,6 +115,14 @@ Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::Ldevol() const
}
template<class ParcelType>
inline Foam::scalar
Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::fh() const
{
return fh_;
}
// * * * * * * * * * * * trackData Member Functions * * * * * * * * * * * * //
template<class ParcelType>

View File

@ -68,7 +68,7 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
}
// scale the mass fractions
const scalarField& YMix = this->YMixture_;
const scalarField& YMix = this->Y_;
YGas_ /= YMix[0] + ROOTVSMALL;
YLiquid_ /= YMix[1] + ROOTVSMALL;
YSolid_ /= YMix[2] + ROOTVSMALL;
@ -129,7 +129,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingMultiphaseParcel<ParcelType>& p = iter();
p.YGas_[j] = YGas[i++]/(p.YMixture()[0] + ROOTVSMALL);
p.YGas_[j] = YGas[i++]/(p.Y()[0] + ROOTVSMALL);
}
}
// Populate YLiquid for each parcel
@ -144,7 +144,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingMultiphaseParcel<ParcelType>& p = iter();
p.YLiquid_[j] = YLiquid[i++]/(p.YMixture()[1] + ROOTVSMALL);
p.YLiquid_[j] = YLiquid[i++]/(p.Y()[1] + ROOTVSMALL);
}
}
// Populate YSolid for each parcel
@ -159,7 +159,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingMultiphaseParcel<ParcelType>& p = iter();
p.YSolid_[j] = YSolid[i++]/(p.YMixture()[2] + ROOTVSMALL);
p.YSolid_[j] = YSolid[i++]/(p.Y()[2] + ROOTVSMALL);
}
}
}
@ -192,7 +192,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
forAllConstIter(typename Cloud<ParcelType>, c, iter)
{
const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
YGas[i++] = p0.YGas()[j]*p0.YMixture()[0];
YGas[i++] = p0.YGas()[j]*p0.Y()[0];
}
YGas.write();
@ -211,7 +211,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
forAllConstIter(typename Cloud<ParcelType>, c, iter)
{
const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
YLiquid[i++] = p0.YLiquid()[j]*p0.YMixture()[1];
YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[1];
}
YLiquid.write();
@ -230,7 +230,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
forAllConstIter(typename Cloud<ParcelType>, c, iter)
{
const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
YSolid[i++] = p0.YSolid()[j]*p0.YMixture()[2];
YSolid[i++] = p0.YSolid()[j]*p0.Y()[2];
}
YSolid.write();
@ -248,9 +248,9 @@ Foam::Ostream& Foam::operator<<
const ReactingMultiphaseParcel<ParcelType>& p
)
{
scalarField YGasLoc = p.YGas()*p.YMixture()[0];
scalarField YLiquidLoc = p.YLiquid()*p.YMixture()[1];
scalarField YSolidLoc = p.YSolid()*p.YMixture()[2];
scalarField YGasLoc = p.YGas()*p.Y()[0];
scalarField YLiquidLoc = p.YLiquid()*p.Y()[1];
scalarField YSolidLoc = p.YSolid()*p.Y()[2];
if (os.format() == IOstream::ASCII)
{
os << static_cast<const ReactingMultiphaseParcel<ParcelType>& >(p)

View File

@ -34,22 +34,38 @@ void Foam::ReactingParcel<ParcelType>::updateCellQuantities
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
)
{
ThermoParcel<ParcelType>::updateCellQuantities(td, dt, celli);
ThermoParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
pc_ = td.pInterp().interpolate(this->position(), celli);
pc_ = td.pInterp().interpolate(this->position(), cellI);
}
template<class ParcelType>
void Foam::ReactingParcel<ParcelType>::updateMassFraction
(
const scalar mass0,
const scalar mass1,
const scalarList& dMass,
scalarField& Y
)
{
forAll(Y, i)
{
Y[i] = (Y[i]*mass0 - dMass[i])/mass1;
}
}
template<class ParcelType>
template<class TrackData>
void Foam::ReactingParcel<ParcelType>::calcCoupled
void Foam::ReactingParcel<ParcelType>::calc
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -61,7 +77,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
const scalar np0 = this->nParticle_;
const scalar T0 = this->T_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
@ -73,7 +88,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
scalar dhTrans = 0.0;
// Mass transfer from particle to carrier phase
// - components exist in particle already
scalarList dMassMT(td.cloud().gases().size(), 0.0);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -87,147 +101,58 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
// Calculate heat transfer - update T
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~
// Calculate phase change
// ~~~~~~~~~~~~~~~~~~~~~~
scalarField X = td.cliud().composition().X(0, YMixture_);
calcPhaseChange(td, dt, T, X, dMassMT);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate phase change - update mass, Y, cp, T, dhTrans
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar dMassPC = calcPhaseChange(td, dt, cellI, T1, dMassMT);
// New total mass
const scalar mass1 = mass0 - sum(dMassMT);
// Update particle mass
scalar mass1 = mass0 - dMassPC;
// Update particle component mass fractions
updateMassFraction(mass0, mass1, dMassMT, Y_);
// New specific heat capacity of
scalar cp1 = td.cloud().composition().cp(0, Y_, pc_, T1);
// Correct temperature due to evaporated components
// TODO: use hl function in liquidMixture???
T1 -= td.constProps().Lvap()*dMassPC/(0.5*(mass0 + mass1)*cp1);
// Correct dhTrans to account for the change in enthalpy due to the
// phase change
scalar HMix = td.cloud().composition().H(0, Y_, pc, 0.5*(T0 + T1));
dhTrans += dMassPC*HMix;
// ~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate source terms
// ~~~~~~~~~~~~~~~~~~~~~~~
// Transfer mass lost from particle to carrier mass source
forAll(dMassMT, i)
if (td.cloud().coupled())
{
td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Update momentum transfer
td.cloud().UTrans()[celli] += np0*dUTrans;
// Accumulate coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[celli] += np0*mass0*Cud;
// Update enthalpy transfer
td.cloud().hTrans()[celli] += np0*dhTrans;
// Accumulate coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mass1 < td.constProps().minParticleMass())
{
td.keepParticle = false;
// Absorb particle(s) into carrier phase
// Transfer mass lost from particle to carrier mass source
forAll(dMassMT, i)
{
td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
td.cloud().rhoTrans(i)[cellI] += np0*dMassMT[i];
}
td.cloud().UTrans()[celli] += np0*mass1*U1;
// 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;
// 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();
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
else
{
this->U_ = U1;
this->T_ = T1;
// this->cp_ = ??? // TODO:
// Update particle density or diameter
if (td.constProps().constantVolume())
{
this->rho_ = mass1/this->volume();
}
else
{
this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
}
}
}
template<class ParcelType>
template<class TrackData>
void Foam::ReactingParcel<ParcelType>::calcUncoupled
(
TrackData& td,
const scalar dt,
const label celli
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Define local properties at beginning of timestep
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar T0 = this->T_;
const scalar mass0 = this->mass();
const scalar cp0 = this->cp_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// Enthalpy transfer from the particle to the carrier phase
scalar dhTrans = 0.0;
// Mass transfer from particle to carrier phase
// - components exist in particle already
scalarList dMassMT(td.cloud().gases().size(), 0.0);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate heat transfer - update T
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~
// Calculate phase change
// ~~~~~~~~~~~~~~~~~~~~~~
scalarField X = td.cloud().composition().X(0, YMixture_);
scalar dMassPC = calcPhaseChange(td, dt, T, X, dMassMT);
T1 -= td.constProps().Lvap()*dMassPC/(0.5*mass0*cp0);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate surface reactions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise enthalpy retention to zero
scalar dhRet = 0.0;
// New total mass
const scalar mass1 = mass0 - sum(dMassMT);
// New specific heat capacity
const scalar cp1 = cp0; // TODO: new cp1
// Add retained enthalpy to particle
T1 += dhRet/(mass0*0.5*(cp0 + cp1));
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Remove the particle when mass falls below minimum threshold
@ -235,6 +160,17 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
if (mass1 < td.constProps().minParticleMass())
{
td.keepParticle = false;
if (td.cloud().coupled())
{
// Absorb particle(s) into carrier phase
forAll(dMassMT, i)
{
td.cloud().rhoTrans(i)[cellI] += np0*mass1*Y_[i];
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
td.cloud().hTrans()[cellI] += np0*mass1*HMix;
}
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
@ -260,39 +196,36 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
template<class ParcelType>
template<class TrackData>
void Foam::ReactingParcel<ParcelType>::calcPhaseChange
Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
(
TrackData& td,
const scalar dt,
const label cellI,
const scalar T,
scalarField& X,
scalarList& dMassMT
)
{
if (!td.cloud().phaseChange().active())
{
return;
return 0.0;
}
// TODO: separate treatment for boiling
scalar dMassTot = td.cloud().phaseChange().calculate
scalar dMass = td.cloud().phaseChange().calculate
(
dt,
cellI,
T,
this->d_,
X,
dMassMT,
this->U_ - this->Uc_,
this->Tc_,
pc_,
this->d_,
this->Tc_,
this->muc_/this->rhoc_,
dt
this->U_ - this->Uc_,
dMassMT
);
td.cloud().addToMassPhaseChange(dMassTot);
// TODO: Re-calculate mass fractions
td.cloud().addToMassPhaseChange(dMass);
return dMass;
}

View File

@ -175,7 +175,7 @@ protected:
scalar mass0_;
//- Mass fractions of mixture []
scalarField YMixture_;
scalarField Y_;
// Cell-based quantities
@ -188,15 +188,24 @@ protected:
//- Calculate Phase change
template<class TrackData>
void calcPhaseChange
scalar calcPhaseChange
(
TrackData& td,
const scalar dt,
const label cellI,
const scalar T,
scalarField& X,
scalarList& dMassMT
);
//- Update mass fraction
void updateMassFraction
(
const scalar mass0,
const scalar mass1,
const scalarList& dMass,
scalarField& Y
);
public:
@ -214,11 +223,11 @@ public:
ReactingCloud<ParcelType>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const scalarField& YMixture0,
const scalarField& Y0,
const constantProperties& constProps
);
@ -245,7 +254,7 @@ public:
inline scalar mass0() const;
//- Return const access to mass fractions of mixture
inline const scalarField& YMixture() const;
inline const scalarField& Y() const;
//- Return the owner cell pressure
inline scalar pc() const;
@ -257,7 +266,7 @@ public:
inline scalar& mass0();
//- Return access to mass fractions of mixture
inline scalarField& YMixture();
inline scalarField& Y();
// Main calculation loop
@ -268,25 +277,16 @@ public:
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
);
//- Coupled calculation with the continuous phase
//- Update parcel properties over the time interval
template<class TrackData>
void calcCoupled
void calc
(
TrackData& td,
const scalar dt,
const label celli
);
//- Uncoupled calculation with the continuous phase
template<class TrackData>
void calcUncoupled
(
TrackData& td,
const scalar dt,
const label
const label cellI
);

View File

@ -77,11 +77,11 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
ReactingCloud<ParcelType>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const scalarField& YMixture0,
const scalarField& Y0,
const constantProperties& constProps
)
:
@ -90,14 +90,14 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
owner,
typeId,
position,
celli,
cellI,
d0,
U0,
nParticle0,
constProps
),
mass0_(0.0),
YMixture_(YMixture0),
Y_(Y0),
pc_(0.0)
{
// Set initial parcel mass
@ -176,9 +176,9 @@ inline Foam::scalar Foam::ReactingParcel<ParcelType>::mass0() const
template<class ParcelType>
inline const Foam::scalarField& Foam::ReactingParcel<ParcelType>::
YMixture() const
Y() const
{
return YMixture_;
return Y_;
}
@ -197,9 +197,9 @@ inline Foam::scalar& Foam::ReactingParcel<ParcelType>::mass0()
template<class ParcelType>
inline Foam::scalarField& Foam::ReactingParcel<ParcelType>::YMixture()
inline Foam::scalarField& Foam::ReactingParcel<ParcelType>::Y()
{
return YMixture_;
return Y_;
}

View File

@ -39,7 +39,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
:
ThermoParcel<ParcelType>(cloud, is, readFields),
mass0_(0.0),
YMixture_(0),
Y_(0),
pc_(0.0)
{
if (readFields)
@ -48,11 +48,11 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
dynamic_cast<const ReactingCloud<ParcelType>& >(cloud);
const label nMixture = cR.composition().phaseTypes().size();
YMixture_.setSize(nMixture);
Y_.setSize(nMixture);
if (is.format() == IOstream::ASCII)
{
is >> mass0_ >> YMixture_;
is >> mass0_ >> Y_;
}
else
{
@ -61,7 +61,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
reinterpret_cast<char*>(&mass0_),
+ sizeof(mass0_)
);
is >> YMixture_;
is >> Y_;
}
}
@ -109,13 +109,13 @@ void Foam::ReactingParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingParcel<ParcelType>& p = iter();
p.YMixture_.setSize(nPhases, 0.0);
p.Y_.setSize(nPhases, 0.0);
}
// Populate YMixture for each parcel
// Populate Y for each parcel
forAll(phaseTypes, j)
{
IOField<scalar> YMixture
IOField<scalar> Y
(
c.fieldIOobject("Y" + phaseTypes[j], IOobject::MUST_READ)
);
@ -124,7 +124,7 @@ void Foam::ReactingParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingParcel<ParcelType>& p = iter();
p.YMixture_[j] = YMixture[i++];
p.Y_[j] = Y[i++];
}
}
}
@ -157,7 +157,7 @@ void Foam::ReactingParcel<ParcelType>::writeFields
const wordList phaseTypes = c.composition().phaseTypes();
forAll(phaseTypes, j)
{
IOField<scalar> YMixture
IOField<scalar> Y
(
c.fieldIOobject("Y" + phaseTypes[j], IOobject::NO_READ),
np
@ -167,10 +167,10 @@ void Foam::ReactingParcel<ParcelType>::writeFields
forAllConstIter(typename Cloud<ParcelType>, c, iter)
{
const ReactingParcel<ParcelType>& p0 = iter();
YMixture[i++] = p0.YMixture()[j];
Y[i++] = p0.Y()[j];
}
YMixture.write();
Y.write();
}
}
}
@ -189,7 +189,7 @@ Foam::Ostream& Foam::operator<<
{
os << static_cast<const ThermoParcel<ParcelType>& >(p)
<< token::SPACE << p.mass0()
<< token::SPACE << p.YMixture();
<< token::SPACE << p.Y();
}
else
{
@ -200,7 +200,7 @@ Foam::Ostream& Foam::operator<<
reinterpret_cast<const char*>(&p.mass0_),
sizeof(p.mass0())
);
os << p.YMixture();
os << p.Y();
}
// Check state of Ostream

View File

@ -34,23 +34,23 @@ void Foam::ThermoParcel<ParcelType>::updateCellQuantities
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
)
{
KinematicParcel<ParcelType>::updateCellQuantities(td, dt, celli);
KinematicParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
Tc_ = td.TInterp().interpolate(this->position(), celli);
cpc_ = td.cpInterp().interpolate(this->position(), celli);
Tc_ = td.TInterp().interpolate(this->position(), cellI);
cpc_ = td.cpInterp().interpolate(this->position(), cellI);
}
template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::calcCoupled
void Foam::ThermoParcel<ParcelType>::calc
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -59,8 +59,6 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
const vector U0 = this->U_;
const scalar mass0 = this->mass();
const scalar np0 = this->nParticle_;
// const scalar T0 = T_;
// const scalar cp0 = cp_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~
@ -85,25 +83,27 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
// Calculate heat transfer - update T
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
const scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
const scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate source terms
// ~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().coupled())
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Update momentum transfer
td.cloud().UTrans()[celli] += np0*dUTrans;
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Accumulate coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[celli] += np0*mass0*Cud;
// Coefficient to be applied in carrier phase momentum coupling
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
// Update enthalpy transfer
td.cloud().hTrans()[celli] += np0*dhTrans;
// Accumulate coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
// 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();
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle properties
@ -113,48 +113,13 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
}
template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::calcUncoupled
(
TrackData& td,
const scalar dt,
const label celli
)
{
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise transfer terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// Enthalpy transfer from the particle to the carrier phase
scalar dhTrans = 0.0;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
this->U_ = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate heat transfer - update T
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
T_ = calcHeatTransfer(td, dt, celli, htc, dhTrans);
}
template<class ParcelType>
template <class TrackData>
Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
(
TrackData& td,
const scalar dt,
const label celli,
const label cellI,
scalar& htc,
scalar& dhTrans
)
@ -193,7 +158,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
const scalar sigma = radiation::sigmaSB.value();
const scalar epsilon = td.constProps().epsilon0();
const scalar epsilonSigmaT3 = epsilon*sigma*pow3(T_);
ap = (htc*Tc_ + 0.25*epsilon*G[celli])/(htc + epsilonSigmaT3);
ap = (htc*Tc_ + 0.25*epsilon*G[cellI])/(htc + epsilonSigmaT3);
bp += epsilonSigmaT3;
}
bp *= 6.0/(this->rho_*this->d_*cp_);

View File

@ -209,7 +209,7 @@ protected:
(
TrackData& td,
const scalar dt,
const label celli,
const label cellI,
scalar& htc,
scalar& dhTrans
);
@ -231,7 +231,7 @@ public:
ThermoCloud<ParcelType>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
@ -281,25 +281,16 @@ public:
(
TrackData& td,
const scalar dt,
const label celli
const label cellI
);
//- Coupled calculation with the continuous phase
//- Update parcel properties over the time interval
template<class TrackData>
void calcCoupled
void calc
(
TrackData& td,
const scalar dt,
const label celli
);
//- Uncoupled calculation with the continuous phase
template<class TrackData>
void calcUncoupled
(
TrackData& td,
const scalar dt,
const label
const label cellI
);

View File

@ -75,7 +75,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
ThermoCloud<ParcelType>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
@ -87,7 +87,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
owner,
typeId,
position,
celli,
cellI,
d0,
U0,
nParticle0,

View File

@ -43,7 +43,7 @@ Foam::basicKinematicParcel::basicKinematicParcel
KinematicCloud<basicKinematicParcel>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
@ -55,7 +55,7 @@ Foam::basicKinematicParcel::basicKinematicParcel
owner,
typeId,
position,
celli,
cellI,
d0,
U0,
nParticle0,

View File

@ -66,7 +66,7 @@ public:
KinematicCloud<basicKinematicParcel>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,

View File

@ -43,14 +43,14 @@ Foam::basicReactingMultiphaseParcel::basicReactingMultiphaseParcel
ReactingMultiphaseCloud<basicReactingMultiphaseParcel>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const scalarField& YGas0,
const scalarField& YLiquid0,
const scalarField& YSolid0,
const scalarField& YMixture0,
const scalarField& Y0,
const constantProperties& constProps
)
:
@ -59,14 +59,14 @@ Foam::basicReactingMultiphaseParcel::basicReactingMultiphaseParcel
owner,
typeId,
position,
celli,
cellI,
d0,
U0,
nParticle0,
YGas0,
YLiquid0,
YSolid0,
YMixture0,
Y0,
constProps
)
{}

View File

@ -66,14 +66,14 @@ public:
ReactingMultiphaseCloud<basicReactingMultiphaseParcel>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const scalarField& YGas0,
const scalarField& YLiquid0,
const scalarField& YSolid0,
const scalarField& YMixture0,
const scalarField& Y0,
const constantProperties& constProps
);

View File

@ -43,11 +43,11 @@ Foam::basicReactingParcel::basicReactingParcel
ReactingCloud<basicReactingParcel>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const scalarField& YMixture0,
const scalarField& Y0,
const constantProperties& constProps
)
:
@ -56,11 +56,11 @@ Foam::basicReactingParcel::basicReactingParcel
owner,
typeId,
position,
celli,
cellI,
d0,
U0,
nParticle0,
YMixture0,
Y0,
constProps
)
{}

View File

@ -66,11 +66,11 @@ public:
ReactingCloud<basicReactingParcel>& owner,
const label typeId,
const vector& position,
const label celli,
const label cellI,
const scalar d0,
const vector& U0,
const scalar nParticle0,
const scalarField& YMixture0,
const scalarField& Y0,
const constantProperties& constProps
);

View File

@ -43,7 +43,7 @@ Foam::basicThermoParcel::basicThermoParcel
ThermoCloud<basicThermoParcel>& owner,
const label typeId,
const vector position,
const label celli,
const label cellI,
const scalar d0,
const vector U0,
const scalar nParticle0,
@ -55,7 +55,7 @@ Foam::basicThermoParcel::basicThermoParcel
owner,
typeId,
position,
celli,
cellI,
d0,
U0,
nParticle0,

View File

@ -65,7 +65,7 @@ public:
ThermoCloud<basicThermoParcel>& owner,
const label typeId,
const vector position,
const label celli,
const label cellI,
const scalar d0,
const vector U0,
const scalar nParticle0,

View File

@ -26,6 +26,13 @@ License
#include "polynomial.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(polynomial, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polynomial::polynomial(const word& entryName, Istream& is)

View File

@ -67,6 +67,26 @@ Foam::label Foam::LiquidEvaporation<CloudType>::carrierSpecieId
}
template<class CloudType>
Foam::scalarField Foam::LiquidEvaporation<CloudType>::calcXc
(
const label cellI
) const
{
scalarField Xc(this->owner().carrierThermo().composition().Y().size());
scalar Winv = 0.0;
forAll(Xc, i)
{
scalar Y = this->owner().carrierThermo().composition().Y()[i][cellI];
Winv += Y/this->owner().gases()[i].W();
Xc[i] = Y/this->owner().gases()[i].W();
}
return Xc/Winv;
}
template <class CloudType>
Foam::scalar Foam::LiquidEvaporation<CloudType>::Sh
(
@ -99,10 +119,10 @@ Foam::LiquidEvaporation<CloudType>::LiquidEvaporation
)
),
Tvap_(readScalar(this->coeffDict().lookup("Tvap"))),
evapProps_(this->coeffDict().lookup("activeLiquids")),
liqToGasMap_(evapProps_.size())
activeLiquids_(this->coeffDict().lookup("activeLiquids")),
liqToGasMap_(activeLiquids_.size())
{
if (evapProps_.size() == 0)
if (activeLiquids_.size() == 0)
{
WarningIn
(
@ -116,9 +136,9 @@ Foam::LiquidEvaporation<CloudType>::LiquidEvaporation
}
// Calculate mapping between liquid and carrier phase species
forAll(evapProps_, i)
forAll(activeLiquids_, i)
{
liqToGasMap_[i] = carrierSpecieId(evapProps_[i].name());
liqToGasMap_[i] = carrierSpecieId(activeLiquids_[i]);
}
}
@ -142,15 +162,15 @@ bool Foam::LiquidEvaporation<CloudType>::active() const
template<class CloudType>
Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
(
const scalar dt,
const label cellI,
const scalar T,
const scalar d,
const scalarField& Xc,
scalarList& dMassMT,
const vector& Ur,
const scalar Tc,
const scalar pc,
const scalar d,
const scalar Tc,
const scalar nuc,
const scalar dt
const vector& Ur,
scalarList& dMassMT
) const
{
scalar dMassTot = 0.0;
@ -162,20 +182,25 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
}
else
{
// droplet area
// Construct carrier phase species volume fractions
scalarField Xc = calcXc(cellI);
// droplet surface area
scalar A = mathematicalConstant::pi*sqr(d);
// Reynolds number
scalar Re = mag(Ur)*d/(nuc + ROOTVSMALL);
// Calculate mass transfer of each specie in liquid
forAll(evapProps_, i)
forAll(activeLiquids_, i)
{
// Diffusion coefficient for species i
scalar Dab = evapProps_[i].Dab();
label gid = liqToGasMap_[i];
// Saturation pressure for species i at temperature T
scalar pSat = evapProps_[i].TvsPSat().value(T);
// Vapour diffusivity [m2/s]
scalar Dab = liquids_->properties()[gid].D(pc, T);
// Saturation pressure for species i [pa]
scalar pSat = liquids_->properties()[gid].pv(pc, T);
// Schmidt number
scalar Sc = nuc/(Dab + ROOTVSMALL);
@ -190,15 +215,14 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
scalar Cs = pSat/(specie::RR*T);
// vapour concentration in bulk gas [kgmol/m3]
scalar Cinf = Xc[i]*pc/(specie::RR*Tc);
scalar Cinf = Xc[gid]*pc/(specie::RR*Tc);
// molar flux of vapour [kgmol/m2/s]
scalar Ni = max(kc*(Cs - Cinf), 0.0);
// mass transfer
label globalLiqId = liqToGasMap_[i];
scalar dm = Ni*A*liquids_->properties()[globalLiqId].W()*dt;
dMassMT[globalLiqId] -= dm;
// mass transfer [kg]
scalar dm = Ni*A*liquids_->properties()[gid].W()*dt;
dMassMT[gid] -= dm;
dMassTot += dm;
}
}

View File

@ -35,9 +35,6 @@ Description
#include "PhaseChangeModel.H"
#include "liquidMixture.H"
#include "Tuple2.H"
#include "evaporationProperties.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -64,8 +61,8 @@ protected:
// allowed
scalar Tvap_;
//- List of evaporation properties
List<evaporationProperties> evapProps_;
//- List of active liquid names
List<word> activeLiquids_;
//- Mapping between liquid and carrier gas species
List<label> liqToGasMap_;
@ -79,6 +76,9 @@ protected:
//- Sherwood number as a function of Reynolds and Schmidt numbers
scalar Sh(const scalar Re, const scalar Sc) const;
//- Calculate the carrier phase comonent volume fractions at cellI
scalarField calcXc(const label cellI) const;
public:
@ -108,15 +108,15 @@ public:
//- Update model
scalar calculate
(
const scalar dt,
const label cellI,
const scalar T,
const scalar d,
const scalarField& Xc,
scalarList& dMassMT,
const vector& Ur,
const scalar Tc,
const scalar pc,
const scalar d,
const scalar Tc,
const scalar nuc,
const scalar dt
const vector& Ur,
scalarList& dMassMT
) const;
};

View File

@ -58,15 +58,15 @@ bool Foam::NoPhaseChange<CloudType>::active() const
template<class CloudType>
Foam::scalar Foam::NoPhaseChange<CloudType>::calculate
(
const scalar,
const label,
const scalar,
const scalar,
const scalar,
const scalar,
const scalar,
const scalarField&,
scalarList&,
const vector&,
const scalar,
const scalar,
const scalar,
const scalar
scalarList&
) const
{
// Nothing to do...

View File

@ -72,15 +72,15 @@ public:
//- Update model
scalar calculate
(
const scalar dt,
const label cellI,
const scalar T,
const scalar d,
const scalarField& Xc,
scalarList& dMassMT,
const vector& Ur,
const scalar Tc,
const scalar pc,
const scalar d,
const scalar Tc,
const scalar nuc,
const scalar dt
const vector& Ur,
scalarList& dMassMT
) const;
};

View File

@ -138,15 +138,15 @@ public:
//- Update model
virtual scalar calculate
(
const scalar dt,
const label cellI,
const scalar T,
const scalar d,
const scalarField& Xc,
scalarList& dMassMT,
const vector& Ur,
const scalar Tc,
const scalar pc,
const scalar d,
const scalar Tc,
const scalar nuc,
const scalar dt
const vector& Ur,
scalarList& dMassMT
) const = 0;
};