Allow for additional body forces.

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);

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

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

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);

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();
 }

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_;
 }

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