Commit 2518e125 authored by Henry Weller's avatar Henry Weller
Browse files

reactingEulerFoam: Support compressibility and mass-transfer independently

Now combinations of incompressible, compressible phases with or without
mass-transfer are supported efficiently.
parent 3487bd49
......@@ -151,6 +151,16 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
bool Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::transfersMass
(
const phaseModel& phase
) const
{
return true;
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
......
......@@ -120,6 +120,9 @@ public:
// Member Functions
//- Return true if there is mass transfer for phase
virtual bool transfersMass(const phaseModel& phase) const;
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
......
......@@ -92,20 +92,7 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
const Foam::phaseModel& phase
) const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
IOobject::groupName("dmdt", phase.name()),
this->mesh().time().timeName(),
this->mesh().time()
),
this->mesh(),
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
return tmp<volScalarField>(NULL);
}
......
......@@ -25,6 +25,7 @@ License
#include "ThermalPhaseChangePhaseSystem.H"
#include "fvCFD.H"
#include "alphatPhaseChangeWallFunctionFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -39,7 +40,61 @@ ThermalPhaseChangePhaseSystem
volatile_(this->lookup("volatile")),
saturationModel_(saturationModel::New(this->subDict("saturationModel"))),
massTransfer_(this->lookup("massTransfer"))
{}
{
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
// Initially assume no mass transfer
iDmdt_.insert
(
pair,
new volScalarField
(
IOobject
(
IOobject::groupName("iDmdt", pair.name()),
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
// Initially assume no mass transfer
wDmdt_.insert
(
pair,
new volScalarField
(
IOobject
(
IOobject::groupName("wDmdt", pair.name()),
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
......@@ -52,6 +107,14 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
const Foam::saturationModel&
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::saturation() const
{
return saturationModel_();
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::massTransferTable>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
......@@ -121,6 +184,9 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
template<class BasePhaseSystem>
void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
{
typedef compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
alphatPhaseChangeWallFunction;
BasePhaseSystem::correctThermo();
forAllConstIter
......@@ -147,6 +213,8 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
const volScalarField& he2(phase2.thermo().he());
volScalarField& dmdt(*this->dmdt_[pair]);
volScalarField& iDmdt(*this->iDmdt_[pair]);
volScalarField& wDmdt(*this->wDmdt_[pair]);
volScalarField& Tf = *this->Tf_[pair];
......@@ -167,25 +235,28 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
Tf = saturationModel_->Tsat(phase1.thermo().p());
dmdt =
(H1*(Tf - T1) + H2*(Tf - T2))
scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
iDmdt =
(1 - iDmdtRelax)*iDmdt
+ iDmdtRelax*(H1*(Tf - T1) + H2*(Tf - T2))
/min
(
(pos(dmdt)*he2 + neg(dmdt)*hef2)
- (neg(dmdt)*he1 + pos(dmdt)*hef1),
(pos(iDmdt)*he2 + neg(iDmdt)*hef2)
- (neg(iDmdt)*he1 + pos(iDmdt)*hef1),
0.3*mag(hef2 - hef1)
);
Info<< "dmdt." << pair.name()
<< ": min = " << min(dmdt.internalField())
<< ", mean = " << average(dmdt.internalField())
<< ", max = " << max(dmdt.internalField())
<< ", integral = " << fvc::domainIntegrate(dmdt).value()
Info<< "iDmdt." << pair.name()
<< ": min = " << min(iDmdt.internalField())
<< ", mean = " << average(iDmdt.internalField())
<< ", max = " << max(iDmdt.internalField())
<< ", integral = " << fvc::domainIntegrate(iDmdt).value()
<< endl;
}
else
{
dmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
iDmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
}
volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
......@@ -200,12 +271,13 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
volScalarField mDotL
(
dmdt*
iDmdt*
(
(pos(dmdt)*he2 + neg(dmdt)*hef2)
- (neg(dmdt)*he1 + pos(dmdt)*hef1)
(pos(iDmdt)*he2 + neg(iDmdt)*hef2)
- (neg(iDmdt)*he1 + pos(iDmdt)*hef1)
)
);
Tf = (H1*T1 + H2*T2 + mDotL)/(H1 + H2);
Info<< "Tf." << pair.name()
......@@ -213,6 +285,68 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
<< ", mean = " << average(Tf.internalField())
<< ", max = " << max(Tf.internalField())
<< endl;
// Accumulate dmdt contributions from boundaries
if
(
phase2.mesh().foundObject<volScalarField>
(
"alphat." + phase2.name()
)
)
{
scalar wDmdtRelax(this->mesh().fieldRelaxationFactor("wDmdt"));
wDmdt *= (1 - wDmdtRelax);
const volScalarField& alphat =
phase2.mesh().lookupObject<volScalarField>
(
"alphat." + phase2.name()
);
const fvPatchList& patches = this->mesh().boundary();
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if
(
isA<alphatPhaseChangeWallFunction>
(
alphat.boundaryField()[patchi]
)
)
{
const scalarField& patchDmdt =
refCast<const alphatPhaseChangeWallFunction>
(
alphat.boundaryField()[patchi]
).dmdt();
forAll(patchDmdt,facei)
{
label faceCelli = currPatch.faceCells()[facei];
wDmdt[faceCelli] += wDmdtRelax*patchDmdt[facei];
}
}
}
Info<< "wDmdt." << pair.name()
<< ": min = " << min(wDmdt.internalField())
<< ", mean = " << average(wDmdt.internalField())
<< ", max = " << max(wDmdt.internalField())
<< ", integral = " << fvc::domainIntegrate(wDmdt).value()
<< endl;
}
dmdt = wDmdt + iDmdt;
Info<< "dmdt." << pair.name()
<< ": min = " << min(dmdt.internalField())
<< ", mean = " << average(dmdt.internalField())
<< ", max = " << max(dmdt.internalField())
<< ", integral = " << fvc::domainIntegrate(dmdt).value()
<< endl;
}
}
......
......@@ -73,6 +73,14 @@ protected:
// Mass transfer enabled
Switch massTransfer_;
//- Interfacial Mass transfer rate
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
iDmdt_;
//- Wall Mass transfer rate
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
wDmdt_;
public:
......@@ -88,6 +96,9 @@ public:
// Member Functions
//- Return the saturationModel
const saturationModel& saturation() const;
//- Return the mass transfer matrices
virtual autoPtr<phaseSystem::massTransferTable> massTransfer() const;
......
......@@ -37,17 +37,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
)
:
BasePhaseModel(fluid, phaseName, index),
divU_
(
IOobject
(
IOobject::groupName("divU", this->name()),
fluid.mesh().time().timeName(),
fluid.mesh()
),
fluid.mesh(),
dimensionedScalar("divU", dimless/dimTime, 0)
),
divU_(NULL),
K_
(
IOobject
......@@ -105,7 +95,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
volScalarField& he = this->thermo_->he();
return
tmp<fvScalarMatrix> tEEqn
(
fvm::ddt(alpha, this->rho(), he) + fvm::div(alphaRhoPhi, he)
- fvm::Sp(contErr, he)
......@@ -119,16 +109,23 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
*fvc::interpolate(alphaEff),
he
)
+ (
he.name() == this->thermo_->phasePropertyName("e")
? fvc::ddt(alpha)*this->thermo().p()
+ fvc::div(alphaPhi, this->thermo().p())
: -alpha*this->fluid().dpdt()
)
==
this->Sh()
);
// Add the appropriate pressure-work term
if (he.name() == this->thermo_->phasePropertyName("e"))
{
tEEqn() +=
fvc::ddt(alpha)*this->thermo().p()
+ fvc::div(alphaPhi, this->thermo().p());
}
else if (this->thermo_->dpdt())
{
tEEqn() -= alpha*this->fluid().dpdt();
}
return tEEqn;
}
......@@ -140,7 +137,7 @@ bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
template<class BasePhaseModel>
const Foam::volScalarField&
const Foam::tmp<Foam::volScalarField>&
Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
{
return divU_;
......@@ -149,7 +146,10 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
template<class BasePhaseModel>
void
Foam::AnisothermalPhaseModel<BasePhaseModel>::divU(const volScalarField& divU)
Foam::AnisothermalPhaseModel<BasePhaseModel>::divU
(
const tmp<volScalarField>& divU
)
{
divU_ = divU;
}
......
......@@ -55,7 +55,7 @@ class AnisothermalPhaseModel
// Private data
//- Dilatation rate
volScalarField divU_;
tmp<volScalarField> divU_;
//- Kinetic energy
volScalarField K_;
......@@ -95,10 +95,10 @@ public:
virtual bool compressible() const;
//- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual const volScalarField& divU() const;
virtual const tmp<volScalarField>& divU() const;
//- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual void divU(const volScalarField& divU);
virtual void divU(const tmp<volScalarField>& divU);
//- Return the phase kinetic energy
virtual const volScalarField& K() const;
......
......@@ -165,16 +165,17 @@ bool Foam::phaseModel::compressible() const
}
const Foam::volScalarField& Foam::phaseModel::divU() const
const Foam::tmp<Foam::volScalarField>& Foam::phaseModel::divU() const
{
notImplemented("Foam::phaseModel::divU()");
return volScalarField::null();
static tmp<Foam::volScalarField> divU_(NULL);
return divU_;
}
void Foam::phaseModel::divU(const volScalarField& divU)
void Foam::phaseModel::divU(const tmp<volScalarField>& divU)
{
WarningIn("phaseModel::divU(const volScalarField& divU)")
WarningIn("phaseModel::divU(const tmp<volScalarField>& divU)")
<< "Attempt to set the dilatation rate of an incompressible phase"
<< endl;
}
......
......@@ -216,10 +216,10 @@ public:
virtual bool compressible() const;
//- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual const volScalarField& divU() const;
virtual const tmp<volScalarField>& divU() const;
//- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual void divU(const volScalarField& divU);
virtual void divU(const tmp<volScalarField>& divU);
//- Return the phase kinetic energy
virtual const volScalarField& K() const;
......
......@@ -120,6 +120,19 @@ constTransport
>
> constRefRhoConstEThermoPhysics;
typedef
constTransport
<
species::thermo
<
hRefConstThermo
<
rhoConst<specie>
>,
sensibleEnthalpy
>
> constRefRhoConstHThermoPhysics;
// pureMixture, sensibleEnthalpy:
......@@ -229,6 +242,18 @@ makeReactionMixtureThermo
);
// multiComponentMixture, sensibleEnthalpy:
makeReactionMixtureThermo
(
rhoThermo,
rhoReactionThermo,
heRhoThermo,
multiComponentMixture,
constRefRhoConstHThermoPhysics
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
......
......@@ -204,7 +204,7 @@ void Foam::multiphaseSystem::solveAlphas()
mesh_
),
mesh_,
dimensionedScalar("Sp", phase.divU().dimensions(), 0.0)
dimensionedScalar("Sp", divU.dimensions(), 0.0)
);
volScalarField::DimensionedInternalField Su
......@@ -220,9 +220,9 @@ void Foam::multiphaseSystem::solveAlphas()
divU*min(alpha, scalar(1))
);
if (phase.compressible())
if (phase.divU().valid())
{
const scalarField& dgdt = phase.divU();
const scalarField& dgdt = phase.divU()();
forAll(dgdt, celli)
{
......@@ -247,9 +247,9 @@ void Foam::multiphaseSystem::solveAlphas()
if (&phase2 == &phase) continue;
if (phase2.compressible())
if (phase2.divU().valid())
{
const scalarField& dgdt2 = phase2.divU();
const scalarField& dgdt2 = phase2.divU()();
forAll(dgdt2, celli)
{
......@@ -515,6 +515,12 @@ Foam::multiphaseSystem::~multiphaseSystem()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::multiphaseSystem::transfersMass(const phaseModel& phase) const
{
return false;
}
Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
(
const phaseModel& phase1
......
......@@ -175,6 +175,9 @@ public:
const PtrList<volScalarField>& rAUs
) const = 0;
//- Return true if there is mass transfer for phase
virtual bool transfersMass(const phaseModel& phase) const;
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const = 0;
......
......@@ -311,7 +311,7 @@ while (pimple.correct())
phasei,
(
(
phase.continuityError() - fluid.dmdt(phase)
phase.continuityError()
- fvc::Sp
(
fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
......@@ -340,7 +340,7 @@ while (pimple.correct())
phasei,