From f12fef0f85ffb1a37bc13e1121271ec1d004c5f0 Mon Sep 17 00:00:00 2001
From: Thomas Lichtenegger <k3b1677@relise.jku.austriangrid.at>
Date: Wed, 20 Oct 2021 09:02:25 +0200
Subject: [PATCH] Some changes to drag and pressure force in the presence of
 fines.

---
 .../BeetstraDragPoly/BeetstraDragPoly.C       | 38 ++++++++++++-------
 .../BeetstraDragPoly/BeetstraDragPoly.H       |  8 ++--
 .../forceModel/gradPForce/gradPForce.C        | 20 ++++++++++
 .../forceModel/gradPForce/gradPForce.H        |  6 +++
 4 files changed, 55 insertions(+), 17 deletions(-)

diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C
index ae6e05e7..2c297636 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C
@@ -52,7 +52,8 @@ BeetstraDragPoly::BeetstraDragPoly
     dFine_(1.0),
     alphaP_(NULL),
     alphaSt_(NULL),
-    dSauter_(NULL)
+    dSauter_(NULL),
+    dSauterMix_(NULL)
 {
     // if fines are present, take mixture dSauter, otherwise normal dSauter
     if (fines_)
@@ -63,13 +64,11 @@ BeetstraDragPoly::BeetstraDragPoly
         volScalarField& alphaSt(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaSt")));
         alphaSt_ = &alphaSt;
         volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauterMix")));
-        dSauter_ = &dSauter;
-    }
-    else
-    {
-        volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauter")));
-        dSauter_ = &dSauter;
+        dSauterMix_ = &dSauter;
     }
+
+    volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauter")));
+    dSauter_ = &dSauter;
 }
 
 
@@ -84,15 +83,19 @@ scalar BeetstraDragPoly::effDiameter(double d, label cellI, label index) const
 {
     scalar dS = (*dSauter_)[cellI];
     scalar effD = d*d / dS + 0.064*d*d*d*d / (dS*dS*dS);
-    
+
     if (fines_)
     {
-        scalar fineCorr = dFine_*dFine_ / dS + 0.064*dFine_*dFine_*dFine_*dFine_ / (dS*dS*dS);
-        fineCorr *= d*d*d / (dFine_*dFine_*dFine_) * (*alphaSt_)[cellI] / (*alphaP_)[cellI];
-        effD += fineCorr;
+        scalar dSMix = (*dSauterMix_)[cellI];
+        scalar aP = (*alphaP_)[cellI];
+        scalar aSt = (*alphaSt_)[cellI];
+
+        effD = d*d*dS + 0.064*d*d*d*d / dS;
+        effD *= 1.0 / (dSMix*dSMix);
+        effD *= (aP + aSt)/aP;
     }
-    
-    if (particleCloud_.getParticleEffVolFactors()) 
+
+    if (particleCloud_.getParticleEffVolFactors())
     {
         scalar effVolFac = particleCloud_.particleEffVolFactor(index);
         effD *= effVolFac;
@@ -103,7 +106,14 @@ scalar BeetstraDragPoly::effDiameter(double d, label cellI, label index) const
 
 scalar BeetstraDragPoly::meanSauterDiameter(double d, label cellI) const
 {
-    return (*dSauter_)[cellI];
+    if (fines_)
+    {
+        return (*dSauterMix_)[cellI];
+    }
+    else
+    {
+        return (*dSauter_)[cellI];
+    }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H
index 8c91b7cf..5b357c95 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H
@@ -43,17 +43,19 @@ class BeetstraDragPoly
 protected:
 
     const bool fines_;
-    
+
     scalar dFine_;
 
     volScalarField* alphaP_;
-    
+
     volScalarField* alphaSt_;
 
     volScalarField* dSauter_;
 
+    volScalarField* dSauterMix_;
+
     scalar effDiameter(double, label, label) const;
-    
+
     scalar meanSauterDiameter(double, label) const;
 
 public:
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C
index ed8481ad..927be44e 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C
@@ -66,6 +66,9 @@ gradPForce::gradPForce
     p_(sm.mesh().lookupObject<volScalarField> (pFieldName_)),
     velocityFieldName_(propsDict_.lookup("velocityFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velocityFieldName_)),
+    alphaP_(NULL),
+    alphaSt_(NULL),
+    fines_(propsDict_.lookupOrDefault<bool>("fines",false)),
     useRho_(false),
     useU_(false),
     addedMassCoeff_(0.0)
@@ -102,6 +105,14 @@ gradPForce::gradPForce
         }
     }
 
+    if (fines_)
+    {
+        volScalarField& alphaP(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaP")));
+        alphaP_ = &alphaP;
+        volScalarField& alphaSt(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaSt")));
+        alphaSt_ = &alphaSt;
+    }
+
     if (propsDict_.found("useU"))
         useU_ = true;
 
@@ -148,6 +159,8 @@ void gradPForce::setForce() const
     vector gradP;
     scalar Vs;
     scalar rho;
+    scalar aP;
+    scalar aSt;
     vector position;
     vector force;
     label cellI;
@@ -179,6 +192,13 @@ void gradPForce::setForce() const
                 Vs = particleCloud_.particleVolume(index);
                 rho = forceSubM(0).rhoField()[cellI];
 
+                if (fines_)
+                {
+                    aP = (*alphaP_)[cellI];
+                    aSt = (*alphaSt_)[cellI];
+                    Vs *= (1+aSt/aP);
+                }
+
                 // calc particle's pressure gradient force
                 if (useRho_)
                     force = -Vs*gradP*rho*(1.0+addedMassCoeff_);
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.H b/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.H
index 68001a69..b8f10fda 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.H
@@ -71,6 +71,12 @@ private:
 
     const volVectorField& U_;
 
+    volScalarField* alphaP_;
+
+    volScalarField* alphaSt_;
+
+    bool fines_;
+
     bool useRho_;
 
     bool useU_;          // if false: substitution p=0.5*rho*U^2