multiphaseEulerFoam::populationBalanceModel: improved dilatation treatment

The population balance model considers dilatation originating from density
change and mass transfer via source terms describing nucleation as well as
"drift" of the size distribution to smaller or larger sizes. Numerically, the
treatment does not necessarily equal the total dilatation, hence a correction is
introduced to ensure boundedness of the size group fractions.

Patch contributed by Institute of Fluid Dynamics,
Helmholtz-Zentrum Dresden - Rossendorf (HZDR)
and VTT Technical Research Centre of Finland Ltd.
This commit is contained in:
Henry Weller
2022-04-29 16:18:03 +01:00
parent 58444464aa
commit 376b51b58b
19 changed files with 183 additions and 206 deletions

View File

@ -46,20 +46,23 @@ PopulationBalancePhaseSystem
const diameterModels::populationBalanceModel& popBal =
populationBalances_[popBali];
forAll(popBal.velocityGroups(), velGrp1i)
{
const diameterModels::velocityGroup& velGrp1 =
popBal.velocityGroups()[velGrp1i];
for
forAllConstIter
(
label velGrp2i = velGrp1i + 1;
velGrp2i < popBal.velocityGroups().size();
++ velGrp2i
HashTable<const diameterModels::velocityGroup*>,
popBal.velocityGroupPtrs(),
iter1
)
{
const diameterModels::velocityGroup& velGrp2 =
popBal.velocityGroups()[velGrp2i];
const diameterModels::velocityGroup& velGrp1 = *iter1();
forAllConstIter
(
HashTable<const diameterModels::velocityGroup*>,
popBal.velocityGroupPtrs(),
iter2
)
{
const diameterModels::velocityGroup& velGrp2 = *iter2();
const phaseInterface interface
(
@ -67,6 +70,13 @@ PopulationBalancePhaseSystem
velGrp2.phase()
);
if
(
&velGrp1 != &velGrp2
&&
!dmdtfs_.found(interface)
)
{
this->template validateMassTransfer
<diameterModels::populationBalanceModel>(interface);
@ -92,6 +102,7 @@ PopulationBalancePhaseSystem
}
}
}
}
}

View File

@ -26,9 +26,9 @@ License
#include "fractal.H"
#include "addToRunTimeSelectionTable.H"
#include "sinteringModel.H"
#include "fvm.H"
#include "fvcDdt.H"
#include "fvcDiv.H"
#include "fvmDdt.H"
#include "fvmDiv.H"
#include "fvmSup.H"
#include "mixedFvPatchField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -122,7 +122,7 @@ Foam::diameterModels::shapeModels::fractal::fractal
if
(
isA<mixedFvPatchScalarField>(kappa_.boundaryFieldRef()[patchi])
isA<const mixedFvPatchScalarField>(kappa_.boundaryField()[patchi])
)
{
mixedFvPatchScalarField& kappa =
@ -205,12 +205,12 @@ void Foam::diameterModels::shapeModels::fractal::correct()
fvScalarMatrix kappaEqn
(
fvm::ddt(alpha, fi, kappa_) + fvm::div(fAlphaPhi, kappa_)
fvm::ddt(alpha, fi, kappa_)
+ fvm::div(fAlphaPhi, kappa_)
==
- sinteringModel_->R()
+ Su_
- fvm::Sp(popBal.Sp(fi.i()())*fi, kappa_)
- fvm::Sp(popBal.Sp(fi.i())*fi, kappa_)
- correction
(
fvm::Sp

View File

@ -75,7 +75,7 @@ Foam::diameterModels::sizeGroup::sizeGroup
if
(
isA<mixedFvPatchScalarField>(this->boundaryFieldRef()[patchi])
isA<const mixedFvPatchScalarField>(this->boundaryField()[patchi])
)
{
mixedFvPatchScalarField& f =
@ -108,7 +108,7 @@ Foam::diameterModels::sizeGroup::clone() const
}
const Foam::autoPtr<Foam::label>& Foam::diameterModels::sizeGroup::i() const
const Foam::label& Foam::diameterModels::sizeGroup::i() const
{
if (!i_.valid())
{
@ -127,7 +127,7 @@ const Foam::autoPtr<Foam::label>& Foam::diameterModels::sizeGroup::i() const
}
}
return i_;
return i_();
}

View File

@ -192,8 +192,8 @@ public:
//- Return const-reference to diameterModel of the phase
inline const autoPtr<shapeModel>& shapeModelPtr() const;
//- Return label of the sizeGroup within the population balance
const autoPtr<label>& i() const;
//- Return index of the size group within the population balance
const label& i() const;
//- Return representative surface area of the sizeGroup
const tmp<volScalarField> a() const;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2017-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2017-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -63,7 +63,7 @@ void Foam::diameterModels::driftModels::constantDrift::addToDriftRate
const label i
)
{
driftRate += rate_;
driftRate += popBal_.sizeGroups()[i].phase()*rate_;
}

View File

@ -70,11 +70,11 @@ void Foam::diameterModels::driftModels::densityChangeDrift::addToDriftRate
const volScalarField& rho = phase.thermo().rho();
driftRate -=
fi.x()/(rho*max(alpha, phase.residualAlpha()))
fi.x()
*(
fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
);
)/rho;
}

View File

@ -97,15 +97,20 @@ void Foam::diameterModels::driftModels::phaseChange::precompute()
forAll(interfaces_, k)
{
forAll(popBal_.velocityGroups(), j)
forAllConstIter
(
HashTable<const velocityGroup*>,
popBal_.velocityGroupPtrs(),
iter
)
{
const velocityGroup& vgj = popBal_.velocityGroups()[j];
const velocityGroup& velGrp = *iter();
if (interfaces_[k].contains(vgj.phase()))
if (interfaces_[k].contains(velGrp.phase()))
{
forAll(vgj.sizeGroups(), i)
forAll(velGrp.sizeGroups(), i)
{
const sizeGroup& fi = vgj.sizeGroups()[i];
const sizeGroup& fi = velGrp.sizeGroups()[i];
if (numberWeighted_)
{
@ -128,11 +133,11 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
const label i
)
{
const velocityGroup& vg = popBal_.sizeGroups()[i].VelocityGroup();
const velocityGroup& velGrp = popBal_.sizeGroups()[i].VelocityGroup();
forAll(interfaces_, k)
{
if (interfaces_[k].contains(vg.phase()))
if (interfaces_[k].contains(velGrp.phase()))
{
const volScalarField& dmidtf =
popBal_.mesh().lookupObject<volScalarField>
@ -149,7 +154,7 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
);
const scalar dmidtfSign =
interfaces_[k].index(vg.phase()) == 0 ? +1 : -1;
interfaces_[k].index(velGrp.phase()) == 0 ? +1 : -1;
const sizeGroup& fi = popBal_.sizeGroups()[i];
@ -163,7 +168,7 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
dDriftRate.ref() *= fi.a();
}
driftRate += dDriftRate;
driftRate += velGrp.phase()*dDriftRate;
}
}
}

View File

@ -31,9 +31,12 @@ License
#include "driftModel.H"
#include "nucleationModel.H"
#include "surfaceTensionModel.H"
#include "fvm.H"
#include "fvmDdt.H"
#include "fvmDiv.H"
#include "fvmSup.H"
#include "fvcDdt.H"
#include "fvcDiv.H"
#include "fvcSup.H"
#include "shapeModel.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -60,12 +63,26 @@ void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
if (velGroup.popBalName() == this->name())
{
velocityGroups_.resize(velocityGroups_.size() + 1);
velocityGroupPtrs_.insert(velGroup.phase().name(), &velGroup);
velocityGroups_.set
dilatationErrors_.insert
(
velocityGroups_.size() - 1,
&const_cast<velocityGroup&>(velGroup)
velGroup.phase().name(),
volScalarField
(
IOobject
(
IOobject::groupName
(
"dilatationError",
velGroup.phase().name()
),
fluid_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar(inv(dimTime), 0)
)
);
forAll(velGroup.sizeGroups(), i)
@ -446,22 +463,21 @@ void Foam::diameterModels::populationBalanceModel::drift
const sizeGroup& fe = sizeGroups()[i+1];
volScalarField& Sue = Sui_;
Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_()*fp.phase();
Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_();
Sue =
fe.x()/(fp.x()*(fe.x() - fp.x()))
*pos(driftRate_())*driftRate_()*fp*fp.phase();
fe.x()/(fp.x()*(fe.x() - fp.x()))*pos(driftRate_())*driftRate_()*fp;
Su_[i+1] += Sue;
const phaseInterface interfaceij(fp.phase(), fe.phase());
const phaseInterface interfacepe(fp.phase(), fe.phase());
if (pDmdt_.found(interfaceij))
if (pDmdt_.found(interfacepe))
{
const scalar dmdtSign =
interfaceij.index(fp.phase()) == 0 ? +1 : -1;
interfacepe.index(fp.phase()) == 0 ? +1 : -1;
*pDmdt_[interfaceij] -= dmdtSign*Sue*fp.phase().rho();
*pDmdt_[interfacepe] -= dmdtSign*Sue*fp.phase().rho();
}
sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model);
@ -469,7 +485,7 @@ void Foam::diameterModels::populationBalanceModel::drift
if (i == sizeGroups().size() - 1)
{
Sp_[i] -= pos(driftRate_())*driftRate_()*fp.phase()/fp.x();
Sp_[i] -= pos(driftRate_())*driftRate_()/fp.x();
}
if (i > 0)
@ -477,22 +493,21 @@ void Foam::diameterModels::populationBalanceModel::drift
const sizeGroup& fw = sizeGroups()[i-1];
volScalarField& Suw = Sui_;
Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_()*fp.phase();
Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_();
Suw =
fw.x()/(fp.x()*(fw.x() - fp.x()))
*neg(driftRate_())*driftRate_()*fp*fp.phase();
fw.x()/(fp.x()*(fw.x() - fp.x()))*neg(driftRate_())*driftRate_()*fp;
Su_[i-1] += Suw;
const phaseInterface interfaceih(fp.phase(), fw.phase());
const phaseInterface interfacepw(fp.phase(), fw.phase());
if (pDmdt_.found(interfaceih))
if (pDmdt_.found(interfacepw))
{
const scalar dmdtSign =
interfaceih.index(fp.phase()) == 0 ? +1 : -1;
interfacepw.index(fp.phase()) == 0 ? +1 : -1;
*pDmdt_[interfaceih] -= dmdtSign*Suw*fp.phase().rho();
*pDmdt_[interfacepw] -= dmdtSign*Suw*fp.phase().rho();
}
sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model);
@ -500,7 +515,7 @@ void Foam::diameterModels::populationBalanceModel::drift
if (i == 0)
{
Sp_[i] -= neg(driftRate_())*driftRate_()*fp.phase()/fp.x();
Sp_[i] -= neg(driftRate_())*driftRate_()/fp.x();
}
}
@ -609,13 +624,43 @@ void Foam::diameterModels::populationBalanceModel::sources()
}
void Foam::diameterModels::populationBalanceModel::correctDilatationError()
{
forAllIter
(
HashTable<volScalarField>,
dilatationErrors_,
iter
)
{
volScalarField& dilatationError = iter();
const word& phaseName = iter.key();
const phaseModel& phase = fluid_.phases()[phaseName];
const velocityGroup& velGroup = *velocityGroupPtrs_[phaseName];
const volScalarField& alpha = phase;
const volScalarField& rho = phase.thermo().rho();
dilatationError =
fvc::ddt(alpha) + fvc::div(phase.alphaPhi())
- (fluid_.fvModels().source(alpha, rho) & rho)/rho;
forAll(velGroup.sizeGroups(), i)
{
const sizeGroup& fi = velGroup.sizeGroups()[i];
dilatationError -= Su_[fi.i()] - fvc::Sp(Sp_[fi.i()], fi);
}
}
}
void Foam::diameterModels::populationBalanceModel::calcAlphas()
{
alphas_() = Zero;
forAll(velocityGroups_, v)
forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
{
const phaseModel& phase = velocityGroups_[v].phase();
const phaseModel& phase = iter()->phase();
alphas_() += max(phase, phase.residualAlpha());
}
@ -637,9 +682,9 @@ Foam::diameterModels::populationBalanceModel::calcDsm()
volScalarField& invDsm = tInvDsm.ref();
forAll(velocityGroups_, v)
forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
{
const phaseModel& phase = velocityGroups_[v].phase();
const phaseModel& phase = iter()->phase();
invDsm += max(phase, phase.residualAlpha())/(phase.d()*alphas_());
}
@ -652,9 +697,9 @@ void Foam::diameterModels::populationBalanceModel::calcVelocity()
{
U_() = Zero;
forAll(velocityGroups_, v)
forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
{
const phaseModel& phase = velocityGroups_[v].phase();
const phaseModel& phase = iter()->phase();
U_() += max(phase, phase.residualAlpha())*phase.U()/alphas_();
}
@ -704,7 +749,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
IOobject::groupName("alpha", dict_.lookup("continuousPhase"))
)
),
velocityGroups_(),
sizeGroups_(),
v_(),
delta_(),
@ -890,7 +934,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
);
}
if (velocityGroups_.size() > 1)
if (velocityGroupPtrs_.size() > 1)
{
alphas_.set
(
@ -1060,6 +1104,8 @@ void Foam::diameterModels::populationBalanceModel::solve()
sources();
}
correctDilatationError();
maxInitialResidual = 0;
forAll(sizeGroups(), i)
@ -1068,15 +1114,18 @@ void Foam::diameterModels::populationBalanceModel::solve()
const phaseModel& phase = fi.phase();
const volScalarField& alpha = phase;
const volScalarField& rho = phase.thermo().rho();
const volScalarField& dilatationError =
dilatationErrors_[phase.name()];
fvScalarMatrix sizeGroupEqn
(
fvm::ddt(alpha, fi) + fvm::div(phase.alphaPhi(), fi)
fvm::ddt(alpha, fi)
+ fvm::div(phase.alphaPhi(), fi)
- fvm::Sp(dilatationError, fi)
==
Su_[i]
fvc::Su(Su_[i], fi)
- fvm::Sp(Sp_[i], fi)
+ fluid_.fvModels().source(alpha, rho, fi)/rho
- correction
(
fvm::Sp
@ -1121,7 +1170,7 @@ void Foam::diameterModels::populationBalanceModel::solve()
void Foam::diameterModels::populationBalanceModel::correct()
{
if (velocityGroups_.size() > 1)
if (velocityGroupPtrs_.size() > 1)
{
calcAlphas();
dsm_() = calcDsm();

View File

@ -257,10 +257,10 @@ private:
//- Continuous phase
const phaseModel& continuousPhase_;
//- velocityGroups belonging to this populationBalance
UPtrList<velocityGroup> velocityGroups_;
//- Velocity groups belonging to this populationBalance
HashTable<const velocityGroup*> velocityGroupPtrs_;
//- sizeGroups belonging to this populationBalance
//- Size groups belonging to this populationBalance
UPtrList<sizeGroup> sizeGroups_;
//- sizeGroup boundaries
@ -275,6 +275,9 @@ private:
//- Implicitly treated sources
PtrList<volScalarField> Sp_;
//- Dilatation errors
HashTable<volScalarField> dilatationErrors_;
//- Field for caching sources
volScalarField Sui_;
@ -357,7 +360,7 @@ private:
void sources();
void dmdt();
void correctDilatationError();
void calcAlphas();
@ -454,15 +457,18 @@ public:
//- Return continuous phase
inline const phaseModel& continuousPhase() const;
//- Return the velocityGroups belonging to this populationBalance
inline const UPtrList<velocityGroup>& velocityGroups() const;
//- Return the velocity groups belonging to this populationBalance
inline const HashTable<const velocityGroup*>& velocityGroupPtrs() const;
//- Return the sizeGroups belonging to this populationBalance
//- Return the size groups belonging to this populationBalance
inline const UPtrList<sizeGroup>& sizeGroups() const;
//- Return implicit source terms
inline const volScalarField& Sp(const label i) const;
//- Return dilatation obtained from source terms
inline const HashTable<volScalarField>& sourceDilatation() const;
//- Return coalescence relevant size group pairs
inline const List<labelPair>& coalescencePairs() const;

View File

@ -82,10 +82,10 @@ Foam::diameterModels::populationBalanceModel::continuousPhase() const
}
inline const Foam::UPtrList<Foam::diameterModels::velocityGroup>&
Foam::diameterModels::populationBalanceModel::velocityGroups() const
inline const Foam::HashTable<const Foam::diameterModels::velocityGroup*>&
Foam::diameterModels::populationBalanceModel::velocityGroupPtrs() const
{
return velocityGroups_;
return velocityGroupPtrs_;
}
@ -120,13 +120,13 @@ Foam::diameterModels::populationBalanceModel::binaryBreakupPairs() const
inline const Foam::volScalarField&
Foam::diameterModels::populationBalanceModel::alphas() const
{
if (velocityGroups_.size() > 1)
if (velocityGroupPtrs_.size() > 1)
{
return alphas_();
}
else
{
return velocityGroups_.first().phase();
return velocityGroupPtrs_[velocityGroupPtrs_.begin().key()]->phase();
}
}
@ -134,13 +134,14 @@ Foam::diameterModels::populationBalanceModel::alphas() const
inline const Foam::volVectorField&
Foam::diameterModels::populationBalanceModel::U() const
{
if (velocityGroups_.size() > 1)
if (velocityGroupPtrs_.size() > 1)
{
return U_();
}
else
{
return velocityGroups_.first().phase().U();
return
velocityGroupPtrs_[velocityGroupPtrs_.begin().key()]->phase().U();
}
}

View File

@ -14,51 +14,6 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
airSource
{
type massSource;
selectionMode all;
massFlowRate 0.02186682;
phase air;
rho thermo:rho.air;
fieldValues
{
f0.air 0;
f1.air 0;
f2.air 0;
f3.air 0;
f4.air 0;
f5.air 0;
f6.air 0;
f7.air 0;
f8.air 0;
f9.air 0;
f10.air 0;
f11.air 0;
f12.air 0;
f13.air 0;
f14.air 0;
f15.air 0;
f16.air 0;
f17.air 0;
f18.air 0;
f19.air 0;
f20.air 0;
f21.air 0;
f22.air 0;
f23.air 0;
f24.air 0;
f25.air 0;
f26.air 0;
f27.air 0;
f28.air 0;
}
}
waterSink
{
type massSource;
@ -71,7 +26,9 @@ waterSink
rho thermo:rho.water;
fieldValues
{}
{
U.water (0 0 0);
}
}
// ************************************************************************* //

View File

@ -20,7 +20,7 @@ thermoType
mixture pureMixture;
transport const;
thermo hConst;
equationOfState rhoConst;
equationOfState perfectGas;
specie specie;
energy sensibleInternalEnergy;
}
@ -31,10 +31,6 @@ mixture
{
molWeight 28.9;
}
equationOfState
{
rho 1;
}
thermodynamics
{
Cp 1007;

View File

@ -60,7 +60,7 @@ PIMPLE
pRefCell 0;
pRefValue 1e5;
flow no;
flow yes;
thermophysics yes;
models yes;
}

View File

@ -12,7 +12,6 @@ gnuplot<<EOF
set xlabel 'd (mm)'
set ylabel 'y (m)'
set yrange [0.05:]
set key bottom center
latestTime = system("foamListTimes -case .. -latestTime")

View File

@ -14,52 +14,7 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
airSource
{
type massSource;
selectionMode all;
massFlowRate -0.030785654;
phase air;
rho thermo:rho.air;
fieldValues
{
f0.air 0;
f1.air 0;
f2.air 0;
f3.air 0;
f4.air 0;
f5.air 0;
f6.air 0;
f7.air 0;
f8.air 0;
f9.air 0;
f10.air 0;
f11.air 0;
f12.air 0;
f13.air 0;
f14.air 0;
f15.air 0;
f16.air 0;
f17.air 0;
f18.air 0;
f19.air 0;
f20.air 0;
f21.air 0;
f22.air 0;
f23.air 0;
f24.air 0;
f25.air 0;
f26.air 0;
f27.air 0;
f28.air 0;
}
}
waterSink
waterSource
{
type massSource;
@ -71,7 +26,9 @@ waterSink
rho thermo:rho.water;
fieldValues
{}
{
U.water (0 0 0);
}
}
// ************************************************************************* //

View File

@ -20,7 +20,7 @@ thermoType
mixture pureMixture;
transport const;
thermo hConst;
equationOfState rhoConst;
equationOfState perfectGas;
specie specie;
energy sensibleInternalEnergy;
}
@ -31,10 +31,6 @@ mixture
{
molWeight 28.9;
}
equationOfState
{
rho 1;
}
thermodynamics
{
Cp 1007;

View File

@ -60,7 +60,7 @@ PIMPLE
pRefCell 0;
pRefValue 1e5;
flow no;
flow yes;
thermophysics yes;
models yes;
}

View File

@ -28,7 +28,7 @@ solvers
tolerance 1e-4;
scale true;
solveOnFinalIterOnly true;
sourceUpdateInterval 1;
sourceUpdateInterval 20;
}
p_rgh

View File

@ -28,7 +28,7 @@ solvers
tolerance 1e-4;
scale true;
solveOnFinalIterOnly true;
sourceUpdateInterval 1;
sourceUpdateInterval 20;
}
p_rgh