ENH: Cohesion force based on energy density between particles.

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,12 +92,14 @@ Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
     alpha_(readScalar(this->coeffDict().lookup("alpha"))),
     b_(readScalar(this->coeffDict().lookup("b"))),
     mu_(readScalar(this->coeffDict().lookup("mu"))),
+    cohesionEnergyDensity_
+    (
+        readScalar(this->coeffDict().lookup("cohesionEnergyDensity"))
+    ),
+    cohesion_(false),
     collisionResolutionSteps_
     (
-        readScalar
-        (
-            this->coeffDict().lookup("collisionResolutionSteps")
-        )
+        readScalar(this->coeffDict().lookup("collisionResolutionSteps"))
     ),
     volumeFactor_(1.0),
     useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
@@ -116,6 +118,8 @@ Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
     scalar G = E/(2.0*(1.0 + nu));
 
     Gstar_ = G/(2.0*(2.0 - nu));
+
+    cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
 }
 
 
@@ -183,13 +187,15 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
         dBEff *= cbrt(pB.nParticle()*volumeFactor_);
     }
 
-    scalar normalOverlapMag = 0.5*(dAEff + dBEff) - mag(r_AB);
+    scalar r_AB_mag = mag(r_AB);
+
+    scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;
 
     if (normalOverlapMag > 0)
     {
         //Particles in collision
 
-        vector rHat_AB = r_AB/(mag(r_AB) + VSMALL);
+        vector rHat_AB = r_AB/(r_AB_mag + VSMALL);
 
         vector U_AB = pA.U() - pB.U();
 
@@ -208,6 +214,15 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
             rHat_AB
            *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
 
+        // Cohesion force
+        if (cohesion_)
+        {
+            fN_AB +=
+                -cohesionEnergyDensity_
+                *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
+                *rHat_AB;
+        }
+
         pA.f() += fN_AB;
         pB.f() += -fN_AB;
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H

@@ -34,6 +34,7 @@ Description
 
 #include "PairModel.H"
 #include "CollisionRecordList.H"
+#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,6 +68,12 @@ class PairSpringSliderDashpot
         //- Coefficient of friction in for tangential sliding
         scalar mu_;
 
+        //- Cohesion energy density [J/m^3]
+        scalar cohesionEnergyDensity_;
+
+        // Switch cohesion on and off
+        bool cohesion_;
+
         //- The number of steps over which to resolve the minimum
         //  harmonic approximation of the collision period
         scalar collisionResolutionSteps_;
@@ -129,6 +136,23 @@ public:
             return volumeFactor_;
         }
 
+        // Return the area of overlap between two spheres of radii rA and rB,
+        // centres separated by a distance rAB.  Assumes rAB < (rA + rB).
+        inline scalar overlapArea(scalar rA, scalar rB, scalar rAB) const
+        {
+            // From:
+            // http://mathworld.wolfram.com/Sphere-SphereIntersection.html
+            return
+                mathematical::pi/4.0
+               /sqr(rAB)
+               *(
+                    (-rAB + rA - rB)
+                   *(-rAB - rA + rB)
+                   *(-rAB + rA + rB)
+                   *( rAB + rA + rA)
+                );
+        }
+
         //- Whether the PairModel has a timestep limit that will
         //  require subCycling
         virtual bool controlsTimestep() const;

tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties

@@ -82,6 +82,7 @@ subModels
 
     pairCollisionCoeffs
     {
+        // Maximum possible particle diameter expected at any time
         maxInteractionDistance  0.006;
 
         writeReferredParticleCloud no;
@@ -94,6 +95,7 @@ subModels
             alpha               0.12;
             b                   1.5;
             mu                  0.52;
+            cohesionEnergyDensity 0;
             collisionResolutionSteps 12;
         };
 

tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties

@@ -91,6 +91,7 @@ subModels
 
     pairCollisionCoeffs
     {
+        // Maximum possible particle diameter expected at any time
         maxInteractionDistance  0.006;
 
         writeReferredParticleCloud no;
@@ -103,6 +104,7 @@ subModels
             alpha               0.12;
             b                   1.5;
             mu                  0.52;
+            cohesionEnergyDensity 0;
             collisionResolutionSteps 12;
         };