mirror of
https://github.com/OpenFOAM/OpenFOAM-6.git
synced 2025-12-08 06:57:46 +00:00
reactingEulerFoam: Multiphase partial elimination and re-organisation
Partial elimination has been implemented for the multiphase Euler-Euler solver. This does a linear solution of the drag system when calculating flux and velocity corrections after the solution of the pressure equation. This can improve the behaviour of the solution in the event that the drag coupling is high. It is controlled by means of a "partialElimination" switch within the PIMPLE control dictionary in fvSolution. A re-organisation has also been done in order to remove the exposure of the sub-modelling from the top-level solver. Rather than looping the drag, virtual mass, lift, etc..., models directly, the solver now calls a set of phase-system methods which group the different force terms. These new methods are documented in MomentumTransferPhaseSystem.H. Many other accessors have been removed as a consequence of this grouping. A bug was also fixed whereby the face-based algorithm was not transferring the momentum associated with a given interfacial mass transfer.
This commit is contained in:
Will Bainbridge
@ -132,10 +132,8 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
template<class BasePhaseSystem>
|
||||
bool Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::transfersMass
|
||||
(
|
||||
const phaseModel& phase
|
||||
) const
|
||||
bool
|
||||
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::transfersMass() const
|
||||
{
|
||||
return true;
|
||||
}
|
||||
@ -202,6 +200,17 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
|
||||
}
|
||||
|
||||
|
||||
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
|
||||
|
||||
@ -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
|
||||
@ -113,8 +113,8 @@ public:
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Return true if there is mass transfer for phase
|
||||
virtual bool transfersMass(const phaseModel& phase) const;
|
||||
//- 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;
|
||||
@ -126,6 +126,10 @@ public:
|
||||
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;
|
||||
|
||||
@ -50,10 +50,7 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::~HeatTransferPhaseSystem()
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
template<class BasePhaseSystem>
|
||||
bool Foam::HeatTransferPhaseSystem<BasePhaseSystem>::transfersMass
|
||||
(
|
||||
const phaseModel& phase
|
||||
) const
|
||||
bool Foam::HeatTransferPhaseSystem<BasePhaseSystem>::transfersMass() const
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
@ -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
|
||||
@ -89,21 +89,17 @@ public:
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Return true if there is mass transfer for phase
|
||||
virtual bool transfersMass(const phaseModel& phase) const;
|
||||
//- 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;
|
||||
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;
|
||||
virtual autoPtr<phaseSystem::heatTransferTable> heatTransfer() const;
|
||||
|
||||
//- Return the mass transfer matrices
|
||||
virtual autoPtr<phaseSystem::massTransferTable> massTransfer() const;
|
||||
|
||||
@ -27,6 +27,40 @@ License
|
||||
#include "interfaceCompositionModel.H"
|
||||
|
||||
|
||||
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
|
||||
|
||||
template<class BasePhaseSystem>
|
||||
void Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
|
||||
addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const
|
||||
{
|
||||
// Source term due to mass transfer
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::phasePairTable,
|
||||
this->phasePairs_,
|
||||
phasePairIter
|
||||
)
|
||||
{
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
template<class BasePhaseSystem>
|
||||
@ -225,33 +259,21 @@ momentumTransfer() const
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
eqnsPtr(BasePhaseSystem::momentumTransfer());
|
||||
|
||||
phaseSystem::momentumTransferTable& eqns = eqnsPtr();
|
||||
addMomentumTransfer(eqnsPtr());
|
||||
|
||||
// Source term due to mass transfer
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::phasePairTable,
|
||||
this->phasePairs_,
|
||||
phasePairIter
|
||||
)
|
||||
{
|
||||
const phasePair& pair(phasePairIter());
|
||||
return eqnsPtr;
|
||||
}
|
||||
|
||||
if (pair.ordered())
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
const volVectorField& U1(pair.phase1().U());
|
||||
const volVectorField& U2(pair.phase2().U());
|
||||
template<class BasePhaseSystem>
|
||||
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
|
||||
Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
|
||||
momentumTransferf() const
|
||||
{
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
eqnsPtr(BasePhaseSystem::momentumTransferf());
|
||||
|
||||
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);
|
||||
}
|
||||
addMomentumTransfer(eqnsPtr());
|
||||
|
||||
return eqnsPtr;
|
||||
}
|
||||
|
||||
@ -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
|
||||
@ -85,6 +85,14 @@ protected:
|
||||
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
|
||||
iDmdtExplicit_;
|
||||
|
||||
// Protected member functions
|
||||
|
||||
//- Add the mass-based momentum transfer to the equations
|
||||
void addMomentumTransfer
|
||||
(
|
||||
phaseSystem::momentumTransferTable& eqns
|
||||
) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
@ -119,6 +127,10 @@ public:
|
||||
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;
|
||||
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -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
|
||||
@ -122,6 +122,9 @@ private:
|
||||
//- Face virtual mass coefficients
|
||||
phaseSystem::VmfTable Vmfs_;
|
||||
|
||||
//- The phase diffusivities divided by the momentum coefficients
|
||||
HashPtrTable<surfaceScalarField> DByAfs_;
|
||||
|
||||
// Sub Models
|
||||
|
||||
//- Drag models
|
||||
@ -139,26 +142,17 @@ private:
|
||||
//- Turbulent dispersion models
|
||||
turbulentDispersionModelTable turbulentDispersionModels_;
|
||||
|
||||
//- Construct element phasei of Fs if not set and return
|
||||
// Used by Fs()
|
||||
volVectorField& setF
|
||||
(
|
||||
PtrList<volVectorField>& Fs, const label phasei
|
||||
) const;
|
||||
|
||||
//- Construct element phasei of Ffs if not set and return
|
||||
// Used by Ffs()
|
||||
surfaceScalarField& setFf
|
||||
(
|
||||
PtrList<surfaceScalarField>& Ffs, const label phasei
|
||||
) const;
|
||||
// Private member functions
|
||||
|
||||
//- Construct element phasei of phiDs if not set and return
|
||||
// Used by phiDs()
|
||||
surfaceScalarField& setPhiD
|
||||
(
|
||||
PtrList<surfaceScalarField>& phiDs, const label phasei
|
||||
) const;
|
||||
//- Return the drag coefficient for the phase pair
|
||||
virtual tmp<volScalarField> Kd(const phasePairKey& key) const;
|
||||
|
||||
//- Return the face drag coefficient for the phase pair
|
||||
virtual tmp<surfaceScalarField> Kdf(const phasePairKey& key) const;
|
||||
|
||||
//- Return the virtual mass coefficient for the phase pair
|
||||
virtual tmp<volScalarField> Vm(const phasePairKey& key) const;
|
||||
|
||||
|
||||
public:
|
||||
@ -175,82 +169,82 @@ public:
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Constant access to drag coefficients
|
||||
virtual const phaseSystem::KdTable& Kds() const
|
||||
{
|
||||
return Kds_;
|
||||
}
|
||||
//- 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;
|
||||
|
||||
//- Constant access to face drag coefficients
|
||||
virtual const phaseSystem::KdfTable& Kdfs() const
|
||||
{
|
||||
return Kdfs_;
|
||||
}
|
||||
//- As momentumTransfer, but for the face-based algorithm
|
||||
virtual autoPtr<phaseSystem::momentumTransferTable>
|
||||
momentumTransferf() const;
|
||||
|
||||
//- Constant access to virtual mass force coefficients
|
||||
virtual const phaseSystem::VmTable& Vms() const
|
||||
{
|
||||
return Vms_;
|
||||
}
|
||||
//- Return implicit force coefficients on the faces, for the face-based
|
||||
// algorithm.
|
||||
virtual Xfer<PtrList<surfaceScalarField>> AFfs() const;
|
||||
|
||||
//- Constant access to virtual mass force coefficients
|
||||
virtual const phaseSystem::VmfTable& Vmfs() const
|
||||
{
|
||||
return Vmfs_;
|
||||
}
|
||||
//- Return the explicit force fluxes for the cell-based algorithm, that
|
||||
// do not depend on phase mass/volume fluxes, and can therefore be
|
||||
// evaluated outside the corrector loop. This includes things like
|
||||
// lift, turbulent dispersion, and wall lubrication.
|
||||
virtual Xfer<PtrList<surfaceScalarField>> phiFs
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs
|
||||
);
|
||||
|
||||
//- Return the drag coefficient
|
||||
virtual tmp<volScalarField> Kd(const phasePairKey& key) const;
|
||||
//- As phiFs, but for the face-based algorithm
|
||||
virtual Xfer<PtrList<surfaceScalarField>> phiFfs
|
||||
(
|
||||
const PtrList<surfaceScalarField>& rAUfs
|
||||
);
|
||||
|
||||
//- Return the face drag coefficient
|
||||
virtual tmp<surfaceScalarField> Kdf(const phasePairKey& key) const;
|
||||
|
||||
//- Return the drag coefficient for phase
|
||||
virtual tmp<volScalarField> Kd(const phaseModel& phase) const;
|
||||
|
||||
//- Return the face drag coefficient for phase
|
||||
virtual tmp<surfaceScalarField> Kdf(const phaseModel& phase) const;
|
||||
|
||||
//- Return the virtual mass coefficient
|
||||
virtual tmp<volScalarField> Vm(const phasePairKey& key) const;
|
||||
|
||||
//- Return the face virtual mass coefficient
|
||||
virtual tmp<surfaceScalarField> Vmf(const phasePairKey& key) const;
|
||||
|
||||
//- Return the face virtual mass force coefficient for phase
|
||||
virtual tmp<surfaceScalarField> Vmf(const phaseModel& phase) const;
|
||||
|
||||
//- Return the combined force (lift + wall-lubrication)
|
||||
virtual tmp<volVectorField> F(const phasePairKey& key) const;
|
||||
|
||||
//- Return the combined force (lift + wall-lubrication)
|
||||
virtual autoPtr<PtrList<volVectorField>> Fs() const;
|
||||
|
||||
//- Return the combined force (lift + wall-lubrication)
|
||||
virtual autoPtr<PtrList<surfaceScalarField>> Ffs() const;
|
||||
|
||||
//- Return the turbulent dispersion force on faces for phase pair
|
||||
virtual autoPtr<PtrList<surfaceScalarField>> phiDs
|
||||
//- Return the explicit drag force fluxes for the cell-based algorithm.
|
||||
// These depend on phase mass/volume fluxes, and must therefore be
|
||||
// evaluated inside the corrector loop.
|
||||
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhis
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs
|
||||
) const;
|
||||
|
||||
//- Return the face based turbulent dispersion force for phase pair
|
||||
virtual autoPtr<PtrList<surfaceScalarField>> phiDfs
|
||||
//- As phiKdPhis, but for the face-based algorithm
|
||||
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhifs
|
||||
(
|
||||
const PtrList<surfaceScalarField>& rAUfs
|
||||
) const;
|
||||
|
||||
//- Return the combined face-force (lift + wall-lubrication)
|
||||
virtual tmp<surfaceScalarField> Ff(const phasePairKey& key) const;
|
||||
//- Return the explicit part of the drag force for the cell-based
|
||||
// algorithm. This is the cell-equivalent of phiKdPhis. These depend on
|
||||
// phase velocities, and must therefore be evaluated inside the
|
||||
// corrector loop.
|
||||
virtual Xfer<PtrList<volVectorField>> KdUByAs
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs
|
||||
) const;
|
||||
|
||||
//- Return the turbulent diffusivity
|
||||
// Multiplies the phase-fraction gradient
|
||||
virtual tmp<volScalarField> D(const phasePairKey& key) const;
|
||||
//- Solve the drag system for the velocities and fluxes
|
||||
virtual void partialElimination
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs
|
||||
);
|
||||
|
||||
//- Return the momentum transfer matrices
|
||||
virtual autoPtr<phaseSystem::momentumTransferTable>
|
||||
momentumTransfer() const;
|
||||
//- As partialElimination, but for the face-based algorithm. Only solves
|
||||
// for the fluxes.
|
||||
virtual void partialEliminationf
|
||||
(
|
||||
const PtrList<surfaceScalarField>& rAUfs
|
||||
);
|
||||
|
||||
//- Return the flux corrections for the cell-based algorithm. These
|
||||
// depend on phase mass/volume fluxes, and must therefore be evaluated
|
||||
// inside the corrector loop.
|
||||
virtual Xfer<PtrList<surfaceScalarField>> ddtCorrByAs
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs,
|
||||
const bool includeVirtualMass = false
|
||||
) const;
|
||||
|
||||
//- Return the phase diffusivities divided by the momentum coefficients
|
||||
virtual const HashPtrTable<surfaceScalarField>& DByAfs() const;
|
||||
|
||||
//- Read base phaseProperties dictionary
|
||||
virtual bool read();
|
||||
|
||||
@ -25,6 +25,41 @@ License
|
||||
|
||||
#include "PopulationBalancePhaseSystem.H"
|
||||
|
||||
|
||||
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
|
||||
|
||||
template<class BasePhaseSystem>
|
||||
void Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::
|
||||
addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const
|
||||
{
|
||||
// Source term due to mass transfer
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::phasePairTable,
|
||||
this->phasePairs_,
|
||||
phasePairIter
|
||||
)
|
||||
{
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
template<class BasePhaseSystem>
|
||||
@ -207,33 +242,20 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransfer() const
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
eqnsPtr(BasePhaseSystem::momentumTransfer());
|
||||
|
||||
phaseSystem::momentumTransferTable& eqns = eqnsPtr();
|
||||
addMomentumTransfer(eqnsPtr());
|
||||
|
||||
// Source term due to mass trasfer
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::phasePairTable,
|
||||
this->phasePairs_,
|
||||
phasePairIter
|
||||
)
|
||||
{
|
||||
const phasePair& pair(phasePairIter());
|
||||
return eqnsPtr;
|
||||
}
|
||||
|
||||
if (pair.ordered())
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
const volVectorField& U1(pair.phase1().U());
|
||||
const volVectorField& U2(pair.phase2().U());
|
||||
template<class BasePhaseSystem>
|
||||
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
|
||||
Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransferf() const
|
||||
{
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
eqnsPtr(BasePhaseSystem::momentumTransferf());
|
||||
|
||||
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);
|
||||
}
|
||||
addMomentumTransfer(eqnsPtr());
|
||||
|
||||
return eqnsPtr;
|
||||
}
|
||||
|
||||
@ -70,6 +70,15 @@ protected:
|
||||
HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash> pDmdt_;
|
||||
|
||||
|
||||
// Protected member functions
|
||||
|
||||
//- Add the mass-based momentum transfer to the equations
|
||||
void addMomentumTransfer
|
||||
(
|
||||
phaseSystem::momentumTransferTable& eqns
|
||||
) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
// Constructors
|
||||
@ -100,6 +109,10 @@ public:
|
||||
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;
|
||||
|
||||
|
||||
@ -28,6 +28,40 @@ License
|
||||
#include "fvcVolumeIntegrate.H"
|
||||
#include "fvmSup.H"
|
||||
|
||||
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
|
||||
|
||||
template<class BasePhaseSystem>
|
||||
void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
|
||||
addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const
|
||||
{
|
||||
// Source term due to mass trasfer
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::phasePairTable,
|
||||
this->phasePairs_,
|
||||
phasePairIter
|
||||
)
|
||||
{
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
template<class BasePhaseSystem>
|
||||
@ -317,33 +351,20 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::momentumTransfer() const
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
eqnsPtr(BasePhaseSystem::momentumTransfer());
|
||||
|
||||
phaseSystem::momentumTransferTable& eqns = eqnsPtr();
|
||||
addMomentumTransfer(eqnsPtr());
|
||||
|
||||
// Source term due to mass trasfer
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::phasePairTable,
|
||||
this->phasePairs_,
|
||||
phasePairIter
|
||||
)
|
||||
{
|
||||
const phasePair& pair(phasePairIter());
|
||||
return eqnsPtr;
|
||||
}
|
||||
|
||||
if (pair.ordered())
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
const volVectorField& U1(pair.phase1().U());
|
||||
const volVectorField& U2(pair.phase2().U());
|
||||
template<class BasePhaseSystem>
|
||||
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
|
||||
Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::momentumTransferf() const
|
||||
{
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
eqnsPtr(BasePhaseSystem::momentumTransferf());
|
||||
|
||||
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);
|
||||
}
|
||||
addMomentumTransfer(eqnsPtr());
|
||||
|
||||
return eqnsPtr;
|
||||
}
|
||||
|
||||
@ -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
|
||||
@ -82,6 +82,14 @@ protected:
|
||||
wMDotL_;
|
||||
|
||||
|
||||
// Protected member functions
|
||||
|
||||
//- Add the mass-based momentum transfer to the equations
|
||||
void addMomentumTransfer
|
||||
(
|
||||
phaseSystem::momentumTransferTable& eqns
|
||||
) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
@ -126,6 +134,10 @@ public:
|
||||
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;
|
||||
|
||||
|
||||
@ -164,17 +164,8 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
|
||||
fluid.mesh(),
|
||||
dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0)
|
||||
),
|
||||
DUDt_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
IOobject::groupName("DUDt", this->name()),
|
||||
fluid.mesh().time().timeName(),
|
||||
fluid.mesh()
|
||||
),
|
||||
fluid.mesh(),
|
||||
dimensionedVector("0", dimAcceleration, Zero)
|
||||
),
|
||||
DUDt_(nullptr),
|
||||
DUDtf_(nullptr),
|
||||
divU_(nullptr),
|
||||
turbulence_
|
||||
(
|
||||
@ -234,7 +225,17 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correctKinematics()
|
||||
{
|
||||
BasePhaseModel::correctKinematics();
|
||||
|
||||
DUDt_ = fvc::ddt(U_) + fvc::div(phi_, U_) - fvc::div(phi_)*U_;
|
||||
if (DUDt_.valid())
|
||||
{
|
||||
DUDt_.clear();
|
||||
DUDt();
|
||||
}
|
||||
|
||||
if (DUDtf_.valid())
|
||||
{
|
||||
DUDtf_.clear();
|
||||
DUDtf();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -260,39 +261,38 @@ template<class BasePhaseModel>
|
||||
Foam::tmp<Foam::fvVectorMatrix>
|
||||
Foam::MovingPhaseModel<BasePhaseModel>::UEqn()
|
||||
{
|
||||
const volScalarField& alpha = *this;
|
||||
volScalarField& rho = this->thermo().rho();
|
||||
|
||||
return
|
||||
(
|
||||
fvm::ddt(*this, this->thermo().rho(), U_)
|
||||
fvm::ddt(alpha, rho, U_)
|
||||
+ fvm::div(alphaRhoPhi_, U_)
|
||||
- fvm::Sp(continuityError_, U_)
|
||||
+ this->fluid().MRF().DDt(*this*this->thermo().rho(), U_)
|
||||
- fvm::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi_), U_)
|
||||
+ fvm::SuSp(this->fluid().fvOptions()(alpha, rho)&rho, U_)
|
||||
+ this->fluid().MRF().DDt(alpha*rho, U_)
|
||||
+ turbulence_->divDevRhoReff(U_)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
template<class BasePhaseModel>
|
||||
const Foam::surfaceScalarField&
|
||||
Foam::MovingPhaseModel<BasePhaseModel>::DbyA() const
|
||||
Foam::tmp<Foam::fvVectorMatrix>
|
||||
Foam::MovingPhaseModel<BasePhaseModel>::UfEqn()
|
||||
{
|
||||
if (DbyA_.valid())
|
||||
{
|
||||
return DbyA_;
|
||||
}
|
||||
else
|
||||
{
|
||||
return surfaceScalarField::null();
|
||||
}
|
||||
}
|
||||
// As the "normal" U-eqn but without the ddt terms
|
||||
|
||||
const volScalarField& alpha = *this;
|
||||
volScalarField& rho = this->thermo().rho();
|
||||
|
||||
template<class BasePhaseModel>
|
||||
void Foam::MovingPhaseModel<BasePhaseModel>::DbyA
|
||||
(
|
||||
const tmp<surfaceScalarField>& DbyA
|
||||
)
|
||||
{
|
||||
DbyA_ = DbyA;
|
||||
return
|
||||
(
|
||||
fvm::div(alphaRhoPhi_, U_)
|
||||
- fvm::Sp(fvc::div(alphaRhoPhi_), U_)
|
||||
+ fvm::SuSp(this->fluid().fvOptions()(alpha, rho)&rho, U_)
|
||||
+ this->fluid().MRF().DDt(alpha*rho, U_)
|
||||
+ turbulence_->divDevRhoReff(U_)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
@ -316,7 +316,25 @@ template<class BasePhaseModel>
|
||||
Foam::tmp<Foam::volVectorField>
|
||||
Foam::MovingPhaseModel<BasePhaseModel>::DUDt() const
|
||||
{
|
||||
return DUDt_;
|
||||
if (!DUDt_.valid())
|
||||
{
|
||||
DUDt_ = fvc::ddt(U_) + fvc::div(phi_, U_) - fvc::div(phi_)*U_;
|
||||
}
|
||||
|
||||
return tmp<volVectorField>(DUDt_());
|
||||
}
|
||||
|
||||
|
||||
template<class BasePhaseModel>
|
||||
Foam::tmp<Foam::surfaceScalarField>
|
||||
Foam::MovingPhaseModel<BasePhaseModel>::DUDtf() const
|
||||
{
|
||||
if (!DUDtf_.valid())
|
||||
{
|
||||
DUDtf_ = byDt(phi_ - phi_.oldTime());
|
||||
}
|
||||
|
||||
return tmp<surfaceScalarField>(DUDtf_());
|
||||
}
|
||||
|
||||
|
||||
@ -329,8 +347,7 @@ Foam::MovingPhaseModel<BasePhaseModel>::divU() const
|
||||
|
||||
|
||||
template<class BasePhaseModel>
|
||||
void
|
||||
Foam::MovingPhaseModel<BasePhaseModel>::divU
|
||||
void Foam::MovingPhaseModel<BasePhaseModel>::divU
|
||||
(
|
||||
const tmp<volScalarField>& divU
|
||||
)
|
||||
@ -475,12 +492,4 @@ Foam::MovingPhaseModel<BasePhaseModel>::pPrime() const
|
||||
}
|
||||
|
||||
|
||||
template<class BasePhaseModel>
|
||||
Foam::tmp<Foam::fvVectorMatrix>
|
||||
Foam::MovingPhaseModel<BasePhaseModel>::divDevRhoReff()
|
||||
{
|
||||
return turbulence_->divDevRhoReff(U_);
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -80,7 +80,10 @@ protected:
|
||||
surfaceScalarField alphaRhoPhi_;
|
||||
|
||||
//- Lagrangian acceleration field (needed for virtual-mass)
|
||||
volVectorField DUDt_;
|
||||
mutable tmp<volVectorField> DUDt_;
|
||||
|
||||
//- Lagrangian acceleration field on the faces (needed for virtual-mass)
|
||||
mutable tmp<surfaceScalarField> DUDtf_;
|
||||
|
||||
//- Dilatation rate
|
||||
tmp<volScalarField> divU_;
|
||||
@ -91,10 +94,6 @@ protected:
|
||||
//- Continuity error
|
||||
volScalarField continuityError_;
|
||||
|
||||
//- Phase diffusivity divided by the momentum coefficient.
|
||||
// Used for implicit treatment of the phase pressure and dispersion
|
||||
tmp<surfaceScalarField> DbyA_;
|
||||
|
||||
|
||||
private:
|
||||
|
||||
@ -137,14 +136,8 @@ public:
|
||||
//- Return the momentum equation
|
||||
virtual tmp<fvVectorMatrix> UEqn();
|
||||
|
||||
|
||||
// Implicit phase pressure and dispersion support
|
||||
|
||||
//- Return the phase diffusivity divided by the momentum coefficient
|
||||
virtual const surfaceScalarField& DbyA() const;
|
||||
|
||||
//- Set the phase diffusivity divided by the momentum coefficient
|
||||
virtual void DbyA(const tmp<surfaceScalarField>& DbyA);
|
||||
//- Return the momentum equation for the face-based algorithm
|
||||
virtual tmp<fvVectorMatrix> UfEqn();
|
||||
|
||||
|
||||
// Momentum
|
||||
@ -158,6 +151,9 @@ public:
|
||||
//- Return the substantive acceleration
|
||||
virtual tmp<volVectorField> DUDt() const;
|
||||
|
||||
//- 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;
|
||||
|
||||
@ -219,9 +215,6 @@ public:
|
||||
//- Return the phase-pressure'
|
||||
// (derivative of phase-pressure w.r.t. phase-fraction)
|
||||
virtual tmp<volScalarField> pPrime() const;
|
||||
|
||||
//- Return the turbulent term for the momentum equation
|
||||
virtual tmp<fvVectorMatrix> divDevRhoReff();
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -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
|
||||
@ -212,18 +212,4 @@ const Foam::volScalarField& Foam::phaseModel::K() const
|
||||
}
|
||||
|
||||
|
||||
const Foam::surfaceScalarField& Foam::phaseModel::DbyA() const
|
||||
{
|
||||
return surfaceScalarField::null();
|
||||
}
|
||||
|
||||
|
||||
void Foam::phaseModel::DbyA(const tmp<surfaceScalarField>& DbyA)
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Attempt to set the dilatation rate of an incompressible phase"
|
||||
<< endl;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -201,6 +201,9 @@ public:
|
||||
//- 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;
|
||||
|
||||
@ -226,15 +229,6 @@ public:
|
||||
virtual const volScalarField& K() const;
|
||||
|
||||
|
||||
// Implicit phase pressure and dispersion support
|
||||
|
||||
//- Return the phase diffusivity divided by the momentum coefficient
|
||||
virtual const surfaceScalarField& DbyA() const;
|
||||
|
||||
//- Set the phase diffusivity divided by the momentum coefficient
|
||||
virtual void DbyA(const tmp<surfaceScalarField>& DbyA);
|
||||
|
||||
|
||||
// Thermo
|
||||
|
||||
//- Return const access to the thermophysical model
|
||||
@ -265,6 +259,9 @@ public:
|
||||
//- Return the substantive acceleration
|
||||
virtual tmp<volVectorField> DUDt() const = 0;
|
||||
|
||||
//- Return the substantive acceleration on the faces
|
||||
virtual tmp<surfaceScalarField> DUDtf() const = 0;
|
||||
|
||||
//- Constant access the continuity error
|
||||
virtual tmp<volScalarField> continuityError() const = 0;
|
||||
|
||||
@ -380,9 +377,6 @@ public:
|
||||
//- Return the phase-pressure'
|
||||
// (derivative of phase-pressure w.r.t. phase-fraction)
|
||||
virtual tmp<volScalarField> pPrime() const = 0;
|
||||
|
||||
//- Return the turbulent term for the momentum equation
|
||||
virtual tmp<fvVectorMatrix> divDevRhoReff() = 0;
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -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
|
||||
@ -273,6 +273,68 @@ protected:
|
||||
>& models
|
||||
);
|
||||
|
||||
//- Add the field to a phase-indexed list, with the given name,
|
||||
// constructing if necessary
|
||||
template<class GeoField>
|
||||
void addField
|
||||
(
|
||||
const phaseModel& phase,
|
||||
const word& fieldName,
|
||||
tmp<GeoField> field,
|
||||
PtrList<GeoField>& fieldList
|
||||
) const;
|
||||
|
||||
//- Add the field to a phase-indexed list, with the given name,
|
||||
// constructing if necessary
|
||||
template<class GeoField>
|
||||
void addField
|
||||
(
|
||||
const phaseModel& phase,
|
||||
const word& fieldName,
|
||||
const GeoField& field,
|
||||
PtrList<GeoField>& fieldList
|
||||
) const;
|
||||
|
||||
//- Add the field to a phase-indexed table, with the given name,
|
||||
// constructing if necessary
|
||||
template<class GeoField>
|
||||
void addField
|
||||
(
|
||||
const phaseModel& phase,
|
||||
const word& fieldName,
|
||||
tmp<GeoField> field,
|
||||
HashPtrTable<GeoField>& fieldTable
|
||||
) const;
|
||||
|
||||
//- Add the field to a phase-indexed table, with the given name,
|
||||
// constructing if necessary
|
||||
template<class GeoField>
|
||||
void addField
|
||||
(
|
||||
const phaseModel& phase,
|
||||
const word& fieldName,
|
||||
const GeoField& field,
|
||||
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:
|
||||
|
||||
|
||||
@ -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
|
||||
@ -198,6 +198,161 @@ void Foam::phaseSystem::generatePairsAndSubModels
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField>
|
||||
void Foam::phaseSystem::addField
|
||||
(
|
||||
const phaseModel& phase,
|
||||
const word& fieldName,
|
||||
tmp<GeoField> field,
|
||||
PtrList<GeoField>& fieldList
|
||||
) const
|
||||
{
|
||||
if (fieldList.set(phase.index()))
|
||||
{
|
||||
fieldList[phase.index()] += field;
|
||||
}
|
||||
else
|
||||
{
|
||||
fieldList.set
|
||||
(
|
||||
phase.index(),
|
||||
new GeoField
|
||||
(
|
||||
IOobject::groupName(fieldName, phase.name()),
|
||||
field
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField>
|
||||
void Foam::phaseSystem::addField
|
||||
(
|
||||
const phaseModel& phase,
|
||||
const word& fieldName,
|
||||
const GeoField& field,
|
||||
PtrList<GeoField>& fieldList
|
||||
) const
|
||||
{
|
||||
addField(phase, fieldName, tmp<GeoField>(field), fieldList);
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField>
|
||||
void Foam::phaseSystem::addField
|
||||
(
|
||||
const phaseModel& phase,
|
||||
const word& fieldName,
|
||||
tmp<GeoField> field,
|
||||
HashPtrTable<GeoField>& fieldTable
|
||||
) const
|
||||
{
|
||||
if (fieldTable.found(phase.name()))
|
||||
{
|
||||
*fieldTable[phase.name()] += field;
|
||||
}
|
||||
else
|
||||
{
|
||||
fieldTable.set
|
||||
(
|
||||
phase.name(),
|
||||
new GeoField
|
||||
(
|
||||
IOobject::groupName(fieldName, phase.name()),
|
||||
field
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField>
|
||||
void Foam::phaseSystem::addField
|
||||
(
|
||||
const phaseModel& phase,
|
||||
const word& fieldName,
|
||||
const GeoField& field,
|
||||
HashPtrTable<GeoField>& fieldTable
|
||||
) const
|
||||
{
|
||||
addField(phase, fieldName, tmp<GeoField>(field), fieldTable);
|
||||
}
|
||||
|
||||
|
||||
template<class Type, template<class> class PatchField, class GeoMesh>
|
||||
void Foam::phaseSystem::fillFields
|
||||
(
|
||||
const word& name,
|
||||
const dimensionSet& dims,
|
||||
PtrList<GeometricField<Type, PatchField, GeoMesh>>& fieldList
|
||||
) const
|
||||
{
|
||||
forAll(this->phaseModels_, phasei)
|
||||
{
|
||||
if (fieldList.set(phasei))
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
const phaseModel& phase = this->phaseModels_[phasei];
|
||||
|
||||
fieldList.set
|
||||
(
|
||||
phasei,
|
||||
new GeometricField<Type, PatchField, GeoMesh>
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
IOobject::groupName(name, phase.name()),
|
||||
this->mesh_.time().timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensioned<Type>("zero", dims, pTraits<Type>::zero)
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class Type, template<class> class PatchField, class GeoMesh>
|
||||
void Foam::phaseSystem::fillFields
|
||||
(
|
||||
const word& name,
|
||||
const dimensionSet& dims,
|
||||
HashPtrTable<GeometricField<Type, PatchField, GeoMesh>>& fieldTable
|
||||
) const
|
||||
{
|
||||
forAll(this->phaseModels_, phasei)
|
||||
{
|
||||
const phaseModel& phase = this->phaseModels_[phasei];
|
||||
|
||||
if (fieldTable.set(phase.name()))
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
fieldTable.set
|
||||
(
|
||||
phase.name(),
|
||||
new GeometricField<Type, PatchField, GeoMesh>
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
IOobject::groupName(name, phase.name()),
|
||||
this->mesh_.time().timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensioned<Type>("zero", dims, pTraits<Type>::zero)
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class modelType>
|
||||
const modelType& Foam::phaseSystem::lookupSubModel(const phasePair& key) const
|
||||
{
|
||||
|
||||
@ -53,13 +53,21 @@ class multiphaseSystem
|
||||
:
|
||||
public phaseSystem
|
||||
{
|
||||
// Private data
|
||||
private:
|
||||
|
||||
volScalarField alphas_;
|
||||
// Private typedefs
|
||||
|
||||
typedef HashTable<scalar, phasePairKey, phasePairKey::hash>
|
||||
cAlphaTable;
|
||||
|
||||
|
||||
// Private data
|
||||
|
||||
//- The indexed phase-fraction field; e.g., 1 for water, 2 for air, 3
|
||||
// for oil, etc...
|
||||
volScalarField alphas_;
|
||||
|
||||
//-
|
||||
cAlphaTable cAlphas_;
|
||||
|
||||
//- Stabilisation for normalisation of the interface normal
|
||||
@ -100,27 +108,26 @@ class multiphaseSystem
|
||||
const phaseModel& alpha2
|
||||
) const;
|
||||
|
||||
|
||||
//- Return the drag coefficient for phase pair
|
||||
virtual tmp<volScalarField> Kd(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the face drag coefficient for phase pair
|
||||
virtual tmp<surfaceScalarField> Kdf(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the virtual mass coefficient for phase pair
|
||||
virtual tmp<volScalarField> Vm(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the face virtual mass coefficient for phase pair
|
||||
virtual tmp<surfaceScalarField> Vmf(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the turbulent diffusivity for phase pair
|
||||
// Multiplies the phase-fraction gradient
|
||||
virtual tmp<volScalarField> D(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the interfacial mass flow rate for phase pair
|
||||
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
// Protected data
|
||||
|
||||
//- Flag to indicate that returned lists of fields are "complete"; i.e.,
|
||||
// that an absence of force is returned as a constructed list of zeros,
|
||||
// rather than a null pointer
|
||||
static const bool fillFields_ = false;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
@ -160,54 +167,69 @@ public:
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Return the drag coefficient for all phase-pairs
|
||||
virtual const phaseSystem::KdTable& Kds() const = 0;
|
||||
using phaseSystem::sigma;
|
||||
using phaseSystem::dmdt;
|
||||
|
||||
//- Return the face drag coefficient for all phase-pairs
|
||||
virtual const phaseSystem::KdfTable& Kdfs() const = 0;
|
||||
//- Return the momentum transfer matrices for the cell-based algorithm
|
||||
virtual autoPtr<momentumTransferTable> momentumTransfer() const = 0;
|
||||
|
||||
//- Return the virtual mass force coefficient for all phase-pairs
|
||||
virtual const phaseSystem::VmTable& Vms() const = 0;
|
||||
//- Return the momentum transfer matrices for the face-based algiorithm
|
||||
virtual autoPtr<momentumTransferTable> momentumTransferf() const = 0;
|
||||
|
||||
//- Return the face based virtual mass force coefficient
|
||||
// for all phase-pairs
|
||||
virtual const phaseSystem::VmfTable& Vmfs() const = 0;
|
||||
//- Return the implicit force coefficients for the face-based algorithm
|
||||
virtual Xfer<PtrList<surfaceScalarField>> AFfs() const = 0;
|
||||
|
||||
//- Return the drag coefficient for phase
|
||||
virtual tmp<volScalarField> Kd(const phaseModel& phase) const = 0;
|
||||
//- Return the force fluxes for the cell-based algorithm
|
||||
virtual Xfer<PtrList<surfaceScalarField>> phiFs
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs
|
||||
) = 0;
|
||||
|
||||
//- Return the face based drag coefficient for phase
|
||||
virtual tmp<surfaceScalarField> Kdf(const phaseModel& phase) const = 0;
|
||||
//- Return the force fluxes for the face-based algorithm
|
||||
virtual Xfer<PtrList<surfaceScalarField>> phiFfs
|
||||
(
|
||||
const PtrList<surfaceScalarField>& rAUfs
|
||||
) = 0;
|
||||
|
||||
//- Return the face based virtual mass force coefficient for phase
|
||||
virtual tmp<surfaceScalarField> Vmf(const phaseModel& phase) const = 0;
|
||||
|
||||
//- Return the combined force (lift + wall-lubrication) for phase pair
|
||||
virtual autoPtr<PtrList<Foam::volVectorField>> Fs() const = 0;
|
||||
|
||||
//- Return the combined face force (lift + wall-lubrication)
|
||||
virtual autoPtr<PtrList<Foam::surfaceScalarField>> Ffs() const = 0;
|
||||
|
||||
//- Return the turbulent dispersion force on faces for phase pair
|
||||
virtual autoPtr<PtrList<Foam::surfaceScalarField>> phiDs
|
||||
//- Return the force fluxes for the cell-based algorithm
|
||||
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhis
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs
|
||||
) const = 0;
|
||||
|
||||
//- Return the face based turbulent dispersion force on faces
|
||||
// for phase pair
|
||||
virtual autoPtr<PtrList<Foam::surfaceScalarField>> phiDfs
|
||||
//- Return the force fluxes for the face-based algorithm
|
||||
virtual Xfer<PtrList<surfaceScalarField>> phiKdPhifs
|
||||
(
|
||||
const PtrList<surfaceScalarField>& rAUfs
|
||||
) const = 0;
|
||||
|
||||
//- Return true if there is mass transfer for phase
|
||||
virtual bool transfersMass(const phaseModel& phase) const = 0;
|
||||
//- Return the explicit part of the drag force
|
||||
virtual Xfer<PtrList<volVectorField>> KdUByAs
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs
|
||||
) const = 0;
|
||||
|
||||
using phaseSystem::dmdt;
|
||||
//- Solve the drag system for the new velocities and fluxes
|
||||
virtual void partialElimination
|
||||
(
|
||||
const PtrList<volScalarField>& rAUs
|
||||
) = 0;
|
||||
|
||||
//- Return the momentum transfer matrices
|
||||
virtual autoPtr<momentumTransferTable> momentumTransfer() const = 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;
|
||||
@ -215,6 +237,10 @@ public:
|
||||
//- Return the mass transfer matrices
|
||||
virtual autoPtr<massTransferTable> massTransfer() const = 0;
|
||||
|
||||
//- Return true if there is mass transfer
|
||||
virtual bool transfersMass() const = 0;
|
||||
|
||||
//- Return the surface tension force
|
||||
tmp<surfaceScalarField> surfaceTension(const phaseModel& phase) const;
|
||||
|
||||
//- Indicator of the proximity of the interface
|
||||
|
||||
@ -2,6 +2,7 @@ PtrList<surfaceScalarField> alphafs(phases.size());
|
||||
PtrList<volScalarField> rAUs(phases.size());
|
||||
PtrList<surfaceScalarField> alpharAUfs(phases.size());
|
||||
|
||||
// Diagonal coefficients
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
@ -33,57 +34,8 @@ forAll(phases, phasei)
|
||||
);
|
||||
}
|
||||
|
||||
// Lift, wall-lubrication and turbulent diffusion fluxes
|
||||
PtrList<surfaceScalarField> phiFs(phases.size());
|
||||
{
|
||||
autoPtr<PtrList<volVectorField>> Fs = fluid.Fs();
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
if (Fs().set(phasei))
|
||||
{
|
||||
phiFs.set
|
||||
(
|
||||
phasei,
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("phiF", phase.name()),
|
||||
fvc::flux(rAUs[phasei]*Fs()[phasei])
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
{
|
||||
autoPtr<PtrList<surfaceScalarField>> phiDs = fluid.phiDs(rAUs);
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
if (phiDs().set(phasei))
|
||||
{
|
||||
if (phiFs.set(phasei))
|
||||
{
|
||||
phiFs[phasei] += phiDs()[phasei];
|
||||
}
|
||||
else
|
||||
{
|
||||
phiFs.set
|
||||
(
|
||||
phasei,
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("phiF", phase.name()),
|
||||
phiDs()[phasei]
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Explicit force fluxes
|
||||
PtrList<surfaceScalarField> phiFs(fluid.phiFs(rAUs));
|
||||
|
||||
// --- Pressure corrector loop
|
||||
while (pimple.correct())
|
||||
@ -93,67 +45,104 @@ while (pimple.correct())
|
||||
// Correct p_rgh for consistency with p and the updated densities
|
||||
p_rgh = p - rho*gh;
|
||||
|
||||
PtrList<volVectorField> HbyAs(phases.size());
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
const volScalarField& alpha = phase;
|
||||
|
||||
// Correct fixed-flux BCs to be consistent with the velocity BCs
|
||||
MRF.correctBoundaryFlux(phase.U(), phase.phi());
|
||||
|
||||
HbyAs.set
|
||||
(
|
||||
phasei,
|
||||
new volVectorField
|
||||
(
|
||||
IOobject::groupName("HbyA", phase.name()),
|
||||
phase.U()
|
||||
)
|
||||
);
|
||||
|
||||
HbyAs[phasei] =
|
||||
rAUs[phasei]
|
||||
*(
|
||||
UEqns[phasei].H()
|
||||
+ byDt(max(phase.residualAlpha() - alpha, scalar(0))*phase.rho())
|
||||
*phase.U().oldTime()
|
||||
);
|
||||
}
|
||||
|
||||
surfaceScalarField ghSnGradRho
|
||||
(
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
|
||||
// Combined buoyancy and force fluxes
|
||||
PtrList<surfaceScalarField> phigFs(phases.size());
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
phigFs.set
|
||||
const surfaceScalarField ghSnGradRho
|
||||
(
|
||||
phasei,
|
||||
(
|
||||
alpharAUfs[phasei]
|
||||
*(
|
||||
ghSnGradRho
|
||||
- (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
|
||||
- fluid.surfaceTension(phase)*mesh.magSf()
|
||||
)
|
||||
).ptr()
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
|
||||
if (phiFs.set(phasei))
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phigFs[phasei] += phiFs[phasei];
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
phigFs.set
|
||||
(
|
||||
phasei,
|
||||
(
|
||||
alpharAUfs[phasei]
|
||||
*(
|
||||
ghSnGradRho
|
||||
- (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
|
||||
- fluid.surfaceTension(phase)*mesh.magSf()
|
||||
)
|
||||
).ptr()
|
||||
);
|
||||
|
||||
if (phiFs.set(phasei))
|
||||
{
|
||||
phigFs[phasei] += phiFs[phasei];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Predicted velocities and fluxes for each phase
|
||||
PtrList<volVectorField> HbyAs(phases.size());
|
||||
PtrList<surfaceScalarField> phiHbyAs(phases.size());
|
||||
{
|
||||
// Correction force fluxes
|
||||
PtrList<surfaceScalarField> ddtCorrByAs(fluid.ddtCorrByAs(rAUs));
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
const volScalarField& alpha = phase;
|
||||
|
||||
// Correct fixed-flux BCs to be consistent with the velocity BCs
|
||||
MRF.correctBoundaryFlux(phase.U(), phase.phi());
|
||||
|
||||
HbyAs.set
|
||||
(
|
||||
phasei,
|
||||
new volVectorField
|
||||
(
|
||||
IOobject::groupName("HbyA", phase.name()),
|
||||
phase.U()
|
||||
)
|
||||
);
|
||||
|
||||
HbyAs[phasei] =
|
||||
rAUs[phasei]
|
||||
*(
|
||||
UEqns[phasei].H()
|
||||
+ byDt
|
||||
(
|
||||
max(phase.residualAlpha() - alpha, scalar(0))
|
||||
*phase.rho()
|
||||
)
|
||||
*phase.U().oldTime()
|
||||
);
|
||||
|
||||
phiHbyAs.set
|
||||
(
|
||||
phasei,
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase.name()),
|
||||
fvc::flux(HbyAs[phasei])
|
||||
- phigFs[phasei]
|
||||
- ddtCorrByAs[phasei]
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
// Add explicit drag forces and fluxes if not doing partial elimination
|
||||
if (!partialElimination)
|
||||
{
|
||||
PtrList<volVectorField> KdUByAs(fluid.KdUByAs(rAUs));
|
||||
PtrList<surfaceScalarField> phiKdPhis(fluid.phiKdPhis(rAUs));
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
HbyAs[phasei] -= KdUByAs[phasei];
|
||||
phiHbyAs[phasei] -= phiKdPhis[phasei];
|
||||
}
|
||||
}
|
||||
|
||||
// Total predicted flux
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
IOobject
|
||||
@ -170,78 +159,20 @@ while (pimple.correct())
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
const volScalarField& alpha = phase;
|
||||
|
||||
// ddtPhiCorr filter -- only apply in pure(ish) phases
|
||||
surfaceScalarField alphafBar
|
||||
(
|
||||
fvc::interpolate(fvc::average(alphafs[phasei]))
|
||||
);
|
||||
surfaceScalarField phiCorrCoeff(pos0(alphafBar - 0.99));
|
||||
|
||||
surfaceScalarField::Boundary& phiCorrCoeffBf =
|
||||
phiCorrCoeff.boundaryFieldRef();
|
||||
|
||||
forAll(mesh.boundary(), patchi)
|
||||
{
|
||||
// Set ddtPhiCorr to 0 on non-coupled boundaries
|
||||
if
|
||||
(
|
||||
!mesh.boundary()[patchi].coupled()
|
||||
|| isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
|
||||
)
|
||||
{
|
||||
phiCorrCoeffBf[patchi] = 0;
|
||||
}
|
||||
}
|
||||
|
||||
phiHbyAs.set
|
||||
(
|
||||
phasei,
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase.name()),
|
||||
fvc::flux(HbyAs[phasei])
|
||||
+ phiCorrCoeff
|
||||
*fvc::interpolate
|
||||
(
|
||||
byDt(alpha.oldTime()*phase.rho()().oldTime()*rAUs[phasei])
|
||||
)
|
||||
*(
|
||||
MRF.absolute(phase.phi().oldTime())
|
||||
- fvc::flux(phase.U().oldTime())
|
||||
)
|
||||
- phigFs[phasei]
|
||||
)
|
||||
);
|
||||
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::KdTable,
|
||||
fluid.Kds(),
|
||||
KdIter
|
||||
)
|
||||
{
|
||||
const volScalarField& K(*KdIter());
|
||||
|
||||
const phasePair& pair(fluid.phasePairs()[KdIter.key()]);
|
||||
|
||||
if (pair.contains(phase))
|
||||
{
|
||||
const phaseModel& otherPhase = pair.otherPhase(phase);
|
||||
|
||||
phiHbyAs[phasei] +=
|
||||
fvc::interpolate(rAUs[phasei]*K)
|
||||
*MRF.absolute(otherPhase.phi());
|
||||
|
||||
HbyAs[phasei] += rAUs[phasei]*K*otherPhase.U();
|
||||
}
|
||||
}
|
||||
|
||||
phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
|
||||
}
|
||||
|
||||
// Add explicit drag fluxes if doing partial elimination
|
||||
if (partialElimination)
|
||||
{
|
||||
PtrList<surfaceScalarField> phiKdPhis(fluid.phiKdPhis(rAUs));
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phiHbyA -= alphafs[phasei]*phiKdPhis[phasei];
|
||||
}
|
||||
}
|
||||
|
||||
MRF.makeRelative(phiHbyA);
|
||||
|
||||
// Construct pressure "diffusivity"
|
||||
@ -261,8 +192,8 @@ while (pimple.correct())
|
||||
{
|
||||
rAUf += alphafs[phasei]*alpharAUfs[phasei];
|
||||
}
|
||||
rAUf = mag(rAUf);
|
||||
|
||||
rAUf = mag(rAUf);
|
||||
|
||||
// Update the fixedFluxPressure BCs to ensure flux consistency
|
||||
{
|
||||
@ -283,6 +214,7 @@ while (pimple.correct())
|
||||
);
|
||||
}
|
||||
|
||||
// Compressible pressure equations
|
||||
PtrList<fvScalarMatrix> pEqnComps(phases.size());
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
@ -355,7 +287,7 @@ while (pimple.correct())
|
||||
);
|
||||
}
|
||||
|
||||
if (fluid.transfersMass(phase))
|
||||
if (fluid.transfersMass())
|
||||
{
|
||||
if (pEqnComps.set(phasei))
|
||||
{
|
||||
@ -439,6 +371,17 @@ while (pimple.correct())
|
||||
alpharAUfs[phasei]*mSfGradp
|
||||
- phigFs[phasei]
|
||||
);
|
||||
}
|
||||
|
||||
if (partialElimination)
|
||||
{
|
||||
fluid.partialElimination(rAUs);
|
||||
}
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
phase.U().correctBoundaryConditions();
|
||||
fvOptions.correct(phase.U());
|
||||
}
|
||||
|
||||
@ -1,7 +0,0 @@
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
ddtPhis[phasei] =
|
||||
(phase.phi() - phase.phi().oldTime())/runTime.deltaT();
|
||||
}
|
||||
@ -3,26 +3,14 @@ Info<< "Constructing face momentum equations" << endl;
|
||||
PtrList<fvVectorMatrix> UEqns(phases.size());
|
||||
|
||||
{
|
||||
// Update of coefficients
|
||||
fluid.momentumTransfer();
|
||||
fluid.momentumTransfer(); // !!! Update coefficients shouldn't be necessary
|
||||
// This should be done on demand
|
||||
|
||||
PtrList<fvVectorMatrix> UgradUs(phases.size());
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
momentumTransferPtr(fluid.momentumTransferf());
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
UgradUs.set
|
||||
(
|
||||
phasei,
|
||||
new fvVectorMatrix
|
||||
(
|
||||
fvm::div(phase.phi(), phase.U())
|
||||
- fvm::Sp(fvc::div(phase.phi()), phase.U())
|
||||
+ MRF.DDt(phase.U())
|
||||
)
|
||||
);
|
||||
}
|
||||
phaseSystem::momentumTransferTable&
|
||||
momentumTransfer(momentumTransferPtr());
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
@ -30,49 +18,20 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
|
||||
|
||||
const volScalarField& alpha = phase;
|
||||
volScalarField& rho = phase.thermo().rho();
|
||||
const surfaceScalarField& alphaRhoPhi = phase.alphaRhoPhi();
|
||||
volVectorField& U = phase.U();
|
||||
|
||||
|
||||
UEqns.set
|
||||
(
|
||||
phasei,
|
||||
new fvVectorMatrix
|
||||
(
|
||||
fvm::div(alphaRhoPhi, U) - fvm::Sp(fvc::div(alphaRhoPhi), U)
|
||||
+ MRF.DDt(alpha*rho, U)
|
||||
+ phase.divDevRhoReff()
|
||||
+ fvm::SuSp(fvOptions(alpha, rho)&rho, U)
|
||||
phase.UfEqn()
|
||||
==
|
||||
fvOptions(alpha, rho, U)
|
||||
*momentumTransfer[phase.name()]
|
||||
+ fvOptions(alpha, rho, U)
|
||||
)
|
||||
);
|
||||
|
||||
// Add gradient part of virtual mass force
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::VmTable,
|
||||
fluid.Vms(),
|
||||
VmIter
|
||||
)
|
||||
{
|
||||
const volScalarField& Vm(*VmIter());
|
||||
|
||||
const phasePair& pair(fluid.phasePairs()[VmIter.key()]);
|
||||
|
||||
if (pair.contains(phase))
|
||||
{
|
||||
const phaseModel& otherPhase = pair.otherPhase(phase);
|
||||
|
||||
UEqns[phasei] +=
|
||||
Vm
|
||||
*(
|
||||
UgradUs[phasei]
|
||||
- (UgradUs[otherPhase.index()] & otherPhase.U())
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
UEqns[phasei].relax();
|
||||
fvOptions.constrain(UEqns[phasei]);
|
||||
U.correctBoundaryConditions();
|
||||
|
||||
@ -1,15 +0,0 @@
|
||||
PtrList<surfaceScalarField> ddtPhis(phases.size());
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
ddtPhis.set
|
||||
(
|
||||
phasei,
|
||||
(
|
||||
(phase.phi() - phase.phi().oldTime())/runTime.deltaT()
|
||||
).ptr()
|
||||
);
|
||||
}
|
||||
@ -3,102 +3,57 @@ PtrList<surfaceScalarField> alphaRho0fs(phases.size());
|
||||
PtrList<surfaceScalarField> rAUfs(phases.size());
|
||||
PtrList<surfaceScalarField> alpharAUfs(phases.size());
|
||||
|
||||
forAll(phases, phasei)
|
||||
// Diagonal coefficients
|
||||
{
|
||||
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()
|
||||
);
|
||||
|
||||
// add implicit part of drag and virtual mass force
|
||||
rAUfs.set
|
||||
(
|
||||
phasei,
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("rAUf", phase.name()),
|
||||
1.0
|
||||
/(
|
||||
byDt(alphaRho0fs[phasei] + fluid.Vmf(phase))
|
||||
+ fvc::interpolate(UEqns[phasei].A())
|
||||
+ fluid.Kdf(phase)
|
||||
)
|
||||
)
|
||||
);
|
||||
|
||||
alpharAUfs.set
|
||||
(
|
||||
phasei,
|
||||
(
|
||||
max(alphafs[phasei], phase.residualAlpha())*rAUfs[phasei]
|
||||
).ptr()
|
||||
);
|
||||
}
|
||||
|
||||
// Lift, wall-lubrication and turbulent diffusion fluxes
|
||||
PtrList<surfaceScalarField> phiFs(phases.size());
|
||||
{
|
||||
autoPtr<PtrList<surfaceScalarField>> Fs = fluid.Ffs();
|
||||
PtrList<surfaceScalarField> AFfs(fluid.AFfs());
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
const volScalarField& alpha = phase;
|
||||
|
||||
if (Fs().set(phasei))
|
||||
{
|
||||
phiFs.set
|
||||
alphafs.set(phasei, fvc::interpolate(alpha).ptr());
|
||||
alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
|
||||
|
||||
alphaRho0fs.set
|
||||
(
|
||||
phasei,
|
||||
(
|
||||
phasei,
|
||||
new surfaceScalarField
|
||||
fvc::interpolate
|
||||
(
|
||||
IOobject::groupName("phiF", phase.name()),
|
||||
rAUfs[phasei]*Fs()[phasei]
|
||||
max(alpha.oldTime(), phase.residualAlpha())
|
||||
*phase.rho()().oldTime()
|
||||
)
|
||||
);
|
||||
}
|
||||
).ptr()
|
||||
);
|
||||
|
||||
rAUfs.set
|
||||
(
|
||||
phasei,
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("rAUf", phase.name()),
|
||||
1.0
|
||||
/(
|
||||
byDt(alphaRho0fs[phasei])
|
||||
+ fvc::interpolate(UEqns[phasei].A())
|
||||
+ AFfs[phasei]
|
||||
)
|
||||
)
|
||||
);
|
||||
|
||||
alpharAUfs.set
|
||||
(
|
||||
phasei,
|
||||
(
|
||||
max(alphafs[phasei], phase.residualAlpha())*rAUfs[phasei]
|
||||
).ptr()
|
||||
);
|
||||
}
|
||||
}
|
||||
{
|
||||
autoPtr<PtrList<surfaceScalarField>> phiDs = fluid.phiDfs(rAUfs);
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
if (phiDs().set(phasei))
|
||||
{
|
||||
if (phiFs.set(phasei))
|
||||
{
|
||||
phiFs[phasei] += phiDs()[phasei];
|
||||
}
|
||||
else
|
||||
{
|
||||
phiFs.set
|
||||
(
|
||||
phasei,
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("phiF", phase.name()),
|
||||
phiDs()[phasei]
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Explicit force fluxes
|
||||
PtrList<surfaceScalarField> phiFfs(fluid.phiFfs(rAUfs));
|
||||
|
||||
// --- Pressure corrector loop
|
||||
while (pimple.correct())
|
||||
@ -116,52 +71,42 @@ while (pimple.correct())
|
||||
MRF.correctBoundaryFlux(phase.U(), phase.phi());
|
||||
}
|
||||
|
||||
surfaceScalarField ghSnGradRho
|
||||
(
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
|
||||
// Combined buoyancy and force fluxes
|
||||
PtrList<surfaceScalarField> phigFs(phases.size());
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
phigFs.set
|
||||
const surfaceScalarField ghSnGradRho
|
||||
(
|
||||
phasei,
|
||||
(
|
||||
alpharAUfs[phasei]
|
||||
*(
|
||||
ghSnGradRho
|
||||
- (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
|
||||
- fluid.surfaceTension(phase)*mesh.magSf()
|
||||
)
|
||||
).ptr()
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
|
||||
if (phiFs.set(phasei))
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phigFs[phasei] += phiFs[phasei];
|
||||
phaseModel& phase = phases[phasei];
|
||||
|
||||
phigFs.set
|
||||
(
|
||||
phasei,
|
||||
(
|
||||
alpharAUfs[phasei]
|
||||
*(
|
||||
ghSnGradRho
|
||||
- (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
|
||||
- fluid.surfaceTension(phase)*mesh.magSf()
|
||||
)
|
||||
).ptr()
|
||||
);
|
||||
|
||||
if (phiFfs.set(phasei))
|
||||
{
|
||||
phigFs[phasei] += phiFfs[phasei];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Predicted fluxes for each phase
|
||||
PtrList<surfaceScalarField> phiHbyAs(phases.size());
|
||||
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"phiHbyA",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh,
|
||||
dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
|
||||
);
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
@ -181,52 +126,53 @@ while (pimple.correct())
|
||||
- phigFs[phasei]
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
// Add explicit part of the drag
|
||||
forAllConstIter
|
||||
(
|
||||
phaseSystem::KdfTable,
|
||||
fluid.Kdfs(),
|
||||
KdfIter
|
||||
)
|
||||
// Add explicit drag forces and fluxes if not doing partial elimination
|
||||
if (!partialElimination)
|
||||
{
|
||||
PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
const surfaceScalarField& Kf(*KdfIter());
|
||||
|
||||
const phasePair& pair(fluid.phasePairs()[KdfIter.key()]);
|
||||
|
||||
if (pair.contains(phase))
|
||||
if (phiKdPhifs.set(phasei))
|
||||
{
|
||||
phiHbyAs[phasei] +=
|
||||
rAUfs[phasei]*Kf*MRF.absolute(pair.otherPhase(phase).phi());
|
||||
phiHbyAs[phasei] -= phiKdPhifs[phasei];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Add explicit part of the virtual mass force
|
||||
forAllConstIter
|
||||
// Total predicted flux
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
phaseSystem::VmfTable,
|
||||
fluid.Vmfs(),
|
||||
VmfIter
|
||||
)
|
||||
{
|
||||
const surfaceScalarField& Vmf(*VmfIter());
|
||||
|
||||
const phasePair& pair(fluid.phasePairs()[VmfIter.key()]);
|
||||
|
||||
if (pair.contains(phase))
|
||||
{
|
||||
phiHbyAs[phasei] +=
|
||||
rAUfs[phasei]
|
||||
*(
|
||||
Vmf*byDt(MRF.absolute(phase.phi().oldTime()))
|
||||
+ Vmf*ddtPhis[pair.otherPhase(phase).index()]
|
||||
);
|
||||
}
|
||||
}
|
||||
"phiHbyA",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh,
|
||||
dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
|
||||
);
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
|
||||
}
|
||||
|
||||
// Add explicit drag fluxes if doing partial elimination
|
||||
if (partialElimination)
|
||||
{
|
||||
PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phiHbyA -= alphafs[phasei]*phiKdPhifs[phasei];
|
||||
}
|
||||
}
|
||||
|
||||
MRF.makeRelative(phiHbyA);
|
||||
|
||||
// Construct pressure "diffusivity"
|
||||
@ -246,8 +192,8 @@ while (pimple.correct())
|
||||
{
|
||||
rAUf += alphafs[phasei]*alpharAUfs[phasei];
|
||||
}
|
||||
rAUf = mag(rAUf);
|
||||
|
||||
rAUf = mag(rAUf);
|
||||
|
||||
// Update the fixedFluxPressure BCs to ensure flux consistency
|
||||
{
|
||||
@ -268,6 +214,7 @@ while (pimple.correct())
|
||||
);
|
||||
}
|
||||
|
||||
// Compressible pressure equations
|
||||
PtrList<fvScalarMatrix> pEqnComps(phases.size());
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
@ -344,7 +291,7 @@ while (pimple.correct())
|
||||
);
|
||||
}
|
||||
|
||||
if (fluid.transfersMass(phase))
|
||||
if (fluid.transfersMass())
|
||||
{
|
||||
if (pEqnComps.set(phasei))
|
||||
{
|
||||
@ -412,6 +359,11 @@ while (pimple.correct())
|
||||
}
|
||||
}
|
||||
|
||||
if (partialElimination)
|
||||
{
|
||||
fluid.partialEliminationf(rAUfs);
|
||||
}
|
||||
|
||||
// Optionally relax pressure for velocity correction
|
||||
p_rgh.relax();
|
||||
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -64,16 +64,10 @@ int main(int argc, char *argv[])
|
||||
(
|
||||
pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
|
||||
);
|
||||
|
||||
// Switch implicitPhasePressure
|
||||
// (
|
||||
// mesh.solverDict(alpha1.name()).lookupOrDefault<Switch>
|
||||
// (
|
||||
// "implicitPhasePressure", false
|
||||
// )
|
||||
// );
|
||||
|
||||
#include "pUf/createDDtU.H"
|
||||
Switch partialElimination
|
||||
(
|
||||
pimple.dict().lookupOrDefault<Switch>("partialElimination", false)
|
||||
);
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -114,7 +108,6 @@ int main(int argc, char *argv[])
|
||||
#include "pUf/UEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pUf/pEqn.H"
|
||||
#include "pUf/DDtU.H"
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
@ -7,12 +7,10 @@ volScalarField& alpha2 = phase2;
|
||||
volVectorField& U1 = phase1.U();
|
||||
surfaceScalarField& phi1 = phase1.phi();
|
||||
surfaceScalarField& alphaPhi1 = phase1.alphaPhi();
|
||||
surfaceScalarField& alphaRhoPhi1 = phase1.alphaRhoPhi();
|
||||
|
||||
volVectorField& U2 = phase2.U();
|
||||
surfaceScalarField& phi2 = phase2.phi();
|
||||
surfaceScalarField& alphaPhi2 = phase2.alphaPhi();
|
||||
surfaceScalarField& alphaRhoPhi2 = phase2.alphaRhoPhi();
|
||||
surfaceScalarField& phi = fluid.phi();
|
||||
|
||||
rhoThermo& thermo1 = phase1.thermo();
|
||||
|
||||
@ -1,24 +1,33 @@
|
||||
const surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1));
|
||||
const surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
|
||||
|
||||
const volScalarField rAU1
|
||||
PtrList<volScalarField> rAUs;
|
||||
rAUs.append
|
||||
(
|
||||
IOobject::groupName("rAU", phase1.name()),
|
||||
1.0
|
||||
/(
|
||||
U1Eqn.A()
|
||||
+ byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
|
||||
new volScalarField
|
||||
(
|
||||
IOobject::groupName("rAU", phase1.name()),
|
||||
1.0
|
||||
/(
|
||||
U1Eqn.A()
|
||||
+ byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
|
||||
)
|
||||
)
|
||||
);
|
||||
const volScalarField rAU2
|
||||
rAUs.append
|
||||
(
|
||||
IOobject::groupName("rAU", phase2.name()),
|
||||
1.0
|
||||
/(
|
||||
U2Eqn.A()
|
||||
+ byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
|
||||
new volScalarField
|
||||
(
|
||||
IOobject::groupName("rAU", phase2.name()),
|
||||
1.0
|
||||
/(
|
||||
U2Eqn.A()
|
||||
+ byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
|
||||
)
|
||||
)
|
||||
);
|
||||
const volScalarField& rAU1 = rAUs[0];
|
||||
const volScalarField& rAU2 = rAUs[1];
|
||||
|
||||
const surfaceScalarField alpharAUf1
|
||||
(
|
||||
@ -29,55 +38,17 @@ const surfaceScalarField alpharAUf2
|
||||
fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2)
|
||||
);
|
||||
|
||||
// Drag coefficients
|
||||
const volScalarField Kd(fluid.Kd());
|
||||
const volScalarField Vm(fluid.Vm());
|
||||
|
||||
// Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
|
||||
tmp<surfaceScalarField> phiF1;
|
||||
tmp<surfaceScalarField> phiF2;
|
||||
{
|
||||
// Turbulent-dispersion diffusivity
|
||||
const volScalarField D(fluid.D());
|
||||
|
||||
// Phase-1 turbulent dispersion and particle-pressure flux
|
||||
tmp<surfaceScalarField> DbyA1
|
||||
(
|
||||
fvc::interpolate
|
||||
(
|
||||
rAU1*(D + phase1.pPrime())
|
||||
)
|
||||
);
|
||||
|
||||
// Phase-2 turbulent dispersion and particle-pressure flux
|
||||
tmp<surfaceScalarField> DbyA2
|
||||
(
|
||||
fvc::interpolate
|
||||
(
|
||||
rAU2*(D + phase2.pPrime())
|
||||
)
|
||||
);
|
||||
|
||||
// Lift and wall-lubrication forces
|
||||
const volVectorField F(fluid.F());
|
||||
|
||||
// Phase-fraction face-gradient
|
||||
const surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
|
||||
|
||||
// Phase-1 dispersion, lift and wall-lubrication flux
|
||||
phiF1 = DbyA1()*snGradAlpha1 + fvc::flux(rAU1*F);
|
||||
|
||||
// Phase-2 dispersion, lift and wall-lubrication flux
|
||||
phiF2 = - DbyA2()*snGradAlpha1 - fvc::flux(rAU2*F);
|
||||
|
||||
// Cache the phase diffusivities for implicit treatment in the
|
||||
// phase-fraction equation
|
||||
if (implicitPhasePressure)
|
||||
{
|
||||
phase1.DbyA(DbyA1);
|
||||
phase2.DbyA(DbyA2);
|
||||
}
|
||||
}
|
||||
const volScalarField rAUKd1(rAU1*Kd);
|
||||
const volScalarField rAUKd2(rAU2*Kd);
|
||||
const surfaceScalarField phiKd1(fvc::interpolate(rAUKd1));
|
||||
const surfaceScalarField phiKd2(fvc::interpolate(rAUKd2));
|
||||
|
||||
// Explicit force fluxes
|
||||
PtrList<surfaceScalarField> phiFs(fluid.phiFs(rAUs));
|
||||
const surfaceScalarField& phiF1 = phiFs[0];
|
||||
const surfaceScalarField& phiF2 = phiFs[1];
|
||||
|
||||
// --- Pressure corrector loop
|
||||
while (pimple.correct())
|
||||
@ -91,6 +62,34 @@ while (pimple.correct())
|
||||
MRF.correctBoundaryFlux(U1, phi1);
|
||||
MRF.correctBoundaryFlux(U2, phi2);
|
||||
|
||||
// Combined buoyancy and force fluxes
|
||||
const surfaceScalarField ghSnGradRho
|
||||
(
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
|
||||
const surfaceScalarField phigF1
|
||||
(
|
||||
alpharAUf1
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
|
||||
)
|
||||
+ phiF1
|
||||
);
|
||||
|
||||
const surfaceScalarField phigF2
|
||||
(
|
||||
alpharAUf2
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
|
||||
)
|
||||
+ phiF2
|
||||
);
|
||||
|
||||
// Predicted velocities
|
||||
volVectorField HbyA1
|
||||
(
|
||||
IOobject::groupName("HbyA", phase1.name()),
|
||||
@ -117,115 +116,40 @@ while (pimple.correct())
|
||||
*U2.oldTime()
|
||||
);
|
||||
|
||||
const surfaceScalarField ghSnGradRho
|
||||
(
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
// Correction force fluxes
|
||||
PtrList<surfaceScalarField> ddtCorrByAs(fluid.ddtCorrByAs(rAUs));
|
||||
|
||||
const surfaceScalarField phig1
|
||||
(
|
||||
alpharAUf1
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
|
||||
)
|
||||
);
|
||||
|
||||
const surfaceScalarField phig2
|
||||
(
|
||||
alpharAUf2
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
|
||||
)
|
||||
);
|
||||
|
||||
|
||||
// ddtPhiCorr filter -- only apply in pure(ish) phases
|
||||
const surfaceScalarField alphaf1Bar
|
||||
(
|
||||
fvc::interpolate(fvc::average(alphaf1))
|
||||
);
|
||||
surfaceScalarField phiCorrCoeff1(pos0(alphaf1Bar - 0.99));
|
||||
surfaceScalarField phiCorrCoeff2(pos0(0.01 - alphaf1Bar));
|
||||
|
||||
{
|
||||
surfaceScalarField::Boundary& phiCorrCoeff1Bf =
|
||||
phiCorrCoeff1.boundaryFieldRef();
|
||||
|
||||
surfaceScalarField::Boundary& phiCorrCoeff2Bf =
|
||||
phiCorrCoeff2.boundaryFieldRef();
|
||||
|
||||
forAll(mesh.boundary(), patchi)
|
||||
{
|
||||
// Set ddtPhiCorr to 0 on non-coupled boundaries
|
||||
if
|
||||
(
|
||||
!mesh.boundary()[patchi].coupled()
|
||||
|| isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
|
||||
)
|
||||
{
|
||||
phiCorrCoeff1Bf[patchi] = 0;
|
||||
phiCorrCoeff2Bf[patchi] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const surfaceScalarField phi1Corr
|
||||
(
|
||||
MRF.absolute(phi1.oldTime()) - fvc::flux(U1.oldTime())
|
||||
);
|
||||
|
||||
const surfaceScalarField phi2Corr
|
||||
(
|
||||
MRF.absolute(phi2.oldTime()) - fvc::flux(U2.oldTime())
|
||||
);
|
||||
|
||||
// Phase-1 predicted flux
|
||||
// Predicted fluxes
|
||||
const surfaceScalarField phiHbyA1
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase1.name()),
|
||||
fvc::flux(HbyA1)
|
||||
+ phiCorrCoeff1
|
||||
*fvc::interpolate(byDt(alpha1.oldTime()*rho1.oldTime()*rAU1))*phi1Corr
|
||||
+ fvc::interpolate(Vm*byDt(rAU1))
|
||||
*(
|
||||
phi1Corr + (MRF.absolute(phi2) - fvc::flux(U2)) - phi2Corr
|
||||
)
|
||||
- phiF1()
|
||||
- phig1
|
||||
fvc::flux(HbyA1) - phigF1 - ddtCorrByAs[0]
|
||||
);
|
||||
|
||||
// Phase-2 predicted flux
|
||||
const surfaceScalarField phiHbyA2
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase2.name()),
|
||||
fvc::flux(HbyA2)
|
||||
+ phiCorrCoeff2
|
||||
*fvc::interpolate(byDt(alpha2.oldTime()*rho2.oldTime()*rAU2))*phi2Corr
|
||||
+ fvc::interpolate(Vm*byDt(rAU2))
|
||||
*(
|
||||
phi2Corr + (MRF.absolute(phi1) - fvc::flux(U1)) - phi1Corr
|
||||
)
|
||||
- phiF2()
|
||||
- phig2
|
||||
fvc::flux(HbyA2) - phigF2 - ddtCorrByAs[1]
|
||||
);
|
||||
|
||||
// Face-drag coefficients
|
||||
const surfaceScalarField rAUKd1(fvc::interpolate(rAU1*Kd));
|
||||
const surfaceScalarField rAUKd2(fvc::interpolate(rAU2*Kd));
|
||||
ddtCorrByAs.clear();
|
||||
|
||||
// Construct the mean predicted flux
|
||||
// including explicit drag contributions based on absolute fluxes
|
||||
// Drag fluxes
|
||||
PtrList<surfaceScalarField> phiKdPhis(fluid.phiKdPhis(rAUs));
|
||||
const surfaceScalarField& phiKdPhi1 = phiKdPhis[0];
|
||||
const surfaceScalarField& phiKdPhi2 = phiKdPhis[1];
|
||||
|
||||
// Total predicted flux
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
"phiHbyA",
|
||||
alphaf1*(phiHbyA1 + rAUKd1*MRF.absolute(phi2))
|
||||
+ alphaf2*(phiHbyA2 + rAUKd2*MRF.absolute(phi1))
|
||||
alphaf1*(phiHbyA1 - phiKdPhi1) + alphaf2*(phiHbyA2 - phiKdPhi2)
|
||||
);
|
||||
|
||||
MRF.makeRelative(phiHbyA);
|
||||
|
||||
phiKdPhis.clear();
|
||||
|
||||
// Construct pressure "diffusivity"
|
||||
const surfaceScalarField rAUf
|
||||
(
|
||||
@ -246,11 +170,8 @@ while (pimple.correct())
|
||||
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
|
||||
);
|
||||
|
||||
tmp<fvScalarMatrix> pEqnComp1;
|
||||
tmp<fvScalarMatrix> pEqnComp2;
|
||||
|
||||
// Construct the compressibility parts of the pressure equation
|
||||
|
||||
tmp<fvScalarMatrix> pEqnComp1, pEqnComp2;
|
||||
if (phase1.compressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
@ -412,8 +333,8 @@ while (pimple.correct())
|
||||
|
||||
surfaceScalarField phir
|
||||
(
|
||||
((phi1s + rAUKd1*phi2s) - (phi2s + rAUKd2*phi1s))
|
||||
/(1 - rAUKd1*rAUKd2)
|
||||
((phi1s + phiKd1*phi2s) - (phi2s + phiKd2*phi1s))
|
||||
/(1 - phiKd1*phiKd2)
|
||||
);
|
||||
|
||||
phi1 = phi + alphaf2*phir;
|
||||
@ -440,25 +361,23 @@ while (pimple.correct())
|
||||
const volVectorField Us1
|
||||
(
|
||||
HbyA1
|
||||
+ fvc::reconstruct(alpharAUf1*mSfGradp - phiF1() - phig1)
|
||||
+ fvc::reconstruct(alpharAUf1*mSfGradp - phigF1)
|
||||
);
|
||||
|
||||
const volVectorField Us2
|
||||
(
|
||||
HbyA2
|
||||
+ fvc::reconstruct(alpharAUf2*mSfGradp - phiF2() - phig2)
|
||||
+ fvc::reconstruct(alpharAUf2*mSfGradp - phigF2)
|
||||
);
|
||||
|
||||
const volScalarField D1(rAU1*Kd);
|
||||
const volScalarField D2(rAU2*Kd);
|
||||
|
||||
const volVectorField U
|
||||
(
|
||||
alpha1*(Us1 + D1*U2) + alpha2*(Us2 + D2*U1)
|
||||
alpha1*(Us1 + rAUKd1*U2) + alpha2*(Us2 + rAUKd2*U1)
|
||||
);
|
||||
|
||||
const volVectorField Ur
|
||||
(
|
||||
((1 - D2)*Us1 - (1 - D1)*Us2)/(1 - D1*D2)
|
||||
((1 - rAUKd2)*Us1 - (1 - rAUKd1)*Us2)/(1 - rAUKd1*rAUKd2)
|
||||
);
|
||||
|
||||
U1 = U + alpha2*Ur;
|
||||
|
||||
@ -1,2 +0,0 @@
|
||||
tddtPhi1 = fvc::ddt(phi1);
|
||||
tddtPhi2 = fvc::ddt(phi2);
|
||||
@ -4,33 +4,22 @@ fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime);
|
||||
fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
|
||||
|
||||
{
|
||||
volScalarField Vm(fluid.Vm());
|
||||
fluid.momentumTransfer(); // !!! Update coefficients shouldn't be necessary
|
||||
// This should be done on demand
|
||||
|
||||
fvVectorMatrix UgradU1
|
||||
(
|
||||
fvm::div(phi1, U1) - fvm::Sp(fvc::div(phi1), U1)
|
||||
+ MRF.DDt(U1)
|
||||
);
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
momentumTransferPtr(fluid.momentumTransferf());
|
||||
|
||||
fvVectorMatrix UgradU2
|
||||
(
|
||||
fvm::div(phi2, U2) - fvm::Sp(fvc::div(phi2), U2)
|
||||
+ MRF.DDt(U2)
|
||||
);
|
||||
|
||||
const volScalarField dmdt21(posPart(fluid.dmdt(phase1)));
|
||||
const volScalarField dmdt12(negPart(fluid.dmdt(phase1)));
|
||||
phaseSystem::momentumTransferTable&
|
||||
momentumTransfer(momentumTransferPtr());
|
||||
|
||||
{
|
||||
U1Eqn =
|
||||
(
|
||||
fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1)
|
||||
+ fvm::Sp(dmdt21, U1) - dmdt21*U2
|
||||
+ MRF.DDt(alpha1*rho1, U1)
|
||||
+ phase1.divDevRhoReff()
|
||||
+ Vm*(UgradU1 - (UgradU2 & U2))
|
||||
- fvOptions(alpha1, rho1, U1)
|
||||
+ fvm::SuSp(fvOptions(alpha1, rho1)&rho1, U1)
|
||||
phase1.UfEqn()
|
||||
==
|
||||
*momentumTransfer[phase1.name()]
|
||||
+ fvOptions(alpha1, rho1, U1)
|
||||
);
|
||||
U1Eqn.relax();
|
||||
fvOptions.constrain(U1Eqn);
|
||||
@ -41,13 +30,10 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
|
||||
{
|
||||
U2Eqn =
|
||||
(
|
||||
fvm::div(alphaRhoPhi2, U2) - fvm::Sp(fvc::div(alphaRhoPhi2), U2)
|
||||
- fvm::Sp(dmdt12, U2) + dmdt12*U1
|
||||
+ MRF.DDt(alpha2*rho2, U2)
|
||||
+ phase2.divDevRhoReff()
|
||||
+ Vm*(UgradU2 - (UgradU1 & U1))
|
||||
- fvOptions(alpha2, rho2, U2)
|
||||
+ fvm::SuSp(fvOptions(alpha2, rho2)&rho2, U2)
|
||||
phase2.UfEqn()
|
||||
==
|
||||
*momentumTransfer[phase2.name()]
|
||||
+ fvOptions(alpha2, rho2, U2)
|
||||
);
|
||||
U2Eqn.relax();
|
||||
fvOptions.constrain(U2Eqn);
|
||||
|
||||
@ -1,7 +0,0 @@
|
||||
tmp<surfaceScalarField> tddtPhi1;
|
||||
tmp<surfaceScalarField> tddtPhi2;
|
||||
|
||||
if (faceMomentum)
|
||||
{
|
||||
#include "pUf/DDtU.H"
|
||||
}
|
||||
@ -22,74 +22,47 @@ surfaceScalarField alphaRhof20
|
||||
);
|
||||
|
||||
// Drag coefficient
|
||||
surfaceScalarField Kdf("Kdf", fluid.Kdf());
|
||||
const surfaceScalarField Kdf("Kdf", fluid.Kdf());
|
||||
|
||||
// Virtual-mass coefficient
|
||||
surfaceScalarField Vmf("Vmf", fluid.Vmf());
|
||||
// Diagonal coefficients
|
||||
PtrList<surfaceScalarField> AFfs(fluid.AFfs());
|
||||
|
||||
surfaceScalarField rAUf1
|
||||
PtrList<surfaceScalarField> rAUfs;
|
||||
rAUfs.append
|
||||
(
|
||||
IOobject::groupName("rAUf", phase1.name()),
|
||||
1.0
|
||||
/(
|
||||
byDt(alphaRhof10 + Vmf)
|
||||
+ fvc::interpolate(U1Eqn.A())
|
||||
+ Kdf
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("rAUf", phase1.name()),
|
||||
1.0
|
||||
/(
|
||||
byDt(alphaRhof10)
|
||||
+ fvc::interpolate(U1Eqn.A())
|
||||
+ AFfs[0]
|
||||
)
|
||||
)
|
||||
);
|
||||
|
||||
surfaceScalarField rAUf2
|
||||
rAUfs.append
|
||||
(
|
||||
IOobject::groupName("rAUf", phase2.name()),
|
||||
1.0
|
||||
/(
|
||||
byDt(alphaRhof20 + Vmf)
|
||||
+ fvc::interpolate(U2Eqn.A())
|
||||
+ Kdf
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("rAUf", phase2.name()),
|
||||
1.0
|
||||
/(
|
||||
byDt(alphaRhof20)
|
||||
+ fvc::interpolate(U2Eqn.A())
|
||||
+ AFfs[0]
|
||||
)
|
||||
)
|
||||
);
|
||||
const surfaceScalarField& rAUf1 = rAUfs[0];
|
||||
const surfaceScalarField& rAUf2 = rAUfs[1];
|
||||
|
||||
AFfs.clear();
|
||||
|
||||
// Turbulent dispersion, particle-pressure, lift and wall-lubrication forces
|
||||
tmp<surfaceScalarField> Ff1;
|
||||
tmp<surfaceScalarField> Ff2;
|
||||
{
|
||||
// Turbulent-dispersion diffusivity
|
||||
volScalarField D(fluid.D());
|
||||
|
||||
// Phase-1 turbulent dispersion and particle-pressure diffusivity
|
||||
surfaceScalarField Df1
|
||||
(
|
||||
fvc::interpolate(D + phase1.pPrime())
|
||||
);
|
||||
|
||||
// Phase-2 turbulent dispersion and particle-pressure diffusivity
|
||||
surfaceScalarField Df2
|
||||
(
|
||||
fvc::interpolate(D + phase2.pPrime())
|
||||
);
|
||||
|
||||
// Cache the phase diffusivities for implicit treatment in the
|
||||
// phase-fraction equation
|
||||
if (implicitPhasePressure)
|
||||
{
|
||||
phase1.DbyA(rAUf1*Df1);
|
||||
phase2.DbyA(rAUf2*Df2);
|
||||
}
|
||||
|
||||
// Lift and wall-lubrication forces
|
||||
surfaceScalarField Ff(fluid.Ff());
|
||||
|
||||
// Phase-fraction face-gradient
|
||||
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
|
||||
|
||||
// Phase-1 dispersion, lift and wall-lubrication force
|
||||
Ff1 = Df1*snGradAlpha1 + Ff;
|
||||
|
||||
// Phase-2 dispersion, lift and wall-lubrication force
|
||||
Ff2 = -Df2*snGradAlpha1 - Ff;
|
||||
}
|
||||
|
||||
// Explicit force fluxes
|
||||
PtrList<surfaceScalarField> phiFfs(fluid.phiFfs(rAUfs));
|
||||
const surfaceScalarField& phiFf1 = phiFfs[0];
|
||||
const surfaceScalarField& phiFf2 = phiFfs[1];
|
||||
|
||||
while (pimple.correct())
|
||||
{
|
||||
@ -105,48 +78,46 @@ while (pimple.correct())
|
||||
MRF.correctBoundaryFlux(U1, phi1);
|
||||
MRF.correctBoundaryFlux(U2, phi2);
|
||||
|
||||
surfaceScalarField alpharAUf1
|
||||
const surfaceScalarField alpharAUf1
|
||||
(
|
||||
IOobject::groupName("alpharAUf", phase1.name()),
|
||||
max(alphaf1, phase1.residualAlpha())*rAUf1
|
||||
);
|
||||
|
||||
surfaceScalarField alpharAUf2
|
||||
const surfaceScalarField alpharAUf2
|
||||
(
|
||||
IOobject::groupName("alpharAUf", phase2.name()),
|
||||
max(alphaf2, phase2.residualAlpha())*rAUf2
|
||||
);
|
||||
|
||||
surfaceScalarField ghSnGradRho
|
||||
// Combined buoyancy and force fluxes
|
||||
const surfaceScalarField ghSnGradRho
|
||||
(
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
|
||||
// Phase-1 buoyancy flux
|
||||
surfaceScalarField phig1
|
||||
const surfaceScalarField phigF1
|
||||
(
|
||||
IOobject::groupName("phig", phase1.name()),
|
||||
alpharAUf1
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf2*(rhof1 - rhof2)*(g & mesh.Sf())
|
||||
)
|
||||
+ phiFf1
|
||||
);
|
||||
|
||||
// Phase-2 buoyancy flux
|
||||
surfaceScalarField phig2
|
||||
const surfaceScalarField phigF2
|
||||
(
|
||||
IOobject::groupName("phig", phase2.name()),
|
||||
alpharAUf2
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf1*(rhof2 - rhof1)*(g & mesh.Sf())
|
||||
)
|
||||
+ phiFf2
|
||||
);
|
||||
|
||||
|
||||
// Phase-1 predicted flux
|
||||
// Predicted fluxes
|
||||
surfaceScalarField phiHbyA1
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase1.name()),
|
||||
@ -156,15 +127,11 @@ while (pimple.correct())
|
||||
phiHbyA1 =
|
||||
rAUf1
|
||||
*(
|
||||
(alphaRhof10 + Vmf)
|
||||
*byDt(MRF.absolute(phi1.oldTime()))
|
||||
alphaRhof10*byDt(MRF.absolute(phi1.oldTime()))
|
||||
+ fvc::flux(U1Eqn.H())
|
||||
+ Vmf*tddtPhi2()
|
||||
+ Kdf*MRF.absolute(phi2)
|
||||
- Ff1()
|
||||
);
|
||||
)
|
||||
- phigF1;
|
||||
|
||||
// Phase-2 predicted flux
|
||||
surfaceScalarField phiHbyA2
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase2.name()),
|
||||
@ -174,26 +141,29 @@ while (pimple.correct())
|
||||
phiHbyA2 =
|
||||
rAUf2
|
||||
*(
|
||||
(alphaRhof20 + Vmf)
|
||||
*byDt(MRF.absolute(phi2.oldTime()))
|
||||
alphaRhof20*byDt(MRF.absolute(phi2.oldTime()))
|
||||
+ fvc::flux(U2Eqn.H())
|
||||
+ Vmf*tddtPhi1()
|
||||
+ Kdf*MRF.absolute(phi1)
|
||||
- Ff2()
|
||||
);
|
||||
)
|
||||
- phigF2;
|
||||
|
||||
// Drag fluxes
|
||||
PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
|
||||
const surfaceScalarField& phiKdPhif1 = phiKdPhifs[0];
|
||||
const surfaceScalarField& phiKdPhif2 = phiKdPhifs[1];
|
||||
|
||||
// Total predicted flux
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
"phiHbyA",
|
||||
alphaf1*(phiHbyA1 - phig1) + alphaf2*(phiHbyA2 - phig2)
|
||||
alphaf1*(phiHbyA1 - phiKdPhif1) + alphaf2*(phiHbyA2 - phiKdPhif2)
|
||||
);
|
||||
|
||||
MRF.makeRelative(phiHbyA);
|
||||
|
||||
phiHbyA1 -= phig1;
|
||||
phiHbyA2 -= phig2;
|
||||
phiKdPhifs.clear();
|
||||
|
||||
surfaceScalarField rAUf
|
||||
// Construct pressure "diffusivity"
|
||||
const surfaceScalarField rAUf
|
||||
(
|
||||
"rAUf",
|
||||
mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2)
|
||||
@ -212,11 +182,8 @@ while (pimple.correct())
|
||||
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
|
||||
);
|
||||
|
||||
tmp<fvScalarMatrix> pEqnComp1;
|
||||
tmp<fvScalarMatrix> pEqnComp2;
|
||||
|
||||
// Construct the compressibility parts of the pressure equation
|
||||
|
||||
tmp<fvScalarMatrix> pEqnComp1, pEqnComp2;
|
||||
if (phase1.compressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
@ -357,18 +324,14 @@ while (pimple.correct())
|
||||
|
||||
phi = phiHbyA + pEqnIncomp.flux();
|
||||
|
||||
surfaceScalarField phi1s
|
||||
const surfaceScalarField phi1s
|
||||
(
|
||||
phiHbyA1
|
||||
+ alpharAUf1*mSfGradp
|
||||
- rAUf1*Kdf*MRF.absolute(phi2)
|
||||
phiHbyA1 + alpharAUf1*mSfGradp
|
||||
);
|
||||
|
||||
surfaceScalarField phi2s
|
||||
const surfaceScalarField phi2s
|
||||
(
|
||||
phiHbyA2
|
||||
+ alpharAUf2*mSfGradp
|
||||
- rAUf2*Kdf*MRF.absolute(phi1)
|
||||
phiHbyA2 + alpharAUf2*mSfGradp
|
||||
);
|
||||
|
||||
surfaceScalarField phir
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -65,16 +65,7 @@ int main(int argc, char *argv[])
|
||||
pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
|
||||
);
|
||||
|
||||
Switch implicitPhasePressure
|
||||
(
|
||||
mesh.solverDict(alpha1.name()).lookupOrDefault<Switch>
|
||||
(
|
||||
"implicitPhasePressure", false
|
||||
)
|
||||
);
|
||||
|
||||
#include "pUf/createRDeltaTf.H"
|
||||
#include "pUf/createDDtU.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -119,7 +110,6 @@ int main(int argc, char *argv[])
|
||||
#include "pUf/UEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pUf/pEqn.H"
|
||||
#include "pUf/DDtU.H"
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
@ -410,7 +410,7 @@ void Foam::RASModels::kineticTheoryModel::correct()
|
||||
// Drag
|
||||
volScalarField beta
|
||||
(
|
||||
refCast<const twoPhaseSystem>(phase_.fluid()).drag(phase_).K()
|
||||
refCast<const twoPhaseSystem>(phase_.fluid()).Kd()
|
||||
);
|
||||
|
||||
// Eq. 3.25, p. 50 Js = J1 - J2
|
||||
|
||||
@ -85,12 +85,6 @@ Foam::twoPhaseSystem::sigma() const
|
||||
}
|
||||
|
||||
|
||||
const Foam::dragModel& Foam::twoPhaseSystem::drag(const phaseModel& phase) const
|
||||
{
|
||||
return lookupSubModel<dragModel>(phase, otherPhase(phase));
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::volScalarField>
|
||||
Foam::twoPhaseSystem::Kd() const
|
||||
{
|
||||
@ -111,13 +105,6 @@ Foam::twoPhaseSystem::Kdf() const
|
||||
}
|
||||
|
||||
|
||||
const Foam::virtualMassModel&
|
||||
Foam::twoPhaseSystem::virtualMass(const phaseModel& phase) const
|
||||
{
|
||||
return lookupSubModel<virtualMassModel>(phase, otherPhase(phase));
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::volScalarField>
|
||||
Foam::twoPhaseSystem::Vm() const
|
||||
{
|
||||
@ -128,52 +115,6 @@ Foam::twoPhaseSystem::Vm() const
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::surfaceScalarField>
|
||||
Foam::twoPhaseSystem::Vmf() const
|
||||
{
|
||||
return Vmf
|
||||
(
|
||||
phasePairKey(phase1().name(), phase2().name())
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::volVectorField>
|
||||
Foam::twoPhaseSystem::F() const
|
||||
{
|
||||
return F
|
||||
(
|
||||
phasePairKey(phase1().name(), phase2().name())
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::surfaceScalarField>
|
||||
Foam::twoPhaseSystem::Ff() const
|
||||
{
|
||||
return Ff
|
||||
(
|
||||
phasePairKey(phase1().name(), phase2().name())
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::volScalarField>
|
||||
Foam::twoPhaseSystem::D() const
|
||||
{
|
||||
return D
|
||||
(
|
||||
phasePairKey(phase1().name(), phase2().name())
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
bool Foam::twoPhaseSystem::transfersMass() const
|
||||
{
|
||||
return transfersMass(phase1());
|
||||
}
|
||||
|
||||
|
||||
void Foam::twoPhaseSystem::solve()
|
||||
{
|
||||
const Time& runTime = mesh_.time();
|
||||
@ -229,10 +170,12 @@ void Foam::twoPhaseSystem::solve()
|
||||
surfaceScalarField phir("phir", phi1 - phi2);
|
||||
|
||||
tmp<surfaceScalarField> alphaDbyA;
|
||||
|
||||
if (notNull(phase1_.DbyA()) && notNull(phase2_.DbyA()))
|
||||
if (DByAfs().found(phase1_.name()) && DByAfs().found(phase2_.name()))
|
||||
{
|
||||
surfaceScalarField DbyA(phase1_.DbyA() + phase2_.DbyA());
|
||||
surfaceScalarField DbyA
|
||||
(
|
||||
*DByAfs()[phase1_.name()] + *DByAfs()[phase2_.name()]
|
||||
);
|
||||
|
||||
alphaDbyA =
|
||||
fvc::interpolate(max(alpha1, scalar(0)))
|
||||
|
||||
@ -53,6 +53,8 @@ class twoPhaseSystem
|
||||
:
|
||||
public phaseSystem
|
||||
{
|
||||
private:
|
||||
|
||||
// Private member functions
|
||||
|
||||
//- Return the drag coefficient for phase pair
|
||||
@ -64,23 +66,6 @@ class twoPhaseSystem
|
||||
//- Return the virtual mass coefficient for phase pair
|
||||
virtual tmp<volScalarField> Vm(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the face virtual mass coefficient for phase pair
|
||||
virtual tmp<surfaceScalarField> Vmf(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the combined force (lift + wall-lubrication) for phase pair
|
||||
virtual tmp<volVectorField> F(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the combined face-force (lift + wall-lubrication)
|
||||
// for phase pair
|
||||
virtual tmp<surfaceScalarField> Ff(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return the turbulent diffusivity for phase pair
|
||||
// Multiplies the phase-fraction gradient
|
||||
virtual tmp<volScalarField> D(const phasePairKey& key) const = 0;
|
||||
|
||||
//- Return true if there is mass transfer for phase
|
||||
virtual bool transfersMass(const phaseModel& phase) const = 0;
|
||||
|
||||
//- Return the interfacial mass flow rate for phase pair
|
||||
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
|
||||
|
||||
@ -89,6 +74,11 @@ protected:
|
||||
|
||||
// Protected data
|
||||
|
||||
//- Flag to indicate that returned lists of fields are "complete"; i.e.,
|
||||
// that an absence of force is returned as a constructed list of zeros,
|
||||
// rather than a null pointer
|
||||
static const bool fillFields_ = true;
|
||||
|
||||
//- Phase model 1
|
||||
phaseModel& phase1_;
|
||||
|
||||
@ -135,6 +125,9 @@ public:
|
||||
|
||||
// Member Functions
|
||||
|
||||
using phaseSystem::sigma;
|
||||
using phaseSystem::dmdt;
|
||||
|
||||
//- Constant access phase model 1
|
||||
const phaseModel& phase1() const;
|
||||
|
||||
@ -148,57 +141,89 @@ public:
|
||||
phaseModel& phase2();
|
||||
|
||||
//- Constant access the phase not given as an argument
|
||||
const phaseModel& otherPhase
|
||||
(
|
||||
const phaseModel& phase
|
||||
) const;
|
||||
const phaseModel& otherPhase(const phaseModel& phase) const;
|
||||
|
||||
//- Return the momentum transfer matrices
|
||||
//- 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;
|
||||
|
||||
using phaseSystem::sigma;
|
||||
|
||||
//- Return the surface tension coefficient
|
||||
tmp<volScalarField> sigma() const;
|
||||
|
||||
//- Return the drag model for the given phase
|
||||
const dragModel& drag(const phaseModel& phase) const;
|
||||
|
||||
//- Return the drag coefficient
|
||||
tmp<volScalarField> Kd() const;
|
||||
|
||||
//- Return the face drag coefficient
|
||||
tmp<surfaceScalarField> Kdf() const;
|
||||
|
||||
//- Return the virtual mass model for the given phase
|
||||
const virtualMassModel& virtualMass(const phaseModel& phase) const;
|
||||
|
||||
//- Return the virtual mass coefficient
|
||||
tmp<volScalarField> Vm() const;
|
||||
|
||||
//- Return the face virtual mass coefficient
|
||||
tmp<surfaceScalarField> Vmf() const;
|
||||
|
||||
//- Return the combined force (lift + wall-lubrication)
|
||||
tmp<volVectorField> F() const;
|
||||
|
||||
//- Return the combined face-force (lift + wall-lubrication)
|
||||
tmp<surfaceScalarField> Ff() const;
|
||||
|
||||
//- Return the turbulent diffusivity
|
||||
// Multiplies the phase-fraction gradient
|
||||
tmp<volScalarField> D() const;
|
||||
|
||||
//- Return true if there is mass transfer
|
||||
bool transfersMass() const;
|
||||
|
||||
using phaseSystem::dmdt;
|
||||
virtual bool transfersMass() const = 0;
|
||||
|
||||
//- Solve for the phase fractions
|
||||
virtual void solve();
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -361,12 +361,12 @@ Foam::RASModels::kineticTheoryModel::divDevRhoReff
|
||||
void Foam::RASModels::kineticTheoryModel::correct()
|
||||
{
|
||||
// Local references
|
||||
const twoPhaseSystem& fluid = refCast<const twoPhaseSystem>(phase_.fluid());
|
||||
volScalarField alpha(max(alpha_, scalar(0)));
|
||||
const volScalarField& rho = phase_.rho();
|
||||
const surfaceScalarField& alphaRhoPhi = alphaRhoPhi_;
|
||||
const volVectorField& U = U_;
|
||||
const volVectorField& Uc_ =
|
||||
refCast<const twoPhaseSystem>(phase_.fluid()).otherPhase(phase_).U();
|
||||
const volVectorField& Uc_ = fluid.otherPhase(phase_).U();
|
||||
|
||||
const scalar sqrtPi = sqrt(constant::mathematical::pi);
|
||||
dimensionedScalar ThetaSmall("ThetaSmall", Theta_.dimensions(), 1.0e-6);
|
||||
@ -410,7 +410,11 @@ void Foam::RASModels::kineticTheoryModel::correct()
|
||||
// Drag
|
||||
volScalarField beta
|
||||
(
|
||||
refCast<const twoPhaseSystem>(phase_.fluid()).drag(phase_).K()
|
||||
fluid.lookupSubModel<dragModel>
|
||||
(
|
||||
phase_,
|
||||
fluid.otherPhase(phase_)
|
||||
).K()
|
||||
);
|
||||
|
||||
// Eq. 3.25, p. 50 Js = J1 - J2
|
||||
|
||||
@ -568,19 +568,6 @@ bool Foam::twoPhaseSystem::read()
|
||||
}
|
||||
|
||||
|
||||
const Foam::dragModel& Foam::twoPhaseSystem::drag(const phaseModel& phase) const
|
||||
{
|
||||
return drag_->phaseModel(phase);
|
||||
}
|
||||
|
||||
|
||||
const Foam::virtualMassModel&
|
||||
Foam::twoPhaseSystem::virtualMass(const phaseModel& phase) const
|
||||
{
|
||||
return virtualMass_->phaseModel(phase);
|
||||
}
|
||||
|
||||
|
||||
const Foam::dimensionedScalar& Foam::twoPhaseSystem::sigma() const
|
||||
{
|
||||
return pair_->sigma();
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -182,11 +182,17 @@ public:
|
||||
|
||||
// Access
|
||||
|
||||
//- Return the drag model for the given phase
|
||||
const dragModel& drag(const phaseModel& phase) const;
|
||||
//- Access a sub model between a phase pair
|
||||
template<class modelType>
|
||||
const modelType& lookupSubModel(const phasePair& key) const;
|
||||
|
||||
//- Return the virtual mass model for the given phase
|
||||
const virtualMassModel& virtualMass(const phaseModel& phase) const;
|
||||
//- Access a sub model between two phases
|
||||
template<class modelType>
|
||||
const modelType& lookupSubModel
|
||||
(
|
||||
const phaseModel& dispersed,
|
||||
const phaseModel& continuous
|
||||
) const;
|
||||
|
||||
//- Return the surface tension coefficient
|
||||
const dimensionedScalar& sigma() const;
|
||||
@ -236,6 +242,12 @@ public:
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#ifdef NoRepository
|
||||
#include "twoPhaseSystemTemplates.C"
|
||||
#endif
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -0,0 +1,77 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / 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 "BlendedInterfacialModel.H"
|
||||
#include "dragModel.H"
|
||||
#include "virtualMassModel.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
template<class modelType>
|
||||
const modelType& twoPhaseSystem::lookupSubModel
|
||||
(
|
||||
const phasePair& key
|
||||
) const
|
||||
{
|
||||
return
|
||||
mesh().lookupObject<modelType>
|
||||
(
|
||||
IOobject::groupName(modelType::typeName, key.name())
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
template<>
|
||||
inline const dragModel& twoPhaseSystem::lookupSubModel<dragModel>
|
||||
(
|
||||
const phaseModel& dispersed,
|
||||
const phaseModel& continuous
|
||||
) const
|
||||
{
|
||||
return drag_->phaseModel(dispersed);
|
||||
}
|
||||
|
||||
|
||||
template<>
|
||||
inline const virtualMassModel& twoPhaseSystem::lookupSubModel<virtualMassModel>
|
||||
(
|
||||
const phaseModel& dispersed,
|
||||
const phaseModel& continuous
|
||||
) const
|
||||
{
|
||||
return virtualMass_->phaseModel(dispersed);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -26,6 +26,7 @@ License
|
||||
#include "NicenoKEqn.H"
|
||||
#include "fvOptions.H"
|
||||
#include "twoPhaseSystem.H"
|
||||
#include "dragModel.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -181,11 +182,13 @@ tmp<volScalarField> NicenoKEqn<BasicTurbulenceModel>::bubbleG() const
|
||||
refCast<const twoPhaseSystem>(liquid.fluid());
|
||||
const transportModel& gas = fluid.otherPhase(liquid);
|
||||
|
||||
const dragModel& drag = fluid.lookupSubModel<dragModel>(gas, liquid);
|
||||
|
||||
volScalarField magUr(mag(this->U_ - gasTurbulence.U()));
|
||||
|
||||
tmp<volScalarField> bubbleG
|
||||
(
|
||||
Cp_*sqr(magUr)*fluid.drag(gas).K()/liquid.rho()
|
||||
Cp_*sqr(magUr)*drag.K()/liquid.rho()
|
||||
);
|
||||
|
||||
return bubbleG;
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -26,6 +26,7 @@ License
|
||||
#include "LaheyKEpsilon.H"
|
||||
#include "fvOptions.H"
|
||||
#include "twoPhaseSystem.H"
|
||||
#include "dragModel.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -191,6 +192,8 @@ tmp<volScalarField> LaheyKEpsilon<BasicTurbulenceModel>::bubbleG() const
|
||||
const twoPhaseSystem& fluid = refCast<const twoPhaseSystem>(liquid.fluid());
|
||||
const transportModel& gas = fluid.otherPhase(liquid);
|
||||
|
||||
const dragModel& drag = fluid.lookupSubModel<dragModel>(gas, liquid);
|
||||
|
||||
volScalarField magUr(mag(this->U_ - gasTurbulence.U()));
|
||||
|
||||
tmp<volScalarField> bubbleG
|
||||
@ -198,7 +201,7 @@ tmp<volScalarField> LaheyKEpsilon<BasicTurbulenceModel>::bubbleG() const
|
||||
Cp_
|
||||
*(
|
||||
pow3(magUr)
|
||||
+ pow(fluid.drag(gas).CdRe()*liquid.nu()/gas.d(), 4.0/3.0)
|
||||
+ pow(drag.CdRe()*liquid.nu()/gas.d(), 4.0/3.0)
|
||||
*pow(magUr, 5.0/3.0)
|
||||
)
|
||||
*gas
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -26,6 +26,7 @@ License
|
||||
#include "continuousGasKEpsilon.H"
|
||||
#include "fvOptions.H"
|
||||
#include "twoPhaseSystem.H"
|
||||
#include "dragModel.H"
|
||||
#include "virtualMassModel.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
@ -122,8 +123,11 @@ void continuousGasKEpsilon<BasicTurbulenceModel>::correctNut()
|
||||
const twoPhaseSystem& fluid = refCast<const twoPhaseSystem>(gas.fluid());
|
||||
const transportModel& liquid = fluid.otherPhase(gas);
|
||||
|
||||
const virtualMassModel& virtualMass =
|
||||
fluid.lookupSubModel<virtualMassModel>(gas, liquid);
|
||||
|
||||
volScalarField thetal(liquidTurbulence.k()/liquidTurbulence.epsilon());
|
||||
volScalarField rhodv(gas.rho() + fluid.virtualMass(gas).Cvm()*liquid.rho());
|
||||
volScalarField rhodv(gas.rho() + virtualMass.Cvm()*liquid.rho());
|
||||
volScalarField thetag((rhodv/(18*liquid.rho()*liquid.nu()))*sqr(gas.d()));
|
||||
volScalarField expThetar
|
||||
(
|
||||
@ -206,12 +210,15 @@ continuousGasKEpsilon<BasicTurbulenceModel>::rhoEff() const
|
||||
const twoPhaseSystem& fluid = refCast<const twoPhaseSystem>(gas.fluid());
|
||||
const transportModel& liquid = fluid.otherPhase(gas);
|
||||
|
||||
const virtualMassModel& virtualMass =
|
||||
fluid.lookupSubModel<virtualMassModel>(gas, liquid);
|
||||
|
||||
return tmp<volScalarField>
|
||||
(
|
||||
new volScalarField
|
||||
(
|
||||
IOobject::groupName("rhoEff", this->alphaRhoPhi_.group()),
|
||||
gas.rho() + (fluid.virtualMass(gas).Cvm() + 3.0/20.0)*liquid.rho()
|
||||
gas.rho() + (virtualMass.Cvm() + 3.0/20.0)*liquid.rho()
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -27,6 +27,7 @@ License
|
||||
#include "fvOptions.H"
|
||||
#include "bound.H"
|
||||
#include "twoPhaseSystem.H"
|
||||
#include "dragModel.H"
|
||||
#include "virtualMassModel.H"
|
||||
#include "fixedValueFvPatchFields.H"
|
||||
#include "inletOutletFvPatchFields.H"
|
||||
@ -389,7 +390,7 @@ tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::Ct2() const
|
||||
volScalarField beta
|
||||
(
|
||||
(6*this->Cmu_/(4*sqrt(3.0/2.0)))
|
||||
*fluid.drag(gas).K()/liquid.rho()
|
||||
*fluid.Kd()/liquid.rho()
|
||||
*(liquidTurbulence.k_/liquidTurbulence.epsilon_)
|
||||
);
|
||||
volScalarField Ct0((3 + beta)/(1 + beta + 2*gas.rho()/liquid.rho()));
|
||||
@ -413,9 +414,9 @@ tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rhogEff() const
|
||||
{
|
||||
const transportModel& gas = this->transport();
|
||||
const twoPhaseSystem& fluid = refCast<const twoPhaseSystem>(gas.fluid());
|
||||
return
|
||||
gas.rho()
|
||||
+ fluid.virtualMass(gas).Cvm()*fluid.otherPhase(gas).rho();
|
||||
const virtualMassModel& virtualMass =
|
||||
fluid.lookupSubModel<virtualMassModel>(gas, fluid.otherPhase(gas));
|
||||
return gas.rho() + virtualMass.Cvm()*fluid.otherPhase(gas).rho();
|
||||
}
|
||||
|
||||
|
||||
@ -491,6 +492,8 @@ tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::bubbleG() const
|
||||
const twoPhaseSystem& fluid = refCast<const twoPhaseSystem>(gas.fluid());
|
||||
const transportModel& liquid = fluid.otherPhase(gas);
|
||||
|
||||
const dragModel& drag = fluid.lookupSubModel<dragModel>(gas, liquid);
|
||||
|
||||
volScalarField magUr(mag(liquidTurbulence.U() - this->U()));
|
||||
|
||||
// Lahey model
|
||||
@ -500,7 +503,7 @@ tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::bubbleG() const
|
||||
*liquid*liquid.rho()
|
||||
*(
|
||||
pow3(magUr)
|
||||
+ pow(fluid.drag(gas).CdRe()*liquid.nu()/gas.d(), 4.0/3.0)
|
||||
+ pow(drag.CdRe()*liquid.nu()/gas.d(), 4.0/3.0)
|
||||
*pow(magUr, 5.0/3.0)
|
||||
)
|
||||
*gas
|
||||
@ -510,7 +513,7 @@ tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::bubbleG() const
|
||||
// Simple model
|
||||
// tmp<volScalarField> bubbleG
|
||||
// (
|
||||
// Cp_*liquid*fluid.drag(gas).K()*sqr(magUr)
|
||||
// Cp_*liquid*drag.K()*sqr(magUr)
|
||||
// );
|
||||
|
||||
return bubbleG;
|
||||
|
||||
@ -17,7 +17,7 @@ FoamFile
|
||||
|
||||
solvers
|
||||
{
|
||||
alpha.air
|
||||
"alpha.*"
|
||||
{
|
||||
nAlphaCorr 1;
|
||||
nAlphaSubCycles 2;
|
||||
|
||||
@ -17,7 +17,7 @@ FoamFile
|
||||
|
||||
solvers
|
||||
{
|
||||
alpha.air
|
||||
"alpha.*"
|
||||
{
|
||||
nAlphaCorr 1;
|
||||
nAlphaSubCycles 2;
|
||||
|
||||
@ -17,7 +17,7 @@ FoamFile
|
||||
|
||||
solvers
|
||||
{
|
||||
alpha.air
|
||||
"alpha.*"
|
||||
{
|
||||
nAlphaCorr 1;
|
||||
nAlphaSubCycles 2;
|
||||
|
||||
@ -17,7 +17,7 @@ FoamFile
|
||||
|
||||
solvers
|
||||
{
|
||||
alpha.air
|
||||
"alpha.*"
|
||||
{
|
||||
nAlphaCorr 1;
|
||||
nAlphaSubCycles 2;
|
||||
|
||||
@ -17,7 +17,7 @@ FoamFile
|
||||
|
||||
solvers
|
||||
{
|
||||
alpha.air
|
||||
"alpha.*"
|
||||
{
|
||||
nAlphaCorr 1;
|
||||
nAlphaSubCycles 2;
|
||||
|
||||
@ -17,7 +17,7 @@ FoamFile
|
||||
|
||||
solvers
|
||||
{
|
||||
alpha.air
|
||||
"alpha.*"
|
||||
{
|
||||
nAlphaCorr 1;
|
||||
nAlphaSubCycles 2;
|
||||
|
||||
@ -17,7 +17,7 @@ FoamFile
|
||||
|
||||
solvers
|
||||
{
|
||||
alpha.air
|
||||
"alpha.*"
|
||||
{
|
||||
nAlphaCorr 1;
|
||||
nAlphaSubCycles 2;
|
||||
|
||||
Reference in New Issue
Block a user