Multi-phase solvers: Improved handling of inflow/outflow BCs in MULES

Avoids slight phase-fraction unboundedness at entertainment BCs and improved
robustness.

Additionally the phase-fractions in the multi-phase (rather than two-phase)
solvers are adjusted to avoid the slow growth of inconsistency ("drift") caused
by solving for all of the phase-fractions rather than deriving one from the
others.

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,11 +64,8 @@ void Foam::multiphaseSystem::solveAlphas()
 
     forAllIter(PtrDictionary<phaseModel>, phases_, iter)
     {
-        phaseModel& phase1 = iter();
+        phaseModel& phase = iter();
-        volScalarField& alpha1 = phase1;
+        volScalarField& alpha1 = phase;
-
-        phase1.alphaPhi() =
-            dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
 
         alphaPhiCorrs.set
         (
@@ -79,7 +76,7 @@ void Foam::multiphaseSystem::solveAlphas()
                 fvc::flux
                 (
                     phi_,
-                    phase1,
+                    phase,
                     "div(phi," + alpha1.name() + ')'
                 )
             )
@@ -92,13 +89,13 @@ void Foam::multiphaseSystem::solveAlphas()
             phaseModel& phase2 = iter2();
             volScalarField& alpha2 = phase2;
 
-            if (&phase2 == &phase1) continue;
+            if (&phase2 == &phase) continue;
 
-            surfaceScalarField phir(phase1.phi() - phase2.phi());
+            surfaceScalarField phir(phase.phi() - phase2.phi());
 
             scalarCoeffSymmTable::const_iterator cAlpha
             (
-                cAlphas_.find(interfacePair(phase1, phase2))
+                cAlphas_.find(interfacePair(phase, phase2))
             );
 
             if (cAlpha != cAlphas_.end())
@@ -108,7 +105,7 @@ void Foam::multiphaseSystem::solveAlphas()
                     (mag(phi_) + mag(phir))/mesh_.magSf()
                 );
 
-                phir += min(cAlpha()*phic, max(phic))*nHatf(phase1, phase2);
+                phir += min(cAlpha()*phic, max(phic))*nHatf(phase, phase2);
             }
 
             word phirScheme
@@ -119,39 +116,18 @@ void Foam::multiphaseSystem::solveAlphas()
             alphaPhiCorr += fvc::flux
             (
                 -fvc::flux(-phir, phase2, phirScheme),
-                phase1,
+                phase,
                 phirScheme
             );
         }
 
-        surfaceScalarField::Boundary& alphaPhiCorrBf =
+        phase.correctInflowOutflow(alphaPhiCorr);
-            alphaPhiCorr.boundaryFieldRef();
-
-        // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhiCorrBf, patchi)
-        {
-            fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];
-
-            if (!alphaPhiCorrp.coupled())
-            {
-                const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
-                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
-                forAll(alphaPhiCorrp, facei)
-                {
-                    if (phi1p[facei] < 0)
-                    {
-                        alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
-                    }
-                }
-            }
-        }
 
         MULES::limit
         (
             1.0/mesh_.time().deltaT().value(),
             geometricOneField(),
-            phase1,
+            phase,
             phi_,
             alphaPhiCorr,
             zeroField(),
@@ -182,29 +158,30 @@ void Foam::multiphaseSystem::solveAlphas()
 
     forAllIter(PtrDictionary<phaseModel>, phases_, iter)
     {
-        phaseModel& phase1 = iter();
+        phaseModel& phase = iter();
 
         surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
-        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase1);
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
+        phase.correctInflowOutflow(alphaPhi);
 
         MULES::explicitSolve
         (
             geometricOneField(),
-            phase1,
+            phase,
             alphaPhi,
             zeroField(),
             zeroField()
         );
 
-        phase1.alphaPhi() += alphaPhi;
+        phase.alphaPhi() = alphaPhi;
 
-        Info<< phase1.name() << " volume fraction, min, max = "
+        Info<< phase.name() << " volume fraction, min, max = "
-            << phase1.weightedAverage(mesh_.V()).value()
+            << phase.weightedAverage(mesh_.V()).value()
-            << ' ' << min(phase1).value()
+            << ' ' << min(phase).value()
-            << ' ' << max(phase1).value()
+            << ' ' << max(phase).value()
             << endl;
 
-        sumAlpha += phase1;
+        sumAlpha += phase;
 
         phasei++;
     }
@@ -215,6 +192,15 @@ void Foam::multiphaseSystem::solveAlphas()
         << ' ' << max(sumAlpha).value()
         << endl;
 
+    // Correct the sum of the phase-fractions to avoid 'drift'
+    volScalarField sumCorr(1.0 - sumAlpha);
+    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+    {
+        phaseModel& phase = iter();
+        volScalarField& alpha = phase;
+        alpha += alpha*sumCorr;
+    }
+
     calcAlphas();
 }
 

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -236,6 +236,24 @@ bool Foam::phaseModel::read(const dictionary& phaseDict)
 }
 
 
+void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
+{
+    surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+    const volScalarField::Boundary& alphaBf = boundaryField();
+    const surfaceScalarField::Boundary& phiBf = phi().boundaryField();
+
+    forAll(alphaPhiBf, patchi)
+    {
+        fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
+
+        if (!alphaPhip.coupled())
+        {
+            alphaPhip = phiBf[patchi]*alphaBf[patchi];
+        }
+    }
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::phaseModel::d() const
 {
     return dPtr_().d();

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -208,6 +208,9 @@ public:
             return alphaPhi_;
         }
 
+        //- Ensure that the flux at inflow/outflow BCs is preserved
+        void correctInflowOutflow(surfaceScalarField& alphaPhi) const;
+
         //- Correct the phase properties
         void correct();
 

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -681,6 +681,14 @@ void Foam::multiphaseMixture::solveAlphas
         << ' ' << max(sumAlpha).value()
         << endl;
 
+    // Correct the sum of the phase-fractions to avoid 'drift'
+    volScalarField sumCorr(1.0 - sumAlpha);
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        phase& alpha = iter();
+        alpha += alpha*sumCorr;
+    }
+
     calcAlphas();
 }
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -165,6 +165,24 @@ bool Foam::phaseModel::compressible() const
 }
 
 
+void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
+{
+    surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+    const volScalarField::Boundary& alphaBf = boundaryField();
+    const surfaceScalarField::Boundary& phiBf = phi()().boundaryField();
+
+    forAll(alphaPhiBf, patchi)
+    {
+        fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
+
+        if (!alphaPhip.coupled())
+        {
+            alphaPhip = phiBf[patchi]*alphaBf[patchi];
+        }
+    }
+}
+
+
 const Foam::tmp<Foam::volScalarField>& Foam::phaseModel::divU() const
 {
     NotImplemented;

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -283,6 +283,9 @@ public:
             //- Access the mass flux of the phase
             virtual surfaceScalarField& alphaRhoPhi() = 0;
 
+            //- Ensure that the flux at inflow/outflow BCs is preserved
+            void correctInflowOutflow(surfaceScalarField& alphaPhi) const;
+
 
         // Transport
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,15 +71,17 @@ void Foam::multiphaseSystem::solveAlphas()
 {
     bool LTS = fv::localEulerDdt::enabled(mesh_);
 
+    forAll(phases(), phasei)
+    {
+        phases()[phasei].correctBoundaryConditions();
+    }
+
     PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
     forAll(phases(), phasei)
     {
         phaseModel& phase = phases()[phasei];
         volScalarField& alpha1 = phase;
 
-        phase.alphaPhi() =
-            dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
-
         alphaPhiCorrs.set
         (
             phasei,
@@ -134,28 +136,7 @@ void Foam::multiphaseSystem::solveAlphas()
             );
         }
 
-        surfaceScalarField::Boundary& alphaPhiCorrBf =
+        phase.correctInflowOutflow(alphaPhiCorr);
-            alphaPhiCorr.boundaryFieldRef();
-
-        // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhiCorr.boundaryField(), patchi)
-        {
-            fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];
-
-            if (!alphaPhiCorrp.coupled())
-            {
-                const scalarField& phi1p = phase.phi().boundaryField()[patchi];
-                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
-                forAll(alphaPhiCorrp, facei)
-                {
-                    if (phi1p[facei] < 0)
-                    {
-                        alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
-                    }
-                }
-            }
-        }
 
         if (LTS)
         {
@@ -215,8 +196,9 @@ void Foam::multiphaseSystem::solveAlphas()
         phaseModel& phase = phases()[phasei];
         volScalarField& alpha = phase;
 
-        surfaceScalarField& alphaPhic = alphaPhiCorrs[phasei];
+        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
-        alphaPhic += upwind<scalar>(mesh_, phi_).flux(phase);
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
+        phase.correctInflowOutflow(alphaPhi);
 
         volScalarField::Internal Sp
         (
@@ -298,12 +280,12 @@ void Foam::multiphaseSystem::solveAlphas()
         (
             geometricOneField(),
             alpha,
-            alphaPhic,
+            alphaPhi,
             Sp,
             Su
         );
 
-        phase.alphaPhi() += alphaPhic;
+        phase.alphaPhi() = alphaPhi;
 
         Info<< phase.name() << " volume fraction, min, max = "
             << phase.weightedAverage(mesh_.V()).value()
@@ -319,6 +301,14 @@ void Foam::multiphaseSystem::solveAlphas()
         << ' ' << min(sumAlpha).value()
         << ' ' << max(sumAlpha).value()
         << endl;
+
+    // Correct the sum of the phase-fractions to avoid 'drift'
+    volScalarField sumCorr(1.0 - sumAlpha);
+    forAll(phases(), phasei)
+    {
+        volScalarField& alpha = phases()[phasei];
+        alpha += alpha*sumCorr;
+    }
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -201,7 +201,6 @@ void Foam::twoPhaseSystem::solve()
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
 
-    const surfaceScalarField& phi = this->phi();
     const surfaceScalarField& phi1 = phase1_.phi();
     const surfaceScalarField& phi2 = phase2_.phi();
 
@@ -236,9 +235,7 @@ void Foam::twoPhaseSystem::solve()
     }
 
     alpha1.correctBoundaryConditions();
-    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
 
-    surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
 
     tmp<surfaceScalarField> alphaDbyA;
@@ -279,7 +276,7 @@ void Foam::twoPhaseSystem::solve()
             ),
             // Divergence term is handled explicitly to be
             // consistent with the explicit transport solution
-            fvc::div(phi)*min(alpha1, scalar(1))
+            fvc::div(phi_)*min(alpha1, scalar(1))
         );
 
         if (tdgdt.valid())
@@ -300,11 +297,11 @@ void Foam::twoPhaseSystem::solve()
             }
         }
 
-        surfaceScalarField alphaPhic1
+        surfaceScalarField alphaPhi1
         (
             fvc::flux
             (
-                phic,
+                phi_,
                 alpha1,
                 alphaScheme
             )
@@ -316,28 +313,7 @@ void Foam::twoPhaseSystem::solve()
             )
         );
 
-        surfaceScalarField::Boundary& alphaPhic1Bf =
+        phase1_.correctInflowOutflow(alphaPhi1);
-            alphaPhic1.boundaryFieldRef();
-
-        // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhic1Bf, patchi)
-        {
-            fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi];
-
-            if (!alphaPhic1p.coupled())
-            {
-                const scalarField& phi1p = phi1.boundaryField()[patchi];
-                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
-                forAll(alphaPhic1p, facei)
-                {
-                    if (phi1p[facei] < 0)
-                    {
-                        alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
-                    }
-                }
-            }
-        }
 
         if (nAlphaSubCycles > 1)
         {
@@ -355,14 +331,14 @@ void Foam::twoPhaseSystem::solve()
                 !(++alphaSubCycle).end();
             )
             {
-                surfaceScalarField alphaPhic10(alphaPhic1);
+                surfaceScalarField alphaPhi10(alphaPhi1);
 
                 MULES::explicitSolve
                 (
                     geometricOneField(),
                     alpha1,
-                    phi,
+                    phi_,
-                    alphaPhic10,
+                    alphaPhi10,
                     (alphaSubCycle.index()*Sp)(),
                     (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
                     phase1_.alphaMax(),
@@ -371,11 +347,11 @@ void Foam::twoPhaseSystem::solve()
 
                 if (alphaSubCycle.index() == 1)
                 {
-                    phase1_.alphaPhi() = alphaPhic10;
+                    phase1_.alphaPhi() = alphaPhi10;
                 }
                 else
                 {
-                    phase1_.alphaPhi() += alphaPhic10;
+                    phase1_.alphaPhi() += alphaPhi10;
                 }
             }
 
@@ -387,15 +363,15 @@ void Foam::twoPhaseSystem::solve()
             (
                 geometricOneField(),
                 alpha1,
-                phi,
+                phi_,
-                alphaPhic1,
+                alphaPhi1,
                 Sp,
                 Su,
                 phase1_.alphaMax(),
                 0
             );
 
-            phase1_.alphaPhi() = alphaPhic1;
+            phase1_.alphaPhi() = alphaPhi1;
         }
 
         if (alphaDbyA.valid())
@@ -415,8 +391,8 @@ void Foam::twoPhaseSystem::solve()
         phase1_.alphaRhoPhi() =
             fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
 
-        phase2_.alphaPhi() = phi - phase1_.alphaPhi();
+        phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
-        alpha2 = scalar(1) - alpha1;
+        phase2_.correctInflowOutflow(phase2_.alphaPhi());
         phase2_.alphaRhoPhi() =
             fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
 
@@ -425,6 +401,13 @@ void Foam::twoPhaseSystem::solve()
             << "  Min(alpha1) = " << min(alpha1).value()
             << "  Max(alpha1) = " << max(alpha1).value()
             << endl;
+
+        // Ensure the phase-fractions are bounded
+        alpha1.max(0);
+        alpha1.min(1);
+
+        // Update the phase-fraction of the other phase
+        alpha2 = scalar(1) - alpha1;
     }
 }
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -248,27 +248,19 @@ bool Foam::phaseModel::read(const dictionary& phaseProperties)
 }
 
 
-void Foam::phaseModel::correctInflowFlux(surfaceScalarField& alphaPhi) const
+void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
 {
     surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+    const volScalarField::Boundary& alphaBf = boundaryField();
+    const surfaceScalarField::Boundary& phiBf = phi().boundaryField();
 
-    // Ensure that the flux at inflow BCs is preserved
     forAll(alphaPhiBf, patchi)
     {
         fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
 
         if (!alphaPhip.coupled())
         {
-            const scalarField& phip = phi().boundaryField()[patchi];
+            alphaPhip = phiBf[patchi]*alphaBf[patchi];
-            const scalarField& alphap = boundaryField()[patchi];
-
-            forAll(alphaPhip, facei)
-            {
-                if (phip[facei] < SMALL)
-                {
-                    alphaPhip[facei] = alphap[facei]*phip[facei];
-                }
-            }
         }
     }
 }

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -319,8 +319,8 @@ public:
             return alphaRhoPhi_;
         }
 
-        //- Ensure that the flux at inflow BCs is preserved
+        //- Ensure that the flux at inflow/outflow BCs is preserved
-        void correctInflowFlux(surfaceScalarField& alphaPhi) const;
+        void correctInflowOutflow(surfaceScalarField& alphaPhi) const;
 
         //- Correct the phase properties
         //  other than the thermodynamics and turbulence

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -444,7 +444,7 @@ void Foam::twoPhaseSystem::solve()
             )
         );
 
-        phase1_.correctInflowFlux(alphaPhic1);
+        phase1_.correctInflowOutflow(alphaPhic1);
 
         if (nAlphaSubCycles > 1)
         {
@@ -515,8 +515,7 @@ void Foam::twoPhaseSystem::solve()
             fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
 
         phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
-        alpha2 = scalar(1) - alpha1;
+        phase2_.correctInflowOutflow(phase2_.alphaPhi());
-        phase2_.correctInflowFlux(phase2_.alphaPhi());
         phase2_.alphaRhoPhi() =
             fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
 
@@ -525,6 +524,12 @@ void Foam::twoPhaseSystem::solve()
             << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
+
+        // Ensure the phase-fractions are bounded
+        alpha1.max(0);
+        alpha1.min(1);
+
+        alpha2 = scalar(1) - alpha1;
     }
 }
 

src/finiteVolume/Make/files

@@ -235,7 +235,6 @@ fvMatrices/fvMatrices.C
 fvMatrices/fvScalarMatrix/fvScalarMatrix.C
 fvMatrices/solvers/MULES/MULES.C
 fvMatrices/solvers/MULES/CMULES.C
-fvMatrices/solvers/MULES/IMULES.C
 fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
 
 interpolation = interpolation/interpolation

src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -151,17 +151,17 @@ void Foam::MULES::limiterCorr
 
     const dictionary& MULEScontrols = mesh.solverDict(psi.name());
 
-    label nLimiterIter
+    const label nLimiterIter
     (
         readLabel(MULEScontrols.lookup("nLimiterIter"))
     );
 
-    scalar smoothLimiter
+    const scalar smoothLimiter
     (
         MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
     );
 
-    scalar extremaCoeff
+    const scalar extremaCoeff
     (
         MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0)
     );
@@ -207,8 +207,8 @@ void Foam::MULES::limiterCorr
 
     forAll(phiCorrIf, facei)
     {
-        label own = owner[facei];
+        const label own = owner[facei];
-        label nei = neighb[facei];
+        const label nei = neighb[facei];
 
         psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
         psiMinn[own] = min(psiMinn[own], psiIf[nei]);
@@ -216,9 +216,9 @@ void Foam::MULES::limiterCorr
         psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
         psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
 
-        scalar phiCorrf = phiCorrIf[facei];
+        const scalar phiCorrf = phiCorrIf[facei];
 
-        if (phiCorrf > 0.0)
+        if (phiCorrf > 0)
         {
             sumPhip[own] += phiCorrf;
             mSumPhim[nei] += phiCorrf;
@@ -253,7 +253,7 @@ void Foam::MULES::limiterCorr
         {
             forAll(phiCorrPf, pFacei)
             {
-                label pfCelli = pFaceCells[pFacei];
+                const label pfCelli = pFaceCells[pFacei];
 
                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
@@ -262,11 +262,11 @@ void Foam::MULES::limiterCorr
 
         forAll(phiCorrPf, pFacei)
         {
-            label pfCelli = pFaceCells[pFacei];
+            const label pfCelli = pFaceCells[pFacei];
 
-            scalar phiCorrf = phiCorrPf[pFacei];
+            const scalar phiCorrf = phiCorrPf[pFacei];
 
-            if (phiCorrf > 0.0)
+            if (phiCorrf > 0)
             {
                 sumPhip[pfCelli] += phiCorrf;
             }
@@ -309,17 +309,17 @@ void Foam::MULES::limiterCorr
 
     for (int j=0; j<nLimiterIter; j++)
     {
-        sumlPhip = 0.0;
+        sumlPhip = 0;
-        mSumlPhim = 0.0;
+        mSumlPhim = 0;
 
         forAll(lambdaIf, facei)
         {
-            label own = owner[facei];
+            const label own = owner[facei];
-            label nei = neighb[facei];
+            const label nei = neighb[facei];
 
-            scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
+            const scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
 
-            if (lambdaPhiCorrf > 0.0)
+            if (lambdaPhiCorrf > 0)
             {
                 sumlPhip[own] += lambdaPhiCorrf;
                 mSumlPhim[nei] += lambdaPhiCorrf;
@@ -344,7 +344,7 @@ void Foam::MULES::limiterCorr
 
                 scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
 
-                if (lambdaPhiCorrf > 0.0)
+                if (lambdaPhiCorrf > 0)
                 {
                     sumlPhip[pfCelli] += lambdaPhiCorrf;
                 }
@@ -379,7 +379,7 @@ void Foam::MULES::limiterCorr
 
         forAll(lambdaIf, facei)
         {
-            if (phiCorrIf[facei] > 0.0)
+            if (phiCorrIf[facei] > 0)
             {
                 lambdaIf[facei] = min
                 (
@@ -415,9 +415,9 @@ void Foam::MULES::limiterCorr
 
                 forAll(lambdaPf, pFacei)
                 {
-                    label pfCelli = pFaceCells[pFacei];
+                    const label pfCelli = pFaceCells[pFacei];
 
-                    if (phiCorrfPf[pFacei] > 0.0)
+                    if (phiCorrfPf[pFacei] > 0)
                     {
                         lambdaPf[pFacei] =
                             min(lambdaPf[pFacei], lambdap[pfCelli]);
@@ -438,11 +438,11 @@ void Foam::MULES::limiterCorr
                 forAll(lambdaPf, pFacei)
                 {
                     // Limit outlet faces only
-                    if (phiPf[pFacei] > SMALL*SMALL)
+                    if ((phiPf[pFacei] + phiCorrfPf[pFacei]) > SMALL*SMALL)
                     {
-                        label pfCelli = pFaceCells[pFacei];
+                        const label pfCelli = pFaceCells[pFacei];
 
-                        if (phiCorrfPf[pFacei] > 0.0)
+                        if (phiCorrfPf[pFacei] > 0)
                         {
                             lambdaPf[pFacei] =
                                 min(lambdaPf[pFacei], lambdap[pfCelli]);

src/finiteVolume/fvMatrices/solvers/MULES/IMULES.C

@@ -1,51 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "IMULES.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::MULES::implicitSolve
-(
-    volScalarField& psi,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiPsi,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    implicitSolve
-    (
-        geometricOneField(),
-        psi,
-        phi,
-        phiPsi,
-        zeroField(), zeroField(),
-        psiMax, psiMin
-    );
-}
-
-
-// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H

@@ -1,94 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Global
-    IMULES
-
-Description
-    IMULES: Multidimensional universal limiter for implicit solution.
-
-    Solve a convective-only transport equation using an explicit universal
-    multi-dimensional limiter applied to an implicit formulation requiring
-    iteration to guarantee boundedness.  The number of iterations required
-    to obtain boundedness increases with the Courant number of the simulation.
-
-    It may be more efficient to use CMULES.
-
-SourceFiles
-    IMULES.C
-    IMULESTemplates.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef IMULES_H
-#define IMULES_H
-
-#include "MULES.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace MULES
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class RhoType, class SpType, class SuType>
-void implicitSolve
-(
-    const RhoType& rho,
-    volScalarField& gamma,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-void implicitSolve
-(
-    volScalarField& gamma,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiCorr,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace MULES
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "IMULESTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C

@@ -1,236 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "IMULES.H"
-#include "gaussConvectionScheme.H"
-#include "surfaceInterpolate.H"
-#include "fvmDdt.H"
-#include "fvmSup.H"
-#include "fvcDiv.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace MULES
-{
-    template<class RhoType>
-    inline tmp<surfaceScalarField> interpolate(const RhoType& rho)
-    {
-        NotImplemented;
-        return tmp<surfaceScalarField>(nullptr);
-    }
-
-    template<>
-    inline tmp<surfaceScalarField> interpolate(const volScalarField& rho)
-    {
-        return fvc::interpolate(rho);
-    }
-}
-}
-
-
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::implicitSolve
-(
-    const RhoType& rho,
-    volScalarField& psi,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiPsi,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    const fvMesh& mesh = psi.mesh();
-
-    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
-
-    label maxIter
-    (
-        readLabel(MULEScontrols.lookup("maxIter"))
-    );
-
-    scalar maxUnboundedness
-    (
-        readScalar(MULEScontrols.lookup("maxUnboundedness"))
-    );
-
-    scalar CoCoeff
-    (
-        readScalar(MULEScontrols.lookup("CoCoeff"))
-    );
-
-    scalarField allCoLambda(mesh.nFaces());
-
-    {
-        slicedSurfaceScalarField CoLambda
-        (
-            IOobject
-            (
-                "CoLambda",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimless,
-            allCoLambda,
-            false   // Use slices for the couples
-        );
-
-        if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
-        {
-            tmp<surfaceScalarField> Cof =
-                mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
-               *mag(phi/interpolate(rho))/mesh.magSf();
-
-            CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
-        }
-        else
-        {
-            tmp<surfaceScalarField> Cof =
-                mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
-               *mag(phi)/mesh.magSf();
-
-            CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
-        }
-    }
-
-    scalarField allLambda(allCoLambda);
-    //scalarField allLambda(mesh.nFaces(), 1.0);
-
-    slicedSurfaceScalarField lambda
-    (
-        IOobject
-        (
-            "lambda",
-            mesh.time().timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh,
-        dimless,
-        allLambda,
-        false   // Use slices for the couples
-    );
-
-    linear<scalar> CDs(mesh);
-    upwind<scalar> UDs(mesh, phi);
-    //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
-
-    fvScalarMatrix psiConvectionDiffusion
-    (
-        fvm::ddt(rho, psi)
-      + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
-        //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
-        //.fvmLaplacian(Dpsif, psi)
-      - fvm::Sp(Sp, psi)
-      - Su
-    );
-
-    surfaceScalarField phiBD(psiConvectionDiffusion.flux());
-
-    surfaceScalarField& phiCorr = phiPsi;
-    phiCorr -= phiBD;
-
-    for (label i=0; i<maxIter; i++)
-    {
-        if (i != 0 && i < 4)
-        {
-            allLambda = allCoLambda;
-        }
-
-        limiter
-        (
-            allLambda,
-            1.0/mesh.time().deltaTValue(),
-            rho,
-            psi,
-            phiBD,
-            phiCorr,
-            Sp,
-            Su,
-            psiMax,
-            psiMin
-        );
-
-        solve
-        (
-            psiConvectionDiffusion + fvc::div(lambda*phiCorr),
-            MULEScontrols
-        );
-
-        scalar maxPsiM1 = gMax(psi.primitiveField()) - 1.0;
-        scalar minPsi = gMin(psi.primitiveField());
-
-        scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
-
-        if (unboundedness < maxUnboundedness)
-        {
-            break;
-        }
-        else
-        {
-            Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
-                << " min(" << psi.name() << ") = " << minPsi << endl;
-
-            phiBD = psiConvectionDiffusion.flux();
-
-            /*
-            word gammaScheme("div(phi,gamma)");
-            word gammarScheme("div(phirb,gamma)");
-
-            const surfaceScalarField& phir =
-                mesh.lookupObject<surfaceScalarField>("phir");
-
-            phiCorr =
-                fvc::flux
-                (
-                    phi,
-                    psi,
-                    gammaScheme
-                )
-              + fvc::flux
-                (
-                    -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
-                    psi,
-                    gammarScheme
-                )
-                - phiBD;
-            */
-        }
-    }
-
-    phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
-}
-
-
-// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -185,12 +185,12 @@ void Foam::MULES::limiter
 
     const dictionary& MULEScontrols = mesh.solverDict(psi.name());
 
-    label nLimiterIter
+    const label nLimiterIter
     (
         MULEScontrols.lookupOrDefault<label>("nLimiterIter", 3)
     );
 
-    scalar smoothLimiter
+    const scalar smoothLimiter
     (
         MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
     );
@@ -228,8 +228,7 @@ void Foam::MULES::limiter
     );
 
     scalarField& lambdaIf = lambda;
-    surfaceScalarField::Boundary& lambdaBf =
+    surfaceScalarField::Boundary& lambdaBf = lambda.boundaryFieldRef();
-        lambda.boundaryFieldRef();
 
     scalarField psiMaxn(psiIf.size(), psiMin);
     scalarField psiMinn(psiIf.size(), psiMax);
@@ -241,8 +240,8 @@ void Foam::MULES::limiter
 
     forAll(phiCorrIf, facei)
     {
-        label own = owner[facei];
+        const label own = owner[facei];
-        label nei = neighb[facei];
+        const label nei = neighb[facei];
 
         psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
         psiMinn[own] = min(psiMinn[own], psiIf[nei]);
@@ -253,9 +252,9 @@ void Foam::MULES::limiter
         sumPhiBD[own] += phiBDIf[facei];
         sumPhiBD[nei] -= phiBDIf[facei];
 
-        scalar phiCorrf = phiCorrIf[facei];
+        const scalar phiCorrf = phiCorrIf[facei];
 
-        if (phiCorrf > 0.0)
+        if (phiCorrf > 0)
         {
             sumPhip[own] += phiCorrf;
             mSumPhim[nei] += phiCorrf;
@@ -281,7 +280,7 @@ void Foam::MULES::limiter
 
             forAll(phiCorrPf, pFacei)
             {
-                label pfCelli = pFaceCells[pFacei];
+                const label pfCelli = pFaceCells[pFacei];
 
                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
@@ -291,7 +290,7 @@ void Foam::MULES::limiter
         {
             forAll(phiCorrPf, pFacei)
             {
-                label pfCelli = pFaceCells[pFacei];
+                const label pfCelli = pFaceCells[pFacei];
 
                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
@@ -300,13 +299,13 @@ void Foam::MULES::limiter
 
         forAll(phiCorrPf, pFacei)
         {
-            label pfCelli = pFaceCells[pFacei];
+            const label pfCelli = pFaceCells[pFacei];
 
             sumPhiBD[pfCelli] += phiBDPf[pFacei];
 
-            scalar phiCorrf = phiCorrPf[pFacei];
+            const scalar phiCorrf = phiCorrPf[pFacei];
 
-            if (phiCorrf > 0.0)
+            if (phiCorrf > 0)
             {
                 sumPhip[pfCelli] += phiCorrf;
             }
@@ -376,17 +375,17 @@ void Foam::MULES::limiter
 
     for (int j=0; j<nLimiterIter; j++)
     {
-        sumlPhip = 0.0;
+        sumlPhip = 0;
-        mSumlPhim = 0.0;
+        mSumlPhim = 0;
 
         forAll(lambdaIf, facei)
         {
-            label own = owner[facei];
+            const label own = owner[facei];
-            label nei = neighb[facei];
+            const label nei = neighb[facei];
 
             scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
 
-            if (lambdaPhiCorrf > 0.0)
+            if (lambdaPhiCorrf > 0)
             {
                 sumlPhip[own] += lambdaPhiCorrf;
                 mSumlPhim[nei] += lambdaPhiCorrf;
@@ -407,11 +406,11 @@ void Foam::MULES::limiter
 
             forAll(lambdaPf, pFacei)
             {
-                label pfCelli = pFaceCells[pFacei];
+                const label pfCelli = pFaceCells[pFacei];
+                const scalar lambdaPhiCorrf =
+                    lambdaPf[pFacei]*phiCorrfPf[pFacei];
 
-                scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
+                if (lambdaPhiCorrf > 0)
-
-                if (lambdaPhiCorrf > 0.0)
                 {
                     sumlPhip[pfCelli] += lambdaPhiCorrf;
                 }
@@ -446,7 +445,7 @@ void Foam::MULES::limiter
 
         forAll(lambdaIf, facei)
         {
-            if (phiCorrIf[facei] > 0.0)
+            if (phiCorrIf[facei] > 0)
             {
                 lambdaIf[facei] = min
                 (
@@ -464,7 +463,6 @@ void Foam::MULES::limiter
             }
         }
 
-
         forAll(lambdaBf, patchi)
         {
             fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
@@ -482,9 +480,9 @@ void Foam::MULES::limiter
 
                 forAll(lambdaPf, pFacei)
                 {
-                    label pfCelli = pFaceCells[pFacei];
+                    const label pfCelli = pFaceCells[pFacei];
 
-                    if (phiCorrfPf[pFacei] > 0.0)
+                    if (phiCorrfPf[pFacei] > 0)
                     {
                         lambdaPf[pFacei] =
                             min(lambdaPf[pFacei], lambdap[pfCelli]);
@@ -496,33 +494,6 @@ void Foam::MULES::limiter
                     }
                 }
             }
-            else
-            {
-                const labelList& pFaceCells =
-                    mesh.boundary()[patchi].faceCells();
-                const scalarField& phiBDPf = phiBDBf[patchi];
-                const scalarField& phiCorrPf = phiCorrBf[patchi];
-
-                forAll(lambdaPf, pFacei)
-                {
-                    // Limit outlet faces only
-                    if ((phiBDPf[pFacei] + phiCorrPf[pFacei]) > SMALL*SMALL)
-                    {
-                        label pfCelli = pFaceCells[pFacei];
-
-                        if (phiCorrfPf[pFacei] > 0.0)
-                        {
-                            lambdaPf[pFacei] =
-                                min(lambdaPf[pFacei], lambdap[pfCelli]);
-                        }
-                        else
-                        {
-                            lambdaPf[pFacei] =
-                                min(lambdaPf[pFacei], lambdam[pfCelli]);
-                        }
-                    }
-                }
-            }
         }
 
         syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
@@ -549,6 +520,19 @@ void Foam::MULES::limit
 
     surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
 
+    surfaceScalarField::Boundary& phiBDBf = phiBD.boundaryFieldRef();
+    const surfaceScalarField::Boundary& phiPsiBf = phiPsi.boundaryField();
+
+    forAll(phiBDBf, patchi)
+    {
+        fvsPatchScalarField& phiBDPf = phiBDBf[patchi];
+
+        if (!phiBDPf.coupled())
+        {
+            phiBDPf = phiPsiBf[patchi];
+        }
+    }
+
     surfaceScalarField& phiCorr = phiPsi;
     phiCorr -= phiBD;
 

tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution

@@ -20,7 +20,7 @@ solvers
     "alpha.*"
     {
         nAlphaSubCycles 4;
-        cAlpha          2;
+        cAlpha          1;
     }
 
     pcorr

tutorials/multiphase/interDyMFoam/RAS/testTubeMixer/system/fvSolution

@@ -21,7 +21,7 @@ solvers
     {
         nAlphaCorr      1;
         nAlphaSubCycles 3;
-        cAlpha          1.5;
+        cAlpha          1;
     }
 
     pcorr

tutorials/multiphase/interFoam/LES/nozzleFlow2D/system/fvSolution

@@ -21,7 +21,7 @@ solvers
     {
         nAlphaCorr      1;
         nAlphaSubCycles 4;
-        cAlpha          2;
+        cAlpha          1;
     }
 
     pcorr

tutorials/multiphase/interFoam/laminar/damBreak/damBreak/system/fvSolution

@@ -24,7 +24,7 @@ solvers
         cAlpha          1;
 
         MULESCorr       yes;
-        nLimiterIter    3;
+        nLimiterIter    5;
 
         solver          smoothSolver;
         smoother        symGaussSeidel;

tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSolution

@@ -20,7 +20,7 @@ solvers
     "alpha.*"
     {
         nAlphaSubCycles 4;
-        cAlpha          2;
+        cAlpha          1;
     }
 
     pcorr

tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSolution

@@ -20,7 +20,7 @@ solvers
     "alpha.*"
     {
         nAlphaSubCycles 4;
-        cAlpha          2;
+        cAlpha          1;
     }
 
     pcorr

tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/system/fvSolution

@@ -21,7 +21,7 @@ solvers
     {
         nAlphaCorr      4;
         nAlphaSubCycles 4;
-        cAlpha          2;
+        cAlpha          1;
     }
 
     pcorr