reactingEulerFoam: Corrected and rationalized pressure-work

In many publications and Euler-Euler codes the pressure-work term in the
total enthalpy is stated and implemented as -alpha*dp/dt rather than the
conservative form derived from the total internal energy equation
-d(alpha*p)/dt.  In order for the enthalpy and internal energy equations
to be consistent this error/simplification propagates to the total
internal energy equation as a spurious additional term p*d(alpha)/dt
which is included in the OpenFOAM Euler-Euler solvers and causes
stability and conservation issues.

I have now re-derived the energy equations for multiphase flow from
first-principles and implemented in the reactingEulerFoam solvers the
correct conservative form of pressure-work in both the internal energy
and enthalpy equations.

Additionally an optional limiter may be applied to the pressure-work
term in either of the energy forms to avoid spurious fluctuations in the
phase temperature in regions where the phase-fraction -> 0.  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.

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

@@ -77,11 +77,39 @@ void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctThermo()
 }
 
 
+template<class BasePhaseModel>
+Foam::tmp<Foam::volScalarField>
+Foam::AnisothermalPhaseModel<BasePhaseModel>::filterPressureWork
+(
+    const tmp<volScalarField>& pressureWork
+) const
+{
+    const volScalarField& alpha = *this;
+
+    scalar pressureWorkAlphaLimit =
+        this->thermo_->lookupOrDefault("pressureWorkAlphaLimit", 0.0);
+
+    if (pressureWorkAlphaLimit > 0)
+    {
+        return
+        (
+            max(alpha - pressureWorkAlphaLimit, scalar(0))
+           /max(alpha - pressureWorkAlphaLimit, pressureWorkAlphaLimit)
+        )*pressureWork;
+    }
+    else
+    {
+        return pressureWork;
+    }
+}
+
+
 template<class BasePhaseModel>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 {
     const volScalarField& alpha = *this;
+    const volVectorField& U = this->U();
     const surfaceScalarField& alphaPhi = this->alphaPhi();
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi();
 
@@ -93,7 +121,8 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 
     tmp<fvScalarMatrix> tEEqn
     (
-        fvm::ddt(alpha, this->rho(), he) + fvm::div(alphaRhoPhi, he)
+        fvm::ddt(alpha, this->rho(), he)
+      + fvm::div(alphaRhoPhi, he)
       - fvm::Sp(contErr, he)
 
       + fvc::ddt(alpha, this->rho(), K_) + fvc::div(alphaRhoPhi, K_)
@@ -112,13 +141,18 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
     // Add the appropriate pressure-work term
     if (he.name() == this->thermo_->phasePropertyName("e"))
     {
-        tEEqn.ref() +=
-            fvc::ddt(alpha)*this->thermo().p()
-          + fvc::div(alphaPhi, this->thermo().p());
+        tEEqn.ref() += filterPressureWork
+        (
+            fvc::div(fvc::absolute(alphaPhi, alpha, U), this->thermo().p())
+        );
     }
     else if (this->thermo_->dpdt())
     {
-        tEEqn.ref() -= alpha*this->fluid().dpdt();
+        tEEqn.ref() -= filterPressureWork
+        (
+           fvc::ddt(alpha, this->thermo().p())
+         + alpha*(this->fluid().dpdt() - fvc::ddt(this->thermo().p()))
+        );
     }
 
     return tEEqn;

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H

@@ -58,6 +58,15 @@ class AnisothermalPhaseModel
         volScalarField K_;
 
 
+    // Private member functions
+
+        //- Optionally filter the pressure work term as the phase-fraction -> 0
+        tmp<volScalarField> filterPressureWork
+        (
+            const tmp<volScalarField>& pressureWork
+        ) const;
+
+
 public:
 
     // Constructors