From 12dc41d0fbc1c6cfd91b02fdf38ab77b2aadf1e6 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Mon, 11 Oct 2021 14:22:29 +0100
Subject: [PATCH] compressibleInterFoam::TEqn: Added internal energy
 formulation option

The temperature equation in compressibleInterFoam is derived from the phase
total internal energy equations including the kinetic energy source terms.
While this formulation handles pressure loss more accurately the kinetic energy
source terms can cause numerical problems and more likely to induce negative
temperatures in regions where the pressure drops rapidly.  For such cases it may
be beneficial to solve a temperature equation derived from the phase internal
energy equations, i.e. without the kinetic energy source terms which is now
selectable by setting the optional totalInternalEnergy entry in
constant/phaseProperties to false:

    totalInternalEnergy false;

When this entry is omitted the total energy form is used providing
backward-compatibility.
---
 .../solvers/multiphase/compressibleInterFoam/TEqn.H   |  4 ++++
 .../compressibleTwoPhaseMixture.C                     | 10 +++++++++-
 .../compressibleTwoPhaseMixture.H                     | 11 +++++++++++
 3 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
index 72eedb6559..7fe9a41e72 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
@@ -4,9 +4,13 @@
         fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T)
       - fvm::laplacian(turbulence.alphaEff(), T)
       + (
+             mixture.totalInternalEnergy()
+          ?
             fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
           + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()()
           - (U()&(fvModels.source(rho, U)&U)()) - contErr*K
+          :
+            p*fvc::div(fvc::absolute(phi, U))()()
         )
        *(
            alpha1()/mixture.thermo1().Cv()()
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.C
index 2494284680..c3b72b4fe8 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.C
@@ -44,6 +44,11 @@ Foam::compressibleTwoPhaseMixture::compressibleTwoPhaseMixture
     twoPhaseMixture(U.mesh()),
     interfaceProperties(alpha1(), alpha2(), U, *this),
 
+    totalInternalEnergy_
+    (
+        lookupOrDefault<Switch>("totalInternalEnergy", true)
+    ),
+
     p_
     (
         IOobject
@@ -204,8 +209,11 @@ Foam::tmp<Foam::scalarField> Foam::compressibleTwoPhaseMixture::nu
 
 bool Foam::compressibleTwoPhaseMixture::read()
 {
-    if (regIOobject::read())
+    if (twoPhaseMixture::read())
     {
+        totalInternalEnergy_ =
+            lookupOrDefault<Switch>("totalInternalEnergy", true);
+
         return interfaceProperties::read();
     }
     else
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.H
index d700882b9d..4b9284c1a7 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.H
@@ -59,6 +59,10 @@ class compressibleTwoPhaseMixture
 {
     // Private Data
 
+        //- Switch to choose between solving for internal energy
+        //  or total internal energy which is the default
+        Switch totalInternalEnergy_;
+
         //- Pressure
         volScalarField p_;
 
@@ -103,6 +107,13 @@ public:
 
     // Member Functions
 
+        //- Return true to solve for total internal energy
+        //  return false to solve for internal energy
+        bool totalInternalEnergy() const
+        {
+            return totalInternalEnergy_;
+        }
+
         //- Return pressure [Pa]
         volScalarField& p()
         {