src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -177,7 +178,11 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
 
         calcLimiter(phi, limiterField);
 
-        return limiterField;
+        return tmp<surfaceScalarField>::New
+        (
+            limiterFieldName,
+            limiterField
+        );
     }
     else
     {

src/phaseSystemModels/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2019 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -704,7 +705,7 @@ void Foam::multiphaseSystem::solve()
         phase.alphaRhoPhiRef() =
             fvc::interpolate(phase.rho())*phase.alphaPhi();
 
-        phase.clip(0, 1);
+        phase.clip(SMALL, 1 - SMALL);
     }
 
     calcAlphas();

src/phaseSystemModels/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -339,7 +340,7 @@ void Foam::twoPhaseSystem::solve()
             << endl;
 
         // Ensure the phase-fractions are bounded
-        alpha1.clip(0, 1);
+        alpha1.clip(SMALL, 1 - SMALL);
 
         // Update the phase-fraction of the other phase
         alpha2 = scalar(1) - alpha1;