Commit 24c7a739 authored by Henry Weller's avatar Henry Weller
Browse files

reactingTwoPhaseEulerFoam: Improved support for boiling/condensation

Includes many contributions from Juho Peltola
parent c9c4e47a
......@@ -47,3 +47,9 @@
}
fluid.correctThermo();
Info<< " phase1.thermo().T(): " << min(phase1.thermo().T()).value()
<< " - " << max(phase1.thermo().T()).value() << endl;
Info<< " phase2.thermo().T(): " << min(phase2.thermo().T()).value()
<< " - " << max(phase2.thermo().T()).value() << endl;
......@@ -60,13 +60,13 @@ Foam::heatTransferModels::RanzMarshall::~RanzMarshall()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
Foam::heatTransferModels::RanzMarshall::K() const
Foam::heatTransferModels::RanzMarshall::K(const scalar residualAlpha) const
{
volScalarField Nu(scalar(2) + 0.6*sqrt(pair_.Re())*cbrt(pair_.Pr()));
return
6.0
*max(pair_.dispersed(), residualAlpha_)
*max(pair_.dispersed(), residualAlpha)
*pair_.continuous().kappa()
*Nu
/sqr(pair_.dispersed().d());
......
......@@ -79,7 +79,7 @@ public:
// Member Functions
//- The heat transfer function K used in the enthalpy equation
tmp<volScalarField> K() const;
tmp<volScalarField> K(const scalar residualAlpha) const;
};
......
......@@ -65,4 +65,13 @@ Foam::heatTransferModel::~heatTransferModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
Foam::heatTransferModel::K() const
{
return K(residualAlpha_.value());
}
// ************************************************************************* //
......@@ -117,7 +117,13 @@ public:
//- The heat transfer function K used in the enthalpy equation
// ddt(alpha1*rho1*ha) + ... = ... K*(Ta - Tb)
// ddt(alpha2*rho2*hb) + ... = ... K*(Tb - Ta)
virtual tmp<volScalarField> K() const = 0;
tmp<volScalarField> K() const;
//- The heat transfer function K used in the enthalpy equation
// ddt(alpha1*rho1*ha) + ... = ... K*(Ta - Tb)
// ddt(alpha2*rho2*hb) + ... = ... K*(Tb - Ta)
// with a specified residual volume fraction
virtual tmp<volScalarField> K(const scalar residualAlpha) const = 0;
};
......
......@@ -65,11 +65,14 @@ Foam::heatTransferModels::sphericalHeatTransfer::~sphericalHeatTransfer()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
Foam::heatTransferModels::sphericalHeatTransfer::K() const
Foam::heatTransferModels::sphericalHeatTransfer::K
(
const scalar residualAlpha
) const
{
return
60.0
*max(pair_.dispersed(), residualAlpha_)
*max(pair_.dispersed(), residualAlpha)
*pair_.continuous().kappa()
/sqr(pair_.dispersed().d());
}
......
......@@ -79,7 +79,7 @@ public:
// Member Functions
//- The heat transfer function K used in the enthalpy equation
tmp<volScalarField> K() const;
tmp<volScalarField> K(const scalar residualAlpha) const;
};
......
......@@ -223,6 +223,66 @@ Foam::BlendedInterfacialModel<ModelType>::K() const
}
template<class ModelType>
Foam::tmp<Foam::volScalarField>
Foam::BlendedInterfacialModel<ModelType>::K(const scalar residualAlpha) const
{
tmp<volScalarField> f1, f2;
if (model_.valid() || model1In2_.valid())
{
f1 = blending_.f1(phase1_, phase2_);
}
if (model_.valid() || model2In1_.valid())
{
f2 = blending_.f2(phase1_, phase2_);
}
tmp<volScalarField> x
(
new volScalarField
(
IOobject
(
ModelType::typeName + ":K",
phase1_.mesh().time().timeName(),
phase1_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
phase1_.mesh(),
dimensionedScalar("zero", ModelType::dimK, 0)
)
);
if (model_.valid())
{
x() += model_->K(residualAlpha)*(scalar(1) - f1() - f2());
}
if (model1In2_.valid())
{
x() += model1In2_->K(residualAlpha)*f1;
}
if (model2In1_.valid())
{
x() += model2In1_->K(residualAlpha)*f2;
}
if
(
correctFixedFluxBCs_
&& (model_.valid() || model1In2_.valid() || model2In1_.valid())
)
{
correctFixedFluxBCs(x());
}
return x;
}
template<class ModelType>
Foam::tmp<Foam::surfaceScalarField>
Foam::BlendedInterfacialModel<ModelType>::Kf() const
......
......@@ -136,6 +136,9 @@ public:
//- Return the blended force coefficient
tmp<volScalarField> K() const;
//- Return the blended force coefficient with a specified residual alpha
tmp<volScalarField> K(const scalar residualAlpha) const;
//- Return the face blended force coefficient
tmp<surfaceScalarField> Kf() const;
......
......@@ -29,4 +29,6 @@ twoPhaseSystem/twoPhaseSystem.C
twoPhaseSystem/newTwoPhaseSystem.C
twoPhaseSystem/twoPhaseSystems.C
reactionThermo/hRefConstThermos.C
LIB = $(FOAM_LIBBIN)/libreactingTwoPhaseSystem
......@@ -192,11 +192,11 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
const volVectorField& U2(pair.phase2().U());
const volScalarField dmdt(this->dmdt(pair));
const volScalarField dmdt12(posPart(dmdt));
const volScalarField dmdt21(negPart(dmdt));
const volScalarField dmdt21(posPart(dmdt));
const volScalarField dmdt12(negPart(dmdt));
*eqns[pair.phase1().name()] += fvm::Sp(dmdt21, U1) - dmdt21*U2;
*eqns[pair.phase2().name()] += dmdt12*U1 - fvm::Sp(dmdt12, U2);
*eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1);
*eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2);
}
return eqnsPtr;
......@@ -283,7 +283,7 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
}
}
// Source term due to mass trasfer
// Source term due to mass transfer
forAllConstIter
(
phaseSystem::phasePairTable,
......@@ -308,17 +308,19 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
const volScalarField& K2(phase2.K());
const volScalarField dmdt(this->dmdt(pair));
const volScalarField dmdt12(posPart(dmdt));
const volScalarField dmdt21(negPart(dmdt));
const volScalarField dmdt21(posPart(dmdt));
const volScalarField dmdt12(negPart(dmdt));
const volScalarField& Tf(*Tf_[pair]);
*eqns[phase1.name()] +=
fvm::Sp(dmdt21, he1) + dmdt21*K1
- dmdt21*(phase2.thermo().he(phase2.thermo().p(), Tf) + K2);
*eqns[phase2.name()] +=
dmdt12*(phase1.thermo().he(phase1.thermo().p(), Tf) + K1)
- fvm::Sp(dmdt12, he2) - dmdt12*K2;
dmdt21*(phase1.thermo().he(phase1.thermo().p(), Tf))
- fvm::Sp(dmdt21, he1)
+ dmdt21*(K2 - K1);
*eqns[phase2.name()] -=
dmdt12*(phase2.thermo().he(phase2.thermo().p(), Tf))
- fvm::Sp(dmdt12, he2)
+ dmdt12*(K1 - K2);
}
return eqnsPtr;
......
......@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "ThermalPhaseChangePhaseSystem.H"
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -36,7 +37,8 @@ ThermalPhaseChangePhaseSystem
:
HeatAndMassTransferPhaseSystem<BasePhaseSystem>(mesh),
volatile_(this->lookup("volatile")),
saturationModel_(saturationModel::New(this->subDict("saturationModel")))
saturationModel_(saturationModel::New(this->subDict("saturationModel"))),
massTransfer_(this->lookup("massTransfer"))
{}
......@@ -143,37 +145,79 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
volScalarField& Tf = *this->Tf_[pair];
Tf = saturationModel_->Tsat(phase1.thermo().p());
Info<< "Tf." << pair.name()
<< ": min = " << min(Tf.internalField())
<< ", mean = " << average(Tf.internalField())
<< ", max = " << max(Tf.internalField())
<< endl;
volScalarField& dmdt(*this->dmdt_[pair]);
volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
const volScalarField& T1(phase1.thermo().T());
const volScalarField& T2(phase2.thermo().T());
const volScalarField& he1(phase1.thermo().he());
const volScalarField& he2(phase2.thermo().he());
volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
volScalarField& dmdt(*this->dmdt_[pair]);
volScalarField& Tf = *this->Tf_[pair];
volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
dmdt =
(H2*(Tf - T1) + H1*(Tf - T2))
/min
if (massTransfer_ )
{
volScalarField H1
(
(pos(dmdt)*he2 + neg(dmdt)*hef2)
- (neg(dmdt)*he1 + pos(dmdt)*hef1),
0.3*(hef1 - hef2)
this->heatTransferModels_[pair][pair.first()]->K(0)
);
volScalarField H2
(
this->heatTransferModels_[pair][pair.second()]->K(0)
);
Tf = saturationModel_->Tsat(phase1.thermo().p());
dmdt =
(H1*(Tf - T1) + H2*(Tf - T2))
/min
(
(pos(dmdt)*he2 + neg(dmdt)*hef2)
- (neg(dmdt)*he1 + pos(dmdt)*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()
<< endl;
}
else
{
dmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
}
volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
// Limit the H[12] boundary field to avoid /0
const scalar HLimit = 1e-4;
H1.boundaryField() =
max(H1.boundaryField(), phase1.boundaryField()*HLimit);
H2.boundaryField() =
max(H2.boundaryField(), phase2.boundaryField()*HLimit);
volScalarField mDotL
(
dmdt*
(
(pos(dmdt)*he2 + neg(dmdt)*hef2)
- (neg(dmdt)*he1 + pos(dmdt)*hef1)
)
);
Tf = (H1*T1 + H2*T2 + mDotL)/(H1 + H2);
Info<< "Tf." << pair.name()
<< ": min = " << min(Tf.internalField())
<< ", mean = " << average(Tf.internalField())
<< ", max = " << max(Tf.internalField())
<< endl;
}
}
......
......@@ -43,6 +43,7 @@ SourceFiles
#include "HeatAndMassTransferPhaseSystem.H"
#include "saturationModel.H"
#include "Switch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -69,6 +70,9 @@ protected:
//- The saturation model used to evaluate Tsat = Tf
autoPtr<saturationModel> saturationModel_;
// Mass transfer enabled
Switch massTransfer_;
public:
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 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 "makeReactionThermo.H"
#include "makeThermo.H"
#include "rhoReactionThermo.H"
#include "heRhoThermo.H"
#include "specie.H"
#include "perfectGas.H"
#include "perfectFluid.H"
#include "rhoConst.H"
#include "incompressiblePerfectGas.H"
#include "hConstThermo.H"
#include "janafThermo.H"
#include "sensibleEnthalpy.H"
#include "absoluteEnthalpy.H"
#include "absoluteInternalEnergy.H"
#include "thermo.H"
#include "constTransport.H"
#include "sutherlandTransport.H"
#include "pureMixture.H"
#include "homogeneousMixture.H"
#include "inhomogeneousMixture.H"
#include "veryInhomogeneousMixture.H"
#include "multiComponentMixture.H"
#include "reactingMixture.H"
#include "singleStepReactingMixture.H"
#include "thermoPhysicsTypes.H"
#include "hRefConstThermo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Thermo type typedefs:
typedef
constTransport
<
species::thermo
<
hRefConstThermo
<
perfectGas<specie>
>,
sensibleEnthalpy
>
> constRefGasHThermoPhysics;
typedef
constTransport
<
species::thermo
<
hRefConstThermo
<
perfectFluid<specie>
>,
sensibleEnthalpy
>
> constRefFluidHThermoPhysics;
typedef
constTransport
<
species::thermo
<
hRefConstThermo
<
perfectGas<specie>
>,
sensibleInternalEnergy
>
> constRefGasEThermoPhysics;
typedef
constTransport
<
species::thermo
<
hRefConstThermo
<
perfectFluid<specie>
>,
sensibleInternalEnergy
>
> constRefFluidEThermoPhysics;
// pureMixture, sensibleEnthalpy:
makeThermo
(
rhoThermo,
heRhoThermo,
pureMixture,
constTransport,
sensibleEnthalpy,
hRefConstThermo,
perfectGas,
specie
);
makeThermo
(
rhoThermo,
heRhoThermo,
pureMixture,
constTransport,
sensibleEnthalpy,
hRefConstThermo,
perfectFluid,
specie
);
makeThermo
(
rhoThermo,
heRhoThermo,
pureMixture,
constTransport,
sensibleEnthalpy,
hRefConstThermo,
rhoConst,
specie
);
// pureMixture, sensibleInternalEnergy:
makeThermo
(
rhoThermo,
heRhoThermo,
pureMixture,
constTransport,
sensibleInternalEnergy,
hRefConstThermo,
perfectGas,
specie
);