Beetstra drag mono can take multiple cg factors, added Beetstra drag poly for polydisp systems, including polydisp parcels and effect of fines phase.

src/lagrangian/cfdemParticle/Make/files

@@ -64,6 +64,7 @@ $(forceModels)/particleCellVolume/particleCellVolume.C
 $(forceModels)/fieldTimeAverage/fieldTimeAverage.C
 $(forceModels)/volWeightedAverage/volWeightedAverage.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C

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

@@ -49,14 +49,26 @@ BeetstraDrag::BeetstraDrag
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    minVoidfraction_(propsDict_.lookupOrDefault<scalar>("minVoidfraction",0.1)),
     UsFieldName_(propsDict_.lookup("granVelFieldName")),
     UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
     scaleDia_(1.),
-    scaleDrag_(1.)
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
+    scaleDrag_(1.),
+    rhoP_(0.),
+    rho_(0.),
+    Lc2_(0.),
+    dPrim_(0.),
+    nuf_(0.),
+    g_(9.81),
+    k_(0.05),
+    useGC_(false),
+    usePC_(false)
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
@@ -78,10 +90,42 @@ BeetstraDrag::BeetstraDrag
     forceSubM(0).readSwitches();
 
     particleCloud_.checkCG(true);
-    if (propsDict_.found("scale"))
+    if (propsDict_.found("scale") && typeCG_.size()==1)
+    {
+        // if "scale" is specified and there's only one single type, use "scale"
         scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        typeCG_[0] = scaleDia_;
+    }
     if (propsDict_.found("scaleDrag"))
+    {
         scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
+    }
+
+    if (typeCG_.size()>1) multiTypes_ = true;
+
+    if (propsDict_.found("useFilteredDragModel"))
+    {
+        useGC_ = true;
+        g_=propsDict_.lookupOrDefault<scalar>("g", 9.81);
+        dPrim_=scalar(readScalar(propsDict_.lookup("dPrim")));
+        rhoP_=scalar(readScalar(propsDict_.lookup("rhoP")));
+        rho_=scalar(readScalar(propsDict_.lookup("rho")));
+        nuf_=scalar(readScalar(propsDict_.lookup("nuf")));
+	    scalar ut = terminalVelocity(1., dPrim_, nuf_, rho_, rhoP_, g_);
+	    scalar Frp = ut*ut/g_/dPrim_;
+        Lc2_ = ut*ut/g_*pow(Frp, -.6666667); // n is hardcoded as -2/3
+        Info << "using grid coarsening correction with Lc2 = " << Lc2_ << " and ut = " << ut << " and Frp = " << Frp<< endl;
+
+        if (propsDict_.found("useParcelSizeDependentFilteredDrag"))
+        {
+            usePC_ = true;
+            if (propsDict_.found("k"))
+                k_=scalar(readScalar(propsDict_.lookup("k")));
+            Info << "using particle coarsening correction with k = " << k_ << endl;
+        }
+    }
+
+
 
 }
 
@@ -96,9 +140,9 @@ BeetstraDrag::~BeetstraDrag()
 
 void BeetstraDrag::setForce() const
 {
-    if (scaleDia_ > 1)
+    if (typeCG_.size()>1 || typeCG_[0] > 1)
     {
-        Info << "Beetstra using scale = " << scaleDia_ << endl;
+        Info << "Beetstra using scale = " << typeCG_ << endl;
     }
     else if (particleCloud_.cg() > 1)
     {
@@ -119,12 +163,17 @@ void BeetstraDrag::setForce() const
     vector Ur(0,0,0);
     scalar ds(0);
     scalar ds_scaled(0);
-    scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_;
+    scalar scaleDia3 = typeCG_[0]*typeCG_[0]*typeCG_[0];
     scalar nuf(0);
     scalar rho(0);
     scalar magUr(0);
     scalar Rep(0);
     scalar localPhiP(0);
+    scalar GCcorr(1.);
+    scalar PCcorr(1.);
+
+    scalar cg = typeCG_[0];
+    label partType = 1;
 
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
@@ -146,27 +195,37 @@ void BeetstraDrag::setForce() const
             if (cellI > -1) // particle found
             {
 
-                if( forceSubM(0).interpolation() )
+                if ( forceSubM(0).interpolation() )
                 {
                     position     = particleCloud_.position(index);
                     voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
                     Ufluid       = UInterpolator_.interpolate(position,cellI);
                     //Ensure interpolated void fraction to be meaningful
                     // Info << " --> voidfraction: " << voidfraction << endl;
-                    if(voidfraction>1.00) voidfraction = 1.0;
+                    if (voidfraction > 1.00) voidfraction = 1.0;
-                    if(voidfraction<0.10) voidfraction = 0.10;
+                    if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
                 }
                 else
                 {
                     voidfraction = voidfraction_[cellI];
                     Ufluid = U_[cellI];
                 }
+                // in case a fines phase is present, void fraction needs to be adapted
+                adaptVoidfraction(voidfraction, cellI);
+
+                if (multiTypes_) 
+                {
+                    partType = particleCloud_.particleType(index);
+                    cg = typeCG_[partType - 1];
+                    scaleDia3 = cg*cg*cg;
+                }
 
                 Us = particleCloud_.velocity(index);
                 Ur = Ufluid-Us;
                 magUr = mag(Ur);
                 ds = 2*particleCloud_.radius(index);
-                ds_scaled = ds/scaleDia_;
+                ds_scaled = ds/cg;
+
                 rho = rhoField[cellI];
                 nuf = nufField[cellI];
 
@@ -175,14 +234,31 @@ void BeetstraDrag::setForce() const
 
                 // calc particle's drag coefficient (i.e., Force per unit slip velocity and Stokes drag)
 
-                Rep=ds_scaled*voidfraction*magUr/nuf+SMALL;
+                Rep=ds_scaled*voidfraction*magUr/nuf + SMALL;
-                dragCoefficient = 10.0*localPhiP/(voidfraction*voidfraction) +
+
-                                  voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP)) +
+                dragCoefficient = F(Rep, voidfraction)
-                                  0.413*Rep/(24*voidfraction*voidfraction)*(1.0/voidfraction+3*voidfraction*localPhiP+8.4*Foam::pow(Rep,-0.343))/
+                                   *3*M_PI*nuf*rho*voidfraction
-                                  (1+Foam::pow(10,3*localPhiP)*Foam::pow(Rep,-0.5*(1+4*localPhiP)));
+                                   *effDiameter(ds_scaled, voidfraction, cellI, index)
+                                   *scaleDia3*scaleDrag_;
+                
+                // calculate filtering corrections
+                if (useGC_)
+                {
+                    GCcorr = 1.-h(localPhiP) 
+                                /(  a(localPhiP) 
+                                    *Lc2_ 
+                                    /Foam::pow(U_.mesh().V()[cellI],.33333333) 
+                                  + 1.
+                                 );
+                    if (usePC_)
+                    {
+                        PCcorr = Foam::exp(k_*(1.-cg));
+                    }
+                }
+
+                // apply filtering corrections
+                dragCoefficient *= GCcorr*PCcorr;
 
-                // calc particle's drag
-                dragCoefficient *= 3*M_PI*ds_scaled*nuf*rho*voidfraction*scaleDia3*scaleDrag_;
                 if (modelType_=="B")
                     dragCoefficient /= voidfraction;
 
@@ -198,11 +274,13 @@ void BeetstraDrag::setForce() const
                     Pout << "Us = " << Us << endl;
                     Pout << "Ur = " << Ur << endl;
                     Pout << "ds = " << ds << endl;
-                    Pout << "ds/scale = " << ds/scaleDia_ << endl;
+                    Pout << "ds/scale = " << ds/cg << endl;
                     Pout << "rho = " << rho << endl;
                     Pout << "nuf = " << nuf << endl;
                     Pout << "voidfraction = " << voidfraction << endl;
                     Pout << "Rep = " << Rep << endl;
+                    Pout << "GCcorr = " << GCcorr << endl;
+		            Pout << "PCcorr = " << PCcorr << endl;
                     Pout << "drag = " << drag << endl;
                 }
 
@@ -223,7 +301,145 @@ void BeetstraDrag::setForce() const
     }
 }
 
+/*********************************************************
+ * "Drag Force of Intermediate Reynolds Number Flow Past *
+ *  Mono- and Bidisperse Arrays of Spheres", eq. 16      *
+ *  R Beetstra, M. A. van der Hoef, JAM Kuipers          *
+ *  AIChE Journal 53(2) (2007)                           *
+ *********************************************************/
+double BeetstraDrag::F(double voidfraction, double Rep) const
+{
+    double localPhiP = max(SMALL,min(1.-SMALL,1.-voidfraction));
+    return   10.0*localPhiP/(voidfraction*voidfraction) 
+           + voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP)) 
+           + 0.413*Rep/(24*voidfraction*voidfraction)
+                       *(1.0/voidfraction
+                          +3*voidfraction*localPhiP
+                          +8.4*Foam::pow(Rep,-0.343)
+                        )
+                        /(1+Foam::pow(10,3*localPhiP)
+                            *Foam::pow(Rep,-0.5*(1+4*localPhiP))
+                         );
 
+}
+ 
+
+/*********************************************************
+ * "A drag model for filtered Euler-Lagange simulations  *
+ *  of clustered gas-particle suspension",               *
+ *  S. Radl, S. Sundaresan,                              *
+ *  Chemical Engineering Science 117 (2014).             *
+ *********************************************************/
+double BeetstraDrag::a(double phiP) const
+{
+    double a0m = 0.;
+    double a1m = 0.;
+    double a2m = 0.;
+    double a3m = 0.;
+    double phipam = 0.;
+    if      (phiP < 0.016)
+    {
+      a0m = 21.51;
+    }
+    else if (phiP < 0.100)
+    {
+      a0m =  1.96; a1m =  29.40; a2m =  164.91; a3m = -1923.;
+    }
+    else if (phiP < 0.180)
+    {
+      a0m =  4.63; a1m =   4.68; a2m = -412.04; a3m =  2254.; phipam = 0.10;
+    }
+    else if (phiP < 0.250)
+    {
+      a0m =  3.52; a1m = -17.99; a2m =  128.80; a3m =  -603.; phipam = 0.18;
+    }
+    else if (phiP < 0.400)
+    {
+      a0m =  2.68; a1m =  -8.20; a2m =    2.18; a3m = 112.33; phipam = 0.25;
+    }
+    else
+    {
+      a0m =  1.79;
+    }
+    return a0m + a1m*(phiP-phipam) + a2m*pow(phiP-phipam,2.) + a3m*pow(phiP-phipam,3.);
+}
+
+double BeetstraDrag::h(double phiP) const
+{
+  double h0m = 0.;
+  double h1m = 0.;
+  double h2m = 0.;
+  double h3m = 0.;
+  double phiphm = 0.;
+
+  if      (phiP < 0.03)
+  {
+                 h1m = 7.97;
+  }
+  else if (phiP < 0.08)
+  {
+    h0m = 0.239; h1m =  4.640; h2m =   -4.410; h3m =  253.630; phiphm = 0.03;
+  }
+  else if (phiP < 0.12)
+  {
+    h0m = 0.492; h1m =  6.100; h2m =   33.630; h3m = -789.600; phiphm = 0.08;
+  }
+  else if (phiP < 0.18)
+  {
+    h0m = 0.739; h1m =  5.010; h2m =  -61.100; h3m =  310.800; phiphm = 0.12;
+  }
+  else if (phiP < 0.34)
+  {
+    h0m = 0.887; h1m =  1.030; h2m =   -5.170; h3m =    5.990; phiphm = 0.18;
+  }
+  else if (phiP < 0.48)
+  {
+    h0m = 0.943; h1m = -0.170; h2m =   -2.290; h3m =   -9.120; phiphm = 0.34;
+  }
+  else if (phiP < 0.55)
+  {
+    h0m = 0.850; h1m = -1.350; h2m =   -6.130; h3m = -132.600; phiphm = 0.48;
+  }
+  else
+  {
+    h0m = 0.680; h1m = -2.340; h2m = -225.200;                 phiphm = 0.55;
+  }
+  return h0m + h1m*(phiP-phiphm) + h2m*pow(phiP-phiphm,2) + h3m*pow(phiP-phiphm,3);
+}
+
+double BeetstraDrag::terminalVelocity(double voidfraction, double dp, double nuf, double rhof, double rhop, double g) const
+{
+    scalar u0(dp*dp*fabs(rhof-rhop)*g/18./rhof/nuf);
+    scalar Re(u0*dp/nuf);
+    scalar res(1.);
+    scalar u(u0);
+    scalar Fi(0);
+    scalar CdSt(0);
+    Info << "uo: " << u0<<endl;
+    int i = 0;
+    while ((res > 1.e-6) && (i<100))
+    {
+        Info << "Iteration " << i;
+        u0 = u;
+        Info << ", u0 = " << u0;
+        CdSt = 24/Re;
+        Info << ", CdSt = " << CdSt;
+        Fi = F(voidfraction, Re);
+        Info << ", F = ";
+        u = sqrt(1.333333333*fabs(rhof-rhop)*g*dp
+                  /(CdSt*voidfraction*Fi*rhof)
+            );
+        Info << ", u = " << u;
+        Re = fabs(u)*dp/nuf*voidfraction;
+        res = fabs((u-u0)/u);
+        Info << "Res: " << res << endl;
+        i++;
+    }
+    if (res >1.e-6)
+        FatalError << "Terminal velocity calculation diverged!" << endl;
+
+    return u;
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -40,9 +40,11 @@ class BeetstraDrag
 :
     public forceModel
 {
-private:
+protected:
     dictionary propsDict_;
 
+    bool multiTypes_;
+
     word velFieldName_;
 
     const volVectorField& U_;
@@ -51,14 +53,50 @@ private:
 
     const volScalarField& voidfraction_;
 
+    const scalar minVoidfraction_;
+
     word UsFieldName_;
 
     const volVectorField& UsField_;
 
     mutable scalar scaleDia_;
 
+    scalarList typeCG_;
+
     mutable scalar scaleDrag_;
 
+    mutable scalar rhoP_;
+
+    mutable scalar rho_;
+
+    mutable scalar Lc2_;
+
+    mutable scalar dPrim_;
+
+    mutable scalar nuf_;
+
+    mutable scalar g_;
+
+    mutable scalar k_;
+
+    bool useGC_;
+
+    bool usePC_;
+
+    virtual void adaptVoidfraction(double&, label) const {}
+
+    virtual scalar effDiameter(double d, double voidfraction, label cellI, label index) const {return d;}
+
+    double F(double, double) const;
+
+    double terminalVelocity(double, double, double, double, double, double) const; 
+
+    double a(double) const;
+
+    double h(double) const;
+
+
+
 public:
 
     //- Runtime type information

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

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "BeetstraDragPoly.H"
+#include "addToRunTimeSelectionTable.H"
+#include "averagingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(BeetstraDragPoly, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    BeetstraDragPoly,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+BeetstraDragPoly::BeetstraDragPoly
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    BeetstraDrag(dict,sm),
+    fines_(propsDict_.lookupOrDefault<bool>("fines",false)),
+    dSauter_(sm.mesh().lookupObject<volScalarField> ("dSauter"))
+{
+    if (fines_)
+    {
+        volScalarField& alphaSt(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaSt")));
+        alphaSt_.set(&alphaSt);
+        volScalarField& dSauterMix(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauterMix")));
+        dSauterMix_.set(&dSauterMix);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+BeetstraDragPoly::~BeetstraDragPoly()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void BeetstraDragPoly::adaptVoidfraction(double& voidfraction, label cellI) const
+{
+    voidfraction -= alphaSt_()[cellI];
+    if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
+}
+
+scalar BeetstraDragPoly::effDiameter(double d, double voidfraction, label cellI, label index) const
+{
+    scalar dS = dSauter_[cellI];
+    scalar effD = d*d/dS;
+    if (fines_)
+    {
+        scalar pureVoidfraction = voidfraction_[cellI];
+        scalar dSmix = dSauterMix_()[cellI];
+        effD *= pureVoidfraction / voidfraction * (1 - voidfraction) / (1 - pureVoidfraction);
+        effD *= dS * dS / (dSmix * dSmix);
+    }
+    if (particleCloud_.getParticleEffVolFactors()) 
+    {
+        scalar effVolFac = particleCloud_.particleEffVolFactor(index);
+        effD *= effVolFac;
+    }
+    return effD;
+    
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+Description
+    drag law for monodisperse systems according to
+    Beetstra et al. AIChE J 53.2 (2007)
+
+SourceFiles
+    BeetstraDragPoly.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef BeetstraDragPoly_H
+#define BeetstraDragPoly_H
+
+#include "BeetstraDrag.H"
+#include "interpolationCellPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class BeetstraDragPoly Declaration
+\*---------------------------------------------------------------------------*/
+
+class BeetstraDragPoly
+:
+    public BeetstraDrag
+{
+protected:
+
+    const bool fines_;
+
+    autoPtr<volScalarField> alphaSt_;
+
+    const volScalarField& dSauter_;
+
+    autoPtr<volScalarField> dSauterMix_;
+
+    void adaptVoidfraction(double&, label) const;
+
+    scalar effDiameter(double, double, label, label) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("BeetstraDragPoly");
+
+
+    // Constructors
+
+        //- Construct from components
+        BeetstraDragPoly
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~BeetstraDragPoly();
+
+
+    // Member Functions
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -56,7 +56,7 @@ dSauter::dSauter
     multiTypes_(false),
     d2_(NULL),
     d3_(NULL),
-    typeCG_(propsDict_.lookupOrDefault<labelList>("coarseGrainingFactors",labelList(1,1))),
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     d2Field_
     (   IOobject
         (
@@ -127,7 +127,7 @@ void dSauter::allocateMyArrays() const
 
 void dSauter::setForce() const
 {
-    if (typeCG_.size()>1 || typeCG_[0] > 1)
+    if (typeCG_.size()>1 || typeCG_[0] > 1.0)
     {
         Info << "dSauter using CG factor(s) = " << typeCG_ << endl;
     }
@@ -135,8 +135,8 @@ void dSauter::setForce() const
     allocateMyArrays();
 
     label cellI = 0;
-    label cg = 1;
     label partType = 1;
+    scalar cg = typeCG_[0];
     scalar ds = 0.0;
     scalar effVolFac = 1.0;
 
@@ -150,8 +150,11 @@ void dSauter::setForce() const
             {
                 effVolFac = particleCloud_.particleEffVolFactor(index);
             }
-            if (multiTypes_) partType = particleCloud_.particleType(index);
+            if (multiTypes_) 
+            {
+                partType = particleCloud_.particleType(index);
                 cg = typeCG_[partType - 1];
+            }
             ds = particleCloud_.d(index);
             d2_[index][0] = ds*ds*effVolFac*cg;
             d3_[index][0] = ds*ds*ds*effVolFac;

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

@@ -49,7 +49,7 @@ private:
     
     mutable double **d3_;
 
-    labelList typeCG_;
+    scalarList typeCG_;
     
     mutable volScalarField d2Field_;
     

src/lagrangian/cfdemParticleComp/Make/files

@@ -60,6 +60,7 @@ $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C
 $(forceModels)/particleCellVolume/particleCellVolume.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties

@@ -33,7 +33,7 @@ couplingInterval 100;
 
 voidFractionModel divided;//centre;//
 
-locateModel engine;//turboEngine;//
+locateModel engine;//turboEngineM2M;//
 
 meshMotionModel noMeshMotion;
 
@@ -49,7 +49,7 @@ averagingModel dense;//dilute;//
 
 clockModel off;//standardClock;//
 
-smoothingModel off;// constDiffSmoothing; //
+smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; //
 
 forceModels
 (
@@ -61,6 +61,7 @@ forceModels
     GidaspowDrag
     //Archimedes
     //volWeightedAverage
+    //totalMomentumExchange
     particleCellVolume
 );
 
@@ -74,6 +75,14 @@ turbulenceModelType "turbulenceProperties";
 //===========================================================================//
 // sub-model properties
 
+localPSizeDiffSmoothingProps
+{
+    lowerLimit 0.1;
+    upperLimit 1e10;
+    dSmoothingLength 1.5e-3;
+    Csmoothing    1.0; 
+}
+
 constDiffSmoothingProps
 {
     lowerLimit 0.1;
@@ -123,7 +132,13 @@ volWeightedAverageProps
     lowerThreshold 0;
     verbose true;
 }
-
+totalMomentumExchangeProps
+{
+    implicitMomExFieldName "Ksl";
+    explicitMomExFieldName "none";
+    fluidVelFieldName "U";
+    granVelFieldName "Us";
+}
 GidaspowDragProps
 {
     verbose true;
@@ -147,12 +162,15 @@ KochHillDragProps
 BeetstraDragProps
 {
     velFieldName "U";
-    gravityFieldName "g";
-    rhoParticle     2000.;
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
     interpolation true;
-    useFilteredDragModel ;
+//    useFilteredDragModel;
-    useParcelSizeDependentFilteredDrag ;
+//    useParcelSizeDependentFilteredDrag;
+    g       9.81;
+    rhoP    7000.;
+    rho     10.;
+    nuf     1.5e-4;
     k       0.05;
     aLimit  0.0;
 //    verbose true;