From f4eee20b9717cc60ae7f5496025dad5ec7a5bfba Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sat, 25 Mar 2023 16:51:07 +0000
Subject: [PATCH] multiphaseEuler: Adjust the drag and heat transfer
 coefficients for the continuous phase

so that the residualAlpha applied to stabilised the dispersed phase does not
affect the continuous phase in the limit of it becoming pure.
---
 .../MomentumTransferPhaseSystem.C             | 41 ++++++++++++++++---
 .../OneResistanceHeatTransferPhaseSystem.C    | 13 ++++--
 2 files changed, 46 insertions(+), 8 deletions(-)

diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 98f14ceaec..37906ed9fe 100644
--- a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -209,7 +209,14 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
             {
                 fvVectorMatrix& eqn = *eqns[iter().name()];
 
-                eqn -= fvm::Sp(K, eqn.psi());
+                const volScalarField phaseK
+                (
+                    iter.otherPhase()
+                   /max(iter.otherPhase(), iter.otherPhase().residualAlpha())
+                   *K
+                );
+
+                eqn -= fvm::Sp(phaseK, eqn.psi());
             }
         }
     }
@@ -481,7 +488,17 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::AFfs() const
 
         forAllConstIter(phaseInterface, interface, iter)
         {
-            addField(iter(), "AFf", Kf, AFfs);
+            addField
+            (
+                iter(),
+                "AFf",
+                fvc::interpolate
+                (
+                    iter.otherPhase()
+                   /max(iter.otherPhase(), iter.otherPhase().residualAlpha())
+                )*Kf,
+                AFfs
+            );
         }
     }
 
@@ -820,7 +837,12 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiKdPhis
             (
                 iter(),
                 "phiKdPhi",
-               -fvc::interpolate(rAUs[iter().index()]*K)
+                fvc::interpolate
+                (
+                   -iter.otherPhase()
+                   /max(iter.otherPhase(), iter.otherPhase().residualAlpha())
+                   *rAUs[iter().index()]*K
+                )
                *fvc::absolute
                 (
                     this->MRF().absolute(iter.otherPhase().phi()),
@@ -866,7 +888,13 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiKdPhifs
             (
                 iter(),
                 "phiKdPhif",
-               -rAUfs[iter().index()]*Kf
+               -rAUfs[iter().index()]
+               *fvc::interpolate
+                (
+                   -iter.otherPhase()
+                   /max(iter.otherPhase(), iter.otherPhase().residualAlpha())
+                )
+               *Kf
                *fvc::absolute
                 (
                     this->MRF().absolute(iter.otherPhase().phi()),
@@ -912,7 +940,10 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdUByAs
             (
                 iter(),
                 "KdUByA",
-               -rAUs[iter().index()]*K*iter.otherPhase().U(),
+               -rAUs[iter().index()]
+               *iter.otherPhase()
+               /max(iter.otherPhase(), iter.otherPhase().residualAlpha())
+               *K*iter.otherPhase().U(),
                 KdUByAs
             );
         }
diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C
index 31dab0bd25..1380a37a52 100644
--- a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C
+++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2020-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2020-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,9 +95,16 @@ heatTransfer() const
             const volScalarField& he = phase.thermo().he();
             const volScalarField Cpv(phase.thermo().Cpv());
 
+            const volScalarField phaseK
+            (
+                iter.otherPhase()
+               /max(iter.otherPhase(), iter.otherPhase().residualAlpha())
+               *K
+            );
+
             *eqns[phase.name()] +=
-                K*(otherPhase.thermo().T() - phase.thermo().T() + he/Cpv)
-              - fvm::Sp(K/Cpv, he);
+                phaseK*(otherPhase.thermo().T() - phase.thermo().T() + he/Cpv)
+              - fvm::Sp(phaseK/Cpv, he);
         }
     }