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

reacting*EulerFoam: Separate continuity error from mass transfer

Browse Source 
continuityError is now just the transport inconsistency. Mass sources,
whether as a result of fvOptions or phase-transfer/change processes, are
not included.
This commit is contained in:
Will Bainbridge
2020-02-27 10:24:32 +00:00
parent 7867ffce1b
commit e3903fdb35
18 changed files with 101 additions and 129 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.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) 2015-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -135,33 +135,19 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::addDmdtUfs
phaseModel& phase1 = this->phases()[pair.phase1().name()];
phaseModel& phase2 = this->phases()[pair.phase2().name()];
// Note that the phase UEqn contains a continuity error term, which
// implicitly adds a mass transfer term of fvm::Sp(dmdt, U). These
// additions remove the part of this term corresponding to the supplied
// mass transfer rate; I.e.:
//
// DDt(alpha1*rho1*U1) + ... = dmdt21*U2 + dmdt12*U1
// DDt(alpha1*rho1*U1) - contErr + ... = dmdt21*U2 + dmdt12*U1 - contErr
// DDt(alpha1*rho1*U1) - contErr + ... = dmdt21*U2 + dmdt12*U1 - dmdt*U1
// DDt(alpha1*rho1*U1) - contErr + ... = dmdt21*U2 - dmdt21*U1
//
// Where contErr is the continuity error associated with this mass
// transfer rate, dmdt.
const volScalarField dmdtf21(posPart(dmdtf));
const volScalarField dmdtf12(negPart(dmdtf));
if (!phase1.stationary())
{
const volScalarField dmdtf21(posPart(dmdtf));
*eqns[phase1.name()] +=
dmdtf21*phase2.U() - fvm::Sp(dmdtf21, phase1.URef());
dmdtf21*phase2.U() + fvm::Sp(dmdtf12, phase1.URef());
}
if (!phase2.stationary())
{
const volScalarField dmdtf12(negPart(dmdtf));
*eqns[phase2.name()] -=
dmdtf12*phase1.U() - fvm::Sp(dmdtf12, phase2.URef());
dmdtf12*phase1.U() + fvm::Sp(dmdtf21, phase2.URef());
}
}
}

26
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.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) 2015-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -55,17 +55,13 @@ void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Note that the phase EEqn contains a continuity error term. See
// MomentumTransferPhaseSystem::addDmdtUfs for an explanation of the
// fvm::Sp terms below.
// Transfer of energy from bulk to bulk
*eqns[phase1.name()] += dmdtf21*he2 - fvm::Sp(dmdtf21, he1);
*eqns[phase2.name()] -= dmdtf12*he1 - fvm::Sp(dmdtf12, he2);
*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 - K1);
*eqns[phase2.name()] -= dmdtf12*(K1 - K2);
*eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2;
}
}
@ -91,10 +87,6 @@ void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Note that the phase EEqn contains a continuity error term. See
// MomentumTransferPhaseSystem::addDmdtUfs for an explanation of the
// fvm::Sp terms below.
forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
{
const word& member = dmidtfJter.key();
@ -135,12 +127,12 @@ void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
}
// Transfer of energy from bulk to bulk
*eqns[phase1.name()] += dmidtf21*hei2 - fvm::Sp(dmidtf21, he1);
*eqns[phase2.name()] -= dmidtf12*hei1 - fvm::Sp(dmidtf12, he2);
*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 - K1);
*eqns[phase2.name()] -= dmidtf12*(K1 - K2);
*eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1;
*eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2;
}
}
}

8
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.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) 2018-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2018-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -92,9 +92,6 @@ void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addDmdtYfs
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
// Note that the phase YiEqn does not contain a continuity error term,
// so the transfers below are complete.
forAll(phase1.Y(), Yi1)
{
const volScalarField& Y1 = phase1.Y()[Yi1];
@ -122,9 +119,6 @@ void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addDmidtYf
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
// Note that the phase YiEqn does not contain a continuity error term,
// so the transfers below are complete.
forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
{
const word& member = dmidtfJter.key();

3
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
View File

@ -447,9 +447,6 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::specieTransfer() const
const volScalarField& Y1 = phase1.Y(volatile_);
const volScalarField& Y2 = phase2.Y(volatile_);
// Note that the phase YiEqn does not contain a continuity error
// term, so these additions represent the entire mass transfer
const volScalarField dmdtf(this->totalDmdtf(pair));
*eqns[Y1.name()] += dmdtf;

42
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.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) 2015-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -80,10 +80,6 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Note that the phase EEqn contains a continuity error term. See
// MomentumTransferPhaseSystem::addDmdtU for an explanation of the
// fvm::Sp terms below.
if (heatTransferModels_.found(key))
{
// Assume a thermally-coupled mass transfer process. Calculate
@ -94,8 +90,8 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
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 - fvm::Sp(dmdtf, he1);
*eqns[phase2.name()] -= dmdtf*hef2 - fvm::Sp(dmdtf, he2);
*eqns[phase1.name()] += dmdtf*hef1;
*eqns[phase2.name()] -= dmdtf*hef2;
// Latent heat contribution
const volScalarField L(hef2 + thermo2.hc() - hef1 - thermo1.hc());
@ -106,8 +102,8 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
*eqns[phase2.name()] += (1 - H1Fac)*dmdtf*L;
// Transfer of kinetic energy
*eqns[phase1.name()] += dmdtf21*(K2 - K1);
*eqns[phase2.name()] -= dmdtf12*(K1 - K2);
*eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2;
}
else
{
@ -116,12 +112,12 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
// change and therefore no interface state or latent heat...
// Transfer of energy from bulk to bulk
*eqns[phase1.name()] += dmdtf21*he2 - fvm::Sp(dmdtf21, he1);
*eqns[phase2.name()] -= dmdtf12*he1 - fvm::Sp(dmdtf12, he2);
*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 - K1);
*eqns[phase2.name()] -= dmdtf12*(K1 - K2);
*eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1;
*eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2;
}
}
}
@ -150,10 +146,6 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
// Note that the phase EEqn contains a continuity error term. See
// MomentumTransferPhaseSystem::addDmdtU for an explanation of the
// fvm::Sp terms below.
if (heatTransferModels_.found(key))
{
// Assume a thermally-coupled mass transfer process. Calculate
@ -233,16 +225,16 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
const volScalarField Li(hefi2 + hci2 - hefi1 - hci1);
// Transfer of energy from the interface into the bulk
*eqns[phase1.name()] += dmidtf*hefi1 - fvm::Sp(dmidtf, he1);
*eqns[phase2.name()] -= dmidtf*hefi2 - fvm::Sp(dmidtf, he2);
*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 - K1);
*eqns[phase2.name()] -= dmidtf12*(K1 - K2);
*eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1;
*eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2;
}
}
else
@ -297,12 +289,12 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
}
// Transfer of energy from bulk to bulk
*eqns[phase1.name()] += dmidtf21*hei2 - fvm::Sp(dmidtf21, he1);
*eqns[phase2.name()] -= dmidtf12*hei1 - fvm::Sp(dmidtf12, he2);
*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 - K1);
*eqns[phase2.name()] -= dmidtf12*(K1 - K2);
*eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1;
*eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2;
}
}
}

17
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C
View File

@ -215,23 +215,6 @@ void Foam::diameterModels::shapeModels::fractal::correct()
+ alpha*rho*fi*fvm::ddt(kappa_)
+ fvm::div(fAlphaRhoPhi, kappa_)
+ fvm::SuSp(-phase.continuityError()*fi, kappa_)
+ fvm::SuSp
(
(
fi.VelocityGroup().dmdt()
// Temporary term pending update of continuityError
- (
phase.fluid().fvOptions()
(
alpha,
const_cast<volScalarField&>(rho)
) & rho
)
)
*fi,
kappa_
)
==
- sinteringModel_->R()
+ fvc::Su(Su_*rho, kappa_)

17
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
View File

@ -203,13 +203,14 @@ Foam::MovingPhaseModel<BasePhaseModel>::~MovingPhaseModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasePhaseModel>
void Foam::MovingPhaseModel<BasePhaseModel>::correctContinuityError()
void Foam::MovingPhaseModel<BasePhaseModel>::correctContinuityError
(
const volScalarField& source
)
{
volScalarField& rho = this->thermoRef().rho();
continuityError_ =
fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_)
- (this->fluid().fvOptions()(*this, rho)&rho);
continuityError_ = fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_) - source;
}
@ -218,7 +219,6 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correct()
{
BasePhaseModel::correct();
this->fluid().MRF().correctBoundaryVelocity(U_);
correctContinuityError();
}
@ -247,13 +247,6 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correctKinematics()
}
template<class BasePhaseModel>
void Foam::MovingPhaseModel<BasePhaseModel>::correctThermo()
{
correctContinuityError();
}
template<class BasePhaseModel>
void Foam::MovingPhaseModel<BasePhaseModel>::correctTurbulence()
{

9
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
View File

@ -103,9 +103,6 @@ private:
//- Calculate and return the flux field
tmp<surfaceScalarField> phi(const volVectorField& U) const;
//- Correct the continuity error
void correctContinuityError();
public:
@ -128,12 +125,12 @@ public:
//- Correct the phase properties other than the thermo and turbulence
virtual void correct();
//- Correct the continuity error
virtual void correctContinuityError(const volScalarField& source);
//- Correct the kinematics
virtual void correctKinematics();
//- Correct the thermodynamics
virtual void correctThermo();
//- Correct the turbulence
virtual void correctTurbulence();

6
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.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) 2015-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -143,6 +143,10 @@ void Foam::phaseModel::correct()
}
void Foam::phaseModel::correctContinuityError(const volScalarField& source)
{}
void Foam::phaseModel::correctKinematics()
{}

3
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
View File

@ -186,6 +186,9 @@ public:
//- Correct the phase properties
virtual void correct();
//- Correct the continuity error
virtual void correctContinuityError(const volScalarField& source);
//- Correct the kinematics
virtual void correctKinematics();

35
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
View File

@ -450,6 +450,41 @@ void Foam::phaseSystem::correct()
}
void Foam::phaseSystem::correctContinuityError()
{
const PtrList<volScalarField> dmdts = this->dmdts();
forAll(movingPhaseModels_, movingPhasei)
{
phaseModel& phase = movingPhaseModels_[movingPhasei];
const volScalarField& alpha = phase;
volScalarField& rho = phase.thermoRef().rho();
volScalarField source
(
volScalarField::New
(
IOobject::groupName("source", phase.name()),
mesh_,
dimensionedScalar(dimDensity/dimTime, 0)
)
);
if (fvOptions().appliesToField(rho.name()))
{
source += fvOptions()(alpha, rho)&rho;
}
if (dmdts.set(phase.index()))
{
source += dmdts[phase.index()];
}
phase.correctContinuityError(source);
}
}
void Foam::phaseSystem::correctKinematics()
{
bool updateDpdt = false;

8
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
View File

@ -187,10 +187,7 @@ protected:
) const;
//- Generate pairs
void generatePairs
(
const dictTable& modelDicts
);
void generatePairs(const dictTable& modelDicts);
//- Generate pairs and sub-model tables
template<class modelType>
@ -523,6 +520,9 @@ public:
//- Correct the fluid properties other than those listed below
virtual void correct();
//- Correct the continuity errors
virtual void correctContinuityError();
//- Correct the kinematics
virtual void correctKinematics();

16
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
View File

@ -1304,21 +1304,7 @@ void Foam::diameterModels::populationBalanceModel::solve()
phase.alphaRhoPhi(),
fi
)
+ fvm::SuSp
(
fi.VelocityGroup().dmdt()
- phase.continuityError()
// Temporary term pending update of continuityError
- (
phase.fluid().fvOptions()
(
alpha,
const_cast<volScalarField&>(rho)
) & rho
),
fi
)
+ fvm::SuSp(-phase.continuityError(), fi)
==
fvc::Su(Su_[i]*rho, fi)
- fvm::SuSp(SuSp_[i]*rho, fi)

1
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H
View File

@ -33,6 +33,7 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
}
fluid.correctThermo();
fluid.correctContinuityError();
}

1
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.C
View File

@ -106,6 +106,7 @@ int main(int argc, char *argv[])
{
fluid.solve(rAUs, rAUfs);
fluid.correct();
fluid.correctContinuityError();
#include "YEqns.H"

1
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H
View File

@ -43,4 +43,5 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
}
fluid.correctThermo();
fluid.correctContinuityError();
}

1
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
View File

@ -102,6 +102,7 @@ int main(int argc, char *argv[])
{
fluid.solve(rAUs, rAUfs);
fluid.correct();
fluid.correctContinuityError();
#include "YEqns.H"

12
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.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) 2013-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -125,8 +125,14 @@ void Foam::diameterModels::IATE::correct()
(
((1.0/3.0)/alphaAv)
*(
fvc::ddt(phase()) + fvc::div(phase().alphaPhi())
- phase().continuityError()/phase().rho()
(
fvc::ddt(phase())
+ fvc::div(phase().alphaPhi())
)
- (
fvc::ddt(phase(), phase().rho()())
+ fvc::div(phase().alphaRhoPhi())
)/phase().rho()
),
kappai_
)