From e352828514abe029065d8032ec78516c980a7a2d Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 22 Mar 2018 11:36:43 +0000 Subject: [PATCH] reactingMultiphaseEulerFoam: Stationary phase Two new phase models have been added as selectable options for reactingMultiphaseEulerFoam; pureStationaryPhaseModel and pureStationaryIsothermalPhaseModel. These phases do not store a velocity and their phase fractions remain constant throughout the simulation. They are intended for use in modelling static particle beds and other forms of porous media by means of the existing Euler-Euler transfer models (drag, heat transfer, etc...). Note that this functionality has not been extended to reactingTwoPhaseEulerFoam, or the non-reacting *EulerFoam solvers. Additional maintenance work has been carried out on the phase model and phase system structure. The system can now loop over subsets of phases with specific functionality (moving, multi-component, etc...) in order to avoid testing for the existence of equations or variables in the top level solver. The mass transfer handling and it's effect on per-phase source terms has been refactored to reduce duplication. Const and non-const access to phase properties has been formalised by renaming non-const accessors with a "Ref" suffix, which is consistent with other recent developments to classes including tmp and GeometricField, among others. More sub-modelling details have been made private in order to reduce the size of interfaces and improve abstraction. This work was supported by Zhen Li, at Evonik --- .../HeatAndMassTransferPhaseSystem.C | 81 ---- .../HeatAndMassTransferPhaseSystem.H | 17 - .../HeatTransferPhaseSystem.C | 61 --- .../HeatTransferPhaseSystem.H | 9 - ...terfaceCompositionPhaseChangePhaseSystem.C | 172 +------ ...terfaceCompositionPhaseChangePhaseSystem.H | 42 +- .../MomentumTransferPhaseSystem.C | 439 +++++++++--------- .../MomentumTransferPhaseSystem.H | 47 +- .../PopulationBalancePhaseSystem.C | 153 ++---- .../PopulationBalancePhaseSystem.H | 38 +- .../ThermalPhaseChangePhaseSystem.C | 260 +++-------- .../ThermalPhaseChangePhaseSystem.H | 63 +-- .../AnisothermalPhaseModel.C | 121 ++--- .../AnisothermalPhaseModel.H | 22 +- .../InertPhaseModel/InertPhaseModel.H | 9 +- .../IsothermalPhaseModel.C | 18 +- .../IsothermalPhaseModel.H | 10 +- .../MovingPhaseModel/MovingPhaseModel.C | 168 ++++--- .../MovingPhaseModel/MovingPhaseModel.H | 98 ++-- .../MultiComponentPhaseModel.C | 72 +-- .../MultiComponentPhaseModel.H | 27 +- .../PurePhaseModel/PurePhaseModel.C | 50 +- .../PurePhaseModel/PurePhaseModel.H | 20 +- .../ReactingPhaseModel/ReactingPhaseModel.H | 9 +- .../StationaryPhaseModel.C | 367 +++++++++++++++ .../StationaryPhaseModel.H | 226 +++++++++ .../ThermoPhaseModel/ThermoPhaseModel.C | 11 +- .../ThermoPhaseModel/ThermoPhaseModel.H | 39 +- .../phaseModel/phaseModel/phaseModel.C | 29 -- .../phaseModel/phaseModel/phaseModel.H | 113 +++-- .../phaseModel/phaseModel/phaseModels.C | 51 +- .../phaseSystems/phaseSystem/phaseSystem.C | 71 ++- .../phaseSystems/phaseSystem/phaseSystem.H | 350 ++++++++------ .../phaseSystems/phaseSystem/phaseSystemI.H | 58 ++- .../phaseSystem/phaseSystemTemplates.C | 2 + .../driftModels/constantDrift/constantDrift.C | 4 +- .../constantNucleation/constantNucleation.C | 6 +- .../wallBoiling/wallBoiling.C | 4 +- .../reactingMultiphaseEulerFoam/EEqns.H | 33 +- .../reactingMultiphaseEulerFoam/YEqns.H | 29 +- .../createFields.H | 2 +- .../multiphaseSystem/multiphaseSystem.C | 156 +++++-- .../multiphaseSystem/multiphaseSystem.H | 75 +-- .../reactingMultiphaseEulerFoam/pU/UEqns.H | 12 +- .../reactingMultiphaseEulerFoam/pU/pEqn.H | 137 ++++-- .../reactingMultiphaseEulerFoam/pUf/UEqns.H | 14 +- .../reactingMultiphaseEulerFoam/pUf/pEqn.H | 169 ++++--- .../reactingTwoPhaseEulerFoam/EEqns.H | 60 ++- .../reactingTwoPhaseEulerFoam/YEqns.H | 31 +- .../createFieldRefs.H | 21 +- .../reactingTwoPhaseEulerFoam/createFields.H | 2 +- .../reactingTwoPhaseEulerFoam/pU/pEqn.H | 75 ++- .../reactingTwoPhaseEulerFoam/pUf/pEqn.H | 75 ++- .../twoPhaseSystem/twoPhaseSystem.C | 19 +- .../twoPhaseSystem/twoPhaseSystem.H | 81 +--- .../laminar/bed/0/T.air | 36 ++ .../laminar/bed/0/T.solid | 33 ++ .../laminar/bed/0/T.water | 36 ++ .../laminar/bed/0/U.air | 35 ++ .../laminar/bed/0/U.water | 35 ++ .../laminar/bed/0/alpha.air.orig | 36 ++ .../laminar/bed/0/alpha.solid.orig | 33 ++ .../laminar/bed/0/alpha.water.orig | 36 ++ .../laminar/bed/0/p | 35 ++ .../laminar/bed/0/p_rgh | 39 ++ .../laminar/bed/Allrun | 11 + .../laminar/bed/constant/fvOptions | 46 ++ .../laminar/bed/constant/g | 22 + .../laminar/bed/constant/phaseProperties | 208 +++++++++ .../bed/constant/thermophysicalProperties.air | 48 ++ .../constant/thermophysicalProperties.solid | 52 +++ .../constant/thermophysicalProperties.water | 53 +++ .../bed/constant/turbulenceProperties.air | 20 + .../bed/constant/turbulenceProperties.water | 20 + .../laminar/bed/system/blockMeshDict | 58 +++ .../laminar/bed/system/controlDict | 54 +++ .../laminar/bed/system/fvSchemes | 61 +++ .../laminar/bed/system/fvSolution | 75 +++ .../laminar/bed/system/setFieldsDict | 40 ++ 79 files changed, 3495 insertions(+), 1955 deletions(-) create mode 100644 applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C create mode 100644 applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.air create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.solid create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.water create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/U.air create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/U.water create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.air.orig create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.solid.orig create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.water.orig create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/p create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/p_rgh create mode 100755 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/Allrun create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/fvOptions create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/g create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/phaseProperties create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.air create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.solid create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.water create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/turbulenceProperties.air create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/turbulenceProperties.water create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/blockMeshDict create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/controlDict create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/fvSchemes create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/fvSolution create mode 100644 tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/setFieldsDict diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C index 015afd01c..f6328ddaa 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C @@ -131,86 +131,6 @@ Foam::HeatAndMassTransferPhaseSystem:: // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -template -bool -Foam::HeatAndMassTransferPhaseSystem::transfersMass() const -{ - return true; -} - - -template -Foam::tmp -Foam::HeatAndMassTransferPhaseSystem::dmdt -( - const phasePairKey& key -) const -{ - return tmp - ( - new volScalarField - ( - IOobject - ( - "dmdt", - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); -} - - -template -Foam::tmp -Foam::HeatAndMassTransferPhaseSystem::dmdt -( - const Foam::phaseModel& phase -) const -{ - tmp tDmdt - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("dmdt", phase.name()), - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); - - return tDmdt; -} - - -template -Foam::autoPtr -Foam::HeatAndMassTransferPhaseSystem::momentumTransfer() const -{ - autoPtr - eqnsPtr(BasePhaseSystem::momentumTransfer()); - - return eqnsPtr; -} - - -template -Foam::autoPtr -Foam::HeatAndMassTransferPhaseSystem::momentumTransferf() const -{ - autoPtr - eqnsPtr(BasePhaseSystem::momentumTransferf()); - - return eqnsPtr; -} - - template Foam::autoPtr Foam::HeatAndMassTransferPhaseSystem::heatTransfer() const @@ -328,7 +248,6 @@ Foam::HeatAndMassTransferPhaseSystem::massTransfer() const template void Foam::HeatAndMassTransferPhaseSystem::correctThermo() { - phaseSystem::correctThermo(); forAllConstIter diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H index 46a16022e..5917ce976 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H @@ -113,23 +113,6 @@ public: // Member Functions - //- Return true if there is mass transfer - virtual bool transfersMass() const; - - //- Return the interfacial mass flow rate - virtual tmp dmdt(const phasePairKey& key) const; - - //- Return the total interfacial mass transfer rate for phase - virtual tmp dmdt(const phaseModel& phase) const; - - //- Return the momentum transfer matrices - virtual autoPtr - momentumTransfer() const; - - //- Return the momentum transfer matrices for the face based algorithm - virtual autoPtr - momentumTransferf() const; - //- Return the heat transfer matrices virtual autoPtr heatTransfer() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C index 6613c0c65..c80c29079 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C @@ -49,67 +49,6 @@ Foam::HeatTransferPhaseSystem::~HeatTransferPhaseSystem() // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -template -bool Foam::HeatTransferPhaseSystem::transfersMass() const -{ - return false; -} - - -template -Foam::tmp -Foam::HeatTransferPhaseSystem::dmdt -( - const phasePairKey& key -) const -{ - return tmp - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName - ( - "dmdt", - this->phasePairs_[key]->name() - ), - this->mesh().time().timeName(), - this->mesh().time() - ), - this->mesh(), - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); -} - - -template -Foam::tmp -Foam::HeatTransferPhaseSystem::dmdt -( - const Foam::phaseModel& phase -) const -{ - tmp tDmdt - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("dmdt", phase.name()), - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); - - return tDmdt; -} - - template Foam::autoPtr Foam::HeatTransferPhaseSystem::heatTransfer() const diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H index da67063c3..664d16f15 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H @@ -89,15 +89,6 @@ public: // Member Functions - //- Return true if there is mass transfer - virtual bool transfersMass() const; - - //- Return the interfacial mass flow rate - virtual tmp dmdt(const phasePairKey& key) const; - - //- Return the total interfacial mass transfer rate for phase - virtual tmp dmdt(const phaseModel& phase) const; - //- Return the heat transfer matrices virtual autoPtr heatTransfer() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C index 98b54af68..2c4b99511 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C @@ -30,37 +30,24 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // template -void Foam::InterfaceCompositionPhaseChangePhaseSystem:: -addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const +Foam::tmp +Foam::InterfaceCompositionPhaseChangePhaseSystem::iDmdt +( + const phasePairKey& key +) const { - // Source term due to mass transfer - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) + if (!iDmdt_.found(key)) { - const phasePair& pair(phasePairIter()); - - if (pair.ordered()) - { - continue; - } - - const volVectorField& U1(pair.phase1().U()); - const volVectorField& U2(pair.phase2().U()); - - const volScalarField dmdt(this->iDmdt(pair)); - const volScalarField dmdt21(posPart(dmdt)); - const volScalarField dmdt12(negPart(dmdt)); - - *eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1); - *eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2); + return phaseSystem::dmdt(key); } + + const scalar dmdtSign(Pair::compare(iDmdt_.find(key).key(), key)); + + return dmdtSign**iDmdt_[key]; } + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template @@ -142,140 +129,31 @@ Foam::InterfaceCompositionPhaseChangePhaseSystem:: template Foam::tmp -Foam::InterfaceCompositionPhaseChangePhaseSystem::iDmdt +Foam::InterfaceCompositionPhaseChangePhaseSystem::dmdt ( const phasePairKey& key ) const { - const scalar dmdtSign(Pair::compare(iDmdt_.find(key).key(), key)); - - return dmdtSign**iDmdt_[key]; + return BasePhaseSystem::dmdt(key) + this->iDmdt(key); } template -Foam::tmp -Foam::InterfaceCompositionPhaseChangePhaseSystem::iDmdt -( - const Foam::phaseModel& phase -) const +Foam::Xfer> +Foam::InterfaceCompositionPhaseChangePhaseSystem::dmdts() const { - tmp tiDmdt - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("iDmdt", phase.name()), - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); + PtrList dmdts(BasePhaseSystem::dmdts()); - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) + forAllConstIter(iDmdtTable, iDmdt_, iDmdtIter) { - const phasePair& pair(phasePairIter()); + const phasePair& pair = this->phasePairs_[iDmdtIter.key()]; + const volScalarField& iDmdt = *iDmdtIter(); - if (pair.ordered()) - { - continue; - } - - if (pair.contains(phase)) - { - tiDmdt.ref() += this->iDmdt - ( - phasePairKey(phase.name(), pair.otherPhase(phase).name(), false) - ); - } + this->addField(pair.phase1(), "dmdt", iDmdt, dmdts); + this->addField(pair.phase2(), "dmdt", - iDmdt, dmdts); } - return tiDmdt; -} - - -template -Foam::tmp -Foam::InterfaceCompositionPhaseChangePhaseSystem::dmdt -( - const Foam::phaseModel& phase -) const -{ - tmp tDmdt - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("dmdt", phase.name()), - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); - - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) - { - const phasePair& pair(phasePairIter()); - - if (pair.ordered()) - { - continue; - } - - if (pair.contains(phase)) - { - tDmdt.ref() += this->dmdt - ( - phasePairKey(phase.name(), pair.otherPhase(phase).name(), false) - ); - } - } - - return tDmdt; -} - - -template -Foam::autoPtr -Foam::InterfaceCompositionPhaseChangePhaseSystem:: -momentumTransfer() const -{ - autoPtr - eqnsPtr(BasePhaseSystem::momentumTransfer()); - - addMomentumTransfer(eqnsPtr()); - - return eqnsPtr; -} - - -template -Foam::autoPtr -Foam::InterfaceCompositionPhaseChangePhaseSystem:: -momentumTransferf() const -{ - autoPtr - eqnsPtr(BasePhaseSystem::momentumTransferf()); - - addMomentumTransfer(eqnsPtr()); - - return eqnsPtr; + return dmdts.xfer(); } @@ -315,8 +193,8 @@ heatTransfer() const const volScalarField& he1(phase1.thermo().he()); const volScalarField& he2(phase2.thermo().he()); - const volScalarField& K1(phase1.K()); - const volScalarField& K2(phase2.K()); + const volScalarField K1(phase1.K()); + const volScalarField K2(phase2.K()); const volScalarField dmdt(this->dmdt(pair)); const volScalarField dmdt21(posPart(dmdt)); diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H index 59ad55674..bc9e9831f 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H @@ -69,6 +69,9 @@ protected: phasePairKey::hash > interfaceCompositionModelTable; + typedef HashPtrTable + iDmdtTable; + // Protected data @@ -78,20 +81,16 @@ protected: interfaceCompositionModelTable interfaceCompositionModels_; //- Interfacial Mass transfer rate - HashPtrTable - iDmdt_; + iDmdtTable iDmdt_; //- Explicit part of the mass transfer rate - HashPtrTable - iDmdtExplicit_; + iDmdtTable iDmdtExplicit_; + // Protected member functions - //- Add the mass-based momentum transfer to the equations - void addMomentumTransfer - ( - phaseSystem::momentumTransferTable& eqns - ) const; + //- Return the interfacial mass transfer rate for a pair for a pair + virtual tmp iDmdt(const phasePairKey& key) const; public: @@ -108,28 +107,11 @@ public: // Member Functions - //- Return the interfacial mass flow rate - virtual tmp iDmdt(const phasePairKey& key) const; + //- Return the mass transfer rate for a pair + virtual tmp dmdt(const phasePairKey& key) const; - //- Return the total interfacial mass transfer rate for phase - virtual tmp iDmdt(const phaseModel& phase) const; - - //- Return the interfacial mass flow rate - virtual tmp dmdt(const phasePairKey& key) const - { - return BasePhaseSystem::dmdt(key) + this->iDmdt(key); - }; - - //- Return the total interfacial mass transfer rate for phase - virtual tmp dmdt(const phaseModel& phase) const; - - //- Return the momentum transfer matrices - virtual autoPtr - momentumTransfer() const; - - //- Return the momentum transfer matrices for the face based algorithm - virtual autoPtr - momentumTransferf() const; + //- Return the mass transfer rates for each phase + virtual Xfer> dmdts() const; //- Return the heat transfer matrices virtual autoPtr heatTransfer() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index 6bc8b0668..cedf2781c 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -44,6 +44,150 @@ License #include "fvcSnGrad.H" #include "fvMatrix.H" + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +template +Foam::tmp +Foam::MomentumTransferPhaseSystem::Kd +( + const phasePairKey& key +) const +{ + if (dragModels_.found(key)) + { + return dragModels_[key]->K(); + } + else + { + return tmp + ( + new volScalarField + ( + IOobject + ( + dragModel::typeName + ":K", + this->mesh_.time().timeName(), + this->mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + this->mesh_, + dimensionedScalar("zero", dragModel::dimK, 0) + ) + ); + } +} + + +template +Foam::tmp +Foam::MomentumTransferPhaseSystem::Kdf +( + const phasePairKey& key +) const +{ + if (dragModels_.found(key)) + { + return dragModels_[key]->Kf(); + } + else + { + return tmp + ( + new surfaceScalarField + ( + IOobject + ( + dragModel::typeName + ":K", + this->mesh_.time().timeName(), + this->mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + this->mesh_, + dimensionedScalar("zero", dragModel::dimK, 0) + ) + ); + } +} + + +template +Foam::tmp +Foam::MomentumTransferPhaseSystem::Vm +( + const phasePairKey& key +) const +{ + if (virtualMassModels_.found(key)) + { + return virtualMassModels_[key]->K(); + } + else + { + return tmp + ( + new volScalarField + ( + IOobject + ( + virtualMassModel::typeName + ":K", + this->mesh_.time().timeName(), + this->mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + this->mesh_, + dimensionedScalar("zero", virtualMassModel::dimK, 0) + ) + ); + } +} + + +template +void Foam::MomentumTransferPhaseSystem:: +addMassTransferMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const +{ + forAllConstIter + ( + phaseSystem::phasePairTable, + this->phasePairs_, + phasePairIter + ) + { + const phasePair& pair(phasePairIter()); + + if (pair.ordered()) + { + continue; + } + + const volScalarField dmdt(this->dmdt(pair)); + + if (!pair.phase1().stationary()) + { + fvVectorMatrix& eqn = *eqns[pair.phase1().name()]; + const volScalarField dmdt21(posPart(dmdt)); + + eqn += dmdt21*pair.phase2().U() - fvm::Sp(dmdt21, eqn.psi()); + } + + if (!pair.phase2().stationary()) + { + fvVectorMatrix& eqn = *eqns[pair.phase2().name()]; + const volScalarField dmdt12(negPart(dmdt)); + + eqn -= dmdt12*pair.phase1().U() - fvm::Sp(dmdt12, eqn.psi()); + } + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template @@ -157,111 +301,9 @@ Foam::MomentumTransferPhaseSystem:: // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -template -Foam::tmp -Foam::MomentumTransferPhaseSystem::Kd -( - const phasePairKey& key -) const -{ - if (dragModels_.found(key)) - { - return dragModels_[key]->K(); - } - else - { - return tmp - ( - new volScalarField - ( - IOobject - ( - dragModel::typeName + ":K", - this->mesh_.time().timeName(), - this->mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - this->mesh_, - dimensionedScalar("zero", dragModel::dimK, 0) - ) - ); - } -} - - -template -Foam::tmp -Foam::MomentumTransferPhaseSystem::Kdf -( - const phasePairKey& key -) const -{ - if (dragModels_.found(key)) - { - return dragModels_[key]->Kf(); - } - else - { - return tmp - ( - new surfaceScalarField - ( - IOobject - ( - dragModel::typeName + ":K", - this->mesh_.time().timeName(), - this->mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - this->mesh_, - dimensionedScalar("zero", dragModel::dimK, 0) - ) - ); - } -} - - -template -Foam::tmp -Foam::MomentumTransferPhaseSystem::Vm -( - const phasePairKey& key -) const -{ - if (virtualMassModels_.found(key)) - { - return virtualMassModels_[key]->K(); - } - else - { - return tmp - ( - new volScalarField - ( - IOobject - ( - virtualMassModel::typeName + ":K", - this->mesh_.time().timeName(), - this->mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - this->mesh_, - dimensionedScalar("zero", virtualMassModel::dimK, 0) - ) - ); - } -} - - template Foam::autoPtr -Foam::MomentumTransferPhaseSystem::momentumTransfer() const +Foam::MomentumTransferPhaseSystem::momentumTransfer() { // Create a momentum transfer matrix for each phase autoPtr eqnsPtr @@ -271,9 +313,9 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const phaseSystem::momentumTransferTable& eqns = eqnsPtr(); - forAll(this->phaseModels_, phasei) + forAll(this->movingPhases(), movingPhasei) { - const phaseModel& phase = this->phaseModels_[phasei]; + const phaseModel& phase = this->movingPhases()[movingPhasei]; eqns.insert ( @@ -295,19 +337,19 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const } // Add the implicit part of the drag force - forAllConstIter - ( - phaseSystem::KdTable, - Kds_, - KdIter - ) + forAllConstIter(KdTable, Kds_, KdIter) { const volScalarField& K(*KdIter()); const phasePair& pair(this->phasePairs_[KdIter.key()]); forAllConstIter(phasePair, pair, iter) { - *eqns[iter().name()] -= fvm::Sp(K, iter().U()()); + if (!iter().stationary()) + { + fvVectorMatrix& eqn = *eqns[iter().name()]; + + eqn -= fvm::Sp(K, eqn.psi()); + } } } @@ -324,12 +366,7 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const } // Add the virtual mass force - forAllConstIter - ( - phaseSystem::VmTable, - Vms_, - VmIter - ) + forAllConstIter(VmTable, Vms_, VmIter) { const volScalarField& Vm(*VmIter()); const phasePair& pair(this->phasePairs_[VmIter.key()]); @@ -339,28 +376,36 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const const phaseModel& phase = iter(); const phaseModel& otherPhase = iter.otherPhase(); - const volVectorField& U = phase.U(); - const surfaceScalarField& phi = phase.phi(); + if (!phase.stationary()) + { + fvVectorMatrix& eqn = *eqns[phase.name()]; - *eqns[phase.name()] -= - Vm - *( - fvm::ddt(U) - + fvm::div(phi, U) - - fvm::Sp(fvc::div(phi), U) - - otherPhase.DUDt() - ) - + this->MRF_.DDt(Vm, U - otherPhase.U()); + const volVectorField& U = eqn.psi(); + const surfaceScalarField& phi = phase.phi(); + + eqn -= + Vm + *( + fvm::ddt(U) + + fvm::div(phi, U) + - fvm::Sp(fvc::div(phi), U) + - otherPhase.DUDt() + ) + + this->MRF_.DDt(Vm, U - otherPhase.U()); + } } } + // Add the source term due to mass trasfer + addMassTransferMomentumTransfer(eqns); + return eqnsPtr; } template Foam::autoPtr -Foam::MomentumTransferPhaseSystem::momentumTransferf() const +Foam::MomentumTransferPhaseSystem::momentumTransferf() { // Create a momentum transfer matrix for each phase autoPtr eqnsPtr @@ -370,9 +415,9 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() const phaseSystem::momentumTransferTable& eqns = eqnsPtr(); - forAll(this->phaseModels_, phasei) + forAll(this->movingPhases(), movingPhasei) { - const phaseModel& phase = this->phaseModels_[phasei]; + const phaseModel& phase = this->movingPhases()[movingPhasei]; eqns.insert ( @@ -401,12 +446,7 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() const } // Add the virtual mass force - forAllConstIter - ( - phaseSystem::VmTable, - Vms_, - VmIter - ) + forAllConstIter(VmTable, Vms_, VmIter) { const volScalarField& Vm(*VmIter()); const phasePair& pair(this->phasePairs_[VmIter.key()]); @@ -416,15 +456,21 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() const const phaseModel& phase = iter(); const phaseModel& otherPhase = iter.otherPhase(); - *eqns[phase.name()] -= - Vm - *( - UgradUs[phase.index()] - - (UgradUs[otherPhase.index()] & otherPhase.U()) - ); + if (!phase.stationary()) + { + *eqns[phase.name()] -= + Vm + *( + UgradUs[phase.index()] + - (UgradUs[otherPhase.index()] & otherPhase.U()) + ); + } } } + // Add the source term due to mass trasfer + addMassTransferMomentumTransfer(eqns); + return eqnsPtr; } @@ -436,12 +482,7 @@ Foam::MomentumTransferPhaseSystem::AFfs() const PtrList AFfs(this->phaseModels_.size()); // Add the implicit part of the drag force - forAllConstIter - ( - phaseSystem::KdfTable, - Kdfs_, - KdfIter - ) + forAllConstIter(KdfTable, Kdfs_, KdfIter) { const surfaceScalarField& Kf(*KdfIter()); const phasePair& pair(this->phasePairs_[KdfIter.key()]); @@ -453,12 +494,7 @@ Foam::MomentumTransferPhaseSystem::AFfs() const } // Add the implicit part of the virtual mass force - forAllConstIter - ( - phaseSystem::VmfTable, - Vmfs_, - VmfIter - ) + forAllConstIter(VmfTable, Vmfs_, VmfIter) { const surfaceScalarField& Vmf(*VmfIter()); const phasePair& pair(this->phasePairs_[VmfIter.key()]); @@ -629,12 +665,7 @@ Foam::MomentumTransferPhaseSystem::phiFfs PtrList phiFfs(this->phaseModels_.size()); // Add the explicit part of the virtual mass force - forAllConstIter - ( - phaseSystem::VmfTable, - Vmfs_, - VmfIter - ) + forAllConstIter(VmfTable, Vmfs_, VmfIter) { const surfaceScalarField& Vmf(*VmfIter()); const phasePair& pair(this->phasePairs_[VmfIter.key()]); @@ -797,12 +828,7 @@ Foam::MomentumTransferPhaseSystem::phiKdPhis PtrList phiKdPhis(this->phaseModels_.size()); // Add the explicit part of the drag force - forAllConstIter - ( - phaseSystem::KdTable, - Kds_, - KdIter - ) + forAllConstIter(KdTable, Kds_, KdIter) { const volScalarField& K(*KdIter()); const phasePair& pair(this->phasePairs_[KdIter.key()]); @@ -844,12 +870,7 @@ Foam::MomentumTransferPhaseSystem::phiKdPhifs PtrList phiKdPhifs(this->phaseModels_.size()); // Add the explicit part of the drag force - forAllConstIter - ( - phaseSystem::KdfTable, - Kdfs_, - KdfIter - ) + forAllConstIter(KdfTable, Kdfs_, KdfIter) { const surfaceScalarField& Kf(*KdfIter()); const phasePair& pair(this->phasePairs_[KdfIter.key()]); @@ -891,12 +912,7 @@ Foam::MomentumTransferPhaseSystem::KdUByAs PtrList KdUByAs(this->phaseModels_.size()); // Add the explicit part of the drag force - forAllConstIter - ( - phaseSystem::KdTable, - Kds_, - KdIter - ) + forAllConstIter(KdTable, Kds_, KdIter) { const volScalarField& K(*KdIter()); const phasePair& pair(this->phasePairs_[KdIter.key()]); @@ -991,12 +1007,7 @@ Foam::MomentumTransferPhaseSystem::ddtCorrByAs // Add virtual mass correction if (includeVirtualMass) { - forAllConstIter - ( - phaseSystem::VmTable, - Vms_, - VmIter - ) + forAllConstIter(VmTable, Vms_, VmIter) { const volScalarField& Vm(*VmIter()); const phasePair& pair(this->phasePairs_[VmIter.key()]); @@ -1056,12 +1067,7 @@ void Foam::MomentumTransferPhaseSystem::partialElimination ); } - forAllConstIter - ( - phaseSystem::KdTable, - Kds_, - KdIter - ) + forAllConstIter(KdTable, Kds_, KdIter) { const volScalarField& K(*KdIter()); const phasePair& pair(this->phasePairs_[KdIter.key()]); @@ -1138,21 +1144,27 @@ void Foam::MomentumTransferPhaseSystem::partialElimination // Solve for the velocities and fluxes for (label i = 1; i < phases.size(); ++ i) { - for (label j = 0; j < i; j ++) + if (!phases[i].stationary()) { - phases[i].U() -= KdByAs[i][j]*phases[j].U(); - phases[i].phi() -= phiKds[i][j]*phases[j].phi(); + for (label j = 0; j < i; j ++) + { + phases[i].URef() -= KdByAs[i][j]*phases[j].U(); + phases[i].phiRef() -= phiKds[i][j]*phases[j].phi(); + } } } for (label i = phases.size() - 1; i >= 0; i --) { - for (label j = phases.size() - 1; j > i; j --) + if (!phases[i].stationary()) { - phases[i].U() -= KdByAs[i][j]*phases[j].U(); - phases[i].phi() -= phiKds[i][j]*phases[j].phi(); + for (label j = phases.size() - 1; j > i; j --) + { + phases[i].URef() -= KdByAs[i][j]*phases[j].U(); + phases[i].phiRef() -= phiKds[i][j]*phases[j].phi(); + } + phases[i].URef() /= KdByAs[i][i]; + phases[i].phiRef() /= phiKds[i][i]; } - phases[i].U() /= KdByAs[i][i]; - phases[i].phi() /= phiKds[i][i]; } } @@ -1179,15 +1191,10 @@ void Foam::MomentumTransferPhaseSystem::partialEliminationf ); } - forAllConstIter - ( - phaseSystem::KdfTable, - Kdfs_, - KdIter - ) + forAllConstIter(KdfTable, Kdfs_, KdfIter) { - const surfaceScalarField& K(*KdIter()); - const phasePair& pair(this->phasePairs_[KdIter.key()]); + const surfaceScalarField& K(*KdfIter()); + const phasePair& pair(this->phasePairs_[KdfIter.key()]); const label phase1i = pair.phase1().index(); const label phase2i = pair.phase2().index(); @@ -1239,18 +1246,24 @@ void Foam::MomentumTransferPhaseSystem::partialEliminationf // Solve for the fluxes for (label i = 1; i < phases.size(); ++ i) { - for (label j = 0; j < i; j ++) + if (!phases[i].stationary()) { - phases[i].phi() -= phiKdfs[i][j]*phases[j].phi(); + for (label j = 0; j < i; j ++) + { + phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi(); + } } } for (label i = phases.size() - 1; i >= 0; i --) { - for (label j = phases.size() - 1; j > i; j --) + if (!phases[i].stationary()) { - phases[i].phi() -= phiKdfs[i][j]*phases[j].phi(); + for (label j = phases.size() - 1; j > i; j --) + { + phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi(); + } + phases[i].phiRef() /= phiKdfs[i][i]; } - phases[i].phi() /= phiKdfs[i][i]; } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H index 8b7206b5d..c22b9f274 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H @@ -70,6 +70,34 @@ protected: // Protected typedefs + typedef HashPtrTable + < + volScalarField, + phasePairKey, + phasePairKey::hash + > KdTable; + + typedef HashPtrTable + < + surfaceScalarField, + phasePairKey, + phasePairKey::hash + > KdfTable; + + typedef HashPtrTable + < + volScalarField, + phasePairKey, + phasePairKey::hash + > VmTable; + + typedef HashPtrTable + < + surfaceScalarField, + phasePairKey, + phasePairKey::hash + > VmfTable; + typedef HashTable < autoPtr>, @@ -111,16 +139,16 @@ private: // Private data //- Drag coefficients - phaseSystem::KdTable Kds_; + KdTable Kds_; //- Face drag coefficients - phaseSystem::KdfTable Kdfs_; + KdfTable Kdfs_; //- Virtual mass coefficients - phaseSystem::VmTable Vms_; + VmTable Vms_; //- Face virtual mass coefficients - phaseSystem::VmfTable Vmfs_; + VmfTable Vmfs_; //- The phase diffusivities divided by the momentum coefficients HashPtrTable DByAfs_; @@ -154,6 +182,11 @@ private: //- Return the virtual mass coefficient for the phase pair virtual tmp Vm(const phasePairKey& key) const; + //- Add the mass-transfer-based momentum transfer to the equations + void addMassTransferMomentumTransfer + ( + phaseSystem::momentumTransferTable& eqns + ) const; public: @@ -172,12 +205,10 @@ public: //- Return the momentum transfer matrices for the cell-based algorithm. // This includes implicit and explicit forces that add into the cell // UEqn in the normal way. - virtual autoPtr - momentumTransfer() const; + virtual autoPtr momentumTransfer(); //- As momentumTransfer, but for the face-based algorithm - virtual autoPtr - momentumTransferf() const; + virtual autoPtr momentumTransferf(); //- Return implicit force coefficients on the faces, for the face-based // algorithm. diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C index 87b45df7b..e84f1fcc8 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C @@ -29,34 +29,20 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // template -void Foam::PopulationBalancePhaseSystem:: -addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const +Foam::tmp +Foam::PopulationBalancePhaseSystem::pDmdt +( + const phasePairKey& key +) const { - // Source term due to mass transfer - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) + if (!pDmdt_.found(key)) { - const phasePair& pair(phasePairIter()); - - if (pair.ordered()) - { - continue; - } - - const volVectorField& U1(pair.phase1().U()); - const volVectorField& U2(pair.phase2().U()); - - const volScalarField dmdt(this->pDmdt(pair)); - const volScalarField dmdt21(posPart(dmdt)); - const volScalarField dmdt12(negPart(dmdt)); - - *eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1); - *eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2); + return phaseSystem::dmdt(key); } + + const scalar pDmdtSign(Pair::compare(pDmdt_.find(key).key(), key)); + + return pDmdtSign**pDmdt_[key]; } @@ -152,112 +138,31 @@ Foam::PopulationBalancePhaseSystem:: template Foam::tmp -Foam::PopulationBalancePhaseSystem::pDmdt +Foam::PopulationBalancePhaseSystem::dmdt ( const phasePairKey& key ) const { - if (pDmdt_.found(key)) - { - const scalar pDmdtSign - ( - Pair::compare(pDmdt_.find(key).key(), key) - ); + return BasePhaseSystem::dmdt(key) + this->pDmdt(key); +} - return pDmdtSign**pDmdt_[key]; + +template +Foam::Xfer> +Foam::PopulationBalancePhaseSystem::dmdts() const +{ + PtrList dmdts(BasePhaseSystem::dmdts()); + + forAllConstIter(pDmdtTable, pDmdt_, pDmdtIter) + { + const phasePair& pair = this->phasePairs_[pDmdtIter.key()]; + const volScalarField& pDmdt = *pDmdtIter(); + + this->addField(pair.phase1(), "dmdt", pDmdt, dmdts); + this->addField(pair.phase2(), "dmdt", - pDmdt, dmdts); } - tmp tpDmdt - ( - new volScalarField - ( - IOobject - ( - "dmdt", - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); - - return tpDmdt; -} - - -template -Foam::tmp -Foam::PopulationBalancePhaseSystem::dmdt -( - const phaseModel& phase -) const -{ - tmp tDmdtFromKey - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("dmdtFromKey", phase.name()), - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); - - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) - { - const phasePair& pair(phasePairIter()); - - if (pair.ordered()) - { - continue; - } - - if (pair.contains(phase)) - { - tDmdtFromKey.ref() += this->dmdt - ( - phasePairKey(phase.name(), pair.otherPhase(phase).name(), false) - ); - } - } - - return tDmdtFromKey; -} - - -template -Foam::autoPtr -Foam::PopulationBalancePhaseSystem::momentumTransfer() const -{ - autoPtr - eqnsPtr(BasePhaseSystem::momentumTransfer()); - - addMomentumTransfer(eqnsPtr()); - - return eqnsPtr; -} - - -template -Foam::autoPtr -Foam::PopulationBalancePhaseSystem::momentumTransferf() const -{ - autoPtr - eqnsPtr(BasePhaseSystem::momentumTransferf()); - - addMomentumTransfer(eqnsPtr()); - - return eqnsPtr; + return dmdts.xfer(); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H index b6eb9658d..f03a860b8 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H @@ -58,25 +58,27 @@ class PopulationBalancePhaseSystem : public BasePhaseSystem { - protected: + // Protected typedefs + + typedef HashPtrTable + pDmdtTable; + + // Protected data //- populationBalanceModels PtrList populationBalances_; //- Interfacial Mass transfer rate - HashPtrTable pDmdt_; + pDmdtTable pDmdt_; // Protected member functions - //- Add the mass-based momentum transfer to the equations - void addMomentumTransfer - ( - phaseSystem::momentumTransferTable& eqns - ) const; + //- Return the population balance mass transfer rate + virtual tmp pDmdt(const phasePairKey& key) const; public: @@ -93,25 +95,11 @@ public: // Member Functions - //- Return the interfacial mass flow rate - virtual tmp pDmdt(const phasePairKey& key) const; + //- Return the mass transfer rate for a pair + virtual tmp dmdt(const phasePairKey& key) const; - //- Return the interfacial mass flow rate - virtual tmp dmdt(const phasePairKey& key) const - { - return BasePhaseSystem::dmdt(key) + this->pDmdt(key); - }; - - //- Return the total interfacial mass transfer rate for phase - virtual tmp dmdt(const phaseModel& phase) const; - - //- Return the momentum transfer matrices - virtual autoPtr - momentumTransfer() const; - - //- Return the momentum transfer matrices for the face based algorithm - virtual autoPtr - momentumTransferf() const; + //- Return the mass transfer rates for each phase + virtual Xfer> dmdts() const; //- Return the heat transfer matrices virtual autoPtr heatTransfer() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C index dacf1bd9b..c0149414e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C @@ -31,34 +31,38 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // template -void Foam::ThermalPhaseChangePhaseSystem:: -addMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const +Foam::tmp +Foam::ThermalPhaseChangePhaseSystem::iDmdt +( + const phasePairKey& key +) const { - // Source term due to mass trasfer - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) + if (!iDmdt_.found(key)) { - const phasePair& pair(phasePairIter()); - - if (pair.ordered()) - { - continue; - } - - const volVectorField& U1(pair.phase1().U()); - const volVectorField& U2(pair.phase2().U()); - - const volScalarField dmdt(this->iDmdt(pair) + this->wDmdt(pair)); - const volScalarField dmdt21(posPart(dmdt)); - const volScalarField dmdt12(negPart(dmdt)); - - *eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1); - *eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2); + return phaseSystem::dmdt(key); } + + const scalar dmdtSign(Pair::compare(iDmdt_.find(key).key(), key)); + + return dmdtSign**iDmdt_[key]; +} + + +template +Foam::tmp +Foam::ThermalPhaseChangePhaseSystem::wDmdt +( + const phasePairKey& key +) const +{ + if (!wDmdt_.found(key)) + { + return phaseSystem::dmdt(key); + } + + const scalar dmdtSign(Pair::compare(wDmdt_.find(key).key(), key)); + + return dmdtSign**wDmdt_[key]; } @@ -172,201 +176,42 @@ Foam::ThermalPhaseChangePhaseSystem::saturation() const } -template -Foam::tmp -Foam::ThermalPhaseChangePhaseSystem::iDmdt -( - const phasePairKey& key -) const -{ - const scalar dmdtSign(Pair::compare(iDmdt_.find(key).key(), key)); - return dmdtSign**iDmdt_[key]; -} - - -template -Foam::tmp -Foam::ThermalPhaseChangePhaseSystem::iDmdt -( - const Foam::phaseModel& phase -) const -{ - tmp tiDmdt - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("iDmdt", phase.name()), - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); - - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) - { - const phasePair& pair(phasePairIter()); - - if (pair.ordered()) - { - continue; - } - - if (pair.contains(phase)) - { - tiDmdt.ref() += this->iDmdt - ( - phasePairKey(phase.name(), pair.otherPhase(phase).name()) - ); - } - } - - return tiDmdt; -} - - -template -Foam::tmp -Foam::ThermalPhaseChangePhaseSystem::wDmdt -( - const phasePairKey& key -) const -{ - const scalar dmdtSign(Pair::compare(wDmdt_.find(key).key(), key)); - - return dmdtSign**wDmdt_[key]; -} - - -template -Foam::tmp -Foam::ThermalPhaseChangePhaseSystem::wDmdt -( - const Foam::phaseModel& phase -) const -{ - tmp twDmdt - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("wDmdt", phase.name()), - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); - - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) - { - const phasePair& pair(phasePairIter()); - - if (pair.ordered()) - { - continue; - } - - if (pair.contains(phase)) - { - twDmdt.ref() += this->wDmdt - ( - phasePairKey(phase.name(), pair.otherPhase(phase).name()) - ); - } - } - - return twDmdt; -} - - template Foam::tmp Foam::ThermalPhaseChangePhaseSystem::dmdt ( - const Foam::phaseModel& phase + const phasePairKey& key ) const { - tmp tDmdt - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("Dmdt", phase.name()), - this->mesh_.time().timeName(), - this->mesh_ - ), - this->mesh_, - dimensionedScalar("zero", dimDensity/dimTime, 0) - ) - ); + return BasePhaseSystem::dmdt(key) + this->iDmdt(key) + this->wDmdt(key); +} - forAllConstIter - ( - phaseSystem::phasePairTable, - this->phasePairs_, - phasePairIter - ) + +template +Foam::Xfer> +Foam::ThermalPhaseChangePhaseSystem::dmdts() const +{ + PtrList dmdts(BasePhaseSystem::dmdts()); + + forAllConstIter(iDmdtTable, iDmdt_, iDmdtIter) { - const phasePair& pair(phasePairIter()); + const phasePair& pair = this->phasePairs_[iDmdtIter.key()]; + const volScalarField& iDmdt = *iDmdtIter(); - if (pair.ordered()) - { - continue; - } - - if (pair.contains(phase)) - { - tDmdt.ref() += this->dmdt - ( - phasePairKey(phase.name(), pair.otherPhase(phase).name()) - ); - } + this->addField(pair.phase1(), "dmdt", iDmdt, dmdts); + this->addField(pair.phase2(), "dmdt", - iDmdt, dmdts); } - return tDmdt; -} + forAllConstIter(wDmdtTable, wDmdt_, wDmdtIter) + { + const phasePair& pair = this->phasePairs_[wDmdtIter.key()]; + const volScalarField& wDmdt = *wDmdtIter(); + this->addField(pair.phase1(), "dmdt", wDmdt, dmdts); + this->addField(pair.phase2(), "dmdt", - wDmdt, dmdts); + } -template -Foam::autoPtr -Foam::ThermalPhaseChangePhaseSystem::momentumTransfer() const -{ - autoPtr - eqnsPtr(BasePhaseSystem::momentumTransfer()); - - addMomentumTransfer(eqnsPtr()); - - return eqnsPtr; -} - - -template -Foam::autoPtr -Foam::ThermalPhaseChangePhaseSystem::momentumTransferf() const -{ - autoPtr - eqnsPtr(BasePhaseSystem::momentumTransferf()); - - addMomentumTransfer(eqnsPtr()); - - return eqnsPtr; + return dmdts.xfer(); } @@ -461,7 +306,6 @@ Foam::ThermalPhaseChangePhaseSystem::massTransfer() const forAll(Yi, i) { - if (Yi[i].member() != volatile_) { continue; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H index d5fef1acb..fef8a2680 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H @@ -58,6 +58,18 @@ class ThermalPhaseChangePhaseSystem protected: + // Protected typedefs + + typedef HashPtrTable + iDmdtTable; + + typedef HashPtrTable + wDmdtTable; + + typedef HashPtrTable + wMDotLTable; + + // Protected data //- Name of the volatile specie @@ -70,25 +82,22 @@ protected: Switch massTransfer_; //- Interfacial Mass transfer rate - HashPtrTable - iDmdt_; + iDmdtTable iDmdt_; //- Boundary Mass transfer rate - HashPtrTable - wDmdt_; + wDmdtTable wDmdt_; //- Boundary thermal energy transfer rate - HashPtrTable - wMDotL_; + wMDotLTable wMDotL_; // Protected member functions - //- Add the mass-based momentum transfer to the equations - void addMomentumTransfer - ( - phaseSystem::momentumTransferTable& eqns - ) const; + //- Return the interfacial mass transfer rate for a pair + tmp iDmdt(const phasePairKey& key) const; + + //- Return the boundary mass transfer rate for a pair + tmp wDmdt(const phasePairKey& key) const; public: @@ -108,35 +117,11 @@ public: //- Return the saturationModel const saturationModel& saturation() const; - //- Return the interfacial mass flow rate - virtual tmp iDmdt(const phasePairKey& key) const; + //- Return the mass transfer rate for a pair + virtual tmp dmdt(const phasePairKey& key) const; - //- Return the total interfacial mass transfer rate for phase - virtual tmp iDmdt(const phaseModel& phase) const; - - //- Return the boundary mass flow rate - virtual tmp wDmdt(const phasePairKey& key) const; - - //- Return the total boundary mass transfer rate for phase - virtual tmp wDmdt(const phaseModel& phase) const; - - //- Return the interfacial mass flow rate - virtual tmp dmdt(const phasePairKey& key) const - { - return BasePhaseSystem::dmdt(key) - + this->iDmdt(key) + this->wDmdt(key); - }; - - //- Return the total interfacial mass transfer rate for phase - virtual tmp dmdt(const phaseModel& phase) const; - - //- Return the momentum transfer matrices - virtual autoPtr - momentumTransfer() const; - - //- Return the momentum transfer matrices for the face based algorithm - virtual autoPtr - momentumTransferf() const; + //- Return the mass transfer rates for each phase + virtual Xfer> dmdts() const; //- Return the heat transfer matrices virtual autoPtr heatTransfer() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C index cd2b1c392..8c1f42077 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C @@ -26,63 +26,7 @@ License #include "AnisothermalPhaseModel.H" #include "phaseSystem.H" -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -Foam::AnisothermalPhaseModel::AnisothermalPhaseModel -( - const phaseSystem& fluid, - const word& phaseName, - const label index -) -: - BasePhaseModel(fluid, phaseName, index), - K_ - ( - IOobject - ( - IOobject::groupName("K", this->name()), - fluid.mesh().time().timeName(), - fluid.mesh() - ), - fluid.mesh(), - dimensionedScalar("K", sqr(dimVelocity), scalar(0)) - ) -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -template -Foam::AnisothermalPhaseModel::~AnisothermalPhaseModel() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -bool Foam::AnisothermalPhaseModel::compressible() const -{ - return !this->thermo().incompressible(); -} - - -template -void Foam::AnisothermalPhaseModel::correctKinematics() -{ - BasePhaseModel::correctKinematics(); - K_ = 0.5*magSqr(this->U()); -} - - -template -void Foam::AnisothermalPhaseModel::correctThermo() -{ - BasePhaseModel::correctThermo(); - - this->thermo_->correct(); -} - +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // template Foam::tmp @@ -111,18 +55,57 @@ Foam::AnisothermalPhaseModel::filterPressureWork } +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::AnisothermalPhaseModel::AnisothermalPhaseModel +( + const phaseSystem& fluid, + const word& phaseName, + const label index +) +: + BasePhaseModel(fluid, phaseName, index) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::AnisothermalPhaseModel::~AnisothermalPhaseModel() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::AnisothermalPhaseModel::correctThermo() +{ + BasePhaseModel::correctThermo(); + + this->thermo_->correct(); +} + + +template +bool Foam::AnisothermalPhaseModel::isothermal() const +{ + return false; +} + + template Foam::tmp Foam::AnisothermalPhaseModel::heEqn() { const volScalarField& alpha = *this; - const volVectorField& U = this->U(); - const surfaceScalarField& alphaPhi = this->alphaPhi(); - const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi(); - const volScalarField& contErr(this->continuityError()); + const volVectorField U(this->U()); + const surfaceScalarField alphaPhi(this->alphaPhi()); + const surfaceScalarField alphaRhoPhi(this->alphaRhoPhi()); - const volScalarField alphaEff(this->alphaEff()); + const volScalarField contErr(this->continuityError()); + const volScalarField K(this->K()); volScalarField& he = this->thermo_->he(); @@ -132,13 +115,13 @@ Foam::AnisothermalPhaseModel::heEqn() + fvm::div(alphaRhoPhi, he) - fvm::Sp(contErr, he) - + fvc::ddt(alpha, this->rho(), K_) + fvc::div(alphaRhoPhi, K_) - - contErr*K_ + + fvc::ddt(alpha, this->rho(), K) + fvc::div(alphaRhoPhi, K) + - contErr*K - fvm::laplacian ( fvc::interpolate(alpha) - *fvc::interpolate(alphaEff), + *fvc::interpolate(this->alphaEff()), he ) == @@ -163,12 +146,4 @@ Foam::AnisothermalPhaseModel::heEqn() } -template -const Foam::volScalarField& -Foam::AnisothermalPhaseModel::K() const -{ - return K_; -} - - // ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H index 95f7437f0..c15fe7107 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -44,7 +44,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class phaseModel Declaration + Class AnisothermalPhaseModel Declaration \*---------------------------------------------------------------------------*/ template @@ -52,12 +52,6 @@ class AnisothermalPhaseModel : public BasePhaseModel { - // Private data - - //- Kinetic energy - volScalarField K_; - - // Private member functions //- Optionally filter the pressure work term as the phase-fraction -> 0 @@ -85,20 +79,14 @@ public: // Member Functions - //- Return true if the phase is compressible otherwise false - virtual bool compressible() const; - - //- Correct the kinematics - virtual void correctKinematics(); - //- Correct the thermodynamics virtual void correctThermo(); + //- Return whether the phase is isothermal + virtual bool isothermal() const; + //- Return the enthalpy equation virtual tmp heEqn(); - - //- Return the phase kinetic energy - virtual const volScalarField& K() const; }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/InertPhaseModel/InertPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/InertPhaseModel/InertPhaseModel.H index c152da488..8404c19a0 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/InertPhaseModel/InertPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/InertPhaseModel/InertPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -44,7 +44,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class phaseModel Declaration + Class InertPhaseModel Declaration \*---------------------------------------------------------------------------*/ template @@ -73,10 +73,7 @@ public: // Thermo //- Return the fuel consumption rate matrix - virtual tmp R - ( - volScalarField& Yi - ) const; + virtual tmp R(volScalarField& Yi) const; //- Return the heat release rate virtual tmp Qdot() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C index 95a9f0c38..b743dfd7f 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -50,21 +50,25 @@ Foam::IsothermalPhaseModel::~IsothermalPhaseModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -bool Foam::IsothermalPhaseModel::compressible() const -{ - return !this->thermo().incompressible(); -} +void Foam::IsothermalPhaseModel::correctThermo() +{} template -void Foam::IsothermalPhaseModel::correctThermo() -{} +bool Foam::IsothermalPhaseModel::isothermal() const +{ + return true; +} template Foam::tmp Foam::IsothermalPhaseModel::heEqn() { + FatalErrorInFunction + << "Cannot construct an energy equation for an isothermal phase" + << exit(FatalError); + return tmp(); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H index f8bb151fe..143fcd076 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -45,7 +45,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class phaseModel Declaration + Class IsothermalPhaseModel Declaration \*---------------------------------------------------------------------------*/ template @@ -71,12 +71,12 @@ public: // Member Functions - //- Return true if the phase is compressible otherwise false - virtual bool compressible() const; - //- Correct the thermodynamics virtual void correctThermo(); + //- Return whether the phase is isothermal + virtual bool isothermal() const; + //- Return the enthalpy equation virtual tmp heEqn(); }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C index 740a66b19..18cf32167 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C @@ -179,19 +179,32 @@ Foam::MovingPhaseModel::MovingPhaseModel *this ) ), - continuityError_ + continuityErrorFlow_ ( IOobject ( - IOobject::groupName("continuityError", this->name()), + IOobject::groupName("continuityErrorFlow", this->name()), fluid.mesh().time().timeName(), fluid.mesh() ), fluid.mesh(), dimensionedScalar("0", dimDensity/dimTime, 0) - ) + ), + continuityErrorSources_ + ( + IOobject + ( + IOobject::groupName("continuityErrorSources", this->name()), + fluid.mesh().time().timeName(), + fluid.mesh() + ), + fluid.mesh(), + dimensionedScalar("0", dimDensity/dimTime, 0) + ), + K_(nullptr) { phi_.writeOpt() = IOobject::AUTO_WRITE; + correctKinematics(); } @@ -212,11 +225,11 @@ void Foam::MovingPhaseModel::correct() this->fluid().MRF().correctBoundaryVelocity(U_); - volScalarField& rho = this->thermo().rho(); + volScalarField& rho = this->thermoRef().rho(); - continuityError_ = - fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_) - - (this->fluid().fvOptions()(*this, rho)&rho); + continuityErrorFlow_ = fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_); + + continuityErrorSources_ = - (this->fluid().fvOptions()(*this, rho)&rho); } @@ -236,6 +249,12 @@ void Foam::MovingPhaseModel::correctKinematics() DUDtf_.clear(); DUDtf(); } + + if (K_.valid()) + { + K_.clear(); + K(); + } } @@ -257,19 +276,25 @@ void Foam::MovingPhaseModel::correctEnergyTransport() } +template +bool Foam::MovingPhaseModel::stationary() const +{ + return false; +} + + template Foam::tmp Foam::MovingPhaseModel::UEqn() { const volScalarField& alpha = *this; - volScalarField& rho = this->thermo().rho(); + const volScalarField& rho = this->thermo().rho(); return ( fvm::ddt(alpha, rho, U_) + fvm::div(alphaRhoPhi_, U_) - - fvm::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi_), U_) - + fvm::SuSp(this->fluid().fvOptions()(alpha, rho)&rho, U_) + + fvm::SuSp(- this->continuityError(), U_) + this->fluid().MRF().DDt(alpha*rho, U_) + turbulence_->divDevRhoReff(U_) ); @@ -283,13 +308,13 @@ Foam::MovingPhaseModel::UfEqn() // As the "normal" U-eqn but without the ddt terms const volScalarField& alpha = *this; - volScalarField& rho = this->thermo().rho(); + const volScalarField& rho = this->thermo().rho(); return ( fvm::div(alphaRhoPhi_, U_) - fvm::Sp(fvc::div(alphaRhoPhi_), U_) - + fvm::SuSp(this->fluid().fvOptions()(alpha, rho)&rho, U_) + + fvm::SuSp(- this->continuityErrorSources(), U_) + this->fluid().MRF().DDt(alpha*rho, U_) + turbulence_->divDevRhoReff(U_) ); @@ -306,12 +331,60 @@ Foam::MovingPhaseModel::U() const template Foam::volVectorField& -Foam::MovingPhaseModel::U() +Foam::MovingPhaseModel::URef() { return U_; } +template +Foam::tmp +Foam::MovingPhaseModel::phi() const +{ + return phi_; +} + + +template +Foam::surfaceScalarField& +Foam::MovingPhaseModel::phiRef() +{ + return phi_; +} + + +template +Foam::tmp +Foam::MovingPhaseModel::alphaPhi() const +{ + return alphaPhi_; +} + + +template +Foam::surfaceScalarField& +Foam::MovingPhaseModel::alphaPhiRef() +{ + return alphaPhi_; +} + + +template +Foam::tmp +Foam::MovingPhaseModel::alphaRhoPhi() const +{ + return alphaRhoPhi_; +} + + +template +Foam::surfaceScalarField& +Foam::MovingPhaseModel::alphaRhoPhiRef() +{ + return alphaRhoPhi_; +} + + template Foam::tmp Foam::MovingPhaseModel::DUDt() const @@ -339,76 +412,59 @@ Foam::MovingPhaseModel::DUDtf() const template -const Foam::tmp& -Foam::MovingPhaseModel::divU() const +Foam::tmp +Foam::MovingPhaseModel::continuityError() const { - return divU_; -} - - -template -void Foam::MovingPhaseModel::divU -( - const tmp& divU -) -{ - divU_ = divU; + return continuityErrorFlow_ + continuityErrorSources_; } template Foam::tmp -Foam::MovingPhaseModel::continuityError() const +Foam::MovingPhaseModel::continuityErrorFlow() const { - return continuityError_; + return continuityErrorFlow_; } template -Foam::tmp -Foam::MovingPhaseModel::phi() const +Foam::tmp +Foam::MovingPhaseModel::continuityErrorSources() const { - return phi_; + return continuityErrorSources_; } template -Foam::surfaceScalarField& -Foam::MovingPhaseModel::phi() +Foam::tmp +Foam::MovingPhaseModel::K() const { - return phi_; + if (!K_.valid()) + { + K_ = + new volScalarField + ( + IOobject::groupName("K", this->name()), + 0.5*magSqr(this->U()) + ); + } + + return tmp(K_()); } template -Foam::tmp -Foam::MovingPhaseModel::alphaPhi() const +Foam::tmp +Foam::MovingPhaseModel::divU() const { - return alphaPhi_; + return divU_.valid() ? tmp(divU_()) : tmp(); } template -Foam::surfaceScalarField& -Foam::MovingPhaseModel::alphaPhi() +void Foam::MovingPhaseModel::divU(tmp divU) { - return alphaPhi_; -} - - -template -Foam::tmp -Foam::MovingPhaseModel::alphaRhoPhi() const -{ - return alphaRhoPhi_; -} - - -template -Foam::surfaceScalarField& -Foam::MovingPhaseModel::alphaRhoPhi() -{ - return alphaRhoPhi_; + divU_ = divU; } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H index 387e485be..e60e601ab 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H @@ -55,7 +55,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class phaseModel Declaration + Class MovingPhaseModel Declaration \*---------------------------------------------------------------------------*/ template @@ -91,8 +91,14 @@ protected: //- Turbulence model autoPtr turbulence_; - //- Continuity error - volScalarField continuityError_; + //- Continuity error due to the flow + volScalarField continuityErrorFlow_; + + //- Continuity error due to any sources + volScalarField continuityErrorSources_; + + //- Kinetic Energy + mutable tmp K_; private: @@ -133,20 +139,41 @@ public: //- Correct the energy transport e.g. alphat virtual void correctEnergyTransport(); - //- Return the momentum equation - virtual tmp UEqn(); - - //- Return the momentum equation for the face-based algorithm - virtual tmp UfEqn(); - // Momentum - //- Constant access the velocity + //- Return whether the phase is stationary + virtual bool stationary() const; + + //- Return the momentum equation + virtual tmp UEqn(); + + //- Return the momentum equation for the face-based algorithm + virtual tmp UfEqn(); + + //- Return the velocity virtual tmp U() const; //- Access the velocity - virtual volVectorField& U(); + virtual volVectorField& URef(); + + //- Return the volumetric flux + virtual tmp phi() const; + + //- Access the volumetric flux + virtual surfaceScalarField& phiRef(); + + //- Return the volumetric flux of the phase + virtual tmp alphaPhi() const; + + //- Access the volumetric flux of the phase + virtual surfaceScalarField& alphaPhiRef(); + + //- Return the mass flux of the phase + virtual tmp alphaRhoPhi() const; + + //- Access the mass flux of the phase + virtual surfaceScalarField& alphaRhoPhiRef(); //- Return the substantive acceleration virtual tmp DUDt() const; @@ -154,35 +181,38 @@ public: //- Return the substantive acceleration on the faces virtual tmp DUDtf() const; - //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi)) - virtual const tmp& divU() const; - - //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi)) - virtual void divU(const tmp& divU); - - //- Constant access the continuity error + //- Return the continuity error virtual tmp continuityError() const; - //- Constant access the volumetric flux - virtual tmp phi() const; + //- Return the continuity error due to the flow field + virtual tmp continuityErrorFlow() const; - //- Access the volumetric flux - virtual surfaceScalarField& phi(); + //- Return the continuity error due to any sources + virtual tmp continuityErrorSources() const; - //- Constant access the volumetric flux of the phase - virtual tmp alphaPhi() const; - - //- Access the volumetric flux of the phase - virtual surfaceScalarField& alphaPhi(); - - //- Constant access the mass flux of the phase - virtual tmp alphaRhoPhi() const; - - //- Access the mass flux of the phase - virtual surfaceScalarField& alphaRhoPhi(); + //- Return the phase kinetic energy + virtual tmp K() const; - // Transport + // Compressibility (variable density) + + //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi)) + virtual tmp divU() const; + + //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi)) + virtual void divU(tmp divU); + + + // Transport (prevents compiler warnings) + + //- Return the effective thermal conductivity + using BasePhaseModel::kappaEff; + + //- Return the effective thermal conductivity for enthalpy + using BasePhaseModel::alphaEff; + + + // Turbulence //- Return the turbulent dynamic viscosity virtual tmp mut() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C index 534ebca40..7a279f75a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C @@ -69,6 +69,18 @@ Foam::MultiComponentPhaseModel::MultiComponentPhaseModel { inertIndex_ = this->thermo_->composition().species()[inertSpecie]; } + + PtrList& Y = this->thermo_->composition().Y(); + + forAll(Y, i) + { + if (i != inertIndex_ && this->thermo_->composition().active(i)) + { + const label j = YActive_.size(); + YActive_.resize(j + 1); + YActive_.set(j, &Y[i]); + } + } } @@ -96,7 +108,7 @@ void Foam::MultiComponentPhaseModel::correctThermo() dimensionedScalar("zero", dimless, 0) ); - PtrList& Yi = Y(); + PtrList& Yi = YRef(); forAll(Yi, i) { @@ -125,39 +137,19 @@ void Foam::MultiComponentPhaseModel::correctThermo() template -Foam::tmp -Foam::MultiComponentPhaseModel::YiEqn -( - volScalarField& Yi -) +bool Foam::MultiComponentPhaseModel::pure() const { - if - ( - (inertIndex_ != -1) - && ( - ( - Yi.name() - == IOobject::groupName - ( - this->thermo_->composition().species()[inertIndex_], - this->name() - ) - ) - || ( - !this->thermo_->composition().active - ( - this->thermo_->composition().species()[Yi.member()] - ) - ) - ) - ) - { - return tmp(); - } + return false; +} + +template +Foam::tmp +Foam::MultiComponentPhaseModel::YiEqn(volScalarField& Yi) +{ const volScalarField& alpha = *this; - const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi(); - const volScalarField& rho = this->rho(); + const surfaceScalarField alphaRhoPhi(this->alphaRhoPhi()); + const volScalarField rho(this->rho()); return ( @@ -190,10 +182,26 @@ Foam::MultiComponentPhaseModel::Y() const template Foam::PtrList& -Foam::MultiComponentPhaseModel::Y() +Foam::MultiComponentPhaseModel::YRef() { return this->thermo_->composition().Y(); } +template +const Foam::UPtrList& +Foam::MultiComponentPhaseModel::YActive() const +{ + return YActive_; +} + + +template +Foam::UPtrList& +Foam::MultiComponentPhaseModel::YActiveRef() +{ + return YActive_; +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.H index 6bbab1cf5..311548016 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -44,7 +44,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class phaseModel Declaration + Class MultiComponentPhaseModel Declaration \*---------------------------------------------------------------------------*/ template @@ -65,6 +65,9 @@ protected: //- Inert species index label inertIndex_; + //- Pointer list to active species + UPtrList YActive_; + public: @@ -87,17 +90,23 @@ public: //- Correct the thermodynamics virtual void correctThermo(); - //- Return the species fraction equation - virtual tmp YiEqn - ( - volScalarField& Yi - ); + //- Return whether the phase is pure (i.e., not multi-component) + virtual bool pure() const; - //- Constant access the species mass fractions + //- Return the species fraction equation + virtual tmp YiEqn(volScalarField& Yi); + + //- Return the species mass fractions virtual const PtrList& Y() const; //- Access the species mass fractions - virtual PtrList& Y(); + virtual PtrList& YRef(); + + //- Return the active species mass fractions + virtual const UPtrList& YActive() const; + + //- Access the active species mass fractions + virtual UPtrList& YActiveRef(); }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.C index 9d94f7627..b0a714ad5 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -50,13 +50,19 @@ Foam::PurePhaseModel::~PurePhaseModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -Foam::tmp -Foam::PurePhaseModel::YiEqn -( - volScalarField& Yi -) +bool Foam::PurePhaseModel::pure() const { - NotImplemented; + return true; +} + + +template +Foam::tmp +Foam::PurePhaseModel::YiEqn(volScalarField& Yi) +{ + FatalErrorInFunction + << "Cannot construct a species fraction equation for a pure phase" + << exit(FatalError); return tmp(); } @@ -66,14 +72,42 @@ template const Foam::PtrList& Foam::PurePhaseModel::Y() const { + // Y_ has never been set, so we are returning an empty list + return Y_; } template Foam::PtrList& -Foam::PurePhaseModel::Y() +Foam::PurePhaseModel::YRef() { + FatalErrorInFunction + << "Cannot access the species fractions of for a pure phase" + << exit(FatalError); + + return Y_; +} + + +template +const Foam::UPtrList& +Foam::PurePhaseModel::YActive() const +{ + // Y_ has never been set, so we are returning an empty list + + return Y_; +} + + +template +Foam::UPtrList& +Foam::PurePhaseModel::YActiveRef() +{ + FatalErrorInFunction + << "Cannot access the species fractions of for a pure phase" + << exit(FatalError); + return Y_; } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.H index ff42fbe71..3fe183297 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/PurePhaseModel/PurePhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -44,7 +44,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class phaseModel Declaration + Class PurePhaseModel Declaration \*---------------------------------------------------------------------------*/ template @@ -78,17 +78,25 @@ public: // Member Functions - //- Return the species fraction equation - virtual tmp YiEqn(volScalarField& Yi); + // Species + //- Return whether the phase is pure (i.e., not multi-component) + virtual bool pure() const; - // Thermo + //- Return the species fraction equation + virtual tmp YiEqn(volScalarField& Yi); //- Return the species mass fractions virtual const PtrList& Y() const; //- Access the species mass fractions - virtual PtrList& Y(); + virtual PtrList& YRef(); + + //- Return the active species mass fractions + virtual const UPtrList& YActive() const; + + //- Access the active species mass fractions + virtual UPtrList& YActiveRef(); }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ReactingPhaseModel/ReactingPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ReactingPhaseModel/ReactingPhaseModel.H index fdda42ca0..6cc8932d8 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ReactingPhaseModel/ReactingPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ReactingPhaseModel/ReactingPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -44,7 +44,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class phaseModel Declaration + Class ReactingPhaseModel Declaration \*---------------------------------------------------------------------------*/ template @@ -82,10 +82,7 @@ public: virtual void correctThermo(); //- Return the species fraction equation - virtual tmp R - ( - volScalarField& Yi - ) const; + virtual tmp R(volScalarField& Yi) const; //- Return heat release rate virtual tmp Qdot() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C new file mode 100644 index 000000000..8ac5d39de --- /dev/null +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C @@ -0,0 +1,367 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "StationaryPhaseModel.H" + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +template +template class PatchField, class GeoMesh> +Foam::tmp> +Foam::StationaryPhaseModel::zeroField +( + const word& name, + const dimensionSet& dims, + const bool cache +) const +{ + return tmp> + ( + new GeometricField + ( + IOobject + ( + IOobject::groupName(name, this->name()), + this->mesh().time().timeName(), + this->mesh() + ), + this->mesh(), + dimensioned("zero", dims, pTraits::zero) + ) + ); +} + + +template +template +Foam::tmp> +Foam::StationaryPhaseModel::zeroVolField +( + const word& name, + const dimensionSet& dims, + const bool cache +) const +{ + return zeroField(name, dims, cache); +} + + +template +template +Foam::tmp> +Foam::StationaryPhaseModel::zeroSurfaceField +( + const word& name, + const dimensionSet& dims, + const bool cache +) const +{ + return zeroField(name, dims, cache); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::StationaryPhaseModel::StationaryPhaseModel +( + const phaseSystem& fluid, + const word& phaseName, + const label index +) +: + BasePhaseModel(fluid, phaseName, index) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::StationaryPhaseModel::~StationaryPhaseModel() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +bool Foam::StationaryPhaseModel::stationary() const +{ + return true; +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::UEqn() +{ + FatalErrorInFunction + << "Cannot construct a momentum equation for a stationary phase" + << exit(FatalError); + + return tmp(); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::UfEqn() +{ + FatalErrorInFunction + << "Cannot construct a momentum equation for a stationary phase" + << exit(FatalError); + + return tmp(); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::U() const +{ + return zeroVolField("U", dimVelocity, true); +} + + +template +Foam::volVectorField& +Foam::StationaryPhaseModel::URef() +{ + FatalErrorInFunction + << "Cannot access the velocity of a stationary phase" + << exit(FatalError); + + return const_cast(volVectorField::null()); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::phi() const +{ + return zeroSurfaceField("phi", dimVolume/dimTime); +} + + +template +Foam::surfaceScalarField& +Foam::StationaryPhaseModel::phiRef() +{ + FatalErrorInFunction + << "Cannot access the flux of a stationary phase" + << exit(FatalError); + + return const_cast(surfaceScalarField::null()); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::alphaPhi() const +{ + return zeroSurfaceField("alphaPhi", dimVolume/dimTime); +} + + +template +Foam::surfaceScalarField& +Foam::StationaryPhaseModel::alphaPhiRef() +{ + FatalErrorInFunction + << "Cannot access the volumetric flux of a stationary phase" + << exit(FatalError); + + return const_cast(surfaceScalarField::null()); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::alphaRhoPhi() const +{ + return zeroSurfaceField("alphaRhoPhi", dimMass/dimTime); +} + + +template +Foam::surfaceScalarField& +Foam::StationaryPhaseModel::alphaRhoPhiRef() +{ + FatalErrorInFunction + << "Cannot access the mass flux of a stationary phase" + << exit(FatalError); + + return const_cast(surfaceScalarField::null()); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::DUDt() const +{ + return zeroVolField("DUDt", dimVelocity/dimTime); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::DUDtf() const +{ + return zeroSurfaceField("DUDtf", dimVelocity*dimArea/dimTime); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::continuityError() const +{ + return zeroVolField("continuityError", dimDensity/dimTime); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::continuityErrorFlow() const +{ + return zeroVolField("continuityErrorFlow", dimDensity/dimTime); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::continuityErrorSources() const +{ + return zeroVolField("continuityErrorSources", dimDensity/dimTime); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::K() const +{ + return zeroVolField("K", sqr(dimVelocity)); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::divU() const +{ + return tmp(); +} + + +template +void Foam::StationaryPhaseModel::divU +( + tmp divU +) +{ + FatalErrorInFunction + << "Cannot set the dilatation rate of a stationary phase" + << exit(FatalError); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::mut() const +{ + return zeroVolField("continuityError", dimDynamicViscosity); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::muEff() const +{ + return this->thermo().mu(); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::nut() const +{ + return zeroVolField("continuityError", dimViscosity); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::nuEff() const +{ + return this->thermo().nu(); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::kappaEff() const +{ + return this->thermo().kappa(); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::kappaEff(const label patchi) const +{ + return this->thermo().kappa(patchi); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::alphaEff() const +{ + return this->thermo().alpha(); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::alphaEff(const label patchi) const +{ + return this->thermo().alpha(patchi); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::k() const +{ + return zeroVolField("k", sqr(dimVelocity)); +} + + +template +Foam::tmp +Foam::StationaryPhaseModel::pPrime() const +{ + return zeroVolField("pPrime", dimPressure); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H new file mode 100644 index 000000000..b2d0822e8 --- /dev/null +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H @@ -0,0 +1,226 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::StationaryPhaseModel + +Description + ... + +SourceFiles + StationaryPhaseModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef StationaryPhaseModel_H +#define StationaryPhaseModel_H + +#include "phaseModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class StationaryPhaseModel Declaration +\*---------------------------------------------------------------------------*/ + +template +class StationaryPhaseModel +: + public BasePhaseModel +{ +private: + + // Private member functions + + //- Create a zero geometric field + template class PatchField, class GeoMesh> + tmp> zeroField + ( + const word& name, + const dimensionSet& dims, + const bool cache = false + ) const; + + //- Create a zero vol field + template + tmp> zeroVolField + ( + const word& name, + const dimensionSet& dims, + const bool cache = false + ) const; + + //- Create a zero surface field + template + tmp> zeroSurfaceField + ( + const word& name, + const dimensionSet& dims, + const bool cache = false + ) const; + + +public: + + // Constructors + + StationaryPhaseModel + ( + const phaseSystem& fluid, + const word& phaseName, + const label index + ); + + + //- Destructor + virtual ~StationaryPhaseModel(); + + + // Member Functions + + // Momentum + + //- Return whether the phase is stationary + virtual bool stationary() const; + + //- Return the momentum equation + virtual tmp UEqn(); + + //- Return the momentum equation for the face-based algorithm + virtual tmp UfEqn(); + + //- Return the velocity + virtual tmp U() const; + + //- Access the velocity + virtual volVectorField& URef(); + + //- Return the volumetric flux + virtual tmp phi() const; + + //- Access the volumetric flux + virtual surfaceScalarField& phiRef(); + + //- Return the volumetric flux of the phase + virtual tmp alphaPhi() const; + + //- Access the volumetric flux of the phase + virtual surfaceScalarField& alphaPhiRef(); + + //- Return the mass flux of the phase + virtual tmp alphaRhoPhi() const; + + //- Access the mass flux of the phase + virtual surfaceScalarField& alphaRhoPhiRef(); + + //- Return the substantive acceleration + virtual tmp DUDt() const; + + //- Return the substantive acceleration on the faces + virtual tmp DUDtf() const; + + //- Return the continuity error + virtual tmp continuityError() const; + + //- Return the continuity error due to the flow field + virtual tmp continuityErrorFlow() const; + + //- Return the continuity error due to any sources + virtual tmp continuityErrorSources() const; + + //- Return the phase kinetic energy + virtual tmp K() const; + + + // Compressibility (variable density) + + //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi)) + virtual tmp divU() const; + + //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi)) + virtual void divU(tmp divU); + + + // Transport (prevents compiler warnings) + + //- Return the effective thermal conductivity + using BasePhaseModel::kappaEff; + + //- Return the effective thermal conductivity for enthalpy + using BasePhaseModel::alphaEff; + + + // Turbulence + + //- Return the turbulent dynamic viscosity + virtual tmp mut() const; + + //- Return the effective dynamic viscosity + virtual tmp muEff() const; + + //- Return the turbulent kinematic viscosity + virtual tmp nut() const; + + //- Return the effective kinematic viscosity + virtual tmp nuEff() const; + + //- Return the effective thermal conductivity + virtual tmp kappaEff() const; + + //- Return the effective thermal conductivity on a patch + virtual tmp kappaEff(const label patchi) const; + + //- Return the effective thermal diffusivity for enthalpy + virtual tmp alphaEff() const; + + //- Return the effective thermal conductivity for enthalpy on a + // patch + virtual tmp alphaEff(const label patchi) const; + + //- Return the turbulent kinetic energy + virtual tmp k() const; + + //- Return the phase-pressure' + // (derivative of phase-pressure w.r.t. phase-fraction) + virtual tmp pPrime() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "StationaryPhaseModel.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C index 7755afec2..18f91ff5d 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -65,6 +65,13 @@ Foam::ThermoPhaseModel::~ThermoPhaseModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +bool Foam::ThermoPhaseModel::compressible() const +{ + return !thermo_().incompressible(); +} + + template const Foam::rhoThermo& Foam::ThermoPhaseModel::thermo() const @@ -75,7 +82,7 @@ Foam::ThermoPhaseModel::thermo() const template Foam::rhoThermo& -Foam::ThermoPhaseModel::thermo() +Foam::ThermoPhaseModel::thermoRef() { return thermo_(); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H index d4e405dd8..6882a0f8a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H @@ -49,7 +49,7 @@ namespace Foam class rhoThermo; /*---------------------------------------------------------------------------*\ - Class phaseModel Declaration + Class ThermoPhaseModel Declaration \*---------------------------------------------------------------------------*/ template @@ -85,12 +85,14 @@ public: // Thermo - //- Return const access to the thermophysical model + //- Return whether the phase is compressible + virtual bool compressible() const; + + //- Return the thermophysical model virtual const rhoThermo& thermo() const; - //- Return non-const access to the thermophysical model - // for correction - virtual rhoThermo& thermo(); + //- Access the thermophysical model + virtual rhoThermo& thermoRef(); //- Return the density field virtual tmp rho() const; @@ -101,28 +103,30 @@ public: //- Return the laminar dynamic viscosity virtual tmp mu() const; - //- Access the laminar dynamic viscosity + //- Return the laminar dynamic viscosity on a patch virtual tmp mu(const label patchi) const; //- Return the laminar kinematic viscosity virtual tmp nu() const; - //- Access the laminar kinematic viscosity + //- Return the laminar kinematic viscosity on a patch virtual tmp nu(const label patchi) const; //- Return the laminar thermal conductivity virtual tmp kappa() const; - //- Access the laminar thermal conductivity + //- Return the laminar thermal conductivity on a patch virtual tmp kappa(const label patchi) const; - //- Return the laminar thermal conductivity + //- Return the laminar thermal conductivity, given the turbulent + // thermal diffusivity virtual tmp kappaEff ( const volScalarField& alphat ) const; - //- Access the laminar thermal conductivity + //- Return the laminar thermal conductivity on a patch, given the + // turbulent thermal diffusivity virtual tmp kappaEff ( const scalarField& alphat, @@ -135,18 +139,29 @@ public: //- Return the thermal diffusivity for enthalpy on a patch virtual tmp alpha(const label patchi) const; - //- Return the thermal diffusivity for enthalpy + //- Return the effective thermal diffusivity for enthalpy, given the + // turbulent thermal diffusivity virtual tmp alphaEff ( const volScalarField& alphat ) const; - //- Return the thermal diffusivity for enthalpy on a patch + //- Return the effective thermal diffusivity for enthalpy on a + // patch, given the turbulent thermal diffusivity virtual tmp alphaEff ( const scalarField& alphat, const label patchi ) const; + + + // Turbulence (prevents compiler warnings) + + //- Return the effective thermal conductivity + using BasePhaseModel::kappaEff; + + //- Return the effective thermal conductivity for enthalpy + using BasePhaseModel::alphaEff; }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C index 3b4dd6d4d..0de3d7782 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C @@ -165,12 +165,6 @@ bool Foam::phaseModel::read() } -bool Foam::phaseModel::compressible() const -{ - return false; -} - - void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const { surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef(); @@ -189,27 +183,4 @@ void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const } -const Foam::tmp& Foam::phaseModel::divU() const -{ - NotImplemented; - static tmp divU_(nullptr); - return divU_; -} - - -void Foam::phaseModel::divU(const tmp& divU) -{ - WarningInFunction - << "Attempt to set the dilatation rate of an incompressible phase" - << endl; -} - - -const Foam::volScalarField& Foam::phaseModel::K() const -{ - NotImplemented; - return volScalarField::null(); -} - - // ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H index 6dd4d3760..dc3a11c5e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H @@ -195,20 +195,11 @@ public: //- Correct the turbulence virtual void correctTurbulence(); - //- Correct the energy transport e.g. alphat + //- Correct the energy transport virtual void correctEnergyTransport(); - //- Return the momentum equation - virtual tmp UEqn() = 0; - - //- Return the momentum equation for the face-based algorithm - virtual tmp UfEqn() = 0; - - //- Return the enthalpy equation - virtual tmp heEqn() = 0; - - //- Return the species fraction equation - virtual tmp YiEqn(volScalarField& Yi) = 0; + //- Ensure that the flux at inflow/outflow BCs is preserved + void correctInflowOutflow(surfaceScalarField& alphaPhi) const; //- Read phase properties dictionary virtual bool read(); @@ -217,44 +208,88 @@ public: // Compressibility (variable density) //- Return true if the phase is compressible otherwise false - virtual bool compressible() const; + virtual bool compressible() const = 0; //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi)) - virtual const tmp& divU() const; + virtual tmp divU() const = 0; //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi)) - virtual void divU(const tmp& divU); - - //- Return the phase kinetic energy - virtual const volScalarField& K() const; + virtual void divU(tmp divU) = 0; // Thermo - //- Return const access to the thermophysical model + //- Return whether the phase is isothermal + virtual bool isothermal() const = 0; + + //- Return the enthalpy equation + virtual tmp heEqn() = 0; + + //- Return the thermophysical model virtual const rhoThermo& thermo() const = 0; - //- Return non-const access to the thermophysical model - // for correction - virtual rhoThermo& thermo() = 0; + //- Access the thermophysical model + virtual rhoThermo& thermoRef() = 0; //- Return the density field virtual tmp rho() const = 0; - //- Constant access the species mass fractions + + // Species + + //- Return whether the phase is pure (i.e., not multi-component) + virtual bool pure() const = 0; + + //- Return the species fraction equation + virtual tmp YiEqn(volScalarField& Yi) = 0; + + //- Return the species mass fractions virtual const PtrList& Y() const = 0; //- Access the species mass fractions - virtual PtrList& Y() = 0; + virtual PtrList& YRef() = 0; + + //- Return the active species mass fractions + virtual const UPtrList& YActive() const = 0; + + //- Access the active species mass fractions + virtual UPtrList& YActiveRef() = 0; // Momentum - //- Constant access the velocity + //- Return whether the phase is stationary + virtual bool stationary() const = 0; + + //- Return the momentum equation + virtual tmp UEqn() = 0; + + //- Return the momentum equation for the face-based algorithm + virtual tmp UfEqn() = 0; + + //- Return the velocity virtual tmp U() const = 0; //- Access the velocity - virtual volVectorField& U() = 0; + virtual volVectorField& URef() = 0; + + //- Return the volumetric flux + virtual tmp phi() const = 0; + + //- Access the volumetric flux + virtual surfaceScalarField& phiRef() = 0; + + //- Return the volumetric flux of the phase + virtual tmp alphaPhi() const = 0; + + //- Access the volumetric flux of the phase + virtual surfaceScalarField& alphaPhiRef() = 0; + + //- Return the mass flux of the phase + virtual tmp alphaRhoPhi() const = 0; + + //- Access the mass flux of the phase + virtual surfaceScalarField& alphaRhoPhiRef() = 0; //- Return the substantive acceleration virtual tmp DUDt() const = 0; @@ -262,29 +297,17 @@ public: //- Return the substantive acceleration on the faces virtual tmp DUDtf() const = 0; - //- Constant access the continuity error + //- Return the continuity error virtual tmp continuityError() const = 0; - //- Constant access the volumetric flux - virtual tmp phi() const = 0; + //- Return the continuity error due to the flow field + virtual tmp continuityErrorFlow() const = 0; - //- Access the volumetric flux - virtual surfaceScalarField& phi() = 0; + //- Return the continuity error due to any sources + virtual tmp continuityErrorSources() const = 0; - //- Constant access the volumetric flux of the phase - virtual tmp alphaPhi() const = 0; - - //- Access the volumetric flux of the phase - virtual surfaceScalarField& alphaPhi() = 0; - - //- Constant access the mass flux of the phase - virtual tmp alphaRhoPhi() const = 0; - - //- Access the mass flux of the phase - virtual surfaceScalarField& alphaRhoPhi() = 0; - - //- Ensure that the flux at inflow/outflow BCs is preserved - void correctInflowOutflow(surfaceScalarField& alphaPhi) const; + //- Return the phase kinetic energy + virtual tmp K() const = 0; // Transport diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModels.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModels.C index 0b2654998..22e2933be 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModels.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModels.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -39,6 +39,7 @@ License #include "InertPhaseModel.H" #include "ReactingPhaseModel.H" #include "MovingPhaseModel.H" +#include "StationaryPhaseModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -68,6 +69,30 @@ namespace Foam purePhaseModel ); + typedef + AnisothermalPhaseModel + < + PurePhaseModel + < + InertPhaseModel + < + StationaryPhaseModel + < + ThermoPhaseModel + > + > + > + > + pureStationaryPhaseModel; + + addNamedToRunTimeSelectionTable + ( + phaseModel, + pureStationaryPhaseModel, + phaseSystem, + pureStationaryPhaseModel + ); + typedef IsothermalPhaseModel < @@ -92,6 +117,30 @@ namespace Foam pureIsothermalPhaseModel ); + typedef + IsothermalPhaseModel + < + PurePhaseModel + < + InertPhaseModel + < + StationaryPhaseModel + < + ThermoPhaseModel + > + > + > + > + pureStationaryIsothermalPhaseModel; + + addNamedToRunTimeSelectionTable + ( + phaseModel, + pureStationaryIsothermalPhaseModel, + phaseSystem, + pureStationaryIsothermalPhaseModel + ); + typedef AnisothermalPhaseModel < diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C index 6f74672be..d1a15f73b 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -155,6 +155,50 @@ Foam::phaseSystem::phaseSystem MRF_(mesh_) { + // Groupings + label movingPhasei = 0; + label stationaryPhasei = 0; + label anisothermalPhasei = 0; + label multiComponentPhasei = 0; + forAll(phaseModels_, phasei) + { + phaseModel& phase = phaseModels_[phasei]; + movingPhasei += !phase.stationary(); + stationaryPhasei += phase.stationary(); + anisothermalPhasei += !phase.isothermal(); + multiComponentPhasei += !phase.pure(); + } + movingPhaseModels_.resize(movingPhasei); + stationaryPhaseModels_.resize(stationaryPhasei); + anisothermalPhaseModels_.resize(anisothermalPhasei); + multiComponentPhaseModels_.resize(multiComponentPhasei); + + movingPhasei = 0; + stationaryPhasei = 0; + anisothermalPhasei = 0; + multiComponentPhasei = 0; + forAll(phaseModels_, phasei) + { + phaseModel& phase = phaseModels_[phasei]; + if (!phase.stationary()) + { + movingPhaseModels_.set(movingPhasei ++, &phase); + } + if (phase.stationary()) + { + stationaryPhaseModels_.set(stationaryPhasei ++, &phase); + } + if (!phase.isothermal()) + { + anisothermalPhaseModels_.set(anisothermalPhasei ++, &phase); + } + if (!phase.pure()) + { + multiComponentPhaseModels_.set(multiComponentPhasei ++, &phase); + } + } + + // Write phi phi_.writeOpt() = IOobject::AUTO_WRITE; // Blending methods @@ -175,6 +219,7 @@ Foam::phaseSystem::phaseSystem generatePairsAndSubModels("surfaceTension", surfaceTensionModels_); generatePairsAndSubModels("aspectRatio", aspectRatioModels_); + // Update motion fields correctKinematics(); } @@ -279,22 +324,18 @@ Foam::phaseSystem::sigma(const phasePairKey& key) const } -void Foam::phaseSystem::solve() -{} - - Foam::tmp Foam::phaseSystem::dmdt ( -const Foam::phaseModel& phase + const phasePairKey& key ) const { - tmp tDmdt + return tmp ( new volScalarField ( IOobject ( - IOobject::groupName("dmdt", phase.name()), + IOobject::groupName("dmdt", phasePairs_[key]->name()), this->mesh_.time().timeName(), this->mesh_ ), @@ -302,11 +343,21 @@ const Foam::phaseModel& phase dimensionedScalar("zero", dimDensity/dimTime, 0) ) ); - - return tDmdt; } +Foam::Xfer> Foam::phaseSystem::dmdts() const +{ + PtrList dmdts(this->phaseModels_.size()); + + return dmdts.xfer(); +} + + +void Foam::phaseSystem::solve() +{} + + void Foam::phaseSystem::correct() { forAll(phaseModels_, phasei) diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H index fe08c1015..19426058f 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H @@ -73,71 +73,16 @@ public: // Public typedefs - typedef - HashPtrTable - < - volScalarField, - phasePairKey, - phasePairKey::hash - > - KdTable; + typedef HashPtrTable momentumTransferTable; - typedef - HashPtrTable - < - surfaceScalarField, - phasePairKey, - phasePairKey::hash - > - KdfTable; + typedef HashPtrTable heatTransferTable; - typedef - HashPtrTable - < - volScalarField, - phasePairKey, - phasePairKey::hash - > - VmTable; - - typedef - HashPtrTable - < - surfaceScalarField, - phasePairKey, - phasePairKey::hash - > - VmfTable; - - typedef - HashPtrTable - < - fvVectorMatrix, - word, - string::hash - > - momentumTransferTable; - - typedef - HashPtrTable - < - fvScalarMatrix, - word, - string::hash - > - heatTransferTable; - - typedef - HashPtrTable - < - fvScalarMatrix, - word, - string::hash - > - massTransferTable; + typedef HashPtrTable massTransferTable; typedef PtrListDictionary phaseModelList; + typedef UPtrList phaseModelPartialList; + protected: @@ -182,6 +127,18 @@ protected: //- Phase models phaseModelList phaseModels_; + //- Moving phase models + phaseModelPartialList movingPhaseModels_; + + //- Stationary phase models + phaseModelPartialList stationaryPhaseModels_; + + //- Anisothermal phase models + phaseModelPartialList anisothermalPhaseModels_; + + //- Multi-component phase models + phaseModelPartialList multiComponentPhaseModels_; + //- Phase pairs phasePairTable phasePairs_; @@ -317,24 +274,6 @@ protected: HashPtrTable& fieldTable ) const; - //- Fill up gaps in a phase-indexed list with zeros - template class PatchField, class GeoMesh> - void fillFields - ( - const word& name, - const dimensionSet& dims, - PtrList>& fieldList - ) const; - - //- Fill up gaps in a phase-indexed table with zeros - template class PatchField, class GeoMesh> - void fillFields - ( - const word& name, - const dimensionSet& dims, - HashPtrTable>& fieldTable - ) const; - public: @@ -357,83 +296,230 @@ public: // Member Functions - //- Constant access the mesh - inline const fvMesh& mesh() const; + // Access - //- Constant access the phase models - inline const phaseModelList& phases() const; + //- Return the mesh + inline const fvMesh& mesh() const; - //- Access the phase models - inline phaseModelList& phases(); + //- Return the phase models + inline const phaseModelList& phases() const; - //- Constant access the phase pairs - inline const phasePairTable& phasePairs() const; + //- Access the phase models + inline phaseModelList& phases(); - //- Return the mixture density - tmp rho() const; + //- Return the models for phases that are moving + inline const phaseModelPartialList& movingPhases() const; - //- Return the mixture velocity - tmp U() const; + //- Access the models for phases that are moving + inline phaseModelPartialList& movingPhases(); - //- Constant access the mixture flux - inline const surfaceScalarField& phi() const; + //- Return the models for phases that are stationary + inline const phaseModelPartialList& stationaryPhases() const; - //- Access the mixture flux - inline surfaceScalarField& phi(); + //- Access the models for phases that are stationary + inline phaseModelPartialList& stationaryPhases(); - //- Constant access the rate of change of the pressure - inline const volScalarField& dpdt() const; + //- Return the models for phases that have variable temperature + inline const phaseModelPartialList& anisothermalPhases() const; - //- Access the rate of change of the pressure - inline volScalarField& dpdt(); + //- Access the models for phases that have variable temperature + inline phaseModelPartialList& anisothermalPhases(); - //- Return the aspect-ratio - tmp E(const phasePairKey& key) const; + //- Return the models for phases that have multiple species + inline const phaseModelPartialList& multiComponentPhases() const; - //- Return the surface tension coefficient - tmp sigma(const phasePairKey& key) const; + //- Access the models for phases that have multiple species + inline phaseModelPartialList& multiComponentPhases(); - //- Return MRF zones - inline const IOMRFZoneList& MRF() const; + //- Return the phase pairs + inline const phasePairTable& phasePairs() const; - //- Optional FV-options - inline fv::options& fvOptions() const; + //- Return the mixture flux + inline const surfaceScalarField& phi() const; - //- Access a sub model between a phase pair - template - const modelType& lookupSubModel(const phasePair& key) const; + //- Access the mixture flux + inline surfaceScalarField& phi(); - //- Access a sub model between two phases - template - const modelType& lookupSubModel - ( - const phaseModel& dispersed, - const phaseModel& continuous - ) const; + //- Return the rate of change of the pressure + inline const volScalarField& dpdt() const; - //- Solve for the phase fractions - virtual void solve(); + //- Access the rate of change of the pressure + inline volScalarField& dpdt(); - //- Return the total interfacial mass transfer rate for a phase - virtual tmp dmdt(const phaseModel& phase) const; + //- Return MRF zones + inline const IOMRFZoneList& MRF() const; - //- Correct the fluid properties other than the thermo and turbulence - virtual void correct(); + //- Access the fvOptions + inline fv::options& fvOptions() const; - //- Correct the kinematics - virtual void correctKinematics(); - //- Correct the thermodynamics - virtual void correctThermo(); + // Sub-model lookup - //- Correct the turbulence - virtual void correctTurbulence(); + //- Return a sub model between a phase pair + template + const modelType& lookupSubModel(const phasePair& key) const; - //- Correct the energy transport e.g. alphat - virtual void correctEnergyTransport(); + //- Return a sub model between two phases + template + const modelType& lookupSubModel + ( + const phaseModel& dispersed, + const phaseModel& continuous + ) const; - //- Read base phaseProperties dictionary - virtual bool read(); + + // Field construction + + //- Fill up gaps in a phase-indexed list of fields with zeros + template + < + class Type, + template class PatchField, + class GeoMesh + > + void fillFields + ( + const word& name, + const dimensionSet& dims, + PtrList>& fieldList + ) const; + + //- Fill up gaps in a phase-indexed table of fields with zeros + template + < + class Type, + template class PatchField, + class GeoMesh + > + void fillFields + ( + const word& name, + const dimensionSet& dims, + HashPtrTable>& + fieldTable + ) const; + + + // Properties + + //- Return the mixture density + tmp rho() const; + + //- Return the mixture velocity + tmp U() const; + + //- Return the aspect-ratio for a pair + tmp E(const phasePairKey& key) const; + + //- Return the surface tension coefficient for a pair + tmp sigma(const phasePairKey& key) const; + + //- Return the mass transfer rate for a pair + virtual tmp dmdt(const phasePairKey& key) const; + + //- Return the mass transfer rates for each phase + virtual Xfer> dmdts() const; + + + // Transfers + + //- Return the momentum transfer matrices for the cell-based + // algorithm + virtual autoPtr momentumTransfer() = 0; + + //- Return the momentum transfer matrices for the face-based + // algiorithm + virtual autoPtr momentumTransferf() = 0; + + //- Return the implicit force coefficients for the face-based + // algorithm + virtual Xfer> AFfs() const = 0; + + //- Return the force fluxes for the cell-based algorithm + virtual Xfer> phiFs + ( + const PtrList& rAUs + ) = 0; + + //- Return the force fluxes for the face-based algorithm + virtual Xfer> phiFfs + ( + const PtrList& rAUfs + ) = 0; + + //- Return the force fluxes for the cell-based algorithm + virtual Xfer> phiKdPhis + ( + const PtrList& rAUs + ) const = 0; + + //- Return the force fluxes for the face-based algorithm + virtual Xfer> phiKdPhifs + ( + const PtrList& rAUfs + ) const = 0; + + //- Return the explicit part of the drag force + virtual Xfer> KdUByAs + ( + const PtrList& rAUs + ) const = 0; + + //- Solve the drag system for the new velocities and fluxes + virtual void partialElimination + ( + const PtrList& rAUs + ) = 0; + + //- Solve the drag system for the new fluxes + virtual void partialEliminationf + ( + const PtrList& rAUfs + ) = 0; + + //- Return the flux corrections for the cell-based algorithm + virtual Xfer> ddtCorrByAs + ( + const PtrList& rAUs, + const bool includeVirtualMass = false + ) const = 0; + + //- Return the phase diffusivities divided by the momentum + // coefficients + virtual const HashPtrTable& DByAfs() const = 0; + + //- Return the heat transfer matrices + virtual autoPtr heatTransfer() const = 0; + + //- Return the mass transfer matrices + virtual autoPtr massTransfer() const = 0; + + + // Evolution + + //- Solve for the phase fractions + virtual void solve(); + + //- Correct the fluid properties other than those listed below + virtual void correct(); + + //- Correct the kinematics + virtual void correctKinematics(); + + //- Correct the thermodynamics + virtual void correctThermo(); + + //- Correct the turbulence + virtual void correctTurbulence(); + + //- Correct the energy transport e.g. alphat + virtual void correctEnergyTransport(); + + + // IO + + //- Read base phaseProperties dictionary + virtual bool read(); }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H index 7f0e73be5..df9ff4539 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -45,6 +45,62 @@ Foam::phaseSystem::phases() } +inline const Foam::phaseSystem::phaseModelPartialList& +Foam::phaseSystem::movingPhases() const +{ + return movingPhaseModels_; +} + + +inline Foam::phaseSystem::phaseModelPartialList& +Foam::phaseSystem::movingPhases() +{ + return movingPhaseModels_; +} + + +inline const Foam::phaseSystem::phaseModelPartialList& +Foam::phaseSystem::stationaryPhases() const +{ + return stationaryPhaseModels_; +} + + +inline Foam::phaseSystem::phaseModelPartialList& +Foam::phaseSystem::stationaryPhases() +{ + return stationaryPhaseModels_; +} + + +inline const Foam::phaseSystem::phaseModelPartialList& +Foam::phaseSystem::anisothermalPhases() const +{ + return anisothermalPhaseModels_; +} + + +inline Foam::phaseSystem::phaseModelPartialList& +Foam::phaseSystem::anisothermalPhases() +{ + return anisothermalPhaseModels_; +} + + +inline const Foam::phaseSystem::phaseModelPartialList& +Foam::phaseSystem::multiComponentPhases() const +{ + return multiComponentPhaseModels_; +} + + +inline Foam::phaseSystem::phaseModelPartialList& +Foam::phaseSystem::multiComponentPhases() +{ + return multiComponentPhaseModels_; +} + + inline const Foam::phaseSystem::phasePairTable& Foam::phaseSystem::phasePairs() const { diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C index 6a67bf0b8..69aebd5ed 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C @@ -281,6 +281,8 @@ void Foam::phaseSystem::addField } +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + template class PatchField, class GeoMesh> void Foam::phaseSystem::fillFields ( diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C index 43c461648..2eb2de358 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C @@ -87,10 +87,8 @@ void Foam::diameterModels::driftModels::constantDrift::driftRate ) { const sizeGroup& fi = *popBal_.sizeGroups()[i]; - phaseModel& phase = const_cast(fi.phase()); - volScalarField& rho = phase.thermo().rho(); - driftRate += (popBal_.fluid().fvOptions()(phase, rho)&rho)/(N_*rho); + driftRate -= fi.phase().continuityErrorSources()/(N_*fi.phase().rho()); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/constantNucleation/constantNucleation.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/constantNucleation/constantNucleation.C index bf5cda759..8e971f215 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/constantNucleation/constantNucleation.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/constantNucleation/constantNucleation.C @@ -105,12 +105,10 @@ nucleationRate ) { const sizeGroup& fi = *popBal_.sizeGroups()[i]; - phaseModel& phase = const_cast(fi.phase()); - volScalarField& rho = phase.thermo().rho(); - nucleationRate += + nucleationRate -= popBal_.gamma(i, velGroup_.formFactor()*pow3(d_)) - *(popBal_.fluid().fvOptions()(phase, rho)&rho)/rho/fi.x(); + *fi.phase().continuityErrorSources()/fi.phase().rho()/fi.x(); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C index 210bf51b3..cd7ae6637 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C @@ -151,8 +151,8 @@ nucleationRate ) { const sizeGroup& fi = *popBal_.sizeGroups()[i]; - phaseModel& phase = const_cast(fi.phase()); - volScalarField& rho = phase.thermo().rho(); + const phaseModel& phase = fi.phase(); + const volScalarField& rho = phase.rho(); const tmp talphat(turbulence_.alphat()); const volScalarField::Boundary& alphatBf = talphat().boundaryField(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H index dc53815cf..22d875dce 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H @@ -7,32 +7,27 @@ for (int Ecorr=0; Ecorr EEqn(phase.heEqn()); + fvScalarMatrix EEqn + ( + phase.heEqn() + == + *heatTransfer[phase.name()] + + alpha*rho*(U&g) + + fvOptions(alpha, rho, phase.thermoRef().he()) + ); - if (EEqn.valid()) - { - EEqn = - ( - EEqn - == - *heatTransfer[phase.name()] - + alpha*rho*(U&g) - + fvOptions(alpha, rho, phase.thermo().he()) - ); - - EEqn->relax(); - fvOptions.constrain(EEqn.ref()); - EEqn->solve(); - fvOptions.correct(phase.thermo().he()); - } + EEqn.relax(); + fvOptions.constrain(EEqn); + EEqn.solve(); + fvOptions.correct(phase.thermoRef().he()); } fluid.correctThermo(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H index 4094e185a..02a7474a9 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H @@ -5,31 +5,26 @@ phaseSystem::massTransferTable& massTransfer(massTransferPtr()); - forAll(phases, phasei) + forAll(fluid.multiComponentPhases(), multiComponentPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.multiComponentPhases()[multiComponentPhasei]; - PtrList& Y = phase.Y(); + UPtrList& Y = phase.YActiveRef(); const volScalarField& alpha = phase; const volScalarField& rho = phase.rho(); forAll(Y, i) { - tmp YiEqn(phase.YiEqn(Y[i])); + fvScalarMatrix YiEqn + ( + phase.YiEqn(Y[i]) + == + *massTransfer[Y[i].name()] + + fvOptions(alpha, rho, Y[i]) + ); - if (YiEqn.valid()) - { - YiEqn = - ( - YiEqn - == - *massTransfer[Y[i].name()] - + fvOptions(alpha, rho, Y[i]) - ); - - YiEqn->relax(); - YiEqn->solve(mesh.solver("Yi")); - } + YiEqn.relax(); + YiEqn.solve(mesh.solver("Yi")); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H index 1fc5b19ca..50f1ed519 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H @@ -20,7 +20,7 @@ dimensionedScalar pMin #include "gh.H" -volScalarField& p = phases[0].thermo().p(); +volScalarField& p = phases[0].thermoRef().p(); Info<< "Reading field p_rgh\n" << endl; volScalarField p_rgh diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 79aa52e03..def9e81c2 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -28,6 +28,7 @@ License #include "MULES.H" #include "subCycle.H" +#include "UniformField.H" #include "fvcDdt.H" #include "fvcDiv.H" @@ -74,32 +75,75 @@ void Foam::multiphaseSystem::solveAlphas() phases()[phasei].correctBoundaryConditions(); } - PtrList alphaPhiCorrs(phases().size()); + // Calculate the void fraction + volScalarField alphaVoid + ( + IOobject + ( + "alphaVoid", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("alphaVoid", dimless, 1) + ); + forAll(stationaryPhases(), stationaryPhasei) + { + alphaVoid -= stationaryPhases()[stationaryPhasei]; + } + + // Generate face-alphas + PtrList alphafs(phases().size()); forAll(phases(), phasei) { phaseModel& phase = phases()[phasei]; - volScalarField& alpha1 = phase; - - alphaPhiCorrs.set + alphafs.set ( phasei, new surfaceScalarField ( - "phi" + alpha1.name() + "Corr", - fvc::flux - ( - phi_, - phase, - "div(phi," + alpha1.name() + ')' - ) + IOobject::groupName("alphaf", phase.name()), + upwind(mesh_, phi_).interpolate(phase) + ) + ); + } + + // Create correction fluxes + PtrList alphaPhiCorrs(phases().size()); + forAll(stationaryPhases(), stationaryPhasei) + { + phaseModel& phase = stationaryPhases()[stationaryPhasei]; + + alphaPhiCorrs.set + ( + phase.index(), + new surfaceScalarField + ( + IOobject::groupName("alphaPhiCorr", phase.name()), + - upwind(mesh_, phi_).flux(phase) + ) + ); + } + forAll(movingPhases(), movingPhasei) + { + phaseModel& phase = movingPhases()[movingPhasei]; + volScalarField& alpha = phase; + + alphaPhiCorrs.set + ( + phase.index(), + new surfaceScalarField + ( + IOobject::groupName("alphaPhiCorr", phase.name()), + fvc::flux(phi_, phase, "div(phi," + alpha.name() + ')') ) ); - surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei]; + surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phase.index()]; - forAll(phases(), phasej) + forAll(phases(), phasei) { - phaseModel& phase2 = phases()[phasej]; + phaseModel& phase2 = phases()[phasei]; volScalarField& alpha2 = phase2; if (&phase2 == &phase) continue; @@ -123,7 +167,7 @@ void Foam::multiphaseSystem::solveAlphas() word phirScheme ( - "div(phir," + alpha2.name() + ',' + alpha1.name() + ')' + "div(phir," + alpha2.name() + ',' + alpha.name() + ')' ); alphaPhiCorr += fvc::flux @@ -144,33 +188,27 @@ void Foam::multiphaseSystem::solveAlphas() alphaPhiCorr, zeroField(), zeroField(), - UniformField(phase.alphaMax()), + min(alphaVoid.primitiveField(), phase.alphaMax())(), zeroField(), true ); } - MULES::limitSum(alphaPhiCorrs); - - volScalarField sumAlpha - ( - IOobject - ( - "sumAlpha", - mesh_.time().timeName(), - mesh_ - ), - mesh_, - dimensionedScalar("sumAlpha", dimless, 0) - ); - - - forAll(phases(), phasei) + // Limit the flux sums, fixing those of the stationary phases + labelHashSet fixedAlphaPhiCorrs; + forAll(stationaryPhases(), stationaryPhasei) { - phaseModel& phase = phases()[phasei]; + fixedAlphaPhiCorrs.insert(stationaryPhases()[stationaryPhasei].index()); + } + MULES::limitSum(alphafs, alphaPhiCorrs, fixedAlphaPhiCorrs); + + // Solve for the moving phase alphas + forAll(movingPhases(), movingPhasei) + { + phaseModel& phase = movingPhases()[movingPhasei]; volScalarField& alpha = phase; - surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei]; + surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()]; alphaPhi += upwind(mesh_, phi_).flux(phase); phase.correctInflowOutflow(alphaPhi); @@ -252,29 +290,47 @@ void Foam::multiphaseSystem::solveAlphas() Su ); - phase.alphaPhi() = alphaPhi; + phase.alphaPhiRef() = alphaPhi; + } - Info<< phase.name() << " volume fraction, min, max = " + // Report the phase fractions and the phase fraction sum + forAll(phases(), phasei) + { + phaseModel& phase = phases()[phasei]; + + Info<< phase.name() << " fraction, min, max = " << phase.weightedAverage(mesh_.V()).value() << ' ' << min(phase).value() << ' ' << max(phase).value() << endl; + } - sumAlpha += phase; + volScalarField sumAlphaMoving + ( + IOobject + ( + "sumAlphaMoving", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("sumAlphaMoving", dimless, 0) + ); + forAll(movingPhases(), movingPhasei) + { + sumAlphaMoving += movingPhases()[movingPhasei]; } Info<< "Phase-sum volume fraction, min, max = " - << sumAlpha.weightedAverage(mesh_.V()).value() - << ' ' << min(sumAlpha).value() - << ' ' << max(sumAlpha).value() + << (sumAlphaMoving + 1 - alphaVoid)().weightedAverage(mesh_.V()).value() + << ' ' << min(sumAlphaMoving + 1 - alphaVoid).value() + << ' ' << max(sumAlphaMoving + 1 - alphaVoid).value() << endl; - // Correct the sum of the phase-fractions to avoid 'drift' - volScalarField sumCorr(1 - sumAlpha); - forAll(phases(), phasei) + // Correct the sum of the phase fractions to avoid drift + forAll(movingPhases(), movingPhasei) { - volScalarField& alpha = phases()[phasei]; - alpha += alpha*sumCorr; + movingPhases()[movingPhasei] *= alphaVoid/sumAlphaMoving; } } @@ -643,9 +699,11 @@ void Foam::multiphaseSystem::solve() forAll(phases(), phasei) { phaseModel& phase = phases()[phasei]; + if (phase.stationary()) continue; + volScalarField& alpha = phase; - phase.alphaPhi() = alphaPhiSums[phasei]/nAlphaSubCycles; + phase.alphaPhiRef() = alphaPhiSums[phasei]/nAlphaSubCycles; // Correct the time index of the field // to correspond to the global time @@ -664,9 +722,11 @@ void Foam::multiphaseSystem::solve() forAll(phases(), phasei) { phaseModel& phase = phases()[phasei]; - phase.alphaRhoPhi() = fvc::interpolate(phase.rho())*phase.alphaPhi(); + if (phase.stationary()) continue; + + phase.alphaRhoPhiRef() = + fvc::interpolate(phase.rho())*phase.alphaPhi(); - // Ensure the phase-fractions are bounded phase.maxMin(0, 1); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H index b133452ea..81b03bb0f 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H @@ -114,9 +114,6 @@ private: //- Return the virtual mass coefficient for phase pair virtual tmp Vm(const phasePairKey& key) const = 0; - //- Return the interfacial mass flow rate for phase pair - virtual tmp dmdt(const phasePairKey& key) const = 0; - protected: @@ -168,77 +165,7 @@ public: // Member Functions using phaseSystem::sigma; - using phaseSystem::dmdt; - - //- Return the momentum transfer matrices for the cell-based algorithm - virtual autoPtr momentumTransfer() const = 0; - - //- Return the momentum transfer matrices for the face-based algiorithm - virtual autoPtr momentumTransferf() const = 0; - - //- Return the implicit force coefficients for the face-based algorithm - virtual Xfer> AFfs() const = 0; - - //- Return the force fluxes for the cell-based algorithm - virtual Xfer> phiFs - ( - const PtrList& rAUs - ) = 0; - - //- Return the force fluxes for the face-based algorithm - virtual Xfer> phiFfs - ( - const PtrList& rAUfs - ) = 0; - - //- Return the force fluxes for the cell-based algorithm - virtual Xfer> phiKdPhis - ( - const PtrList& rAUs - ) const = 0; - - //- Return the force fluxes for the face-based algorithm - virtual Xfer> phiKdPhifs - ( - const PtrList& rAUfs - ) const = 0; - - //- Return the explicit part of the drag force - virtual Xfer> KdUByAs - ( - const PtrList& rAUs - ) const = 0; - - //- Solve the drag system for the new velocities and fluxes - virtual void partialElimination - ( - const PtrList& rAUs - ) = 0; - - //- Solve the drag system for the new fluxes - virtual void partialEliminationf - ( - const PtrList& rAUfs - ) = 0; - - //- Return the flux corrections for the cell-based algorithm - virtual Xfer> ddtCorrByAs - ( - const PtrList& rAUs, - const bool includeVirtualMass = false - ) const = 0; - - //- Return the phase diffusivities divided by the momentum coefficients - virtual const HashPtrTable& DByAfs() const = 0; - - //- Return the heat transfer matrices - virtual autoPtr heatTransfer() const = 0; - - //- Return the mass transfer matrices - virtual autoPtr massTransfer() const = 0; - - //- Return true if there is mass transfer - virtual bool transfersMass() const = 0; + using phaseSystem::dmdts; //- Return the surface tension force tmp surfaceTension(const phaseModel& phase) const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H index cf70a99e1..bb439af80 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H @@ -9,17 +9,17 @@ PtrList UEqns(phases.size()); phaseSystem::momentumTransferTable& momentumTransfer(momentumTransferPtr()); - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; const volScalarField& alpha = phase; const volScalarField& rho = phase.rho(); - volVectorField& U = phase.U(); + volVectorField& U = phase.URef(); UEqns.set ( - phasei, + phase.index(), new fvVectorMatrix ( phase.UEqn() @@ -29,8 +29,8 @@ PtrList UEqns(phases.size()); ) ); - UEqns[phasei].relax(); - fvOptions.constrain(UEqns[phasei]); + UEqns[phase.index()].relax(); + fvOptions.constrain(UEqns[phase.index()]); fvOptions.correct(U); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H index 243e4b9a4..5e488110b 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H @@ -1,8 +1,5 @@ +// Face volume fractions PtrList alphafs(phases.size()); -PtrList rAUs(phases.size()); -PtrList alpharAUfs(phases.size()); - -// Diagonal coefficients forAll(phases, phasei) { phaseModel& phase = phases[phasei]; @@ -10,20 +7,37 @@ forAll(phases, phasei) alphafs.set(phasei, fvc::interpolate(alpha).ptr()); alphafs[phasei].rename("pEqn" + alphafs[phasei].name()); +} + +// Diagonal coefficients +PtrList rAUs(phases.size()); +forAll(fluid.movingPhases(), movingPhasei) +{ + phaseModel& phase = fluid.movingPhases()[movingPhasei]; + const volScalarField& alpha = phase; rAUs.set ( - phasei, + phase.index(), new volScalarField ( IOobject::groupName("rAU", phase.name()), 1.0 /( - UEqns[phasei].A() + UEqns[phase.index()].A() + byDt(max(phase.residualAlpha() - alpha, scalar(0))*phase.rho()) ) ) ); +} +fluid.fillFields("rAU", dimTime/dimDensity, rAUs); + +// Phase diagonal coefficients +PtrList alpharAUfs(phases.size()); +forAll(phases, phasei) +{ + phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; alpharAUfs.set ( @@ -45,6 +59,13 @@ while (pimple.correct()) // Correct p_rgh for consistency with p and the updated densities p_rgh = p - rho*gh; + // Correct fixed-flux BCs to be consistent with the velocity BCs + forAll(fluid.movingPhases(), movingPhasei) + { + phaseModel& phase = fluid.movingPhases()[movingPhasei]; + MRF.correctBoundaryFlux(phase.U(), phase.phiRef()); + } + // Combined buoyancy and force fluxes PtrList phigFs(phases.size()); { @@ -85,17 +106,14 @@ while (pimple.correct()) // Correction force fluxes PtrList ddtCorrByAs(fluid.ddtCorrByAs(rAUs)); - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; const volScalarField& alpha = phase; - // Correct fixed-flux BCs to be consistent with the velocity BCs - MRF.correctBoundaryFlux(phase.U(), phase.phi()); - HbyAs.set ( - phasei, + phase.index(), new volVectorField ( IOobject::groupName("HbyA", phase.name()), @@ -103,31 +121,33 @@ while (pimple.correct()) ) ); - HbyAs[phasei] = - rAUs[phasei] + HbyAs[phase.index()] = + rAUs[phase.index()] *( - UEqns[phasei].H() + UEqns[phase.index()].H() + byDt ( max(phase.residualAlpha() - alpha, scalar(0)) *phase.rho() ) - *phase.U().oldTime() + *phase.U()().oldTime() ); phiHbyAs.set ( - phasei, + phase.index(), new surfaceScalarField ( IOobject::groupName("phiHbyA", phase.name()), - fvc::flux(HbyAs[phasei]) - - phigFs[phasei] - - ddtCorrByAs[phasei] + fvc::flux(HbyAs[phase.index()]) + - phigFs[phase.index()] + - ddtCorrByAs[phase.index()] ) ); } } + fluid.fillFields("HbyA", dimVelocity, HbyAs); + fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs); // Add explicit drag forces and fluxes if not doing partial elimination if (!partialElimination) @@ -211,7 +231,8 @@ while (pimple.correct()) forAll(phases, phasei) { phaseModel& phase = phases[phasei]; - phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField(); + phib += + alphafs[phasei].boundaryField()*phase.phi()().boundaryField(); } setSnGrad @@ -225,11 +246,12 @@ while (pimple.correct()) // Compressible pressure equations PtrList pEqnComps(phases.size()); + PtrList dmdts(fluid.dmdts()); forAll(phases, phasei) { phaseModel& phase = phases[phasei]; const volScalarField& alpha = phase; - volScalarField& rho = phase.thermo().rho(); + volScalarField& rho = phase.thermoRef().rho(); if (phase.compressible()) { @@ -246,7 +268,7 @@ while (pimple.correct()) phasei, ( ( - phase.continuityError() + fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi()) - fvc::Sp ( fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), @@ -274,7 +296,7 @@ while (pimple.correct()) phasei, ( ( - phase.continuityError() + fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi()) - fvc::Sp ( (fvc::ddt(alpha) + fvc::div(phase.alphaPhi())), @@ -287,27 +309,38 @@ while (pimple.correct()) ); } } - else - { - pEqnComps.set - ( - phasei, - fvm::Su(-(fvOptions(alpha, rho)&rho)/rho, p_rgh).ptr() - ); - } - if (fluid.transfersMass()) + // Add option sources + if (fvOptions.appliesToField(rho.name())) { + tmp optEqn = fvOptions(alpha, rho); if (pEqnComps.set(phasei)) { - pEqnComps[phasei] -= fluid.dmdt(phase)/rho; + pEqnComps[phasei] -= (optEqn&rho)/rho; } else { pEqnComps.set ( phasei, - fvm::Su(-fluid.dmdt(phase)/rho, p_rgh) + fvm::Su(- (optEqn&rho)/rho, p_rgh).ptr() + ); + } + } + + // Add mass transfer + if (dmdts.set(phasei)) + { + if (pEqnComps.set(phasei)) + { + pEqnComps[phasei] -= dmdts[phasei]/rho; + } + else + { + pEqnComps.set + ( + phasei, + fvm::Su(- dmdts[phasei]/rho, p_rgh) ); } } @@ -351,16 +384,18 @@ while (pimple.correct()) surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf); - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; - phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp; + phase.phiRef() = + phiHbyAs[phase.index()] + + alpharAUfs[phase.index()]*mSfGradp; // Set the phase dilatation rates - if (pEqnComps.set(phasei)) + if (pEqnComps.set(phase.index())) { - phase.divU(-pEqnComps[phasei] & p_rgh); + phase.divU(-pEqnComps[phase.index()] & p_rgh); } } @@ -369,16 +404,16 @@ while (pimple.correct()) mSfGradp = pEqnIncomp.flux()/rAUf; - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; - phase.U() = - HbyAs[phasei] + phase.URef() = + HbyAs[phase.index()] + fvc::reconstruct ( - alpharAUfs[phasei]*mSfGradp - - phigFs[phasei] + alpharAUfs[phase.index()]*mSfGradp + - phigFs[phase.index()] ); } @@ -387,12 +422,12 @@ while (pimple.correct()) fluid.partialElimination(rAUs); } - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; - phase.U().correctBoundaryConditions(); - fvOptions.correct(phase.U()); + phase.URef().correctBoundaryConditions(); + fvOptions.correct(phase.URef()); } } } @@ -407,7 +442,7 @@ while (pimple.correct()) forAll(phases, phasei) { phaseModel& phase = phases[phasei]; - phase.thermo().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0); + phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0); } // Correct p_rgh for consistency with p and the updated densities diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H index 04b88edb5..cf72498e9 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H @@ -12,17 +12,17 @@ PtrList UEqns(phases.size()); phaseSystem::momentumTransferTable& momentumTransfer(momentumTransferPtr()); - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; const volScalarField& alpha = phase; - volScalarField& rho = phase.thermo().rho(); - volVectorField& U = phase.U(); + const volScalarField& rho = phase.rho(); + volVectorField& U = phase.URef(); UEqns.set ( - phasei, + phase.index(), new fvVectorMatrix ( phase.UfEqn() @@ -32,8 +32,8 @@ PtrList UEqns(phases.size()); ) ); - UEqns[phasei].relax(); - fvOptions.constrain(UEqns[phasei]); + UEqns[phase.index()].relax(); + fvOptions.constrain(UEqns[phase.index()]); U.correctBoundaryConditions(); fvOptions.correct(U); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H index cf3ed440e..98da1ded9 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H @@ -1,56 +1,68 @@ +// Face volume fractions PtrList alphafs(phases.size()); PtrList alphaRho0fs(phases.size()); -PtrList rAUfs(phases.size()); -PtrList alpharAUfs(phases.size()); +forAll(phases, phasei) +{ + phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + + alphafs.set(phasei, fvc::interpolate(alpha).ptr()); + alphafs[phasei].rename("pEqn" + alphafs[phasei].name()); + + alphaRho0fs.set + ( + phasei, + ( + fvc::interpolate + ( + max(alpha.oldTime(), phase.residualAlpha()) + *phase.rho()().oldTime() + ) + ).ptr() + ); +} // Diagonal coefficients +PtrList rAUfs(phases.size()); { PtrList AFfs(fluid.AFfs()); - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; - const volScalarField& alpha = phase; - - alphafs.set(phasei, fvc::interpolate(alpha).ptr()); - alphafs[phasei].rename("pEqn" + alphafs[phasei].name()); - - alphaRho0fs.set - ( - phasei, - ( - fvc::interpolate - ( - max(alpha.oldTime(), phase.residualAlpha()) - *phase.rho()().oldTime() - ) - ).ptr() - ); + phaseModel& phase = fluid.movingPhases()[movingPhasei]; rAUfs.set ( - phasei, + phase.index(), new surfaceScalarField ( IOobject::groupName("rAUf", phase.name()), 1.0 /( - byDt(alphaRho0fs[phasei]) - + fvc::interpolate(UEqns[phasei].A()) - + AFfs[phasei] + byDt(alphaRho0fs[phase.index()]) + + fvc::interpolate(UEqns[phase.index()].A()) + + AFfs[phase.index()] ) ) ); - - alpharAUfs.set - ( - phasei, - ( - max(alphafs[phasei], phase.residualAlpha())*rAUfs[phasei] - ).ptr() - ); } } +fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs); + +// Phase diagonal coefficients +PtrList alpharAUfs(phases.size()); +forAll(phases, phasei) +{ + phaseModel& phase = phases[phasei]; + alpharAUfs.set + ( + phase.index(), + ( + max(alphafs[phase.index()], phase.residualAlpha()) + *rAUfs[phase.index()] + ).ptr() + ); +} // Explicit force fluxes PtrList phiFfs(fluid.phiFfs(rAUfs)); @@ -63,12 +75,11 @@ while (pimple.correct()) // Correct p_rgh for consistency with p and the updated densities p_rgh = p - rho*gh; - forAll(phases, phasei) + // Correct fixed-flux BCs to be consistent with the velocity BCs + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; - - // Correct fixed-flux BCs to be consistent with the velocity BCs - MRF.correctBoundaryFlux(phase.U(), phase.phi()); + phaseModel& phase = fluid.movingPhases()[movingPhasei]; + MRF.correctBoundaryFlux(phase.U(), phase.phiRef()); } // Combined buoyancy and force fluxes @@ -106,27 +117,27 @@ while (pimple.correct()) // Predicted fluxes for each phase PtrList phiHbyAs(phases.size()); - - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; phiHbyAs.set ( - phasei, + phase.index(), new surfaceScalarField ( IOobject::groupName("phiHbyA", phase.name()), - rAUfs[phasei] + rAUfs[phase.index()] *( - fvc::flux(UEqns[phasei].H()) - + alphaRho0fs[phasei] - *byDt(MRF.absolute(phase.phi().oldTime())) + fvc::flux(UEqns[phase.index()].H()) + + alphaRho0fs[phase.index()] + *byDt(MRF.absolute(phase.phi()().oldTime())) ) - - phigFs[phasei] + - phigFs[phase.index()] ) ); } + fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs); // Add explicit drag forces and fluxes if not doing partial elimination if (!partialElimination) @@ -205,7 +216,8 @@ while (pimple.correct()) forAll(phases, phasei) { phaseModel& phase = phases[phasei]; - phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField(); + phib += + alphafs[phasei].boundaryField()*phase.phi()().boundaryField(); } setSnGrad @@ -219,11 +231,12 @@ while (pimple.correct()) // Compressible pressure equations PtrList pEqnComps(phases.size()); + PtrList dmdts(fluid.dmdts()); forAll(phases, phasei) { phaseModel& phase = phases[phasei]; const volScalarField& alpha = phase; - volScalarField& rho = phase.thermo().rho(); + volScalarField& rho = phase.thermoRef().rho(); if (phase.compressible()) { @@ -240,7 +253,7 @@ while (pimple.correct()) phasei, ( ( - phase.continuityError() + fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi()) - fvc::Sp ( fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), @@ -263,6 +276,7 @@ while (pimple.correct()) ( pEqnComps[phasei].faceFluxCorrectionPtr() ); + pEqnComps[phasei].relax(); } else @@ -272,7 +286,7 @@ while (pimple.correct()) phasei, ( ( - phase.continuityError() + fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi()) - fvc::Sp ( (fvc::ddt(alpha) + fvc::div(phase.alphaPhi())), @@ -285,27 +299,36 @@ while (pimple.correct()) ); } } - else - { - pEqnComps.set - ( - phasei, - fvm::Su(-(fvOptions(alpha, rho)&rho)/rho, p_rgh).ptr() - ); - } - if (fluid.transfersMass()) + if (fvOptions.appliesToField(rho.name())) { + tmp optEqn = fvOptions(alpha, rho); if (pEqnComps.set(phasei)) { - pEqnComps[phasei] -= fluid.dmdt(phase)/rho; + pEqnComps[phasei] -= (optEqn&rho)/rho; } else { pEqnComps.set ( phasei, - fvm::Su(-fluid.dmdt(phase)/rho, p_rgh) + fvm::Su(- (optEqn&rho)/rho, p_rgh).ptr() + ); + } + } + + if (dmdts.set(phasei)) + { + if (pEqnComps.set(phasei)) + { + pEqnComps[phasei] -= dmdts[phasei]/rho; + } + else + { + pEqnComps.set + ( + phasei, + fvm::Su(- dmdts[phasei]/rho, p_rgh) ); } } @@ -349,16 +372,18 @@ while (pimple.correct()) surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf); - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; - phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp; + phase.phiRef() = + phiHbyAs[phase.index()] + + alpharAUfs[phase.index()]*mSfGradp; // Set the phase dilatation rates - if (pEqnComps.set(phasei)) + if (pEqnComps.set(phase.index())) { - phase.divU(-pEqnComps[phasei] & p_rgh); + phase.divU(-pEqnComps[phase.index()] & p_rgh); } } @@ -372,13 +397,13 @@ while (pimple.correct()) mSfGradp = pEqnIncomp.flux()/rAUf; - forAll(phases, phasei) + forAll(fluid.movingPhases(), movingPhasei) { - phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid.movingPhases()[movingPhasei]; - phase.U() = fvc::reconstruct(MRF.absolute(phase.phi())); - phase.U().correctBoundaryConditions(); - fvOptions.correct(phase.U()); + phase.URef() = fvc::reconstruct(MRF.absolute(phase.phi())); + phase.URef().correctBoundaryConditions(); + fvOptions.correct(phase.URef()); } } } @@ -393,7 +418,7 @@ while (pimple.correct()) forAll(phases, phasei) { phaseModel& phase = phases[phasei]; - phase.thermo().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0); + phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0); } // Correct p_rgh for consistency with p and the updated densities diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H index 246bea1fb..5dfc536aa 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H @@ -8,46 +8,38 @@ for (int Ecorr=0; Ecorr E1Eqn(phase1.heEqn()); + fvScalarMatrix E1Eqn + ( + phase1.heEqn() + == + *heatTransfer[phase1.name()] + + alpha1*rho1*(U1&g) + + fvOptions(alpha1, rho1, thermo1.he()) + ); - if (E1Eqn.valid()) - { - E1Eqn = - ( - E1Eqn - == - *heatTransfer[phase1.name()] - + alpha1*rho1*(U1&g) - + fvOptions(alpha1, rho1, phase1.thermo().he()) - ); - - E1Eqn->relax(); - fvOptions.constrain(E1Eqn.ref()); - E1Eqn->solve(); - fvOptions.correct(phase1.thermo().he()); - } + E1Eqn.relax(); + fvOptions.constrain(E1Eqn); + E1Eqn.solve(); + fvOptions.correct(thermo1.he()); } + if (!phase2.isothermal()) { - tmp E2Eqn(phase2.heEqn()); + fvScalarMatrix E2Eqn + ( + phase2.heEqn() + == + *heatTransfer[phase2.name()] + + alpha2*rho2*(U2&g) + + fvOptions(alpha2, rho2, phase2.thermoRef().he()) + ); - if (E2Eqn.valid()) - { - E2Eqn = - ( - E2Eqn - == - *heatTransfer[phase2.name()] - + alpha2*rho2*(U2&g) - + fvOptions(alpha2, rho2, phase2.thermo().he()) - ); - - E2Eqn->relax(); - fvOptions.constrain(E2Eqn.ref()); - E2Eqn->solve(); - fvOptions.correct(phase2.thermo().he()); - } + E2Eqn.relax(); + fvOptions.constrain(E2Eqn); + E2Eqn.solve(); + fvOptions.correct(thermo2.he()); } fluid.correctThermo(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/YEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/YEqns.H index 176f491f3..4c9ecb94a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/YEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/YEqns.H @@ -5,44 +5,41 @@ phaseSystem::massTransferTable& massTransfer(massTransferPtr()); - PtrList& Y1 = phase1.Y(); - PtrList& Y2 = phase2.Y(); - - forAll(Y1, i) + if (!phase1.pure()) { - tmp Y1iEqn(phase1.YiEqn(Y1[i])); + UPtrList& Y1 = phase1.YActiveRef(); - if (Y1iEqn.valid()) + forAll(Y1, i) { - Y1iEqn = + fvScalarMatrix Y1iEqn ( - Y1iEqn + phase1.YiEqn(Y1[i]) == *massTransfer[Y1[i].name()] + fvOptions(alpha1, rho1, Y1[i]) ); - Y1iEqn->relax(); - Y1iEqn->solve(mesh.solver("Yi")); + Y1iEqn.relax(); + Y1iEqn.solve(mesh.solver("Yi")); } } - forAll(Y2, i) + if (!phase2.pure()) { - tmp Y2iEqn(phase2.YiEqn(Y2[i])); + UPtrList& Y2 = phase2.YActiveRef(); - if (Y2iEqn.valid()) + forAll(Y2, i) { - Y2iEqn = + fvScalarMatrix Y2iEqn ( - Y2iEqn + phase2.YiEqn(Y2[i]) == *massTransfer[Y2[i].name()] + fvOptions(alpha2, rho2, Y2[i]) ); - Y2iEqn->relax(); - Y2iEqn->solve(mesh.solver("Yi")); + Y2iEqn.relax(); + Y2iEqn.solve(mesh.solver("Yi")); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFieldRefs.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFieldRefs.H index 193e5388d..edad54fa1 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFieldRefs.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFieldRefs.H @@ -1,20 +1,21 @@ phaseModel& phase1 = fluid.phase1(); phaseModel& phase2 = fluid.phase2(); -volScalarField& alpha1 = phase1; -volScalarField& alpha2 = phase2; +const volScalarField& alpha1 = phase1; +const volScalarField& alpha2 = phase2; -volVectorField& U1 = phase1.U(); -surfaceScalarField& phi1 = phase1.phi(); -surfaceScalarField& alphaPhi1 = phase1.alphaPhi(); +volVectorField& U1 = phase1.URef(); +surfaceScalarField& phi1 = phase1.phiRef(); +const surfaceScalarField& alphaPhi1 = phase1.alphaPhi(); + +volVectorField& U2 = phase2.URef(); +surfaceScalarField& phi2 = phase2.phiRef(); +const surfaceScalarField& alphaPhi2 = phase2.alphaPhi(); -volVectorField& U2 = phase2.U(); -surfaceScalarField& phi2 = phase2.phi(); -surfaceScalarField& alphaPhi2 = phase2.alphaPhi(); surfaceScalarField& phi = fluid.phi(); -rhoThermo& thermo1 = phase1.thermo(); -rhoThermo& thermo2 = phase2.thermo(); +rhoThermo& thermo1 = phase1.thermoRef(); +rhoThermo& thermo2 = phase2.thermoRef(); volScalarField& rho1 = thermo1.rho(); const volScalarField& psi1 = thermo1.psi(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFields.H index 697269ff4..df7c20a79 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFields.H @@ -19,7 +19,7 @@ dimensionedScalar pMin #include "gh.H" -volScalarField& p = fluid.phase1().thermo().p(); +volScalarField& p = fluid.phase1().thermoRef().p(); Info<< "Reading field p_rgh\n" << endl; volScalarField p_rgh diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H index 783f1bcc7..d94526e01 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H @@ -184,7 +184,7 @@ while (pimple.correct()) pEqnComp1 = ( - phase1.continuityError() + fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi()) - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) )/rho1 + correction @@ -202,17 +202,12 @@ while (pimple.correct()) { pEqnComp1 = ( - phase1.continuityError() + fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi()) - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) )/rho1 + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)); } } - else - { - pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh); - } - if (phase2.compressible()) { if (pimple.transonic()) @@ -225,7 +220,7 @@ while (pimple.correct()) pEqnComp2 = ( - phase2.continuityError() + fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi()) - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) )/rho2 + correction @@ -243,35 +238,65 @@ while (pimple.correct()) { pEqnComp2 = ( - phase2.continuityError() + fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi()) - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) )/rho2 + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)); } } - else + + // Add option sources { - pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh); + if (fvOptions.appliesToField(rho1.name())) + { + tmp optEqn1 = fvOptions(alpha1, rho1); + if (pEqnComp1.valid()) + { + pEqnComp1.ref() -= (optEqn1 & rho1)/rho1; + } + else + { + pEqnComp1 = fvm::Su(- (optEqn1 & rho1)/rho1, p_rgh); + } + } + if (fvOptions.appliesToField(rho2.name())) + { + tmp optEqn2 = fvOptions(alpha2, rho2); + if (pEqnComp2.valid()) + { + pEqnComp2.ref() -= (optEqn2 & rho2)/rho2; + } + else + { + pEqnComp2 = fvm::Su(- (optEqn2 & rho2)/rho2, p_rgh); + } + } } - if (fluid.transfersMass()) + // Add mass transfer { - if (pEqnComp1.valid()) + PtrList dmdts(fluid.dmdts()); + if (dmdts.set(0)) { - pEqnComp1.ref() -= fluid.dmdt(phase1)/rho1; + if (pEqnComp1.valid()) + { + pEqnComp1.ref() -= dmdts[0]/rho1; + } + else + { + pEqnComp1 = fvm::Su(- dmdts[0]/rho1, p_rgh); + } } - else + if (dmdts.set(1)) { - pEqnComp1 = fvm::Su(-fluid.dmdt(phase1)/rho1, p_rgh); - } - - if (pEqnComp2.valid()) - { - pEqnComp2.ref() -= fluid.dmdt(phase2)/rho2; - } - else - { - pEqnComp2 = fvm::Su(-fluid.dmdt(phase2)/rho2, p_rgh); + if (pEqnComp2.valid()) + { + pEqnComp2.ref() -= dmdts[1]/rho2; + } + else + { + pEqnComp2 = fvm::Su(- dmdts[1]/rho2, p_rgh); + } } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H index 11032d514..c5d0e7869 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H @@ -196,7 +196,7 @@ while (pimple.correct()) pEqnComp1 = ( - phase1.continuityError() + fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi()) - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) )/rho1 + correction @@ -214,17 +214,12 @@ while (pimple.correct()) { pEqnComp1 = ( - phase1.continuityError() + fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi()) - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) )/rho1 + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)); } } - else - { - pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh); - } - if (phase2.compressible()) { if (pimple.transonic()) @@ -237,7 +232,7 @@ while (pimple.correct()) pEqnComp2 = ( - phase2.continuityError() + fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi()) - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) )/rho2 + correction @@ -255,35 +250,65 @@ while (pimple.correct()) { pEqnComp2 = ( - phase2.continuityError() + fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi()) - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) )/rho2 + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)); } } - else + + // Add option sources { - pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh); + if (fvOptions.appliesToField(rho1.name())) + { + tmp optEqn1 = fvOptions(alpha1, rho1); + if (pEqnComp1.valid()) + { + pEqnComp1.ref() -= (optEqn1 & rho1)/rho1; + } + else + { + pEqnComp1 = fvm::Su(- (optEqn1 & rho1)/rho1, p_rgh); + } + } + if (fvOptions.appliesToField(rho2.name())) + { + tmp optEqn2 = fvOptions(alpha2, rho2); + if (pEqnComp2.valid()) + { + pEqnComp2.ref() -= (optEqn2 & rho2)/rho2; + } + else + { + pEqnComp2 = fvm::Su(- (optEqn2 & rho2)/rho2, p_rgh); + } + } } - if (fluid.transfersMass()) + // Add mass transfer { - if (pEqnComp1.valid()) + PtrList dmdts(fluid.dmdts()); + if (dmdts.set(0)) { - pEqnComp1.ref() -= fluid.dmdt(phase1)/rho1; + if (pEqnComp1.valid()) + { + pEqnComp1.ref() -= dmdts[0]/rho1; + } + else + { + pEqnComp1 = fvm::Su(- dmdts[0]/rho1, p_rgh); + } } - else + if (dmdts.set(1)) { - pEqnComp1 = fvm::Su(-fluid.dmdt(phase1)/rho1, p_rgh); - } - - if (pEqnComp2.valid()) - { - pEqnComp2.ref() -= fluid.dmdt(phase2)/rho2; - } - else - { - pEqnComp2 = fvm::Su(-fluid.dmdt(phase2)/rho2, p_rgh); + if (pEqnComp2.valid()) + { + pEqnComp2.ref() -= dmdts[1]/rho2; + } + else + { + pEqnComp2 = fvm::Su(- dmdts[1]/rho2, p_rgh); + } } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index 796df17cc..bbb2b75bb 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -29,6 +29,7 @@ License #include "MULES.H" #include "subCycle.H" +#include "UniformField.H" #include "fvcDdt.H" #include "fvcDiv.H" @@ -280,15 +281,15 @@ void Foam::twoPhaseSystem::solve() if (alphaSubCycle.index() == 1) { - phase1_.alphaPhi() = alphaPhi10; + phase1_.alphaPhiRef() = alphaPhi10; } else { - phase1_.alphaPhi() += alphaPhi10; + phase1_.alphaPhiRef() += alphaPhi10; } } - phase1_.alphaPhi() /= nAlphaSubCycles; + phase1_.alphaPhiRef() /= nAlphaSubCycles; } else { @@ -304,7 +305,7 @@ void Foam::twoPhaseSystem::solve() zeroField() ); - phase1_.alphaPhi() = alphaPhi1; + phase1_.alphaPhiRef() = alphaPhi1; } if (alphaDbyA.valid()) @@ -318,15 +319,15 @@ void Foam::twoPhaseSystem::solve() alpha1Eqn.relax(); alpha1Eqn.solve(); - phase1_.alphaPhi() += alpha1Eqn.flux(); + phase1_.alphaPhiRef() += alpha1Eqn.flux(); } - phase1_.alphaRhoPhi() = + phase1_.alphaRhoPhiRef() = fvc::interpolate(phase1_.rho())*phase1_.alphaPhi(); - phase2_.alphaPhi() = phi_ - phase1_.alphaPhi(); - phase2_.correctInflowOutflow(phase2_.alphaPhi()); - phase2_.alphaRhoPhi() = + phase2_.alphaPhiRef() = phi_ - phase1_.alphaPhi(); + phase2_.correctInflowOutflow(phase2_.alphaPhiRef()); + phase2_.alphaRhoPhiRef() = fvc::interpolate(phase2_.rho())*phase2_.alphaPhi(); Info<< alpha1.name() << " volume fraction = " diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H index aaade2c13..eb4315a63 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H @@ -66,9 +66,6 @@ private: //- Return the virtual mass coefficient for phase pair virtual tmp Vm(const phasePairKey& key) const = 0; - //- Return the interfacial mass flow rate for phase pair - virtual tmp dmdt(const phasePairKey& key) const = 0; - protected: @@ -126,90 +123,23 @@ public: // Member Functions using phaseSystem::sigma; - using phaseSystem::dmdt; + using phaseSystem::dmdts; - //- Constant access phase model 1 + //- Return phase model 1 const phaseModel& phase1() const; //- Access phase model 1 phaseModel& phase1(); - //- Constant access phase model 2 + //- Return phase model 2 const phaseModel& phase2() const; //- Access phase model 2 phaseModel& phase2(); - //- Constant access the phase not given as an argument + //- Return the phase not given as an argument const phaseModel& otherPhase(const phaseModel& phase) const; - //- Return the momentum transfer matrices for the cell-based algorithm - virtual autoPtr momentumTransfer() const = 0; - - //- Return the momentum transfer matrices for the face-based algorithm - virtual autoPtr momentumTransferf() const = 0; - - //- Return the implicit force coefficients for the face-based algorithm - virtual Xfer> AFfs() const = 0; - - //- Return the force fluxes for the cell-based algorithm - virtual Xfer> phiFs - ( - const PtrList& rAUs - ) = 0; - - //- Return the force fluxes for the face-based algorithm - virtual Xfer> phiFfs - ( - const PtrList& rAUfs - ) = 0; - - //- Return the force fluxes for the cell-based algorithm - virtual Xfer> phiKdPhis - ( - const PtrList& rAUs - ) const = 0; - - //- Return the force fluxes for the face-based algorithm - virtual Xfer> phiKdPhifs - ( - const PtrList& rAUfs - ) const = 0; - - //- Return the explicit part of the drag force - virtual Xfer> KdUByAs - ( - const PtrList& rAUs - ) const = 0; - - //- Solve the drag system for the new velocities and fluxes - virtual void partialElimination - ( - const PtrList& rAUs - ) = 0; - - //- Solve the drag system for the new fluxes - virtual void partialEliminationf - ( - const PtrList& rAUfs - ) = 0; - - //- Return the flux corrections for the cell-based algorithm - virtual Xfer> ddtCorrByAs - ( - const PtrList& rAUs, - const bool includeVirtualMass = false - ) const = 0; - - //- Return the phase diffusivities divided by the momentum coefficients - virtual const HashPtrTable& DByAfs() const = 0; - - //- Return the heat transfer matrices - virtual autoPtr heatTransfer() const = 0; - - //- Return the mass transfer matrices - virtual autoPtr massTransfer() const = 0; - //- Return the surface tension coefficient tmp sigma() const; @@ -222,9 +152,6 @@ public: //- Return the virtual mass coefficient tmp Vm() const; - //- Return true if there is mass transfer - virtual bool transfersMass() const = 0; - //- Solve for the phase fractions virtual void solve(); }; diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.air b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.air new file mode 100644 index 000000000..19718b298 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.air @@ -0,0 +1,36 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object T.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 300; + +boundaryField +{ + top + { + type inletOutlet; + phi phi.air; + inletValue $internalField; + value $internalField; + } + walls + { + type zeroGradient; + } +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.solid b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.solid new file mode 100644 index 000000000..a768f5bd4 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.solid @@ -0,0 +1,33 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object T.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 300; + +boundaryField +{ + top + { + type zeroGradient; + } + walls + { + type zeroGradient; + } +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.water b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.water new file mode 100644 index 000000000..ca4a14aab --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/T.water @@ -0,0 +1,36 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object T.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 300; + +boundaryField +{ + top + { + type inletOutlet; + phi phi.water; + inletValue $internalField; + value $internalField; + } + walls + { + type zeroGradient; + } +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/U.air b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/U.air new file mode 100644 index 000000000..42f7b6477 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/U.air @@ -0,0 +1,35 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format binary; + class volVectorField; + object U.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + top + { + type pressureInletOutletVelocity; + phi phi.air; + value $internalField; + } + walls + { + type noSlip; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/U.water b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/U.water new file mode 100644 index 000000000..8c9544420 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/U.water @@ -0,0 +1,35 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format binary; + class volVectorField; + object U.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + top + { + type pressureInletOutletVelocity; + phi phi.water; + value $internalField; + } + walls + { + type noSlip; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.air.orig b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.air.orig new file mode 100644 index 000000000..96eb2aa71 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.air.orig @@ -0,0 +1,36 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object alpha.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + top + { + type inletOutlet; + phi phi.air; + inletValue $internalField; + value $internalField; + } + walls + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.solid.orig b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.solid.orig new file mode 100644 index 000000000..5d8f456c7 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.solid.orig @@ -0,0 +1,33 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object alpha.solid; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + top + { + type zeroGradient; + } + walls + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.water.orig b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.water.orig new file mode 100644 index 000000000..d0a28e019 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/alpha.water.orig @@ -0,0 +1,36 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object alpha.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + top + { + type inletOutlet; + phi phi.water; + value $internalField; + inletValue $internalField; + } + walls + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/p b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/p new file mode 100644 index 000000000..0430fdf9c --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/p @@ -0,0 +1,35 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform 1e5; + +boundaryField +{ + top + { + type calculated; + value $internalField; + } + walls + { + type calculated; + value $internalField; + } +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/p_rgh b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/p_rgh new file mode 100644 index 000000000..229f70b87 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/0/p_rgh @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p_rgh; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform 1e5; + +boundaryField +{ + top + { + type prghTotalPressure; + p0 $internalField; + U U.air; + phi phi.air; + rho thermo:rho.air; + value $internalField; + } + walls + { + type fixedFluxPressure; + value $internalField; + } +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/Allrun b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/Allrun new file mode 100755 index 000000000..ed61ac543 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/Allrun @@ -0,0 +1,11 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions + +runApplication blockMesh +runApplication setFields +runApplication $(getApplication) + +#------------------------------------------------------------------------------ diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/fvOptions b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/fvOptions new file mode 100644 index 000000000..276b05d1f --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/fvOptions @@ -0,0 +1,46 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object fvOptions; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +injector1 +{ + timeStart 0; + duration 1e6; + selectionMode points; + points + ( + (0.075 0.925 0.05) + ); +} + +options +{ + massSource1 + { + type scalarSemiImplicitSource; + + $injector1; + + volumeMode absolute; + injectionRateSuSp + { + thermo:rho.water (4.0e-1 0); // kg/s + } + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/g b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/g new file mode 100644 index 000000000..0cc222ca3 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/g @@ -0,0 +1,22 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value (0 -9.81 0); + + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/phaseProperties b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/phaseProperties new file mode 100644 index 000000000..301afa9ae --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/phaseProperties @@ -0,0 +1,208 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object phaseProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +type heatAndMomentumTransferMultiphaseSystem; + +phases (air water solid); + +air +{ + type pureIsothermalPhaseModel; + diameterModel isothermal; + isothermalCoeffs + { + d0 3e-3; + p0 1e5; + } + + residualAlpha 1e-6; +} + +water +{ + type pureIsothermalPhaseModel; + diameterModel constant; + constantCoeffs + { + d 1e-4; + } + + residualAlpha 1e-6; +} + +solid +{ + type pureStationaryIsothermalPhaseModel; + diameterModel constant; + constantCoeffs + { + d 1e-3; + } + + residualAlpha 1e-6; +} + +blending +{ + default + { + type linear; + minFullyContinuousAlpha.air 0.7; + minPartlyContinuousAlpha.air 0.3; + minFullyContinuousAlpha.water 0.7; + minPartlyContinuousAlpha.water 0.3; + minFullyContinuousAlpha.solid 0; + minPartlyContinuousAlpha.solid 0; + } +} + +surfaceTension +( + (air and water) + { + type constant; + sigma 0.07; + } + + (air and solid) + { + type constant; + sigma 0; + } + + (solid and water) + { + type constant; + sigma 0; + } +); + +aspectRatio +( + (air in water) + { + type constant; + E0 1.0; + } + + (water in air) + { + type constant; + E0 1.0; + } + + (air in solid) + { + type constant; + E0 1.0; + } + + (solid in air) + { + type constant; + E0 1.0; + } + + (water in solid) + { + type constant; + E0 1.0; + } + + (solid in water) + { + type constant; + E0 1.0; + } +); + +drag +( + (air in water) + { + type SchillerNaumann; + residualRe 1e-3; + swarmCorrection + { + type none; + } + } + + (water in air) + { + type SchillerNaumann; + residualRe 1e-3; + swarmCorrection + { + type none; + } + } + + (solid in air) + { + type Ergun; + residualRe 1e-3; + swarmCorrection + { + type none; + } + } + + (solid in water) + { + type Ergun; + residualRe 1e-3; + swarmCorrection + { + type none; + } + } +); + +virtualMass +( + (air in water) + { + type constantCoefficient; + Cvm 0.5; + } + + (water in air) + { + type constantCoefficient; + Cvm 0.5; + } +); + +heatTransfer +(); + +lift +(); + +wallLubrication +(); + +turbulentDispersion +(); + +interfaceCompression +(); + +pMin 10000; + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.air b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.air new file mode 100644 index 000000000..eb9bb8ab7 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.air @@ -0,0 +1,48 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object thermophysicalProperties.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState perfectGas; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + molWeight 28.9; + } + thermodynamics + { + Cp 1007; + Hf 0; + } + transport + { + mu 1.84e-05; + Pr 0.7; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.solid b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.solid new file mode 100644 index 000000000..19a2516e6 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.solid @@ -0,0 +1,52 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object thermophysicalProperties.solid; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState rhoConst; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + molWeight 100; + } + equationOfState + { + rho 2500; + } + thermodynamics + { + Cp 6000; + Hf 0; + } + transport + { + mu 0; + Pr 1; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.water b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.water new file mode 100644 index 000000000..31d8cf058 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/thermophysicalProperties.water @@ -0,0 +1,53 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object thermophysicalProperties.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState perfectFluid; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + molWeight 18; + } + equationOfState + { + R 3000; + rho0 1027; + } + thermodynamics + { + Cp 4195; + Hf 0; + } + transport + { + mu 3.645e-4; + Pr 2.289; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/turbulenceProperties.air b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/turbulenceProperties.air new file mode 100644 index 000000000..1296429b7 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/turbulenceProperties.air @@ -0,0 +1,20 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/turbulenceProperties.water b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/turbulenceProperties.water new file mode 100644 index 000000000..7f0d75a82 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/constant/turbulenceProperties.water @@ -0,0 +1,20 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/blockMeshDict b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/blockMeshDict new file mode 100644 index 000000000..1ef1e3835 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/blockMeshDict @@ -0,0 +1,58 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 1; + +vertices +( + (0 0 0) + (0.15 0 0) + (0.15 1 0) + (0 1 0) + (0 0 0.1) + (0.15 0 0.1) + (0.15 1 0.1) + (0 1 0.1) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (25 75 1) simpleGrading (1 1 1) +); + +edges +( +); + +patches +( + patch top + ( + (3 7 6 2) + ) + wall walls + ( + (1 5 4 0) + (0 4 7 3) + (2 6 5 1) + ) +); + +mergePatchPairs +( +); + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/controlDict b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/controlDict new file mode 100644 index 000000000..e94cda13c --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/controlDict @@ -0,0 +1,54 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application reactingMultiphaseEulerFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 9; + +deltaT 0.002; + +writeControl runTime; + +writeInterval 0.1; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression uncompressed; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +adjustTimeStep no; + +maxCo 0.5; + +maxDeltaT 1; + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/fvSchemes b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/fvSchemes new file mode 100644 index 000000000..3356dda13 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/fvSchemes @@ -0,0 +1,61 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; +} + +divSchemes +{ + default none; + + "div\(phi,alpha.*\)" Gauss vanLeer; + "div\(phir,alpha.*,alpha.*\)" Gauss vanLeer; + + "div\(alphaRhoPhi.*,U.*\)" Gauss limitedLinearV 1; + "div\(phi.*,U.*\)" Gauss limitedLinearV 1; + + "div\(alphaRhoPhi.*,(h|e).*\)" Gauss limitedLinear 1; + "div\(alphaRhoPhi.*,K.*\)" Gauss limitedLinear 1; + "div\(alphaPhi.*,p\)" Gauss limitedLinear 1; + + "div\(\(\(\(alpha.*\*thermo:rho.*\)*nuEff.*\)\*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear; +} + +laplacianSchemes +{ + default Gauss linear uncorrected; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default uncorrected; +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/fvSolution b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/fvSolution new file mode 100644 index 000000000..5ea62ead3 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/fvSolution @@ -0,0 +1,75 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + "alpha.*" + { + nAlphaCorr 1; + nAlphaSubCycles 2; + } + + p_rgh + { + solver GAMG; + smoother DIC; + tolerance 1e-8; + relTol 0; + } + + p_rghFinal + { + $p_rgh; + relTol 0; + } + + "U.*" + { + solver smoothSolver; + smoother symGaussSeidel; + tolerance 1e-5; + relTol 0; + minIter 1; + } + + "e.*" + { + solver smoothSolver; + smoother symGaussSeidel; + tolerance 1e-8; + relTol 0; + minIter 1; + } +} + +PIMPLE +{ + nOuterCorrectors 3; + nCorrectors 1; + nNonOrthogonalCorrectors 0; +} + +relaxationFactors +{ + equations + { + ".*" 1; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/setFieldsDict b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/setFieldsDict new file mode 100644 index 000000000..d6f271df0 --- /dev/null +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bed/system/setFieldsDict @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object setFieldsDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +defaultFieldValues +( + volScalarFieldValue alpha.air 1 + volScalarFieldValue alpha.water 0 + volScalarFieldValue alpha.solid 0 +); + +regions +( + boxToCell + { + box (0 0.3333 -0.1) (0.15 0.6501 0.1); + fieldValues + ( + volScalarFieldValue alpha.air 0.4 + volScalarFieldValue alpha.water 0 + volScalarFieldValue alpha.solid 0.6 + ); + } +); + + +// ************************************************************************* //