zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of github.com-OpenFOAM:OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
Henry Weller
2021-06-08 13:25:38 +01:00
parent fd9a97c42e bc844c383c
commit e6ce04f070
19 changed files with 297 additions and 396 deletions
Download Patch File Download Diff File Expand all files Collapse all files

34
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
View File

@ -32,16 +32,9 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::dmdtfTable>
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
totalDmdtfs() const
void Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
correctDmdtfs()
{
autoPtr<phaseSystem::dmdtfTable> totalDmdtfsPtr
(
new phaseSystem::dmdtfTable
);
phaseSystem::dmdtfTable& totalDmdtfs = totalDmdtfsPtr();
forAllConstIter
(
interfaceCompositionModelTable,
@ -52,6 +45,8 @@ totalDmdtfs() const
const phasePair& pair =
this->phasePairs_[interfaceCompositionModelIter.key()];
*dmdtfs_[pair] = Zero;
forAllConstIter(phasePair, pair, pairIter)
{
const autoPtr<interfaceCompositionModel>& compositionModelPtr =
@ -73,28 +68,15 @@ totalDmdtfs() const
{
const word& specie = *specieIter;
tmp<volScalarField> dmidtf
(
*dmdtfs_[pair] +=
(pairIter.index() == 0 ? +1 : -1)
*(
*(*dmidtfSus_[pair])[specie]
+ *(*dmidtfSps_[pair])[specie]*phase.Y(specie)
)
);
if (totalDmdtfs.found(pair))
{
*totalDmdtfs[pair] += dmidtf;
}
else
{
totalDmdtfs.insert(pair, dmidtf.ptr());
}
);
}
}
}
return totalDmdtfsPtr;
}
@ -610,7 +592,7 @@ correct()
}
}
dmdtfs_ = totalDmdtfs();
correctDmdtfs();
}
@ -620,7 +602,7 @@ correctSpecies()
{
BasePhaseSystem::correctSpecies();
dmdtfs_ = totalDmdtfs();
correctDmdtfs();
}

4
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H
View File

@ -108,8 +108,8 @@ class InterfaceCompositionPhaseChangePhaseSystem
// Private Member Functions
//- Return mass transfers across each interface
autoPtr<phaseSystem::dmdtfTable> totalDmdtfs() const;
//- Correct mass transfers across each interface
void correctDmdtfs();
//- Return species mass transfers across each interface
autoPtr<phaseSystem::dmidtfTable> totalDmidtfs() const;

24
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -184,8 +184,8 @@ Foam::diameterModels::shapeModels::fractal::dColl() const
volScalarField& dColl = tDColl.ref();
dColl =
6.0/kappa_
*pow(sizeGroup_.x()*pow3(kappa_)/(36.0*pi*alphaC_), 1.0/Df_);
6/kappa_
*pow(sizeGroup_.x()*pow3(kappa_)/(36*pi*alphaC_), 1/Df_);
return tDColl;
}
@ -216,7 +216,7 @@ void Foam::diameterModels::shapeModels::fractal::correct()
==
- sinteringModel_->R()
+ Su_
- fvm::SuSp(popBal.SuSp(fi.i()())*fi, kappa_)
- fvm::Sp(popBal.Sp(fi.i()())*fi, kappa_)
+ fvc::ddt(fi.phase().residualAlpha(), kappa_)
- fvm::ddt(fi.phase().residualAlpha(), kappa_)
);
@ -230,8 +230,8 @@ void Foam::diameterModels::shapeModels::fractal::correct()
kappa_ =
min
(
max(kappa_, 6.0/sizeGroup_.dSph()),
6.0/popBal.sizeGroups().first().dSph()
max(kappa_, 6/sizeGroup_.dSph()),
6/popBal.sizeGroups().first().dSph()
);
kappa_.correctBoundaryConditions();
@ -283,7 +283,7 @@ void Foam::diameterModels::shapeModels::fractal::addDrift
fu.shapeModelPtr()()
);
volScalarField dp(6.0/sourceKappa);
volScalarField dp(6/sourceKappa);
const volScalarField a(sourceKappa*fu.x());
const dimensionedScalar dv(sizeGroup_.x() - fu.x());
@ -292,18 +292,18 @@ void Foam::diameterModels::shapeModels::fractal::addDrift
(2.0/3.0)*dv
*(
sourceKappa
+ sourceShape.Df_*(1.0/sourceShape.d() - 1.0/dp)
+ sourceShape.Df_*(1/sourceShape.d() - 1/dp)
)
);
dp += 6.0*(dv*a - fu.x()*da1)/sqr(a);
dp += 6*(dv*a - fu.x()*da1)/sqr(a);
const volScalarField np(6.0*sizeGroup_.x()/pi/pow3(dp));
const volScalarField dc(dp*pow(np/alphaC_, 1.0/Df_));
const volScalarField np(6*sizeGroup_.x()/pi/pow3(dp));
const volScalarField dc(dp*pow(np/alphaC_, 1/Df_));
const volScalarField da2
(
dv*(4.0/dp + 2.0*Df_/3.0*(1.0/dc - 1.0/dp))
dv*(4/dp + 2*Df_/3*(1/dc - 1/dp))
);
Su_ += (a + 0.5*da1 + 0.5*da2)/sizeGroup_.x()*Su;

6
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/sinteringModels/KochFriedlander/KochFriedlander.C
View File

@ -82,7 +82,7 @@ Foam::diameterModels::shapeModels::sinteringModels::KochFriedlander::tau() const
(
"tau",
fractal_.SizeGroup().mesh(),
dimensionedScalar(dimTime, 0.0)
dimensionedScalar(dimTime, Zero)
)
);
@ -120,7 +120,7 @@ Foam::diameterModels::shapeModels::sinteringModels::KochFriedlander::R() const
fi.mesh()
),
fi.mesh(),
dimensionedScalar(inv(dimTime), 0)
dimensionedScalar(inv(dimTime), Zero)
);
volScalarField::Internal tau(this->tau());
@ -130,7 +130,7 @@ Foam::diameterModels::shapeModels::sinteringModels::KochFriedlander::R() const
R[celli] = fi[celli]*alpha[celli]/tau[celli];
}
return fvm::Sp(R, kappai) - 6.0/fi.dSph()*R;
return fvm::Sp(R, kappai) - 6/fi.dSph()*R;
}

4
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/spherical/spherical.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -79,7 +79,7 @@ Foam::diameterModels::shapeModels::spherical::a() const
(
"a",
sizeGroup_.mesh(),
6.0/sizeGroup_.dSph()*sizeGroup_.x()
6/sizeGroup_.dSph()*sizeGroup_.x()
)
);
}

2
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C
View File

@ -67,7 +67,7 @@ Foam::diameterModels::sizeGroup::sizeGroup
phase_(phase),
velocityGroup_(velocityGroup),
dSph_("dSph", dimLength, dict),
x_("x", pi/6.0*pow3(dSph_)),
x_("x", pi/6*pow3(dSph_)),
value_(dict.lookup<scalar>("value"))
{
// Adjust refValue at mixedFvPatchField boundaries

2
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C
View File

@ -64,7 +64,7 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::velocityGroup::dsm() const
invDsm += fi.a()*fi/fi.x();
}
return 6.0/tInvDsm;
return 6/tInvDsm;
}

4
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LehrMilliesMewes/LehrMilliesMewes.C
View File

@ -102,9 +102,9 @@ addToBinaryBreakupRate
binaryBreakupRate +=
0.5*pow(fj.dSph()/L, 5.0/3.0)
*exp(-sqrt(2.0)/pow3(fj.dSph()/L))
*6.0/pow(pi, 1.5)/pow3(fi.dSph()/L)
*6/pow(pi, 1.5)/pow3(fi.dSph()/L)
*exp(-9.0/4.0*sqr(log(pow(2.0, 0.4)*fi.dSph()/L)))
/max(1.0 + erf(1.5*log(pow(2.0, 1.0/15.0)*fj.dSph()/L)), small)
/max(1 + erf(1.5*log(pow(2.0, 1.0/15.0)*fj.dSph()/L)), small)
/(T*pow3(L));
}

6
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C
View File

@ -180,7 +180,7 @@ Foam::diameterModels::binaryBreakupModels::LuoSvendsen::addToBinaryBreakupRate
const volScalarField b
(
12.0*cf*popBal_.sigmaWithContinuousPhase(fi.phase())
12*cf*popBal_.sigmaWithContinuousPhase(fi.phase())
/(
beta_*continuousPhase.rho()*pow(fj.dSph(), 5.0/3.0)
*pow(popBal_.continuousTurbulence().epsilon(), 2.0/3.0)
@ -191,12 +191,12 @@ Foam::diameterModels::binaryBreakupModels::LuoSvendsen::addToBinaryBreakupRate
const volScalarField tMin(b/pow(xiMin, 11.0/3.0));
volScalarField integral(3.0/(11.0*pow(b, 8.0/11.0)));
volScalarField integral(3/(11*pow(b, 8.0/11.0)));
forAll(integral, celli)
{
integral[celli] *=
2.0*pow(b[celli], 3.0/11.0)*tgamma(5.0/11.0)
2*pow(b[celli], 3.0/11.0)*tgamma(5.0/11.0)
*(
gammaUpperReg5by11_->value(b[celli])
- gammaUpperReg5by11_->value(tMin[celli])

4
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/powerLawUniformBinary/powerLawUniformBinary.C
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-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -74,7 +74,7 @@ addToBinaryBreakupRate
const sizeGroup& fj = popBal_.sizeGroups()[j];
binaryBreakupRate.primitiveFieldRef() +=
pow(fj.x().value(), power_)*2.0/fj.x().value();
pow(fj.x().value(), power_)*2/fj.x().value();
}

4
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/LehrMilliesMewesCoalescence/LehrMilliesMewesCoalescence.C
View File

@ -95,11 +95,11 @@ addToCoalescenceRate
);
coalescenceRate +=
pi/4.0*sqr(fi.dSph() + fj.dSph())*min(uChar, uCrit_)
pi/4*sqr(fi.dSph() + fj.dSph())*min(uChar, uCrit_)
*exp
(
- sqr(cbrt(alphaMax_)
/cbrt(max(popBal_.alphas(), fi.phase().residualAlpha())) - 1.0)
/cbrt(max(popBal_.alphas(), fi.phase().residualAlpha())) - 1)
);
}

10
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/Luo/Luo.C
View File

@ -62,7 +62,7 @@ Luo
:
coalescenceModel(popBal, dict),
beta_(dimensionedScalar::lookupOrDefault("beta", dict, dimless, 2.05)),
C1_(dimensionedScalar::lookupOrDefault("C1", dict, dimless, 1.0))
C1_(dimensionedScalar::lookupOrDefault("C1", dict, dimless, 1))
{}
@ -102,18 +102,18 @@ addToCoalescenceRate
(
sqrt(beta_)
*cbrt(popBal_.continuousTurbulence().epsilon()*fi.dSph())
*sqrt(1.0 + pow(xi, -2.0/3.0))
*sqrt(1 + pow(xi, -2.0/3.0))
);
coalescenceRate +=
pi/4.0*sqr(fi.dSph() + fj.dSph())*uij
pi/4*sqr(fi.dSph() + fj.dSph())*uij
*exp
(
- C1_
*sqrt(0.75*(1.0 + sqr(xi))*(1.0 + pow3(xi)))
*sqrt(0.75*(1 + sqr(xi))*(1 + pow3(xi)))
/(
sqrt(fi.phase().rho()/continuousPhase.rho()
+ vm.Cvm())*pow3(1.0 + xi)
+ vm.Cvm())*pow3(1 + xi)
)
*sqrt
(

2
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/LaakkonenDaughterSizeDistribution/LaakkonenDaughterSizeDistribution.C
View File

@ -115,7 +115,7 @@ LaakkonenDaughterSizeDistribution::calcNik
if (k == 0)
{
return 1.0;
return 1;
}
return

2
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C
View File

@ -85,7 +85,7 @@ Foam::diameterModels::daughterSizeDistributionModels::uniformBinary::calcNik
{
if (k == 0)
{
return 1.0;
return 1;
}
return (sizeGroups[i+1].x() - xi)/(xk - x0);

4
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/reactionDriven/reactionDriven.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -139,7 +139,7 @@ Foam::diameterModels::nucleationModels::reactionDriven::addToNucleationRate
velGroup_.phase().name() == pair_.first() ? +1 : -1;
nucleationRate +=
popBal_.eta(i, pi/6.0*pow3(dNuc_))*dmidtfSign*dmidtf/rho/fi.x();
popBal_.eta(i, pi/6*pow3(dNuc_))*dmidtfSign*dmidtf/rho/fi.x();
}

2
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C
View File

@ -174,7 +174,7 @@ Foam::diameterModels::nucleationModels::wallBoiling::addToNucleationRate
const labelList& faceCells = alphatw.patch().faceCells();
dimensionedScalar unitLength("unitLength", dimLength, 1.0);
dimensionedScalar unitLength("unitLength", dimLength, 1);
forAll(alphatw, facei)
{

462
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
View File

@ -88,9 +88,9 @@ void Foam::diameterModels::populationBalanceModel::registerSizeGroups
{
if
(
sizeGroups_.size() != 0
sizeGroups().size() != 0
&&
group.x().value() <= sizeGroups_.last().x().value()
group.x().value() <= sizeGroups().last().x().value()
)
{
FatalErrorInFunction
@ -99,53 +99,44 @@ void Foam::diameterModels::populationBalanceModel::registerSizeGroups
<< exit(FatalError);
}
sizeGroups_.resize(sizeGroups_.size() + 1);
sizeGroups_.set(sizeGroups_.size() - 1, &group);
sizeGroups_.resize(sizeGroups().size() + 1);
sizeGroups_.set(sizeGroups().size() - 1, &group);
// Grid generation over property space
if (sizeGroups_.size() == 1)
if (sizeGroups().size() == 1)
{
// Set the first sizeGroup boundary to the representative value of
// the first sizeGroup
v_.append
(
new dimensionedScalar
(
"v",
sizeGroups_.last().x()
sizeGroups().last().x()
)
);
// Set the last sizeGroup boundary to the representative size of the
// last sizeGroup
v_.append
(
new dimensionedScalar
(
"v",
sizeGroups_.last().x()
sizeGroups().last().x()
)
);
}
else
{
// Overwrite the next-to-last sizeGroup boundary to lie halfway between
// the last two sizeGroups
v_.last() =
0.5
*(
sizeGroups_[sizeGroups_.size()-2].x()
+ sizeGroups_.last().x()
sizeGroups()[sizeGroups().size()-2].x()
+ sizeGroups().last().x()
);
// Set the last sizeGroup boundary to the representative size of the
// last sizeGroup
v_.append
(
new dimensionedScalar
(
"v",
sizeGroups_.last().x()
sizeGroups().last().x()
)
);
}
@ -167,13 +158,13 @@ void Foam::diameterModels::populationBalanceModel::registerSizeGroups
)
);
SuSp_.append
Sp_.append
(
new volScalarField
(
IOobject
(
"SuSp",
"Sp",
fluid_.time().timeName(),
mesh_
),
@ -226,23 +217,21 @@ void Foam::diameterModels::populationBalanceModel::createPhasePairs()
void Foam::diameterModels::populationBalanceModel::precompute()
{
calcDeltas();
forAll(coalescence_, model)
forAll(coalescenceModels_, model)
{
coalescence_[model].precompute();
coalescenceModels_[model].precompute();
}
forAll(breakup_, model)
forAll(breakupModels_, model)
{
breakup_[model].precompute();
breakupModels_[model].precompute();
breakup_[model].dsdPtr()->precompute();
breakupModels_[model].dsdPtr()->precompute();
}
forAll(binaryBreakup_, model)
forAll(binaryBreakupModels_, model)
{
binaryBreakup_[model].precompute();
binaryBreakupModels_[model].precompute();
}
forAll(drift_, model)
@ -257,42 +246,37 @@ void Foam::diameterModels::populationBalanceModel::precompute()
}
void Foam::diameterModels::populationBalanceModel::
birthByCoalescence
void Foam::diameterModels::populationBalanceModel::birthByCoalescence
(
const label j,
const label k
)
{
const sizeGroup& fj = sizeGroups_[j];
const sizeGroup& fk = sizeGroups_[k];
const sizeGroup& fj = sizeGroups()[j];
const sizeGroup& fk = sizeGroups()[k];
dimensionedScalar Eta;
dimensionedScalar v = fj.x() + fk.x();
for (label i = j; i < sizeGroups_.size(); i++)
for (label i = j; i < sizeGroups().size(); i++)
{
// Calculate fraction for intra-interval events
Eta = eta(i, v);
if (Eta.value() == 0) continue;
const sizeGroup& fi = sizeGroups_[i];
const sizeGroup& fi = sizeGroups()[i];
// Avoid double counting of events
if (j == k)
{
Sui_ =
0.5*fi.x()*coalescenceRate_()*Eta
*fj*fj.phase()/fj.x()
*fk*fk.phase()/fk.x();
0.5*fi.x()/(fj.x()*fk.x())*Eta
*coalescenceRate_()*fj*fj.phase()*fk*fk.phase();
}
else
{
Sui_ =
fi.x()*coalescenceRate_()*Eta
*fj*fj.phase()/fj.x()
*fk*fk.phase()/fk.x();
fi.x()/(fj.x()*fk.x())*Eta
*coalescenceRate_()*fj*fj.phase()*fk*fk.phase();
}
Su_[i] += Sui_;
@ -310,7 +294,7 @@ birthByCoalescence
Pair<word>::compare(pDmdt_.find(pairij).key(), pairij)
);
*pDmdt_[pairij] += dmdtSign*fj.x()/v*Sui_*fi.phase().rho();
*pDmdt_[pairij] += dmdtSign*fj.x()/v*Sui_*fj.phase().rho();
}
const phasePairKey pairik
@ -326,7 +310,7 @@ birthByCoalescence
Pair<word>::compare(pDmdt_.find(pairik).key(), pairik)
);
*pDmdt_[pairik] += dmdtSign*fk.x()/v*Sui_*fi.phase().rho();
*pDmdt_[pairik] += dmdtSign*fk.x()/v*Sui_*fk.phase().rho();
}
sizeGroups_[i].shapeModelPtr()->addCoalescence(Sui_, fj, fk);
@ -334,41 +318,39 @@ birthByCoalescence
}
void Foam::diameterModels::populationBalanceModel::
deathByCoalescence
void Foam::diameterModels::populationBalanceModel::deathByCoalescence
(
const label i,
const label j
)
{
const sizeGroup& fi = sizeGroups_[i];
const sizeGroup& fj = sizeGroups_[j];
const sizeGroup& fi = sizeGroups()[i];
const sizeGroup& fj = sizeGroups()[j];
SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
Sp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
if (i != j)
{
SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
Sp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
}
}
void Foam::diameterModels::populationBalanceModel::
birthByBreakup
void Foam::diameterModels::populationBalanceModel::birthByBreakup
(
const label k,
const label model
)
{
const sizeGroup& fk = sizeGroups_[k];
const sizeGroup& fk = sizeGroups()[k];
for (label i = 0; i <= k; i++)
{
const sizeGroup& fi = sizeGroups_[i];
const sizeGroup& fi = sizeGroups()[i];
Sui_ =
fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i, k)
*fk*fk.phase()/fk.x();
fi.x()*breakupModels_[model].dsdPtr()().nik(i, k)/fk.x()
*breakupRate_()*fk*fk.phase();
Su_[i] += Sui_;
@ -385,7 +367,7 @@ birthByBreakup
Pair<word>::compare(pDmdt_.find(pair).key(), pair)
);
*pDmdt_[pair] += dmdtSign*Sui_*fi.phase().rho();
*pDmdt_[pair] += dmdtSign*Sui_*fk.phase().rho();
}
sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fk);
@ -395,20 +377,20 @@ birthByBreakup
void Foam::diameterModels::populationBalanceModel::deathByBreakup(const label i)
{
const sizeGroup& fi = sizeGroups_[i];
SuSp_[i] += breakupRate_()*fi.phase();
Sp_[i] += breakupRate_()*sizeGroups()[i].phase();
}
void Foam::diameterModels::populationBalanceModel::calcDeltas()
{
forAll(sizeGroups_, i)
forAll(sizeGroups(), i)
{
if (delta_[i].empty())
{
for (label j = 0; j <= sizeGroups_.size() - 1; j++)
for (label j = 0; j <= sizeGroups().size() - 1; j++)
{
const sizeGroup& fj = sizeGroups()[j];
delta_[i].append
(
new dimensionedScalar
@ -421,14 +403,14 @@ void Foam::diameterModels::populationBalanceModel::calcDeltas()
if
(
v_[i].value() < 0.5*sizeGroups_[j].x().value()
v_[i].value() < 0.5*fj.x().value()
&&
0.5*sizeGroups_[j].x().value() < v_[i+1].value()
0.5*fj.x().value() < v_[i+1].value()
)
{
delta_[i][j] = mag(0.5*sizeGroups_[j].x() - v_[i]);
delta_[i][j] = mag(0.5*fj.x() - v_[i]);
}
else if (0.5*sizeGroups_[j].x().value() <= v_[i].value())
else if (0.5*fj.x().value() <= v_[i].value())
{
delta_[i][j].value() = 0;
}
@ -438,17 +420,18 @@ void Foam::diameterModels::populationBalanceModel::calcDeltas()
}
void Foam::diameterModels::populationBalanceModel::
birthByBinaryBreakup
void Foam::diameterModels::populationBalanceModel::birthByBinaryBreakup
(
const label i,
const label j
)
{
const sizeGroup& fj = sizeGroups_[j];
const sizeGroup& fi = sizeGroups_[i];
const sizeGroup& fi = sizeGroups()[i];
const sizeGroup& fj = sizeGroups()[j];
Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*fj.phase()/fj.x();
const volScalarField Su(binaryBreakupRate_()*fj*fj.phase());
Sui_ = fi.x()*delta_[i][j]/fj.x()*Su;
Su_[i] += Sui_;
@ -467,7 +450,7 @@ birthByBinaryBreakup
Pair<word>::compare(pDmdt_.find(pairij).key(), pairij)
);
*pDmdt_[pairij] += dmdtSign*Sui_*fi.phase().rho();
*pDmdt_[pairij] += dmdtSign*Sui_*fj.phase().rho();
}
dimensionedScalar Eta;
@ -475,18 +458,15 @@ birthByBinaryBreakup
for (label k = 0; k <= j; k++)
{
// Calculate fraction for intra-interval events
Eta = eta(k, v);
if (Eta.value() == 0) continue;
const sizeGroup& fk = sizeGroups_[k];
const sizeGroup& fk = sizeGroups()[k];
volScalarField& Suk = Sui_;
Suk =
sizeGroups_[k].x()*binaryBreakupRate_()*delta_[i][j]*Eta
*fj*fj.phase()/fj.x();
Suk = fk.x()*delta_[i][j]*Eta/fj.x()*Su;
Su_[k] += Suk;
@ -507,24 +487,21 @@ birthByBinaryBreakup
)
);
*pDmdt_[pairkj] += dmdtSign*Suk*fi.phase().rho();
*pDmdt_[pairkj] += dmdtSign*Suk*fj.phase().rho();
}
sizeGroups_[i].shapeModelPtr()->addBreakup(Suk, fj);
sizeGroups_[k].shapeModelPtr()->addBreakup(Suk, fj);
}
}
void Foam::diameterModels::populationBalanceModel::
deathByBinaryBreakup
void Foam::diameterModels::populationBalanceModel::deathByBinaryBreakup
(
const label j,
const label i
)
{
const volScalarField& alphai = sizeGroups_[i].phase();
SuSp_[i] += alphai*binaryBreakupRate_()*delta_[j][i];
Sp_[i] += sizeGroups()[i].phase()*binaryBreakupRate_()*delta_[j][i];
}
@ -536,77 +513,18 @@ void Foam::diameterModels::populationBalanceModel::drift
{
model.addToDriftRate(driftRate_(), i);
const sizeGroup& fp = sizeGroups_[i];
const sizeGroup& fp = sizeGroups()[i];
if (i == 0)
if (i < sizeGroups().size() - 1)
{
rx_() = pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
}
else if (i == sizeGroups_.size() - 1)
{
rx_() = neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
}
else
{
rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x()
+ neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
}
SuSp_[i] +=
(neg(1 - rx_()) + neg(rx_() - rx_()/(1 - rx_())))*driftRate_()
*fp.phase()/((rx_() - 1)*fp.x());
rx_() = Zero;
rdx_() = Zero;
if (i == sizeGroups_.size() - 2)
{
rx_() = pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
rdx_() =
pos(driftRate_())
*(sizeGroups_[i+1].x() - sizeGroups_[i].x())
/(sizeGroups_[i].x() - sizeGroups_[i-1].x());
}
else if (i < sizeGroups_.size() - 2)
{
rx_() = pos(driftRate_())*sizeGroups_[i+2].x()/sizeGroups_[i+1].x();
rdx_() =
pos(driftRate_())
*(sizeGroups_[i+2].x() - sizeGroups_[i+1].x())
/(sizeGroups_[i+1].x() - sizeGroups_[i].x());
}
if (i == 1)
{
rx_() += neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
rdx_() +=
neg(driftRate_())
*(sizeGroups_[i].x() - sizeGroups_[i-1].x())
/(sizeGroups_[i+1].x() - sizeGroups_[i].x());
}
else if (i > 1)
{
rx_() += neg(driftRate_())*sizeGroups_[i-2].x()/sizeGroups_[i-1].x();
rdx_() +=
neg(driftRate_())
*(sizeGroups_[i-1].x() - sizeGroups_[i-2].x())
/(sizeGroups_[i].x() - sizeGroups_[i-1].x());
}
if (i != sizeGroups_.size() - 1)
{
const sizeGroup& fe = sizeGroups_[i+1];
const sizeGroup& fe = sizeGroups()[i+1];
volScalarField& Sue = Sui_;
Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_()*fp.phase();
Sue =
pos(driftRate_())*driftRate_()*rdx_()
*fp*fp.phase()/fp.x()
/(rx_() - 1);
fe.x()/(fp.x()*(fe.x() - fp.x()))
*pos(driftRate_())*driftRate_()*fp*fp.phase();
Su_[i+1] += Sue;
@ -629,15 +547,21 @@ void Foam::diameterModels::populationBalanceModel::drift
sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model);
}
if (i != 0)
if (i == sizeGroups().size() - 1)
{
const sizeGroup& fw = sizeGroups_[i-1];
Sp_[i] -= pos(driftRate_())*driftRate_()*fp.phase()/fp.x();
}
if (i > 0)
{
const sizeGroup& fw = sizeGroups()[i-1];
volScalarField& Suw = Sui_;
Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_()*fp.phase();
Suw =
neg(driftRate_())*driftRate_()*rdx_()
*fp*fp.phase()/fp.x()
/(rx_() - 1);
fw.x()/(fp.x()*(fw.x() - fp.x()))
*neg(driftRate_())*driftRate_()*fp*fp.phase();
Su_[i-1] += Suw;
@ -659,21 +583,25 @@ void Foam::diameterModels::populationBalanceModel::drift
sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model);
}
if (i == 0)
{
Sp_[i] -= neg(driftRate_())*driftRate_()*fp.phase()/fp.x();
}
}
void Foam::diameterModels::populationBalanceModel::
nucleation
void Foam::diameterModels::populationBalanceModel::nucleation
(
const label i,
nucleationModel& model
)
{
const sizeGroup& fi = sizeGroups_[i];
const sizeGroup& fi = sizeGroups()[i];
model.addToNucleationRate(nucleationRate_(), i);
Sui_ = sizeGroups_[i].x()*nucleationRate_();
Sui_ = fi.x()*nucleationRate_();
Su_[i] += Sui_;
@ -683,11 +611,11 @@ nucleation
void Foam::diameterModels::populationBalanceModel::sources()
{
forAll(sizeGroups_, i)
forAll(sizeGroups(), i)
{
sizeGroups_[i].shapeModelPtr()->reset();
Su_[i] = Zero;
SuSp_[i] = Zero;
Sp_[i] = Zero;
}
forAllConstIter
@ -700,84 +628,73 @@ void Foam::diameterModels::populationBalanceModel::sources()
*pDmdt_(phasePairIter()) = Zero;
}
forAll(sizeGroups_, i)
forAll(coalescencePairs_, coalescencePairi)
{
if (coalescence_.size() != 0)
label i = coalescencePairs_[coalescencePairi].first();
label j = coalescencePairs_[coalescencePairi].second();
coalescenceRate_() = Zero;
forAll(coalescenceModels_, model)
{
for (label j = 0; j <= i; j++)
{
coalescenceRate_() = Zero;
forAll(coalescence_, model)
{
coalescence_[model].addToCoalescenceRate
(
coalescenceRate_(),
i,
j
);
}
birthByCoalescence(i, j);
deathByCoalescence(i, j);
}
coalescenceModels_[model].addToCoalescenceRate
(
coalescenceRate_(),
i,
j
);
}
if (breakup_.size() != 0)
birthByCoalescence(i, j);
deathByCoalescence(i, j);
}
forAll(binaryBreakupPairs_, binaryBreakupPairi)
{
label i = binaryBreakupPairs_[binaryBreakupPairi].first();
label j = binaryBreakupPairs_[binaryBreakupPairi].second();
binaryBreakupRate_() = Zero;
forAll(binaryBreakupModels_, model)
{
forAll(breakup_, model)
{
breakup_[model].setBreakupRate(breakupRate_(), i);
birthByBreakup(i, model);
deathByBreakup(i);
}
binaryBreakupModels_[model].addToBinaryBreakupRate
(
binaryBreakupRate_(),
j,
i
);
}
if (binaryBreakup_.size() != 0)
birthByBinaryBreakup(j, i);
deathByBinaryBreakup(j, i);
}
forAll(sizeGroups(), i)
{
forAll(breakupModels_, model)
{
label j = 0;
breakupModels_[model].setBreakupRate(breakupRate_(), i);
while (delta_[j][i].value() != 0)
{
binaryBreakupRate_() = Zero;
birthByBreakup(i, model);
forAll(binaryBreakup_, model)
{
binaryBreakup_[model].addToBinaryBreakupRate
(
binaryBreakupRate_(),
j,
i
);
}
birthByBinaryBreakup(j, i);
deathByBinaryBreakup(j, i);
j++;
}
deathByBreakup(i);
}
if (drift_.size() != 0)
forAll(drift_, model)
{
forAll(drift_, model)
{
driftRate_() = Zero;
drift(i, drift_[model]);
}
driftRate_() = Zero;
drift(i, drift_[model]);
}
if (nucleation_.size() != 0)
forAll(nucleation_, model)
{
forAll(nucleation_, model)
{
nucleationRate_() = Zero;
nucleation(i, nucleation_[model]);
}
nucleationRate_() = Zero;
nucleation(i, nucleation_[model]);
}
}
}
@ -818,7 +735,7 @@ Foam::diameterModels::populationBalanceModel::calcDsm()
invDsm += max(phase, phase.residualAlpha())/(phase.d()*alphas_());
}
return 1.0/tInvDsm;
return 1/tInvDsm;
}
@ -883,7 +800,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
v_(),
delta_(),
Su_(),
SuSp_(),
Sp_(),
Sui_
(
IOobject
@ -895,32 +812,32 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
mesh_,
dimensionedScalar(inv(dimTime), Zero)
),
coalescence_
coalescenceModels_
(
dict_.lookup("coalescenceModels"),
coalescenceModel::iNew(*this)
),
coalescenceRate_(),
breakup_
coalescencePairs_(),
breakupModels_
(
dict_.lookup("breakupModels"),
breakupModel::iNew(*this)
),
breakupRate_(),
binaryBreakup_
binaryBreakupModels_
(
dict_.lookup("binaryBreakupModels"),
binaryBreakupModel::iNew(*this)
),
binaryBreakupRate_(),
binaryBreakupPairs_(),
drift_
(
dict_.lookup("driftModels"),
driftModel::iNew(*this)
),
driftRate_(),
rx_(),
rdx_(),
nucleation_
(
dict_.lookup("nucleationModels"),
@ -936,7 +853,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
this->createPhasePairs();
if (sizeGroups_.size() < 3)
if (sizeGroups().size() < 3)
{
FatalErrorInFunction
<< "The populationBalance " << name_
@ -945,7 +862,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
}
if (coalescence_.size() != 0)
if (coalescenceModels_.size() != 0)
{
coalescenceRate_.set
(
@ -961,9 +878,17 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
dimensionedScalar(dimVolume/dimTime, Zero)
)
);
forAll(sizeGroups(), i)
{
for (label j = 0; j <= i; j++)
{
coalescencePairs_.append(labelPair(i, j));
}
}
}
if (breakup_.size() != 0)
if (breakupModels_.size() != 0)
{
breakupRate_.set
(
@ -981,7 +906,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
);
}
if (binaryBreakup_.size() != 0)
if (binaryBreakupModels_.size() != 0)
{
binaryBreakupRate_.set
(
@ -1002,6 +927,19 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
)
)
);
calcDeltas();
forAll(sizeGroups(), i)
{
label j = 0;
while (delta_[j][i].value() != 0)
{
binaryBreakupPairs_.append(labelPair(i, j));
j++;
}
}
}
if (drift_.size() != 0)
@ -1020,36 +958,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
dimensionedScalar(dimVolume/dimTime, Zero)
)
);
rx_.set
(
new volScalarField
(
IOobject
(
"r",
fluid_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar(dimless, Zero)
)
);
rdx_.set
(
new volScalarField
(
IOobject
(
"r",
fluid_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar(dimless, Zero)
)
);
}
if (nucleation_.size() != 0)
@ -1158,17 +1066,17 @@ Foam::diameterModels::populationBalanceModel::eta
const dimensionedScalar& v
) const
{
const dimensionedScalar& x0 = sizeGroups_[0].x();
const dimensionedScalar& xi = sizeGroups_[i].x();
const dimensionedScalar& xm = sizeGroups_.last().x();
const dimensionedScalar& x0 = sizeGroups()[0].x();
const dimensionedScalar& xi = sizeGroups()[i].x();
const dimensionedScalar& xm = sizeGroups().last().x();
dimensionedScalar lowerBoundary(x0);
dimensionedScalar upperBoundary(xm);
if (i != 0) lowerBoundary = sizeGroups_[i-1].x();
if (i != 0) lowerBoundary = sizeGroups()[i-1].x();
if (i != sizeGroups_.size() - 1) upperBoundary = sizeGroups_[i+1].x();
if (i != sizeGroups().size() - 1) upperBoundary = sizeGroups()[i+1].x();
if ((i == 0 && v < x0) || (i == sizeGroups_.size() - 1 && v > xm))
if ((i == 0 && v < x0) || (i == sizeGroups().size() - 1 && v > xm))
{
return v/xi;
}
@ -1178,7 +1086,7 @@ Foam::diameterModels::populationBalanceModel::eta
}
else if (v.value() == xi.value())
{
return 1.0;
return 1;
}
else if (v > xi)
{
@ -1255,7 +1163,7 @@ void Foam::diameterModels::populationBalanceModel::solve()
maxInitialResidual = 0;
forAll(sizeGroups_, i)
forAll(sizeGroups(), i)
{
sizeGroup& fi = sizeGroups_[i];
const phaseModel& phase = fi.phase();
@ -1269,7 +1177,7 @@ void Foam::diameterModels::populationBalanceModel::solve()
+ fvm::div(phase.alphaPhi(), fi)
==
Su_[i]
- fvm::SuSp(SuSp_[i], fi)
- fvm::Sp(Sp_[i], fi)
+ fluid_.fvModels().source(alpha, rho, fi)/rho
+ fvc::ddt(residualAlpha, fi)
- fvm::ddt(residualAlpha, fi)
@ -1290,12 +1198,12 @@ void Foam::diameterModels::populationBalanceModel::solve()
volScalarField fAlpha0
(
sizeGroups_.first()*sizeGroups_.first().phase()
sizeGroups().first()*sizeGroups().first().phase()
);
volScalarField fAlphaN
(
sizeGroups_.last()*sizeGroups_.last().phase()
sizeGroups().last()*sizeGroups().last().phase()
);
Info<< this->name() << " sizeGroup phase fraction first, last = "

46
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
View File

@ -53,7 +53,7 @@ Description
\verbatim
Coalescence and breakup term formulation:
Kumar, S., & Ramkrishna, D. (1996).
On the solution of population balance equations by discretisation-I. A
On the solution of population balance equations by discretization-I. A
fixed pivot technique.
Chemical Engineering Science, 51(8), 1311-1332.
\endverbatim
@ -154,6 +154,7 @@ SourceFiles
#include "pimpleControl.H"
#include "phaseDynamicMomentumTransportModel.H"
#include "HashPtrTable.H"
#include "Pair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -233,42 +234,42 @@ private:
//- Explicitly treated sources
PtrList<volScalarField> Su_;
//- Sources treated implicitly or explicitly depending on sign
PtrList<volScalarField> SuSp_;
//- Implicitly treated sources
PtrList<volScalarField> Sp_;
//- Field for caching sources
volScalarField Sui_;
//- Coalescence models
PtrList<coalescenceModel> coalescence_;
PtrList<coalescenceModel> coalescenceModels_;
//- Coalescence rate
autoPtr<volScalarField> coalescenceRate_;
//- BreakupModels
PtrList<breakupModel> breakup_;
//- Coalescence relevant size group pairs
List<labelPair> coalescencePairs_;
//- Breakup models
PtrList<breakupModel> breakupModels_;
//- Breakup rate
autoPtr<volScalarField> breakupRate_;
//- Binary breakup models
PtrList<binaryBreakupModel> binaryBreakup_;
PtrList<binaryBreakupModel> binaryBreakupModels_;
//- Binary breakup rate
autoPtr<volScalarField> binaryBreakupRate_;
//- Binary breakup relevant size group pairs
List<labelPair> binaryBreakupPairs_;
//- Drift models
PtrList<driftModel> drift_;
//- Drift rate
autoPtr<volScalarField> driftRate_;
//- Ratio between successive representative volumes
autoPtr<volScalarField> rx_;
//- Ratio between successive class widths
autoPtr<volScalarField> rdx_;
//- Zeroeth order models
PtrList<nucleationModel> nucleation_;
@ -416,18 +417,17 @@ public:
//- Return the sizeGroups belonging to this populationBalance
inline const UPtrList<sizeGroup>& sizeGroups() const;
//- Return semi-implicit source terms
inline const volScalarField& SuSp(const label i) const;
//- Return list of unordered phasePairs in this populationBalance
inline const phasePairTable& phasePairs() const;
//- Returns true if both phases are velocity groups and
// belong to this populationBalance
inline bool isVelocityGroupPair(const phasePair& pair) const;
//- Return implicit source terms
inline const volScalarField& Sp(const label i) const;
//- Return the sizeGroup boundaries
inline const PtrList<dimensionedScalar>& v() const;
//- Return coalescence relevant size group pairs
inline const List<labelPair>& coalescencePairs() const;
//- Return binary breakup relevant size group pairs
inline const List<labelPair>& binaryBreakupPairs() const;
//- Return total void of phases belonging to this populationBalance
inline const volScalarField& alphas() const;
@ -435,6 +435,10 @@ public:
//- Return average velocity
inline const volVectorField& U() const;
//- Returns true if both phases are velocity groups and
// belong to this populationBalance
inline bool isVelocityGroupPair(const phasePair& pair) const;
//- Return allocation coefficient
const dimensionedScalar eta
(

71
applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
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-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -84,13 +84,6 @@ Foam::diameterModels::populationBalanceModel::sizeGroups() const
}
inline const Foam::volScalarField&
Foam::diameterModels::populationBalanceModel::SuSp(const label i) const
{
return SuSp_[i];
}
inline const Foam::diameterModels::populationBalanceModel::phasePairTable&
Foam::diameterModels::populationBalanceModel::phasePairs() const
{
@ -98,35 +91,24 @@ Foam::diameterModels::populationBalanceModel::phasePairs() const
}
inline bool Foam::diameterModels::populationBalanceModel::isVelocityGroupPair
(
const phasePair& pair
) const
inline const Foam::volScalarField&
Foam::diameterModels::populationBalanceModel::Sp(const label i) const
{
if
(
isA<velocityGroup>(pair.phase1().dPtr()())
&&
isA<velocityGroup>(pair.phase2().dPtr()())
)
{
const velocityGroup& velGroup1 =
refCast<const velocityGroup>(pair.phase1().dPtr()());
const velocityGroup& velGroup2 =
refCast<const velocityGroup>(pair.phase2().dPtr()());
return velGroup1.popBalName() == velGroup2.popBalName();
}
return false;
return Sp_[i];
}
inline const Foam::PtrList<Foam::dimensionedScalar>&
Foam::diameterModels::populationBalanceModel::v() const
inline const Foam::List<Foam::Pair<Foam::label>>&
Foam::diameterModels::populationBalanceModel::coalescencePairs() const
{
return v_;
return coalescencePairs_;
}
inline const Foam::List<Foam::Pair<Foam::label>>&
Foam::diameterModels::populationBalanceModel::binaryBreakupPairs() const
{
return binaryBreakupPairs_;
}
@ -158,4 +140,29 @@ Foam::diameterModels::populationBalanceModel::U() const
}
inline bool Foam::diameterModels::populationBalanceModel::isVelocityGroupPair
(
const phasePair& pair
) const
{
if
(
isA<velocityGroup>(pair.phase1().dPtr()())
&&
isA<velocityGroup>(pair.phase2().dPtr()())
)
{
const velocityGroup& velGroup1 =
refCast<const velocityGroup>(pair.phase1().dPtr()());
const velocityGroup& velGroup2 =
refCast<const velocityGroup>(pair.phase2().dPtr()());
return velGroup1.popBalName() == velGroup2.popBalName();
}
return false;
}
// ************************************************************************* //