diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
index 5cb2092818c443bc1a63a302f9ead652a74a5dd4..ab2b45799150685d36f3b48334a96bae5717e82a 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
+++ b/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
- phiGByA
- (
- "phiGByA",
- deltaT*(g & mesh.Sf())*fvc::interpolate(1.0 - rhoc/rho)
- );
-
- alphaEqn += fvm::div(phiGByA, alpha_);
+ 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>
(
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.H b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.H
index 9f90dabcbf58d1d42235614549bae875ef3c08cd..dc3117000400b5e616d38ae1699404e162bf136c 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.H
+++ b/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_;
diff --git a/tutorials/lagrangian/MPPICFoam/Goldschmidt/constant/kinematicCloudProperties b/tutorials/lagrangian/MPPICFoam/Goldschmidt/constant/kinematicCloudProperties
index fdb104bc2d338ac0908c4493ffc1c35c8666e0f6..9f46325f606942dd6216d6d8198688aefc568fc0 100644
--- a/tutorials/lagrangian/MPPICFoam/Goldschmidt/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/MPPICFoam/Goldschmidt/constant/kinematicCloudProperties
@@ -148,6 +148,7 @@ subModels
{
alphaMin 0.0001;
rhoMin 1.0;
+ applyLimiting true;
applyGravity false;
particleStressModel
{
diff --git a/tutorials/lagrangian/MPPICFoam/column/constant/kinematicCloudProperties b/tutorials/lagrangian/MPPICFoam/column/constant/kinematicCloudProperties
index f086e389ec1305cffae1301d6cb537f26bee6704..00d205137b5a6fd8432bdb6e7cc888f11dd06be4 100644
--- a/tutorials/lagrangian/MPPICFoam/column/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/MPPICFoam/column/constant/kinematicCloudProperties
@@ -140,6 +140,7 @@ subModels
{
alphaMin 0.0001;
rhoMin 1.0;
+ applyLimiting true;
applyGravity true;
particleStressModel
{
diff --git a/tutorials/lagrangian/MPPICFoam/cyclone/constant/kinematicCloudProperties b/tutorials/lagrangian/MPPICFoam/cyclone/constant/kinematicCloudProperties
index 78bac15385f98a382ee98194643ad9cfd49c8d9b..d72fd6316304611fc0b37626202f448fd97adad7 100644
--- a/tutorials/lagrangian/MPPICFoam/cyclone/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/MPPICFoam/cyclone/constant/kinematicCloudProperties
@@ -145,6 +145,7 @@ subModels
{
alphaMin 0.0001;
rhoMin 1.0;
+ applyLimiting true;
applyGravity false;
particleStressModel
{