Merge pull request lammps#4409 from willzunker/mdr-rebase2

pair_style granular - MDR contact model

doc/src/fix_wall_gran.rst

@@ -222,10 +222,10 @@ restart file, so that the operation of the fix continues in an
 uninterrupted fashion.
 
 If the :code:`contacts` option is used, this fix generates a per-atom array
-with 8 columns as output, containing the contact information for owned
+with at least 8 columns as output, containing the contact information for owned
 particles (nlocal on each processor). All columns in this per-atom array will
-be zero if no contact has occurred.  The values of these columns are listed in
-the following table:
+be zero if no contact has occurred.  The first 8 values of these columns are
+listed in the following table.
 
 +-------+----------------------------------------------------+----------------+
 | Index | Value                                              | Units          |
@@ -248,6 +248,14 @@ the following table:
 |     8 | Radius :math:`r` of atom                           | distance units |
 +-------+----------------------------------------------------+----------------+
 
+If a granular submodel calculates additional contact information (e.g. the
+heat submodels calculate the amount of heat exchanged), these quantities
+are appended to the end of this array. First, any extra values from the
+normal submodel are appended followed by the damping, tangential, rolling,
+twisting, then heat models. See the descriptions of granular submodels in
+the :doc:`pair granular <pair_granular>` page for information on any extra
+quantities.
+
 None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
 No parameter of this fix can be used with the *start/stop* keywords of the
 :doc:`run <run>` command. This fix is not invoked during :doc:`energy

doc/src/fix_wall_gran_region.rst

@@ -243,10 +243,10 @@ uninterrupted fashion.
    with a different region ID.
 
 If the :code:`contacts` option is used, this fix generates a per-atom array
-with 8 columns as output, containing the contact information for owned
+with at least 8 columns as output, containing the contact information for owned
 particles (nlocal on each processor). All columns in this per-atom array will
-be zero if no contact has occurred. The values of these columns are listed in
-the following table:
+be zero if no contact has occurred.  The first 8 values of these columns are
+listed in the following table.
 
 +-------+----------------------------------------------------+----------------+
 | Index | Value                                              | Units          |
@@ -269,6 +269,14 @@ the following table:
 |     8 | Radius :math:`r` of atom                           | distance units |
 +-------+----------------------------------------------------+----------------+
 
+If a granular submodel calculates additional contact information (e.g. the
+heat submodels calculate the amount of heat exchanged), these quantities
+are appended to the end of this array. First, any extra values from the
+normal submodel are appended followed by the damping, tangential, rolling,
+twisting, then heat models. See the descriptions of granular submodels in
+the :doc:`pair granular <pair_granular>` page for information on any extra
+quantities.
+
 None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
 No parameter of this fix can be used with the *start/stop* keywords of the
 :doc:`run <run>` command. This fix is not invoked during :doc:`energy

doc/src/pair_granular.rst

@@ -40,6 +40,9 @@ Examples
    pair_style granular
    pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 heat area 0.1
 
+   pair_style granular
+   pair_coeff * * mdr 5e6 0.4 1.9e5 2.0 0.5 0.5 tangential linear_history 940.0 0.0 0.7 rolling sds 2.7e5 0.0 0.6 damping none
+
 Description
 """""""""""
 
@@ -82,6 +85,7 @@ and their required arguments are:
 3. *hertz/material* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`
 4. *dmt* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`, :math:`\gamma`
 5. *jkr* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`, :math:`\gamma`
+6. *mdr* : :math:`E`, :math:`\nu`, :math:`Y`, :math:`\Delta\gamma`, :math:`\psi_b`, :math:`e`
 
 Here, :math:`k_n` is spring stiffness (with units that depend on model
 choice, see below); :math:`\eta_{n0}` is a damping prefactor (or, in its
@@ -162,6 +166,143 @@ initially will not experience force until they come into contact
 experience a tensile force up to :math:`3\pi\gamma R`, at which point they
 lose contact.
 
+The *mdr* model is a mechanically-derived contact model designed to capture the
+contact response between adhesive elastic-plastic particles into large deformation.
+The theoretical foundations of the *mdr* model are detailed in the
+two-part series :ref:`Zunker and Kamrin Part I <Zunker2024I>` and
+:ref:`Zunker and Kamrin Part II <Zunker2024II>`. Further development
+and demonstrations of its application to industrially relevant powder
+compaction processes are presented in :ref:`Zunker et al. <Zunker2025>`.
+
+The model requires the following inputs:
+
+   1. *Young's modulus* :math:`E > 0` : The Young's modulus is commonly reported
+   for various powders.
+
+   2. *Poisson's ratio* :math:`0 \le \nu \le 0.5` : The Poisson's ratio is commonly
+   reported for various powders.
+
+   3. *Yield stress* :math:`Y \ge 0` : The yield stress is often known for powders
+   composed of materials such as metals but may be unreported for ductile organic
+   materials, in which case it can be treated as a free parameter.
+
+   4. *Effective surface energy* :math:`\Delta\gamma \ge 0` : The effective surface
+   energy for powder compaction applications is most easily determined through its
+   relation to the more commonly reported critical stress intensity factor
+   :math:`K_{Ic} = \sqrt{2\Delta\gamma E/(1-\nu^2)}`.
+
+   5. *Critical confinement ratio* :math:`0 \le \psi_b \le 1` : The critical confinement
+   ratio is a tunable parameter that determines when the bulk elastic response is
+   triggered. Lower values of :math:`\psi_b` delay the onset of the bulk elastic
+   response.
+
+   6. *Coefficient of restitution* :math:`0 \le e \le 1` : The coefficient of
+   restitution is a tunable parameter that controls damping in the normal direction.
+
+.. note::
+
+   The values for :math:`E`, :math:`\nu`, :math:`Y`, and :math:`\Delta\gamma` (i.e.,
+   :math:`K_{Ic}`) should be selected for zero porosity to reflect the intrinsic
+   material property rather than the bulk powder property.
+
+The *mdr* model produces a nonlinear force-displacement response, therefore the
+critical timestep :math:`\Delta t` depends on the inputs and level of
+deformation. As a conservative starting point the timestep can be assumed to be
+dictated by the bulk elastic response such that
+:math:`\Delta t = 0.35\sqrt{m/k_\textrm{bulk}}`, where :math:`m` is the mass of
+the smallest particle and :math:`k_\textrm{bulk} = \kappa R_\textrm{min}` is an
+effective stiffness related to the bulk elastic response.
+Here, :math:`\kappa = E/(3(1-2\nu))` is the bulk modulus and
+:math:`R_\textrm{min}` is the radius of the smallest particle.
+
+.. note::
+
+   The *mdr* model requires some specific settings to function properly,
+   please read the following text carefully to ensure all requirements are
+   followed.
+
+The *atom_style* must be set to *sphere 1* to enable dynamic particle
+radii. The *mdr* model is designed to respect the incompressibility of
+plastic deformation and inherently tracks free surface displacements
+induced by all particle contacts. In practice, this means that all particles
+begin with an initial radius, however as compaction occurs and plastic
+deformation is accumulated, a new enlarged apparent radius is defined to
+ensure that that volume change due to plastic deformation is not lost.
+This apparent radius is stored as the *atom radius* meaning it is used
+for subsequent neighbor list builds and contact detection checks. The
+advantage of this is that multi-neighbor dependent effects such as
+formation of secondary contacts caused by radial expansion are captured
+by the *mdr* model. Setting *atom_style sphere 1* ensures that updates to
+the particle radii are properly reflected throughout the simulation.
+
+.. code-block:: LAMMPS
+
+   atom_style sphere 1
+
+Newton's third law must be set to *off*. This ensures that the neighbor lists
+are constructed properly for the topological penalty algorithm used to screen
+for non-physical contacts occurring through obstructing particles, an issue
+prevalent under large deformation conditions. For more information on this
+algorithm see :ref:`Zunker et al. <Zunker2025>`.
+
+.. code-block:: LAMMPS
+
+   newton off
+
+The damping model must be set to *none*. The *mdr* model already has a built
+in damping model.
+
+.. code-block:: LAMMPS
+
+   pair_coeff * * mdr 5e6 0.4 1.9e5 2 0.5 0.5 damping none
+
+The definition of multiple *mdr* models in the *pair_style* is currently not
+supported. Similarly, the *mdr* model cannot be combined with a different normal
+model in the *pair_style*. Physically this means that only one homogenous
+collection of particles governed by a single *mdr* model is allowed.
+
+The *mdr* model currently only supports *fix wall/gran/region*, not
+*fix wall/gran*. If the *mdr* model is specified for the *pair_style*
+any *fix wall/gran/region* commands must also use the *mdr* model.
+Additionally, the following *mdr* inputs must match between the
+*pair_style* and *fix wall/gran/region* definitions: :math:`E`,
+:math:`\nu`, :math:`Y`, :math:`\psi_b`, and :math:`e`. The exception
+is :math:`\Delta\gamma`, which may vary, permitting different
+adhesive behaviors between particle-particle and particle-wall interactions.
+
+.. note::
+
+   The *mdr* model has a number of custom *property/atom* and *pair/local* definitions that
+   can be called in the input file. The useful properties for visualization
+   and analysis are described below.
+
+In addition to contact forces the *mdr* model also tracks the following
+quantities for each particle: elastic volume change, average normal
+stress components, total surface area involved in
+contact, and individual contact areas. In the input script, these quantities are
+initialized by calling *run 0* and can then be accessed using subsequent *compute*
+commands. The last *compute* command uses *pair/local p13* to calculate the pairwise
+contact areas for each active contact in the *group-ID*. Due to the use of an apparent
+radius in the *mdr* model, the keyword/arg pair *cutoff radius* must be specified for
+*pair/local* to properly detect existing contacts.
+
+.. code-block:: LAMMPS
+
+   run 0
+   compute ID group-ID property/atom d_Velas
+   compute ID group-ID property/atom d_sigmaxx
+   compute ID group-ID property/atom d_sigmayy
+   compute ID group-ID property/atom d_sigmazz
+   compute ID group-ID property/atom d_Acon1
+   compute ID group-ID pair/local p13 cutoff radius
+
+.. note::
+
+   The *mdr* model has two example input scripts within the
+   *examples/granular* directory. The first is a die compaction
+   simulation involving 200 particles named *in.tableting.200*.
+   The second is a triaxial compaction simulation involving 12
+   particles named *in.triaxial.compaction.12*.
 ----------
 
 In addition, the normal force is augmented by a damping term of the
@@ -674,7 +815,10 @@ supported are:
 2. *radius* : :math:`k_{s}`
 3. *area* : :math:`h_{s}`
 
-If the *heat* keyword is not specified, the model defaults to *none*.
+If the *heat* keyword is not specified, the model defaults to *none*. All
+heat models calculate an additional pairwise quantity accessible by the
+single() function (described below) which is the heat conducted between the
+two particles.
 
 For *heat* *radius*, the heat
 :math:`Q` conducted between two particles is given by
@@ -789,17 +933,25 @@ The single() function of these pair styles returns 0.0 for the energy
 of a pairwise interaction, since energy is not conserved in these
 dissipative potentials.  It also returns only the normal component of
 the pairwise interaction force.  However, the single() function also
-calculates 13 extra pairwise quantities.  The first 3 are the
+calculates at least 13 extra pairwise quantities.  The first 3 are the
 components of the tangential force between particles I and J, acting
 on particle I.  The fourth is the magnitude of this tangential force.
 The next 3 (5-7) are the components of the rolling torque acting on
 particle I. The next entry (8) is the magnitude of the rolling torque.
 The next entry (9) is the magnitude of the twisting torque acting
 about the vector connecting the two particle centers.
 The next 3 (10-12) are the components of the vector connecting
-the centers of the two particles (x_I - x_J). The last quantity (13)
-is the heat flow between the two particles, set to 0 if no heat model
-is active.
+the centers of the two particles (x_I - x_J). If a granular submodel
+calculates additional contact information (e.g. the heat submodels
+calculate the amount of heat exchanged), these quantities are appended
+to the end of this list. First, any extra values from the normal submodel
+are appended followed by the damping, tangential, rolling, twisting, then
+heat models. See the descriptions of specific granular submodels above
+for information on any extra quantities. If two or more models are
+defined by pair coefficients, the size of the array is set by the
+maximum number of extra quantities in a model but the order of quantities
+is determined by each model's specific set of submodels. Any unused
+quantities are zeroed.
 
 These extra quantities can be accessed by the :doc:`compute pair/local <compute_pair_local>` command, as *p1*, *p2*, ...,
 *p12*\ .
@@ -870,10 +1022,32 @@ solids. Proc. R. Soc. Lond. A, 324(1558), 301-313.
 
 .. _DMT1975:
 
-**Derjaguin et al, 1975)** Derjaguin, B. V., Muller, V. M., & Toporov,
+**(Derjaguin et al, 1975)** Derjaguin, B. V., Muller, V. M., & Toporov,
 Y. P. (1975). Effect of contact deformations on the adhesion of
 particles. Journal of Colloid and interface science, 53(2), 314-326.
 
+.. _Zunker2024I:
+
+**(Zunker and Kamrin, 2024)** Zunker, W., & Kamrin, K. (2024).
+A mechanically-derived contact model for adhesive elastic-perfectly
+plastic particles, Part I: Utilizing the method of dimensionality
+reduction. Journal of the Mechanics and Physics of Solids, 183, 105492.
+
+.. _Zunker2024II:
+
+**(Zunker and Kamrin, 2024)** Zunker, W., & Kamrin, K. (2024).
+A mechanically-derived contact model for adhesive elastic-perfectly
+plastic particles, Part II: Contact under high compaction—modeling
+a bulk elastic response. Journal of the Mechanics and Physics of Solids,
+183, 105493.
+
+.. _Zunker2025:
+
+**(Zunker et al, 2025)** Zunker, W., Dunatunga, S., Thakur, S.,
+Tang, P., & Kamrin, K. (2025). Experimentally validated DEM for large
+deformation powder compaction: mechanically-derived contact model and
+screening of non-physical contacts. engrXiv.
+
 .. _Luding2008:
 
 **(Luding, 2008)** Luding, S. (2008). Cohesive, frictional powders: