From ad476af9b33b995ffbb063adb0421d13da0c2db4 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Wed, 9 Nov 2016 11:14:26 +0000
Subject: [PATCH] reactingEulerFoam, twoPhaseEulerFoam: Reinstated interfacial
 pressure-work

Added the interfacial pressure-work terms according to:

Ishii, M., Hibiki, T.,
Thermo-fluid dynamics of two-phase flow,
ISBN-10: 0-387-28321-8, 2006

While this is the most common approach to handling the interfacial
pressure-work it introduces numerical stability issues in regions of low
phase-fraction and rapid flow deformation.  To alleviate this problem an
optional limiter may be applied to the pressure-work term in either of
the energy forms.  This may specified in the
"thermophysicalProperties.<phase>" file, e.g.

pressureWorkAlphaLimit 1e-3;

which sets the pressure work term to 0 for phase-fractions below 1e-3.

For particularly unstable cases a limit of 1e-2 may be necessary.
---
 .../AnisothermalPhaseModel/AnisothermalPhaseModel.C        | 7 ++-----
 applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H  | 6 ++++--
 2 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index 4285a40c0c..b5dd4ab37a 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -144,15 +144,12 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
         tEEqn.ref() += filterPressureWork
         (
             fvc::div(fvc::absolute(alphaPhi, alpha, U), this->thermo().p())
+          + this->thermo().p()*fvc::ddt(alpha)
         );
     }
     else if (this->thermo_->dpdt())
     {
-        tEEqn.ref() -= filterPressureWork
-        (
-           fvc::ddt(alpha, this->thermo().p())
-         + alpha*(this->fluid().dpdt() - fvc::ddt(this->thermo().p()))
-        );
+        tEEqn.ref() -= filterPressureWork(alpha*this->fluid().dpdt());
     }
 
     return tEEqn;
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
index 25b97c17a1..308e7589f9 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
@@ -16,7 +16,8 @@
       - contErr1*K1
       + (
             he1.name() == thermo1.phasePropertyName("e")
-          ? fvc::ddt(alpha1)*p + fvc::div(alphaPhi1, p)
+          ? fvc::div(fvc::absolute(alphaPhi1, alpha1, U1), p)
+          + p*fvc::ddt(alpha1)
           : -alpha1*dpdt
         )
 
@@ -48,7 +49,8 @@
       - contErr2*K2
       + (
             he2.name() == thermo2.phasePropertyName("e")
-          ? fvc::ddt(alpha2)*p + fvc::div(alphaPhi2, p)
+          ? fvc::div(fvc::absolute(alphaPhi2, alpha2, U2), p)
+          + p*fvc::ddt(alpha1)
           : -alpha2*dpdt
         )