TwoResistanceHeatTransferPhaseSystem: Changed stabilisation term

The stabilisation term for the heat transfer with the interface
temperature has been changed to mirror the local heat transfer in the
phase, rather than an effective heat transfer across both phases. This
makes the stabilisation term match the actual temperature-based transfer
terms more accurately. The difference is particularly significant when
the mass transfer rate is high, and cases of this type gain a
significant stability benefit from this change as a result.

Patch contributed by Juho Peltola, VTT.

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -173,26 +173,17 @@ heatTransfer() const
             heatTransferModelIter().second()->K()
         );
 
-        const volScalarField KEff
-        (
-            Ks.first()()*Ks.second()()
-           /max
-            (
-                Ks.first()() + Ks.second()(),
-                dimensionedScalar(heatTransferModel::dimK, small)
-            )
-        );
-
         forAllConstIter(phasePair, pair, iter)
         {
             const phaseModel& phase = iter();
 
             const volScalarField& he(phase.thermo().he());
             const volScalarField Cpv(phase.thermo().Cpv());
+            const volScalarField& K(Ks[iter.index()]);
 
             *eqns[phase.name()] +=
-                Ks[iter.index()]*(Tf - phase.thermo().T())
-              + KEff/Cpv*he - fvm::Sp(KEff/Cpv, he);
+                K*(Tf - phase.thermo().T())
+              + K/Cpv*he - fvm::Sp(K/Cpv, he);
         }
     }