diff --git a/applications/solvers/cfdemSolverRhoPimple/UEqn.H b/applications/solvers/cfdemSolverRhoPimple/UEqn.H
index 59f47b8a..ed1c117b 100644
--- a/applications/solvers/cfdemSolverRhoPimple/UEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimple/UEqn.H
@@ -1,4 +1,5 @@
 // Solve the Momentum equation
+particleCloud.otherForces(fOther);
 
 tmp<fvVectorMatrix> UEqn
 
@@ -16,14 +17,14 @@ UEqn().relax();
 fvOptions.constrain(UEqn());
 if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
 {
-    solve(UEqn() == -fvc::grad(p)+ Ksl*Us);
+    solve(UEqn() == -fvc::grad(p)+ Ksl*Us + fOther);
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
 }
 else if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -voidfraction*fvc::grad(p)+ Ksl*Us);
+    solve(UEqn() == -voidfraction*fvc::grad(p)+ Ksl*Us + fOther);
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
diff --git a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimpleChem/createFields.H b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimpleChem/createFields.H
index d6bc4f58..286b136c 100644
--- a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimpleChem/createFields.H
+++ b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimpleChem/createFields.H
@@ -128,6 +128,21 @@
         mesh,
         dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
     );
+    
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
 
     #ifndef compressibleCreatePhi_H
     #define compressibleCreatePhi_H
diff --git a/applications/solvers/cfdemSolverRhoPimple/createFields.H b/applications/solvers/cfdemSolverRhoPimple/createFields.H
index 1643d95e..1c69db2b 100644
--- a/applications/solvers/cfdemSolverRhoPimple/createFields.H
+++ b/applications/solvers/cfdemSolverRhoPimple/createFields.H
@@ -6,11 +6,11 @@ Info<< "Reading thermophysical properties\n" << endl;
     );
     psiThermo& thermo = pThermo();
     thermo.validate(args.executable(), "h", "e");
-
     volScalarField& p = thermo.p();
     volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
+    Info<< "Reading field rho\n" << endl;
     volScalarField rho
     (
         IOobject
@@ -25,7 +25,7 @@ Info<< "Reading thermophysical properties\n" << endl;
     );
 
 
-
+    
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -98,6 +98,21 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh,
         dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
     );
+    
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
 
     #ifndef compressibleCreatePhi_H
     #define compressibleCreatePhi_H
diff --git a/applications/solvers/cfdemSolverRhoPimple/pEqn.H b/applications/solvers/cfdemSolverRhoPimple/pEqn.H
index 35a3b50b..967dc5d9 100644
--- a/applications/solvers/cfdemSolverRhoPimple/pEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimple/pEqn.H
@@ -10,6 +10,7 @@ volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
 surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
+surfaceScalarField phiOther("phiOther", fvc::interpolate(rhoeps*rAU*fOther)& mesh.Sf());
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -77,6 +78,7 @@ else
           + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p)
 	  + fvc::div(phiUs)
+	  + fvc::div(phiOther)
           ==
             fvOptions(psi, p, rho.name())
         );
@@ -87,7 +89,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA + pEqn.flux()+phiUs;
+            phi = phiHbyA + pEqn.flux()+phiUs + phiOther;
         }
     }
 }
@@ -107,7 +109,7 @@ rho.relax();
 Info<< "rho max/min : " << max(rho).value()
     << " " << min(rho).value() << endl;
 
-U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);
+U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us - fOther);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
index 9f5baffe..e6307385 100644
--- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
+++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
@@ -753,7 +753,8 @@ void cfdemCloud::resetArray(double**& array,int length,int width,double resetVal
 
 void cfdemCloud::otherForces(volVectorField& forcefield)
 {
-  forcefield = vector::zero;
+  forcefield.internalField() = vector::zero;
+  forcefield.boundaryField() = vector::zero;
   for (int i=0;i<otherForceModels_.size();i++)
       forcefield += otherForceModel_[i]().exportForceField();
 }
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C
index 3fc03a59..01dae3a0 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C
@@ -456,7 +456,29 @@ void FinesFields::updateUDyn()
     volScalarField denom = FanningCoeff_ + DragCoeff_;
 
     uDyn_ = num / denom;
-
+    
+    // limit uDyn for stability reasons
+    forAll(uDyn_,cellI)
+    {
+        scalar mU(mag(U_[cellI]));
+	scalar muDyn(mag(uDyn_[cellI]));
+        if(muDyn > mU && muDyn > SMALL)
+	{
+	    uDyn_[cellI] *= mU / muDyn;
+	}
+    }
+    
+    forAll(uDyn_.boundaryField(), patchI)
+        forAll(uDyn_.boundaryField()[patchI], faceI)
+        {
+	    scalar mU(mag(U_.boundaryField()[patchI][faceI]));
+	    scalar muDyn(mag(uDyn_.boundaryField()[patchI][faceI]));
+            if(muDyn > mU && muDyn > SMALL)
+	    {
+	        uDyn_.boundaryField()[patchI][faceI] *= mU / muDyn;
+	    }
+        }
+    
     alphaDyn_.correctBoundaryConditions();
     massFluxDyn_ = rhoFine_ * alphaDyn_ * uDyn_;
 }
diff --git a/src/lagrangian/cfdemParticleComp/Make/files b/src/lagrangian/cfdemParticleComp/Make/files
index 0c21b7c8..d5eb919f 100644
--- a/src/lagrangian/cfdemParticleComp/Make/files
+++ b/src/lagrangian/cfdemParticleComp/Make/files
@@ -62,7 +62,7 @@ $(forceModels)/dSauter/dSauter.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C
 $(forceModels)/Fines/FanningDynFines.C
-$(forceModels)/Fines/FinesStatErgun.C
+$(forceModels)/Fines/ErgunStatFines.C
 
 $(forceSubModels)/forceSubModel/newForceSubModel.C
 $(forceSubModels)/forceSubModel/forceSubModel.C