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 + ); + } +); + + +// ************************************************************************* //