multiphaseEulerFoam: Mass/heat transfer consistency and linearisation

All heat transfers that result from mass-transfer are now implemented in
terms of sensible enthalpy, so that they are consistent regardless of
which form of energy is being solved for. This has removed some spurious
temperature anomalies from a number of cases involving mass-transfer.

All heat transfers that result from mass-transfer are now linearised. In
the case of multi-specie systems this requires the specification of a
residual mass fraction, which is given by a new "residualY" keyword in
the constant/phaseProperties dictionary. If this entry is omitted for
multi-specie systems then linearisation is deactivated.

**** Details for developers ****

Methods have been added to the base heat transfer phase systems to
permit energy transfer as a result of phase change, without coupling to
a diffusive heat transfer model. These functions require a "weight" to
be specified in the call to define how the latent heat is divided
between either side of the interface. A weight of 0 indicates that the
latent heat is dissipated entirely in the upwind phase, and 1 means it
is entirely in the downwind phase.

The forms of latent heat calculation and transfer have been standardised
between the various phase systems. There are now two methods of
calculating the latent heat, and two methods of applying the transfer
(see below for details). These options are currently hard-coded into the
systems that use them, but they could be made user modifiable
per-mass-transfer in future.

Interface temperatures are now stored by the derived phase systems
alongside their corresponding mass transfer rates. These temperatures
are passed by argument to the phase-change heat transfer methods
provided by the base heat transfer systems. This allows multiple
mechanisms of mass transfer each involving different interface state to
occur across the same interface.

These changes have allowed all phase systems to use the same set of
base energy-transfer functionality.

**** Even more details for developers ****

The two forms of latent heat scheme available are:

    symmetric: The latent heat is calculated as the difference between
               the interface enthalpies on either side of an interface.
               This is the simplest form.

       upwind: The latent heat is calculated as the difference between
               the bulk enthalpy on the side of the interface that mass
               is being transferred from and the interface enthalpy on
               the side of the interface that mass is transferring to.
               This form may confer some stability benefits.

The two format of latent heat transfer are:

         heat: The latent heat is applied by transferring heat unequally
               on either side of an interface using the difference
               between the bulk phase temperatures and the interface
               temperature. No explicit latent heat source is required.
               This method has a stability advantage over the "mass"
               option, but the transfer is not energy conservative
               unless the interface temperature is exactly correct.

         mass: The latent heat is applied as an explicit mass transfer
               source to both sides of an interface. The ratio between
               the heat transfer coefficients on either side determines
               what proportion of the latent heat source ends up in each
               phase. Heat transfer is calculated equally on both sides
               of an interface using bulk phase temperatures and is not
               coupled to the thermal effect of phase change. This
               method has the advantage of being energy conservative
               even if the interface temperature is not exact, but it is
               less stable than the "heat" option at extreme conditions.
This commit is contained in:
Will Bainbridge
2020-03-20 09:49:18 +00:00
parent 1d91ec3a6b
commit 0efc492a77
30 changed files with 1891 additions and 1080 deletions

View File

@ -97,6 +97,16 @@ Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModel::dY
}
Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModel::dYfPrime
(
const word& speciesName,
const volScalarField& Tf
) const
{
return YfPrime(speciesName, Tf);
}
Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModel::D
(
const word& speciesName
@ -117,53 +127,4 @@ Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModel::D
}
Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModel::L
(
const word& speciesName,
const volScalarField& Tf
) const
{
const label speciei = composition().species()[speciesName];
const volScalarField& p(thermo_.p());
volScalarField Ha(composition().Ha(speciei, p, Tf));
const volScalarField& otherP(otherThermo_.p());
tmp<volScalarField> otherHa(nullptr);
if (otherHasComposition())
{
const label otherSpeciei = otherComposition().species()[speciesName];
otherHa = otherComposition().Ha(otherSpeciei, otherP, Tf);
}
else
{
otherHa = otherThermo_.ha(otherP, Tf);
}
return
volScalarField::New
(
IOobject::groupName("L", pair_.name()),
otherHa - Ha
);
}
void Foam::interfaceCompositionModel::addDmdtL
(
const volScalarField& K,
const volScalarField& Tf,
volScalarField& dmdtL,
volScalarField& dmdtLPrime
) const
{
forAllConstIter(hashedWordList, species_, iter)
{
const volScalarField rhoKDL(thermo_.rho()*K*D(*iter)*L(*iter, Tf));
dmdtL += rhoKDL*dY(*iter, Tf);
dmdtLPrime += rhoKDL*YfPrime(*iter, Tf);
}
}
// ************************************************************************* //

View File

@ -168,34 +168,26 @@ public:
) const = 0;
//- Mass fraction difference between the interface and the field
virtual tmp<volScalarField> dY
tmp<volScalarField> dY
(
const word& speciesName,
const volScalarField& Tf
) const;
//- Mass fraction difference between the interface and the field
// derivative w.r.t. temperature
tmp<volScalarField> dYfPrime
(
const word& speciesName,
const volScalarField& Tf
) const;
//- Mass diffusivity
virtual tmp<volScalarField> D
tmp<volScalarField> D
(
const word& speciesName
) const;
//- Latent heat
virtual tmp<volScalarField> L
(
const word& speciesName,
const volScalarField& Tf
) const;
//- Add latent heat flow rate to total
virtual void addDmdtL
(
const volScalarField& K,
const volScalarField& Tf,
volScalarField& dmdtL,
volScalarField& dmdtLPrime
) const;
//- Update the composition
virtual void update(const volScalarField& Tf) = 0;

View File

@ -44,8 +44,7 @@ alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF),
fixedDmdtf_(0),
L_(0)
fixedDmdtf_(0)
{}
@ -58,8 +57,7 @@ alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF, dict),
fixedDmdtf_(dict.lookupOrDefault<scalar>("fixedDmdtf", 0)),
L_(dict.lookupOrDefault<scalar>("L", 0))
fixedDmdtf_(dict.lookupOrDefault<scalar>("fixedDmdtf", 0))
{}
@ -79,8 +77,7 @@ alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
iF,
mapper
),
fixedDmdtf_(psf.fixedDmdtf_),
L_(psf.L_)
fixedDmdtf_(psf.fixedDmdtf_)
{}
@ -91,8 +88,7 @@ alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(psf),
fixedDmdtf_(psf.fixedDmdtf_),
L_(psf.L_)
fixedDmdtf_(psf.fixedDmdtf_)
{}
@ -104,8 +100,7 @@ alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(psf, iF),
fixedDmdtf_(psf.fixedDmdtf_),
L_(psf.L_)
fixedDmdtf_(psf.fixedDmdtf_)
{}
@ -119,7 +114,6 @@ void alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
}
dmdtf_ = (1 - relax_)*dmdtf_ + relax_*fixedDmdtf_;
dmdtLf_ = dmdtf_*L_;
operator==(calcAlphat(*this));
@ -135,7 +129,6 @@ void alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::write
alphatPhaseChangeWallFunctionFvPatchScalarField::write(os);
writeEntry(os, "fixedDmdtf", fixedDmdtf_);
writeEntry(os, "L", L_);
}

View File

@ -62,9 +62,6 @@ class alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
//- Volumetric phase-change mass flux in near wall cells
scalar fixedDmdtf_;
//- Latent heat
scalar L_;
public:

View File

@ -51,8 +51,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
alphatPhaseJayatillekeWallFunctionFvPatchScalarField(p, iF),
otherPhaseName_(word::null),
relax_(1),
dmdtf_(p.size(), 0),
dmdtLf_(p.size(), 0)
dmdtf_(p.size(), 0)
{}
@ -67,8 +66,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
alphatPhaseJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
otherPhaseName_(dict.lookup("otherPhase")),
relax_(dict.lookupOrDefault<scalar>("relax", 1)),
dmdtf_(p.size(), 0),
dmdtLf_(p.size(), 0)
dmdtf_(p.size(), 0)
{
// Check that otherPhaseName != this phase
if (internalField().group() == otherPhaseName_)
@ -85,11 +83,6 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
{
dmdtf_ = scalarField("dmdtf", dict, p.size());
}
if (dict.found("dmdtLf"))
{
dmdtLf_ = scalarField("dmdtLf", dict, p.size());
}
}
@ -105,8 +98,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
alphatPhaseJayatillekeWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
otherPhaseName_(ptf.otherPhaseName_),
relax_(ptf.relax_),
dmdtf_(mapper(ptf.dmdtf_)),
dmdtLf_(mapper(ptf.dmdtLf_))
dmdtf_(mapper(ptf.dmdtf_))
{}
@ -119,8 +111,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
alphatPhaseJayatillekeWallFunctionFvPatchScalarField(awfpsf),
otherPhaseName_(awfpsf.otherPhaseName_),
relax_(awfpsf.relax_),
dmdtf_(awfpsf.dmdtf_),
dmdtLf_(awfpsf.dmdtLf_)
dmdtf_(awfpsf.dmdtf_)
{}
@ -134,8 +125,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
alphatPhaseJayatillekeWallFunctionFvPatchScalarField(awfpsf, iF),
otherPhaseName_(awfpsf.otherPhaseName_),
relax_(awfpsf.relax_),
dmdtf_(awfpsf.dmdtf_),
dmdtLf_(awfpsf.dmdtLf_)
dmdtf_(awfpsf.dmdtf_)
{}
@ -175,32 +165,7 @@ dmdtf(const phasePairKey& phasePair) const
<< " dmdtf requested for invalid phasePair!"
<< abort(FatalError);
return dmdtLf_;
}
}
const scalarField&
alphatPhaseChangeWallFunctionFvPatchScalarField::dmdtLf() const
{
return dmdtLf_;
}
const scalarField& alphatPhaseChangeWallFunctionFvPatchScalarField::
dmdtLf(const phasePairKey& phasePair) const
{
if (activePhasePair(phasePair))
{
return dmdtLf_;
}
else
{
FatalErrorInFunction
<< " dmdtLf requested for invalid phasePair!"
<< abort(FatalError);
return dmdtLf_;
return dmdtf_;
}
}
@ -213,7 +178,6 @@ void alphatPhaseChangeWallFunctionFvPatchScalarField::autoMap
alphatPhaseJayatillekeWallFunctionFvPatchScalarField::autoMap(m);
m(dmdtf_, dmdtf_);
m(dmdtLf_, dmdtLf_);
}
@ -229,7 +193,6 @@ void alphatPhaseChangeWallFunctionFvPatchScalarField::rmap
refCast<const alphatPhaseChangeWallFunctionFvPatchScalarField>(ptf);
dmdtf_.rmap(tiptf.dmdtf_, addr);
dmdtLf_.rmap(tiptf.dmdtLf_, addr);
}
@ -240,7 +203,6 @@ void alphatPhaseChangeWallFunctionFvPatchScalarField::write(Ostream& os) const
writeEntry(os, "otherPhase", otherPhaseName_);
writeEntry(os, "relax", relax_);
writeEntry(os, "dmdtf", dmdtf_);
writeEntry(os, "dmdtLf", dmdtLf_);
}

View File

@ -69,9 +69,6 @@ protected:
//- Rate of phase-change
scalarField dmdtf_;
//- Latent heat of the phase-change
scalarField dmdtLf_;
public:
@ -132,12 +129,6 @@ public:
//- Return the rate of phase-change for specific phase pair
const scalarField& dmdtf(const phasePairKey&) const;
//- Return the enthalpy source due to phase-change
const scalarField& dmdtLf() const;
//- Return the rate of phase-change for specific phase pair
const scalarField& dmdtLf(const phasePairKey&) const;
// Mapping functions

View File

@ -25,6 +25,7 @@ License
#include "alphatWallBoilingWallFunctionFvPatchScalarField.H"
#include "phaseSystem.H"
#include "heatTransferPhaseSystem.H"
#include "compressibleMomentumTransportModel.H"
#include "phaseCompressibleMomentumTransportModel.H"
#include "saturationModel.H"
@ -373,7 +374,6 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
const scalarField& muw = tmuw();
const rhoThermo& lThermo = liquid.thermo();
const rhoThermo& vThermo = vapor.thermo();
const scalarField& alphaw = lThermo.alpha(patchi);
@ -395,9 +395,6 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
const fvPatchScalarField& hew =
lThermo.he().boundaryField()[patchi];
const fvPatchScalarField& pw =
lThermo.p().boundaryField()[patchi];
const fvPatchScalarField& Tw =
lThermo.T().boundaryField()[patchi];
@ -430,39 +427,29 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
const fvPatchScalarField& Tsatw(Tsat.boundaryField()[patchi]);
// Latent heat
scalarField liquidHaw(lThermo.ha(Tc, patchi));
scalarField vaporHaw(vThermo.ha(Tsatw, patchi));
if (volatileSpecie != "none" && isA<rhoReactionThermo>(lThermo))
{
const basicSpecieMixture& composition =
refCast<const rhoReactionThermo>(lThermo).composition();
liquidHaw =
composition.Ha
(
composition.species()[volatileSpecie],
pw,
Tc
);
}
if (volatileSpecie != "none" && isA<rhoReactionThermo>(vThermo))
{
const basicSpecieMixture& composition =
refCast<const rhoReactionThermo>(vThermo).composition();
vaporHaw =
composition.Ha
(
composition.species()[volatileSpecie],
pw,
Tsatw
);
}
const scalarField L(vaporHaw - liquidHaw);
const scalarField L
(
volatileSpecie != "none"
? -refCast<const heatTransferPhaseSystem>(fluid)
.Li
(
pair,
volatileSpecie,
dmdtf_,
Tsat,
patch().faceCells(),
heatTransferPhaseSystem::latentHeatScheme::upwind
)
: -refCast<const heatTransferPhaseSystem>(fluid)
.L
(
pair,
dmdtf_,
Tsat,
patch().faceCells(),
heatTransferPhaseSystem::latentHeatScheme::upwind
)
);
// Liquid phase fraction at the wall
const scalarField liquidw(liquid.boundaryField()[patchi]);
@ -571,10 +558,6 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
(1 - relax_)*dmdtf_
+ relax_*(1.0/6.0)*A2E*dDep_*rhoVaporw*fDep*AbyV_;
// Volumetric source in the near wall cell due to the wall
// boiling
dmdtLf_ = dmdtf_*L;
// Quenching heat transfer coefficient
const scalarField hQ
(
@ -590,6 +573,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
(1 - relax_)*qq_
+ relax_*(A2*hQ*max(Tw - Tl, scalar(0)));
// Evaporation heat flux
const scalarField qe(dmdtf_*L/AbyV_);
// Effective thermal diffusivity that corresponds to the
// calculated convective, quenching and evaporative heat
// fluxes
@ -598,7 +584,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
(
(
A1*alphatConv_
+ (qq_ + qe())/max(hew.snGrad(), scalar(1e-16))
+ (qq_ + qe)/max(hew.snGrad(), scalar(1e-16))
)
/max(liquidw, scalar(1e-8))
);
@ -641,10 +627,10 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
<< gMax(dmdtf_) << endl;
Info<< " qc: " << gMin(qc) << " - " << gMax(qc)
<< endl;
Info<< " qq: " << gMin(fLiquid*qq()) << " - "
<< gMax(fLiquid*qq()) << endl;
Info<< " qe: " << gMin(fLiquid*qe()) << " - "
<< gMax(fLiquid*qe()) << endl;
Info<< " qq: " << gMin(fLiquid*qq_) << " - "
<< gMax(fLiquid*qq_) << endl;
Info<< " qe: " << gMin(fLiquid*qe) << " - "
<< gMax(fLiquid*qe) << endl;
Info<< " qEff: " << gMin(qEff) << " - "
<< gMax(qEff) << endl;
Info<< " alphat: " << gMin(*this) << " - "

View File

@ -282,18 +282,6 @@ public:
return dDep_;
}
//- Return the quenching surface heat flux [W/m^2]
const scalarField& qq() const
{
return qq_;
}
//- Return the evaporation surface heat flux [W/m^2]
tmp<scalarField> qe() const
{
return dmdtLf_/AbyV_;
}
// Mapping functions

View File

@ -10,6 +10,9 @@ phaseSystem/phaseSystem.C
phaseSystem/phaseSystemNew.C
phaseSystem/phaseSystemSolve.C
PhaseSystems/HeatTransferPhaseSystem/heatTransferPhaseSystem.C
PhaseSystems/TwoResistanceHeatTransferPhaseSystem/twoResistanceHeatTransferPhaseSystem.C
diameterModels/diameterModel/diameterModel.C
diameterModels/diameterModel/diameterModelNew.C
diameterModels/sphericalDiameter/sphericalDiameter.C

View File

@ -0,0 +1,768 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2020 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "HeatTransferPhaseSystem.H"
#include "fvmSup.H"
#include "rhoReactionThermo.H"
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class BasePhaseSystem>
void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
(
const phaseSystem::dmdtfTable& dmdtfs,
phaseSystem::heatTransferTable& eqns
) const
{
// Loop the pairs
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
{
const phasePairKey& key = dmdtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
const volScalarField dmdtf21(posPart(dmdtf));
const volScalarField dmdtf12(negPart(dmdtf));
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const volScalarField& he1 = thermo1.he();
const volScalarField& he2 = thermo2.he();
const volScalarField hs1(thermo1.hs());
const volScalarField hs2(thermo2.hs());
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Transfer of sensible enthalpy within the phases
*eqns[phase1.name()] +=
dmdtf*hs1 + fvm::Sp(dmdtf12, he1) - dmdtf12*he1;
*eqns[phase2.name()] -=
dmdtf*hs2 + fvm::Sp(dmdtf21, he2) - dmdtf21*he2;
// Transfer of sensible enthalpy between the phases
*eqns[phase1.name()] += dmdtf21*(hs2 - hs1);
*eqns[phase2.name()] -= dmdtf12*(hs1 - hs2);
// Transfer of kinetic energy
*eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2;
}
}
template<class BasePhaseSystem>
void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs
(
const phaseSystem::dmidtfTable& dmidtfs,
phaseSystem::heatTransferTable& eqns
) const
{
static const dimensionedScalar one(dimless, 1);
// Loop the pairs
forAllConstIter(phaseSystem::dmidtfTable, dmidtfs, dmidtfIter)
{
const phasePairKey& key = dmidtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const basicSpecieMixture* compositionPtr1 =
isA<rhoReactionThermo>(thermo1)
? &refCast<const rhoReactionThermo>(thermo1).composition()
: static_cast<const basicSpecieMixture*>(nullptr);
const basicSpecieMixture* compositionPtr2 =
isA<rhoReactionThermo>(thermo2)
? &refCast<const rhoReactionThermo>(thermo2).composition()
: static_cast<const basicSpecieMixture*>(nullptr);
const volScalarField& he1 = thermo1.he();
const volScalarField& he2 = thermo2.he();
const volScalarField hs1(thermo1.hs());
const volScalarField hs2(thermo2.hs());
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Loop the species
forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
{
const word& specie = dmidtfJter.key();
// Mass transfer rates
const volScalarField dmidtf
(
Pair<word>::compare(pair, key)**dmidtfJter()
);
const volScalarField dmidtf21(posPart(dmidtf));
const volScalarField dmidtf12(negPart(dmidtf));
// Specie indices
const label speciei1 =
compositionPtr1 ? compositionPtr1->species()[specie] : -1;
const label speciei2 =
compositionPtr2 ? compositionPtr2->species()[specie] : -1;
// Enthalpies
const volScalarField hsi1
(
compositionPtr1
? compositionPtr1->Hs(speciei1, thermo1.p(), thermo1.T())
: tmp<volScalarField>(hs1)
);
const volScalarField hsi2
(
compositionPtr2
? compositionPtr2->Hs(speciei2, thermo2.p(), thermo2.T())
: tmp<volScalarField>(hs2)
);
// Limited mass fractions
tmp<volScalarField> tYi1, tYi2;
if (residualY_ > 0)
{
tYi1 =
compositionPtr1
? max(compositionPtr1->Y(speciei1), residualY_)
: volScalarField::New("Yi1", this->mesh(), one);
tYi2 =
compositionPtr2
? max(compositionPtr2->Y(speciei2), residualY_)
: volScalarField::New("Yi2", this->mesh(), one);
}
// Transfer of sensible enthalpy within the phases
*eqns[phase1.name()] += dmidtf*hsi1;
*eqns[phase2.name()] -= dmidtf*hsi2;
if (residualY_ > 0)
{
*eqns[phase1.name()] +=
fvm::Sp(dmidtf12/tYi1(), he1) - dmidtf12/tYi1()*he1;
*eqns[phase2.name()] -=
fvm::Sp(dmidtf21/tYi2(), he2) - dmidtf21/tYi2()*he2;
}
// Transfer of sensible enthalpy between the phases
*eqns[phase1.name()] += dmidtf21*(hsi2 - hsi1);
*eqns[phase2.name()] -= dmidtf12*(hsi1 - hsi2);
// Transfer of kinetic energy
*eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1;
*eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2;
}
}
}
template<class BasePhaseSystem>
void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefsWithoutL
(
const phaseSystem::dmdtfTable& dmdtfs,
const phaseSystem::dmdtfTable& Tfs,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const
{
// Loop the pairs
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
{
const phasePairKey& key = dmdtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
const volScalarField dmdtf21(posPart(dmdtf));
const volScalarField dmdtf12(negPart(dmdtf));
const volScalarField& Tf = *Tfs[key];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const volScalarField& he1 = thermo1.he();
const volScalarField& he2 = thermo2.he();
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Interface enthalpies
const volScalarField hsf1(thermo1.hs(thermo1.p(), Tf));
const volScalarField hsf2(thermo2.hs(thermo1.p(), Tf));
// Transfer of energy from the interface into the bulk
switch (scheme)
{
case latentHeatScheme::symmetric:
{
*eqns[phase1.name()] += dmdtf*hsf1;
*eqns[phase2.name()] -= dmdtf*hsf2;
break;
}
case latentHeatScheme::upwind:
{
// Bulk enthalpies
const volScalarField hs1(thermo1.hs());
const volScalarField hs2(thermo2.hs());
*eqns[phase1.name()] += dmdtf21*hsf1 + dmdtf12*hs1;
*eqns[phase2.name()] -= dmdtf12*hsf2 + dmdtf21*hs2;
break;
}
}
*eqns[phase1.name()] += fvm::Sp(dmdtf12, he1) - dmdtf12*he1;
*eqns[phase2.name()] -= fvm::Sp(dmdtf21, he2) - dmdtf21*he2;
// Transfer of kinetic energy
*eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2;
}
}
template<class BasePhaseSystem>
void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmdtL
(
const phaseSystem::dmdtfTable& dmdtfs,
const phaseSystem::dmdtfTable& Tfs,
const scalar weight,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const
{
// Loop the pairs
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
{
const phasePairKey& key = dmdtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
const volScalarField dmdtf21(posPart(dmdtf));
const volScalarField dmdtf12(negPart(dmdtf));
const volScalarField& Tf = *Tfs[key];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
// Latent heat contribution
const volScalarField L(this->L(pair, dmdtf, Tf, scheme));
*eqns[phase1.name()] += ((1 - weight)*dmdtf12 + weight*dmdtf21)*L;
*eqns[phase2.name()] += ((1 - weight)*dmdtf21 + weight*dmdtf12)*L;
}
}
template<class BasePhaseSystem>
void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
(
const phaseSystem::dmdtfTable& dmdtfs,
const phaseSystem::dmdtfTable& Tfs,
const scalar weight,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const
{
addDmdtHefsWithoutL(dmdtfs, Tfs, scheme, eqns);
addDmdtL(dmdtfs, Tfs, weight, scheme, eqns);
}
template<class BasePhaseSystem>
void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefsWithoutL
(
const phaseSystem::dmidtfTable& dmidtfs,
const phaseSystem::dmdtfTable& Tfs,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const
{
static const dimensionedScalar one(dimless, 1);
// Loop the pairs
forAllConstIter(phaseSystem::dmidtfTable, dmidtfs, dmidtfIter)
{
const phasePairKey& key = dmidtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const volScalarField& Tf = *Tfs[key];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const basicSpecieMixture* compositionPtr1 =
isA<rhoReactionThermo>(thermo1)
? &refCast<const rhoReactionThermo>(thermo1).composition()
: static_cast<const basicSpecieMixture*>(nullptr);
const basicSpecieMixture* compositionPtr2 =
isA<rhoReactionThermo>(thermo2)
? &refCast<const rhoReactionThermo>(thermo2).composition()
: static_cast<const basicSpecieMixture*>(nullptr);
const volScalarField& he1 = thermo1.he();
const volScalarField& he2 = thermo2.he();
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Interface enthalpies
const volScalarField hsf1(thermo1.hs(thermo1.p(), Tf));
const volScalarField hsf2(thermo2.hs(thermo2.p(), Tf));
// Loop the species
forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
{
const word& specie = dmidtfJter.key();
// Mass transfer rates
const volScalarField dmidtf
(
Pair<word>::compare(pair, key)**dmidtfJter()
);
const volScalarField dmidtf21(posPart(dmidtf));
const volScalarField dmidtf12(negPart(dmidtf));
// Specie indices
const label speciei1 =
compositionPtr1 ? compositionPtr1->species()[specie] : -1;
const label speciei2 =
compositionPtr2 ? compositionPtr2->species()[specie] : -1;
// Interface enthalpies
const volScalarField hsfi1
(
compositionPtr1
? compositionPtr1->Hs(speciei1, thermo1.p(), Tf)
: tmp<volScalarField>(hsf1)
);
const volScalarField hsfi2
(
compositionPtr2
? compositionPtr2->Hs(speciei2, thermo2.p(), Tf)
: tmp<volScalarField>(hsf2)
);
// Limited mass fractions
tmp<volScalarField> tYi1, tYi2;
if (this->residualY_ > 0)
{
tYi1 =
compositionPtr1
? max(compositionPtr1->Y(speciei1), this->residualY_)
: volScalarField::New("Yi1", this->mesh(), one);
tYi2 =
compositionPtr2
? max(compositionPtr2->Y(speciei2), this->residualY_)
: volScalarField::New("Yi2", this->mesh(), one);
}
// Transfer of energy from the interface into the bulk
switch (scheme)
{
case latentHeatScheme::symmetric:
{
*eqns[phase1.name()] += dmidtf*hsfi1;
*eqns[phase2.name()] -= dmidtf*hsfi2;
break;
}
case latentHeatScheme::upwind:
{
// Bulk enthalpies
const volScalarField hsi1
(
compositionPtr1
? compositionPtr1->Hs(speciei1, thermo1.p(), thermo1.T())
: thermo1.hs()
);
const volScalarField hsi2
(
compositionPtr2
? compositionPtr2->Hs(speciei2, thermo2.p(), thermo2.T())
: thermo2.hs()
);
*eqns[phase1.name()] += dmidtf21*hsfi1 + dmidtf12*hsi1;
*eqns[phase2.name()] -= dmidtf12*hsfi2 + dmidtf21*hsi2;
break;
}
}
if (this->residualY_ > 0)
{
*eqns[phase1.name()] +=
fvm::Sp(dmidtf12/tYi1(), he1) - dmidtf12/tYi1()*he1;
}
if (this->residualY_ > 0)
{
*eqns[phase2.name()] -=
fvm::Sp(dmidtf21/tYi2(), he2) - dmidtf21/tYi2()*he2;
}
// Transfer of kinetic energy
*eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1;
*eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2;
}
}
}
template<class BasePhaseSystem>
void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtL
(
const phaseSystem::dmidtfTable& dmidtfs,
const phaseSystem::dmdtfTable& Tfs,
const scalar weight,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const
{
// Loop the pairs
forAllConstIter(phaseSystem::dmidtfTable, dmidtfs, dmidtfIter)
{
const phasePairKey& key = dmidtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const volScalarField& Tf = *Tfs[key];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
// Loop the species
forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
{
const word& specie = dmidtfJter.key();
// Mass transfer rates
const volScalarField dmidtf
(
Pair<word>::compare(pair, key)**dmidtfJter()
);
const volScalarField dmidtf21(posPart(dmidtf));
const volScalarField dmidtf12(negPart(dmidtf));
// Latent heat contribution
const volScalarField Li(this->Li(pair, specie, dmidtf, Tf, scheme));
*eqns[phase1.name()] +=
((1 - weight)*dmidtf12 + weight*dmidtf21)*Li;
*eqns[phase2.name()] +=
((1 - weight)*dmidtf21 + weight*dmidtf12)*Li;
}
}
}
template<class BasePhaseSystem>
void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs
(
const phaseSystem::dmidtfTable& dmidtfs,
const phaseSystem::dmdtfTable& Tfs,
const scalar weight,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const
{
addDmidtHefsWithoutL(dmidtfs, Tfs, scheme, eqns);
addDmidtL(dmidtfs, Tfs, weight, scheme, eqns);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::HeatTransferPhaseSystem
(
const fvMesh& mesh
)
:
heatTransferPhaseSystem(),
BasePhaseSystem(mesh),
residualY_(this->template lookupOrDefault<scalar>("residualY", -1))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::~HeatTransferPhaseSystem()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::L
(
const phasePair& pair,
const volScalarField& dmdtf,
const volScalarField& Tf,
const latentHeatScheme scheme
) const
{
const rhoThermo& thermo1 = pair.phase1().thermo();
const rhoThermo& thermo2 = pair.phase2().thermo();
// Interface enthalpies
const volScalarField haf1(thermo1.ha(thermo1.p(), Tf));
const volScalarField haf2(thermo2.ha(thermo2.p(), Tf));
switch (scheme)
{
case latentHeatScheme::symmetric:
{
return haf2 - haf1;
}
case latentHeatScheme::upwind:
{
// Bulk enthalpies
const volScalarField ha1(thermo1.ha());
const volScalarField ha2(thermo2.ha());
return
neg0(dmdtf)*haf2 + pos(dmdtf)*ha2
- pos0(dmdtf)*haf1 - neg(dmdtf)*ha1;
}
}
return tmp<volScalarField>(nullptr);
}
template<class BasePhaseSystem>
Foam::tmp<Foam::scalarField>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::L
(
const phasePair& pair,
const scalarField& dmdtf,
const scalarField& Tf,
const labelUList& cells,
const latentHeatScheme scheme
) const
{
const rhoThermo& thermo1 = pair.phase1().thermo();
const rhoThermo& thermo2 = pair.phase2().thermo();
// Interface enthalpies
const scalarField haf1(thermo1.ha(Tf, cells));
const scalarField haf2(thermo2.ha(Tf, cells));
switch (scheme)
{
case latentHeatScheme::symmetric:
{
return haf2 - haf1;
}
case latentHeatScheme::upwind:
{
const scalarField T1(UIndirectList<scalar>(thermo1.T(), cells));
const scalarField T2(UIndirectList<scalar>(thermo2.T(), cells));
// Bulk enthalpies
const scalarField ha1(thermo1.ha(T1, cells));
const scalarField ha2(thermo2.ha(T2, cells));
return
neg0(dmdtf)*haf2 + pos(dmdtf)*ha2
- pos0(dmdtf)*haf1 - neg(dmdtf)*ha1;
}
}
return tmp<scalarField>(nullptr);
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
(
const phasePair& pair,
const word& specie,
const volScalarField& dmdtf,
const volScalarField& Tf,
const latentHeatScheme scheme
) const
{
const rhoThermo& thermo1 = pair.phase1().thermo();
const rhoThermo& thermo2 = pair.phase2().thermo();
const basicSpecieMixture* compositionPtr1 =
isA<rhoReactionThermo>(thermo1)
? &refCast<const rhoReactionThermo>(thermo1).composition()
: static_cast<const basicSpecieMixture*>(nullptr);
const basicSpecieMixture* compositionPtr2 =
isA<rhoReactionThermo>(thermo2)
? &refCast<const rhoReactionThermo>(thermo2).composition()
: static_cast<const basicSpecieMixture*>(nullptr);
const label speciei1 =
compositionPtr1 ? compositionPtr1->species()[specie] : -1;
const label speciei2 =
compositionPtr2 ? compositionPtr2->species()[specie] : -1;
// Interface enthalpies
const volScalarField hafi1
(
compositionPtr1
? compositionPtr1->Ha(speciei1, thermo1.p(), Tf)
: thermo1.ha(thermo1.p(), Tf)
);
const volScalarField hafi2
(
compositionPtr2
? compositionPtr2->Ha(speciei2, thermo2.p(), Tf)
: thermo2.ha(thermo1.p(), Tf)
);
switch (scheme)
{
case latentHeatScheme::symmetric:
{
return hafi2 - hafi1;
}
case latentHeatScheme::upwind:
{
// Bulk enthalpies
const volScalarField hai1
(
compositionPtr1
? compositionPtr1->Ha(speciei1, thermo1.p(), thermo1.T())
: thermo1.ha()
);
const volScalarField hai2
(
compositionPtr2
? compositionPtr2->Ha(speciei2, thermo2.p(), thermo2.T())
: thermo2.ha()
);
return
neg0(dmdtf)*hafi2 + pos(dmdtf)*hai2
- pos0(dmdtf)*hafi1 - neg(dmdtf)*hai1;
}
}
return tmp<volScalarField>(nullptr);
}
template<class BasePhaseSystem>
Foam::tmp<Foam::scalarField>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
(
const phasePair& pair,
const word& specie,
const scalarField& dmdtf,
const scalarField& Tf,
const labelUList& cells,
const latentHeatScheme scheme
) const
{
const rhoThermo& thermo1 = pair.phase1().thermo();
const rhoThermo& thermo2 = pair.phase2().thermo();
const basicSpecieMixture* compositionPtr1 =
isA<rhoReactionThermo>(thermo1)
? &refCast<const rhoReactionThermo>(thermo1).composition()
: static_cast<const basicSpecieMixture*>(nullptr);
const basicSpecieMixture* compositionPtr2 =
isA<rhoReactionThermo>(thermo2)
? &refCast<const rhoReactionThermo>(thermo2).composition()
: static_cast<const basicSpecieMixture*>(nullptr);
const label speciei1 =
compositionPtr1 ? compositionPtr1->species()[specie] : -1;
const label speciei2 =
compositionPtr2 ? compositionPtr2->species()[specie] : -1;
const scalarField p1(UIndirectList<scalar>(thermo1.p(), cells));
const scalarField p2(UIndirectList<scalar>(thermo2.p(), cells));
// Interface enthalpies
const scalarField hafi1
(
compositionPtr1
? compositionPtr1->Ha(speciei1, p1, Tf)
: thermo1.ha(Tf, cells)
);
const scalarField hafi2
(
compositionPtr2
? compositionPtr2->Ha(speciei2, p2, Tf)
: thermo2.ha(Tf, cells)
);
switch (scheme)
{
case latentHeatScheme::symmetric:
{
return hafi2 - hafi1;
}
case latentHeatScheme::upwind:
{
const scalarField T1(UIndirectList<scalar>(thermo1.T(), cells));
const scalarField T2(UIndirectList<scalar>(thermo2.T(), cells));
// Bulk enthalpies
const scalarField hai1
(
compositionPtr1
? compositionPtr1->Ha(speciei1, p1, T1)
: thermo1.ha(T1, cells)
);
const scalarField hai2
(
compositionPtr2
? compositionPtr2->Ha(speciei2, p2, T2)
: thermo2.ha(T2, cells)
);
return
neg0(dmdtf)*hafi2 + pos(dmdtf)*hai2
- pos0(dmdtf)*hafi1 - neg(dmdtf)*hai1;
}
}
return tmp<scalarField>(nullptr);
}
template<class BasePhaseSystem>
bool Foam::HeatTransferPhaseSystem<BasePhaseSystem>::read()
{
if (BasePhaseSystem::read())
{
bool readOK = true;
// Models ...
return readOK;
}
else
{
return false;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,226 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::HeatTransferPhaseSystem
Description
...
SourceFiles
HeatTransferPhaseSystem.C
\*---------------------------------------------------------------------------*/
#ifndef HeatTransferPhaseSystem_H
#define HeatTransferPhaseSystem_H
#include "heatTransferPhaseSystem.H"
#include "phaseSystem.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class HeatTransferPhaseSystem Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseSystem>
class HeatTransferPhaseSystem
:
public heatTransferPhaseSystem,
public BasePhaseSystem
{
protected:
// Protected Member Functions
//- Add energy transfer terms which result from bulk mass transfers
void addDmdtHefs
(
const phaseSystem::dmdtfTable& dmdtfs,
phaseSystem::heatTransferTable& eqns
) const;
//- Add energy transfer terms which result from specie mass transfers
void addDmidtHefs
(
const phaseSystem::dmidtfTable& dmidtfs,
phaseSystem::heatTransferTable& eqns
) const;
//- Add energy transfer terms which result from bulk phase changes,
// without the latent heat contribution
void addDmdtHefsWithoutL
(
const phaseSystem::dmdtfTable& dmdtfs,
const phaseSystem::dmdtfTable& Tfs,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const;
//- Add latent heat terms which result from bulk phase changes.
// Weight is the proportion of the latent heat contribution applied to
// the downwind side of the mass transfer.
void addDmdtL
(
const phaseSystem::dmdtfTable& dmdtfs,
const phaseSystem::dmdtfTable& Tfs,
const scalar weight,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const;
//- Add energy transfer terms which result from bulk phase changes.
// Weight is the proportion of the latent heat contribution applied to
// the downwind side of the mass transfer.
void addDmdtHefs
(
const phaseSystem::dmdtfTable& dmdtfs,
const phaseSystem::dmdtfTable& Tfs,
const scalar weight,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const;
//- Add energy transfer terms which result from specie phase changes,
// without the latent heat contribution
void addDmidtHefsWithoutL
(
const phaseSystem::dmidtfTable& dmidtfs,
const phaseSystem::dmdtfTable& Tfs,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const;
//- Add latent heat terms which result from specie phase changes.
// Weight is the proportion of the latent heat contribution applied to
// the downwind side of the mass transfer.
void addDmidtL
(
const phaseSystem::dmidtfTable& dmidtfs,
const phaseSystem::dmdtfTable& Tfs,
const scalar weight,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const;
//- Add energy transfer terms which result from specie phase changes
// Weight is the proportion of the latent heat contribution applied to
// the downwind side of the mass transfer.
void addDmidtHefs
(
const phaseSystem::dmidtfTable& dmidtfs,
const phaseSystem::dmdtfTable& Tfs,
const scalar weight,
const latentHeatScheme scheme,
phaseSystem::heatTransferTable& eqns
) const;
// Protected data
//- Residual mass fraction used for linearisation of heat transfer
// terms that result from mass transfers of individual species
scalar residualY_;
public:
// Constructors
//- Construct from fvMesh
HeatTransferPhaseSystem(const fvMesh&);
//- Destructor
virtual ~HeatTransferPhaseSystem();
// Member Functions
//- Return the latent heat for a given pair, mass transfer rate (used
// only for it's sign), and interface temperature
virtual tmp<volScalarField> L
(
const phasePair& pair,
const volScalarField& dmdtf,
const volScalarField& Tf,
const latentHeatScheme scheme
) const;
//- As above, but for a cell-set
virtual tmp<scalarField> L
(
const phasePair& pair,
const scalarField& dmdtf,
const scalarField& Tf,
const labelUList& cells,
const latentHeatScheme scheme
) const;
//- Return the latent heat for a given pair, specie, mass transfer rate
// (used only for it's sign), and interface temperature
virtual tmp<volScalarField> Li
(
const phasePair& pair,
const word& member,
const volScalarField& dmidtf,
const volScalarField& Tf,
const latentHeatScheme scheme
) const;
//- As above, but for a cell-set
virtual tmp<scalarField> Li
(
const phasePair& pair,
const word& member,
const scalarField& dmidtf,
const scalarField& Tf,
const labelUList& cells,
const latentHeatScheme scheme
) const;
//- Read base phaseProperties dictionary
virtual bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "HeatTransferPhaseSystem.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,67 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2020 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "heatTransferPhaseSystem.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
template<>
const char* NamedEnum
<
Foam::heatTransferPhaseSystem::latentHeatScheme,
2
>::names[] = {"symmetric", "upwind"};
}
const Foam::NamedEnum
<
Foam::heatTransferPhaseSystem::latentHeatScheme,
2
>
Foam::heatTransferPhaseSystem::latentHeatSchemeNames_;
namespace Foam
{
defineTypeNameAndDebug(heatTransferPhaseSystem, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::heatTransferPhaseSystem::heatTransferPhaseSystem()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::heatTransferPhaseSystem::~heatTransferPhaseSystem()
{}
// ************************************************************************* //

View File

@ -0,0 +1,133 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2020 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::heatTransferPhaseSystem
SourceFiles
heatTransferPhaseSystem.C
\*---------------------------------------------------------------------------*/
#ifndef heatTransferPhaseSystem_H
#define heatTransferPhaseSystem_H
#include "NamedEnum.H"
#include "primitiveFieldsFwd.H"
#include "volFieldsFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class phasePair;
/*---------------------------------------------------------------------------*\
Class heatTransferPhaseSystem Declaration
\*---------------------------------------------------------------------------*/
class heatTransferPhaseSystem
{
public:
//- Enumeration for the latent heat scheme
enum class latentHeatScheme
{
symmetric,
upwind
};
//- Names of the latent heat schemes
static const NamedEnum<latentHeatScheme, 2> latentHeatSchemeNames_;
//- Runtime type information
TypeName("heatTransferPhaseSystem");
// Constructors
//- Default constructor
heatTransferPhaseSystem();
//- Destructor
virtual ~heatTransferPhaseSystem();
// Member Functions
//- Return the latent heat for a given pair, mass transfer rate (used
// only for it's sign), and interface temperature
virtual tmp<volScalarField> L
(
const phasePair& pair,
const volScalarField& dmdtf,
const volScalarField& Tf,
const latentHeatScheme scheme
) const = 0;
//- As above, but for a cell-set
virtual tmp<scalarField> L
(
const phasePair& pair,
const scalarField& dmdtf,
const scalarField& Tf,
const labelUList& cells,
const latentHeatScheme scheme
) const = 0;
//- Return the latent heat for a given pair, specie, mass transfer rate
// (used only for it's sign), and interface temperature
virtual tmp<volScalarField> Li
(
const phasePair& pair,
const word& member,
const volScalarField& dmidtf,
const volScalarField& Tf,
const latentHeatScheme scheme
) const = 0;
//- As above, but for a cell-set
virtual tmp<scalarField> Li
(
const phasePair& pair,
const word& member,
const scalarField& dmidtf,
const scalarField& Tf,
const labelUList& cells,
const latentHeatScheme scheme
) const = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -257,8 +257,30 @@ InterfaceCompositionPhaseChangePhaseSystem
this->phasePairs_[interfaceCompositionModelIter.key()];
dmidtfSus_.insert(pair, new HashPtrTable<volScalarField>());
dmidtfSps_.insert(pair, new HashPtrTable<volScalarField>());
Tfs_.insert
(
pair,
new volScalarField
(
IOobject
(
IOobject::groupName
(
"interfaceCompositionPhaseChange:Tf",
pair.name()
),
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
(pair.phase1().thermo().T() + pair.phase2().thermo().T())/2
)
);
forAllConstIter(phasePair, pair, pairIter)
{
const autoPtr<interfaceCompositionModel>& compositionModelPtr =
@ -454,7 +476,14 @@ heatTransfer() const
phaseSystem::heatTransferTable& eqns = eqnsPtr();
this->addDmidtHef(dmidtfs(), eqns);
this->addDmidtHefs
(
dmidtfs(),
Tfs_,
latentHeatScheme::symmetric,
latentHeatTransfer::mass,
eqns
);
return eqnsPtr;
}
@ -545,7 +574,7 @@ correct()
const phasePair& pair =
this->phasePairs_[interfaceCompositionModelIter.key()];
const volScalarField& Tf(*this->Tf_[pair]);
const volScalarField& Tf(*this->Tfs_[pair]);
forAllConstIter(phasePair, pair, pairIter)
{
@ -624,7 +653,13 @@ correctInterfaceThermo()
const volScalarField H2(this->heatTransferModels_[pair].second()->K());
const dimensionedScalar HSmall("small", heatTransferModel::dimK, small);
volScalarField& Tf = *this->Tf_[pair];
const typename diffusiveMassTransferModelTable::value_type&
diffusiveMassTransfers = this->diffusiveMassTransferModels_[pair];
const typename interfaceCompositionModelTable::value_type &
interfaceCompositions = this->interfaceCompositionModels_[pair];
volScalarField& Tf = *this->Tfs_[pair];
for (label i = 0; i < nInterfaceCorrectors_; ++ i)
{
@ -644,25 +679,58 @@ correctInterfaceThermo()
);
// Add latent heats from forward and backward models
if (this->interfaceCompositionModels_[pair].first().valid())
forAllConstIter(phasePair, pair, pairIter)
{
this->interfaceCompositionModels_[pair].first()->addDmdtL
(
diffusiveMassTransferModels_[pair].first()->K(),
Tf,
dmdtLf.ref(),
dmdtLfPrime.ref()
);
}
if (this->interfaceCompositionModels_[pair].second().valid())
{
this->interfaceCompositionModels_[pair].second()->addDmdtL
(
- diffusiveMassTransferModels_[pair].second()->K(),
Tf,
dmdtLf.ref(),
dmdtLfPrime.ref()
);
if (interfaceCompositions[pairIter.index()].valid())
{
const BlendedInterfacialModel<diffusiveMassTransferModel>&
diffusiveMassTransfer =
diffusiveMassTransfers[pairIter.index()];
const interfaceCompositionModel&
interfaceComposition =
interfaceCompositions[pairIter.index()];
const label sign = pairIter.index() == 0 ? 1 : -1;
forAllConstIter
(
hashedWordList,
interfaceComposition.species(),
specieIter
)
{
const word& specie = *specieIter;
const volScalarField dY
(
interfaceComposition.dY(specie, Tf)
);
const volScalarField dYfPrime
(
interfaceComposition.dYfPrime(specie, Tf)
);
const volScalarField rhoKDL
(
pairIter().thermo().rho()
*diffusiveMassTransfer.K()
*interfaceComposition.D(specie)
*this->Li
(
pair,
specie,
dY,
Tf,
latentHeatScheme::symmetric
)
);
dmdtLf.ref() += sign*rhoKDL*dY;
dmdtLfPrime.ref() += sign*rhoKDL*dYfPrime;
}
}
}
// Update the interface temperature by applying one step of newton's

View File

@ -75,6 +75,10 @@ class InterfaceCompositionPhaseChangePhaseSystem
phasePairKey::hash
> diffusiveMassTransferModelTable;
using latentHeatScheme = typename BasePhaseSystem::latentHeatScheme;
using latentHeatTransfer = typename BasePhaseSystem::latentHeatTransfer;
// Private data
@ -95,6 +99,9 @@ class InterfaceCompositionPhaseChangePhaseSystem
//- The implicit part of the interfacial mass transfer rates
phaseSystem::dmidtfTable dmidtfSps_;
//- Interface temperatures
phaseSystem::dmdtfTable Tfs_;
// Private Member Functions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -27,116 +27,6 @@ License
#include "BlendedInterfacialModel.H"
#include "heatTransferModel.H"
#include "fvmSup.H"
#include "rhoReactionThermo.H"
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class BasePhaseSystem>
void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
(
const phaseSystem::dmdtfTable& dmdtfs,
phaseSystem::heatTransferTable& eqns
) const
{
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
{
const phasePairKey& key = dmdtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
const volScalarField dmdtf21(posPart(dmdtf));
const volScalarField dmdtf12(negPart(dmdtf));
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const volScalarField& he1(thermo1.he());
const volScalarField& he2(thermo2.he());
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Transfer of energy from bulk to bulk
*eqns[phase1.name()] += dmdtf21*he2 + fvm::Sp(dmdtf12, he1);
*eqns[phase2.name()] -= dmdtf12*he1 + fvm::Sp(dmdtf21, he2);
// Transfer of kinetic energy
*eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2;
}
}
template<class BasePhaseSystem>
void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
(
const phaseSystem::dmidtfTable& dmidtfs,
phaseSystem::heatTransferTable& eqns
) const
{
forAllConstIter(phaseSystem::dmidtfTable, dmidtfs, dmidtfIter)
{
const phasePairKey& key = dmidtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const volScalarField& he1(thermo1.he());
const volScalarField& he2(thermo2.he());
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
{
const word& member = dmidtfJter.key();
const volScalarField dmidtf
(
Pair<word>::compare(pair, key)**dmidtfJter()
);
const volScalarField dmidtf21(posPart(dmidtf));
const volScalarField dmidtf12(negPart(dmidtf));
// Create the energies for the transferring specie
volScalarField hei1(he1);
if (isA<rhoReactionThermo>(thermo1))
{
const basicSpecieMixture& composition1 =
refCast<const rhoReactionThermo>(thermo1).composition();
hei1 =
composition1.HE
(
composition1.species()[member],
thermo1.p(),
thermo1.T()
);
}
volScalarField hei2(he2);
if (isA<rhoReactionThermo>(thermo2))
{
const basicSpecieMixture& composition2 =
refCast<const rhoReactionThermo>(thermo2).composition();
hei2 =
composition2.HE
(
composition2.species()[member],
thermo2.p(),
thermo2.T()
);
}
// Transfer of energy from bulk to bulk
*eqns[phase1.name()] += dmidtf21*hei2 + fvm::Sp(dmidtf12, he1);
*eqns[phase2.name()] -= dmidtf12*hei1 + fvm::Sp(dmidtf21, he2);
// Transfer of kinetic energy
*eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1;
*eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -147,7 +37,7 @@ OneResistanceHeatTransferPhaseSystem
const fvMesh& mesh
)
:
BasePhaseSystem(mesh)
HeatTransferPhaseSystem<BasePhaseSystem>(mesh)
{
this->generatePairsAndSubModels
(
@ -208,8 +98,8 @@ heatTransfer() const
const phaseModel& phase = iter();
const phaseModel& otherPhase = iter.otherPhase();
const volScalarField& he(phase.thermo().he());
volScalarField Cpv(phase.thermo().Cpv());
const volScalarField& he = phase.thermo().he();
const volScalarField Cpv(phase.thermo().Cpv());
*eqns[phase.name()] +=
K*(otherPhase.thermo().T() - phase.thermo().T() + he/Cpv)
@ -224,7 +114,7 @@ heatTransfer() const
template<class BasePhaseSystem>
bool Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::read()
{
if (BasePhaseSystem::read())
if (HeatTransferPhaseSystem<BasePhaseSystem>::read())
{
bool readOK = true;

View File

@ -39,7 +39,7 @@ SourceFiles
#ifndef OneResistanceHeatTransferPhaseSystem_H
#define OneResistanceHeatTransferPhaseSystem_H
#include "phaseSystem.H"
#include "HeatTransferPhaseSystem.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -57,7 +57,7 @@ class heatTransferModel;
template<class BasePhaseSystem>
class OneResistanceHeatTransferPhaseSystem
:
public BasePhaseSystem
public HeatTransferPhaseSystem<BasePhaseSystem>
{
protected:
@ -79,23 +79,6 @@ protected:
heatTransferModelTable heatTransferModels_;
// Protected Member Functions
//- Add energy transfer terms which result from bulk mass transfers
void addDmdtHefs
(
const phaseSystem::dmdtfTable& dmdtfs,
phaseSystem::heatTransferTable& eqns
) const;
//- Add energy transfer terms which result from specie mass transfers
void addDmidtHef
(
const phaseSystem::dmidtfTable& dmidtfs,
phaseSystem::heatTransferTable& eqns
) const;
public:
// Constructors

View File

@ -332,7 +332,7 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
phaseSystem::heatTransferTable& eqns = eqnsPtr();
this->addDmdtHefs(dmdtfs_, eqns);
this->addDmidtHef(dmidtfs_, eqns);
this->addDmidtHefs(dmidtfs_, eqns);
return eqnsPtr;
}

View File

@ -32,29 +32,28 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::totalDmdtf
void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::addDmdts
(
const phasePairKey& key
PtrList<volScalarField>& dmdts
) const
{
tmp<volScalarField> tTotalDmdtf = phaseSystem::dmdtf(key);
volScalarField& totalDmdtf = tTotalDmdtf.ref();
const scalar tTotalDmdtfSign =
Pair<word>::compare(this->phasePairs_[key], key);
if (dmdtfs_.found(key))
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
{
totalDmdtf += *dmdtfs_[key];
const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
const volScalarField& dmdtf = *dmdtfIter();
addField(pair.phase1(), "dmdt", dmdtf, dmdts);
addField(pair.phase2(), "dmdt", - dmdtf, dmdts);
}
if (nDmdtfs_.found(key))
forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
{
totalDmdtf += *nDmdtfs_[key];
}
const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
const volScalarField& nDmdtf = *nDmdtfIter();
return tTotalDmdtfSign*tTotalDmdtf;
addField(pair.phase1(), "dmdt", nDmdtf, dmdts);
addField(pair.phase2(), "dmdt", - nDmdtf, dmdts);
}
}
@ -77,21 +76,43 @@ ThermalPhaseChangePhaseSystem
saturationModels_
);
// Check that models have been specified in the correct combinations
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
saturationModelTable,
saturationModels_,
saturationModelIter
)
{
const phasePair& pair(phasePairIter());
const phasePair& pair =
this->phasePairs_[saturationModelIter.key()];
if (pair.ordered())
if
(
!this->heatTransferModels_.found(pair)
|| !this->heatTransferModels_[pair].first().valid()
|| !this->heatTransferModels_[pair].second().valid()
)
{
continue;
FatalErrorInFunction
<< "A heat transfer model for both sides of the " << pair
<< "pair is not specified. This is required by the "
<< "corresponding saturation model"
<< exit(FatalError);
}
}
// Generate interfacial mass transfer fields, initially assumed to be zero
forAllConstIter
(
saturationModelTable,
saturationModels_,
saturationModelIter
)
{
const phasePair& pair =
this->phasePairs_[saturationModelIter.key()];
// Initially assume no mass transfer
dmdtfs_.insert
(
pair,
@ -114,7 +135,27 @@ ThermalPhaseChangePhaseSystem
)
);
// Initially assume no mass transfer
Tfs_.insert
(
pair,
new volScalarField
(
IOobject
(
IOobject::groupName
(
"thermalPhaseChange:Tf",
pair.name()
),
this->mesh().time().timeName(),
this->mesh(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
(pair.phase1().thermo().T() + pair.phase2().thermo().T())/2
)
);
nDmdtfs_.insert
(
pair,
@ -136,29 +177,6 @@ ThermalPhaseChangePhaseSystem
dimensionedScalar(dimDensity/dimTime, 0)
)
);
// Initially assume no mass transfer
nDmdtLfs_.insert
(
pair,
new volScalarField
(
IOobject
(
IOobject::groupName
(
"thermalPhaseChange:nucleation:dmdtLf",
pair.name()
),
this->mesh().time().timeName(),
this->mesh(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedScalar(dimEnergy/dimTime/dimVolume, 0)
)
);
}
}
@ -191,7 +209,22 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::dmdtf
const phasePairKey& key
) const
{
return BasePhaseSystem::dmdtf(key) + this->totalDmdtf(key);
tmp<volScalarField> tDmdtf = BasePhaseSystem::dmdtf(key);
const scalar dmdtfSign =
Pair<word>::compare(this->phasePairs_[key], key);
if (dmdtfs_.found(key))
{
tDmdtf.ref() += dmdtfSign**dmdtfs_[key];
}
if (nDmdtfs_.found(key))
{
tDmdtf.ref() += dmdtfSign**nDmdtfs_[key];
}
return tDmdtf;
}
@ -201,23 +234,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::dmdts() const
{
PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
{
const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
const volScalarField& dmdtf = *dmdtfIter();
addField(pair.phase1(), "dmdt", dmdtf, dmdts);
addField(pair.phase2(), "dmdt", - dmdtf, dmdts);
}
forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
{
const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
const volScalarField& nDmdtf = *nDmdtfIter();
addField(pair.phase1(), "dmdt", nDmdtf, dmdts);
addField(pair.phase2(), "dmdt", - nDmdtf, dmdts);
}
addDmdts(dmdts);
return dmdts;
}
@ -266,156 +283,100 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
phaseSystem::heatTransferTable& eqns = eqnsPtr();
forAllConstIter
(
typename BasePhaseSystem::heatTransferModelTable,
this->heatTransferModels_,
heatTransferModelIter
)
// Create temperatures at which to evaluate nucleation mass transfers
phaseSystem::dmdtfTable Tns;
forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
{
const phasePair& pair
const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
const saturationModel& satModel = this->saturation(nDmdtfIter.key());
Tns.insert(pair, satModel.Tsat(pair.phase1().thermo().p()).ptr());
}
// Mass transfer terms
if (volatile_ != "none")
{
{
phaseSystem::dmidtfTable dmidtfs;
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
{
const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
const volScalarField& dmdtf = *dmdtfIter();
dmidtfs.insert(pair, new HashPtrTable<volScalarField>());
dmidtfs[pair]->insert(volatile_, new volScalarField(dmdtf));
}
this->addDmidtHefs
(
dmidtfs,
Tfs_,
latentHeatScheme::upwind,
latentHeatTransfer::heat,
eqns
);
}
{
phaseSystem::dmidtfTable nDmidtfs;
forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
{
const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
const volScalarField& nDmdtf = *nDmdtfIter();
nDmidtfs.insert(pair, new HashPtrTable<volScalarField>());
nDmidtfs[pair]->insert(volatile_, new volScalarField(nDmdtf));
}
this->addDmidtHefs
(
nDmidtfs,
Tns,
0,
latentHeatScheme::upwind,
eqns
);
}
}
else
{
this->addDmdtHefs
(
this->phasePairs_[heatTransferModelIter.key()]
dmdtfs_,
Tfs_,
latentHeatScheme::upwind,
latentHeatTransfer::heat,
eqns
);
this->addDmdtHefs
(
nDmdtfs_,
Tns,
0,
latentHeatScheme::upwind,
eqns
);
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const volScalarField& Tf(*this->Tf_[pair]);
const volScalarField H1(heatTransferModelIter().first()->K());
const volScalarField H2(heatTransferModelIter().second()->K());
const volScalarField HEff(H1*H2/(H1 + H2));
*eqns[phase1.name()] +=
H1*(Tf - phase1.thermo().T())
- HEff*(phase2.thermo().T() - phase1.thermo().T());
*eqns[phase2.name()] +=
H2*(Tf - phase2.thermo().T())
- HEff*(phase1.thermo().T() - phase2.thermo().T());
}
forAllConstIter
(
phaseSystem::phaseModelList,
this->phases(),
phaseIter
)
// Lagging
{
const phaseModel& phase(phaseIter());
PtrList<volScalarField> dmdts(this->phases().size());
if (dmdt0s_.set(phase.index()))
addDmdts(dmdts);
forAllConstIter(phaseSystem::phaseModelList, this->phases(), phaseIter)
{
*eqns[phase.name()] +=
fvm::Sp(dmdt0s_[phase.index()],phase.thermo().he());
}
}
const phaseModel& phase = phaseIter();
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const volScalarField& he1(thermo1.he());
const volScalarField& he2(thermo2.he());
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
const volScalarField dmdtf(this->totalDmdtf(pair));
const volScalarField dmdtf21(posPart(dmdtf));
const volScalarField dmdtf12(negPart(dmdtf));
*eqns[phase1.name()] +=
- fvm::Sp(dmdtf21, he1)
+ dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -=
- fvm::Sp(dmdtf12, he2)
+ dmdtf12*K1 + dmdtf21*K2;
if (this->saturationModels_.found(phasePairIter.key()))
{
const volScalarField& Tf(*this->Tf_[pair]);
if (volatile_ != "none" && isA<rhoReactionThermo>(thermo1))
if (dmdt0s_.set(phase.index()))
{
const basicSpecieMixture& composition1 =
refCast<const rhoReactionThermo>(thermo1).composition();
*eqns[phase1.name()] +=
dmdtf21
*composition1.HE
*eqns[phase.name()] +=
fvm::Sp
(
composition1.species()[volatile_],
thermo1.p(),
Tf
dmdt0s_[phase.index()] - dmdts[phase.index()],
phase.thermo().he()
);
}
else
{
*eqns[phase1.name()] +=
dmdtf21*thermo1.he(thermo1.p(), Tf);
}
if (volatile_ != "none" && isA<rhoReactionThermo>(thermo2))
{
const basicSpecieMixture& composition2 =
refCast<const rhoReactionThermo>(thermo2).composition();
*eqns[phase2.name()] -=
dmdtf12
*composition2.HE
(
composition2.species()[volatile_],
thermo2.p(),
Tf
);
}
else
{
*eqns[phase2.name()] -=
dmdtf12*thermo2.he(thermo2.p(), Tf);
}
}
else
{
*eqns[phase1.name()] += dmdtf21*he2;
*eqns[phase2.name()] -= dmdtf12*he1;
}
if (this->nDmdtLfs_.found(phasePairIter.key()))
{
*eqns[phase1.name()] += negPart(*this->nDmdtLfs_[pair]);
*eqns[phase2.name()] -= posPart(*this->nDmdtLfs_[pair]);
}
if (phase1.thermo().he().member() == "e")
{
*eqns[phase1.name()] +=
phase1.thermo().p()*dmdtf/phase1.thermo().rho();
}
if (phase2.thermo().he().member() == "e")
{
*eqns[phase2.name()] -=
phase2.thermo().p()*dmdtf/phase2.thermo().rho();
}
}
@ -430,45 +391,44 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::specieTransfer() const
autoPtr<phaseSystem::specieTransferTable> eqnsPtr =
BasePhaseSystem::specieTransfer();
if (volatile_ == "none")
{
return eqnsPtr;
}
phaseSystem::specieTransferTable& eqns = eqnsPtr();
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
if (volatile_ != "none")
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
phaseSystem::dmidtfTable dmidtfs;
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
{
const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
const volScalarField& dmdtf = *dmdtfIter();
dmidtfs.insert(pair, new HashPtrTable<volScalarField>());
dmidtfs[pair]->insert(volatile_, new volScalarField(dmdtf));
}
this->addDmidtYf(dmidtfs, eqns);
}
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
if (phase1.pure() || phase2.pure())
{
FatalErrorInFunction
<< "Volatile specie was given, but at least one phase in pair "
<< pair << " is pure."
<< exit(FatalError);
phaseSystem::dmidtfTable nDmidtfs;
forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
{
const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
const volScalarField& nDmdtf = *nDmdtfIter();
nDmidtfs.insert(pair, new HashPtrTable<volScalarField>());
nDmidtfs[pair]->insert(volatile_, new volScalarField(nDmdtf));
}
this->addDmidtYf(nDmidtfs, eqns);
}
const volScalarField& Y1 = phase1.Y(volatile_);
const volScalarField& Y2 = phase2.Y(volatile_);
const volScalarField dmdtf(this->totalDmdtf(pair));
*eqns[Y1.name()] += dmdtf;
*eqns[Y2.name()] -= dmdtf;
}
else
{
this->addDmdtYfs(dmdtfs_, eqns);
this->addDmdtYfs(nDmdtfs_, eqns);
}
return eqnsPtr;
@ -481,23 +441,7 @@ correctContinuityError()
{
dmdt0s_ = PtrList<volScalarField>(this->phases().size());
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
{
const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
const volScalarField& dmdtf = *dmdtfIter();
addField(pair.phase1(), "dmdt", dmdtf, dmdt0s_);
addField(pair.phase2(), "dmdt", - dmdtf, dmdt0s_);
}
forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
{
const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
const volScalarField& nDmdtf = *nDmdtfIter();
addField(pair.phase1(), "dmdt", nDmdtf, dmdt0s_);
addField(pair.phase2(), "dmdt", - nDmdtf, dmdt0s_);
}
addDmdts(dmdt0s_);
BasePhaseSystem::correctContinuityError();
}
@ -512,97 +456,38 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
forAllConstIter
(
typename BasePhaseSystem::heatTransferModelTable,
this->heatTransferModels_,
heatTransferModelIter
saturationModelTable,
saturationModels_,
saturationModelIter
)
{
const phasePair& pair
(
this->phasePairs_[heatTransferModelIter.key()]
);
const phasePair& pair = this->phasePairs_[saturationModelIter.key()];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const volScalarField& T1(thermo1.T());
const volScalarField& T2(thermo2.T());
const volScalarField& p(thermo1.p());
volScalarField& dmdtf(*this->dmdtfs_[pair]);
volScalarField& Tf(*this->Tf_[pair]);
volScalarField dmdtfNew(dmdtf);
if (saturationModels_.found(heatTransferModelIter.key()))
// Interfacial mass transfer update
{
const phasePairKey& key = heatTransferModelIter.key();
const volScalarField Tsat(saturation(key).Tsat(thermo1.p()));
volScalarField& dmdtf(*this->dmdtfs_[pair]);
volScalarField& Tf(*this->Tfs_[pair]);
volScalarField ha1(thermo1.ha());
volScalarField haf1(thermo1.ha(p, Tsat));
const volScalarField Tsat(saturationModelIter()->Tsat(thermo1.p()));
if (volatile_ != "none" && isA<rhoReactionThermo>(thermo1))
{
const basicSpecieMixture& composition1 =
refCast<const rhoReactionThermo>(thermo1).composition();
ha1 =
composition1.Ha
(
composition1.species()[volatile_],
p,
T1
);
haf1 =
composition1.Ha
(
composition1.species()[volatile_],
p,
Tsat
);
}
volScalarField ha2(thermo2.ha());
volScalarField haf2(thermo2.ha(p, Tsat));
if (volatile_ != "none" && isA<rhoReactionThermo>(thermo2))
{
const basicSpecieMixture& composition2 =
refCast<const rhoReactionThermo>(thermo2).composition();
ha2 =
composition2.Ha
(
composition2.species()[volatile_],
p,
T2
);
haf2 =
composition2.Ha
(
composition2.species()[volatile_],
p,
Tsat
);
}
volScalarField L
const volScalarField L
(
(neg0(dmdtf)*haf2 + pos(dmdtf)*ha2)
- (pos0(dmdtf)*haf1 + neg(dmdtf)*ha1)
volatile_ != "none"
? this->Li(pair, volatile_, dmdtf, Tsat, latentHeatScheme::upwind)
: this->L(pair, dmdtf, Tsat, latentHeatScheme::upwind)
);
volScalarField H1(heatTransferModelIter().first()->K(0));
volScalarField H2(heatTransferModelIter().second()->K(0));
volScalarField H1(this->heatTransferModels_[pair].first()->K(0));
volScalarField H2(this->heatTransferModels_[pair].second()->K(0));
dmdtfNew = (H1*(Tsat - T1) + H2*(Tsat - T2))/L;
volScalarField dmdtfNew((H1*(Tsat - T1) + H2*(Tsat - T2))/L);
if (volatile_ != "none")
{
@ -611,36 +496,26 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
+ pos(dmdtfNew)*phase2.Y(volatile_);
}
H1 = heatTransferModelIter().first()->K();
H2 = heatTransferModelIter().second()->K();
H1 = this->heatTransferModels_[pair].first()->K();
H2 = this->heatTransferModels_[pair].second()->K();
// Limit the H[12] to avoid /0
H1.max(small);
H2.max(small);
Tf = (H1*T1 + H2*T2 + dmdtfNew*L)/(H1 + H2);
}
else
{
dmdtfNew = Zero;
volScalarField H1(heatTransferModelIter().first()->K());
volScalarField H2(heatTransferModelIter().second()->K());
Tf = (H1*T1 + H2*T2)/(H1 + H2);
}
Info<< Tf.name()
<< ": min = " << min(Tf.primitiveField())
<< ", mean = " << average(Tf.primitiveField())
<< ", max = " << max(Tf.primitiveField())
<< endl;
Info<< "Tf." << pair.name()
<< ": min = " << min(Tf.primitiveField())
<< ", mean = " << average(Tf.primitiveField())
<< ", max = " << max(Tf.primitiveField())
<< endl;
const scalar dmdtfRelax =
this->mesh().fieldRelaxationFactor(dmdtf.member());
scalar dmdtfRelax =
this->mesh().fieldRelaxationFactor(dmdtf.member());
dmdtf = (1 - dmdtfRelax)*dmdtf + dmdtfRelax*dmdtfNew;
dmdtf = (1 - dmdtfRelax)*dmdtf + dmdtfRelax*dmdtfNew;
if (saturationModels_.found(heatTransferModelIter.key()))
{
Info<< dmdtf.name()
<< ": min = " << min(dmdtf.primitiveField())
<< ", mean = " << average(dmdtf.primitiveField())
@ -649,88 +524,70 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
<< endl;
}
volScalarField& nDmdtf(*this->nDmdtfs_[pair]);
volScalarField& nDmdtLf(*this->nDmdtLfs_[pair]);
nDmdtf = Zero;
nDmdtLf = Zero;
bool wallBoilingActive = false;
forAllConstIter(phasePair, pair, iter)
// Nucleation mass transfer update
{
const phaseModel& phase = iter();
const phaseModel& otherPhase = iter.otherPhase();
volScalarField& nDmdtf(*this->nDmdtfs_[pair]);
if
(
phase.mesh().foundObject<volScalarField>
(
IOobject::groupName("alphat", phase.name())
) &&
saturationModels_.found(heatTransferModelIter.key())
)
nDmdtf = Zero;
bool wallBoilingActive = false;
forAllConstIter(phasePair, pair, pairIter)
{
const volScalarField& alphat =
phase.mesh().lookupObject<volScalarField>
(
IOobject::groupName("alphat", phase.name())
);
const phaseModel& phase = pairIter();
const fvPatchList& patches = this->mesh().boundary();
forAll(patches, patchi)
const word alphatName =
IOobject::groupName("alphat", phase.name());
if (phase.mesh().foundObject<volScalarField>(alphatName))
{
const fvPatch& currPatch = patches[patchi];
const volScalarField& alphat =
phase.mesh().lookupObject<volScalarField>(alphatName);
if
(
isA<alphatPhaseChangeWallFunction>
(
alphat.boundaryField()[patchi]
)
)
forAll(alphat.boundaryField(), patchi)
{
const alphatPhaseChangeWallFunction& PCpatch =
refCast<const alphatPhaseChangeWallFunction>
(
alphat.boundaryField()[patchi]
);
const fvPatchScalarField& alphatp =
alphat.boundaryField()[patchi];
phasePairKey key(phase.name(), otherPhase.name());
if (PCpatch.activePhasePair(key))
if (isA<alphatPhaseChangeWallFunction>(alphatp))
{
wallBoilingActive = true;
const alphatPhaseChangeWallFunction& alphatw =
refCast<const alphatPhaseChangeWallFunction>
(
alphatp
);
const scalarField& patchDmdtf = PCpatch.dmdtf(key);
const scalarField& patchDmdtLf =
PCpatch.dmdtLf(key);
const scalar sign
(
Pair<word>::compare(pair, key)
);
forAll(patchDmdtf, facei)
if (alphatw.activePhasePair(pair))
{
const label faceCelli =
currPatch.faceCells()[facei];
nDmdtf[faceCelli] -= sign*patchDmdtf[facei];
nDmdtLf[faceCelli] -= sign*patchDmdtLf[facei];
wallBoilingActive = true;
const scalarField& patchDmdtf =
alphatw.dmdtf(pair);
const scalar sign =
pairIter.index() == 0 ? +1 : -1;
forAll(patchDmdtf, facei)
{
const label celli =
alphatw.patch().faceCells()[facei];
nDmdtf[celli] -= sign*patchDmdtf[facei];
}
}
}
}
}
}
}
if (wallBoilingActive)
{
Info<< nDmdtf.name()
<< ": min = " << min(nDmdtf.primitiveField())
<< ", mean = " << average(nDmdtf.primitiveField())
<< ", max = " << max(nDmdtf.primitiveField())
<< ", integral = " << fvc::domainIntegrate(nDmdtf).value()
<< endl;
if (wallBoilingActive)
{
Info<< nDmdtf.name()
<< ": min = " << min(nDmdtf.primitiveField())
<< ", mean = " << average(nDmdtf.primitiveField())
<< ", max = " << max(nDmdtf.primitiveField())
<< ", integral = " << fvc::domainIntegrate(nDmdtf).value()
<< endl;
}
}
}
}

View File

@ -85,6 +85,10 @@ class ThermalPhaseChangePhaseSystem
phasePairKey::hash
> saturationModelTable;
using latentHeatScheme = typename BasePhaseSystem::latentHeatScheme;
using latentHeatTransfer = typename BasePhaseSystem::latentHeatTransfer;
// Private data
@ -97,6 +101,9 @@ class ThermalPhaseChangePhaseSystem
//- Mass transfer rates
phaseSystem::dmdtfTable dmdtfs_;
//- Interface temperatures
phaseSystem::dmdtfTable Tfs_;
//- Nucleate Mass transfer rates
phaseSystem::dmdtfTable nDmdtfs_;
@ -110,8 +117,8 @@ class ThermalPhaseChangePhaseSystem
// Private Member Functions
//- Return the total mass transfer rate for a pair
tmp<volScalarField> totalDmdtf(const phasePairKey& key) const;
//- Sum the mass transfer rates for each phase into the given list
void addDmdts(PtrList<volScalarField>&) const;
public:

View File

@ -35,266 +35,158 @@ template<class BasePhaseSystem>
void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
(
const phaseSystem::dmdtfTable& dmdtfs,
const phaseSystem::dmdtfTable& Tfs,
const latentHeatScheme scheme,
const latentHeatTransfer transfer,
phaseSystem::heatTransferTable& eqns
) const
{
// In the presence of thermally-coupled mass transfer, the relation between
// heat transfers across the interface between phases 1 and 2 is:
//
// Q1 + Q2 = mDot*L
// H1*(Tf - T1) + H2*(Tf - T1) = mDot*L
//
// Where Q1 and Q2 are the net transfer into phases 1 and 2 respectively,
// H1 and H2 are the heat transfer coefficients on either side, Tf is the
// temperature at the interface, mDot is the mass transfer rate from phase
// 2 to phase 1, and L is the latent heat of phase 2 minus phase 1.
//
// Rearranging for Tf:
//
// Tf = (H1*T1 + H2*T2 + mDot*L)/(H1 + H2)
//
// And for Q1 and Q2:
//
// Q1 = HEff*(T2 - T1) + H1/(H1 + H2)*mDot*L
// Q2 = HEff*(T1 - T2) + H2/(H1 + H2)*mDot*L
//
// The HEff terms are those implemented in the heatTransfer method. The
// latent heat mDot*L terms are added below (as well as other standard
// terms representing sensible energy and kinetic energy differences).
HeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefsWithoutL
(
dmdtfs,
Tfs,
scheme,
eqns
);
// Loop the pairs
forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
{
const phasePairKey& key = dmdtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
const volScalarField dmdtf21(posPart(dmdtf));
const volScalarField dmdtf12(negPart(dmdtf));
const volScalarField& Tf = *Tfs[key];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const volScalarField& he1(thermo1.he());
const volScalarField& he2(thermo2.he());
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
if (heatTransferModels_.found(key))
// Transfer coefficients
const volScalarField H1(heatTransferModels_[key].first()->K());
const volScalarField H2(heatTransferModels_[key].second()->K());
const volScalarField H1Fac(H1/(H1 + H2));
const volScalarField HEff(H1Fac*H2);
// Latent heat contribution
switch (transfer)
{
// Assume a thermally-coupled mass transfer process. Calculate
// heat transfers by evaluating properties at the interface
// temperature...
case latentHeatTransfer::heat:
{
*eqns[phase1.name()] +=
- HEff*(thermo2.T() - thermo1.T()) + H1*(Tf - thermo1.T());
// Transfer of energy from the interface into the bulk
const volScalarField& Tf(*Tf_[key]);
const volScalarField hef1(thermo1.he(thermo1.p(), Tf));
const volScalarField hef2(thermo2.he(thermo2.p(), Tf));
*eqns[phase1.name()] += dmdtf*hef1;
*eqns[phase2.name()] -= dmdtf*hef2;
*eqns[phase2.name()] +=
- HEff*(thermo1.T() - thermo2.T()) + H2*(Tf - thermo2.T());
// Latent heat contribution
const volScalarField L(hef2 + thermo2.hc() - hef1 - thermo1.hc());
const volScalarField H1(heatTransferModels_[key].first()->K());
const volScalarField H2(heatTransferModels_[key].second()->K());
const volScalarField H1Fac(H1/(H1 + H2));
*eqns[phase1.name()] += H1Fac*dmdtf*L;
*eqns[phase2.name()] += (1 - H1Fac)*dmdtf*L;
break;
}
case latentHeatTransfer::mass:
{
const volScalarField L(this->L(pair, dmdtf, Tf, scheme));
// Transfer of kinetic energy
*eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2;
}
else
{
// In the absence of a heat transfer model, assume a
// non-thermally-coupled mass transfer process; i.e., no phase
// change and therefore no interface state or latent heat...
*eqns[phase1.name()] += H1Fac*dmdtf*L;
*eqns[phase2.name()] += (1 - H1Fac)*dmdtf*L;
// Transfer of energy from bulk to bulk
*eqns[phase1.name()] += dmdtf21*he2 + fvm::Sp(dmdtf12, he1);
*eqns[phase2.name()] -= dmdtf12*he1 + fvm::Sp(dmdtf21, he2);
// Transfer of kinetic energy
*eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2;
break;
}
}
}
}
template<class BasePhaseSystem>
void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs
(
const phaseSystem::dmidtfTable& dmidtfs,
const phaseSystem::dmdtfTable& Tfs,
const latentHeatScheme scheme,
const latentHeatTransfer transfer,
phaseSystem::heatTransferTable& eqns
) const
{
// See the addDmdtHe method for a description of the latent heat terms
HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefsWithoutL
(
dmidtfs,
Tfs,
scheme,
eqns
);
// Loop the pairs
forAllConstIter(phaseSystem::dmidtfTable, dmidtfs, dmidtfIter)
{
const phasePairKey& key = dmidtfIter.key();
const phasePair& pair(this->phasePairs_[key]);
const volScalarField& Tf = *Tfs[key];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const rhoThermo& thermo1 = phase1.thermo();
const rhoThermo& thermo2 = phase2.thermo();
const volScalarField& he1(thermo1.he());
const volScalarField& he2(thermo2.he());
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
if (heatTransferModels_.found(key))
// Transfer coefficients
const volScalarField H1(heatTransferModels_[key].first()->K());
const volScalarField H2(heatTransferModels_[key].second()->K());
const volScalarField H1Fac(H1/(H1 + H2));
const volScalarField HEff(H1Fac*H2);
// Loop the species
forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
{
// Assume a thermally-coupled mass transfer process. Calculate
// heat transfers by evaluating properties at the interface
// temperature...
const word& specie = dmidtfJter.key();
// Transfer coefficients
const volScalarField H1(heatTransferModels_[key].first()->K());
const volScalarField H2(heatTransferModels_[key].second()->K());
const volScalarField H1Fac(H1/(H1 + H2));
// Interface properties
const volScalarField& Tf(*Tf_[key]);
const volScalarField hef1(thermo1.he(thermo1.p(), Tf));
const volScalarField hef2(thermo2.he(thermo2.p(), Tf));
const volScalarField hc1(thermo1.hc());
const volScalarField hc2(thermo2.hc());
// Loop the species
forAllConstIter
// Mass transfer rates
const volScalarField dmidtf
(
HashPtrTable<volScalarField>,
*dmidtfIter(),
dmidtfJter
)
Pair<word>::compare(pair, key)**dmidtfJter()
);
// Latent heat contribution
switch (transfer)
{
const word& member = dmidtfJter.key();
const volScalarField dmidtf
(
Pair<word>::compare(pair, key)**dmidtfJter()
);
const volScalarField dmidtf21(posPart(dmidtf));
const volScalarField dmidtf12(negPart(dmidtf));
// Create the interface energies for the transferring specie
volScalarField hefi1(hef1), hci1(hc1);
if (isA<rhoReactionThermo>(thermo1))
case latentHeatTransfer::heat:
{
const basicSpecieMixture& composition1 =
refCast<const rhoReactionThermo>(thermo1).composition();
hefi1 =
composition1.HE
(
composition1.species()[member],
thermo1.p(),
Tf
);
hci1 =
dimensionedScalar
(
dimEnergy/dimMass,
composition1.Hf(composition1.species()[member])
);
// Do nothing. This term is handled outside the specie loop.
break;
}
volScalarField hefi2(hef2), hci2(hc2);
if (isA<rhoReactionThermo>(thermo2))
case latentHeatTransfer::mass:
{
const basicSpecieMixture& composition2 =
refCast<const rhoReactionThermo>(thermo2).composition();
hefi2 =
composition2.HE
(
composition2.species()[member],
thermo2.p(),
Tf
);
hci2 =
dimensionedScalar
(
dimEnergy/dimMass,
composition2.Hf(composition2.species()[member])
);
const volScalarField Li
(
this->Li(pair, specie, dmidtf, Tf, scheme)
);
*eqns[phase1.name()] += H1Fac*dmidtf*Li;
*eqns[phase2.name()] += (1 - H1Fac)*dmidtf*Li;
break;
}
// Create the latent heat for the transferring specie
const volScalarField Li(hefi2 + hci2 - hefi1 - hci1);
// Transfer of energy from the interface into the bulk
*eqns[phase1.name()] += dmidtf*hefi1;
*eqns[phase2.name()] -= dmidtf*hefi2;
// Latent heat contribution
*eqns[phase1.name()] += H1Fac*dmidtf*Li;
*eqns[phase2.name()] += (1 - H1Fac)*dmidtf*Li;
// Transfer of kinetic energy
*eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1;
*eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2;
}
}
else
// Latent heat contribution
switch (transfer)
{
// In the absence of a heat transfer model, assume a
// non-thermally-coupled mass transfer process; i.e., no phase
// change and therefore no interface state or latent heat...
// Loop the species
forAllConstIter
(
HashPtrTable<volScalarField>,
*dmidtfIter(),
dmidtfJter
)
case latentHeatTransfer::heat:
{
const word& member = dmidtfJter.key();
*eqns[phase1.name()] +=
- HEff*(thermo2.T() - thermo1.T()) + H1*(Tf - thermo1.T());
const volScalarField dmidtf
(
Pair<word>::compare(pair, key)**dmidtfJter()
);
const volScalarField dmidtf21(posPart(dmidtf));
const volScalarField dmidtf12(negPart(dmidtf));
*eqns[phase2.name()] +=
- HEff*(thermo1.T() - thermo2.T()) + H2*(Tf - thermo2.T());
// Create the energies for the transferring specie
volScalarField hei1(he1);
if (isA<rhoReactionThermo>(thermo1))
{
const basicSpecieMixture& composition1 =
refCast<const rhoReactionThermo>(thermo1).composition();
hei1 =
composition1.HE
(
composition1.species()[member],
thermo1.p(),
thermo1.T()
);
}
volScalarField hei2(he2);
if (isA<rhoReactionThermo>(thermo2))
{
const basicSpecieMixture& composition2 =
refCast<const rhoReactionThermo>(thermo2).composition();
hei2 =
composition2.HE
(
composition2.species()[member],
thermo2.p(),
thermo2.T()
);
}
break;
}
case latentHeatTransfer::mass:
{
// Do nothing. This term is handled inside the specie loop.
// Transfer of energy from bulk to bulk
*eqns[phase1.name()] += dmidtf21*hei2 + fvm::Sp(dmidtf12, he1);
*eqns[phase2.name()] -= dmidtf12*hei1 + fvm::Sp(dmidtf21, he2);
// Transfer of kinetic energy
*eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1;
*eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2;
break;
}
}
}
@ -310,7 +202,7 @@ TwoResistanceHeatTransferPhaseSystem
const fvMesh& mesh
)
:
BasePhaseSystem(mesh)
HeatTransferPhaseSystem<BasePhaseSystem>(mesh)
{
this->generatePairsAndSubModels
(
@ -344,37 +236,6 @@ TwoResistanceHeatTransferPhaseSystem
<< exit(FatalError);
}
}
// Calculate initial Tf-s as the mean between the two phases
forAllConstIter
(
heatTransferModelTable,
heatTransferModels_,
heatTransferModelIter
)
{
const phasePair& pair = this->phasePairs_[heatTransferModelIter.key()];
const volScalarField& T1(pair.phase1().thermo().T());
const volScalarField& T2(pair.phase2().thermo().T());
Tf_.insert
(
pair,
new volScalarField
(
IOobject
(
IOobject::groupName("Tf", pair.name()),
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
(T1 + T2)/2
)
);
}
}
@ -411,32 +272,6 @@ heatTransfer() const
);
}
// In the absence of mass transfer, the relation between heat transfers
// across the interface between phases 1 and 2 is:
//
// Q1 + Q2 = 0
// H1*(Tf - T1) + H2*(Tf - T1) = 0
//
// Where Q1 and Q2 are the net transfer into phases 1 and 2 respectively,
// H1 and H2 are the heat transfer coefficients on either side, and Tf is
// the temperature at the interface.
//
// Rearranging for Tf:
//
// Tf = (H1*T1 + H2*T2)/(H1 + H2)
//
// Which gives the following for Q1 and Q2:
//
// Q1 = H1*H2/(H1 + H2)*(T2 - T1)
// Q2 = H1*H2/(H1 + H2)*(T1 - T2)
//
// The transfer is therefore the same as for the single-resistance system
// with a single effective heat transfer coefficient:
//
// HEff = H1*H2/(H1 + H2)
//
// This is implemented below
forAllConstIter
(
heatTransferModelTable,
@ -444,24 +279,19 @@ heatTransfer() const
heatTransferModelIter
)
{
const phasePair& pair
(
this->phasePairs_[heatTransferModelIter.key()]
);
const phasePair& pair = this->phasePairs_[heatTransferModelIter.key()];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const volScalarField& he1(phase1.thermo().he());
const volScalarField& he2(phase2.thermo().he());
const volScalarField& he1 = phase1.thermo().he();
const volScalarField& he2 = phase2.thermo().he();
const volScalarField Cpv1(phase1.thermo().Cpv());
const volScalarField Cpv2(phase2.thermo().Cpv());
// Heat transfer between phases in the absence of mass transfer
const volScalarField H1(heatTransferModelIter().first()->K());
const volScalarField H2(heatTransferModelIter().second()->K());
const volScalarField HEff(H1*H2/(H1 + H2));
*eqns[phase1.name()] +=
HEff*(phase2.thermo().T() - phase1.thermo().T())
+ H1/Cpv1*he1 - fvm::Sp(H1/Cpv1, he1);
@ -484,42 +314,6 @@ correctEnergyTransport()
}
template<class BasePhaseSystem>
void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::
correctInterfaceThermo()
{
forAllConstIter
(
heatTransferModelTable,
heatTransferModels_,
heatTransferModelIter
)
{
const phasePair& pair = this->phasePairs_[heatTransferModelIter.key()];
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const volScalarField& T1(phase1.thermo().T());
const volScalarField& T2(phase2.thermo().T());
const volScalarField H1(heatTransferModels_[pair].first()->K());
const volScalarField H2(heatTransferModels_[pair].second()->K());
const dimensionedScalar HSmall(H1.dimensions(), small);
volScalarField& Tf = *Tf_[pair];
Tf = (H1*T1 + H2*T2)/max(H1 + H2, HSmall);
Info<< "Tf." << pair.name()
<< ": min = " << min(Tf.primitiveField())
<< ", mean = " << average(Tf.primitiveField())
<< ", max = " << max(Tf.primitiveField())
<< endl;
}
}
template<class BasePhaseSystem>
bool Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::read()
{

View File

@ -42,7 +42,8 @@ SourceFiles
#ifndef TwoResistanceHeatTransferPhaseSystem_H
#define TwoResistanceHeatTransferPhaseSystem_H
#include "phaseSystem.H"
#include "twoResistanceHeatTransferPhaseSystem.H"
#include "HeatTransferPhaseSystem.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -62,7 +63,8 @@ class heatTransferModel;
template<class BasePhaseSystem>
class TwoResistanceHeatTransferPhaseSystem
:
public BasePhaseSystem
public twoResistanceHeatTransferPhaseSystem,
public HeatTransferPhaseSystem<BasePhaseSystem>
{
protected:
@ -75,12 +77,13 @@ protected:
phasePairKey::hash
> heatTransferModelTable;
typedef
typename HeatTransferPhaseSystem<BasePhaseSystem>::latentHeatScheme
latentHeatScheme;
// Protected data
//- Interface temperatures
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash> Tf_;
// Sub Models
//- Heat transfer models
@ -90,16 +93,32 @@ protected:
// Protected Member Functions
//- Add energy transfer terms which result from bulk mass transfers
// and phase changes
using HeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs;
//- Add energy transfer terms which result from bulk phase changes
// that are coupled to the two-resistance heat transfer model
void addDmdtHefs
(
const phaseSystem::dmdtfTable& dmdtfs,
const phaseSystem::dmdtfTable& Tfs,
const latentHeatScheme scheme,
const latentHeatTransfer transfer,
phaseSystem::heatTransferTable& eqns
) const;
//- Add energy transfer terms which result from specie mass transfers
void addDmidtHef
// and phase changes
using HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs;
//- Add energy transfer terms which result from specie phase changes
// that are coupled to the two-resistance heat transfer model
void addDmidtHefs
(
const phaseSystem::dmidtfTable& dmidtfs,
const phaseSystem::dmdtfTable& Tfs,
const latentHeatScheme scheme,
const latentHeatTransfer transfer,
phaseSystem::heatTransferTable& eqns
) const;
@ -125,7 +144,7 @@ public:
virtual void correctEnergyTransport();
//- Correct the interface thermodynamics
virtual void correctInterfaceThermo();
virtual void correctInterfaceThermo() = 0;
//- Read base phaseProperties dictionary
virtual bool read();

View File

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2020 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "twoResistanceHeatTransferPhaseSystem.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
template<>
const char* NamedEnum
<
Foam::twoResistanceHeatTransferPhaseSystem::latentHeatTransfer,
2
>::names[] = {"heat", "mass"};
}
const Foam::NamedEnum
<
Foam::twoResistanceHeatTransferPhaseSystem::latentHeatTransfer,
2
>
Foam::twoResistanceHeatTransferPhaseSystem::latentHeatTransferNames_;
// ************************************************************************* //

View File

@ -0,0 +1,70 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2020 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::twoResistanceHeatTransferPhaseSystem
SourceFiles
twoResistanceHeatTransferPhaseSystem.C
\*---------------------------------------------------------------------------*/
#ifndef twoResistanceHeatTransferPhaseSystem_H
#define twoResistanceHeatTransferPhaseSystem_H
#include "NamedEnum.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class twoResistanceHeatTransferPhaseSystem Declaration
\*---------------------------------------------------------------------------*/
class twoResistanceHeatTransferPhaseSystem
{
public:
//- Enumeration for the form of the latent heat transfer
enum class latentHeatTransfer
{
heat,
mass
};
//- Names of the forms of the latent heat transfer
static const NamedEnum<latentHeatTransfer, 2> latentHeatTransferNames_;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -72,9 +72,9 @@ plot \
'postProcessing_'.gY0.'_'.lY0.'/plot/0/T.liquid' w l lc 2 t 'Liquid', \
for [gY in gYs] for [lY in lYs] \
'postProcessing_'.gY.'_'.lY.'/plot/0/T.liquid' w l lc 2 notitle, \
'postProcessing_'.gY0.'_'.lY0.'/plot/0/Tf.gasAndLiquid' w l lc 3 t 'Interface', \
'postProcessing_'.gY0.'_'.lY0.'/plot/0/interfaceCompositionPhaseChange:Tf.gasAndLiquid' w l lc 3 t 'Interface', \
for [gY in gYs] for [lY in lYs] \
'postProcessing_'.gY.'_'.lY.'/plot/0/Tf.gasAndLiquid' w l lc 3 notitle
'postProcessing_'.gY.'_'.lY.'/plot/0/interfaceCompositionPhaseChange:Tf.gasAndLiquid' w l lc 3 notitle
set ytics nomirror
set y2tics

View File

@ -25,7 +25,7 @@ fields
p
T.gas
T.liquid
Tf.gasAndLiquid
interfaceCompositionPhaseChange:Tf.gasAndLiquid
dMass.gas
dMass.liquid
dMass

View File

@ -108,9 +108,9 @@ plot \
'postProcessing_'.gY0.'_'.gHe0.'_'.lHe0.'/plot/0/T.liquid' w l lc 2 t 'Liquid', \
for [gY in gYs] for [ gHe in gHes ] for [ lHe in lHes ] \
'postProcessing_'.gY.'_'.gHe.'_'.lHe.'/plot/0/T.liquid' w l lc 2 notitle, \
'postProcessing_'.gY0.'_'.gHe0.'_'.lHe0.'/plot/0/Tf.gasAndLiquid' w l lc 3 t 'Interface', \
'postProcessing_'.gY0.'_'.gHe0.'_'.lHe0.'/plot/0/interfaceCompositionPhaseChange:Tf.gasAndLiquid' w l lc 3 t 'Interface', \
for [gY in gYs] for [ gHe in gHes ] for [ lHe in lHes ] \
'postProcessing_'.gY.'_'.gHe.'_'.lHe.'/plot/0/Tf.gasAndLiquid' w l lc 3 notitle
'postProcessing_'.gY.'_'.gHe.'_'.lHe.'/plot/0/interfaceCompositionPhaseChange:Tf.gasAndLiquid' w l lc 3 notitle
set ytics nomirror
set y2tics

View File

@ -22,7 +22,7 @@ fields
p
T.gas
T.liquid
Tf.gasAndLiquid
interfaceCompositionPhaseChange:Tf.gasAndLiquid
dMass.gas
dMass.liquid
dMass

View File

@ -81,9 +81,9 @@ plot \
'postProcessing_'.gHe0.'_'.lHe0.'/plot/0/T.water' w l lc 2 t 'Water', \
for [gHe in gHes] for [lHe in lHes] \
'postProcessing_'.gHe.'_'.lHe.'/plot/0/T.water' w l lc 2 notitle, \
'postProcessing_'.gHe0.'_'.lHe0.'/plot/0/Tf.steamAndWater' w l lc 3 t 'Interface', \
'postProcessing_'.gHe0.'_'.lHe0.'/plot/0/thermalPhaseChange:Tf.steamAndWater' w l lc 3 t 'Interface', \
for [gHe in gHes] for [lHe in lHes] \
'postProcessing_'.gHe.'_'.lHe.'/plot/0/Tf.steamAndWater' w l lc 3 notitle
'postProcessing_'.gHe.'_'.lHe.'/plot/0/thermalPhaseChange:Tf.steamAndWater' w l lc 3 notitle
set ytics nomirror
set y2tics

View File

@ -22,7 +22,7 @@ fields
p
T.steam
T.water
Tf.steamAndWater
thermalPhaseChange:Tf.steamAndWater
dMass.steam
dMass.water
dMass