From 047e9640b65bbe7456cf8a2f686dcc75dca0361a Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 27 Feb 2020 10:25:02 +0000
Subject: [PATCH] reacting*EulerFoam/.../ThermalPhaseChange: Update to new
 definition of continuity error

The ThermalPhaseChangePhaseSystem stores the thermal phase change dmdt
used in the previous continuity error update and uses that to stabilize
the interfacial heat transfer calculations when phase fractions approach
zero.

Patch contributed by Juho Peltola, VTT.
---
 .../ThermalPhaseChangePhaseSystem.C           | 44 ++++++++++++++++---
 .../ThermalPhaseChangePhaseSystem.H           | 14 +++---
 2 files changed, 45 insertions(+), 13 deletions(-)

diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index 21e52385e9..77a86260d5 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -68,7 +68,8 @@ ThermalPhaseChangePhaseSystem
 )
 :
     BasePhaseSystem(mesh),
-    volatile_(this->template lookupOrDefault<word>("volatile", "none"))
+    volatile_(this->template lookupOrDefault<word>("volatile", "none")),
+    dmdt0s_(this->phases().size())
 {
     this->generatePairsAndSubModels
     (
@@ -325,9 +326,15 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
         const volScalarField dmdtf21(posPart(dmdtf));
         const volScalarField dmdtf12(negPart(dmdtf));
 
-        *eqns[phase1.name()] += - fvm::Sp(dmdtf21, he1) + dmdtf21*(K2 - K1);
+        *eqns[phase1.name()] +=
+            fvm::Sp(dmdt0s_[phase1.index()],he1)
+          - fvm::Sp(dmdtf21, he1)
+          + dmdtf21*K2 + dmdtf12*K1;
 
-        *eqns[phase2.name()] -= - fvm::Sp(dmdtf12, he2) + dmdtf12*(K1 - K2);
+        *eqns[phase2.name()] -=
+          - fvm::Sp(dmdt0s_[phase2.index()],he2)
+          - fvm::Sp(dmdtf12, he2)
+          + dmdtf12*K1 + dmdtf21*K2;
 
         if (this->saturationModels_.found(phasePairIter.key()))
         {
@@ -376,14 +383,12 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
         else
         {
             *eqns[phase1.name()] += dmdtf21*he2;
-
             *eqns[phase2.name()] -= dmdtf12*he1;
         }
 
         if (this->nDmdtLfs_.found(phasePairIter.key()))
         {
             *eqns[phase1.name()] += negPart(*this->nDmdtLfs_[pair]);
-
             *eqns[phase2.name()] -= posPart(*this->nDmdtLfs_[pair]);
         }
 
@@ -398,7 +403,6 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
             *eqns[phase2.name()] -=
                 phase2.thermo().p()*dmdtf/phase2.thermo().rho();
         }
-
     }
 
     return eqnsPtr;
@@ -457,6 +461,34 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::specieTransfer() const
 }
 
 
+template<class BasePhaseSystem>
+void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
+correctContinuityError()
+{
+    dmdt0s_ = PtrList<volScalarField>(this->phases().size());
+
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
+    {
+        const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
+        const volScalarField& dmdtf = *dmdtfIter();
+
+        addField(pair.phase1(), "dmdt", dmdtf, dmdt0s_);
+        addField(pair.phase2(), "dmdt", - dmdtf, dmdt0s_);
+    }
+
+    forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
+    {
+        const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
+        const volScalarField& nDmdtf = *nDmdtfIter();
+
+        addField(pair.phase1(), "dmdt", nDmdtf, dmdt0s_);
+        addField(pair.phase2(), "dmdt", - nDmdtf, dmdt0s_);
+    }
+
+    BasePhaseSystem::correctContinuityError();
+}
+
+
 template<class BasePhaseSystem>
 void
 Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
index e40a63d85e..756bffe50b 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
@@ -91,16 +91,9 @@ class ThermalPhaseChangePhaseSystem
         //- Name of the volatile specie
         word volatile_;
 
-        //- The saturation model used to evaluate Tsat = Tf
-        //autoPtr<saturationModel> saturationModel_;
-
         //- Saturation models used to evaluate Tsat = Tf
         saturationModelTable saturationModels_;
 
-        //- PhasePairs between which phaseChange occurs, e.g.,
-        //  "((gasI and liquid) (gasII and liquid))"
-        //List<phasePairKey> phaseChangePairKeys_;
-
         //- Mass transfer rates
         phaseSystem::dmdtfTable dmdtfs_;
 
@@ -110,6 +103,10 @@ class ThermalPhaseChangePhaseSystem
         //- Nucleate thermal energy transfer rates
         phaseSystem::dmdtfTable nDmdtLfs_;
 
+        //- Previous continuity error update phase dmdts for the heat transfer
+        //  function
+        PtrList<volScalarField> dmdt0s_;
+
 
     // Private member functions
 
@@ -153,6 +150,9 @@ public:
         virtual autoPtr<phaseSystem::specieTransferTable>
             specieTransfer() const;
 
+        //- Store phase dmdts at the during the continuity error update
+        virtual void correctContinuityError();
+
         //- Correct the interface thermodynamics
         virtual void correctInterfaceThermo();