git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2354 f3b2605a-c512-4ea7-a41b-209d697bcdaa

doc/Eqs/pair_gran_hertz.tex

@@ -7,7 +7,7 @@ $$
      \sqrt{\delta} \sqrt{\frac{R_i R_j}{R_i + R_j}} 
      \Big[ (k_n \delta \mathbf{n}_{ij} -  
        m_{\mbox{\scriptsize{eff}}} \: \gamma_n \mathbf{ v}_n) -
-       (k_t \delta\mathbf{ \Delta s}_t +
+       (k_t \mathbf{ \Delta s}_t +
        m_{\mbox{\scriptsize{eff}}} \: \gamma_t \mathbf{v}_t) \Big]
 $$
 

doc/Eqs/pair_gran_hooke.tex

@@ -5,7 +5,7 @@
 $$
    F_{hk} = (k_n \delta \mathbf{n}_{ij} -  
    m_{\mbox{\scriptsize{eff}}} \gamma_n\mathbf{ v}_n) - 
-   (k_t \delta\mathbf{ \Delta s}_t +
+   (k_t \mathbf{ \Delta s}_t +
    m_{\mbox{\scriptsize{eff}}} \gamma_t \mathbf{v}_t)
 $$
 

doc/Section_commands.html

@@ -366,7 +366,7 @@ full description:
 <TR ALIGN="center"><TD ><A HREF = "pair_coul.html">coul/cut</A></TD><TD ><A HREF = "pair_coul.html">coul/debye</A></TD><TD ><A HREF = "pair_coul.html">coul/long</A></TD><TD ><A HREF = "pair_dipole.html">dipole/cut</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_dpd.html">dpd</A></TD><TD ><A HREF = "pair_eam.html">eam</A></TD><TD ><A HREF = "pair_eam.html">eam/opt</A></TD><TD ><A HREF = "pair_eam.html">eam/alloy</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_eam.html">eam/alloy/opt</A></TD><TD ><A HREF = "pair_eam.html">eam/fs</A></TD><TD ><A HREF = "pair_eam.html">eam/fs/opt</A></TD><TD ><A HREF = "pair_gayberne.html">gayberne</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_gran.html">gran/hertzian</A></TD><TD ><A HREF = "pair_gran.html">gran/history</A></TD><TD ><A HREF = "pair_gran.html">gran/no_history</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_gran_hertz_history.html">gran/hertz/history</A></TD><TD ><A HREF = "pair_gran_hooke.html">gran/hooke</A></TD><TD ><A HREF = "pair_gran_hooke_history.html">gran/hooke/history</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm/implicit</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long/opt</A></TD><TD ><A HREF = "pair_class2.html">lj/class2</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_class2.html">lj/class2/coul/cut</A></TD><TD ><A HREF = "pair_class2.html">lj/class2/coul/long</A></TD><TD ><A HREF = "pair_lj.html">lj/cut</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/opt</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/coul/cut</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/debye</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long/tip4p</A></TD></TR>

doc/Section_commands.txt

@@ -504,9 +504,9 @@ full description:
 "eam/fs"_pair_eam.html,
 "eam/fs/opt"_pair_eam.html,
 "gayberne"_pair_gayberne.html,
-"gran/hertzian"_pair_gran.html,
+"gran/hertz/history"_pair_gran_hertz_history.html,
-"gran/history"_pair_gran.html,
+"gran/hooke"_pair_gran_hooke.html,
-"gran/no_history"_pair_gran.html,
+"gran/hooke/history"_pair_gran_hooke_history.html,
 "lj/charmm/coul/charmm"_pair_charmm.html,
 "lj/charmm/coul/charmm/implicit"_pair_charmm.html,
 "lj/charmm/coul/long"_pair_charmm.html,

doc/fix_wall_gran.html

@@ -13,26 +13,33 @@
 </H3>
 <P><B>Syntax:</B>
 </P>
-<PRE>fix ID group-ID wall/gran wallstyle args keyword values ... 
+<PRE>fix ID group-ID wall/gran Kn Kt gamma_n gamma_t xmu dampflag wallstyle args keyword values ... 
 </PRE>
 <UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command 
 
 <LI>wall/gran = style name of this fix command 
 
-<LI>style = <I>xplane</I> or <I>yplane</I> or <I>zplane</I> or <I>zcylinder</I> 
+<LI>Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) 
+
+<LI>Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) 
+
+<LI>gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) 
+
+<LI>gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) 
+
+<LI>xmu = static yield criterion (unitless fraction between 0.0 and 1.0) 
+
+<LI>dampflag = 0 or 1 if tangential damping force is excluded or included 
+
+<LI>wallstyle = <I>xplane</I> or <I>yplane</I> or <I>zplane</I> or <I>zcylinder</I> 
 
 <LI>args = list of arguments for a particular style 
 
-<PRE>  <I>xplane</I> or <I>yplane</I> or <I>zplane</I> args = lo hi gamma xmu
+<LI>  <I>xplane</I> or <I>yplane</I> or <I>zplane</I> args = lo hi
     lo,hi = position of lower and upper plane (distance units), either can be NULL)
-    gamman = damping coeff for normal direction collisions with wall
+  <I>zcylinder</I> args = radius
-    xmu = friction coeff for the wall
-  <I>zcylinder</I> args = radius gamma xmu
     radius = cylinder radius (distance units)
-    gamman = damping coeff for normal direction collisions with wall
+zero or more keyword/value pairs may be appended to args 
-    xmu = friction coeff for the wall 
-</PRE>
-<LI>zero or more keyword/value pairs may be appended to args 
 
 <LI>keyword = <I>wiggle</I> or <I>shear</I> 
 
@@ -48,9 +55,9 @@
 </UL>
 <P><B>Examples:</B>
 </P>
-<PRE>fix 1 all wall/gran xplane -10.0 10.0 50.0 0.5
+<PRE>fix 1 all wall/gran 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 
-fix 2 all wall/gran zcylinder 15.0 50.0 0.5 wiggle z 3.0 2.0
+fix 1 all wall/gran 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL 
-fix 1 all wall/gran zplane 0.0 NULL 100.0 0.5 
+fix 2 all wall/gran 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -58,13 +65,48 @@ fix 1 all wall/gran zplane 0.0 NULL 100.0 0.5
 wall.  All particles in the group interact with the wall when they are
 close enough to touch it.
 </P>
+<P>The first set of parameters (Kn, Kt, gamma_n, gamma_t, xmu, and
+dampflag) have the same meaning as those specified with the
+<A HREF = "pair_gran.html">pair_style granular</A> force fields.  This means a NULL
+can be used for either Kt or gamma_t as described on that page.  If a
+NULL is used for Kt, then a default value is used where Kt = 2/7 Kn.
+If a NULL is used for gamma_t, then a default value is used where
+gamma_t = 1/2 gamma_n.
+</P>
+<P>The nature of the wall/particle interactions are determined by which
+pair_style is used in your input script: <I>hooke</I>, <I>hooke/history</I>, or
+<I>hertz/history</I>.  The equation for the force between the wall and
+particles touching it is the same as the corresponding equation on the
+<A HREF = "pair_gran.html">pair_style granular</A> doc page, in the limit of one of
+the two particles going to infinite radius and mass (flat wall).
+I.e. delta = radius - r = overlap of particle with wall, m_eff = mass
+of particle, and sqrt(RiRj/Ri+Rj) becomes sqrt(radius of particle).
+The units for Kn, Kt, gamma_n, and gamma_t are as described on that
+doc page.  The meaning of xmu and dampflag are also as described on
+that page.  Note that you can choose different values for these 6
+wall/particle coefficients than for particle/particle interactions, if
+you wish your wall to interact differently with the particles, e.g. if
+the wall is a different material.
+</P>
+<P>IMPORTANT NOTE: As discussed on the doc page for <A HREF = "pair_gran.html">pair_style
+granular</A>, versions of LAMMPS before 9Jan09 used a
+different equation for Hertzian interactions.  This means Hertizian
+wall/particle interactions have also changed.  They now include a
+sqrt(radius) term which was not present before.  Also the previous
+versions used Kn and Kt from the pairwise interaction and hardwired
+dampflag to 1, rather than letting them be specified directly.  This
+means you can set the values of the wall/particle coefficients
+appropriately in the current code to reproduce the results of a
+prevoius Hertzian monodisperse calculation.  For example, for the
+common case of a monodisperse system with particles of diameter 1, Kn,
+Kt, gamma_n, and gamma_s should be set sqrt(2.0) larger than they were
+previously.
+</P>
 <P>The <I>wallstyle</I> can be planar or cylindrical.  The 3 planar options
 specify a pair of walls in a dimension.  Wall positions are given by
 <I>lo</I> and <I>hi</I>.  Either of the values can be specified as NULL if a
 single wall is desired.  For a <I>zcylinder</I> wallstyle, the cylinder's
 axis is at x = y = 0.0, and the radius of the cylinder is specified.
-For all wallstyles, a damping and friction coefficient for
-particle-wall interactions are also specified.
 </P>
 <P>Optionally, the wall can be moving, if the <I>wiggle</I> or <I>shear</I>
 keywords are appended.  Both keywords cannot be used together.
@@ -126,7 +168,7 @@ LAMMPS</A> section for more info.
 </P>
 <P><B>Related commands:</B>
 </P>
-<P><A HREF = "fix_wiggle.html">fix_wiggle</A>
+<P><A HREF = "fix_wiggle.html<A HREF = "pair_gran.html">>fix_wiggle</A>, pair_style granular</A>
 </P>
 <P><B>Default:</B> none
 </P>

doc/fix_wall_gran.txt

@@ -10,20 +10,22 @@ fix wall/gran command :h3
 
 [Syntax:]
 
-fix ID group-ID wall/gran wallstyle args keyword values ... :pre
+fix ID group-ID wall/gran Kn Kt gamma_n gamma_t xmu dampflag wallstyle args keyword values ... :pre
 
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
 wall/gran = style name of this fix command :l
-style = {xplane} or {yplane} or {zplane} or {zcylinder} :l
+Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) :l
+Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) :l
+gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) :l
+gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) :l
+xmu = static yield criterion (unitless fraction between 0.0 and 1.0) :l
+dampflag = 0 or 1 if tangential damping force is excluded or included :l
+wallstyle = {xplane} or {yplane} or {zplane} or {zcylinder} :l
 args = list of arguments for a particular style :l
-  {xplane} or {yplane} or {zplane} args = lo hi gamma xmu
+  {xplane} or {yplane} or {zplane} args = lo hi
     lo,hi = position of lower and upper plane (distance units), either can be NULL)
-    gamman = damping coeff for normal direction collisions with wall
+  {zcylinder} args = radius
-    xmu = friction coeff for the wall
-  {zcylinder} args = radius gamma xmu
     radius = cylinder radius (distance units)
-    gamman = damping coeff for normal direction collisions with wall
-    xmu = friction coeff for the wall :pre
 zero or more keyword/value pairs may be appended to args :l
 keyword = {wiggle} or {shear} :l
   {wiggle} values = dim amplitude period
@@ -37,9 +39,9 @@ keyword = {wiggle} or {shear} :l
 
 [Examples:]
 
-fix 1 all wall/gran xplane -10.0 10.0 50.0 0.5
+fix 1 all wall/gran 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 
-fix 2 all wall/gran zcylinder 15.0 50.0 0.5 wiggle z 3.0 2.0
+fix 1 all wall/gran 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL 
-fix 1 all wall/gran zplane 0.0 NULL 100.0 0.5 :pre
+fix 2 all wall/gran 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 :pre
 
 [Description:]
 
@@ -47,13 +49,48 @@ Bound the simulation domain of a granular system with a frictional
 wall.  All particles in the group interact with the wall when they are
 close enough to touch it.
 
+The first set of parameters (Kn, Kt, gamma_n, gamma_t, xmu, and
+dampflag) have the same meaning as those specified with the
+"pair_style granular"_pair_gran.html force fields.  This means a NULL
+can be used for either Kt or gamma_t as described on that page.  If a
+NULL is used for Kt, then a default value is used where Kt = 2/7 Kn.
+If a NULL is used for gamma_t, then a default value is used where
+gamma_t = 1/2 gamma_n.
+
+The nature of the wall/particle interactions are determined by which
+pair_style is used in your input script: {hooke}, {hooke/history}, or
+{hertz/history}.  The equation for the force between the wall and
+particles touching it is the same as the corresponding equation on the
+"pair_style granular"_pair_gran.html doc page, in the limit of one of
+the two particles going to infinite radius and mass (flat wall).
+I.e. delta = radius - r = overlap of particle with wall, m_eff = mass
+of particle, and sqrt(RiRj/Ri+Rj) becomes sqrt(radius of particle).
+The units for Kn, Kt, gamma_n, and gamma_t are as described on that
+doc page.  The meaning of xmu and dampflag are also as described on
+that page.  Note that you can choose different values for these 6
+wall/particle coefficients than for particle/particle interactions, if
+you wish your wall to interact differently with the particles, e.g. if
+the wall is a different material.
+
+IMPORTANT NOTE: As discussed on the doc page for "pair_style
+granular"_pair_gran.html, versions of LAMMPS before 9Jan09 used a
+different equation for Hertzian interactions.  This means Hertizian
+wall/particle interactions have also changed.  They now include a
+sqrt(radius) term which was not present before.  Also the previous
+versions used Kn and Kt from the pairwise interaction and hardwired
+dampflag to 1, rather than letting them be specified directly.  This
+means you can set the values of the wall/particle coefficients
+appropriately in the current code to reproduce the results of a
+prevoius Hertzian monodisperse calculation.  For example, for the
+common case of a monodisperse system with particles of diameter 1, Kn,
+Kt, gamma_n, and gamma_s should be set sqrt(2.0) larger than they were
+previously.
+
 The {wallstyle} can be planar or cylindrical.  The 3 planar options
 specify a pair of walls in a dimension.  Wall positions are given by
 {lo} and {hi}.  Either of the values can be specified as NULL if a
 single wall is desired.  For a {zcylinder} wallstyle, the cylinder's
 axis is at x = y = 0.0, and the radius of the cylinder is specified.
-For all wallstyles, a damping and friction coefficient for
-particle-wall interactions are also specified.
 
 Optionally, the wall can be moving, if the {wiggle} or {shear}
 keywords are appended.  Both keywords cannot be used together.
@@ -115,6 +152,6 @@ Any dimension (xyz) that has a granular wall must be non-periodic.
 
 [Related commands:]
 
-"fix_wiggle"_fix_wiggle.html
+"fix_wiggle"_fix_wiggle.html, pair_style granular"_pair_gran.html
 
 [Default:] none

doc/pair_gran.html

@@ -17,21 +17,35 @@
 </H3>
 <P><B>Syntax:</B>
 </P>
-<PRE>pair_style style Kn gamma_n xmu dampflag 
+<PRE>pair_style style Kn Kt gamma_n gamma_t xmu dampflag 
 </PRE>
 <UL><LI>style = <I>gran/hooke</I> or <I>gran/hooke/history</I> or <I>gran/hertz/history</I> 
 
-<LI>Kn = spring constant for particle repulsion (see units below) 
+<LI>Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) 
 
-<LI>gamma_n = damping coefficient for collisions in normal direction (see units below) 
+<LI>Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) 
+
+<LI>gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) 
+
+<LI>gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) 
 
 <LI>xmu = static yield criterion (unitless fraction between 0.0 and 1.0) 
 
 <LI>dampflag = 0 or 1 if tangential damping force is excluded or included 
 </UL>
+<P>IMPORTANT NOTE: Versions of LAMMPS before 9Jan09 had different style
+names for granular force fields.  This is to emphasize the fact that
+the Hertzian equation has changed to model polydispersity more
+accurately.  A side effect of the change is that the Kn, Kt, gamma_n,
+and gamma_t coefficients in the pair_style command must be specified
+with different values in order to reproduce calculations made with
+earlier versions of LAMMPS, even for monodisperse systems.  See the
+NOTE below for details.
+</P>
 <P><B>Examples:</B>
 </P>
-<PRE>pair_style gran/history 200000.0 50.0 0.5 1 
+<PRE>pair_style gran/hooke/history 200000.0 NULL 50.0 NULL 0.5 1
+pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -65,51 +79,67 @@ if <I>dampflag</I> is set to 0.
 </P>
 <UL><LI>delta = d - r = overlap distance of 2 particles
 <LI>Kn = elastic constant for normal contact
-<LI>Kt = elastic constant for tangential contact = 2/7 of Kn
+<LI>Kt = elastic constant for tangential contact
 <LI>gamma_n = viscoelastic damping constant for normal contact
-<LI>gamma_t = viscoelastic damping constant for tangential contact = 1/2 of gamma_n
+<LI>gamma_t = viscoelastic damping constant for tangential contact
 <LI>m_eff = Mi Mj / (Mi + Mj) = effective mass of 2 particles of mass Mi and Mj
 <LI>Delta St = tangential displacement vector between 2 spherical particles       which is truncated to satisfy a frictional yield criterion
 <LI>n_ij = unit vector along the line connecting the centers of the 2 particles
 <LI>Vn = normal component of the relative velocity of the 2 particles
 <LI>Vt = tangential component of the relative velocity of the 2 particles 
 </UL>
-<P>The Kn and gamma_n coefficients are specified as parameters to the
+<P>The Kn, Kt, gamma_n, and gamma_t coefficients are specified as
-pair_style command.  The interpretation and units for these
+parameters to the pair_style command.  If a NULL is used for Kt, then
-coefficients are different in the Hookean versus Hertzian formulas.
+a default value is used where Kt = 2/7 Kn.  If a NULL is used for
+gamma_t, then a default value is used where gamma_t = 1/2 gamma_n.
+</P>
+<P>The interpretation and units for these 4 coefficients are different in
+the Hookean versus Hertzian equations.
 </P>
 <P>The Hookean model is one where the normal push-back force for two
 overlapping particles is a linear function of the overlap distance.
 Thus the specified Kn is in units of (force/distance).  Note that this
 push-back force is independent of absolute particle size (in the
-monodisperse case) or the relative sizes of the two particles (in the
+monodisperse case) and of the relative sizes of the two particles (in
-polydisperse case).  This model also applies to the other terms in the
+the polydisperse case).  This model also applies to the other terms in
-force equation so that the specified gamma_n is in units of (1/time),
+the force equation so that the specified gamma_n is in units of
-Kt is in units of (force/distance), and gamma_t is in units of
+(1/time), Kt is in units of (force/distance), and gamma_t is in units
-(1/time).
+of (1/time).
 </P>
 <P>The Hertzian model is one where the normal push-back force for two
 overlapping particles is proportional to the area of overlap of the
 two particles, and is thus a non-linear function of overlap distance.
-Thus the specified Kn is in units of (force/area).  The effects of
+Thus Kn has units of force per area and is thus specified in units of
-absolute particle size (monodispersity) and relative size
+(pressure).  The effects of absolute particle size (monodispersity)
-(polydispersity) are captured in the radii-dependent pre-factors.
+and relative size (polydispersity) are captured in the radii-dependent
-When these pre-factors are carried through to the other terms in the
+pre-factors.  When these pre-factors are carried through to the other
-force equation it means that the specified gamma_n is in units of
+terms in the force equation it means that the specified gamma_n is in
-(1/time-distance), Kt is in units of (force/area), and gamma_t is in
+units of (1/(time*distance)), Kt is in units of (pressure), and
-units of (1/time-distance).
+gamma_t is in units of (1/(time*distance)).
 </P>
-<P>Note that in the Hookean case, Kn can be thought of as a spring
+<P>Note that in the Hookean case, Kn can be thought of as a linear spring
 constant with units of force/distance.  In the Hertzian case, Kn is
-like a non-linear spring constant with units of force/area, and as
+like a non-linear spring constant with units of force/area or
-shown in the <A HREF = "#Zhang">(Zhang)</A> paper, Kn = 4G / (3(1-nu)) where nu =
+pressure, and as shown in the <A HREF = "#Zhang">(Zhang)</A> paper, Kn = 4G /
-the Poisson ratio, G = shear modulus = E / (1(1+nu)), and E = Young's
+(3(1-nu)) where nu = the Poisson ratio, G = shear modulus = E /
-modulus.  Thus in the Hertzian case Kn can be set to a value that
+(1(1+nu)), and E = Young's modulus.  Similarly, Kt = 8G / (2-nu).
+Thus in the Hertzian case Kn and Kt can be set to values that
 corresponds to properties of the material being modeled.  This is also
 true in the Hookean case, except that a spring constant must be chosen
-that is appropriate for the size of particles in the model.  Since
+that is appropriate for the absolute size of particles in the model.
-relative particle sizes are not accounted for, the Hookean styles may
+Since relative particle sizes are not accounted for, the Hookean
-not be a suitable model for polydisperse systems.
+styles may not be a suitable model for polydisperse systems.
+</P>
+<P>IMPORTANT NOTE: In versions of LAMMPS before 9Jan09, the equation for
+Hertzian interactions did not include the sqrt(RiRj/Ri+Rj) term and
+thus was not as accurate for polydisperse systems.  For monodisperse
+systems, sqrt(RiRj/Ri+Rj) is a constant factor that effectively scales
+all 4 coefficients: Kn, Kt, gamma_n, gamma_t.  Thus you can set the
+values of these 4 coefficients appropriately in the current code to
+reproduce the results of a previous Hertzian monodisperse calculation.
+For example, for the common case of a monodisperse system with
+particles of diameter 1, all 4 of these coefficients should now be set
+2x larger than they were previously.
 </P>
 <P>Xmu is also specified in the pair_style command and is the upper limit
 of the tangential force through the Coulomb criterion Ft = xmu*Fn,

doc/pair_gran.txt

@@ -12,17 +12,29 @@ pair_style gran/hertz/history command :h3
 
 [Syntax:]
 
-pair_style style Kn gamma_n xmu dampflag :pre
+pair_style style Kn Kt gamma_n gamma_t xmu dampflag :pre
 
 style = {gran/hooke} or {gran/hooke/history} or {gran/hertz/history} :ulb,l
-Kn = spring constant for particle repulsion (see units below) :l
+Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) :l
-gamma_n = damping coefficient for collisions in normal direction (see units below) :l
+Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) :l
+gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) :l
+gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) :l
 xmu = static yield criterion (unitless fraction between 0.0 and 1.0) :l
 dampflag = 0 or 1 if tangential damping force is excluded or included :l,ule
 
+IMPORTANT NOTE: Versions of LAMMPS before 9Jan09 had different style
+names for granular force fields.  This is to emphasize the fact that
+the Hertzian equation has changed to model polydispersity more
+accurately.  A side effect of the change is that the Kn, Kt, gamma_n,
+and gamma_t coefficients in the pair_style command must be specified
+with different values in order to reproduce calculations made with
+earlier versions of LAMMPS, even for monodisperse systems.  See the
+NOTE below for details.
+
 [Examples:]
 
-pair_style gran/history 200000.0 50.0 0.5 1 :pre
+pair_style gran/hooke/history 200000.0 NULL 50.0 NULL 0.5 1
+pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 :pre
 
 [Description:]
 
@@ -56,9 +68,9 @@ The other quantities in the equations are as follows:
 
 delta = d - r = overlap distance of 2 particles
 Kn = elastic constant for normal contact
-Kt = elastic constant for tangential contact = 2/7 of Kn
+Kt = elastic constant for tangential contact
 gamma_n = viscoelastic damping constant for normal contact
-gamma_t = viscoelastic damping constant for tangential contact = 1/2 of gamma_n
+gamma_t = viscoelastic damping constant for tangential contact
 m_eff = Mi Mj / (Mi + Mj) = effective mass of 2 particles of mass Mi and Mj
 Delta St = tangential displacement vector between 2 spherical particles \
       which is truncated to satisfy a frictional yield criterion
@@ -66,42 +78,58 @@ n_ij = unit vector along the line connecting the centers of the 2 particles
 Vn = normal component of the relative velocity of the 2 particles
 Vt = tangential component of the relative velocity of the 2 particles :ul
 
-The Kn and gamma_n coefficients are specified as parameters to the
+The Kn, Kt, gamma_n, and gamma_t coefficients are specified as
-pair_style command.  The interpretation and units for these
+parameters to the pair_style command.  If a NULL is used for Kt, then
-coefficients are different in the Hookean versus Hertzian formulas.
+a default value is used where Kt = 2/7 Kn.  If a NULL is used for
+gamma_t, then a default value is used where gamma_t = 1/2 gamma_n.
+
+The interpretation and units for these 4 coefficients are different in
+the Hookean versus Hertzian equations.
 
 The Hookean model is one where the normal push-back force for two
 overlapping particles is a linear function of the overlap distance.
 Thus the specified Kn is in units of (force/distance).  Note that this
 push-back force is independent of absolute particle size (in the
-monodisperse case) or the relative sizes of the two particles (in the
+monodisperse case) and of the relative sizes of the two particles (in
-polydisperse case).  This model also applies to the other terms in the
+the polydisperse case).  This model also applies to the other terms in
-force equation so that the specified gamma_n is in units of (1/time),
+the force equation so that the specified gamma_n is in units of
-Kt is in units of (force/distance), and gamma_t is in units of
+(1/time), Kt is in units of (force/distance), and gamma_t is in units
-(1/time).
+of (1/time).
 
 The Hertzian model is one where the normal push-back force for two
 overlapping particles is proportional to the area of overlap of the
 two particles, and is thus a non-linear function of overlap distance.
-Thus the specified Kn is in units of (force/area).  The effects of
+Thus Kn has units of force per area and is thus specified in units of
-absolute particle size (monodispersity) and relative size
+(pressure).  The effects of absolute particle size (monodispersity)
-(polydispersity) are captured in the radii-dependent pre-factors.
+and relative size (polydispersity) are captured in the radii-dependent
-When these pre-factors are carried through to the other terms in the
+pre-factors.  When these pre-factors are carried through to the other
-force equation it means that the specified gamma_n is in units of
+terms in the force equation it means that the specified gamma_n is in
-(1/time-distance), Kt is in units of (force/area), and gamma_t is in
+units of (1/(time*distance)), Kt is in units of (pressure), and
-units of (1/time-distance).
+gamma_t is in units of (1/(time*distance)).
 
-Note that in the Hookean case, Kn can be thought of as a spring
+Note that in the Hookean case, Kn can be thought of as a linear spring
 constant with units of force/distance.  In the Hertzian case, Kn is
-like a non-linear spring constant with units of force/area, and as
+like a non-linear spring constant with units of force/area or
-shown in the "(Zhang)"_#Zhang paper, Kn = 4G / (3(1-nu)) where nu =
+pressure, and as shown in the "(Zhang)"_#Zhang paper, Kn = 4G /
-the Poisson ratio, G = shear modulus = E / (1(1+nu)), and E = Young's
+(3(1-nu)) where nu = the Poisson ratio, G = shear modulus = E /
-modulus.  Thus in the Hertzian case Kn can be set to a value that
+(1(1+nu)), and E = Young's modulus.  Similarly, Kt = 8G / (2-nu).
+Thus in the Hertzian case Kn and Kt can be set to values that
 corresponds to properties of the material being modeled.  This is also
 true in the Hookean case, except that a spring constant must be chosen
-that is appropriate for the size of particles in the model.  Since
+that is appropriate for the absolute size of particles in the model.
-relative particle sizes are not accounted for, the Hookean styles may
+Since relative particle sizes are not accounted for, the Hookean
-not be a suitable model for polydisperse systems.
+styles may not be a suitable model for polydisperse systems.
+
+IMPORTANT NOTE: In versions of LAMMPS before 9Jan09, the equation for
+Hertzian interactions did not include the sqrt(RiRj/Ri+Rj) term and
+thus was not as accurate for polydisperse systems.  For monodisperse
+systems, sqrt(RiRj/Ri+Rj) is a constant factor that effectively scales
+all 4 coefficients: Kn, Kt, gamma_n, gamma_t.  Thus you can set the
+values of these 4 coefficients appropriately in the current code to
+reproduce the results of a previous Hertzian monodisperse calculation.
+For example, for the common case of a monodisperse system with
+particles of diameter 1, all 4 of these coefficients should now be set
+2x larger than they were previously.
 
 Xmu is also specified in the pair_style command and is the upper limit
 of the tangential force through the Coulomb criterion Ft = xmu*Fn,

doc/region.html

@@ -59,7 +59,7 @@
 </P>
 <PRE>region 1 block -3.0 5.0 INF 10.0 INF INF
 region 2 sphere 0.0 0.0 0.0 5 side out
-region void cylinder y 2 3 5 -5.0 INF units box
+region void cylinder y 2 3 5 -5.0 EDGE units box
 region 1 prism 0 10 0 10 0 10 2 0 0
 region outside union 4 side1 side2 side3 side4 
 </PRE>

doc/region.txt

@@ -50,7 +50,7 @@ keyword = {side} or {units} :l
 
 region 1 block -3.0 5.0 INF 10.0 INF INF
 region 2 sphere 0.0 0.0 0.0 5 side out
-region void cylinder y 2 3 5 -5.0 INF units box
+region void cylinder y 2 3 5 -5.0 EDGE units box
 region 1 prism 0 10 0 10 0 10 2 0 0
 region outside union 4 side1 side2 side3 side4 :pre