Corrected the two-way coupling of the immersed boundary method.

Corrected the torque calculation for the cfdemSolverIBTorque and
modified the ShirgaonkarIBTorque calculations.

Included fvOptions.H in the solvers for activation the usage of
fvOptions.

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -53,6 +53,8 @@ Contributions
 
 #include "cellSet.H"
 
+#include "fvOptions.H"  // added the fvOptions library
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])

applications/solvers/cfdemSolverIBTorque/cfdemSolverIBTorque.C

@@ -53,6 +53,8 @@ Contributions
 
 #include "cellSet.H"
 
+#include "fvOptions.H"      // added fvOptions
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])

etc/solver-list.txt

@@ -3,3 +3,4 @@ cfdemSolverPiso/dir
 cfdemSolverRhoPimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
+cfdemSolverIBTorque/dir

src/lagrangian/cfdemParticle/Make/files

@@ -25,6 +25,7 @@ $(cfdemCloud)/cfdemCloudIO.C
 derived/cfdemCloudIB/cfdemCloudIB.C
 derived/cfdemCloudMS/cfdemCloudMS.C
 derived/cfdemCloudEnergy/cfdemCloudEnergy.C
+derived/cfdemCloudIBmodified/cfdemCloudIBmodified.C
 
 $(cfdTools)/global.C
 $(cfdTools)/newGlobal.C
@@ -69,6 +70,7 @@ $(forceModels)/Fines/FinesFields.C
 $(forceModels)/Fines/FanningDynFines.C
 $(forceModels)/Fines/ErgunStatFines.C
 $(forceModels)/granKineticEnergy/granKineticEnergy.C
+$(forceModels)/ShirgaonkarIBTorque/ShirgaonkarIBTorque.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C

src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C

@@ -192,6 +192,7 @@ void cfdemCloudIB::calcVelocityCorrection
                     // impose field velocity
                     U[cellI]=(1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI];
                 }
+
             }
         //}
     }

src/lagrangian/cfdemParticle/derived/cfdemCloudIBmodified/cfdemCloudIBmodified.H

@@ -62,6 +62,7 @@ protected:
     label pRefCell_;
     scalar pRefValue_;
 
+    // Modifications
     mutable double **DEMTorques_;
 
     mutable bool haveEvolvedOnce_;
@@ -89,6 +90,7 @@ public:
 
         bool reAllocArrays();
 
+    // Modification
         inline double ** DEMTorques() const;
 
         bool evolve();
@@ -102,7 +104,6 @@ public:
         {
             return angularVelocities_;
         }
-
 };
 
 

src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C

@@ -65,6 +65,7 @@ ShirgaonkarIB::ShirgaonkarIB
     propsDict_(dict.subDict(typeName + "Props")),
     verbose_(false),
     twoDimensional_(false),
+    multisphere_(false), //  drag for a multisphere particle
     depth_(1),
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
@@ -85,6 +86,9 @@ ShirgaonkarIB::ShirgaonkarIB
         Info << "depth of domain is assumed to be :" << depth_ << endl;
     }
 
+    // Switch to initiate the multisphere calculation
+    if(propsDict_.found("multisphere")) multisphere_=true;
+
     // init force sub model
     setForceSubModels(propsDict_);
 
@@ -151,8 +155,58 @@ void ShirgaonkarIB::setForce() const
             if(verbose_) Info << "impForces = " << impForces()[index][0]<<","<<impForces()[index][1]<<","<<impForces()[index][2] << endl;
         //}
     }
+
+    // Calculate the force if the particle is multisphere template
+    if(multisphere_) calcForce();
 }
 
+void ShirgaonkarIB::calcForce() const
+{
+    label cellID;
+    vector dragMS;
+
+    volVectorField h=forceSubM(0).IBDragPerV(U_,p_);
+
+    dragMS = vector::zero;
+
+    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
+    {
+        int prev;
+
+        for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
+        {
+            cellID = particleCloud_.cellIDs()[index][subCell];
+
+            prev = 0;
+
+            if (cellID > -1 && index == 0) // Force on 1st particle
+            {
+                dragMS += h[cellID]*h.mesh().V()[cellID];
+
+            }
+
+            else if(cellID > -1 && index > 0)
+            {
+
+                // Check for cellID in previous particles
+                for(int i=index-1; i>= 0; i--)
+                {
+                    for(int j=0; j< particleCloud_.voidFractionM().cellsPerParticle()[i][0];j++)
+                    {
+                        label cell = particleCloud_.cellIDs()[i][j];
+                        if(cellID == cell){prev++;}
+                    }
+                }
+                if(prev == 0){ dragMS += h[cellID]*h.mesh().V()[cellID];}
+            }
+        }
+    }
+
+    Info << "Drag force on particle clump = " << dragMS[0]<<","<<dragMS[1]<<","<<dragMS[2] << endl;
+}
+
+
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.H

@@ -66,6 +66,8 @@ private:
 
     bool depth_;
 
+    bool multisphere_;
+
     word velFieldName_;
 
     const volVectorField& U_;
@@ -96,6 +98,7 @@ public:
 
     // Member Functions
         void setForce() const;
+        void calcForce() const;
 };
 
 

src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIBTorque/ShirgaonkarIBTorque.C

@@ -157,7 +157,7 @@ void ShirgaonkarIBTorque::setForce() const
             if(verbose_) Info << "impForces = " << impForces()[index][0]<<","<<impForces()[index][1]<<","<<impForces()[index][2] << endl;
             if(useTorque_)
             {
-                    for(int j=0;j<3;j++) particleCloud_.DEMTorques()[index][j] = torque[j]; // or is it particleCloud_.DEMTorques()[index][j] += torque[j];
+                    for(int j=0;j<3;j++) particleCloud_.DEMTorques()[index][j] = torque[j]; // Adding to the particle torque;
                     Info << "DEMTorques = " << particleCloud_.DEMTorques()[index][0] << "," << particleCloud_.DEMTorques()[index][1] << "," << particleCloud_.DEMTorques()[index][2] << endl;
             }
     }

src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C

@@ -169,7 +169,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
                             scalar lambda = segmentParticleIntersection(index, minPeriodicParticlePos, vertexPosition, cellCentrePosition);
                             voidfractionNext_[particleCenterCellID] -= ratio * lambda;
                         }
-                    }
+                     }
                 } //end particle partially overlapping with cell
 
                 //generating list with cell and subcells