MULES: nLimiterIter and smoothLimiter are now user-input via the corresponding fvSolution sub-dict

nLimiterIter: Number of iterations during limiter construction
    3 (default) is sufficient for 3D simulations with a Courant number 0.5 or so
    For larger Courant numbers larger values may be needed but this is
    only relevant for IMULES and CMULES

smoothLimiter: Coefficient to smooth the limiter to avoid "diamond"
    staggering patters seen in regions of low particle phase-fraction in
    fluidised-bed simulations.

    The default is 0 as it is not needed for all simulations.
    A value of 0.1 is appropriate for fluidised-bed simulations.
    The useful range is 0 -> 0.5.
    Values larger than 0.5 may cause excessive smearing of the solution.

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C

@@ -996,7 +996,6 @@ void Foam::multiphaseMixtureThermo::solveAlphas
             zeroField(),
             1,
             0,
-            3,
             true
         );
 

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H

@@ -83,8 +83,7 @@
             zeroField(),
             zeroField(),
             1,
-            0,
+            0
-            3
         );
 
         // Create the complete flux for alpha2
@@ -125,8 +124,7 @@
             zeroField(),
             zeroField(),
             1,
-            0,
+            0
-            3
         );
 
         // Construct the limited fluxes

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

@@ -158,7 +158,6 @@ void Foam::multiphaseSystem::solveAlphas()
             zeroField(),
             1,
             0,
-            3,
             true
         );
 

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

@@ -129,8 +129,7 @@ void limiterCorr
     const SpType& Sp,
     const SuType& Su,
     const scalar psiMax,
-    const scalar psiMin,
+    const scalar psiMin
-    const label nLimiterIter
 );
 
 template<class RdeltaTType, class RhoType, class SpType, class SuType>
@@ -144,8 +143,7 @@ void limitCorr
     const SpType& Sp,
     const SuType& Su,
     const scalar psiMax,
-    const scalar psiMin,
+    const scalar psiMin
-    const label nLimiterIter
 );
 
 

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

@@ -89,13 +89,6 @@ void Foam::MULES::correct
     const fvMesh& mesh = psi.mesh();
     const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
 
-    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
-
-    label nLimiterIter
-    (
-        readLabel(MULEScontrols.lookup("nLimiterIter"))
-    );
-
     limitCorr
     (
         rDeltaT,
@@ -103,9 +96,10 @@ void Foam::MULES::correct
         psi,
         phi,
         phiCorr,
-        Sp, Su,
+        Sp,
-        psiMax, psiMin,
+        Su,
-        nLimiterIter
+        psiMax,
+        psiMin
     );
 
     correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
@@ -130,13 +124,6 @@ void Foam::MULES::LTScorrect
     const volScalarField& rDeltaT =
         mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
 
-    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
-
-    label nLimiterIter
-    (
-        readLabel(MULEScontrols.lookup("nLimiterIter"))
-    );
-
     limitCorr
     (
         rDeltaT,
@@ -144,9 +131,10 @@ void Foam::MULES::LTScorrect
         psi,
         phi,
         phiCorr,
-        Sp, Su,
+        Sp,
-        psiMax, psiMin,
+        Su,
-        nLimiterIter
+        psiMax,
+        psiMin
     );
     correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
 }
@@ -164,8 +152,7 @@ void Foam::MULES::limiterCorr
     const SpType& Sp,
     const SuType& Su,
     const scalar psiMax,
-    const scalar psiMin,
+    const scalar psiMin
-    const label nLimiterIter
 )
 {
     const scalarField& psiIf = psi;
@@ -175,9 +162,19 @@ void Foam::MULES::limiterCorr
 
     const dictionary& MULEScontrols = mesh.solverDict(psi.name());
 
+    label nLimiterIter
+    (
+        readLabel(MULEScontrols.lookup("nLimiterIter"))
+    );
+
+    scalar smoothLimiter
+    (
+        MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
+    );
+
     scalar extremaCoeff
     (
-        MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0.0)
+        MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0)
     );
 
     const labelUList& owner = mesh.owner();
@@ -294,9 +291,13 @@ void Foam::MULES::limiterCorr
     psiMaxn = min(psiMaxn + extremaCoeff*(psiMax - psiMin), psiMax);
     psiMinn = max(psiMinn - extremaCoeff*(psiMax - psiMin), psiMin);
 
-    // scalar smooth = 0.5;
+    if (smoothLimiter > SMALL)
-    // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
+    {
-    // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
+        psiMaxn =
+            min(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMaxn, psiMax);
+        psiMinn =
+            max(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMinn, psiMin);
+    }
 
     psiMaxn =
         V
@@ -483,8 +484,7 @@ void Foam::MULES::limitCorr
     const SpType& Sp,
     const SuType& Su,
     const scalar psiMax,
-    const scalar psiMin,
+    const scalar psiMin
-    const label nLimiterIter
 )
 {
     const fvMesh& mesh = psi.mesh();
@@ -519,8 +519,7 @@ void Foam::MULES::limitCorr
         Sp,
         Su,
         psiMax,
-        psiMin,
+        psiMin
-        nLimiterIter
     );
 
     phiCorr *= lambda;

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

@@ -77,11 +77,6 @@ void Foam::MULES::implicitSolve
         readLabel(MULEScontrols.lookup("maxIter"))
     );
 
-    label nLimiterIter
-    (
-        readLabel(MULEScontrols.lookup("nLimiterIter"))
-    );
-
     scalar maxUnboundedness
     (
         readScalar(MULEScontrols.lookup("maxUnboundedness"))
@@ -187,8 +182,7 @@ void Foam::MULES::implicitSolve
             Sp,
             Su,
             psiMax,
-            psiMin,
+            psiMin
-            nLimiterIter
         );
 
         solve

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

@@ -138,8 +138,7 @@ void limiter
     const SpType& Sp,
     const SuType& Su,
     const scalar psiMax,
-    const scalar psiMin,
+    const scalar psiMin
-    const label nLimiterIter
 );
 
 template<class RdeltaTType, class RhoType, class SpType, class SuType>
@@ -154,7 +153,6 @@ void limit
     const SuType& Su,
     const scalar psiMax,
     const scalar psiMin,
-    const label nLimiterIter,
     const bool returnCorr
 );
 

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

@@ -108,8 +108,23 @@ void Foam::MULES::explicitSolve
 {
     const fvMesh& mesh = psi.mesh();
     const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
+
     psi.correctBoundaryConditions();
-    limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
+
+    limit
+    (
+        rDeltaT,
+        rho,
+        psi,
+        phi,
+        phiPsi,
+        Sp,
+        Su,
+        psiMax,
+        psiMin,
+        false
+    );
+
     explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
 }
 
@@ -133,7 +148,21 @@ void Foam::MULES::explicitLTSSolve
         mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
 
     psi.correctBoundaryConditions();
-    limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
+
+    limit
+    (
+        rDeltaT,
+        rho,
+        psi,
+        phi,
+        phiPsi,
+        Sp,
+        Su,
+        psiMax,
+        psiMin,
+        false
+    );
+
     explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
 }
 
@@ -150,17 +179,28 @@ void Foam::MULES::limiter
     const SpType& Sp,
     const SuType& Su,
     const scalar psiMax,
-    const scalar psiMin,
+    const scalar psiMin
-    const label nLimiterIter
 )
 {
     const scalarField& psiIf = psi;
     const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
 
-    const scalarField& psi0 = psi.oldTime();
-
     const fvMesh& mesh = psi.mesh();
 
+    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+    label nLimiterIter
+    (
+        MULEScontrols.lookupOrDefault<label>("nLimiterIter", 3)
+    );
+
+    scalar smoothLimiter
+    (
+        MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
+    );
+
+    const scalarField& psi0 = psi.oldTime();
+
     const labelUList& owner = mesh.owner();
     const labelUList& neighb = mesh.neighbour();
     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
@@ -284,9 +324,13 @@ void Foam::MULES::limiter
     psiMaxn = min(psiMaxn, psiMax);
     psiMinn = max(psiMinn, psiMin);
 
-    //scalar smooth = 0.5;
+    if (smoothLimiter > SMALL)
-    //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
+    {
-    //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
+        psiMaxn =
+            min(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMaxn, psiMax);
+        psiMinn =
+            max(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMinn, psiMin);
+    }
 
     if (mesh.moving())
     {
@@ -502,7 +546,6 @@ void Foam::MULES::limit
     const SuType& Su,
     const scalar psiMax,
     const scalar psiMin,
-    const label nLimiterIter,
     const bool returnCorr
 )
 {
@@ -543,8 +586,7 @@ void Foam::MULES::limit
         Sp,
         Su,
         psiMax,
-        psiMin,
+        psiMin
-        nLimiterIter
     );
 
     if (returnCorr)

tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    "alpha.*"
     {
         nAlphaSubCycles 2;
     }

tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    "alpha.*"
     {
         nAlphaSubCycles 3;
     }

tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    "alpha.*"
     {
         nAlphaSubCycles 3;
     }

tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    "alpha.*"
     {
         nAlphaSubCycles 2;
     }

tutorials/multiphase/multiphaseInterDyMFoam/laminar/mixerVesselAMI2D/system/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    "alpha.*"
     {
         nAlphaSubCycles 1;
         cAlpha          1;

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

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    "alpha.*"
     {
         nAlphaSubCycles 4;
         cAlpha          2;

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

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    "alpha.*"
     {
         nAlphaSubCycles 4;
         cAlpha          2;

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

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    "alpha.*"
     {
         nAlphaCorr      4;
         nAlphaSubCycles 4;

tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution

@@ -21,8 +21,10 @@ solvers
     {
         nAlphaCorr      1;
         nAlphaSubCycles 2;
-        implicitPhasePressure yes;
 
+        smoothLimiter   0.1;
+
+        implicitPhasePressure yes;
         solver          smoothSolver;
         smoother        symGaussSeidel;
         tolerance       1e-9;

tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution

@@ -21,8 +21,10 @@ solvers
     {
         nAlphaCorr      1;
         nAlphaSubCycles 2;
-        implicitPhasePressure yes;
 
+        smoothLimiter   0.1;
+
+        implicitPhasePressure yes;
         solver          smoothSolver;
         smoother        symGaussSeidel;
         tolerance       1e-9;