Renaming SpringSliderDashpot to PairSpringSliderDashpot to allow a

WallSpringSliderDashpot model.

Correctly identifying non-flat - sharp - wall interactions, i.e. with
concave points and edges.  Rigoruous methods for exclusion distances
of existing wall interaction sites.

Separated off wall force calculation into WallSpringSliderDashpot
model.

src/lagrangian/intermediate/parcels/include/makeParcelCollisionModels.H

@@ -34,7 +34,9 @@ License
 #include "NoCollision.H"
 #include "PairCollision.H"
 
-#include "SpringSliderDashpot.H"
+#include "PairSpringSliderDashpot.H"
+
+#include "WallSpringSliderDashpot.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,7 +62,16 @@ License
                                                                               \
     makePairModelType                                                         \
     (                                                                         \
-        SpringSliderDashpot,                                                  \
+        PairSpringSliderDashpot,                                              \
+        InteractingKinematicCloud,                                            \
+        ParcelType                                                            \
+    );                                                                        \
+                                                                              \
+    makeWallModel(InteractingKinematicCloud<ParcelType>);                     \
+                                                                              \
+    makeWallModelType                                                         \
+    (                                                                         \
+        WallSpringSliderDashpot,                                              \
         InteractingKinematicCloud,                                            \
         ParcelType                                                            \
     );

src/lagrangian/intermediate/parcels/include/makeReactingParcelCollisionModels.H

@@ -35,7 +35,9 @@ License
 #include "NoCollision.H"
 #include "PairCollision.H"
 
-#include "SpringSliderDashpot.H"
+#include "PairSpringSliderDashpot.H"
+
+#include "WallSpringSliderDashpot.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -84,7 +86,17 @@ License
                                                                               \
     makePairModelThermoType                                                   \
     (                                                                         \
-        SpringSliderDashpot,                                                  \
+        PairSpringSliderDashpot,                                              \
+        InteractingKinematicCloud,                                            \
+        ParcelType,                                                           \
+        ThermoType                                                            \
+    );                                                                        \
+                                                                              \
+    makeWallModel(InteractingKinematicCloud<ParcelType<ThermoType> >);        \
+                                                                              \
+    makeWallModelThermoType                                                   \
+    (                                                                         \
+        WallSpringSliderDashpot,                                              \
         InteractingKinematicCloud,                                            \
         ParcelType,                                                           \
         ThermoType                                                            \

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

@@ -26,6 +26,17 @@ License
 
 #include "PairCollision.H"
 #include "PairModel.H"
+#include "WallModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::PairCollision<CloudType>::cosPhiMinFlatWall = 1 - 1e-6;
+
+template<class CloudType>
+Foam::scalar Foam::PairCollision<CloudType>::flatWallDuplicateExclusion =
+    sqrt(3e-6);
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -144,11 +155,15 @@ void Foam::PairCollision<CloudType>::wallInteraction()
     const DirectInteractionList<typename CloudType::parcelType>& dil =
         il_.dil();
 
-    ReferredCellList<typename CloudType::parcelType>& ril = il_.ril();
+    const ReferredCellList<typename CloudType::parcelType>& ril = il_.ril();
 
-    DynamicList<point> allWallInteractionSites;
+    // Storage for the wall interaction sites
-    DynamicList<point> flatWallInteractionSites;
+    DynamicList<point> flatSites;
-    DynamicList<point> sharpWallInteractionSites;
+    DynamicList<scalar> flatSiteExclusionDistancesSqr;
+    DynamicList<point> otherSites;
+    DynamicList<scalar> otherSiteDistances;
+    DynamicList<point> sharpSites;
+    DynamicList<scalar> sharpSiteExclusionDistancesSqr;
 
     forAll(dil, realCellI)
     {
@@ -162,14 +177,17 @@ void Foam::PairCollision<CloudType>::wallInteraction()
         // Loop over all Parcels in cell
         forAll(cellOccupancy_[realCellI], cellParticleI)
         {
-            allWallInteractionSites.clear();
+            flatSites.clear();
-            flatWallInteractionSites.clear();
+            flatSiteExclusionDistancesSqr.clear();
-            sharpWallInteractionSites.clear();
+            otherSites.clear();
+            otherSiteDistances.clear();
+            sharpSites.clear();
+            sharpSiteExclusionDistancesSqr.clear();
 
             typename CloudType::parcelType& p =
                 *cellOccupancy_[realCellI][cellParticleI];
 
-            const point& pt = p.position();
+            const point& pos = p.position();
 
             // real wallFace interactions
 
@@ -179,7 +197,7 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                 pointHit nearest = mesh.faces()[realFaceI].nearestPoint
                 (
-                    pt,
+                    pos,
                     mesh.points()
                 );
 
@@ -189,15 +207,41 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                     normal /= mag(normal);
 
-                    vector pW = nearest.rawPoint() - pt;
+                    const vector& nearPt = nearest.rawPoint();
+
+                    vector pW = nearPt - pos;
 
                     scalar normalAlignment = normal & pW/mag(pW);
 
-                    allWallInteractionSites.append(nearest.rawPoint());
+                    if (normalAlignment > cosPhiMinFlatWall)
-
-                    if (normalAlignment > 1 - SMALL)
                     {
-                        flatWallInteractionSites.append(nearest.rawPoint());
+                        // Guard against a flat interaction being
+                        // present on the boundary of two or more
+                        // faces, which would create duplicate contact
+                        // points. Duplicates are discarded.
+                        if
+                        (
+                            !duplicatePointInList
+                            (
+                                flatSites,
+                                nearPt,
+                                sqr(p.r()*flatWallDuplicateExclusion)
+                            )
+                        )
+                        {
+                            flatSites.append(nearPt);
+
+                            flatSiteExclusionDistancesSqr.append
+                            (
+                                sqr(p.r()) - sqr(nearest.distance())
+                            );
+                        }
+                    }
+                    else
+                    {
+                        otherSites.append(nearPt);
+
+                        otherSiteDistances.append(nearest.distance());
                     }
                 }
             }
@@ -206,7 +250,7 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
             forAll(referredCellsInRange, refCellInRangeI)
             {
-                ReferredCell<typename CloudType::parcelType>& refCell =
+                const ReferredCell<typename CloudType::parcelType>& refCell =
                     ril[referredCellsInRange[refCellInRangeI]];
 
                 const labelList& refWallFaces = refCell.wallFaces();
@@ -217,7 +261,7 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                     pointHit nearest = refCell.faces()[refFaceI].nearestPoint
                     (
-                        pt,
+                        pos,
                         refCell.points()
                     );
 
@@ -227,53 +271,141 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                         normal /= mag(normal);
 
-                        vector pW = nearest.rawPoint() - pt;
+                        const vector& nearPt = nearest.rawPoint();
+
+                        vector pW = nearPt - pos;
 
                         scalar normalAlignment = normal & pW/mag(pW);
 
-                        allWallInteractionSites.append(nearest.rawPoint());
+                        if (normalAlignment > cosPhiMinFlatWall)
-
-                        if (normalAlignment > 1 - SMALL)
                         {
-                            flatWallInteractionSites.append(nearest.rawPoint());
+                            // Guard against a flat interaction being
-                        }
+                            // present on the boundary of two or more
-                    }
+                            // faces, which would create duplicate contact
-                }
+                            // points. Duplicates are discarded.
-            }
+                            if
-
+                            (
-            Pout<< flatWallInteractionSites << endl;
+                                !duplicatePointInList
-
+                                (
-            forAll(flatWallInteractionSites, siteI)
+                                    flatSites,
+                                    nearPt,
+                                    sqr(p.r()*flatWallDuplicateExclusion)
+                                )
+                            )
                             {
+                                flatSites.append(nearPt);
 
-                scalar nu = this->owner().constProps().poissonsRatio();
+                                flatSiteExclusionDistancesSqr.append
+                                (
+                                    sqr(p.r()) - sqr(nearest.distance())
+                                );
+                            }
+                        }
+                        else
+                        {
+                            otherSites.append(nearPt);
 
-                scalar E = this->owner().constProps().youngsModulus();
+                            otherSiteDistances.append(nearest.distance());
-
-                scalar b = 1.5;
-
-                scalar alpha = 0.52;
-
-                vector r_PW = pt - flatWallInteractionSites[siteI];
-
-                scalar normalOverlapMag = p.r() - mag(r_PW);
-
-                vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
-
-                scalar kN = (4.0/3.0)*sqrt(p.r())*E/(2.0*(1.0 - sqr(nu)));
-
-                scalar etaN = alpha*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
-
-                vector fN_PW =
-                    rHat_PW
-                   *(kN*pow(normalOverlapMag, b) - etaN*(p.U() & rHat_PW));
-
-                p.f() += fN_PW;
-
-                Pout<< "Wall force " << fN_PW << endl;
                         }
                     }
                 }
+            }
+
+            // All flat interaction sites found, now classify the
+            // other sites as being in range of a flat interaction, or
+            // a sharp interaction, being aware of not duplicating the
+            // sharp interaction sites.
+
+            // The "other" sites need to evaluated in order of
+            // ascending distance to their nearest point so that
+            // grouping occurs around the closest in any group
+
+            labelList sortedOtherSiteIndices;
+
+            sortedOrder(otherSiteDistances, sortedOtherSiteIndices);
+
+            forAll(sortedOtherSiteIndices, siteI)
+            {
+                label orderedIndex = sortedOtherSiteIndices[siteI];
+
+                const point& otherPt = otherSites[orderedIndex];
+
+                if
+                (
+                    !duplicatePointInList
+                    (
+                        flatSites,
+                        otherPt,
+                        flatSiteExclusionDistancesSqr
+                    )
+                )
+                {
+                    // Not in range of a flat interaction, must be a
+                    // sharp interaction.
+
+                    if
+                    (
+                        !duplicatePointInList
+                        (
+                            sharpSites,
+                            otherPt,
+                            sharpSiteExclusionDistancesSqr
+                        )
+                    )
+                    {
+                        sharpSites.append(otherPt);
+
+                        sharpSiteExclusionDistancesSqr.append
+                        (
+                            sqr(p.r()) - sqr(otherSiteDistances[orderedIndex])
+                        );
+                    }
+                }
+            }
+
+            evaluateWall(p, flatSites, sharpSites);
+        }
+    }
+}
+
+
+template<class CloudType>
+bool Foam::PairCollision<CloudType>::duplicatePointInList
+(
+    const DynamicList<point>& existingPoints,
+    const point& pointToTest,
+    scalar duplicateRangeSqr
+) const
+{
+    forAll(existingPoints, i)
+    {
+        if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr)
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+template<class CloudType>
+bool Foam::PairCollision<CloudType>::duplicatePointInList
+(
+    const DynamicList<point>& existingPoints,
+    const point& pointToTest,
+    const scalarList& duplicateRangeSqr
+) const
+{
+    forAll(existingPoints, i)
+    {
+        if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr[i])
+        {
+            return true;
+        }
+    }
+
+    return false;
 }
 
 
@@ -323,6 +455,18 @@ void Foam::PairCollision<CloudType>::evaluatePair
 }
 
 
+template<class CloudType>
+void Foam::PairCollision<CloudType>::evaluateWall
+(
+    typename CloudType::parcelType& p,
+    const List<point>& flatSites,
+    const List<point>& sharpSites
+) const
+{
+    wallModel_->evaluateWall(p, flatSites, sharpSites);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -342,6 +486,14 @@ Foam::PairCollision<CloudType>::PairCollision
             this->owner()
         )
     ),
+    wallModel_
+    (
+        WallModel<CloudType>::New
+        (
+            this->coeffDict(),
+            this->owner()
+        )
+    ),
     il_
     (
         owner.mesh(),

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

@@ -48,6 +48,9 @@ namespace Foam
 template<class CloudType>
 class PairModel;
 
+template<class CloudType>
+class WallModel;
+
 
 /*---------------------------------------------------------------------------*\
                      Class PairCollision Declaration
@@ -58,6 +61,19 @@ class PairCollision
 :
     public CollisionModel<CloudType>
 {
+    // Static data
+
+        //- Tolerance to determine flat wall interactions
+        static scalar cosPhiMinFlatWall;
+
+        //- Distance to multiply the distance to a wall with within
+        //  which a duplicate flat wall interaction is defined.  If
+        //  cosPhiFlatWall = 1 - e, this = sqrt(2e - e^2) < sqrt(3e)
+        //  for all e < 1, giving a safe tolerance (i.e. no double
+        //  interactions next to each other.)
+        static scalar flatWallDuplicateExclusion;
+
+
     // Private data
 
         //- Cell occupancy information for each parcel
@@ -66,6 +82,9 @@ class PairCollision
         //- PairModel to calculate the interaction between two parcels
         autoPtr<PairModel<CloudType> > pairModel_;
 
+        //- WallModel to calculate the interaction between the parcel and walls
+        autoPtr<WallModel<CloudType> > wallModel_;
+
         //- Interactions lists determining which cells are in
         //  interaction range of each other
         InteractionLists<typename CloudType::parcelType> il_;
@@ -85,6 +104,20 @@ class PairCollision
         //- Interactions with walls
         void wallInteraction();
 
+        bool duplicatePointInList
+        (
+            const DynamicList<point>& existingPoints,
+            const point& pointToTest,
+            scalar duplicateRangeSqr
+        ) const;
+
+        bool duplicatePointInList
+        (
+            const DynamicList<point>& existingPoints,
+            const point& pointToTest,
+            const scalarList& duplicateRangeSqr
+        ) const;
+
         //- Post collision tasks
         void postInteraction();
 
@@ -98,6 +131,14 @@ class PairCollision
             typename CloudType::parcelType& pB
         ) const;
 
+        //- Calculate the wall forces on a parcel
+        void evaluateWall
+        (
+            typename CloudType::parcelType& p,
+            const List<point>& flatSites,
+            const List<point>& sharpSites
+        ) const;
+
 
 public:
 

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

@@ -24,12 +24,12 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "SpringSliderDashpot.H"
+#include "PairSpringSliderDashpot.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template <class CloudType>
-void Foam::SpringSliderDashpot<CloudType>::findMinMaxProperties
+void Foam::PairSpringSliderDashpot<CloudType>::findMinMaxProperties
 (
     scalar& RMin,
     scalar& rhoMax,
@@ -45,6 +45,7 @@ void Foam::SpringSliderDashpot<CloudType>::findMinMaxProperties
         const typename CloudType::parcelType& p = iter();
 
         // Finding minimum diameter to avoid excessive arithmetic
+
         RMin = min(p.d(), RMin);
 
         rhoMax = max(p.rho(), rhoMax);
@@ -61,9 +62,11 @@ void Foam::SpringSliderDashpot<CloudType>::findMinMaxProperties
     // then rMin into minimum R,
     //     1/RMin = 1/rMin + 1/rMin,
     //     RMin = rMin/2 = dMin/4
+
     RMin /= 4.0;
 
     // Multiply by two to create the worst-case relative velocity
+
     UMagMax = 2*UMagMax;
 }
 
@@ -71,7 +74,7 @@ void Foam::SpringSliderDashpot<CloudType>::findMinMaxProperties
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template <class CloudType>
-Foam::SpringSliderDashpot<CloudType>::SpringSliderDashpot
+Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
 (
     const dictionary& dict,
     CloudType& cloud
@@ -106,21 +109,21 @@ Foam::SpringSliderDashpot<CloudType>::SpringSliderDashpot
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template <class CloudType>
-Foam::SpringSliderDashpot<CloudType>::~SpringSliderDashpot()
+Foam::PairSpringSliderDashpot<CloudType>::~PairSpringSliderDashpot()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-bool Foam::SpringSliderDashpot<CloudType>::controlsTimestep() const
+bool Foam::PairSpringSliderDashpot<CloudType>::controlsTimestep() const
 {
     return true;
 }
 
 
 template<class CloudType>
-Foam::label Foam::SpringSliderDashpot<CloudType>::nSubCycles() const
+Foam::label Foam::PairSpringSliderDashpot<CloudType>::nSubCycles() const
 {
     if (!(this->owner().size()))
     {
@@ -145,7 +148,7 @@ Foam::label Foam::SpringSliderDashpot<CloudType>::nSubCycles() const
 
 
 template<class CloudType>
-void Foam::SpringSliderDashpot<CloudType>::evaluatePair
+void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
 (
     typename CloudType::parcelType& pA,
     typename CloudType::parcelType& pB

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

@@ -23,15 +23,15 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::SpringSliderDashpot
+    Foam::PairSpringSliderDashpot
 
 Description
     Pair forces between particles colliding with a spring, slider, damper model
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpringSliderDashpot_H
+#ifndef PairSpringSliderDashpot_H
-#define SpringSliderDashpot_H
+#define PairSpringSliderDashpot_H
 
 #include "PairModel.H"
 #include "CollisionRecord.H"
@@ -41,11 +41,11 @@ Description
 namespace Foam
 {
 /*---------------------------------------------------------------------------*\
-                      Class SpringSliderDashpot Declaration
+                    Class PairSpringSliderDashpot Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CloudType>
-class SpringSliderDashpot
+class PairSpringSliderDashpot
 :
     public PairModel<CloudType>
 {
@@ -90,17 +90,17 @@ class SpringSliderDashpot
 public:
 
     //- Runtime type information
-    TypeName("SpringSliderDashpot");
+    TypeName("PairSpringSliderDashpot");
 
 
     // Constructors
 
         //- Construct from dictionary
-        SpringSliderDashpot(const dictionary& dict, CloudType& cloud);
+        PairSpringSliderDashpot(const dictionary& dict, CloudType& cloud);
 
 
     //- Destructor
-    virtual ~SpringSliderDashpot();
+    virtual ~PairSpringSliderDashpot();
 
 
     // Member Functions
@@ -130,7 +130,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "SpringSliderDashpot.C"
+#   include "PairSpringSliderDashpot.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/NewWallModel.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "WallModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::WallModel<CloudType> >
+Foam::WallModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word WallModelType(dict.lookup("WallModel"));
+
+    Info<< "Selecting WallModel " << WallModelType
+        << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(WallModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "WallModel<CloudType>::New"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "Unknown WallModelType type "
+            << WallModelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid WallModel types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError);
+    }
+
+    return autoPtr<WallModel<CloudType> >(cstrIter()(dict, owner));
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/WallModel.C

@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "WallModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::WallModel<CloudType>::WallModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    dict_(dict),
+    owner_(owner),
+    coeffDict_(dict.subDict(type + "Coeffs"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::WallModel<CloudType>::~WallModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+const CloudType&
+Foam::WallModel<CloudType>::owner() const
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary& Foam::WallModel<CloudType>::dict() const
+{
+    return dict_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary&
+Foam::WallModel<CloudType>::coeffDict() const
+{
+    return coeffDict_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NewWallModel.C"
+
+// ************************************************************************* //

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

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::WallModel
+
+Description
+    Templated wall interaction class
+
+SourceFiles
+    WallModel.C
+    NewWallModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef WallModel_H
+#define WallModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class WallModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class WallModel
+{
+    // Private data
+
+        //- The cloud dictionary
+        const dictionary& dict_;
+
+        //- Reference to the owner cloud class
+        CloudType& owner_;
+
+        //- The coefficients dictionary
+        const dictionary coeffDict_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("WallModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        WallModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct from components
+        WallModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+
+    //- Destructor
+    virtual ~WallModel();
+
+
+    //- Selector
+    static autoPtr<WallModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Access
+
+        //- Return the owner cloud object
+        const CloudType& owner() const;
+
+        //- Return the dictionary
+        const dictionary& dict() const;
+
+        //- Return the coefficients dictionary
+        const dictionary& coeffDict() const;
+
+
+    // Member Functions
+
+        //- Whether the WallModel has a timestep limit that will
+        //  require subCycling
+        virtual bool controlsTimestep() const = 0;
+
+        //- For WallModels that control the timestep, calculate the
+        //  number of subCycles needed to satisfy the minimum
+        //  allowable timestep
+        virtual label nSubCycles() const = 0;
+
+        //- Calculate the wall interaction between parcels
+        virtual void evaluateWall
+        (
+            typename CloudType::parcelType& p,
+            const List<point>& flatSites,
+            const List<point>& sharpSites
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeWallModel(CloudType)                                              \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(WallModel<CloudType>, 0);             \
+                                                                              \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        WallModel<CloudType>,                                                 \
+        dictionary                                                            \
+    );
+
+
+#define makeWallModelType(SS, CloudType, ParcelType)                          \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(SS<CloudType<ParcelType> >, 0);       \
+                                                                              \
+    WallModel<CloudType<ParcelType> >::                                       \
+        adddictionaryConstructorToTable<SS<CloudType<ParcelType> > >          \
+            add##SS##CloudType##ParcelType##ConstructorToTable_;
+
+
+#define makeWallModelThermoType(SS, CloudType, ParcelType, ThermoType)        \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        SS<CloudType<ParcelType<ThermoType> > >,                              \
+        0                                                                     \
+    );                                                                        \
+                                                                              \
+    WallModel<CloudType<ParcelType<ThermoType> > >::                          \
+        adddictionaryConstructorToTable                                       \
+            <SS<CloudType<ParcelType<ThermoType> > > >                        \
+            add##SS##CloudType##ParcelType##ThermoType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "WallModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "WallSpringSliderDashpot.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template <class CloudType>
+void Foam::WallSpringSliderDashpot<CloudType>::findMinMaxProperties
+(
+    scalar& RMin,
+    scalar& rhoMax,
+    scalar& UMagMax
+) const
+{
+    RMin = VGREAT;
+    rhoMax = -VGREAT;
+    UMagMax = -VGREAT;
+
+    forAllConstIter(typename CloudType, this->owner(), iter)
+    {
+        const typename CloudType::parcelType& p = iter();
+
+        // Finding minimum diameter to avoid excessive arithmetic
+        RMin = min(p.d(), RMin);
+
+        rhoMax = max(p.rho(), rhoMax);
+
+        UMagMax = max
+        (
+            mag(p.U()) + mag(p.omega())*p.r(),
+            UMagMax
+        );
+    }
+
+    // Transform the minimum diameter into minimum radius
+    //     rMin = dMin/2
+    // then rMin into minimum R,
+    //     1/RMin = 1/rMin + 1/rMin,
+    //     RMin = rMin/2 = dMin/4
+
+    RMin /= 4.0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    WallModel<CloudType>(dict, cloud, typeName),
+    E_(dimensionedScalar(this->coeffDict().lookup("youngsModulus")).value()),
+    nu_(dimensionedScalar(this->coeffDict().lookup("poissonsRatio")).value()),
+    alpha_(dimensionedScalar(this->coeffDict().lookup("alpha")).value()),
+    b_(dimensionedScalar(this->coeffDict().lookup("b")).value()),
+    mu_(dimensionedScalar(this->coeffDict().lookup("mu")).value()),
+    collisionResolutionSteps_
+    (
+        readScalar
+        (
+            this->coeffDict().lookup("collisionResolutionSteps")
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::WallSpringSliderDashpot<CloudType>::~WallSpringSliderDashpot()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::WallSpringSliderDashpot<CloudType>::controlsTimestep() const
+{
+    return true;
+}
+
+
+template<class CloudType>
+Foam::label Foam::WallSpringSliderDashpot<CloudType>::nSubCycles() const
+{
+    // if (!(this->owner().size()))
+    // {
+    //     return 1;
+    // }
+
+    // scalar RMin;
+    // scalar rhoMax;
+    // scalar UMagMax;
+
+    // findMinMaxProperties(RMin, rhoMax, UMagMax);
+
+    // // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
+    // scalar minCollisionDeltaT =
+    //     5.429675
+    //    *RMin
+    //    *pow(rhoMax/(Estar_*sqrt(UMagMax) + VSMALL), 0.4)
+    //    /collisionResolutionSteps_;
+
+    // return ceil(this->owner().time().deltaT().value()/minCollisionDeltaT);
+
+    return 1;
+}
+
+
+template<class CloudType>
+void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
+(
+    typename CloudType::parcelType& p,
+    const List<point>& flatSites,
+    const List<point>& sharpSites
+) const
+{
+    scalar pNu = this->owner().constProps().poissonsRatio();
+
+    scalar pE = this->owner().constProps().youngsModulus();
+
+    scalar Estar = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu_))/E_);
+
+    scalar kN = (4.0/3.0)*sqrt(p.r())*Estar;
+
+    forAll(flatSites, siteI)
+    {
+        vector r_PW = p.position() - flatSites[siteI];
+
+        scalar normalOverlapMag = p.r() - mag(r_PW);
+
+        vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
+
+        scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
+
+        vector fN_PW =
+            rHat_PW
+           *(kN*pow(normalOverlapMag, b_) - etaN*(p.U() & rHat_PW));
+
+        p.f() += fN_PW;
+    }
+
+    // Treating sharp sites like flat sites
+
+    forAll(sharpSites, siteI)
+    {
+        vector r_PW = p.position() - sharpSites[siteI];
+
+        scalar normalOverlapMag = p.r() - mag(r_PW);
+
+        vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
+
+        scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
+
+        vector fN_PW =
+            rHat_PW
+           *(kN*pow(normalOverlapMag, b_) - etaN*(p.U() & rHat_PW));
+
+        p.f() += fN_PW;
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::WallSpringSliderDashpot
+
+Description
+    Forces between particles and walls, interacting with a spring,
+    slider, damper model
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef WallSpringSliderDashpot_H
+#define WallSpringSliderDashpot_H
+
+#include "WallModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                    Class WallSpringSliderDashpot Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class WallSpringSliderDashpot
+:
+    public WallModel<CloudType>
+{
+    // Private data
+
+        //- Young's modulus of the wall
+        scalar E_;
+
+        //- Poisson's ratio of the wall
+        scalar nu_;
+
+        //- alpha-coefficient, related to coefficient of restitution
+        scalar alpha_;
+
+        //- Spring power (b = 1 for linear, b = 3/2 for Hertzian)
+        scalar b_;
+
+        //- Coefficient of friction in for tangential sliding
+        scalar mu_;
+
+        //- The number of steps over which to resolve the minimum
+        //  harmonic approximation of the collision period
+        scalar collisionResolutionSteps_;
+
+
+    // Private Member Functions
+
+        //- Find the appropriate properties for determining the minimum
+        //- allowable timestep
+        void findMinMaxProperties
+        (
+            scalar& RMin,
+            scalar& rhoMax,
+            scalar& vMagMax
+        ) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("WallSpringSliderDashpot");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        WallSpringSliderDashpot(const dictionary& dict, CloudType& cloud);
+
+
+    //- Destructor
+    virtual ~WallSpringSliderDashpot();
+
+
+    // Member Functions
+
+        //- Whether the WallModel has a timestep limit that will
+        //  require subCycling
+        virtual bool controlsTimestep() const;
+
+        //- For WallModels that control the timestep, calculate the
+        //  number of subCycles needed to satisfy the minimum
+        //  allowable timestep
+        virtual label nSubCycles() const;
+
+        //- Calculate the wall interaction between parcels
+        virtual void evaluateWall
+        (
+            typename CloudType::parcelType& p,
+            const List<point>& flatSites,
+            const List<point>& sharpSites
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "WallSpringSliderDashpot.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //