From 03bfaaa4fa98317d9af848da5701a8acdcfb60f1 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Fri, 21 Jun 2024 16:40:13 +0100
Subject: [PATCH] fvModels: phaseChange: Corrections

The form of the sensible energy/enthalpy transfer and the sign of the
latent heat contribution have both been corrected
---
 .../phaseChange/multicomponentPhaseChange.C   | 30 +++++++++--------
 .../phaseChange/singleComponentPhaseChange.C  | 32 +++++++++++--------
 2 files changed, 34 insertions(+), 28 deletions(-)

diff --git a/src/fvModels/derived/phaseChange/multicomponentPhaseChange.C b/src/fvModels/derived/phaseChange/multicomponentPhaseChange.C
index 48d8305c41..581729c8c8 100644
--- a/src/fvModels/derived/phaseChange/multicomponentPhaseChange.C
+++ b/src/fvModels/derived/phaseChange/multicomponentPhaseChange.C
@@ -126,7 +126,19 @@ void Foam::fv::multicomponentPhaseChange::addSup
             const label speciei =
                 specieThermos()[i].species()[species()[mDoti]];
 
-            const volScalarField::Internal mDotIn(s*mDot(mDoti));
+            const volScalarField::Internal mDot(this->mDot(mDoti));
+
+            // Direct transfer of energy due to mass transfer
+            tmp<volScalarField::Internal> hs =
+                vfToVif(specieThermos()[i].hsi(speciei, p, Tchange));
+            if (energySemiImplicit_)
+            {
+                eqn += -fvm::SuSp(-s*mDot, heOrYi) + s*mDot*(hs - heOrYi);
+            }
+            else
+            {
+                eqn += s*mDot*hs;
+            }
 
             // Absolute enthalpies at the interface
             Pair<tmp<volScalarField::Internal>> has;
@@ -135,20 +147,10 @@ void Foam::fv::multicomponentPhaseChange::addSup
                 has[j] = vfToVif(specieThermos()[j].hai(speciei, p, Tchange));
             }
 
-            // Direct transfer of energy due to mass transfer
-            if (energySemiImplicit_)
-            {
-                eqn += -fvm::SuSp(-mDotIn, heOrYi) + mDotIn*(has[i]() - heOrYi);
-            }
-            else
-            {
-                eqn += mDotIn*has[i]();
-            }
-
             // Latent heat of phase change
-            eqn +=
-                (i == 0 ? Lfraction() - 1 : Lfraction())
-               *mDotIn*(has.second() - has.first());
+            eqn -=
+                (i == 0 ? 1 - Lfraction() : Lfraction())
+               *mDot*(has.second() - has.first());
         }
     }
     // Mass fraction equation
diff --git a/src/fvModels/derived/phaseChange/singleComponentPhaseChange.C b/src/fvModels/derived/phaseChange/singleComponentPhaseChange.C
index f999b9713b..0be4bd672c 100644
--- a/src/fvModels/derived/phaseChange/singleComponentPhaseChange.C
+++ b/src/fvModels/derived/phaseChange/singleComponentPhaseChange.C
@@ -120,7 +120,21 @@ void Foam::fv::singleComponentPhaseChange::addSup
         const volScalarField& p = this->p();
         const volScalarField Tchange(vifToVf(this->Tchange()));
 
-        const volScalarField::Internal mDotIn(s*mDot());
+        const volScalarField::Internal mDot(this->mDot());
+
+        // Direct transfer of energy due to mass transfer
+        tmp<volScalarField::Internal> hs =
+            specieThermos().valid()[i]
+          ? vfToVif(specieThermos()[i].hsi(specieis_[i], p, Tchange))
+          : vfToVif(thermos()[i].hs(p, Tchange));
+        if (energySemiImplicit_)
+        {
+            eqn += -fvm::SuSp(-s*mDot, heOrYi) + s*mDot*(hs - heOrYi);
+        }
+        else
+        {
+            eqn += s*mDot*hs;
+        }
 
         // Absolute enthalpies at the interface
         Pair<tmp<volScalarField::Internal>> has;
@@ -132,20 +146,10 @@ void Foam::fv::singleComponentPhaseChange::addSup
               : vfToVif(thermos()[j].ha(p, Tchange));
         }
 
-        // Direct transfer of energy due to mass transfer
-        if (energySemiImplicit_)
-        {
-            eqn += -fvm::SuSp(-mDotIn, heOrYi) + mDotIn*(has[i]() - heOrYi);
-        }
-        else
-        {
-            eqn += mDotIn*has[i]();
-        }
-
         // Latent heat of phase change
-        eqn +=
-            (i == 0 ? Lfraction() - 1 : Lfraction())
-           *mDotIn*(has.second() - has.first());
+        eqn -=
+            (i == 0 ? 1 - Lfraction() : Lfraction())
+           *mDot*(has.second() - has.first());
     }
     // Mass fraction equation
     else if