\(\renewcommand{\AA}{\text{Å}}\)
3.16. Granular SubModel styles
In granular models, particles are spheres with a finite radius and rotational degrees of freedom as further described in the Howto granular page. Interactions between pair of particles or particles and walls may therefore depend on many different modes of motion as described in pair granular and fix wall/gran. In both cases, the exchange of forces, torques, and heat flow between two types of bodies is defined using a GranularModel class. The GranularModel class organizes the details of an interaction using a series of granular submodels each of which describe a particular interaction mode (e.g. normal forces or rolling friction). From a parent GranSubMod class, several types of submodel classes are derived:
GranSubModNormal: normal force submodel
GranSubModDamping: normal damping submodel
GranSubModTangential: tangential forces and sliding friction submodel
GranSubModRolling: rolling friction submodel
GranSubModTwisting: twisting friction submodel
GranSubModHeat: heat conduction submodel
For each type of submodel, more classes are further derived, each describing a specific implementation. For instance, from the GranSubModNormal class the GranSubModNormalHooke, GranSubModNormalHertz, and GranSubModNormalJKR classes are derived which calculate Hookean, Hertzian, or JKR normal forces, respectively. This modular structure simplifies the addition of new granular contact models as one only needs to create a new GranSubMod class without having to modify the more complex PairGranular, FixGranWall, and GranularModel classes. Most GranSubMod methods are also already defined by the parent classes so new contact models typically only require edits to a few relevant methods (e.g. methods that define coefficients and calculate forces).
Each GranSubMod class has a pointer to both the LAMMPS class and the GranularModel
class which owns it, lmp
and gm
, respectively. The GranularModel class
includes several public variables that describe the geometry/dynamics of the
contact such as

Positions of the two contacting bodies 

Velocities of the two contacting bodies 

Angular velocities of the two contacting bodies 

The displacement and normalized displacement vectors 

The distance, distance squared, and inverse distance 

The sum of particle radii 

The relative velocity vector and its normal and tangential components 

The relative rotational velocity 
These quantities, among others, are calculated in the GranularModel>check_contact()
and GranularModel>calculate_forces()
methods which can be referred to for more
details.
To create a new GranSubMod class, it is recommended that one first looks at similar GranSubMod classes. All GranSubMod classes share several general methods which may need to be defined

Optional method to define how coefficients are mixed for different atom types. By default, coefficients are mixed using a geometric mean. 

Parses coefficients to define local variables. Run once at model construction. 

Optional method to define local variables after other GranSubMod types were created. For instance, this method may be used by a tangential model that derives parameters from the normal model. 
There are also several typespecific methods

Optional method to test when particles are in contact. By default, this is when particles overlap. 

Optional method to return the distance at which particles stop interacting. By default, this is when particles no longer overlap. 

Optional method to return the radius of the contact. By default, this returns the radius of the geometric cross section. 

Optional method that defines the critical force to break the contact used by some tangential, rolling, and twisting submodels. By default, this is the current total normal force including damping. 

Required method that returns the normal contact force 

Required method that returns the normal damping force 

Required method that calculates tangential forces/torques 

Required method that calculates twisting friction forces/torques 

Required method that calculates rolling friction forces/torques 

Required method that returns the rate of heat flow 
As an example, say one wanted to create a new normal force option that consisted
of a Hookean force with a piecewise stiffness. This could be done by adding a new
set of files gran_sub_mod_custom.h
:
#ifdef GranSubMod_CLASS
// clangformat off
GranSubModStyle(hooke/piecewise,GranSubModNormalHookePiecewise,NORMAL);
// clangformat on
#else
#ifndef GRAN_SUB_MOD_CUSTOM_H_
#define GRAN_SUB_MOD_CUSTOM_H_
#include "gran_sub_mod.h"
#include "gran_sub_mod_normal.h"
namespace LAMMPS_NS {
namespace Granular_NS {
class GranSubModNormalHookePiecewise : public GranSubModNormal {
public:
GranSubModNormalHookePiecewise(class GranularModel *, class LAMMPS *);
void coeffs_to_local() override;
double calculate_forces() override;
protected:
double k1, k2, delta_switch;
};
} // namespace Granular_NS
} // namespace LAMMPS_NS
#endif /*GRAN_SUB_MOD_CUSTOM_H_ */
#endif /*GRAN_SUB_MOD_CLASS_H_ */
and gran_sub_mod_custom.cpp
#include "gran_sub_mod_custom.h"
#include "gran_sub_mod_normal.h"
#include "granular_model.h"
using namespace LAMMPS_NS;
using namespace Granular_NS;
GranSubModNormalHookePiecewise::GranSubModNormalHookePiecewise(GranularModel *gm, LAMMPS *lmp) :
GranSubModNormal(gm, lmp)
{
num_coeffs = 4;
}
/*  */
void GranSubModNormalHookePiecewise::coeffs_to_local()
{
k1 = coeffs[0];
k2 = coeffs[1];
damp = coeffs[2];
delta_switch = coeffs[3];
}
/*  */
double GranSubModNormalHookePiecewise::calculate_forces()
{
double Fne;
if (gm>delta >= delta_switch) {
Fne = k1 * delta_switch + k2 * (gm>delta  delta_switch);
} else {
Fne = k1 * gm>delta;
}
return Fne;
}