MPPICFoam PackingModels: Implicit: Added limiting to the calculation of the correction flux.

Vastly reduces the scattering and churning behaviour of packed beds.

Development provided by Will Bainbridge <github.com/will-bainbridge>

See also http://www.openfoam.org/mantisbt/view.php?id=1994

src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C

@@ -48,6 +48,7 @@ Foam::PackingModels::Implicit<CloudType>::Implicit
     ),
     phiCorrect_(NULL),
     uCorrect_(NULL),
+    applyLimiting_(this->coeffDict().lookup("applyLimiting")),
     applyGravity_(this->coeffDict().lookup("applyGravity")),
     alphaMin_(readScalar(this->coeffDict().lookup("alphaMin"))),
     rhoMin_(readScalar(this->coeffDict().lookup("rhoMin")))
@@ -66,6 +67,7 @@ Foam::PackingModels::Implicit<CloudType>::Implicit
     alpha_(cm.alpha_),
     phiCorrect_(cm.phiCorrect_()),
     uCorrect_(cm.uCorrect_()),
+    applyLimiting_(cm.applyLimiting_),
     applyGravity_(cm.applyGravity_),
     alphaMin_(cm.alphaMin_),
     rhoMin_(cm.rhoMin_)
@@ -102,6 +104,11 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
             (
                 cloudName + ":rhoAverage"
             );
+        const AveragingMethod<vector>& uAverage =
+            mesh.lookupObject<AveragingMethod<vector> >
+            (
+                cloudName + ":uAverage"
+            );
         const AveragingMethod<scalar>& uSqrAverage =
             mesh.lookupObject<AveragingMethod<scalar>>
             (
@@ -135,13 +142,6 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
         rho.internalField() = max(rhoAverage.internalField(), rhoMin_);
         rho.correctBoundaryConditions();
 
-        //Info << "       x: " << mesh.C().internalField().component(2) << endl;
-        //Info << "   alpha: " << alpha_.internalField() << endl;
-        //Info << "alphaOld: " << alpha_.oldTime().internalField() << endl;
-        //Info << "     rho: " << rho.internalField() << endl;
-        //Info << endl;
-
-
         // Stress field
         // ~~~~~~~~~~~~
 
@@ -172,6 +172,24 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
         tauPrime.correctBoundaryConditions();
 
 
+        // Gravity flux
+        // ~~~~~~~~~~~~
+
+        tmp<surfaceScalarField> phiGByA;
+
+        if (applyGravity_)
+        (
+            phiGByA = tmp<surfaceScalarField>
+            (
+                new surfaceScalarField
+                (
+                    "phiGByA",
+                    deltaT*(g & mesh.Sf())*fvc::interpolate(1.0 - rhoc/rho)
+                )
+            )
+        );
+
+
         // Implicit solution for the volume fraction
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -191,14 +209,7 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
 
         if (applyGravity_)
         {
-            surfaceScalarField
+            alphaEqn += fvm::div(phiGByA(), alpha_);
-                phiGByA
-                (
-                    "phiGByA",
-                    deltaT*(g & mesh.Sf())*fvc::interpolate(1.0 - rhoc/rho)
-                );
-
-            alphaEqn += fvm::div(phiGByA, alpha_);
         }
 
         alphaEqn.solve();
@@ -217,6 +228,67 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
             )
         );
 
+        // limit the correction flux
+        if (applyLimiting_)
+        {
+            volVectorField U
+            (
+                IOobject
+                (
+                    cloudName + ":U",
+                    this->owner().db().time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedVector("zero", dimVelocity, vector::zero),
+                fixedValueFvPatchField<vector>::typeName
+            );
+            U.internalField() = uAverage.internalField();
+            U.correctBoundaryConditions();
+
+            surfaceScalarField phi
+            (
+                cloudName + ":phi",
+                linearInterpolate(U) & mesh.Sf()
+            );
+
+            if (applyGravity_)
+            {
+                phiCorrect_() -= phiGByA();
+            }
+
+            forAll(phiCorrect_(), faceI)
+            {
+                // Current and correction fluxes
+                const scalar phiCurr = phi[faceI];
+                scalar& phiCorr = phiCorrect_()[faceI];
+
+                // Don't limit if the correction is in the opposite direction to
+                // the flux. We need all the help we can get in this state.
+                if (phiCurr*phiCorr < 0)
+                {}
+
+                // If the correction and the flux are in the same direction then
+                // don't apply any more correction than is already present in
+                // the flux.
+                else if (phiCorr > 0)
+                {
+                    phiCorr = max(phiCorr - phiCurr, 0);
+                }
+                else
+                {
+                    phiCorr = min(phiCorr - phiCurr, 0);
+                }
+            }
+
+            if (applyGravity_)
+            {
+                phiCorrect_() += phiGByA();
+            }
+        }
+
         // correction velocity
         uCorrect_ = tmp<volVectorField>
         (

src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.H

@@ -58,8 +58,6 @@ class Implicit
 :
     public PackingModel<CloudType>
 {
-private:
-
     //- Private data
 
         //- Volume fraction field
@@ -71,6 +69,9 @@ private:
         //- Correction cell-centred velocity
         tmp<volVectorField> uCorrect_;
 
+        //- Flag to indicate whether implicit limiting is applied
+        Switch applyLimiting_;
+
         //- Flag to indicate whether gravity is applied
         Switch applyGravity_;
 

tutorials/lagrangian/MPPICFoam/Goldschmidt/constant/kinematicCloudProperties

@@ -148,6 +148,7 @@ subModels
     {
         alphaMin 0.0001;
         rhoMin 1.0;
+        applyLimiting true;
         applyGravity false;
         particleStressModel
         {

tutorials/lagrangian/MPPICFoam/column/constant/kinematicCloudProperties

@@ -140,6 +140,7 @@ subModels
     {
         alphaMin 0.0001;
         rhoMin 1.0;
+        applyLimiting true;
         applyGravity true;
         particleStressModel
         {

tutorials/lagrangian/MPPICFoam/cyclone/constant/kinematicCloudProperties

@@ -145,6 +145,7 @@ subModels
     {
         alphaMin 0.0001;
         rhoMin 1.0;
+        applyLimiting true;
         applyGravity false;
         particleStressModel
         {