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.

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;

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
         )