From a5eaa2bafd42b2f9f474ba59614a10315a69b674 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 5 Dec 2011 10:17:15 +0000
Subject: [PATCH] multiphaseEulerFoam: Change the specification of the residual
 drag and smooth the drag phase-fraction pre-factor

* The residual drag is now specified in terms of a residual phase-fraction and slip-velocity.
* The drag phase-fraction pre-factor is smoothed so that the computational molecule matches that
  of the transport terms.
---
 .../dragModels/dragModel/dragModel.C           | 14 ++++++++++----
 .../dragModels/dragModel/dragModel.H           | 12 +++++++++---
 .../multiphaseSystem/multiphaseSystem.C        | 18 +++++++++++++++---
 3 files changed, 34 insertions(+), 10 deletions(-)

diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C
index f3ec039684..acab521e20 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C
@@ -46,11 +46,17 @@ Foam::dragModel::dragModel
     interfaceDict_(interfaceDict),
     phase1_(phase1),
     phase2_(phase2),
-    residualDrag_
+    residualPhaseFraction_
     (
-        "residualDrag",
-        dimensionSet(1, -3, -1, 0, 0),
-        interfaceDict.lookup("residualDrag")
+        "residualPhaseFraction",
+        dimless,
+        interfaceDict.lookup("residualPhaseFraction")
+    ),
+    residualSlip_
+    (
+        "residualSlip",
+        dimVelocity,
+        interfaceDict.lookup("residualSlip")
     )
 {}
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H
index 37413d83e4..a775d91e45 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H
@@ -57,7 +57,8 @@ protected:
         const dictionary& interfaceDict_;
         const phaseModel& phase1_;
         const phaseModel& phase2_;
-        dimensionedScalar residualDrag_;
+        dimensionedScalar residualPhaseFraction_;
+        dimensionedScalar residualSlip_;
 
 public:
 
@@ -117,9 +118,14 @@ public:
             return phase2_;
         }
 
-        const dimensionedScalar& residualDrag() const
+        const dimensionedScalar& residualPhaseFraction() const
         {
-            return residualDrag_;
+            return residualPhaseFraction_;
+        }
+
+        const dimensionedScalar& residualSlip() const
+        {
+            return residualSlip_;
         }
 
         //- the dragfunction K used in the momentum eq.
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 2403597d9c..6ae60e42dc 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -30,6 +30,7 @@ License
 #include "MULES.H"
 #include "fvcSnGrad.H"
 #include "fvcFlux.H"
+#include "fvcAverage.H"
 
 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
 
@@ -626,9 +627,20 @@ Foam::multiphaseSystem::dragCoeffs() const
         (
             iter.key(),
             (
-                dm.phase1()*dm.phase2()
-               *dm.K(mag(dm.phase1().U() - dm.phase2().U()))
-              + dm.residualDrag()
+                max
+                (
+                    fvc::average(dm.phase1())*fvc::average(dm.phase2()),
+                    //dm.phase1()*dm.phase2(),
+                    dm.residualPhaseFraction()
+                )
+               *dm.K
+                (
+                    max
+                    (
+                        mag(dm.phase1().U() - dm.phase2().U()),
+                        dm.residualSlip()
+                    )
+                )
             ).ptr()
         );
     }