ENH: Adding wall cohesion forces.

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

@@ -214,7 +214,8 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
             rHat_AB
            *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
 
-        // Cohesion force
+        // Cohesion force, energy density multiplied by the area of
+        // particle-particle overlap
         if (cohesion_)
         {
             fN_AB +=

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.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
@@ -76,7 +76,8 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
     typename CloudType::parcelType& p,
     const point& site,
     const WallSiteData<vector>& data,
-    scalar pREff
+    scalar pREff,
+    bool cohesion
 ) const
 {
     // wall patch index
@@ -88,14 +89,18 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
     scalar alpha = alpha_[wPI];
     scalar b = b_[wPI];
     scalar mu = mu_[wPI];
+    scalar cohesionEnergyDensity = cohesionEnergyDensity_[wPI];
+    cohesion = cohesion && cohesion_[wPI];
 
     vector r_PW = p.position() - site;
 
     vector U_PW = p.U() - data.wallData();
 
-    scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
+    scalar r_PW_mag = mag(r_PW);
 
-    vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
+    scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
+
+    vector rHat_PW = r_PW/(r_PW_mag + VSMALL);
 
     scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
 
@@ -105,6 +110,16 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
         rHat_PW
        *(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW));
 
+    // Cohesion force, energy density multiplied by the area of wall/particle
+    // overlap
+    if (cohesion)
+    {
+        fN_PW +=
+           -cohesionEnergyDensity
+           *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
+           *rHat_PW;
+    }
+
     p.f() += fN_PW;
 
     vector USlip_PW =
@@ -168,6 +183,8 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
     alpha_(),
     b_(),
     mu_(),
+    cohesionEnergyDensity_(),
+    cohesion_(),
     patchMap_(),
     maxEstarIndex_(-1),
     collisionResolutionSteps_
@@ -212,6 +229,8 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
     alpha_.setSize(nWallPatches);
     b_.setSize(nWallPatches);
     mu_.setSize(nWallPatches);
+    cohesionEnergyDensity_.setSize(nWallPatches);
+    cohesion_.setSize(nWallPatches);
 
     scalar maxEstar = -GREAT;
 
@@ -238,6 +257,13 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
 
         mu_[wPI] = readScalar(patchCoeffDict.lookup("mu"));
 
+        cohesionEnergyDensity_[wPI] = readScalar
+        (
+            patchCoeffDict.lookup("cohesionEnergyDensity")
+        );
+
+        cohesion_[wPI] = (mag(cohesionEnergyDensity_[wPI]) > VSMALL);
+
         if (Estar_[wPI] > maxEstar)
         {
             maxEstarIndex_ = wPI;
@@ -325,20 +351,22 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
             p,
             flatSitePoints[siteI],
             flatSiteData[siteI],
-            pREff
+            pREff,
+            true
         );
     }
 
     forAll(sharpSitePoints, siteI)
     {
-        // Treating sharp sites like flat sites
+        // Treating sharp sites like flat sites, except suppress cohesion
 
         evaluateWall
         (
             p,
             sharpSitePoints[siteI],
             sharpSiteData[siteI],
-            pREff
+            pREff,
+            false
         );
     }
 }

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H

@@ -65,6 +65,12 @@ class WallLocalSpringSliderDashpot
         //- Coefficient of friction in for tangential sliding
         scalarList mu_;
 
+        //- Cohesion energy density [J/m^3]
+        scalarList cohesionEnergyDensity_;
+
+        // Switch cohesion on and off
+        boolList cohesion_;
+
         //- Mapping the patch index to the model data
         labelList patchMap_;
 
@@ -115,7 +121,8 @@ class WallLocalSpringSliderDashpot
             typename CloudType::parcelType& p,
             const point& site,
             const WallSiteData<vector>& data,
-            scalar pREff
+            scalar pREff,
+            bool cohesion
         ) const;
 
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.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
@@ -77,16 +77,19 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
     const point& site,
     const WallSiteData<vector>& data,
     scalar pREff,
-    scalar kN
+    scalar kN,
+    bool cohesion
 ) const
 {
     vector r_PW = p.position() - site;
 
     vector U_PW = p.U() - data.wallData();
 
-    scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
+    scalar r_PW_mag = mag(r_PW);
 
-    vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
+    scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
+
+    vector rHat_PW = r_PW/(r_PW_mag + VSMALL);
 
     scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
 
@@ -94,6 +97,16 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
         rHat_PW
        *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
 
+    // Cohesion force, energy density multiplied by the area of wall/particle
+    // overlap
+    if (cohesion)
+    {
+        fN_PW +=
+           -cohesionEnergyDensity_
+           *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
+           *rHat_PW;
+    }
+
     p.f() += fN_PW;
 
     vector USlip_PW =
@@ -157,6 +170,11 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
     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
@@ -183,6 +201,8 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
     Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
 
     Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
+
+    cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
 }
 
 
@@ -266,13 +286,14 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
             flatSitePoints[siteI],
             flatSiteData[siteI],
             pREff,
-            kN
+            kN,
+            cohesion_
         );
     }
 
     forAll(sharpSitePoints, siteI)
     {
-        // Treating sharp sites like flat sites
+        // Treating sharp sites like flat sites, except suppress cohesion
 
         evaluateWall
         (
@@ -280,7 +301,8 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
             sharpSitePoints[siteI],
             sharpSiteData[siteI],
             pREff,
-            kN
+            kN,
+            false
         );
     }
 }

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H

@@ -65,6 +65,12 @@ class WallSpringSliderDashpot
         //- 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_;
@@ -110,7 +116,8 @@ class WallSpringSliderDashpot
             const point& site,
             const WallSiteData<vector>& data,
             scalar pREff,
-            scalar kN
+            scalar kN,
+            bool cohesion
         ) const;
 
 

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

@@ -112,6 +112,7 @@ subModels
                 alpha           0.12;
                 b               1.5;
                 mu              0.43;
+                cohesionEnergyDensity 0;
             }
             frontAndBack
             {
@@ -120,6 +121,7 @@ subModels
                 alpha           0.12;
                 b               1.5;
                 mu              0.1;
+                cohesionEnergyDensity 0;
             }
         };
     }

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

@@ -121,6 +121,7 @@ subModels
                 alpha           0.12;
                 b               1.5;
                 mu              0.43;
+                cohesionEnergyDensity 0;
             }
             frontAndBack
             {
@@ -129,6 +130,7 @@ subModels
                 alpha           0.12;
                 b               1.5;
                 mu              0.1;
+                cohesionEnergyDensity 0;
             }
         };
     }