MMULES: new development of the MULES limiter to support limiting the sum of phase-fractions

multiphaseInterFoam: Upgraded to use the new MMULES algorithm

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2006-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,4 +80,60 @@ void Foam::MULES::implicitSolve
 }
 
 
+void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
+{
+    forAll(phiPsiCorrs[0], facei)
+    {
+        scalar sumPos = 0;
+        scalar sumNeg = 0;
+
+        for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
+        {
+            if (phiPsiCorrs[phasei][facei] > 0)
+            {
+                sumPos += phiPsiCorrs[phasei][facei];
+            }
+            else
+            {
+                sumNeg += phiPsiCorrs[phasei][facei];
+            }
+        }
+
+        scalar sum = sumPos + sumNeg;
+
+        if (sum > 0 && sumPos > VSMALL)
+        {
+            scalar lambda = -sumNeg/sumPos;
+
+            for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
+            {
+                if (phiPsiCorrs[phasei][facei] > 0)
+                {
+                    phiPsiCorrs[phasei][facei] *= lambda;
+                }
+            }
+        }
+        else if (sum < 0 && sumNeg < -VSMALL)
+        {
+            scalar lambda = -sumPos/sumNeg;
+
+            for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
+            {
+                if (phiPsiCorrs[phasei][facei] < 0)
+                {
+                    phiPsiCorrs[phasei][facei] *= lambda;
+                }
+            }
+        }
+        else
+        {
+            for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
+            {
+                phiPsiCorrs[phasei][facei] = 0;
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2006-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,18 +45,27 @@ SourceFiles
 #include "volFields.H"
 #include "surfaceFieldsFwd.H"
 #include "primitiveFieldsFwd.H"
-#include "zero.H"
 #include "geometricOneField.H"
+#include "zero.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+namespace MULES
+{
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace MULES
-{
+template<class RhoType, class SpType, class SuType>
+void explicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su
+);
 
 template<class RhoType, class SpType, class SuType>
 void explicitSolve
@@ -117,10 +126,29 @@ void limiter
     const label nLimiterIter
 );
 
-} // End namespace MULES
+template<class RhoType, class SpType, class SuType>
+void limit
+(
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter,
+    const bool returnCorr
+);
+
+void limitSum(UPtrList<scalarField>& phiPsiCorrs);
+
+template<class SurfaceScalarFieldList>
+void limitSum(SurfaceScalarFieldList& phiPsiCorrs);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace MULES
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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) 2004-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2006-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,58 +43,14 @@ void Foam::MULES::explicitSolve
 (
     const RhoType& rho,
     volScalarField& psi,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiPsi,
+    const surfaceScalarField& phiPsi,
     const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const SuType& Su
 )
 {
     Info<< "MULES: Solving for " << psi.name() << endl;
 
     const fvMesh& mesh = psi.mesh();
-    psi.correctBoundaryConditions();
-
-    surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
-
-    surfaceScalarField& phiCorr = phiPsi;
-    phiCorr -= phiBD;
-
-    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
-    );
-
-    limiter
-    (
-        allLambda,
-        rho,
-        psi,
-        phiBD,
-        phiCorr,
-        Sp,
-        Su,
-        psiMax,
-        psiMin,
-        3
-    );
-
-    phiPsi = phiBD + lambda*phiCorr;
 
     scalarField& psiIf = psi;
     const scalarField& psi0 = psi.oldTime();
@@ -127,6 +83,25 @@ void Foam::MULES::explicitSolve
 }
 
 
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::explicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    psi.correctBoundaryConditions();
+    limit(rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
+    explicitSolve(rho, psi, phiPsi, Sp, Su);
+}
+
+
 template<class RhoType, class SpType, class SuType>
 void Foam::MULES::implicitSolve
 (
@@ -614,4 +589,100 @@ void Foam::MULES::limiter
 }
 
 
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::limit
+(
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter,
+    const bool returnCorr
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
+
+    surfaceScalarField& phiCorr = phiPsi;
+    phiCorr -= phiBD;
+
+    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
+    );
+
+    limiter
+    (
+        allLambda,
+        rho,
+        psi,
+        phiBD,
+        phiCorr,
+        Sp,
+        Su,
+        psiMax,
+        psiMin,
+        nLimiterIter
+    );
+
+    if (returnCorr)
+    {
+        phiCorr *= lambda;
+    }
+    else
+    {
+        phiPsi = phiBD + lambda*phiCorr;
+    }
+}
+
+
+template<class SurfaceScalarFieldList>
+void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
+{
+    {
+        UPtrList<scalarField> phiPsiCorrsInternal(phiPsiCorrs.size());
+        forAll(phiPsiCorrs, phasei)
+        {
+            phiPsiCorrsInternal.set(phasei, &phiPsiCorrs[phasei]);
+        }
+
+        limitSum(phiPsiCorrsInternal);
+    }
+
+    forAll(phiPsiCorrs[0].boundaryField(), patchi)
+    {
+        UPtrList<scalarField> phiPsiCorrsPatch(phiPsiCorrs.size());
+        forAll(phiPsiCorrs, phasei)
+        {
+            phiPsiCorrsPatch.set
+            (
+                phasei,
+                &phiPsiCorrs[phasei].boundaryField()[patchi]
+            );
+        }
+
+        limitSum(phiPsiCorrsPatch);
+    }
+}
+
+
 // ************************************************************************* //