twoPhaseEulerFoam: Improve stability in the limit alpha -> 0

applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H

@@ -14,7 +14,6 @@
 
       + fvc::ddt(alpha1, rho1, K1) + fvc::div(alphaRhoPhi1, K1)
       - contErr1*K1
-
       + (
             he1.name() == thermo1.phasePropertyName("e")
           ? fvc::ddt(alpha1)*p + fvc::div(alphaPhi1, p)
@@ -27,9 +26,12 @@
            *fvc::interpolate(thermo1.alphaEff(phase1.turbulence().mut())),
             he1
         )
+    );
 
-     ==
+    he1Eqn.relax();
 
+    he1Eqn -=
+    (
         heatTransferCoeff*(thermo2.T() - thermo1.T())
       + heatTransferCoeff*he1/Cpv1
       - fvm::Sp(heatTransferCoeff/Cpv1, he1)
@@ -43,7 +45,6 @@
 
       + fvc::ddt(alpha2, rho2, K2) + fvc::div(alphaRhoPhi2, K2)
       - contErr2*K2
-
       + (
             he2.name() == thermo2.phasePropertyName("e")
           ? fvc::ddt(alpha2)*p + fvc::div(alphaPhi2, p)
@@ -56,23 +57,29 @@
            *fvc::interpolate(thermo2.alphaEff(phase2.turbulence().mut())),
             he2
         )
+    );
 
-     ==
+    he2Eqn.relax();
 
+    he2Eqn -=
+    (
         heatTransferCoeff*(thermo1.T() - thermo2.T())
       + heatTransferCoeff*he2/Cpv2
       - fvm::Sp(heatTransferCoeff/Cpv2, he2)
       + fvOptions(alpha2, rho2, he2)
     );
 
-    he1Eqn.relax();
     fvOptions.constrain(he1Eqn);
     he1Eqn.solve();
 
-    he2Eqn.relax();
     fvOptions.constrain(he2Eqn);
     he2Eqn.solve();
 
     thermo1.correct();
+    Info<< "min " << thermo1.T().name()
+        << " " << min(thermo1.T()).value() << endl;
+
     thermo2.correct();
+    Info<< "min " << thermo2.T().name()
+        << " " << min(thermo2.T()).value() << endl;
 }

applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H

@@ -16,13 +16,11 @@ volScalarField dragCoeff(fluid.dragCoeff());
     {
         U1Eqn =
         (
-            fvm::ddt(alpha1, rho1, U1)
-          + fvm::div(alphaRhoPhi1, U1)
+            fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
           - fvm::Sp(contErr1, U1)
           + mrfZones(alpha1*rho1 + virtualMassCoeff, U1)
           + phase1.turbulence().divDevRhoReff(U1)
          ==
-          - fvm::Sp(dragCoeff, U1)
           - liftForce
           - wallLubricationForce
           - turbulentDispersionForce
@@ -36,20 +34,19 @@ volScalarField dragCoeff(fluid.dragCoeff());
           + fvOptions(alpha1, rho1, U1)
         );
         U1Eqn.relax();
+        U1Eqn += fvm::Sp(dragCoeff, U1);
         fvOptions.constrain(U1Eqn);
     }
 
     {
         U2Eqn =
         (
-            fvm::ddt(alpha2, rho2, U2)
-          + fvm::div(alphaRhoPhi2, U2)
+            fvm::ddt(alpha2, rho2, U2) + fvm::div(alphaRhoPhi2, U2)
           - fvm::Sp(contErr2, U2)
           + mrfZones(alpha2*rho2 + virtualMassCoeff, U2)
           + phase2.turbulence().divDevRhoReff(U2)
          ==
-          - fvm::Sp(dragCoeff, U2)
-          + liftForce
+            liftForce
           + wallLubricationForce
           + turbulentDispersionForce
           - virtualMassCoeff
@@ -62,6 +59,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
           + fvOptions(alpha2, rho2, U2)
         );
         U2Eqn.relax();
+        U2Eqn += fvm::Sp(dragCoeff, U2);
         fvOptions.constrain(U2Eqn);
     }
 }