zhyang-dev/OpenFOAM-6
1
0
Fork 0
mirror of https://github.com/OpenFOAM/OpenFOAM-6.git synced 2025-12-08 06:57:46 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

reactingMultiphaseEulerFoam: Stationary phase

Browse Source 
Two new phase models have been added as selectable options for
reactingMultiphaseEulerFoam; pureStationaryPhaseModel and
pureStationaryIsothermalPhaseModel. These phases do not store a
velocity and their phase fractions remain constant throughout the
simulation. They are intended for use in modelling static particle beds
and other forms of porous media by means of the existing Euler-Euler
transfer models (drag, heat transfer, etc...).

Note that this functionality has not been extended to
reactingTwoPhaseEulerFoam, or the non-reacting *EulerFoam solvers.

Additional maintenance work has been carried out on the phase model
and phase system structure. The system can now loop over subsets of
phases with specific functionality (moving, multi-component, etc...) in
order to avoid testing for the existence of equations or variables in
the top level solver. The mass transfer handling and it's effect on
per-phase source terms has been refactored to reduce duplication. Const
and non-const access to phase properties has been formalised by renaming
non-const accessors with a "Ref" suffix, which is consistent with other
recent developments to classes including tmp and GeometricField, among
others. More sub-modelling details have been made private in order to
reduce the size of interfaces and improve abstraction.

This work was supported by Zhen Li, at Evonik
This commit is contained in:
Will Bainbridge
2018-03-22 11:36:43 +00:00
parent f2cc03bf8d
commit e352828514
79 changed files with 3495 additions and 1955 deletions
Download Patch File Download Diff File Expand all files Collapse all files

81
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
View File

@ -131,86 +131,6 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
bool
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::transfersMass() const
{
return true;
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
(
const phasePairKey& key
) const
{
return tmp<Foam::volScalarField>
(
new volScalarField
(
IOobject
(
"dmdt",
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
(
const Foam::phaseModel& phase
) const
{
tmp<volScalarField> tDmdt
(
new volScalarField
(
IOobject
(
IOobject::groupName("dmdt", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
return tDmdt;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
{
autoPtr<phaseSystem::momentumTransferTable>
eqnsPtr(BasePhaseSystem::momentumTransfer());
return eqnsPtr;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransferf() const
{
autoPtr<phaseSystem::momentumTransferTable>
eqnsPtr(BasePhaseSystem::momentumTransferf());
return eqnsPtr;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::heatTransferTable>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
@ -328,7 +248,6 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::massTransfer() const
template<class BasePhaseSystem>
void Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::correctThermo()
{
phaseSystem::correctThermo();
forAllConstIter

17
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
View File

@ -113,23 +113,6 @@ public:
// Member Functions
//- Return true if there is mass transfer
virtual bool transfersMass() const;
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
//- Return the momentum transfer matrices
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransfer() const;
//- Return the momentum transfer matrices for the face based algorithm
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransferf() const;
//- Return the heat transfer matrices
virtual autoPtr<phaseSystem::heatTransferTable>
heatTransfer() const;

61
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
View File

@ -49,67 +49,6 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::~HeatTransferPhaseSystem()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
bool Foam::HeatTransferPhaseSystem<BasePhaseSystem>::transfersMass() const
{
return false;
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
(
const phasePairKey& key
) const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
IOobject::groupName
(
"dmdt",
this->phasePairs_[key]->name()
),
this->mesh().time().timeName(),
this->mesh().time()
),
this->mesh(),
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
(
const Foam::phaseModel& phase
) const
{
tmp<volScalarField> tDmdt
(
new volScalarField
(
IOobject
(
IOobject::groupName("dmdt", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
return tDmdt;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::heatTransferTable>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const

9
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H
View File

@ -89,15 +89,6 @@ public:
// Member Functions
//- Return true if there is mass transfer
virtual bool transfersMass() const;
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
//- Return the heat transfer matrices
virtual autoPtr<phaseSystem::heatTransferTable> heatTransfer() const;

172
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
View File

@ -30,37 +30,24 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class BasePhaseSystem>
void Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const
Foam::tmp<Foam::volScalarField>
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::iDmdt
(
const phasePairKey& key
) const
{
// Source term due to mass transfer
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
if (!iDmdt_.found(key))
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
const volVectorField& U1(pair.phase1().U());
const volVectorField& U2(pair.phase2().U());
const volScalarField dmdt(this->iDmdt(pair));
const volScalarField dmdt21(posPart(dmdt));
const volScalarField dmdt12(negPart(dmdt));
*eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1);
*eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2);
return phaseSystem::dmdt(key);
}
const scalar dmdtSign(Pair<word>::compare(iDmdt_.find(key).key(), key));
return dmdtSign**iDmdt_[key];
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
@ -142,140 +129,31 @@ Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::iDmdt
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::dmdt
(
const phasePairKey& key
) const
{
const scalar dmdtSign(Pair<word>::compare(iDmdt_.find(key).key(), key));
return dmdtSign**iDmdt_[key];
return BasePhaseSystem::dmdt(key) + this->iDmdt(key);
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::iDmdt
(
const Foam::phaseModel& phase
) const
Foam::Xfer<Foam::PtrList<Foam::volScalarField>>
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::dmdts() const
{
tmp<volScalarField> tiDmdt
(
new volScalarField
(
IOobject
(
IOobject::groupName("iDmdt", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
forAllConstIter(iDmdtTable, iDmdt_, iDmdtIter)
{
const phasePair& pair(phasePairIter());
const phasePair& pair = this->phasePairs_[iDmdtIter.key()];
const volScalarField& iDmdt = *iDmdtIter();
if (pair.ordered())
{
continue;
}
if (pair.contains(phase))
{
tiDmdt.ref() += this->iDmdt
(
phasePairKey(phase.name(), pair.otherPhase(phase).name(), false)
);
}
this->addField(pair.phase1(), "dmdt", iDmdt, dmdts);
this->addField(pair.phase2(), "dmdt", - iDmdt, dmdts);
}
return tiDmdt;
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::dmdt
(
const Foam::phaseModel& phase
) const
{
tmp<volScalarField> tDmdt
(
new volScalarField
(
IOobject
(
IOobject::groupName("dmdt", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
if (pair.contains(phase))
{
tDmdt.ref() += this->dmdt
(
phasePairKey(phase.name(), pair.otherPhase(phase).name(), false)
);
}
}
return tDmdt;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
momentumTransfer() const
{
autoPtr<phaseSystem::momentumTransferTable>
eqnsPtr(BasePhaseSystem::momentumTransfer());
addMomentumTransfer(eqnsPtr());
return eqnsPtr;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
momentumTransferf() const
{
autoPtr<phaseSystem::momentumTransferTable>
eqnsPtr(BasePhaseSystem::momentumTransferf());
addMomentumTransfer(eqnsPtr());
return eqnsPtr;
return dmdts.xfer();
}
@ -315,8 +193,8 @@ heatTransfer() const
const volScalarField& he1(phase1.thermo().he());
const volScalarField& he2(phase2.thermo().he());
const volScalarField& K1(phase1.K());
const volScalarField& K2(phase2.K());
const volScalarField K1(phase1.K());
const volScalarField K2(phase2.K());
const volScalarField dmdt(this->dmdt(pair));
const volScalarField dmdt21(posPart(dmdt));

42
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H
View File

@ -69,6 +69,9 @@ protected:
phasePairKey::hash
> interfaceCompositionModelTable;
typedef HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
iDmdtTable;
// Protected data
@ -78,20 +81,16 @@ protected:
interfaceCompositionModelTable interfaceCompositionModels_;
//- Interfacial Mass transfer rate
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
iDmdt_;
iDmdtTable iDmdt_;
//- Explicit part of the mass transfer rate
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
iDmdtExplicit_;
iDmdtTable iDmdtExplicit_;
// Protected member functions
//- Add the mass-based momentum transfer to the equations
void addMomentumTransfer
(
phaseSystem::momentumTransferTable& eqns
) const;
//- Return the interfacial mass transfer rate for a pair for a pair
virtual tmp<volScalarField> iDmdt(const phasePairKey& key) const;
public:
@ -108,28 +107,11 @@ public:
// Member Functions
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> iDmdt(const phasePairKey& key) const;
//- Return the mass transfer rate for a pair
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> iDmdt(const phaseModel& phase) const;
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const
{
return BasePhaseSystem::dmdt(key) + this->iDmdt(key);
};
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
//- Return the momentum transfer matrices
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransfer() const;
//- Return the momentum transfer matrices for the face based algorithm
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransferf() const;
//- Return the mass transfer rates for each phase
virtual Xfer<PtrList<volScalarField>> dmdts() const;
//- Return the heat transfer matrices
virtual autoPtr<phaseSystem::heatTransferTable> heatTransfer() const;

439
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
View File

@ -44,6 +44,150 @@ License
#include "fvcSnGrad.H"
#include "fvMatrix.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kd
(
const phasePairKey& key
) const
{
if (dragModels_.found(key))
{
return dragModels_[key]->K();
}
else
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
dragModel::typeName + ":K",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar("zero", dragModel::dimK, 0)
)
);
}
}
template<class BasePhaseSystem>
Foam::tmp<Foam::surfaceScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kdf
(
const phasePairKey& key
) const
{
if (dragModels_.found(key))
{
return dragModels_[key]->Kf();
}
else
{
return tmp<surfaceScalarField>
(
new surfaceScalarField
(
IOobject
(
dragModel::typeName + ":K",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar("zero", dragModel::dimK, 0)
)
);
}
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vm
(
const phasePairKey& key
) const
{
if (virtualMassModels_.found(key))
{
return virtualMassModels_[key]->K();
}
else
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
virtualMassModel::typeName + ":K",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar("zero", virtualMassModel::dimK, 0)
)
);
}
}
template<class BasePhaseSystem>
void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::
addMassTransferMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const
{
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
const volScalarField dmdt(this->dmdt(pair));
if (!pair.phase1().stationary())
{
fvVectorMatrix& eqn = *eqns[pair.phase1().name()];
const volScalarField dmdt21(posPart(dmdt));
eqn += dmdt21*pair.phase2().U() - fvm::Sp(dmdt21, eqn.psi());
}
if (!pair.phase2().stationary())
{
fvVectorMatrix& eqn = *eqns[pair.phase2().name()];
const volScalarField dmdt12(negPart(dmdt));
eqn -= dmdt12*pair.phase1().U() - fvm::Sp(dmdt12, eqn.psi());
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
@ -157,111 +301,9 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kd
(
const phasePairKey& key
) const
{
if (dragModels_.found(key))
{
return dragModels_[key]->K();
}
else
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
dragModel::typeName + ":K",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar("zero", dragModel::dimK, 0)
)
);
}
}
template<class BasePhaseSystem>
Foam::tmp<Foam::surfaceScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kdf
(
const phasePairKey& key
) const
{
if (dragModels_.found(key))
{
return dragModels_[key]->Kf();
}
else
{
return tmp<surfaceScalarField>
(
new surfaceScalarField
(
IOobject
(
dragModel::typeName + ":K",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar("zero", dragModel::dimK, 0)
)
);
}
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vm
(
const phasePairKey& key
) const
{
if (virtualMassModels_.found(key))
{
return virtualMassModels_[key]->K();
}
else
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
virtualMassModel::typeName + ":K",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar("zero", virtualMassModel::dimK, 0)
)
);
}
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
{
// Create a momentum transfer matrix for each phase
autoPtr<phaseSystem::momentumTransferTable> eqnsPtr
@ -271,9 +313,9 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
phaseSystem::momentumTransferTable& eqns = eqnsPtr();
forAll(this->phaseModels_, phasei)
forAll(this->movingPhases(), movingPhasei)
{
const phaseModel& phase = this->phaseModels_[phasei];
const phaseModel& phase = this->movingPhases()[movingPhasei];
eqns.insert
(
@ -295,19 +337,19 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
}
// Add the implicit part of the drag force
forAllConstIter
(
phaseSystem::KdTable,
Kds_,
KdIter
)
forAllConstIter(KdTable, Kds_, KdIter)
{
const volScalarField& K(*KdIter());
const phasePair& pair(this->phasePairs_[KdIter.key()]);
forAllConstIter(phasePair, pair, iter)
{
*eqns[iter().name()] -= fvm::Sp(K, iter().U()());
if (!iter().stationary())
{
fvVectorMatrix& eqn = *eqns[iter().name()];
eqn -= fvm::Sp(K, eqn.psi());
}
}
}
@ -324,12 +366,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
}
// Add the virtual mass force
forAllConstIter
(
phaseSystem::VmTable,
Vms_,
VmIter
)
forAllConstIter(VmTable, Vms_, VmIter)
{
const volScalarField& Vm(*VmIter());
const phasePair& pair(this->phasePairs_[VmIter.key()]);
@ -339,28 +376,36 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
const phaseModel& phase = iter();
const phaseModel& otherPhase = iter.otherPhase();
const volVectorField& U = phase.U();
const surfaceScalarField& phi = phase.phi();
if (!phase.stationary())
{
fvVectorMatrix& eqn = *eqns[phase.name()];
*eqns[phase.name()] -=
Vm
*(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
- otherPhase.DUDt()
)
+ this->MRF_.DDt(Vm, U - otherPhase.U());
const volVectorField& U = eqn.psi();
const surfaceScalarField& phi = phase.phi();
eqn -=
Vm
*(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
- otherPhase.DUDt()
)
+ this->MRF_.DDt(Vm, U - otherPhase.U());
}
}
}
// Add the source term due to mass trasfer
addMassTransferMomentumTransfer(eqns);
return eqnsPtr;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf() const
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf()
{
// Create a momentum transfer matrix for each phase
autoPtr<phaseSystem::momentumTransferTable> eqnsPtr
@ -370,9 +415,9 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf() const
phaseSystem::momentumTransferTable& eqns = eqnsPtr();
forAll(this->phaseModels_, phasei)
forAll(this->movingPhases(), movingPhasei)
{
const phaseModel& phase = this->phaseModels_[phasei];
const phaseModel& phase = this->movingPhases()[movingPhasei];
eqns.insert
(
@ -401,12 +446,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf() const
}
// Add the virtual mass force
forAllConstIter
(
phaseSystem::VmTable,
Vms_,
VmIter
)
forAllConstIter(VmTable, Vms_, VmIter)
{
const volScalarField& Vm(*VmIter());
const phasePair& pair(this->phasePairs_[VmIter.key()]);
@ -416,15 +456,21 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf() const
const phaseModel& phase = iter();
const phaseModel& otherPhase = iter.otherPhase();
*eqns[phase.name()] -=
Vm
*(
UgradUs[phase.index()]
- (UgradUs[otherPhase.index()] & otherPhase.U())
);
if (!phase.stationary())
{
*eqns[phase.name()] -=
Vm
*(
UgradUs[phase.index()]
- (UgradUs[otherPhase.index()] & otherPhase.U())
);
}
}
}
// Add the source term due to mass trasfer
addMassTransferMomentumTransfer(eqns);
return eqnsPtr;
}
@ -436,12 +482,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::AFfs() const
PtrList<surfaceScalarField> AFfs(this->phaseModels_.size());
// Add the implicit part of the drag force
forAllConstIter
(
phaseSystem::KdfTable,
Kdfs_,
KdfIter
)
forAllConstIter(KdfTable, Kdfs_, KdfIter)
{
const surfaceScalarField& Kf(*KdfIter());
const phasePair& pair(this->phasePairs_[KdfIter.key()]);
@ -453,12 +494,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::AFfs() const
}
// Add the implicit part of the virtual mass force
forAllConstIter
(
phaseSystem::VmfTable,
Vmfs_,
VmfIter
)
forAllConstIter(VmfTable, Vmfs_, VmfIter)
{
const surfaceScalarField& Vmf(*VmfIter());
const phasePair& pair(this->phasePairs_[VmfIter.key()]);
@ -629,12 +665,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiFfs
PtrList<surfaceScalarField> phiFfs(this->phaseModels_.size());
// Add the explicit part of the virtual mass force
forAllConstIter
(
phaseSystem::VmfTable,
Vmfs_,
VmfIter
)
forAllConstIter(VmfTable, Vmfs_, VmfIter)
{
const surfaceScalarField& Vmf(*VmfIter());
const phasePair& pair(this->phasePairs_[VmfIter.key()]);
@ -797,12 +828,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiKdPhis
PtrList<surfaceScalarField> phiKdPhis(this->phaseModels_.size());
// Add the explicit part of the drag force
forAllConstIter
(
phaseSystem::KdTable,
Kds_,
KdIter
)
forAllConstIter(KdTable, Kds_, KdIter)
{
const volScalarField& K(*KdIter());
const phasePair& pair(this->phasePairs_[KdIter.key()]);
@ -844,12 +870,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiKdPhifs
PtrList<surfaceScalarField> phiKdPhifs(this->phaseModels_.size());
// Add the explicit part of the drag force
forAllConstIter
(
phaseSystem::KdfTable,
Kdfs_,
KdfIter
)
forAllConstIter(KdfTable, Kdfs_, KdfIter)
{
const surfaceScalarField& Kf(*KdfIter());
const phasePair& pair(this->phasePairs_[KdfIter.key()]);
@ -891,12 +912,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdUByAs
PtrList<volVectorField> KdUByAs(this->phaseModels_.size());
// Add the explicit part of the drag force
forAllConstIter
(
phaseSystem::KdTable,
Kds_,
KdIter
)
forAllConstIter(KdTable, Kds_, KdIter)
{
const volScalarField& K(*KdIter());
const phasePair& pair(this->phasePairs_[KdIter.key()]);
@ -991,12 +1007,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::ddtCorrByAs
// Add virtual mass correction
if (includeVirtualMass)
{
forAllConstIter
(
phaseSystem::VmTable,
Vms_,
VmIter
)
forAllConstIter(VmTable, Vms_, VmIter)
{
const volScalarField& Vm(*VmIter());
const phasePair& pair(this->phasePairs_[VmIter.key()]);
@ -1056,12 +1067,7 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::partialElimination
);
}
forAllConstIter
(
phaseSystem::KdTable,
Kds_,
KdIter
)
forAllConstIter(KdTable, Kds_, KdIter)
{
const volScalarField& K(*KdIter());
const phasePair& pair(this->phasePairs_[KdIter.key()]);
@ -1138,21 +1144,27 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::partialElimination
// Solve for the velocities and fluxes
for (label i = 1; i < phases.size(); ++ i)
{
for (label j = 0; j < i; j ++)
if (!phases[i].stationary())
{
phases[i].U() -= KdByAs[i][j]*phases[j].U();
phases[i].phi() -= phiKds[i][j]*phases[j].phi();
for (label j = 0; j < i; j ++)
{
phases[i].URef() -= KdByAs[i][j]*phases[j].U();
phases[i].phiRef() -= phiKds[i][j]*phases[j].phi();
}
}
}
for (label i = phases.size() - 1; i >= 0; i --)
{
for (label j = phases.size() - 1; j > i; j --)
if (!phases[i].stationary())
{
phases[i].U() -= KdByAs[i][j]*phases[j].U();
phases[i].phi() -= phiKds[i][j]*phases[j].phi();
for (label j = phases.size() - 1; j > i; j --)
{
phases[i].URef() -= KdByAs[i][j]*phases[j].U();
phases[i].phiRef() -= phiKds[i][j]*phases[j].phi();
}
phases[i].URef() /= KdByAs[i][i];
phases[i].phiRef() /= phiKds[i][i];
}
phases[i].U() /= KdByAs[i][i];
phases[i].phi() /= phiKds[i][i];
}
}
@ -1179,15 +1191,10 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::partialEliminationf
);
}
forAllConstIter
(
phaseSystem::KdfTable,
Kdfs_,
KdIter
)
forAllConstIter(KdfTable, Kdfs_, KdfIter)
{
const surfaceScalarField& K(*KdIter());
const phasePair& pair(this->phasePairs_[KdIter.key()]);
const surfaceScalarField& K(*KdfIter());
const phasePair& pair(this->phasePairs_[KdfIter.key()]);
const label phase1i = pair.phase1().index();
const label phase2i = pair.phase2().index();
@ -1239,18 +1246,24 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::partialEliminationf
// Solve for the fluxes
for (label i = 1; i < phases.size(); ++ i)
{
for (label j = 0; j < i; j ++)
if (!phases[i].stationary())
{
phases[i].phi() -= phiKdfs[i][j]*phases[j].phi();
for (label j = 0; j < i; j ++)
{
phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi();
}
}
}
for (label i = phases.size() - 1; i >= 0; i --)
{
for (label j = phases.size() - 1; j > i; j --)
if (!phases[i].stationary())
{
phases[i].phi() -= phiKdfs[i][j]*phases[j].phi();
for (label j = phases.size() - 1; j > i; j --)
{
phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi();
}
phases[i].phiRef() /= phiKdfs[i][i];
}
phases[i].phi() /= phiKdfs[i][i];
}
}

47
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
View File

@ -70,6 +70,34 @@ protected:
// Protected typedefs
typedef HashPtrTable
<
volScalarField,
phasePairKey,
phasePairKey::hash
> KdTable;
typedef HashPtrTable
<
surfaceScalarField,
phasePairKey,
phasePairKey::hash
> KdfTable;
typedef HashPtrTable
<
volScalarField,
phasePairKey,
phasePairKey::hash
> VmTable;
typedef HashPtrTable
<
surfaceScalarField,
phasePairKey,
phasePairKey::hash
> VmfTable;
typedef HashTable
<
autoPtr<BlendedInterfacialModel<dragModel>>,
@ -111,16 +139,16 @@ private:
// Private data
//- Drag coefficients
phaseSystem::KdTable Kds_;
KdTable Kds_;
//- Face drag coefficients
phaseSystem::KdfTable Kdfs_;
KdfTable Kdfs_;
//- Virtual mass coefficients
phaseSystem::VmTable Vms_;
VmTable Vms_;
//- Face virtual mass coefficients
phaseSystem::VmfTable Vmfs_;
VmfTable Vmfs_;
//- The phase diffusivities divided by the momentum coefficients
HashPtrTable<surfaceScalarField> DByAfs_;
@ -154,6 +182,11 @@ private:
//- Return the virtual mass coefficient for the phase pair
virtual tmp<volScalarField> Vm(const phasePairKey& key) const;
//- Add the mass-transfer-based momentum transfer to the equations
void addMassTransferMomentumTransfer
(
phaseSystem::momentumTransferTable& eqns
) const;
public:
@ -172,12 +205,10 @@ public:
//- Return the momentum transfer matrices for the cell-based algorithm.
// This includes implicit and explicit forces that add into the cell
// UEqn in the normal way.
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransfer() const;
virtual autoPtr<phaseSystem::momentumTransferTable> momentumTransfer();
//- As momentumTransfer, but for the face-based algorithm
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransferf() const;
virtual autoPtr<phaseSystem::momentumTransferTable> momentumTransferf();
//- Return implicit force coefficients on the faces, for the face-based
// algorithm.

153
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
View File

@ -29,34 +29,20 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class BasePhaseSystem>
void Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::
addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const
Foam::tmp<Foam::volScalarField>
Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::pDmdt
(
const phasePairKey& key
) const
{
// Source term due to mass transfer
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
if (!pDmdt_.found(key))
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
const volVectorField& U1(pair.phase1().U());
const volVectorField& U2(pair.phase2().U());
const volScalarField dmdt(this->pDmdt(pair));
const volScalarField dmdt21(posPart(dmdt));
const volScalarField dmdt12(negPart(dmdt));
*eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1);
*eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2);
return phaseSystem::dmdt(key);
}
const scalar pDmdtSign(Pair<word>::compare(pDmdt_.find(key).key(), key));
return pDmdtSign**pDmdt_[key];
}
@ -152,112 +138,31 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::pDmdt
Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdt
(
const phasePairKey& key
) const
{
if (pDmdt_.found(key))
{
const scalar pDmdtSign
(
Pair<word>::compare(pDmdt_.find(key).key(), key)
);
return BasePhaseSystem::dmdt(key) + this->pDmdt(key);
}
return pDmdtSign**pDmdt_[key];
template<class BasePhaseSystem>
Foam::Xfer<Foam::PtrList<Foam::volScalarField>>
Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdts() const
{
PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
forAllConstIter(pDmdtTable, pDmdt_, pDmdtIter)
{
const phasePair& pair = this->phasePairs_[pDmdtIter.key()];
const volScalarField& pDmdt = *pDmdtIter();
this->addField(pair.phase1(), "dmdt", pDmdt, dmdts);
this->addField(pair.phase2(), "dmdt", - pDmdt, dmdts);
}
tmp<volScalarField> tpDmdt
(
new volScalarField
(
IOobject
(
"dmdt",
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
return tpDmdt;
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdt
(
const phaseModel& phase
) const
{
tmp<volScalarField> tDmdtFromKey
(
new volScalarField
(
IOobject
(
IOobject::groupName("dmdtFromKey", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
if (pair.contains(phase))
{
tDmdtFromKey.ref() += this->dmdt
(
phasePairKey(phase.name(), pair.otherPhase(phase).name(), false)
);
}
}
return tDmdtFromKey;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransfer() const
{
autoPtr<phaseSystem::momentumTransferTable>
eqnsPtr(BasePhaseSystem::momentumTransfer());
addMomentumTransfer(eqnsPtr());
return eqnsPtr;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransferf() const
{
autoPtr<phaseSystem::momentumTransferTable>
eqnsPtr(BasePhaseSystem::momentumTransferf());
addMomentumTransfer(eqnsPtr());
return eqnsPtr;
return dmdts.xfer();
}

38
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H
View File

@ -58,25 +58,27 @@ class PopulationBalancePhaseSystem
:
public BasePhaseSystem
{
protected:
// Protected typedefs
typedef HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
pDmdtTable;
// Protected data
//- populationBalanceModels
PtrList<diameterModels::populationBalanceModel> populationBalances_;
//- Interfacial Mass transfer rate
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash> pDmdt_;
pDmdtTable pDmdt_;
// Protected member functions
//- Add the mass-based momentum transfer to the equations
void addMomentumTransfer
(
phaseSystem::momentumTransferTable& eqns
) const;
//- Return the population balance mass transfer rate
virtual tmp<volScalarField> pDmdt(const phasePairKey& key) const;
public:
@ -93,25 +95,11 @@ public:
// Member Functions
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> pDmdt(const phasePairKey& key) const;
//- Return the mass transfer rate for a pair
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const
{
return BasePhaseSystem::dmdt(key) + this->pDmdt(key);
};
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
//- Return the momentum transfer matrices
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransfer() const;
//- Return the momentum transfer matrices for the face based algorithm
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransferf() const;
//- Return the mass transfer rates for each phase
virtual Xfer<PtrList<volScalarField>> dmdts() const;
//- Return the heat transfer matrices
virtual autoPtr<phaseSystem::heatTransferTable> heatTransfer() const;

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

@ -31,34 +31,38 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class BasePhaseSystem>
void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const
Foam::tmp<Foam::volScalarField>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::iDmdt
(
const phasePairKey& key
) const
{
// Source term due to mass trasfer
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
if (!iDmdt_.found(key))
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
const volVectorField& U1(pair.phase1().U());
const volVectorField& U2(pair.phase2().U());
const volScalarField dmdt(this->iDmdt(pair) + this->wDmdt(pair));
const volScalarField dmdt21(posPart(dmdt));
const volScalarField dmdt12(negPart(dmdt));
*eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1);
*eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2);
return phaseSystem::dmdt(key);
}
const scalar dmdtSign(Pair<word>::compare(iDmdt_.find(key).key(), key));
return dmdtSign**iDmdt_[key];
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::wDmdt
(
const phasePairKey& key
) const
{
if (!wDmdt_.found(key))
{
return phaseSystem::dmdt(key);
}
const scalar dmdtSign(Pair<word>::compare(wDmdt_.find(key).key(), key));
return dmdtSign**wDmdt_[key];
}
@ -172,201 +176,42 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::saturation() const
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::iDmdt
(
const phasePairKey& key
) const
{
const scalar dmdtSign(Pair<word>::compare(iDmdt_.find(key).key(), key));
return dmdtSign**iDmdt_[key];
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::iDmdt
(
const Foam::phaseModel& phase
) const
{
tmp<volScalarField> tiDmdt
(
new volScalarField
(
IOobject
(
IOobject::groupName("iDmdt", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
if (pair.contains(phase))
{
tiDmdt.ref() += this->iDmdt
(
phasePairKey(phase.name(), pair.otherPhase(phase).name())
);
}
}
return tiDmdt;
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::wDmdt
(
const phasePairKey& key
) const
{
const scalar dmdtSign(Pair<word>::compare(wDmdt_.find(key).key(), key));
return dmdtSign**wDmdt_[key];
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::wDmdt
(
const Foam::phaseModel& phase
) const
{
tmp<volScalarField> twDmdt
(
new volScalarField
(
IOobject
(
IOobject::groupName("wDmdt", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
if (pair.contains(phase))
{
twDmdt.ref() += this->wDmdt
(
phasePairKey(phase.name(), pair.otherPhase(phase).name())
);
}
}
return twDmdt;
}
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::dmdt
(
const Foam::phaseModel& phase
const phasePairKey& key
) const
{
tmp<volScalarField> tDmdt
(
new volScalarField
(
IOobject
(
IOobject::groupName("Dmdt", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
return BasePhaseSystem::dmdt(key) + this->iDmdt(key) + this->wDmdt(key);
}
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
template<class BasePhaseSystem>
Foam::Xfer<Foam::PtrList<Foam::volScalarField>>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::dmdts() const
{
PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
forAllConstIter(iDmdtTable, iDmdt_, iDmdtIter)
{
const phasePair& pair(phasePairIter());
const phasePair& pair = this->phasePairs_[iDmdtIter.key()];
const volScalarField& iDmdt = *iDmdtIter();
if (pair.ordered())
{
continue;
}
if (pair.contains(phase))
{
tDmdt.ref() += this->dmdt
(
phasePairKey(phase.name(), pair.otherPhase(phase).name())
);
}
this->addField(pair.phase1(), "dmdt", iDmdt, dmdts);
this->addField(pair.phase2(), "dmdt", - iDmdt, dmdts);
}
return tDmdt;
}
forAllConstIter(wDmdtTable, wDmdt_, wDmdtIter)
{
const phasePair& pair = this->phasePairs_[wDmdtIter.key()];
const volScalarField& wDmdt = *wDmdtIter();
this->addField(pair.phase1(), "dmdt", wDmdt, dmdts);
this->addField(pair.phase2(), "dmdt", - wDmdt, dmdts);
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::momentumTransfer() const
{
autoPtr<phaseSystem::momentumTransferTable>
eqnsPtr(BasePhaseSystem::momentumTransfer());
addMomentumTransfer(eqnsPtr());
return eqnsPtr;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::momentumTransferf() const
{
autoPtr<phaseSystem::momentumTransferTable>
eqnsPtr(BasePhaseSystem::momentumTransferf());
addMomentumTransfer(eqnsPtr());
return eqnsPtr;
return dmdts.xfer();
}
@ -461,7 +306,6 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
forAll(Yi, i)
{
if (Yi[i].member() != volatile_)
{
continue;

63
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
View File

@ -58,6 +58,18 @@ class ThermalPhaseChangePhaseSystem
protected:
// Protected typedefs
typedef HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
iDmdtTable;
typedef HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
wDmdtTable;
typedef HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
wMDotLTable;
// Protected data
//- Name of the volatile specie
@ -70,25 +82,22 @@ protected:
Switch massTransfer_;
//- Interfacial Mass transfer rate
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
iDmdt_;
iDmdtTable iDmdt_;
//- Boundary Mass transfer rate
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
wDmdt_;
wDmdtTable wDmdt_;
//- Boundary thermal energy transfer rate
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
wMDotL_;
wMDotLTable wMDotL_;
// Protected member functions
//- Add the mass-based momentum transfer to the equations
void addMomentumTransfer
(
phaseSystem::momentumTransferTable& eqns
) const;
//- Return the interfacial mass transfer rate for a pair
tmp<volScalarField> iDmdt(const phasePairKey& key) const;
//- Return the boundary mass transfer rate for a pair
tmp<volScalarField> wDmdt(const phasePairKey& key) const;
public:
@ -108,35 +117,11 @@ public:
//- Return the saturationModel
const saturationModel& saturation() const;
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> iDmdt(const phasePairKey& key) const;
//- Return the mass transfer rate for a pair
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> iDmdt(const phaseModel& phase) const;
//- Return the boundary mass flow rate
virtual tmp<volScalarField> wDmdt(const phasePairKey& key) const;
//- Return the total boundary mass transfer rate for phase
virtual tmp<volScalarField> wDmdt(const phaseModel& phase) const;
//- Return the interfacial mass flow rate
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const
{
return BasePhaseSystem::dmdt(key)
+ this->iDmdt(key) + this->wDmdt(key);
};
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
//- Return the momentum transfer matrices
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransfer() const;
//- Return the momentum transfer matrices for the face based algorithm
virtual autoPtr<phaseSystem::momentumTransferTable>
momentumTransferf() const;
//- Return the mass transfer rates for each phase
virtual Xfer<PtrList<volScalarField>> dmdts() const;
//- Return the heat transfer matrices
virtual autoPtr<phaseSystem::heatTransferTable> heatTransfer() const;

121
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
View File

@ -26,63 +26,7 @@ License
#include "AnisothermalPhaseModel.H"
#include "phaseSystem.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasePhaseModel>
Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
(
const phaseSystem& fluid,
const word& phaseName,
const label index
)
:
BasePhaseModel(fluid, phaseName, index),
K_
(
IOobject
(
IOobject::groupName("K", this->name()),
fluid.mesh().time().timeName(),
fluid.mesh()
),
fluid.mesh(),
dimensionedScalar("K", sqr(dimVelocity), scalar(0))
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class BasePhaseModel>
Foam::AnisothermalPhaseModel<BasePhaseModel>::~AnisothermalPhaseModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasePhaseModel>
bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
{
return !this->thermo().incompressible();
}
template<class BasePhaseModel>
void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctKinematics()
{
BasePhaseModel::correctKinematics();
K_ = 0.5*magSqr(this->U());
}
template<class BasePhaseModel>
void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctThermo()
{
BasePhaseModel::correctThermo();
this->thermo_->correct();
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
@ -111,18 +55,57 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::filterPressureWork
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasePhaseModel>
Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
(
const phaseSystem& fluid,
const word& phaseName,
const label index
)
:
BasePhaseModel(fluid, phaseName, index)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class BasePhaseModel>
Foam::AnisothermalPhaseModel<BasePhaseModel>::~AnisothermalPhaseModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasePhaseModel>
void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctThermo()
{
BasePhaseModel::correctThermo();
this->thermo_->correct();
}
template<class BasePhaseModel>
bool Foam::AnisothermalPhaseModel<BasePhaseModel>::isothermal() const
{
return false;
}
template<class BasePhaseModel>
Foam::tmp<Foam::fvScalarMatrix>
Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
{
const volScalarField& alpha = *this;
const volVectorField& U = this->U();
const surfaceScalarField& alphaPhi = this->alphaPhi();
const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi();
const volScalarField& contErr(this->continuityError());
const volVectorField U(this->U());
const surfaceScalarField alphaPhi(this->alphaPhi());
const surfaceScalarField alphaRhoPhi(this->alphaRhoPhi());
const volScalarField alphaEff(this->alphaEff());
const volScalarField contErr(this->continuityError());
const volScalarField K(this->K());
volScalarField& he = this->thermo_->he();
@ -132,13 +115,13 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
+ fvm::div(alphaRhoPhi, he)
- fvm::Sp(contErr, he)
+ fvc::ddt(alpha, this->rho(), K_) + fvc::div(alphaRhoPhi, K_)
- contErr*K_
+ fvc::ddt(alpha, this->rho(), K) + fvc::div(alphaRhoPhi, K)
- contErr*K
- fvm::laplacian
(
fvc::interpolate(alpha)
*fvc::interpolate(alphaEff),
*fvc::interpolate(this->alphaEff()),
he
)
==
@ -163,12 +146,4 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
}
template<class BasePhaseModel>
const Foam::volScalarField&
Foam::AnisothermalPhaseModel<BasePhaseModel>::K() const
{
return K_;
}
// ************************************************************************* //

22
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -44,7 +44,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phaseModel Declaration
Class AnisothermalPhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel>
@ -52,12 +52,6 @@ class AnisothermalPhaseModel
:
public BasePhaseModel
{
// Private data
//- Kinetic energy
volScalarField K_;
// Private member functions
//- Optionally filter the pressure work term as the phase-fraction -> 0
@ -85,20 +79,14 @@ public:
// Member Functions
//- Return true if the phase is compressible otherwise false
virtual bool compressible() const;
//- Correct the kinematics
virtual void correctKinematics();
//- Correct the thermodynamics
virtual void correctThermo();
//- Return whether the phase is isothermal
virtual bool isothermal() const;
//- Return the enthalpy equation
virtual tmp<fvScalarMatrix> heEqn();
//- Return the phase kinetic energy
virtual const volScalarField& K() const;
};

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -44,7 +44,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phaseModel Declaration
Class InertPhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel>
@ -73,10 +73,7 @@ public:
// Thermo
//- Return the fuel consumption rate matrix
virtual tmp<fvScalarMatrix> R
(
volScalarField& Yi
) const;
virtual tmp<fvScalarMatrix> R(volScalarField& Yi) const;
//- Return the heat release rate
virtual tmp<volScalarField> Qdot() const;

18
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -50,21 +50,25 @@ Foam::IsothermalPhaseModel<BasePhaseModel>::~IsothermalPhaseModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasePhaseModel>
bool Foam::IsothermalPhaseModel<BasePhaseModel>::compressible() const
{
return !this->thermo().incompressible();
}
void Foam::IsothermalPhaseModel<BasePhaseModel>::correctThermo()
{}
template<class BasePhaseModel>
void Foam::IsothermalPhaseModel<BasePhaseModel>::correctThermo()
{}
bool Foam::IsothermalPhaseModel<BasePhaseModel>::isothermal() const
{
return true;
}
template<class BasePhaseModel>
Foam::tmp<Foam::fvScalarMatrix>
Foam::IsothermalPhaseModel<BasePhaseModel>::heEqn()
{
FatalErrorInFunction
<< "Cannot construct an energy equation for an isothermal phase"
<< exit(FatalError);
return tmp<fvScalarMatrix>();
}

10
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -45,7 +45,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phaseModel Declaration
Class IsothermalPhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel>
@ -71,12 +71,12 @@ public:
// Member Functions
//- Return true if the phase is compressible otherwise false
virtual bool compressible() const;
//- Correct the thermodynamics
virtual void correctThermo();
//- Return whether the phase is isothermal
virtual bool isothermal() const;
//- Return the enthalpy equation
virtual tmp<fvScalarMatrix> heEqn();
};

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

@ -179,19 +179,32 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
*this
)
),
continuityError_
continuityErrorFlow_
(
IOobject
(
IOobject::groupName("continuityError", this->name()),
IOobject::groupName("continuityErrorFlow", this->name()),
fluid.mesh().time().timeName(),
fluid.mesh()
),
fluid.mesh(),
dimensionedScalar("0", dimDensity/dimTime, 0)
)
),
continuityErrorSources_
(
IOobject
(
IOobject::groupName("continuityErrorSources", this->name()),
fluid.mesh().time().timeName(),
fluid.mesh()
),
fluid.mesh(),
dimensionedScalar("0", dimDensity/dimTime, 0)
),
K_(nullptr)
{
phi_.writeOpt() = IOobject::AUTO_WRITE;
correctKinematics();
}
@ -212,11 +225,11 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correct()
this->fluid().MRF().correctBoundaryVelocity(U_);
volScalarField& rho = this->thermo().rho();
volScalarField& rho = this->thermoRef().rho();
continuityError_ =
fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_)
- (this->fluid().fvOptions()(*this, rho)&rho);
continuityErrorFlow_ = fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_);
continuityErrorSources_ = - (this->fluid().fvOptions()(*this, rho)&rho);
}
@ -236,6 +249,12 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correctKinematics()
DUDtf_.clear();
DUDtf();
}
if (K_.valid())
{
K_.clear();
K();
}
}
@ -257,19 +276,25 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correctEnergyTransport()
}
template<class BasePhaseModel>
bool Foam::MovingPhaseModel<BasePhaseModel>::stationary() const
{
return false;
}
template<class BasePhaseModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::MovingPhaseModel<BasePhaseModel>::UEqn()
{
const volScalarField& alpha = *this;
volScalarField& rho = this->thermo().rho();
const volScalarField& rho = this->thermo().rho();
return
(
fvm::ddt(alpha, rho, U_)
+ fvm::div(alphaRhoPhi_, U_)
- fvm::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi_), U_)
+ fvm::SuSp(this->fluid().fvOptions()(alpha, rho)&rho, U_)
+ fvm::SuSp(- this->continuityError(), U_)
+ this->fluid().MRF().DDt(alpha*rho, U_)
+ turbulence_->divDevRhoReff(U_)
);
@ -283,13 +308,13 @@ Foam::MovingPhaseModel<BasePhaseModel>::UfEqn()
// As the "normal" U-eqn but without the ddt terms
const volScalarField& alpha = *this;
volScalarField& rho = this->thermo().rho();
const volScalarField& rho = this->thermo().rho();
return
(
fvm::div(alphaRhoPhi_, U_)
- fvm::Sp(fvc::div(alphaRhoPhi_), U_)
+ fvm::SuSp(this->fluid().fvOptions()(alpha, rho)&rho, U_)
+ fvm::SuSp(- this->continuityErrorSources(), U_)
+ this->fluid().MRF().DDt(alpha*rho, U_)
+ turbulence_->divDevRhoReff(U_)
);
@ -306,12 +331,60 @@ Foam::MovingPhaseModel<BasePhaseModel>::U() const
template<class BasePhaseModel>
Foam::volVectorField&
Foam::MovingPhaseModel<BasePhaseModel>::U()
Foam::MovingPhaseModel<BasePhaseModel>::URef()
{
return U_;
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::phi() const
{
return phi_;
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::MovingPhaseModel<BasePhaseModel>::phiRef()
{
return phi_;
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::alphaPhi() const
{
return alphaPhi_;
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::MovingPhaseModel<BasePhaseModel>::alphaPhiRef()
{
return alphaPhi_;
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::alphaRhoPhi() const
{
return alphaRhoPhi_;
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::MovingPhaseModel<BasePhaseModel>::alphaRhoPhiRef()
{
return alphaRhoPhi_;
}
template<class BasePhaseModel>
Foam::tmp<Foam::volVectorField>
Foam::MovingPhaseModel<BasePhaseModel>::DUDt() const
@ -339,76 +412,59 @@ Foam::MovingPhaseModel<BasePhaseModel>::DUDtf() const
template<class BasePhaseModel>
const Foam::tmp<Foam::volScalarField>&
Foam::MovingPhaseModel<BasePhaseModel>::divU() const
Foam::tmp<Foam::volScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::continuityError() const
{
return divU_;
}
template<class BasePhaseModel>
void Foam::MovingPhaseModel<BasePhaseModel>::divU
(
const tmp<volScalarField>& divU
)
{
divU_ = divU;
return continuityErrorFlow_ + continuityErrorSources_;
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::continuityError() const
Foam::MovingPhaseModel<BasePhaseModel>::continuityErrorFlow() const
{
return continuityError_;
return continuityErrorFlow_;
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::phi() const
Foam::tmp<Foam::volScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::continuityErrorSources() const
{
return phi_;
return continuityErrorSources_;
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::MovingPhaseModel<BasePhaseModel>::phi()
Foam::tmp<Foam::volScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::K() const
{
return phi_;
if (!K_.valid())
{
K_ =
new volScalarField
(
IOobject::groupName("K", this->name()),
0.5*magSqr(this->U())
);
}
return tmp<volScalarField>(K_());
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::alphaPhi() const
Foam::tmp<Foam::volScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::divU() const
{
return alphaPhi_;
return divU_.valid() ? tmp<volScalarField>(divU_()) : tmp<volScalarField>();
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::MovingPhaseModel<BasePhaseModel>::alphaPhi()
void Foam::MovingPhaseModel<BasePhaseModel>::divU(tmp<volScalarField> divU)
{
return alphaPhi_;
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::MovingPhaseModel<BasePhaseModel>::alphaRhoPhi() const
{
return alphaRhoPhi_;
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::MovingPhaseModel<BasePhaseModel>::alphaRhoPhi()
{
return alphaRhoPhi_;
divU_ = divU;
}

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

@ -55,7 +55,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phaseModel Declaration
Class MovingPhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel>
@ -91,8 +91,14 @@ protected:
//- Turbulence model
autoPtr<phaseCompressibleTurbulenceModel> turbulence_;
//- Continuity error
volScalarField continuityError_;
//- Continuity error due to the flow
volScalarField continuityErrorFlow_;
//- Continuity error due to any sources
volScalarField continuityErrorSources_;
//- Kinetic Energy
mutable tmp<volScalarField> K_;
private:
@ -133,20 +139,41 @@ public:
//- Correct the energy transport e.g. alphat
virtual void correctEnergyTransport();
//- Return the momentum equation
virtual tmp<fvVectorMatrix> UEqn();
//- Return the momentum equation for the face-based algorithm
virtual tmp<fvVectorMatrix> UfEqn();
// Momentum
//- Constant access the velocity
//- Return whether the phase is stationary
virtual bool stationary() const;
//- Return the momentum equation
virtual tmp<fvVectorMatrix> UEqn();
//- Return the momentum equation for the face-based algorithm
virtual tmp<fvVectorMatrix> UfEqn();
//- Return the velocity
virtual tmp<volVectorField> U() const;
//- Access the velocity
virtual volVectorField& U();
virtual volVectorField& URef();
//- Return the volumetric flux
virtual tmp<surfaceScalarField> phi() const;
//- Access the volumetric flux
virtual surfaceScalarField& phiRef();
//- Return the volumetric flux of the phase
virtual tmp<surfaceScalarField> alphaPhi() const;
//- Access the volumetric flux of the phase
virtual surfaceScalarField& alphaPhiRef();
//- Return the mass flux of the phase
virtual tmp<surfaceScalarField> alphaRhoPhi() const;
//- Access the mass flux of the phase
virtual surfaceScalarField& alphaRhoPhiRef();
//- Return the substantive acceleration
virtual tmp<volVectorField> DUDt() const;
@ -154,35 +181,38 @@ public:
//- Return the substantive acceleration on the faces
virtual tmp<surfaceScalarField> DUDtf() const;
//- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual const tmp<volScalarField>& divU() const;
//- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual void divU(const tmp<volScalarField>& divU);
//- Constant access the continuity error
//- Return the continuity error
virtual tmp<volScalarField> continuityError() const;
//- Constant access the volumetric flux
virtual tmp<surfaceScalarField> phi() const;
//- Return the continuity error due to the flow field
virtual tmp<volScalarField> continuityErrorFlow() const;
//- Access the volumetric flux
virtual surfaceScalarField& phi();
//- Return the continuity error due to any sources
virtual tmp<volScalarField> continuityErrorSources() const;
//- Constant access the volumetric flux of the phase
virtual tmp<surfaceScalarField> alphaPhi() const;
//- Access the volumetric flux of the phase
virtual surfaceScalarField& alphaPhi();
//- Constant access the mass flux of the phase
virtual tmp<surfaceScalarField> alphaRhoPhi() const;
//- Access the mass flux of the phase
virtual surfaceScalarField& alphaRhoPhi();
//- Return the phase kinetic energy
virtual tmp<volScalarField> K() const;
// Transport
// Compressibility (variable density)
//- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual tmp<volScalarField> divU() const;
//- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual void divU(tmp<volScalarField> divU);
// Transport (prevents compiler warnings)
//- Return the effective thermal conductivity
using BasePhaseModel::kappaEff;
//- Return the effective thermal conductivity for enthalpy
using BasePhaseModel::alphaEff;
// Turbulence
//- Return the turbulent dynamic viscosity
virtual tmp<volScalarField> mut() const;

72
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
View File

@ -69,6 +69,18 @@ Foam::MultiComponentPhaseModel<BasePhaseModel>::MultiComponentPhaseModel
{
inertIndex_ = this->thermo_->composition().species()[inertSpecie];
}
PtrList<volScalarField>& Y = this->thermo_->composition().Y();
forAll(Y, i)
{
if (i != inertIndex_ && this->thermo_->composition().active(i))
{
const label j = YActive_.size();
YActive_.resize(j + 1);
YActive_.set(j, &Y[i]);
}
}
}
@ -96,7 +108,7 @@ void Foam::MultiComponentPhaseModel<BasePhaseModel>::correctThermo()
dimensionedScalar("zero", dimless, 0)
);
PtrList<volScalarField>& Yi = Y();
PtrList<volScalarField>& Yi = YRef();
forAll(Yi, i)
{
@ -125,39 +137,19 @@ void Foam::MultiComponentPhaseModel<BasePhaseModel>::correctThermo()
template<class BasePhaseModel>
Foam::tmp<Foam::fvScalarMatrix>
Foam::MultiComponentPhaseModel<BasePhaseModel>::YiEqn
(
volScalarField& Yi
)
bool Foam::MultiComponentPhaseModel<BasePhaseModel>::pure() const
{
if
(
(inertIndex_ != -1)
&& (
(
Yi.name()
== IOobject::groupName
(
this->thermo_->composition().species()[inertIndex_],
this->name()
)
)
|| (
!this->thermo_->composition().active
(
this->thermo_->composition().species()[Yi.member()]
)
)
)
)
{
return tmp<fvScalarMatrix>();
}
return false;
}
template<class BasePhaseModel>
Foam::tmp<Foam::fvScalarMatrix>
Foam::MultiComponentPhaseModel<BasePhaseModel>::YiEqn(volScalarField& Yi)
{
const volScalarField& alpha = *this;
const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi();
const volScalarField& rho = this->rho();
const surfaceScalarField alphaRhoPhi(this->alphaRhoPhi());
const volScalarField rho(this->rho());
return
(
@ -190,10 +182,26 @@ Foam::MultiComponentPhaseModel<BasePhaseModel>::Y() const
template<class BasePhaseModel>
Foam::PtrList<Foam::volScalarField>&
Foam::MultiComponentPhaseModel<BasePhaseModel>::Y()
Foam::MultiComponentPhaseModel<BasePhaseModel>::YRef()
{
return this->thermo_->composition().Y();
}
template<class BasePhaseModel>
const Foam::UPtrList<Foam::volScalarField>&
Foam::MultiComponentPhaseModel<BasePhaseModel>::YActive() const
{
return YActive_;
}
template<class BasePhaseModel>
Foam::UPtrList<Foam::volScalarField>&
Foam::MultiComponentPhaseModel<BasePhaseModel>::YActiveRef()
{
return YActive_;
}
// ************************************************************************* //

27
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -44,7 +44,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phaseModel Declaration
Class MultiComponentPhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel>
@ -65,6 +65,9 @@ protected:
//- Inert species index
label inertIndex_;
//- Pointer list to active species
UPtrList<volScalarField> YActive_;
public:
@ -87,17 +90,23 @@ public:
//- Correct the thermodynamics
virtual void correctThermo();
//- Return the species fraction equation
virtual tmp<fvScalarMatrix> YiEqn
(
volScalarField& Yi
);
//- Return whether the phase is pure (i.e., not multi-component)
virtual bool pure() const;
//- Constant access the species mass fractions
//- Return the species fraction equation
virtual tmp<fvScalarMatrix> YiEqn(volScalarField& Yi);
//- Return the species mass fractions
virtual const PtrList<volScalarField>& Y() const;
//- Access the species mass fractions
virtual PtrList<volScalarField>& Y();
virtual PtrList<volScalarField>& YRef();
//- Return the active species mass fractions
virtual const UPtrList<volScalarField>& YActive() const;
//- Access the active species mass fractions
virtual UPtrList<volScalarField>& YActiveRef();
};

50
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -50,13 +50,19 @@ Foam::PurePhaseModel<BasePhaseModel>::~PurePhaseModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasePhaseModel>
Foam::tmp<Foam::fvScalarMatrix>
Foam::PurePhaseModel<BasePhaseModel>::YiEqn
(
volScalarField& Yi
)
bool Foam::PurePhaseModel<BasePhaseModel>::pure() const
{
NotImplemented;
return true;
}
template<class BasePhaseModel>
Foam::tmp<Foam::fvScalarMatrix>
Foam::PurePhaseModel<BasePhaseModel>::YiEqn(volScalarField& Yi)
{
FatalErrorInFunction
<< "Cannot construct a species fraction equation for a pure phase"
<< exit(FatalError);
return tmp<fvScalarMatrix>();
}
@ -66,14 +72,42 @@ template<class BasePhaseModel>
const Foam::PtrList<Foam::volScalarField>&
Foam::PurePhaseModel<BasePhaseModel>::Y() const
{
// Y_ has never been set, so we are returning an empty list
return Y_;
}
template<class BasePhaseModel>
Foam::PtrList<Foam::volScalarField>&
Foam::PurePhaseModel<BasePhaseModel>::Y()
Foam::PurePhaseModel<BasePhaseModel>::YRef()
{
FatalErrorInFunction
<< "Cannot access the species fractions of for a pure phase"
<< exit(FatalError);
return Y_;
}
template<class BasePhaseModel>
const Foam::UPtrList<Foam::volScalarField>&
Foam::PurePhaseModel<BasePhaseModel>::YActive() const
{
// Y_ has never been set, so we are returning an empty list
return Y_;
}
template<class BasePhaseModel>
Foam::UPtrList<Foam::volScalarField>&
Foam::PurePhaseModel<BasePhaseModel>::YActiveRef()
{
FatalErrorInFunction
<< "Cannot access the species fractions of for a pure phase"
<< exit(FatalError);
return Y_;
}

20
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -44,7 +44,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phaseModel Declaration
Class PurePhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel>
@ -78,17 +78,25 @@ public:
// Member Functions
//- Return the species fraction equation
virtual tmp<fvScalarMatrix> YiEqn(volScalarField& Yi);
// Species
//- Return whether the phase is pure (i.e., not multi-component)
virtual bool pure() const;
// Thermo
//- Return the species fraction equation
virtual tmp<fvScalarMatrix> YiEqn(volScalarField& Yi);
//- Return the species mass fractions
virtual const PtrList<volScalarField>& Y() const;
//- Access the species mass fractions
virtual PtrList<volScalarField>& Y();
virtual PtrList<volScalarField>& YRef();
//- Return the active species mass fractions
virtual const UPtrList<volScalarField>& YActive() const;
//- Access the active species mass fractions
virtual UPtrList<volScalarField>& YActiveRef();
};

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -44,7 +44,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phaseModel Declaration
Class ReactingPhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel, class ReactionType>
@ -82,10 +82,7 @@ public:
virtual void correctThermo();
//- Return the species fraction equation
virtual tmp<fvScalarMatrix> R
(
volScalarField& Yi
) const;
virtual tmp<fvScalarMatrix> R(volScalarField& Yi) const;
//- Return heat release rate
virtual tmp<volScalarField> Qdot() const;

367
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C Normal file
View File

@ -0,0 +1,367 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2018 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 "StationaryPhaseModel.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class BasePhaseModel>
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
Foam::StationaryPhaseModel<BasePhaseModel>::zeroField
(
const word& name,
const dimensionSet& dims,
const bool cache
) const
{
return tmp<GeometricField<Type, PatchField, GeoMesh>>
(
new GeometricField<Type, PatchField, GeoMesh>
(
IOobject
(
IOobject::groupName(name, this->name()),
this->mesh().time().timeName(),
this->mesh()
),
this->mesh(),
dimensioned<Type>("zero", dims, pTraits<Type>::zero)
)
);
}
template<class BasePhaseModel>
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
Foam::StationaryPhaseModel<BasePhaseModel>::zeroVolField
(
const word& name,
const dimensionSet& dims,
const bool cache
) const
{
return zeroField<Type, fvPatchField, volMesh>(name, dims, cache);
}
template<class BasePhaseModel>
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
Foam::StationaryPhaseModel<BasePhaseModel>::zeroSurfaceField
(
const word& name,
const dimensionSet& dims,
const bool cache
) const
{
return zeroField<Type, fvsPatchField, surfaceMesh>(name, dims, cache);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasePhaseModel>
Foam::StationaryPhaseModel<BasePhaseModel>::StationaryPhaseModel
(
const phaseSystem& fluid,
const word& phaseName,
const label index
)
:
BasePhaseModel(fluid, phaseName, index)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class BasePhaseModel>
Foam::StationaryPhaseModel<BasePhaseModel>::~StationaryPhaseModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasePhaseModel>
bool Foam::StationaryPhaseModel<BasePhaseModel>::stationary() const
{
return true;
}
template<class BasePhaseModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::StationaryPhaseModel<BasePhaseModel>::UEqn()
{
FatalErrorInFunction
<< "Cannot construct a momentum equation for a stationary phase"
<< exit(FatalError);
return tmp<fvVectorMatrix>();
}
template<class BasePhaseModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::StationaryPhaseModel<BasePhaseModel>::UfEqn()
{
FatalErrorInFunction
<< "Cannot construct a momentum equation for a stationary phase"
<< exit(FatalError);
return tmp<fvVectorMatrix>();
}
template<class BasePhaseModel>
Foam::tmp<Foam::volVectorField>
Foam::StationaryPhaseModel<BasePhaseModel>::U() const
{
return zeroVolField<vector>("U", dimVelocity, true);
}
template<class BasePhaseModel>
Foam::volVectorField&
Foam::StationaryPhaseModel<BasePhaseModel>::URef()
{
FatalErrorInFunction
<< "Cannot access the velocity of a stationary phase"
<< exit(FatalError);
return const_cast<volVectorField&>(volVectorField::null());
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::phi() const
{
return zeroSurfaceField<scalar>("phi", dimVolume/dimTime);
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::StationaryPhaseModel<BasePhaseModel>::phiRef()
{
FatalErrorInFunction
<< "Cannot access the flux of a stationary phase"
<< exit(FatalError);
return const_cast<surfaceScalarField&>(surfaceScalarField::null());
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::alphaPhi() const
{
return zeroSurfaceField<scalar>("alphaPhi", dimVolume/dimTime);
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::StationaryPhaseModel<BasePhaseModel>::alphaPhiRef()
{
FatalErrorInFunction
<< "Cannot access the volumetric flux of a stationary phase"
<< exit(FatalError);
return const_cast<surfaceScalarField&>(surfaceScalarField::null());
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::alphaRhoPhi() const
{
return zeroSurfaceField<scalar>("alphaRhoPhi", dimMass/dimTime);
}
template<class BasePhaseModel>
Foam::surfaceScalarField&
Foam::StationaryPhaseModel<BasePhaseModel>::alphaRhoPhiRef()
{
FatalErrorInFunction
<< "Cannot access the mass flux of a stationary phase"
<< exit(FatalError);
return const_cast<surfaceScalarField&>(surfaceScalarField::null());
}
template<class BasePhaseModel>
Foam::tmp<Foam::volVectorField>
Foam::StationaryPhaseModel<BasePhaseModel>::DUDt() const
{
return zeroVolField<vector>("DUDt", dimVelocity/dimTime);
}
template<class BasePhaseModel>
Foam::tmp<Foam::surfaceScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::DUDtf() const
{
return zeroSurfaceField<scalar>("DUDtf", dimVelocity*dimArea/dimTime);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::continuityError() const
{
return zeroVolField<scalar>("continuityError", dimDensity/dimTime);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::continuityErrorFlow() const
{
return zeroVolField<scalar>("continuityErrorFlow", dimDensity/dimTime);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::continuityErrorSources() const
{
return zeroVolField<scalar>("continuityErrorSources", dimDensity/dimTime);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::K() const
{
return zeroVolField<scalar>("K", sqr(dimVelocity));
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::divU() const
{
return tmp<volScalarField>();
}
template<class BasePhaseModel>
void Foam::StationaryPhaseModel<BasePhaseModel>::divU
(
tmp<volScalarField> divU
)
{
FatalErrorInFunction
<< "Cannot set the dilatation rate of a stationary phase"
<< exit(FatalError);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::mut() const
{
return zeroVolField<scalar>("continuityError", dimDynamicViscosity);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::muEff() const
{
return this->thermo().mu();
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::nut() const
{
return zeroVolField<scalar>("continuityError", dimViscosity);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::nuEff() const
{
return this->thermo().nu();
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::kappaEff() const
{
return this->thermo().kappa();
}
template<class BasePhaseModel>
Foam::tmp<Foam::scalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::kappaEff(const label patchi) const
{
return this->thermo().kappa(patchi);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::alphaEff() const
{
return this->thermo().alpha();
}
template<class BasePhaseModel>
Foam::tmp<Foam::scalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::alphaEff(const label patchi) const
{
return this->thermo().alpha(patchi);
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::k() const
{
return zeroVolField<scalar>("k", sqr(dimVelocity));
}
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
Foam::StationaryPhaseModel<BasePhaseModel>::pPrime() const
{
return zeroVolField<scalar>("pPrime", dimPressure);
}
// ************************************************************************* //

226
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H Normal file
View File

@ -0,0 +1,226 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2018 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/>.
Class
Foam::StationaryPhaseModel
Description
...
SourceFiles
StationaryPhaseModel.C
\*---------------------------------------------------------------------------*/
#ifndef StationaryPhaseModel_H
#define StationaryPhaseModel_H
#include "phaseModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class StationaryPhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel>
class StationaryPhaseModel
:
public BasePhaseModel
{
private:
// Private member functions
//- Create a zero geometric field
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<Type, PatchField, GeoMesh>> zeroField
(
const word& name,
const dimensionSet& dims,
const bool cache = false
) const;
//- Create a zero vol field
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh>> zeroVolField
(
const word& name,
const dimensionSet& dims,
const bool cache = false
) const;
//- Create a zero surface field
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> zeroSurfaceField
(
const word& name,
const dimensionSet& dims,
const bool cache = false
) const;
public:
// Constructors
StationaryPhaseModel
(
const phaseSystem& fluid,
const word& phaseName,
const label index
);
//- Destructor
virtual ~StationaryPhaseModel();
// Member Functions
// Momentum
//- Return whether the phase is stationary
virtual bool stationary() const;
//- Return the momentum equation
virtual tmp<fvVectorMatrix> UEqn();
//- Return the momentum equation for the face-based algorithm
virtual tmp<fvVectorMatrix> UfEqn();
//- Return the velocity
virtual tmp<volVectorField> U() const;
//- Access the velocity
virtual volVectorField& URef();
//- Return the volumetric flux
virtual tmp<surfaceScalarField> phi() const;
//- Access the volumetric flux
virtual surfaceScalarField& phiRef();
//- Return the volumetric flux of the phase
virtual tmp<surfaceScalarField> alphaPhi() const;
//- Access the volumetric flux of the phase
virtual surfaceScalarField& alphaPhiRef();
//- Return the mass flux of the phase
virtual tmp<surfaceScalarField> alphaRhoPhi() const;
//- Access the mass flux of the phase
virtual surfaceScalarField& alphaRhoPhiRef();
//- Return the substantive acceleration
virtual tmp<volVectorField> DUDt() const;
//- Return the substantive acceleration on the faces
virtual tmp<surfaceScalarField> DUDtf() const;
//- Return the continuity error
virtual tmp<volScalarField> continuityError() const;
//- Return the continuity error due to the flow field
virtual tmp<volScalarField> continuityErrorFlow() const;
//- Return the continuity error due to any sources
virtual tmp<volScalarField> continuityErrorSources() const;
//- Return the phase kinetic energy
virtual tmp<volScalarField> K() const;
// Compressibility (variable density)
//- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual tmp<volScalarField> divU() const;
//- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual void divU(tmp<volScalarField> divU);
// Transport (prevents compiler warnings)
//- Return the effective thermal conductivity
using BasePhaseModel::kappaEff;
//- Return the effective thermal conductivity for enthalpy
using BasePhaseModel::alphaEff;
// Turbulence
//- Return the turbulent dynamic viscosity
virtual tmp<volScalarField> mut() const;
//- Return the effective dynamic viscosity
virtual tmp<volScalarField> muEff() const;
//- Return the turbulent kinematic viscosity
virtual tmp<volScalarField> nut() const;
//- Return the effective kinematic viscosity
virtual tmp<volScalarField> nuEff() const;
//- Return the effective thermal conductivity
virtual tmp<volScalarField> kappaEff() const;
//- Return the effective thermal conductivity on a patch
virtual tmp<scalarField> kappaEff(const label patchi) const;
//- Return the effective thermal diffusivity for enthalpy
virtual tmp<volScalarField> alphaEff() const;
//- Return the effective thermal conductivity for enthalpy on a
// patch
virtual tmp<scalarField> alphaEff(const label patchi) const;
//- Return the turbulent kinetic energy
virtual tmp<volScalarField> k() const;
//- Return the phase-pressure'
// (derivative of phase-pressure w.r.t. phase-fraction)
virtual tmp<volScalarField> pPrime() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "StationaryPhaseModel.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

11
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -65,6 +65,13 @@ Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::~ThermoPhaseModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasePhaseModel, class ThermoType>
bool Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::compressible() const
{
return !thermo_().incompressible();
}
template<class BasePhaseModel, class ThermoType>
const Foam::rhoThermo&
Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::thermo() const
@ -75,7 +82,7 @@ Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::thermo() const
template<class BasePhaseModel, class ThermoType>
Foam::rhoThermo&
Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::thermo()
Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::thermoRef()
{
return thermo_();
}

39
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H
View File

@ -49,7 +49,7 @@ namespace Foam
class rhoThermo;
/*---------------------------------------------------------------------------*\
Class phaseModel Declaration
Class ThermoPhaseModel Declaration
\*---------------------------------------------------------------------------*/
template<class BasePhaseModel, class ThermoType>
@ -85,12 +85,14 @@ public:
// Thermo
//- Return const access to the thermophysical model
//- Return whether the phase is compressible
virtual bool compressible() const;
//- Return the thermophysical model
virtual const rhoThermo& thermo() const;
//- Return non-const access to the thermophysical model
// for correction
virtual rhoThermo& thermo();
//- Access the thermophysical model
virtual rhoThermo& thermoRef();
//- Return the density field
virtual tmp<volScalarField> rho() const;
@ -101,28 +103,30 @@ public:
//- Return the laminar dynamic viscosity
virtual tmp<volScalarField> mu() const;
//- Access the laminar dynamic viscosity
//- Return the laminar dynamic viscosity on a patch
virtual tmp<scalarField> mu(const label patchi) const;
//- Return the laminar kinematic viscosity
virtual tmp<volScalarField> nu() const;
//- Access the laminar kinematic viscosity
//- Return the laminar kinematic viscosity on a patch
virtual tmp<scalarField> nu(const label patchi) const;
//- Return the laminar thermal conductivity
virtual tmp<volScalarField> kappa() const;
//- Access the laminar thermal conductivity
//- Return the laminar thermal conductivity on a patch
virtual tmp<scalarField> kappa(const label patchi) const;
//- Return the laminar thermal conductivity
//- Return the laminar thermal conductivity, given the turbulent
// thermal diffusivity
virtual tmp<volScalarField> kappaEff
(
const volScalarField& alphat
) const;
//- Access the laminar thermal conductivity
//- Return the laminar thermal conductivity on a patch, given the
// turbulent thermal diffusivity
virtual tmp<scalarField> kappaEff
(
const scalarField& alphat,
@ -135,18 +139,29 @@ public:
//- Return the thermal diffusivity for enthalpy on a patch
virtual tmp<scalarField> alpha(const label patchi) const;
//- Return the thermal diffusivity for enthalpy
//- Return the effective thermal diffusivity for enthalpy, given the
// turbulent thermal diffusivity
virtual tmp<volScalarField> alphaEff
(
const volScalarField& alphat
) const;
//- Return the thermal diffusivity for enthalpy on a patch
//- Return the effective thermal diffusivity for enthalpy on a
// patch, given the turbulent thermal diffusivity
virtual tmp<scalarField> alphaEff
(
const scalarField& alphat,
const label patchi
) const;
// Turbulence (prevents compiler warnings)
//- Return the effective thermal conductivity
using BasePhaseModel::kappaEff;
//- Return the effective thermal conductivity for enthalpy
using BasePhaseModel::alphaEff;
};

29
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
View File

@ -165,12 +165,6 @@ bool Foam::phaseModel::read()
}
bool Foam::phaseModel::compressible() const
{
return false;
}
void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
{
surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
@ -189,27 +183,4 @@ void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
}
const Foam::tmp<Foam::volScalarField>& Foam::phaseModel::divU() const
{
NotImplemented;
static tmp<Foam::volScalarField> divU_(nullptr);
return divU_;
}
void Foam::phaseModel::divU(const tmp<volScalarField>& divU)
{
WarningInFunction
<< "Attempt to set the dilatation rate of an incompressible phase"
<< endl;
}
const Foam::volScalarField& Foam::phaseModel::K() const
{
NotImplemented;
return volScalarField::null();
}
// ************************************************************************* //

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

@ -195,20 +195,11 @@ public:
//- Correct the turbulence
virtual void correctTurbulence();
//- Correct the energy transport e.g. alphat
//- Correct the energy transport
virtual void correctEnergyTransport();
//- Return the momentum equation
virtual tmp<fvVectorMatrix> UEqn() = 0;
//- Return the momentum equation for the face-based algorithm
virtual tmp<fvVectorMatrix> UfEqn() = 0;
//- Return the enthalpy equation
virtual tmp<fvScalarMatrix> heEqn() = 0;
//- Return the species fraction equation
virtual tmp<fvScalarMatrix> YiEqn(volScalarField& Yi) = 0;
//- Ensure that the flux at inflow/outflow BCs is preserved
void correctInflowOutflow(surfaceScalarField& alphaPhi) const;
//- Read phase properties dictionary
virtual bool read();
@ -217,44 +208,88 @@ public:
// Compressibility (variable density)
//- Return true if the phase is compressible otherwise false
virtual bool compressible() const;
virtual bool compressible() const = 0;
//- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual const tmp<volScalarField>& divU() const;
virtual tmp<volScalarField> divU() const = 0;
//- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual void divU(const tmp<volScalarField>& divU);
//- Return the phase kinetic energy
virtual const volScalarField& K() const;
virtual void divU(tmp<volScalarField> divU) = 0;
// Thermo
//- Return const access to the thermophysical model
//- Return whether the phase is isothermal
virtual bool isothermal() const = 0;
//- Return the enthalpy equation
virtual tmp<fvScalarMatrix> heEqn() = 0;
//- Return the thermophysical model
virtual const rhoThermo& thermo() const = 0;
//- Return non-const access to the thermophysical model
// for correction
virtual rhoThermo& thermo() = 0;
//- Access the thermophysical model
virtual rhoThermo& thermoRef() = 0;
//- Return the density field
virtual tmp<volScalarField> rho() const = 0;
//- Constant access the species mass fractions
// Species
//- Return whether the phase is pure (i.e., not multi-component)
virtual bool pure() const = 0;
//- Return the species fraction equation
virtual tmp<fvScalarMatrix> YiEqn(volScalarField& Yi) = 0;
//- Return the species mass fractions
virtual const PtrList<volScalarField>& Y() const = 0;
//- Access the species mass fractions
virtual PtrList<volScalarField>& Y() = 0;
virtual PtrList<volScalarField>& YRef() = 0;
//- Return the active species mass fractions
virtual const UPtrList<volScalarField>& YActive() const = 0;
//- Access the active species mass fractions
virtual UPtrList<volScalarField>& YActiveRef() = 0;
// Momentum
//- Constant access the velocity
//- Return whether the phase is stationary
virtual bool stationary() const = 0;
//- Return the momentum equation
virtual tmp<fvVectorMatrix> UEqn() = 0;
//- Return the momentum equation for the face-based algorithm
virtual tmp<fvVectorMatrix> UfEqn() = 0;
//- Return the velocity
virtual tmp<volVectorField> U() const = 0;
//- Access the velocity
virtual volVectorField& U() = 0;
virtual volVectorField& URef() = 0;
//- Return the volumetric flux
virtual tmp<surfaceScalarField> phi() const = 0;
//- Access the volumetric flux
virtual surfaceScalarField& phiRef() = 0;
//- Return the volumetric flux of the phase
virtual tmp<surfaceScalarField> alphaPhi() const = 0;
//- Access the volumetric flux of the phase
virtual surfaceScalarField& alphaPhiRef() = 0;
//- Return the mass flux of the phase
virtual tmp<surfaceScalarField> alphaRhoPhi() const = 0;
//- Access the mass flux of the phase
virtual surfaceScalarField& alphaRhoPhiRef() = 0;
//- Return the substantive acceleration
virtual tmp<volVectorField> DUDt() const = 0;
@ -262,29 +297,17 @@ public:
//- Return the substantive acceleration on the faces
virtual tmp<surfaceScalarField> DUDtf() const = 0;
//- Constant access the continuity error
//- Return the continuity error
virtual tmp<volScalarField> continuityError() const = 0;
//- Constant access the volumetric flux
virtual tmp<surfaceScalarField> phi() const = 0;
//- Return the continuity error due to the flow field
virtual tmp<volScalarField> continuityErrorFlow() const = 0;
//- Access the volumetric flux
virtual surfaceScalarField& phi() = 0;
//- Return the continuity error due to any sources
virtual tmp<volScalarField> continuityErrorSources() const = 0;
//- Constant access the volumetric flux of the phase
virtual tmp<surfaceScalarField> alphaPhi() const = 0;
//- Access the volumetric flux of the phase
virtual surfaceScalarField& alphaPhi() = 0;
//- Constant access the mass flux of the phase
virtual tmp<surfaceScalarField> alphaRhoPhi() const = 0;
//- Access the mass flux of the phase
virtual surfaceScalarField& alphaRhoPhi() = 0;
//- Ensure that the flux at inflow/outflow BCs is preserved
void correctInflowOutflow(surfaceScalarField& alphaPhi) const;
//- Return the phase kinetic energy
virtual tmp<volScalarField> K() const = 0;
// Transport

51
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModels.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -39,6 +39,7 @@ License
#include "InertPhaseModel.H"
#include "ReactingPhaseModel.H"
#include "MovingPhaseModel.H"
#include "StationaryPhaseModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -68,6 +69,30 @@ namespace Foam
purePhaseModel
);
typedef
AnisothermalPhaseModel
<
PurePhaseModel
<
InertPhaseModel
<
StationaryPhaseModel
<
ThermoPhaseModel<phaseModel, rhoThermo>
>
>
>
>
pureStationaryPhaseModel;
addNamedToRunTimeSelectionTable
(
phaseModel,
pureStationaryPhaseModel,
phaseSystem,
pureStationaryPhaseModel
);
typedef
IsothermalPhaseModel
<
@ -92,6 +117,30 @@ namespace Foam
pureIsothermalPhaseModel
);
typedef
IsothermalPhaseModel
<
PurePhaseModel
<
InertPhaseModel
<
StationaryPhaseModel
<
ThermoPhaseModel<phaseModel, rhoThermo>
>
>
>
>
pureStationaryIsothermalPhaseModel;
addNamedToRunTimeSelectionTable
(
phaseModel,
pureStationaryIsothermalPhaseModel,
phaseSystem,
pureStationaryIsothermalPhaseModel
);
typedef
AnisothermalPhaseModel
<

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -155,6 +155,50 @@ Foam::phaseSystem::phaseSystem
MRF_(mesh_)
{
// Groupings
label movingPhasei = 0;
label stationaryPhasei = 0;
label anisothermalPhasei = 0;
label multiComponentPhasei = 0;
forAll(phaseModels_, phasei)
{
phaseModel& phase = phaseModels_[phasei];
movingPhasei += !phase.stationary();
stationaryPhasei += phase.stationary();
anisothermalPhasei += !phase.isothermal();
multiComponentPhasei += !phase.pure();
}
movingPhaseModels_.resize(movingPhasei);
stationaryPhaseModels_.resize(stationaryPhasei);
anisothermalPhaseModels_.resize(anisothermalPhasei);
multiComponentPhaseModels_.resize(multiComponentPhasei);
movingPhasei = 0;
stationaryPhasei = 0;
anisothermalPhasei = 0;
multiComponentPhasei = 0;
forAll(phaseModels_, phasei)
{
phaseModel& phase = phaseModels_[phasei];
if (!phase.stationary())
{
movingPhaseModels_.set(movingPhasei ++, &phase);
}
if (phase.stationary())
{
stationaryPhaseModels_.set(stationaryPhasei ++, &phase);
}
if (!phase.isothermal())
{
anisothermalPhaseModels_.set(anisothermalPhasei ++, &phase);
}
if (!phase.pure())
{
multiComponentPhaseModels_.set(multiComponentPhasei ++, &phase);
}
}
// Write phi
phi_.writeOpt() = IOobject::AUTO_WRITE;
// Blending methods
@ -175,6 +219,7 @@ Foam::phaseSystem::phaseSystem
generatePairsAndSubModels("surfaceTension", surfaceTensionModels_);
generatePairsAndSubModels("aspectRatio", aspectRatioModels_);
// Update motion fields
correctKinematics();
}
@ -279,22 +324,18 @@ Foam::phaseSystem::sigma(const phasePairKey& key) const
}
void Foam::phaseSystem::solve()
{}
Foam::tmp<Foam::volScalarField> Foam::phaseSystem::dmdt
(
const Foam::phaseModel& phase
const phasePairKey& key
) const
{
tmp<volScalarField> tDmdt
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
IOobject::groupName("dmdt", phase.name()),
IOobject::groupName("dmdt", phasePairs_[key]->name()),
this->mesh_.time().timeName(),
this->mesh_
),
@ -302,11 +343,21 @@ const Foam::phaseModel& phase
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
return tDmdt;
}
Foam::Xfer<Foam::PtrList<Foam::volScalarField>> Foam::phaseSystem::dmdts() const
{
PtrList<volScalarField> dmdts(this->phaseModels_.size());
return dmdts.xfer();
}
void Foam::phaseSystem::solve()
{}
void Foam::phaseSystem::correct()
{
forAll(phaseModels_, phasei)

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

@ -73,71 +73,16 @@ public:
// Public typedefs
typedef
HashPtrTable
<
volScalarField,
phasePairKey,
phasePairKey::hash
>
KdTable;
typedef HashPtrTable<fvVectorMatrix> momentumTransferTable;
typedef
HashPtrTable
<
surfaceScalarField,
phasePairKey,
phasePairKey::hash
>
KdfTable;
typedef HashPtrTable<fvScalarMatrix> heatTransferTable;
typedef
HashPtrTable
<
volScalarField,
phasePairKey,
phasePairKey::hash
>
VmTable;
typedef
HashPtrTable
<
surfaceScalarField,
phasePairKey,
phasePairKey::hash
>
VmfTable;
typedef
HashPtrTable
<
fvVectorMatrix,
word,
string::hash
>
momentumTransferTable;
typedef
HashPtrTable
<
fvScalarMatrix,
word,
string::hash
>
heatTransferTable;
typedef
HashPtrTable
<
fvScalarMatrix,
word,
string::hash
>
massTransferTable;
typedef HashPtrTable<fvScalarMatrix> massTransferTable;
typedef PtrListDictionary<phaseModel> phaseModelList;
typedef UPtrList<phaseModel> phaseModelPartialList;
protected:
@ -182,6 +127,18 @@ protected:
//- Phase models
phaseModelList phaseModels_;
//- Moving phase models
phaseModelPartialList movingPhaseModels_;
//- Stationary phase models
phaseModelPartialList stationaryPhaseModels_;
//- Anisothermal phase models
phaseModelPartialList anisothermalPhaseModels_;
//- Multi-component phase models
phaseModelPartialList multiComponentPhaseModels_;
//- Phase pairs
phasePairTable phasePairs_;
@ -317,24 +274,6 @@ protected:
HashPtrTable<GeoField>& fieldTable
) const;
//- Fill up gaps in a phase-indexed list with zeros
template<class Type, template<class> class PatchField, class GeoMesh>
void fillFields
(
const word& name,
const dimensionSet& dims,
PtrList<GeometricField<Type, PatchField, GeoMesh>>& fieldList
) const;
//- Fill up gaps in a phase-indexed table with zeros
template<class Type, template<class> class PatchField, class GeoMesh>
void fillFields
(
const word& name,
const dimensionSet& dims,
HashPtrTable<GeometricField<Type, PatchField, GeoMesh>>& fieldTable
) const;
public:
@ -357,83 +296,230 @@ public:
// Member Functions
//- Constant access the mesh
inline const fvMesh& mesh() const;
// Access
//- Constant access the phase models
inline const phaseModelList& phases() const;
//- Return the mesh
inline const fvMesh& mesh() const;
//- Access the phase models
inline phaseModelList& phases();
//- Return the phase models
inline const phaseModelList& phases() const;
//- Constant access the phase pairs
inline const phasePairTable& phasePairs() const;
//- Access the phase models
inline phaseModelList& phases();
//- Return the mixture density
tmp<volScalarField> rho() const;
//- Return the models for phases that are moving
inline const phaseModelPartialList& movingPhases() const;
//- Return the mixture velocity
tmp<volVectorField> U() const;
//- Access the models for phases that are moving
inline phaseModelPartialList& movingPhases();
//- Constant access the mixture flux
inline const surfaceScalarField& phi() const;
//- Return the models for phases that are stationary
inline const phaseModelPartialList& stationaryPhases() const;
//- Access the mixture flux
inline surfaceScalarField& phi();
//- Access the models for phases that are stationary
inline phaseModelPartialList& stationaryPhases();
//- Constant access the rate of change of the pressure
inline const volScalarField& dpdt() const;
//- Return the models for phases that have variable temperature
inline const phaseModelPartialList& anisothermalPhases() const;
//- Access the rate of change of the pressure
inline volScalarField& dpdt();
//- Access the models for phases that have variable temperature
inline phaseModelPartialList& anisothermalPhases();
//- Return the aspect-ratio
tmp<volScalarField> E(const phasePairKey& key) const;
//- Return the models for phases that have multiple species
inline const phaseModelPartialList& multiComponentPhases() const;
//- Return the surface tension coefficient
tmp<volScalarField> sigma(const phasePairKey& key) const;
//- Access the models for phases that have multiple species
inline phaseModelPartialList& multiComponentPhases();
//- Return MRF zones
inline const IOMRFZoneList& MRF() const;
//- Return the phase pairs
inline const phasePairTable& phasePairs() const;
//- Optional FV-options
inline fv::options& fvOptions() const;
//- Return the mixture flux
inline const surfaceScalarField& phi() const;
//- Access a sub model between a phase pair
template<class modelType>
const modelType& lookupSubModel(const phasePair& key) const;
//- Access the mixture flux
inline surfaceScalarField& phi();
//- Access a sub model between two phases
template<class modelType>
const modelType& lookupSubModel
(
const phaseModel& dispersed,
const phaseModel& continuous
) const;
//- Return the rate of change of the pressure
inline const volScalarField& dpdt() const;
//- Solve for the phase fractions
virtual void solve();
//- Access the rate of change of the pressure
inline volScalarField& dpdt();
//- Return the total interfacial mass transfer rate for a phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
//- Return MRF zones
inline const IOMRFZoneList& MRF() const;
//- Correct the fluid properties other than the thermo and turbulence
virtual void correct();
//- Access the fvOptions
inline fv::options& fvOptions() const;
//- Correct the kinematics
virtual void correctKinematics();
//- Correct the thermodynamics
virtual void correctThermo();
// Sub-model lookup
//- Correct the turbulence
virtual void correctTurbulence();
//- Return a sub model between a phase pair
template<class modelType>
const modelType& lookupSubModel(const phasePair& key) const;
//- Correct the energy transport e.g. alphat
virtual void correctEnergyTransport();
//- Return a sub model between two phases
template<class modelType>
const modelType& lookupSubModel
(
const phaseModel& dispersed,
const phaseModel& continuous
) const;
//- Read base phaseProperties dictionary
virtual bool read();
// Field construction
//- Fill up gaps in a phase-indexed list of fields with zeros
template
<
class Type,
template<class> class PatchField,
class GeoMesh
>
void fillFields
(
const word& name,
const dimensionSet& dims,
PtrList<GeometricField<Type, PatchField, GeoMesh>>& fieldList
) const;
//- Fill up gaps in a phase-indexed table of fields with zeros
template
<
class Type,
template<class> class PatchField,
class GeoMesh
>
void fillFields
(
const word& name,
const dimensionSet& dims,
HashPtrTable<GeometricField<Type, PatchField, GeoMesh>>&
fieldTable
) const;
// Properties
//- Return the mixture density
tmp<volScalarField> rho() const;
//- Return the mixture velocity
tmp<volVectorField> U() const;
//- Return the aspect-ratio for a pair
tmp<volScalarField> E(const phasePairKey& key) const;
//- Return the surface tension coefficient for a pair
tmp<volScalarField> sigma(const phasePairKey& key) const;
//- Return the mass transfer rate for a pair
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
//- Return the mass transfer rates for each phase
virtual Xfer<PtrList<volScalarField>> dmdts() const;
// Transfers
//- Return the momentum transfer matrices for the cell-based
// algorithm
virtual autoPtr<momentumTransferTable> momentumTransfer() = 0;
//- Return the momentum transfer matrices for the face-based
// algiorithm
virtual autoPtr<momentumTransferTable> momentumTransferf() = 0;
//- Return the implicit force coefficients for the face-based
// algorithm
virtual Xfer<PtrList<surfaceScalarField>> AFfs() const = 0;
//- Return the force fluxes for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiFs
(
const PtrList<volScalarField>& rAUs
) = 0;
//- Return the force fluxes for the face-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiFfs
(
const PtrList<surfaceScalarField>& rAUfs
) = 0;
//- Return the force fluxes for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhis
(
const PtrList<volScalarField>& rAUs
) const = 0;
//- Return the force fluxes for the face-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhifs
(
const PtrList<surfaceScalarField>& rAUfs
) const = 0;
//- Return the explicit part of the drag force
virtual Xfer<PtrList<volVectorField>> KdUByAs
(
const PtrList<volScalarField>& rAUs
) const = 0;
//- Solve the drag system for the new velocities and fluxes
virtual void partialElimination
(
const PtrList<volScalarField>& rAUs
) = 0;
//- Solve the drag system for the new fluxes
virtual void partialEliminationf
(
const PtrList<surfaceScalarField>& rAUfs
) = 0;
//- Return the flux corrections for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> ddtCorrByAs
(
const PtrList<volScalarField>& rAUs,
const bool includeVirtualMass = false
) const = 0;
//- Return the phase diffusivities divided by the momentum
// coefficients
virtual const HashPtrTable<surfaceScalarField>& DByAfs() const = 0;
//- Return the heat transfer matrices
virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
//- Return the mass transfer matrices
virtual autoPtr<massTransferTable> massTransfer() const = 0;
// Evolution
//- Solve for the phase fractions
virtual void solve();
//- Correct the fluid properties other than those listed below
virtual void correct();
//- Correct the kinematics
virtual void correctKinematics();
//- Correct the thermodynamics
virtual void correctThermo();
//- Correct the turbulence
virtual void correctTurbulence();
//- Correct the energy transport e.g. alphat
virtual void correctEnergyTransport();
// IO
//- Read base phaseProperties dictionary
virtual bool read();
};

58
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2014-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -45,6 +45,62 @@ Foam::phaseSystem::phases()
}
inline const Foam::phaseSystem::phaseModelPartialList&
Foam::phaseSystem::movingPhases() const
{
return movingPhaseModels_;
}
inline Foam::phaseSystem::phaseModelPartialList&
Foam::phaseSystem::movingPhases()
{
return movingPhaseModels_;
}
inline const Foam::phaseSystem::phaseModelPartialList&
Foam::phaseSystem::stationaryPhases() const
{
return stationaryPhaseModels_;
}
inline Foam::phaseSystem::phaseModelPartialList&
Foam::phaseSystem::stationaryPhases()
{
return stationaryPhaseModels_;
}
inline const Foam::phaseSystem::phaseModelPartialList&
Foam::phaseSystem::anisothermalPhases() const
{
return anisothermalPhaseModels_;
}
inline Foam::phaseSystem::phaseModelPartialList&
Foam::phaseSystem::anisothermalPhases()
{
return anisothermalPhaseModels_;
}
inline const Foam::phaseSystem::phaseModelPartialList&
Foam::phaseSystem::multiComponentPhases() const
{
return multiComponentPhaseModels_;
}
inline Foam::phaseSystem::phaseModelPartialList&
Foam::phaseSystem::multiComponentPhases()
{
return multiComponentPhaseModels_;
}
inline const Foam::phaseSystem::phasePairTable&
Foam::phaseSystem::phasePairs() const
{

2
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C
View File

@ -281,6 +281,8 @@ void Foam::phaseSystem::addField
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::phaseSystem::fillFields
(

4
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C
View File

@ -87,10 +87,8 @@ void Foam::diameterModels::driftModels::constantDrift::driftRate
)
{
const sizeGroup& fi = *popBal_.sizeGroups()[i];
phaseModel& phase = const_cast<phaseModel&>(fi.phase());
volScalarField& rho = phase.thermo().rho();
driftRate += (popBal_.fluid().fvOptions()(phase, rho)&rho)/(N_*rho);
driftRate -= fi.phase().continuityErrorSources()/(N_*fi.phase().rho());
}

6
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/constantNucleation/constantNucleation.C
View File

@ -105,12 +105,10 @@ nucleationRate
)
{
const sizeGroup& fi = *popBal_.sizeGroups()[i];
phaseModel& phase = const_cast<phaseModel&>(fi.phase());
volScalarField& rho = phase.thermo().rho();
nucleationRate +=
nucleationRate -=
popBal_.gamma(i, velGroup_.formFactor()*pow3(d_))
*(popBal_.fluid().fvOptions()(phase, rho)&rho)/rho/fi.x();
*fi.phase().continuityErrorSources()/fi.phase().rho()/fi.x();
}

4
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C
View File

@ -151,8 +151,8 @@ nucleationRate
)
{
const sizeGroup& fi = *popBal_.sizeGroups()[i];
phaseModel& phase = const_cast<phaseModel&>(fi.phase());
volScalarField& rho = phase.thermo().rho();
const phaseModel& phase = fi.phase();
const volScalarField& rho = phase.rho();
const tmp<volScalarField> talphat(turbulence_.alphat());
const volScalarField::Boundary& alphatBf = talphat().boundaryField();

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

@ -7,32 +7,27 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
phaseSystem::heatTransferTable& heatTransfer = heatTransferPtr();
forAll(phases, phasei)
forAll(fluid.anisothermalPhases(), anisothermalPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.anisothermalPhases()[anisothermalPhasei];
const volScalarField& alpha = phase;
const volScalarField& rho = phase.rho();
const volVectorField& U = phase.U();
tmp<fvScalarMatrix> EEqn(phase.heEqn());
fvScalarMatrix EEqn
(
phase.heEqn()
==
*heatTransfer[phase.name()]
+ alpha*rho*(U&g)
+ fvOptions(alpha, rho, phase.thermoRef().he())
);
if (EEqn.valid())
{
EEqn =
(
EEqn
==
*heatTransfer[phase.name()]
+ alpha*rho*(U&g)
+ fvOptions(alpha, rho, phase.thermo().he())
);
EEqn->relax();
fvOptions.constrain(EEqn.ref());
EEqn->solve();
fvOptions.correct(phase.thermo().he());
}
EEqn.relax();
fvOptions.constrain(EEqn);
EEqn.solve();
fvOptions.correct(phase.thermoRef().he());
}
fluid.correctThermo();

29
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H
View File

@ -5,31 +5,26 @@
phaseSystem::massTransferTable&
massTransfer(massTransferPtr());
forAll(phases, phasei)
forAll(fluid.multiComponentPhases(), multiComponentPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.multiComponentPhases()[multiComponentPhasei];
PtrList<volScalarField>& Y = phase.Y();
UPtrList<volScalarField>& Y = phase.YActiveRef();
const volScalarField& alpha = phase;
const volScalarField& rho = phase.rho();
forAll(Y, i)
{
tmp<fvScalarMatrix> YiEqn(phase.YiEqn(Y[i]));
fvScalarMatrix YiEqn
(
phase.YiEqn(Y[i])
==
*massTransfer[Y[i].name()]
+ fvOptions(alpha, rho, Y[i])
);
if (YiEqn.valid())
{
YiEqn =
(
YiEqn
==
*massTransfer[Y[i].name()]
+ fvOptions(alpha, rho, Y[i])
);
YiEqn->relax();
YiEqn->solve(mesh.solver("Yi"));
}
YiEqn.relax();
YiEqn.solve(mesh.solver("Yi"));
}
}

2
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H
View File

@ -20,7 +20,7 @@ dimensionedScalar pMin
#include "gh.H"
volScalarField& p = phases[0].thermo().p();
volScalarField& p = phases[0].thermoRef().p();
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh

156
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
View File

@ -28,6 +28,7 @@ License
#include "MULES.H"
#include "subCycle.H"
#include "UniformField.H"
#include "fvcDdt.H"
#include "fvcDiv.H"
@ -74,32 +75,75 @@ void Foam::multiphaseSystem::solveAlphas()
phases()[phasei].correctBoundaryConditions();
}
PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
// Calculate the void fraction
volScalarField alphaVoid
(
IOobject
(
"alphaVoid",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar("alphaVoid", dimless, 1)
);
forAll(stationaryPhases(), stationaryPhasei)
{
alphaVoid -= stationaryPhases()[stationaryPhasei];
}
// Generate face-alphas
PtrList<surfaceScalarField> alphafs(phases().size());
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
volScalarField& alpha1 = phase;
alphaPhiCorrs.set
alphafs.set
(
phasei,
new surfaceScalarField
(
"phi" + alpha1.name() + "Corr",
fvc::flux
(
phi_,
phase,
"div(phi," + alpha1.name() + ')'
)
IOobject::groupName("alphaf", phase.name()),
upwind<scalar>(mesh_, phi_).interpolate(phase)
)
);
}
// Create correction fluxes
PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
forAll(stationaryPhases(), stationaryPhasei)
{
phaseModel& phase = stationaryPhases()[stationaryPhasei];
alphaPhiCorrs.set
(
phase.index(),
new surfaceScalarField
(
IOobject::groupName("alphaPhiCorr", phase.name()),
- upwind<scalar>(mesh_, phi_).flux(phase)
)
);
}
forAll(movingPhases(), movingPhasei)
{
phaseModel& phase = movingPhases()[movingPhasei];
volScalarField& alpha = phase;
alphaPhiCorrs.set
(
phase.index(),
new surfaceScalarField
(
IOobject::groupName("alphaPhiCorr", phase.name()),
fvc::flux(phi_, phase, "div(phi," + alpha.name() + ')')
)
);
surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phase.index()];
forAll(phases(), phasej)
forAll(phases(), phasei)
{
phaseModel& phase2 = phases()[phasej];
phaseModel& phase2 = phases()[phasei];
volScalarField& alpha2 = phase2;
if (&phase2 == &phase) continue;
@ -123,7 +167,7 @@ void Foam::multiphaseSystem::solveAlphas()
word phirScheme
(
"div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
"div(phir," + alpha2.name() + ',' + alpha.name() + ')'
);
alphaPhiCorr += fvc::flux
@ -144,33 +188,27 @@ void Foam::multiphaseSystem::solveAlphas()
alphaPhiCorr,
zeroField(),
zeroField(),
UniformField<scalar>(phase.alphaMax()),
min(alphaVoid.primitiveField(), phase.alphaMax())(),
zeroField(),
true
);
}
MULES::limitSum(alphaPhiCorrs);
volScalarField sumAlpha
(
IOobject
(
"sumAlpha",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar("sumAlpha", dimless, 0)
);
forAll(phases(), phasei)
// Limit the flux sums, fixing those of the stationary phases
labelHashSet fixedAlphaPhiCorrs;
forAll(stationaryPhases(), stationaryPhasei)
{
phaseModel& phase = phases()[phasei];
fixedAlphaPhiCorrs.insert(stationaryPhases()[stationaryPhasei].index());
}
MULES::limitSum(alphafs, alphaPhiCorrs, fixedAlphaPhiCorrs);
// Solve for the moving phase alphas
forAll(movingPhases(), movingPhasei)
{
phaseModel& phase = movingPhases()[movingPhasei];
volScalarField& alpha = phase;
surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()];
alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
phase.correctInflowOutflow(alphaPhi);
@ -252,29 +290,47 @@ void Foam::multiphaseSystem::solveAlphas()
Su
);
phase.alphaPhi() = alphaPhi;
phase.alphaPhiRef() = alphaPhi;
}
Info<< phase.name() << " volume fraction, min, max = "
// Report the phase fractions and the phase fraction sum
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
Info<< phase.name() << " fraction, min, max = "
<< phase.weightedAverage(mesh_.V()).value()
<< ' ' << min(phase).value()
<< ' ' << max(phase).value()
<< endl;
}
sumAlpha += phase;
volScalarField sumAlphaMoving
(
IOobject
(
"sumAlphaMoving",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar("sumAlphaMoving", dimless, 0)
);
forAll(movingPhases(), movingPhasei)
{
sumAlphaMoving += movingPhases()[movingPhasei];
}
Info<< "Phase-sum volume fraction, min, max = "
<< sumAlpha.weightedAverage(mesh_.V()).value()
<< ' ' << min(sumAlpha).value()
<< ' ' << max(sumAlpha).value()
<< (sumAlphaMoving + 1 - alphaVoid)().weightedAverage(mesh_.V()).value()
<< ' ' << min(sumAlphaMoving + 1 - alphaVoid).value()
<< ' ' << max(sumAlphaMoving + 1 - alphaVoid).value()
<< endl;
// Correct the sum of the phase-fractions to avoid 'drift'
volScalarField sumCorr(1 - sumAlpha);
forAll(phases(), phasei)
// Correct the sum of the phase fractions to avoid drift
forAll(movingPhases(), movingPhasei)
{
volScalarField& alpha = phases()[phasei];
alpha += alpha*sumCorr;
movingPhases()[movingPhasei] *= alphaVoid/sumAlphaMoving;
}
}
@ -643,9 +699,11 @@ void Foam::multiphaseSystem::solve()
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
if (phase.stationary()) continue;
volScalarField& alpha = phase;
phase.alphaPhi() = alphaPhiSums[phasei]/nAlphaSubCycles;
phase.alphaPhiRef() = alphaPhiSums[phasei]/nAlphaSubCycles;
// Correct the time index of the field
// to correspond to the global time
@ -664,9 +722,11 @@ void Foam::multiphaseSystem::solve()
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
phase.alphaRhoPhi() = fvc::interpolate(phase.rho())*phase.alphaPhi();
if (phase.stationary()) continue;
phase.alphaRhoPhiRef() =
fvc::interpolate(phase.rho())*phase.alphaPhi();
// Ensure the phase-fractions are bounded
phase.maxMin(0, 1);
}

75
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
View File

@ -114,9 +114,6 @@ private:
//- Return the virtual mass coefficient for phase pair
virtual tmp<volScalarField> Vm(const phasePairKey& key) const = 0;
//- Return the interfacial mass flow rate for phase pair
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
protected:
@ -168,77 +165,7 @@ public:
// Member Functions
using phaseSystem::sigma;
using phaseSystem::dmdt;
//- Return the momentum transfer matrices for the cell-based algorithm
virtual autoPtr<momentumTransferTable> momentumTransfer() const = 0;
//- Return the momentum transfer matrices for the face-based algiorithm
virtual autoPtr<momentumTransferTable> momentumTransferf() const = 0;
//- Return the implicit force coefficients for the face-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> AFfs() const = 0;
//- Return the force fluxes for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiFs
(
const PtrList<volScalarField>& rAUs
) = 0;
//- Return the force fluxes for the face-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiFfs
(
const PtrList<surfaceScalarField>& rAUfs
) = 0;
//- Return the force fluxes for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhis
(
const PtrList<volScalarField>& rAUs
) const = 0;
//- Return the force fluxes for the face-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhifs
(
const PtrList<surfaceScalarField>& rAUfs
) const = 0;
//- Return the explicit part of the drag force
virtual Xfer<PtrList<volVectorField>> KdUByAs
(
const PtrList<volScalarField>& rAUs
) const = 0;
//- Solve the drag system for the new velocities and fluxes
virtual void partialElimination
(
const PtrList<volScalarField>& rAUs
) = 0;
//- Solve the drag system for the new fluxes
virtual void partialEliminationf
(
const PtrList<surfaceScalarField>& rAUfs
) = 0;
//- Return the flux corrections for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> ddtCorrByAs
(
const PtrList<volScalarField>& rAUs,
const bool includeVirtualMass = false
) const = 0;
//- Return the phase diffusivities divided by the momentum coefficients
virtual const HashPtrTable<surfaceScalarField>& DByAfs() const = 0;
//- Return the heat transfer matrices
virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
//- Return the mass transfer matrices
virtual autoPtr<massTransferTable> massTransfer() const = 0;
//- Return true if there is mass transfer
virtual bool transfersMass() const = 0;
using phaseSystem::dmdts;
//- Return the surface tension force
tmp<surfaceScalarField> surfaceTension(const phaseModel& phase) const;

12
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H
View File

@ -9,17 +9,17 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
phaseSystem::momentumTransferTable&
momentumTransfer(momentumTransferPtr());
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
const volScalarField& alpha = phase;
const volScalarField& rho = phase.rho();
volVectorField& U = phase.U();
volVectorField& U = phase.URef();
UEqns.set
(
phasei,
phase.index(),
new fvVectorMatrix
(
phase.UEqn()
@ -29,8 +29,8 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
)
);
UEqns[phasei].relax();
fvOptions.constrain(UEqns[phasei]);
UEqns[phase.index()].relax();
fvOptions.constrain(UEqns[phase.index()]);
fvOptions.correct(U);
}
}

137
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
View File

@ -1,8 +1,5 @@
// Face volume fractions
PtrList<surfaceScalarField> alphafs(phases.size());
PtrList<volScalarField> rAUs(phases.size());
PtrList<surfaceScalarField> alpharAUfs(phases.size());
// Diagonal coefficients
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
@ -10,20 +7,37 @@ forAll(phases, phasei)
alphafs.set(phasei, fvc::interpolate(alpha).ptr());
alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
}
// Diagonal coefficients
PtrList<volScalarField> rAUs(phases.size());
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = fluid.movingPhases()[movingPhasei];
const volScalarField& alpha = phase;
rAUs.set
(
phasei,
phase.index(),
new volScalarField
(
IOobject::groupName("rAU", phase.name()),
1.0
/(
UEqns[phasei].A()
UEqns[phase.index()].A()
+ byDt(max(phase.residualAlpha() - alpha, scalar(0))*phase.rho())
)
)
);
}
fluid.fillFields("rAU", dimTime/dimDensity, rAUs);
// Phase diagonal coefficients
PtrList<surfaceScalarField> alpharAUfs(phases.size());
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
const volScalarField& alpha = phase;
alpharAUfs.set
(
@ -45,6 +59,13 @@ while (pimple.correct())
// Correct p_rgh for consistency with p and the updated densities
p_rgh = p - rho*gh;
// Correct fixed-flux BCs to be consistent with the velocity BCs
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = fluid.movingPhases()[movingPhasei];
MRF.correctBoundaryFlux(phase.U(), phase.phiRef());
}
// Combined buoyancy and force fluxes
PtrList<surfaceScalarField> phigFs(phases.size());
{
@ -85,17 +106,14 @@ while (pimple.correct())
// Correction force fluxes
PtrList<surfaceScalarField> ddtCorrByAs(fluid.ddtCorrByAs(rAUs));
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
const volScalarField& alpha = phase;
// Correct fixed-flux BCs to be consistent with the velocity BCs
MRF.correctBoundaryFlux(phase.U(), phase.phi());
HbyAs.set
(
phasei,
phase.index(),
new volVectorField
(
IOobject::groupName("HbyA", phase.name()),
@ -103,31 +121,33 @@ while (pimple.correct())
)
);
HbyAs[phasei] =
rAUs[phasei]
HbyAs[phase.index()] =
rAUs[phase.index()]
*(
UEqns[phasei].H()
UEqns[phase.index()].H()
+ byDt
(
max(phase.residualAlpha() - alpha, scalar(0))
*phase.rho()
)
*phase.U().oldTime()
*phase.U()().oldTime()
);
phiHbyAs.set
(
phasei,
phase.index(),
new surfaceScalarField
(
IOobject::groupName("phiHbyA", phase.name()),
fvc::flux(HbyAs[phasei])
- phigFs[phasei]
- ddtCorrByAs[phasei]
fvc::flux(HbyAs[phase.index()])
- phigFs[phase.index()]
- ddtCorrByAs[phase.index()]
)
);
}
}
fluid.fillFields("HbyA", dimVelocity, HbyAs);
fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs);
// Add explicit drag forces and fluxes if not doing partial elimination
if (!partialElimination)
@ -211,7 +231,8 @@ while (pimple.correct())
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField();
phib +=
alphafs[phasei].boundaryField()*phase.phi()().boundaryField();
}
setSnGrad<fixedFluxPressureFvPatchScalarField>
@ -225,11 +246,12 @@ while (pimple.correct())
// Compressible pressure equations
PtrList<fvScalarMatrix> pEqnComps(phases.size());
PtrList<volScalarField> dmdts(fluid.dmdts());
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
const volScalarField& alpha = phase;
volScalarField& rho = phase.thermo().rho();
volScalarField& rho = phase.thermoRef().rho();
if (phase.compressible())
{
@ -246,7 +268,7 @@ while (pimple.correct())
phasei,
(
(
phase.continuityError()
fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
- fvc::Sp
(
fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
@ -274,7 +296,7 @@ while (pimple.correct())
phasei,
(
(
phase.continuityError()
fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
- fvc::Sp
(
(fvc::ddt(alpha) + fvc::div(phase.alphaPhi())),
@ -287,27 +309,38 @@ while (pimple.correct())
);
}
}
else
{
pEqnComps.set
(
phasei,
fvm::Su(-(fvOptions(alpha, rho)&rho)/rho, p_rgh).ptr()
);
}
if (fluid.transfersMass())
// Add option sources
if (fvOptions.appliesToField(rho.name()))
{
tmp<fvScalarMatrix> optEqn = fvOptions(alpha, rho);
if (pEqnComps.set(phasei))
{
pEqnComps[phasei] -= fluid.dmdt(phase)/rho;
pEqnComps[phasei] -= (optEqn&rho)/rho;
}
else
{
pEqnComps.set
(
phasei,
fvm::Su(-fluid.dmdt(phase)/rho, p_rgh)
fvm::Su(- (optEqn&rho)/rho, p_rgh).ptr()
);
}
}
// Add mass transfer
if (dmdts.set(phasei))
{
if (pEqnComps.set(phasei))
{
pEqnComps[phasei] -= dmdts[phasei]/rho;
}
else
{
pEqnComps.set
(
phasei,
fvm::Su(- dmdts[phasei]/rho, p_rgh)
);
}
}
@ -351,16 +384,18 @@ while (pimple.correct())
surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp;
phase.phiRef() =
phiHbyAs[phase.index()]
+ alpharAUfs[phase.index()]*mSfGradp;
// Set the phase dilatation rates
if (pEqnComps.set(phasei))
if (pEqnComps.set(phase.index()))
{
phase.divU(-pEqnComps[phasei] & p_rgh);
phase.divU(-pEqnComps[phase.index()] & p_rgh);
}
}
@ -369,16 +404,16 @@ while (pimple.correct())
mSfGradp = pEqnIncomp.flux()/rAUf;
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
phase.U() =
HbyAs[phasei]
phase.URef() =
HbyAs[phase.index()]
+ fvc::reconstruct
(
alpharAUfs[phasei]*mSfGradp
- phigFs[phasei]
alpharAUfs[phase.index()]*mSfGradp
- phigFs[phase.index()]
);
}
@ -387,12 +422,12 @@ while (pimple.correct())
fluid.partialElimination(rAUs);
}
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
phase.U().correctBoundaryConditions();
fvOptions.correct(phase.U());
phase.URef().correctBoundaryConditions();
fvOptions.correct(phase.URef());
}
}
}
@ -407,7 +442,7 @@ while (pimple.correct())
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
phase.thermo().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
}
// Correct p_rgh for consistency with p and the updated densities

14
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H
View File

@ -12,17 +12,17 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
phaseSystem::momentumTransferTable&
momentumTransfer(momentumTransferPtr());
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
const volScalarField& alpha = phase;
volScalarField& rho = phase.thermo().rho();
volVectorField& U = phase.U();
const volScalarField& rho = phase.rho();
volVectorField& U = phase.URef();
UEqns.set
(
phasei,
phase.index(),
new fvVectorMatrix
(
phase.UfEqn()
@ -32,8 +32,8 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
)
);
UEqns[phasei].relax();
fvOptions.constrain(UEqns[phasei]);
UEqns[phase.index()].relax();
fvOptions.constrain(UEqns[phase.index()]);
U.correctBoundaryConditions();
fvOptions.correct(U);
}

169
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H
View File

@ -1,56 +1,68 @@
// Face volume fractions
PtrList<surfaceScalarField> alphafs(phases.size());
PtrList<surfaceScalarField> alphaRho0fs(phases.size());
PtrList<surfaceScalarField> rAUfs(phases.size());
PtrList<surfaceScalarField> alpharAUfs(phases.size());
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
const volScalarField& alpha = phase;
alphafs.set(phasei, fvc::interpolate(alpha).ptr());
alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
alphaRho0fs.set
(
phasei,
(
fvc::interpolate
(
max(alpha.oldTime(), phase.residualAlpha())
*phase.rho()().oldTime()
)
).ptr()
);
}
// Diagonal coefficients
PtrList<surfaceScalarField> rAUfs(phases.size());
{
PtrList<surfaceScalarField> AFfs(fluid.AFfs());
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
const volScalarField& alpha = phase;
alphafs.set(phasei, fvc::interpolate(alpha).ptr());
alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
alphaRho0fs.set
(
phasei,
(
fvc::interpolate
(
max(alpha.oldTime(), phase.residualAlpha())
*phase.rho()().oldTime()
)
).ptr()
);
phaseModel& phase = fluid.movingPhases()[movingPhasei];
rAUfs.set
(
phasei,
phase.index(),
new surfaceScalarField
(
IOobject::groupName("rAUf", phase.name()),
1.0
/(
byDt(alphaRho0fs[phasei])
+ fvc::interpolate(UEqns[phasei].A())
+ AFfs[phasei]
byDt(alphaRho0fs[phase.index()])
+ fvc::interpolate(UEqns[phase.index()].A())
+ AFfs[phase.index()]
)
)
);
alpharAUfs.set
(
phasei,
(
max(alphafs[phasei], phase.residualAlpha())*rAUfs[phasei]
).ptr()
);
}
}
fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs);
// Phase diagonal coefficients
PtrList<surfaceScalarField> alpharAUfs(phases.size());
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
alpharAUfs.set
(
phase.index(),
(
max(alphafs[phase.index()], phase.residualAlpha())
*rAUfs[phase.index()]
).ptr()
);
}
// Explicit force fluxes
PtrList<surfaceScalarField> phiFfs(fluid.phiFfs(rAUfs));
@ -63,12 +75,11 @@ while (pimple.correct())
// Correct p_rgh for consistency with p and the updated densities
p_rgh = p - rho*gh;
forAll(phases, phasei)
// Correct fixed-flux BCs to be consistent with the velocity BCs
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
// Correct fixed-flux BCs to be consistent with the velocity BCs
MRF.correctBoundaryFlux(phase.U(), phase.phi());
phaseModel& phase = fluid.movingPhases()[movingPhasei];
MRF.correctBoundaryFlux(phase.U(), phase.phiRef());
}
// Combined buoyancy and force fluxes
@ -106,27 +117,27 @@ while (pimple.correct())
// Predicted fluxes for each phase
PtrList<surfaceScalarField> phiHbyAs(phases.size());
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
phiHbyAs.set
(
phasei,
phase.index(),
new surfaceScalarField
(
IOobject::groupName("phiHbyA", phase.name()),
rAUfs[phasei]
rAUfs[phase.index()]
*(
fvc::flux(UEqns[phasei].H())
+ alphaRho0fs[phasei]
*byDt(MRF.absolute(phase.phi().oldTime()))
fvc::flux(UEqns[phase.index()].H())
+ alphaRho0fs[phase.index()]
*byDt(MRF.absolute(phase.phi()().oldTime()))
)
- phigFs[phasei]
- phigFs[phase.index()]
)
);
}
fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs);
// Add explicit drag forces and fluxes if not doing partial elimination
if (!partialElimination)
@ -205,7 +216,8 @@ while (pimple.correct())
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField();
phib +=
alphafs[phasei].boundaryField()*phase.phi()().boundaryField();
}
setSnGrad<fixedFluxPressureFvPatchScalarField>
@ -219,11 +231,12 @@ while (pimple.correct())
// Compressible pressure equations
PtrList<fvScalarMatrix> pEqnComps(phases.size());
PtrList<volScalarField> dmdts(fluid.dmdts());
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
const volScalarField& alpha = phase;
volScalarField& rho = phase.thermo().rho();
volScalarField& rho = phase.thermoRef().rho();
if (phase.compressible())
{
@ -240,7 +253,7 @@ while (pimple.correct())
phasei,
(
(
phase.continuityError()
fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
- fvc::Sp
(
fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
@ -263,6 +276,7 @@ while (pimple.correct())
(
pEqnComps[phasei].faceFluxCorrectionPtr()
);
pEqnComps[phasei].relax();
}
else
@ -272,7 +286,7 @@ while (pimple.correct())
phasei,
(
(
phase.continuityError()
fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
- fvc::Sp
(
(fvc::ddt(alpha) + fvc::div(phase.alphaPhi())),
@ -285,27 +299,36 @@ while (pimple.correct())
);
}
}
else
{
pEqnComps.set
(
phasei,
fvm::Su(-(fvOptions(alpha, rho)&rho)/rho, p_rgh).ptr()
);
}
if (fluid.transfersMass())
if (fvOptions.appliesToField(rho.name()))
{
tmp<fvScalarMatrix> optEqn = fvOptions(alpha, rho);
if (pEqnComps.set(phasei))
{
pEqnComps[phasei] -= fluid.dmdt(phase)/rho;
pEqnComps[phasei] -= (optEqn&rho)/rho;
}
else
{
pEqnComps.set
(
phasei,
fvm::Su(-fluid.dmdt(phase)/rho, p_rgh)
fvm::Su(- (optEqn&rho)/rho, p_rgh).ptr()
);
}
}
if (dmdts.set(phasei))
{
if (pEqnComps.set(phasei))
{
pEqnComps[phasei] -= dmdts[phasei]/rho;
}
else
{
pEqnComps.set
(
phasei,
fvm::Su(- dmdts[phasei]/rho, p_rgh)
);
}
}
@ -349,16 +372,18 @@ while (pimple.correct())
surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp;
phase.phiRef() =
phiHbyAs[phase.index()]
+ alpharAUfs[phase.index()]*mSfGradp;
// Set the phase dilatation rates
if (pEqnComps.set(phasei))
if (pEqnComps.set(phase.index()))
{
phase.divU(-pEqnComps[phasei] & p_rgh);
phase.divU(-pEqnComps[phase.index()] & p_rgh);
}
}
@ -372,13 +397,13 @@ while (pimple.correct())
mSfGradp = pEqnIncomp.flux()/rAUf;
forAll(phases, phasei)
forAll(fluid.movingPhases(), movingPhasei)
{
phaseModel& phase = phases[phasei];
phaseModel& phase = fluid.movingPhases()[movingPhasei];
phase.U() = fvc::reconstruct(MRF.absolute(phase.phi()));
phase.U().correctBoundaryConditions();
fvOptions.correct(phase.U());
phase.URef() = fvc::reconstruct(MRF.absolute(phase.phi()));
phase.URef().correctBoundaryConditions();
fvOptions.correct(phase.URef());
}
}
}
@ -393,7 +418,7 @@ while (pimple.correct())
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
phase.thermo().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
}
// Correct p_rgh for consistency with p and the updated densities

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

@ -8,46 +8,38 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
phaseSystem::heatTransferTable&
heatTransfer = heatTransferPtr();
if (!phase1.isothermal())
{
tmp<fvScalarMatrix> E1Eqn(phase1.heEqn());
fvScalarMatrix E1Eqn
(
phase1.heEqn()
==
*heatTransfer[phase1.name()]
+ alpha1*rho1*(U1&g)
+ fvOptions(alpha1, rho1, thermo1.he())
);
if (E1Eqn.valid())
{
E1Eqn =
(
E1Eqn
==
*heatTransfer[phase1.name()]
+ alpha1*rho1*(U1&g)
+ fvOptions(alpha1, rho1, phase1.thermo().he())
);
E1Eqn->relax();
fvOptions.constrain(E1Eqn.ref());
E1Eqn->solve();
fvOptions.correct(phase1.thermo().he());
}
E1Eqn.relax();
fvOptions.constrain(E1Eqn);
E1Eqn.solve();
fvOptions.correct(thermo1.he());
}
if (!phase2.isothermal())
{
tmp<fvScalarMatrix> E2Eqn(phase2.heEqn());
fvScalarMatrix E2Eqn
(
phase2.heEqn()
==
*heatTransfer[phase2.name()]
+ alpha2*rho2*(U2&g)
+ fvOptions(alpha2, rho2, phase2.thermoRef().he())
);
if (E2Eqn.valid())
{
E2Eqn =
(
E2Eqn
==
*heatTransfer[phase2.name()]
+ alpha2*rho2*(U2&g)
+ fvOptions(alpha2, rho2, phase2.thermo().he())
);
E2Eqn->relax();
fvOptions.constrain(E2Eqn.ref());
E2Eqn->solve();
fvOptions.correct(phase2.thermo().he());
}
E2Eqn.relax();
fvOptions.constrain(E2Eqn);
E2Eqn.solve();
fvOptions.correct(thermo2.he());
}
fluid.correctThermo();

31
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/YEqns.H
View File

@ -5,44 +5,41 @@
phaseSystem::massTransferTable&
massTransfer(massTransferPtr());
PtrList<volScalarField>& Y1 = phase1.Y();
PtrList<volScalarField>& Y2 = phase2.Y();
forAll(Y1, i)
if (!phase1.pure())
{
tmp<fvScalarMatrix> Y1iEqn(phase1.YiEqn(Y1[i]));
UPtrList<volScalarField>& Y1 = phase1.YActiveRef();
if (Y1iEqn.valid())
forAll(Y1, i)
{
Y1iEqn =
fvScalarMatrix Y1iEqn
(
Y1iEqn
phase1.YiEqn(Y1[i])
==
*massTransfer[Y1[i].name()]
+ fvOptions(alpha1, rho1, Y1[i])
);
Y1iEqn->relax();
Y1iEqn->solve(mesh.solver("Yi"));
Y1iEqn.relax();
Y1iEqn.solve(mesh.solver("Yi"));
}
}
forAll(Y2, i)
if (!phase2.pure())
{
tmp<fvScalarMatrix> Y2iEqn(phase2.YiEqn(Y2[i]));
UPtrList<volScalarField>& Y2 = phase2.YActiveRef();
if (Y2iEqn.valid())
forAll(Y2, i)
{
Y2iEqn =
fvScalarMatrix Y2iEqn
(
Y2iEqn
phase2.YiEqn(Y2[i])
==
*massTransfer[Y2[i].name()]
+ fvOptions(alpha2, rho2, Y2[i])
);
Y2iEqn->relax();
Y2iEqn->solve(mesh.solver("Yi"));
Y2iEqn.relax();
Y2iEqn.solve(mesh.solver("Yi"));
}
}

21
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFieldRefs.H
View File

@ -1,20 +1,21 @@
phaseModel& phase1 = fluid.phase1();
phaseModel& phase2 = fluid.phase2();
volScalarField& alpha1 = phase1;
volScalarField& alpha2 = phase2;
const volScalarField& alpha1 = phase1;
const volScalarField& alpha2 = phase2;
volVectorField& U1 = phase1.U();
surfaceScalarField& phi1 = phase1.phi();
surfaceScalarField& alphaPhi1 = phase1.alphaPhi();
volVectorField& U1 = phase1.URef();
surfaceScalarField& phi1 = phase1.phiRef();
const surfaceScalarField& alphaPhi1 = phase1.alphaPhi();
volVectorField& U2 = phase2.URef();
surfaceScalarField& phi2 = phase2.phiRef();
const surfaceScalarField& alphaPhi2 = phase2.alphaPhi();
volVectorField& U2 = phase2.U();
surfaceScalarField& phi2 = phase2.phi();
surfaceScalarField& alphaPhi2 = phase2.alphaPhi();
surfaceScalarField& phi = fluid.phi();
rhoThermo& thermo1 = phase1.thermo();
rhoThermo& thermo2 = phase2.thermo();
rhoThermo& thermo1 = phase1.thermoRef();
rhoThermo& thermo2 = phase2.thermoRef();
volScalarField& rho1 = thermo1.rho();
const volScalarField& psi1 = thermo1.psi();

2
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFields.H
View File

@ -19,7 +19,7 @@ dimensionedScalar pMin
#include "gh.H"
volScalarField& p = fluid.phase1().thermo().p();
volScalarField& p = fluid.phase1().thermoRef().p();
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh

75
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
View File

@ -184,7 +184,7 @@ while (pimple.correct())
pEqnComp1 =
(
phase1.continuityError()
fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ correction
@ -202,17 +202,12 @@ while (pimple.correct())
{
pEqnComp1 =
(
phase1.continuityError()
fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
}
}
else
{
pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh);
}
if (phase2.compressible())
{
if (pimple.transonic())
@ -225,7 +220,7 @@ while (pimple.correct())
pEqnComp2 =
(
phase2.continuityError()
fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ correction
@ -243,35 +238,65 @@ while (pimple.correct())
{
pEqnComp2 =
(
phase2.continuityError()
fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
}
}
else
// Add option sources
{
pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh);
if (fvOptions.appliesToField(rho1.name()))
{
tmp<fvScalarMatrix> optEqn1 = fvOptions(alpha1, rho1);
if (pEqnComp1.valid())
{
pEqnComp1.ref() -= (optEqn1 & rho1)/rho1;
}
else
{
pEqnComp1 = fvm::Su(- (optEqn1 & rho1)/rho1, p_rgh);
}
}
if (fvOptions.appliesToField(rho2.name()))
{
tmp<fvScalarMatrix> optEqn2 = fvOptions(alpha2, rho2);
if (pEqnComp2.valid())
{
pEqnComp2.ref() -= (optEqn2 & rho2)/rho2;
}
else
{
pEqnComp2 = fvm::Su(- (optEqn2 & rho2)/rho2, p_rgh);
}
}
}
if (fluid.transfersMass())
// Add mass transfer
{
if (pEqnComp1.valid())
PtrList<volScalarField> dmdts(fluid.dmdts());
if (dmdts.set(0))
{
pEqnComp1.ref() -= fluid.dmdt(phase1)/rho1;
if (pEqnComp1.valid())
{
pEqnComp1.ref() -= dmdts[0]/rho1;
}
else
{
pEqnComp1 = fvm::Su(- dmdts[0]/rho1, p_rgh);
}
}
else
if (dmdts.set(1))
{
pEqnComp1 = fvm::Su(-fluid.dmdt(phase1)/rho1, p_rgh);
}
if (pEqnComp2.valid())
{
pEqnComp2.ref() -= fluid.dmdt(phase2)/rho2;
}
else
{
pEqnComp2 = fvm::Su(-fluid.dmdt(phase2)/rho2, p_rgh);
if (pEqnComp2.valid())
{
pEqnComp2.ref() -= dmdts[1]/rho2;
}
else
{
pEqnComp2 = fvm::Su(- dmdts[1]/rho2, p_rgh);
}
}
}

75
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
View File

@ -196,7 +196,7 @@ while (pimple.correct())
pEqnComp1 =
(
phase1.continuityError()
fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ correction
@ -214,17 +214,12 @@ while (pimple.correct())
{
pEqnComp1 =
(
phase1.continuityError()
fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
}
}
else
{
pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh);
}
if (phase2.compressible())
{
if (pimple.transonic())
@ -237,7 +232,7 @@ while (pimple.correct())
pEqnComp2 =
(
phase2.continuityError()
fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ correction
@ -255,35 +250,65 @@ while (pimple.correct())
{
pEqnComp2 =
(
phase2.continuityError()
fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
}
}
else
// Add option sources
{
pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh);
if (fvOptions.appliesToField(rho1.name()))
{
tmp<fvScalarMatrix> optEqn1 = fvOptions(alpha1, rho1);
if (pEqnComp1.valid())
{
pEqnComp1.ref() -= (optEqn1 & rho1)/rho1;
}
else
{
pEqnComp1 = fvm::Su(- (optEqn1 & rho1)/rho1, p_rgh);
}
}
if (fvOptions.appliesToField(rho2.name()))
{
tmp<fvScalarMatrix> optEqn2 = fvOptions(alpha2, rho2);
if (pEqnComp2.valid())
{
pEqnComp2.ref() -= (optEqn2 & rho2)/rho2;
}
else
{
pEqnComp2 = fvm::Su(- (optEqn2 & rho2)/rho2, p_rgh);
}
}
}
if (fluid.transfersMass())
// Add mass transfer
{
if (pEqnComp1.valid())
PtrList<volScalarField> dmdts(fluid.dmdts());
if (dmdts.set(0))
{
pEqnComp1.ref() -= fluid.dmdt(phase1)/rho1;
if (pEqnComp1.valid())
{
pEqnComp1.ref() -= dmdts[0]/rho1;
}
else
{
pEqnComp1 = fvm::Su(- dmdts[0]/rho1, p_rgh);
}
}
else
if (dmdts.set(1))
{
pEqnComp1 = fvm::Su(-fluid.dmdt(phase1)/rho1, p_rgh);
}
if (pEqnComp2.valid())
{
pEqnComp2.ref() -= fluid.dmdt(phase2)/rho2;
}
else
{
pEqnComp2 = fvm::Su(-fluid.dmdt(phase2)/rho2, p_rgh);
if (pEqnComp2.valid())
{
pEqnComp2.ref() -= dmdts[1]/rho2;
}
else
{
pEqnComp2 = fvm::Su(- dmdts[1]/rho2, p_rgh);
}
}
}

19
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
View File

@ -29,6 +29,7 @@ License
#include "MULES.H"
#include "subCycle.H"
#include "UniformField.H"
#include "fvcDdt.H"
#include "fvcDiv.H"
@ -280,15 +281,15 @@ void Foam::twoPhaseSystem::solve()
if (alphaSubCycle.index() == 1)
{
phase1_.alphaPhi() = alphaPhi10;
phase1_.alphaPhiRef() = alphaPhi10;
}
else
{
phase1_.alphaPhi() += alphaPhi10;
phase1_.alphaPhiRef() += alphaPhi10;
}
}
phase1_.alphaPhi() /= nAlphaSubCycles;
phase1_.alphaPhiRef() /= nAlphaSubCycles;
}
else
{
@ -304,7 +305,7 @@ void Foam::twoPhaseSystem::solve()
zeroField()
);
phase1_.alphaPhi() = alphaPhi1;
phase1_.alphaPhiRef() = alphaPhi1;
}
if (alphaDbyA.valid())
@ -318,15 +319,15 @@ void Foam::twoPhaseSystem::solve()
alpha1Eqn.relax();
alpha1Eqn.solve();
phase1_.alphaPhi() += alpha1Eqn.flux();
phase1_.alphaPhiRef() += alpha1Eqn.flux();
}
phase1_.alphaRhoPhi() =
phase1_.alphaRhoPhiRef() =
fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
phase2_.correctInflowOutflow(phase2_.alphaPhi());
phase2_.alphaRhoPhi() =
phase2_.alphaPhiRef() = phi_ - phase1_.alphaPhi();
phase2_.correctInflowOutflow(phase2_.alphaPhiRef());
phase2_.alphaRhoPhiRef() =
fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
Info<< alpha1.name() << " volume fraction = "

81
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
View File

@ -66,9 +66,6 @@ private:
//- Return the virtual mass coefficient for phase pair
virtual tmp<volScalarField> Vm(const phasePairKey& key) const = 0;
//- Return the interfacial mass flow rate for phase pair
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
protected:
@ -126,90 +123,23 @@ public:
// Member Functions
using phaseSystem::sigma;
using phaseSystem::dmdt;
using phaseSystem::dmdts;
//- Constant access phase model 1
//- Return phase model 1
const phaseModel& phase1() const;
//- Access phase model 1
phaseModel& phase1();
//- Constant access phase model 2
//- Return phase model 2
const phaseModel& phase2() const;
//- Access phase model 2
phaseModel& phase2();
//- Constant access the phase not given as an argument
//- Return the phase not given as an argument
const phaseModel& otherPhase(const phaseModel& phase) const;
//- Return the momentum transfer matrices for the cell-based algorithm
virtual autoPtr<momentumTransferTable> momentumTransfer() const = 0;
//- Return the momentum transfer matrices for the face-based algorithm
virtual autoPtr<momentumTransferTable> momentumTransferf() const = 0;
//- Return the implicit force coefficients for the face-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> AFfs() const = 0;
//- Return the force fluxes for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiFs
(
const PtrList<volScalarField>& rAUs
) = 0;
//- Return the force fluxes for the face-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiFfs
(
const PtrList<surfaceScalarField>& rAUfs
) = 0;
//- Return the force fluxes for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhis
(
const PtrList<volScalarField>& rAUs
) const = 0;
//- Return the force fluxes for the face-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhifs
(
const PtrList<surfaceScalarField>& rAUfs
) const = 0;
//- Return the explicit part of the drag force
virtual Xfer<PtrList<volVectorField>> KdUByAs
(
const PtrList<volScalarField>& rAUs
) const = 0;
//- Solve the drag system for the new velocities and fluxes
virtual void partialElimination
(
const PtrList<volScalarField>& rAUs
) = 0;
//- Solve the drag system for the new fluxes
virtual void partialEliminationf
(
const PtrList<surfaceScalarField>& rAUfs
) = 0;
//- Return the flux corrections for the cell-based algorithm
virtual Xfer<PtrList<surfaceScalarField>> ddtCorrByAs
(
const PtrList<volScalarField>& rAUs,
const bool includeVirtualMass = false
) const = 0;
//- Return the phase diffusivities divided by the momentum coefficients
virtual const HashPtrTable<surfaceScalarField>& DByAfs() const = 0;
//- Return the heat transfer matrices
virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
//- Return the mass transfer matrices
virtual autoPtr<massTransferTable> massTransfer() const = 0;
//- Return the surface tension coefficient
tmp<volScalarField> sigma() const;
@ -222,9 +152,6 @@ public:
//- Return the virtual mass coefficient
tmp<volScalarField> Vm() const;
//- Return true if there is mass transfer
virtual bool transfersMass() const = 0;
//- Solve for the phase fractions
virtual void solve();
};