reactingEulerFoam: Rationalize the phaseSystem base-class

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C

@@ -155,6 +155,63 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kdf
 }
 
 
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kd
+(
+    const Foam::phaseModel& phase
+) const
+{
+    tmp<volScalarField> tKd
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("Kd", phase.name()),
+                this->mesh_.time().timeName(),
+                this->mesh_
+            ),
+            this->mesh_,
+            dimensionedScalar
+            (
+                IOobject::groupName("Kd", phase.name()),
+                dimensionSet(1, -3, -1, 0, 0),
+                0
+            )
+        )
+    );
+
+    // Add the implicit part of the drag force
+    forAllConstIter
+    (
+        phaseSystem::KdTable,
+        Kds_,
+        KdIter
+    )
+    {
+        const volScalarField& K(*KdIter());
+
+        const phasePair& pair(this->phasePairs_[KdIter.key()]);
+
+        const phaseModel* phase1 = &pair.phase1();
+        const phaseModel* phase2 = &pair.phase2();
+
+        forAllConstIter(phasePair, pair, iter)
+        {
+            if (phase1 == &phase)
+            {
+                tKd() += K;
+            }
+
+            Swap(phase1, phase2);
+        }
+    }
+
+    return tKd;
+}
+
+
 template<class BasePhaseSystem>
 Foam::tmp<Foam::volScalarField>
 Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vm
@@ -385,6 +442,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
     }
 
     // Add the implicit part of the drag force
+    /* ***HGW Currently this is handled in the pEqn
     forAllConstIter
     (
         phaseSystem::KdTable,
@@ -401,13 +459,14 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
 
         forAllConstIter(phasePair, pair, iter)
         {
-            const volVectorField& U(phase->U());
+            const volVectorField& U = phase->U();
 
             *eqns[phase->name()] -= fvm::Sp(K, U);
 
             Swap(phase, otherPhase);
         }
     }
+    */
 
     // Update the virtual mass coefficients
     forAllConstIter
@@ -437,8 +496,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
 
         forAllConstIter(phasePair, pair, iter)
         {
-            const volVectorField& U(phase->U());
+            const volVectorField& U = phase->U();
-            const surfaceScalarField& phi(phase->phi());
+            const surfaceScalarField& phi = phase->phi();
 
             *eqns[phase->name()] +=
               - Vm

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H

@@ -148,12 +148,21 @@ public:
 
     // Member Functions
 
+        //- Constant access to drag coefficients
+        virtual const phaseSystem::KdTable& Kds() const
+        {
+            return Kds_;
+        }
+
         //- Return the drag coefficient
         virtual tmp<volScalarField> Kd(const phasePairKey& key) const;
 
         //- Return the face drag coefficient
         virtual tmp<surfaceScalarField> Kdf(const phasePairKey& key) const;
 
+        //- Return the drag coefficient for phase
+        virtual tmp<volScalarField> Kd(const phaseModel& phase) const;
+
         //- Return the virtual mass coefficient
         virtual tmp<volScalarField> Vm(const phasePairKey& key) const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C

@@ -260,6 +260,7 @@ Foam::MovingPhaseModel<BasePhaseModel>::UEqn()
         fvm::ddt(*this, this->thermo().rho(), U_)
       + fvm::div(alphaRhoPhi_, U_)
       - fvm::Sp(continuityError_, U_)
+      + this->fluid().MRF().DDt(*this*this->thermo().rho(), U_)
       + turbulence_->divDevRhoReff(U_)
     );
 }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C

@@ -282,8 +282,6 @@ void Foam::phaseSystem::correctKinematics()
         updateDpdt = updateDpdt || phaseModelIter().thermo().dpdt();
     }
 
-    //phaseModelTable::iterator iter = phaseModels_.begin();
-
     // Update the pressure time-derivative if required
     if (updateDpdt)
     {

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H

@@ -118,10 +118,6 @@ public:
             >
             massTransferTable;
 
-        // typedef
-        //     HashTable<autoPtr<phaseModel>, word, word::hash>
-        //     phaseModelTable;
-
         typedef PtrDictionary<phaseModel> phaseModelTable;
 
 
@@ -284,48 +280,56 @@ public:
 
     // Member Functions
 
+        //- Constant access the mesh
+        inline const fvMesh& mesh() const;
+
+        //- Constant access the phase models
+        inline const phaseModelTable& phases() const;
+
+        //- Access the phase models
+        inline phaseModelTable& phases();
+
+        //- Constant access the phase pairs
+        inline const phasePairTable& phasePairs() const;
+
         //- Return the mixture density
         tmp<volScalarField> rho() const;
 
         //- Return the mixture velocity
         tmp<volVectorField> U() const;
 
+        //- Constant access the mixture flux
+        inline const surfaceScalarField& phi() const;
+
+        //- Access the mixture flux
+        inline surfaceScalarField& phi();
+
+        //- Constant access the rate of change of the pressure
+        inline const volScalarField& dpdt() const;
+
+        //- Access the rate of change of the pressure
+        inline volScalarField& dpdt();
+
         //- Return the surface tension coefficient
         tmp<volScalarField> sigma(const phasePairKey& key) const;
 
-        //- Return the drag coefficient
+        //- Return MRF zones
-        virtual tmp<volScalarField> Kd(const phasePairKey& key) const = 0;
+        inline const IOMRFZoneList& MRF() const;
 
-        //- Return the face drag coefficient
+        //- Optional FV-options
-        virtual tmp<surfaceScalarField> Kdf(const phasePairKey& key) const = 0;
+        inline fv::IOoptionList& fvOptions() const;
 
-        //- Return the virtual mass coefficient
+        //- Access a sub model between a phase pair
-        virtual tmp<volScalarField> Vm(const phasePairKey& key) const = 0;
+        template <class modelType>
+        const modelType& lookupSubModel(const phasePair& key) const;
 
-        //- Return the face virtual mass coefficient
+        //- Access a sub model between two phases
-        virtual tmp<surfaceScalarField> Vmf(const phasePairKey& key) const = 0;
+        template <class modelType>
-
+        const modelType& lookupSubModel
-        //- Return the combined force (lift + wall-lubrication)
+        (
-        virtual tmp<volVectorField> F(const phasePairKey& key) const = 0;
+            const phaseModel& dispersed,
-
+            const phaseModel& continuous
-        //- Return the combined face-force (lift + wall-lubrication)
+        ) const;
-        virtual tmp<surfaceScalarField> Ff(const phasePairKey& key) const = 0;
-
-        //- Return the turbulent diffusivity
-        //  Multiplies the phase-fraction gradient
-        virtual tmp<volScalarField> D(const phasePairKey& key) const = 0;
-
-        //- Return the interfacial mass flow rate
-        virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
-
-        //- Return the momentum transfer matrices
-        virtual autoPtr<momentumTransferTable> momentumTransfer() const = 0;
-
-        //- Return the heat transfer matrices
-        virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
-
-        //- Return the mass transfer matrices
-        virtual autoPtr<massTransferTable> massTransfer() const = 0;
 
         //- Solve for the phase fractions
         virtual void solve();
@@ -347,48 +351,6 @@ public:
 
         //- Read base phaseProperties dictionary
         virtual bool read();
-
-
-        // Access
-
-            //- Constant access the mesh
-            inline const fvMesh& mesh() const;
-
-            //- Constant access the phase models
-            inline const phaseModelTable& phases() const;
-
-            //- Access the phase models
-            inline phaseModelTable& phases();
-
-            //- Constant access the mixture flux
-            inline const surfaceScalarField& phi() const;
-
-            //- Access the mixture flux
-            inline surfaceScalarField& phi();
-
-            //- Constant access the rate of change of the pressure
-            inline const volScalarField& dpdt() const;
-
-            //- Access the rate of change of the pressure
-            inline volScalarField& dpdt();
-
-            //- Return MRF zones
-            inline const IOMRFZoneList& MRF() const;
-
-            //- Optional FV-options
-            inline fv::IOoptionList& fvOptions() const;
-
-            //- Access a sub model between a phase pair
-            template <class modelType>
-            const modelType& lookupSubModel(const phasePair& key) const;
-
-            //- Access a sub model between two phases
-            template <class modelType>
-            const modelType& lookupSubModel
-            (
-                const phaseModel& dispersed,
-                const phaseModel& continuous
-            ) const;
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H

@@ -45,6 +45,13 @@ Foam::phaseSystem::phases()
 }
 
 
+inline const Foam::phaseSystem::phasePairTable&
+Foam::phaseSystem::phasePairs() const
+{
+    return phasePairs_;
+}
+
+
 inline const Foam::surfaceScalarField& Foam::phaseSystem::phi() const
 {
     return phi_;

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H

@@ -1,8 +1,5 @@
 Info<< "Constructing momentum equations" << endl;
 
-MRF.correctBoundaryVelocity(U1);
-MRF.correctBoundaryVelocity(U2);
-
 fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime);
 fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H

@@ -1,8 +1,5 @@
 Info<< "Constructing face momentum equations" << endl;
 
-MRF.correctBoundaryVelocity(U1);
-MRF.correctBoundaryVelocity(U2);
-
 fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime);
 fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -85,6 +85,12 @@ Foam::twoPhaseSystem::sigma() const
 }
 
 
+const Foam::dragModel& Foam::twoPhaseSystem::drag(const phaseModel& phase) const
+{
+    return lookupSubModel<dragModel>(phase, otherPhase(phase));
+}
+
+
 Foam::tmp<Foam::volScalarField>
 Foam::twoPhaseSystem::Kd() const
 {
@@ -105,6 +111,13 @@ Foam::twoPhaseSystem::Kdf() const
 }
 
 
+const Foam::virtualMassModel&
+Foam::twoPhaseSystem::virtualMass(const phaseModel& phase) const
+{
+    return lookupSubModel<virtualMassModel>(phase, otherPhase(phase));
+}
+
+
 Foam::tmp<Foam::volScalarField>
 Foam::twoPhaseSystem::Vm() const
 {
@@ -409,17 +422,4 @@ void Foam::twoPhaseSystem::solve()
 }
 
 
-const Foam::dragModel& Foam::twoPhaseSystem::drag(const phaseModel& phase) const
-{
-    return lookupSubModel<dragModel>(phase, otherPhase(phase));
-}
-
-
-const Foam::virtualMassModel&
-Foam::twoPhaseSystem::virtualMass(const phaseModel& phase) const
-{
-    return lookupSubModel<virtualMassModel>(phase, otherPhase(phase));
-}
-
-
 // ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H

@@ -53,6 +53,35 @@ class twoPhaseSystem
 :
     public phaseSystem
 {
+    // Private member functions
+
+        //- Return the drag coefficient for phase pair
+        virtual tmp<volScalarField> Kd(const phasePairKey& key) const = 0;
+
+        //- Return the face drag coefficient for phase pair
+        virtual tmp<surfaceScalarField> Kdf(const phasePairKey& key) const = 0;
+
+        //- Return the virtual mass coefficient for phase pair
+        virtual tmp<volScalarField> Vm(const phasePairKey& key) const = 0;
+
+        //- Return the face virtual mass coefficient for phase pair
+        virtual tmp<surfaceScalarField> Vmf(const phasePairKey& key) const = 0;
+
+        //- Return the combined force (lift + wall-lubrication) for phase pair
+        virtual tmp<volVectorField> F(const phasePairKey& key) const = 0;
+
+        //- Return the combined face-force (lift + wall-lubrication)
+        //  for phase pair
+        virtual tmp<surfaceScalarField> Ff(const phasePairKey& key) const = 0;
+
+        //- Return the turbulent diffusivity for phase pair
+        //  Multiplies the phase-fraction gradient
+        virtual tmp<volScalarField> D(const phasePairKey& key) const = 0;
+
+        //- Return the interfacial mass flow rate for phase pair
+        virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
+
+
 protected:
 
     // Protected data
@@ -103,28 +132,50 @@ public:
 
     // Member Functions
 
-        //- Solve for the phase fractions
+        //- Constant access phase model 1
-        virtual void solve();
+        const phaseModel& phase1() const;
+
+        //- Access phase model 1
+        phaseModel& phase1();
+
+        //- Constant access phase model 2
+        const phaseModel& phase2() const;
+
+        //- Access phase model 2
+        phaseModel& phase2();
+
+        //- Constant access the phase not given as an argument
+        const phaseModel& otherPhase
+        (
+            const phaseModel& phase
+        ) const;
+
+        //- Return the momentum transfer matrices
+        virtual autoPtr<momentumTransferTable> momentumTransfer() const = 0;
+
+        //- Return the heat transfer matrices
+        virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
+
+        //- Return the mass transfer matrices
+        virtual autoPtr<massTransferTable> massTransfer() const = 0;
 
         using phaseSystem::sigma;
-        using phaseSystem::Kd;
-        using phaseSystem::Kdf;
-        using phaseSystem::Vm;
-        using phaseSystem::Vmf;
-        using phaseSystem::F;
-        using phaseSystem::Ff;
-        using phaseSystem::D;
-        using phaseSystem::dmdt;
 
         //- Return the surface tension coefficient
         tmp<volScalarField> sigma() const;
 
+        //- Return the drag model for the given phase
+        const dragModel& drag(const phaseModel& phase) const;
+
         //- Return the drag coefficient
         tmp<volScalarField> Kd() const;
 
         //- Return the face drag coefficient
         tmp<surfaceScalarField> Kdf() const;
 
+        //- Return the virtual mass model for the given phase
+        const virtualMassModel& virtualMass(const phaseModel& phase) const;
+
         //- Return the virtual mass coefficient
         tmp<volScalarField> Vm() const;
 
@@ -144,32 +195,8 @@ public:
         //- Return the interfacial mass flow rate
         tmp<volScalarField> dmdt() const;
 
-
+        //- Solve for the phase fractions
-        // Access
+        virtual void solve();
-
-            //- Constant access phase model 1
-            const phaseModel& phase1() const;
-
-            //- Access phase model 1
-            phaseModel& phase1();
-
-            //- Constant access phase model 2
-            const phaseModel& phase2() const;
-
-            //- Access phase model 2
-            phaseModel& phase2();
-
-            //- Constant access the phase not given as an argument
-            const phaseModel& otherPhase
-            (
-                const phaseModel& phase
-            ) const;
-
-            //- Return the drag model for the given phase
-            const dragModel& drag(const phaseModel& phase) const;
-
-            //- Return the virtual mass model for the given phase
-            const virtualMassModel& virtualMass(const phaseModel& phase) const;
 };