From dd3f68b1fac19093d64ae9f7beb8508493fd9188 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 4 Sep 2015 17:01:31 +0100
Subject: [PATCH] reactingEulerFoam: Rationalize the phaseSystem base-class

---
 .../MomentumTransferPhaseSystem.C             |  65 ++++++++++-
 .../MomentumTransferPhaseSystem.H             |   9 ++
 .../MovingPhaseModel/MovingPhaseModel.C       |   1 +
 .../phaseSystems/phaseSystem/phaseSystem.C    |   2 -
 .../phaseSystems/phaseSystem/phaseSystem.H    | 110 ++++++------------
 .../phaseSystems/phaseSystem/phaseSystemI.H   |   7 ++
 .../reactingTwoPhaseEulerFoam/pU/UEqns.H      |   3 -
 .../reactingTwoPhaseEulerFoam/pUf/UEqns.H     |   3 -
 .../twoPhaseSystem/twoPhaseSystem.C           |  26 ++---
 .../twoPhaseSystem/twoPhaseSystem.H           |  99 ++++++++++------
 10 files changed, 191 insertions(+), 134 deletions(-)

diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 9646a6a80c1..a209933b095 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/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 surfaceScalarField& phi(phase->phi());
+            const volVectorField& U = phase->U();
+            const surfaceScalarField& phi = phase->phi();
 
             *eqns[phase->name()] +=
               - Vm
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
index 5c9f057bedc..b6aeadd34ca 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
+++ b/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;
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index d64a5c862ef..090d29fafe0 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/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_)
     );
 }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index 3c5dafdfa6d..dee40846b63 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/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)
     {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index 84c69aa655c..bcb268bdc5f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/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;
 
-        //- Return the surface tension coefficient
-        tmp<volScalarField> sigma(const phasePairKey& key) const;
-
-        //- Return the drag coefficient
-        virtual tmp<volScalarField> Kd(const phasePairKey& key) const = 0;
-
-        //- Return the face drag coefficient
-        virtual tmp<surfaceScalarField> Kdf(const phasePairKey& key) const = 0;
+        //- Constant access the mixture flux
+        inline const surfaceScalarField& phi() const;
 
-        //- Return the virtual mass coefficient
-        virtual tmp<volScalarField> Vm(const phasePairKey& key) const = 0;
+        //- Access the mixture flux
+        inline surfaceScalarField& phi();
 
-        //- Return the face virtual mass coefficient
-        virtual tmp<surfaceScalarField> Vmf(const phasePairKey& key) const = 0;
+        //- Constant access the rate of change of the pressure
+        inline const volScalarField& dpdt() const;
 
-        //- Return the combined force (lift + wall-lubrication)
-        virtual tmp<volVectorField> F(const phasePairKey& key) const = 0;
+        //- Access the rate of change of the pressure
+        inline volScalarField& dpdt();
 
-        //- Return the combined face-force (lift + wall-lubrication)
-        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 surface tension coefficient
+        tmp<volScalarField> sigma(const phasePairKey& key) const;
 
-        //- Return the interfacial mass flow rate
-        virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
+        //- Return MRF zones
+        inline const IOMRFZoneList& MRF() const;
 
-        //- Return the momentum transfer matrices
-        virtual autoPtr<momentumTransferTable> momentumTransfer() const = 0;
+        //- Optional FV-options
+        inline fv::IOoptionList& fvOptions() const;
 
-        //- Return the heat transfer matrices
-        virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
+        //- Access a sub model between a phase pair
+        template <class modelType>
+        const modelType& lookupSubModel(const phasePair& key) const;
 
-        //- Return the mass transfer matrices
-        virtual autoPtr<massTransferTable> massTransfer() const = 0;
+        //- Access a sub model between two phases
+        template <class modelType>
+        const modelType& lookupSubModel
+        (
+            const phaseModel& dispersed,
+            const phaseModel& continuous
+        ) const;
 
         //- 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;
 };
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
index 35b61ab4516..343d31f49b7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
+++ b/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_;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H
index 31bd20ff7a6..92526304290 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H
+++ b/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);
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H
index eb1456bc158..83c3a9ce7f7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H
+++ b/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);
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 50b79b5672b..13a98b885f9 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/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));
-}
-
-
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
index 0ac971b790e..deca55fa6e5 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
+++ b/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
-        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 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;
 
-
-        // Access
-
-            //- 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;
+        //- Solve for the phase fractions
+        virtual void solve();
 };
 
 
-- 
GitLab